Two-dimensional spectroscopy on a THz quantum cascade structure

1 Introduction

The terahertz (THz) spectral region has a long standing interest for a wide range of applications such as imaging [1], telecommunication [2], [3], spectroscopy [4], as well as thickness measurements [5]. It is also indispensable in current fundamental research, e.g., in the study of superconductivity [6], [7], quasiparticles such as magnons [8], heavy fermions [9] and in dynamically induced phase transitions [10]. However, a lack of compact, powerful and efficient THz radiation sources has hindered a wide-spread deployment of such applications. A promising candidate is the THz quantum cascade laser (QCL) [11], a unipolar device based on intersubband transitions within the conduction band of an electrically pumped semiconductor heterostructure. The state-of-the-art devices have emission bandwidths of 1–5 THz [12], [13], Watt-level power on a single mode [14] and operate up to 210 K in pulsed mode [15]. Recent developments in the field of optical sensing by means of frequency combs in the THz [16] and mid-infrared [17] spectral regions has gained a lot of attraction. This is because such a spectroscopic technique (dual-comb spectroscopy [18]) provides high resolution and fast acquisition speeds without moving optical or optomechanical parts. One important aspect for spectroscopy applications is to increase the spectral coverage of the frequency combs and to control their noise properties. This requires understanding of the frequency comb formation mechanism. Currently, there are two proposed mechanisms which have different origin. The first relies on four-wave mixing (FWM) which arises from the third-order susceptibility (χ⁽³⁾) of the gain material itself [17], [19]. The second one is associated with spatial hole burning (SHB) in a Fabry–Perot (FP) cavity. SHB is essentially a dynamical gain grating formed by gain saturation within the laser cavity, which enables mode coupling and can be seen as an additional, effective χ⁽³⁾ nonlinearity within the cavity. Recent models [20], [21], [22], that take SHB into account, reproduced spontaneous comb formation with the experimentally observed [23] linear frequency chirp in mid-infrared QCLs. However, in a ring cavity laser, SHB can be completely eliminated, and hence only FWM from the gain medium is present and responsible for comb formation [24]. The interplay of SHB and FWM makes comb formation in FP cavities a highly nonlinear phenomenon and offer a rich playground for cavity and laser structure designs. Since both mechanisms rely on χ⁽³⁾ processes, it is essential to investigate the fundamental nonlinear material properties. In addition, devices such as quantum cascade detectors [25], saturable absorbers [26], and nonlinear conversion media [27] with giant nonlinearities rely on the same operation principle of intersubband transitions in multi-quantum-well structures.

To our knowledge, no investigation of nonlinear properties has been carried out on quantum cascade structures (QCSs) at THz frequencies with two-dimensional spectroscopy, although single and double quantum wells has been studied extensively [28], [29], [30], [31]. Since several systems rely on cascaded intersubband transitions, it is of vital importance to analyze and explore the nonlinear optical properties and quantum-coherent evolution in these type of structures, since modelocking of QCLs due to self-induced transparency has been predicted [32].

The focus of our investigation is a THz QCS and its nonlinear intersubband properties, as well as its coherent quantum evolution. We perform coherent collinear two-dimensional THz (2D-THz) spectroscopy on a passive QCS with the aim to determine the ultrafast nonlinear response and relaxation processes originating from the intersubband states. The measurements in this work are performed on an unbiased QCS because this structure has also the potential to be used as a saturable absorber, which can be monolithically integrated into the QCL cavity. Since this multi-quantum-well system is far from a trivial system and involves many coupled electron states, an advanced modeling is required to fully understand and interpret the experimental results. Therefore, we go beyond the commonly employed models based on classical [33] or few-level Maxwell–Bloch equations with phenomenological scattering and dephasing rates [26], [30], [34], [35]. We model our experiment using a time-resolved density matrix model formalism which takes into account all the relevant intersubband states, as well as the relevant scattering mechanisms (except electron–electron interactions) treated with a master equation in Lindblad form [36]. This allows for ab initio dynamical simulations of the quantum evolution of the full system density matrix, containing all interactions between coherences and populations mediated by the various scattering mechanisms. Explicitly including the time-dependent electric field pulses on a semiclassical level, allows for calculating the nonlinear optical response nonperturbatively.

2 Materials and methods

2.1 Sample and measurements

In order to gain a detailed insight into the rich nonlinear features of a multilevel intersubband system, we perform 2D-THz spectroscopy on an unbiased THz QCL structure [37] which is designed to lase at 1.2 THz. Additional reason to perform the 2D spectroscopy on an unbiased structure is to keep the system complexity and validate the developed model. Figure 1(a) shows the schematic of the experiment (performed at 10 K), where two phase-locked THz pulses, A and B, are generated by optical rectification in ZnTe and GaP, respectively. These pulses arrive at the sample with a relative time delay τ. The investigated active region is 10 µm in height, 800 µm long, 400 µm wide and is grown on a 500 µm thick n-doped GaAs substrate. We lap the substrate down to 190 µm thickness in order to increase the efficiency of the light–matter interaction, polish the substrate to optical quality and coat it with 200 gold for better wave guiding. In this way, the active region and 190-µm-thick substrate are embedded in to a 200-µm-thick single-plasmon wave guide. The optical rectification in the two crystals is realized with 120 fs infrared pulses originating from an amplified Ti:Sapphire laser operating at 1 kHz repetition rate with the center wavelength of 800 nm. In order to obey the intersubband selection rules, the THz pulses are polarized along the growth direction of the QCS. The transmitted THz pulses, which carry the imprint of the induced nonlinear sample response, are detected by a phase-resolved linear detector based on free-space electro-optic (EO) sampling, hence measuring the electric field as a function of detection time t. Figure 1(b) shows the net electric field of two incident THz pulses in air separated by a delay time τ = 2.5 ps. The corresponding spectrum of a single excitation pulse is shown in Figure 1(c) and covers a bandwidth of ∼2.5 THz. In order to separate the nonlinear sample response, which arises due to the presence of the combined pulses from the single excitation pulses, a sequence of measurements is performed for each delay τ where the sample interacts with the electric field of pulse A (E_A) and pulse B (E_B) separately, as well as with the two pulses superimposed (E_AB). The nonlinear system response is extracted as

(1)ENL(t,τ)=EAB(t,τ)−EA(t,τ)−EB(t).

Figure 1:

(a) Schematic of the experimental setup showing a collinear pulse geometry used to perform a 2D-THz spectroscopy on quantum cascade structure (QCS). (b) THz electric filed amplitude of pulses A and B separated by a time delay τ = 2.5 ps before coupling in to the quantum cascade wave guide. Black curve represent the experimental pulses while the red curve shows the pulses used in the simulation. (c) Spectra of a single experimental and simulation pulse. (d) Nonlinear signal extracted from the pulse sequence measurements: ENL(t,τ)=EAB(t,τ)−EA(t,τ)−EB(t). The green and orange dashed lines indicate the position of the driving fields A and B, respectively. The dashed horizontal line is a guide to the eye for τ = 0 ps.

This results in a two-dimensional (2D) map of E_NL (t, τ) as shown in Figure 1(d). Since the measurement is performed in a collinear geometry, all nonlinear signals are simultaneously present in the same measured 2D map. Different orders of the nonlinear sample response can be separated by Fourier transforming into the frequency domain, along the delay axis τ and sampling axis t. The corresponding 2D frequency map as a function of the excitation frequency ν_τ and detection frequency ν_t is shown in Figure 2(a). The collinear geometry of the experiment makes the frequency vectors of the two pulses correspond to vectors with specific directions in 2D frequency space [38], [39]; the frequency vector ν_B of pulse B is parallel to the detection frequency axis, whereas the frequency vector ν_A of pulse A is along the diagonal. This is due to the fact that pulse A is moving with respect to the pulse B (which is kept fixed in time), and thus the wave fronts are tilted by 45° with regards to the sampling axis. Thus, the location of the nonlinear signals in the 2D frequency map can be expressed as a linear combination of frequency vectors ν_A and ν_B, whose lengths are given by the frequencies of the two pulses. The experimental results in Figure 2(a) show multiple nonlinear resonances, marked by circles and labeled as PP₁, PP₂ and PP₃. These correspond to pump-probe (PP) signals and are generated by the following interaction sequence: Two field interactions of pulse A, creating a nonequilibrium level occupation, followed by a time-delayed single field interaction of pulse B, emitting the nonlinear signal. The nonlinear signals are then detected with EO sampling. Since the detected signal is generated by three light–matter interactions, the origin of this signal is due to the strong χ⁽³⁾ of the electronic intersubband system, and the location in the 2D frequency map is given by the frequency vectors ν_AB = ν_A − ν_A + ν_B. We note that any higher odd order could contribute to this signal, but for the experimental field strength of ∼3 kV/cm we estimate the χ⁽⁵⁾ contribution to be at least three orders of magnitude smaller. The signals along the diagonal in Figure 2(a) are also PP signals, but where the order of pulses A and B is reversed. In this way, there are two interactions with pulse B and one with pulse A. These PP signals are located along the diagonal given by ν_BA = ν_B − ν_B + ν_A. Since the electric field strength of pulse A is larger than that of pulse B, the PP signals are smeared out along the diagonal due to the strong system perturbation caused by pulse A. The elongated shape of the B–B–A pulse sequence makes the PP signals almost impossible to distinguish individually. Hence, in the following analysis, we will focus only on the PP signals at zero excitation frequency resulting from A–A–B pulse sequences.

Figure 2:

(a) 2D spectrum of the quantum cascade structure (QCS) obtained by performing a 2D Fourier transform of Figure 1(d) along the sampling and delay axis. Signals along the detection frequency axis ν_t at zero excitation frequency are pump-probe (A_pumpB_probe) signals. Pump-probe (B_pumpA_probe) signals along the diagonal are smeared out due to the strong system saturation with pulse A. (b) Simulated E_NL induced by the time-dependent nonlinear current density j_NL evaluated in the density matrix model. All relevant states of QCS are taken into account. Color-coded circles in (a) and (b) mark pump-probe signals which originates from the engineered states. The corresponding energy levels contributing to these A_pumpB_probe signals are shown in Figure 3 with similarly color-coded arrows.

2.2 Simulation model

For simulating the optical response of the structure, we employ a time-resolved density matrix model extended from a study by Kiršanskas et al. [36] to include arbitrary pulse shapes, mean-field potential, as well as a phenomenological electron–electron dephasing time [40]. The elements ρ_ij of the full density matrix are evolved according to the Liouville–von-Neumann equation in Lindblad form

(2)ρ˙ij=1iℏ[HS,ρ]ij+∑α,kl(LikαρklLkjα−12LikαLklαρlj−12ρikLklαLljα)≡1iℏ[HS,ρ]ij+∑α,klRijklαρkl

where HS=H0+HEM(t) denotes the Hamiltonian of the intersubband system (H₀) including electromagnetic field (E(t))

(3)HEM(t)=−|e|zˆE(t)

in Lorenz gauge, Lijα denote Lindblad jump operators as defined in [36], and the equivalent scattering tensors Rijklα contain couplings between the intersubband system and multiple baths labeled by α, computed within Fermi’s golden rule. Due to the computational load of resolving explicitly the in-plane momenta of the electrons, the density matrix indices and the rates Rijklα are thermal averages with the electron temperature set to the lattice one. For details, see in the study by Kiršanskas et al. [36]. In our simulations, we include longitudinal optical (LO) phonon scattering and pure dephasing implemented as in the study by Gordon and Majer [40], as well as scattering with ionized impurities and interface roughness. The coherent system evolution is governed by the first term H_S, which includes the interaction with the electric field pulses. Since the model is solved in the time domain, the nonlinear response is included to infinite order, and it is thus able to capture any nonlinear effects that may occur in the experiment. The second term Rijklα is responsible for decoherence and dephasing. In addition, this full treatment of scattering on Lindblad form (without the secular approximation), allows the baths to generate coherences as well, allowing for a full treatment of coherent effects such as scattering assisted tunneling, important when the levels are closely spaced as in the present structure. As a basis, we use the Wannier–Stark states of the unbiased structure (a negligible bias of 0.1 meV/period is necessary to break the periodic symmetry of the structure). Due to the computational cost of calculating the scattering tensors Rijklα, these are calculated without time-dependent fields, and assumed to remain constant. This will make a quantitative difference to the exact scattering rates in the system while the pulses act but will not affect the relaxation behavior once the pulses have subsided or affect significantly the coherent evolution imposed by the pulse, since its duration is short compared to the subband lifetimes. The first simulation step is to reach the self-consistent equilibrium state by iterating the Poisson equation and the evolution of Eq. (2) without external fields. This results in the basis states and potential profile (including the mean-field from the electrons and ionized dopants) shown in Figure 3. The electromagnetic field pulses (parametrized to be similar to the experimental ones and shown in Figure 1(b)) are then applied with varying delay according to the experimental parameters. This results in a response of the density matrix ρ_ij(t), and thus the current density according to

(4)〈j〉=Ne〈pˆ〉m*=Nem*Tr{ρpˆ}=Nem*∑klρklplk

where N is the 3D electron density, m* the effective mass and pˆ is the momentum operator. Since p_ii = 0, the current is given solely by the off-diagonal elements of the density matrix, i.e., the coherences. In order to allow for a direct comparison to the experiment, which measures the electric field emitted from the structure rather than the induced polarization current, the induced electric field is evaluated in the frequency domain as

(5)E(ω)=−iμ0j(ω)ω,

where μ₀ denotes the vacuum permeability. The nonlinear contribution to the electric field E_NL(ω) can be evaluated via the Fourier transform of

(6)jNL(t,τ)=jAB(t,τ)−jA(t−τ)−jB(t)+jbg(t)

where j_AB(t, τ) is the current response of both pulse B (with fixed delay) and A (with delay τ), and j_A/B(t) are the single-pulse current responses. Here, we also subtract each of the current contributions on the right hand side with the background current j_bg(t), which removes any small transients and DC currents coming from the small added bias and the finite number of periods included in the simulation.

Figure 3:

(a) Calculated band diagram of the quantum cascade structure at zero bias that shows the conduction band profile E_C and the six relevant energy states which contribute to the nonlinear signals. (b) The transitions which give rise to the same pump-probe signal. (c) Three of four levels (ladder-type system) which have a dominant contribution to the PP₁ signal and are responsible for nontrivial relaxation decay after the pump drives the system far from equilibrium.

By varying the delay between the pulses over the same range as used in the experiment, we obtain a similar 2D map of jNL(ντ,νt), which is Fourier transformed to ENL(ντ,νt) using Eq. (5). As seen in Figure 2(b), the simulation replicate qualitatively all experimental features. Both experiment and simulation show multiple resonances along the detection frequency axis νt at zero excitation frequency and a strong smeared-out PP signals along the diagonal. We verified, with simulations, that by reducing the A pulse intensity the PP signals along the diagonal are more localized in the 2D frequency space. By increasing the barrier AlAs fraction by ∼3% (from 0.10 to 0.103) from its nominal value, which results in a change in band offset which is small compared to the discrepancy between literature values [41], [42], a good agreement for the absolute frequency positions of the PP resonances PP₁, PP₂ and PP₃ is found with respect to the experiment. This confirms the computed band structure shown in Figure 3, and allows us to identify the transitions responsible for each PP peak as indicated by the arrows in Figure 3.

3 Results

The origin of the PP signals at detection frequencies of 0.2–0.5, 0.7–0.9 and 1.2–1.5 THz can be deduced from the energy band diagram of the QCS shown in Figure 3(a). Due to the low doping (2.6 × 10¹⁶ cm⁻³), only the lowest three energy levels are significantly populated at 10 K (in equilibrium). Three groups of intersubband transitions with energies close to the different PP signals are identified and marked in Figure 3(a) as black, blue, and green arrows, respectively. For the sake of clarity we have tabulated these transitions in Figure 3(b). Since the calculated band diagram at zero bias is highly sensitive to a number of factors, such as the mean-field potential, barrier height, and effective masses, the slight frequency deviation of a few tens of GHz of the PP signals between the experiment and the simulation is expected. The spectral amplitude of the pulses has a minimum at the frequency of the PP₃ signal, as shown in Figure 1(c). The higher energy states (E₄, E₅ and E₆) which contributes to this signal are insignificantly populated. Therefore the PP₃ signal strength is one order of magnitude lower than the other PP signal strengths, even though the oscillator strengths in the PP₃ transitions are in general higher.

Next, we study the time-resolved dynamics of individual nonlinear PP signals. This is realized by multiplying the 2D spectrum with a 2D Gaussian filter and performing the inverse Fourier transforms of the individual frequency components back into the time domain. The 2D Gaussian filter has a spectral bandwidth with full width half maximum as shown in Figure 2 (color-coded circles). The results are shown in Figure 4 (a), (b) and (c). As multiple states contribute to the same PP signal, it is difficult to extract the incoherent state life time of individual energy levels directly from the time domain. We therefore compute a spectrogram by performing a Fourier transform along the sampling time axis t for each delay time τ. The corresponding spectrograms are shown in Figure 4 (a_i), (b_i) and (c_i), which describe the time evolution of the frequency-resolved (nonlinear) state populations of the respective PP signals. In order to extract the population decay time of individual PP signals, cuts in the spectrograms along the delay time axis τ (dashed lines in Figure 4 (a_i), (b_i), and (c_i)) are extracted and shown in Figure 4 (a_ii), (b_ii), and (c_ii). The spectrally resolved PP signals show a general decay trend with an additional nontrivial oscillatory behavior. The overall decay times of the PP₁, PP₂ and PP₃ signals are around 4, 8 and 6 ps, respectively. Note that, since several transitions contribute to each PP signal, this method does not provide decay times of individual transitions. Rather, it provides the average rate of decay of the nonlinear response through all “incoherent channels” with transition energies corresponding to the respective PP signal frequency. By incoherent channel, we mean essentially the component of the population decay which is independent of the delay time τ (thus showing no significant frequency components at the filtered ν_τ = 0). For clarity, the average spectral decay is computed and is visible as a black curve in Figure 4 (a_ii), (b_ii), and (c_ii). Additionally, we have verified the presence of coherent oscillations by computing the Fourier transform of average signals in Figure 4 (a_ii), (b_ii) and (c_ii). The corresponding oscillations appear at frequencies of 0.4, 0.33 and 0.22 THz for PP₁, PP₂ and PP₃, respectively.

Figure 4:

(a), (b) and (c) are inverse Fourier transforms of pump-probe signals PP₁, PP₂ and PP₃ from Figure 2(a) and show the time evolution of each nonlinear signal. Since multiple states contribute to the same pump-probe signal, it is more convenient to study the system dynamics by looking at the time-resolved spectrum along the system evolution axis. This is realized by calculation of the Fourier transform of (a), (b) and (c) along the sampling time axis t for every time delay τ. In this way calculated spectrograms are shown in (a_i), (b_i) and (c_i). Extracted cuts at fixed detection frequencies along delay time τ display spectral time evolution and are shown in (a_ii), (b_ii) and (c_ii). They are also indicated with color-coded dashed lines in (a_i), (b_i) and (c_i). The black curve in (a_ii),(b_ii) and (c_ii) shows the average detection frequency evolution spectrum.

In contrast to previous studies on similar, but simpler, quantum well samples, our system shows a population decay as well as non-trivial oscillation (see Figure 4). As already mentioned, Raab et al. [26] investigated a multi-quantum-well system which contains only two levels. Such a system will possess only two PP signals (A_pumpB_probe and B_pumpA_probe). In that two-level system, the expected monotonic, exponential decay of the exited state population was observed. Systems with more than two levels, such as that studied by Kuehn et al. [31], can possess several PP signals at zero excitation frequency, which can be separated in frequency domain provided the energy level separation is larger than the frequency resolution and the level broadening. In these cases, several transitions also contribute to the same PP signals. However, due to the fast dephasing times in the structure in a study by Kuehn et al. [31], with level separations above the LO phonon energy, pulsations in population differences, such as those we observe in Figure 4, could not be seen. We attribute the observed nontrivial oscillatory decay to the coherent carrier exchange among the states. This is enabled by the existence of long coherence times as well as several populated intersubband states. Such a coherent population exchange is driven by the coherent evolution of coherences sharing common state, which are excited by the strong THz pulses. The fact that we can observe it experimentally means that it is associated with an electrical dipole, which requires a non-negligible dipole matrix element between the involved subbands. The existence of such nontrivial incoherent decay is supported by our simulations. Figure 5 show the same analysis on computed 2D spectrum (PP signals from Figure 2(b)). We can identify qualitatively oscillations of the average spectrum in Figure 4 (a_ii), (b_ii) and (c_ii) as observed in the experiment. A detailed discussion on observed coherent oscillations in experiment and theory is elaborated in the following section.

Figure 5:

(a), (b) and (c) are inverse Fourier transforms of pump-probe signals PP₁, PP₂ and PP₃ from Figure 2(b) and show the time evolution of each nonlinear signal. Since multiple states contribute to the same pump-probe signal, it is more convenient to study the system dynamics by looking at the time-resolved spectrum along the system evolution axis. This is realized by calculation of the Fourier transform of (a), (b) and (c) along the sampling time axis t for every time delay τ. In this way calculated spectrograms are shown in (a_i), (b_i) and (c_i). Extracted cuts at fixed detection frequencies along delay time τ display spectral time evolution and are shown in (a_ii), (b_ii) and (c_ii). They are also indicated with color-coded dashed lines in (a_i), (b_i) and (c_i).The black curve in (a_ii), (b_ii) and (c_ii) shows the average detection frequency evolution spectrum.

4 Discussion

In order to find a microscopic explanation to the nonmonotonic decay patterns seen in Figure 4 (a_ii–c_ii) and Figure 5 (a_ii–c_ii), we have investigated the nonlinear contributions to the density matrix elements contributing to the respective PP peaks in Figure 2(b). The nonlinear density matrix of the system provides access to the full time evolution of the nonlinear system. The experimental and computed signals at zero excitation frequency include only the indirect contributions to the decay, i.e. not involving the coherences attributed to the respective PP signals. By including the signals around zero excitation frequency, coherent oscillation effects can be taken into account, which influence the overall population decay. In our model, the nonlinear contribution to the diagonal density matrix elements, which describe the populations, show an oscillating behavior after the excitation pulse has acted. Furthermore, the population oscillations in levels 2 and 3 are out of phase, which indicates a coherent population exchange between the levels. This phenomenon is known from atomic physics, where coherent light fields interact with an atomic system, commonly with three levels. Such a three level system can be either of Λ, V or ladder type [43]. Depending on the system parameters and the system drive, a quantum interference between the excitation pathways will control the optical response of the system and can result in effects such as electromagnetically induced transparency or coherent population trapping (transfer). In our case, the first subband plays a negligible role for the coherences of the system. Thus, the first three excited states are almost decoupled from the ground state and forms a ladder-type system where the population can be coherently exchanged among the levels, as schematically shown in Figure 3(c). Thus, since all PP peaks involve transitions with at least one of levels 2 and 3, they all show similar oscillation patterns in the delay time. The exact behavior also depends on the scattering rates, which are underestimated in the simulation model since electron–electron scattering is neglected. We have checked the calculations by adding an extra pure dephasing channel [40], which does not alter the general trends as it does not increase the transition rates between populations.

Besides the microscopic origin of the observed coherent oscillations, we would like to discuss whether the usage of the single-plasmon wave guide influences the THz pulses (dispersion and nonlinear generation along the entire wave guide) and hence the experimental results. We have computed the losses in the wave guide which are around 44 cm⁻¹ and mainly are due to the highly doped substrate. Such relatively high loss results in rapid absorption of the in-coupled THz pulses. In this way, the measured nonlinear signal originates from the first ∼200 μm of the wave guide. From the measured transmission, we can observe a slight broadening of the THz pulse which can be explained by the group velocity dispersion along the wave guide.

Naturally, the simulation model does not include wave guide dispersion or loss, which is equivalent to assuming a thin wave guide slab. However, since the nonlinear signal mainly originates from the facet of the wave guide as discussed above, the qualitative comparison between the experiment and simulation is still valid. Thus, we do not expect that including a more realistic wave guide, by e.g. performing several simulations with decreasing THz pulse amplitudes and obtaining an average 2D map, would alter the quality of the observed phenomenon (coherent oscillation).

5 Conclusion

In conclusion, we have conducted the first experimental and theoretical 2D-THz spectroscopic study on a QCS at THz frequencies. We observe multiple nonlinear (PP) signals originating from the engineered intersubband states, from which we determine the average incoherent intersubband lifetimes. Moreover, we observe nontrivial system dynamics attributed to field-induced coherences, the so-called coherent population exchange. We find incoherent lifetimes to be several times smaller than the typical gain recovery times in THz QCLs [44], [45], [46], as well as saturation times to be shorter than the QCL gain saturation time. Under these circumstances, one can consider such a structure as a monolithically integrated saturable absorber into a QCL cavity in order to achieve a passive mode locking. However, this is complicated due to the observed nontrivial dynamics, as the interaction between all density matrix elements must be considered. Further investigation of this QCS with respect to bias, field strength and temperature dependence is needed in order to explore its full potential as saturable absorber, ultrafast pulse shaper, frequency converter, as well as its FWM efficiency.