1. Introduction

A central feature of the quantum cascade laser (QCL) is that the optical properties of the nanostructured gain medium can be custom-tailored by band structure engineering. This is not only relevant for obtaining the desired gain characteristics, but can also be used to design the nonlinear optical properties. The resulting giant artificial nonlinearities have recently been exploited for frequency conversion structures, where highly efficient mid-infrared (MIR) QCLs are utilized as a pump to generate radiation in hardly accessible frequency regimes by using nonlinear optical processes. One of the most promising applications is terahertz (THz) difference frequency generation (DFG). This approach enables the generation of coherent THz radiation at room temperature based on compact semiconductor devices [1–9], avoiding the limitations of conventional THz QCLs which are currently restricted to operating temperatures below 200K [10]. Pulsed operation with up to 1.4mW THz peak power has been demonstrated [8], and recently, also continuous wave (cw) room temperature operation with 3μW output power has been achieved [8]. Furthermore, wide wavelength tunability has been demonstrated [4, 5, 9], as needed for various applications in sensing and imaging. A maximum tuning range of 1.2 − 5.9 THz has been obtained, with a peak output power between 5 and 90μW depending on the frequency [9]. This remarkable improvement in THz output power and tuning range has mainly been achieved by optimizing the active region design and using Cherenkov waveguide structures for reduced THz absorption and more efficient outcoupling [3, 5, 8]. Still, regarding the power levels these sources are clearly outperformed by conventional THz QCLs, where cw output powers above 100mW have been demonstrated, albeit at cryogenic temperatures [11].

For most applications, THz powers of at least several mW at room temperature are desirable, preferably tunable over an extended frequency range. Regarding the active region, a design priority for further enhancing the output power is to achieve a large nonlinear susceptibility [5], which should not degrade with temperature and be broadband for wide tunability. Hence, for a further systematic design improvement, a detailed understanding of the artificial nonlinear susceptibility and its dependence on operating parameters is essential. In this context, self-consistent carrier transport simulations are especially valuable, since they do not rely on free fitting parameters and can provide insight into the optical nonlinearity on a microscopic level, hardly accessible to experiment [12, 13]. Recently, we have adapted the ensemble Monte Carlo (EMC) method [14], which is widely used for carrier transport simulations of MIR [15–19] as well as THz [20–27] QCLs, to the modeling of THz DFG structures [12]. Based on this approach, we investigate the dependence of the nonlinear susceptibility on the temperature, bias voltage and frequency detuning of the MIR pump beams for experimental THz DFG structures. The simulation yields good agreement with available experimental data, enables interpretation of the experimental results, and provides valuable insight for the further optimization of QCL-based THz DFG sources.

2. Theoretical model

In THz DFG structures, two MIR QCL active regions operating at frequencies ω₁ = 2π f₁ and ω₂ = 2π f₂ are commonly sandwiched into a single waveguide, generating THz radiation at a frequency ω₃ = ω₁ −ω₂. In our modeling approach [12], the carrier transport is simulated employing the EMC method, and the subband eigenenergies and wavefunctions are found using a Schrödinger-Poisson solver [14, 28]. The nonlinear frequency conversion process is in each active region described by the corresponding effective second order susceptibility. Employing a microscopic approach, the corresponding expression is obtained from a density matrix description of the light-matter interaction by applying time-dependent perturbation theory [29],

Assuming periodicity of the QCL active region, it is appropriate to consider the states ℓ of a single central period, along with all available states m,n (including those in neighboring periods). L_P then corresponds to the length of a single QCL period. The constants e, ε₀ and denote the elementary charge, vacuum permittivity and reduced Planck constant, respectively. The occupation probability f_ℓ (ε) of each subband ℓ is directly obtained from the EMC simulation as a function of the kinetic electron energy ε. Nonparabolicity effects are in our modeling approach considered by assigning individual effective masses $m_{\ell }^{*}$ to each subband, which are extracted from the Schrödinger-Poisson solver [14]. The resonance frequency of an optical transition between two subbands i = m,n with corresponding eigenenergies E_i is then given by ω_mn (ε) = (E_m +ε_m − E_n −ε_n)/ħ, where ${\epsilon }_i=\epsilon m_{\ell }^{*}/m_i^{*}$. The associated optical linewidth is self-consistently determined from the EMC simulation assuming lifetime broadening, i.e., γ_mn (ε) = [γ_m (ε_m) +γ_n (ε_n)]/2 with the intersubband outscattering rates γ_i (ε) [30]. Furthermore, −ez_mn is the transition dipole matrix element, calculated from the wavefunctions provided by the Schrödinger-Poisson solver.

The generated THz power generally scales with |χ⁽²⁾|². Specifically, for outcoupling through the front facet, the available THz output power is given by [12, 31]

P_i and n_i are the optical power and mode refractive index at frequency ω_i, c is the vacuum speed of light, and Δ_k = (ω₃n₃ − ω₁n₁ + ω₂n₂)/c is the phase mismatch. Furthermore, a_w and T_f are the THz waveguide loss coefficient and facet transmittance at frequency ω₃, and S_eff denotes the effective interaction area. For Cherenkov emission schemes, detailed modeling of the waveguide structure is needed to determine the outcoupled power [5, 32].

Our simulation approach does not contain adjustable fitting parameters and allows for a self-consistent evaluation of the carrier transport and nonlinear susceptibility, Eq. (1) [12]. The structure studied in Sections 3 and 4 contains a bound-to-continuum (BTC) and a double-resonant-phonon (DRP) active region, designed for f₁ = 32THz and f₂ = 28THz, respectively [1, 2]. For each active region, a separate carrier transport simulation is performed based on EMC, coupled to a Schrödinger-Poisson solver which provides the subband eigenenergies and wave-functions [14, 28]. Furthermore, the interaction with the optical cavity field is self-consistently included [33, 34] based on the experimentally specified waveguide parameters [2]. The simulation has to account for the fact that the current density is identical in both active regions. Furthermore, the optical intensity at f₁ generated in the BTC active region also affects the carrier dynamics in the DRP section, and the BTC region also interacts with the radiation at f₂. This is taken into account by performing alternate simulations of both active regions: For the DRP section, the intensity at f₁ is fixed to the value extracted from the simulation of the BTC section at the same current density, and vice versa. Both active regions are simulated iteratively until convergence is reached [12]. In Section 5, a widely tunable THz DFG structure is investigated [5, 9], which involves simulations at various frequencies f₁ to evaluate the spectral dependence of χ⁽²⁾. To reduce the numerical load of EMC, the MIR pump powers are here set to the experimentally measured values at the corresponding frequencies f₁ [5], thus avoiding the need to perform alternate simulations of both active regions.

Since the generated THz power is relatively weak, its influence on the carrier transport is negligible. Thus, only the MIR modes have to be considered in the EMC simulations, and the THz DFG process can be subsequently evaluated based on the obtained EMC results. To calculate the available THz power for facet outcoupling, Eq. (2), the quantities T_f, S_eff, and a_w of the fundamental THz mode are determined from waveguide simulations based on the effective index approximation [14, 35]. Here, the complex permittivities of the active regions are directly extracted from the EMC simulation [34], while the Drude model is used for the bulk layers [36].

3. Temperature dependence of difference frequency generation

We have used our simulation approach to analyze the nonlinear susceptibility in an experimental THz DFG structure with a distributed feedback grating, where the THz radiation is outcoupled through the front facet [2]. The MIR and THz powers are self-consistently obtained by coupled carrier transport and optical waveguide simulations, as described in [12] and in Section 2.

In Fig. 1, the physical quantities governing the THz output power P_THz, Eq. (2), are displayed as a function of the bias applied across the BTC structure for different temperatures. |χ⁽²⁾|, shown in Fig. 1(a), exhibits a moderate bias dependence, and the peak value degrades only slightly with increasing temperature. Also the waveguide loss coefficient a_w, displayed in Fig. 1(b), only changes slightly with bias and temperature. Thus, the temperature degradation of P_THz, shown in Fig. 1(d), is dominated by the decreasing MIR power product P₁P₂ displayed in Fig. 1(c), while the conversion efficiency η = P_THz/(P₁P₂) changes much less with temperature, see Fig. 1(e). Since the term ${\left({\mathrm{\Delta }}_k^2+a_{\mathrm{w}}^2/4\right)}^{-1}$ in Eq. (2) is governed by a_w [2], η is approximately proportional to the figure of merit ${\left|{\chi }^{(2)}\right|}^2{a}_{\mathrm{w}}^{-2}$ plotted in Fig. 1(f). By comparison with Fig. 1(e) it becomes clear that the bias and temperature dependence of η is largely determined by ${\left|{\chi }^{(2)}\right|}^2{a}_{\mathrm{w}}^{-2}$, while the other quantities such as effective interaction area and refractive indices do not change noticeably. For optimizing the conversion efficiency in an active region design it is thus essential to optimize the figure of merit ${\left|{\chi }^{(2)}\right|}^2{a}_{\mathrm{w}}^{-2}$, rather than only aiming at a maximum |χ⁽²⁾| [37]. In this context, it is important to consider that the THz waveguide loss a_w depends not only on the active region design, but also on the other waveguide layers including the substrate [1]. Notably, the absorption coefficients of both active regions decrease with increasing temperature and exhibit a more pronounced bias dependence than a_w. For an optimization of the THz output power, it is furthermore important to note that the figure of merit and the MIR power product usually reach their respective maxima at different biases, e.g., 38kV/cm and 40kV/cm at room temperature. In this case, the maximum THz output power is then obtained at 39kV/cm where |χ⁽²⁾| = 45nm/V and a_w = 150cm⁻¹, which agrees well with the experimentally estimated values of |χ⁽²⁾| ≈ 40nm/V and a_w ≈ 150cm⁻¹ [2].

 

Fig. 1 Bias and temperature dependence of the physical quantities governing the THz output power, Eq. (2): (a) Nonlinear susceptibility; (b) waveguide loss coefficient; (c) MIR power product; (d) resulting THz output power; (e) conversion efficiency; (f) figure of merit for difference frequency generation.

Download Full Size | PDF

4. Bias dependence of the optical nonlinearity and validity of simplified models

In the following, we want to take a closer look at the underlying reasons for the moderate bias dependence of |χ⁽²⁾|, as shown in Fig. 1(a). In this context, also the validity of simplified models for the nonlinear susceptibility is investigated. For a rough estimation of the nonlinear susceptibility in experimental setups, simplified expressions are commonly used [1, 5, 31]. In particular, the intrasubband kinetic electron distribution is not considered, and nonparabolicity effects are neglected, which corresponds to setting all subband effective masses in Eq. (1) to the value at the conduction band bottom, $m_{\ell }^{*}=m^{*}$. Under these assumptions, Eq. (1a) simplifies with $n_{\ell }^{\mathrm{s}}=m^{*}{\displaystyle {\int }_0^{\infty }{f}_{\ell }\mathrm{d}\epsilon }/\left(\pi {\hslash }^2\right)$ to

The bias dependence of |χ⁽²⁾| in Fig. 1(a) results from the change of the wavefunctions and eigenenergies, implicating altered resonance frequencies ω_mn and, more importantly, matrix elements z_mn in Eq. (1). In addition, the carrier transport is affected. The resulting changes in occupation probabilities f_ℓ also contribute significantly to the bias dependence of |χ⁽²⁾|, while linewidth changes play only a secondary role in this context. In Fig. 2(a), the approximate |χ⁽²⁾| computed from Eq. (3) is shown as a function of bias at room temperature. For comparison, the corresponding result from Eq. (1) is also displayed. While the obtained values agree reasonably well, the bias dependence is smoother for Eq. (1) which includes the intrasubband kinetic electron distributions and nonparabolicity effects. This is mainly due to the fact that with nonparabolicity included, the resonance frequencies ω_mn (ε) in Eq. (1) are not discrete, but depend on the kinetic electron energy ε. Furthermore, the contributions of the subband triplets 9-10-3’, 10-2’-3’ and 9-10-2’ to Eq. (3) are shown, where the prime denotes subbands belonging to the next-higher period. These are the main contributions at a bias of 39kV/cm where the maximum THz power is obtained, see Fig. 1(d). The bias dependence of the subband energies E_i, expressed in terms of frequency E_i/(2πħ), is displayed in Fig. 2(b), the conduction band profile and probability densities at 39kV/cm are shown in Fig. 2(c). The MIR transition frequencies are f_3′9 = 35.7THz, f_2′9 = 33.8THz, f_3′10 = 32.1THz, and f_2′10 = 30.3THz, which are not exactly in resonance with the MIR pump frequencies f₁ = 32THz and f₂ = 28THz. Also for other biases, χ⁽²⁾ is typically built up from different not perfectly resonant contributions, rather than being associated with a single resonant triplet. This explains the only moderate bias dependence of the nonlinear susceptibility as shown in Fig. 1(a), and is also consistent with the experimentally observed broad optical bandwidth of the nonlinearity [4, 5, 9].

 

Fig. 2 (a) Nonlinear susceptibility vs. applied bias, as obtained with Eq. (1) (solid) and Eq. (3) (dashed), and contributions of individual subband triplets to Eq. (3) (dotted). (b) Subband energies, expressed in terms of frequency, vs. bias. (c) Conduction band profile and probability densities at 39kV/cm.

Download Full Size | PDF

The carrier transport simulation allows us to determine the microscopic quantities relevant for the optical nonlinearity, such as subband occupations and transition linewidths which are difficult to obtain from experimental measurements. The linewidths in Eq. (3) are extracted from the EMC simulations assuming lifetime broadening, i.e., γ_mn = (γ_m +γ_n)/2 with the intersubband outscattering rates γ_i [30]. In Fig. 3(a), the relative occupation $p_i=n_i^s/n^s$, where n^s is the total sheet electron density per period, and level broadening γ_i/(2π) is shown for the subbands i of two adjacent QCL periods as a function of energy, expressed in terms of frequency. The level occupations of the closely spaced subbands in a period resemble a thermal distribution. Specifically, we obtain p₉ = 0.071, p₁₀ = 0.068, p_2′ = 0.132, and p_3′ = 0.095. The corresponding level broadenings are 2.02THz, 1.79THz, 2.05THz, and 1.40THz, resulting in linewidths of around 2THz or 8meV which is consistent with experimental estimates [1, 5]. Thus, the Lorentzian transition linewidths are comparable to the energy separation between the closely spaced subbands in the bound-to-continuum active region, implicating that several subbands contribute to the optical nonlinearity as demonstrated in Fig. 2(a).

 

Fig. 3 (a) Simulated relative level occupation (crosses) and level broadening (points) for two QCL periods at 39kV/cm. (b) and (c) Diagram illustrating near-resonant DFG between QCL subbands with (b) ω₁ ≈ ω₁₃, ω₂ ≈ ω₁₂, ω₃ ≈ ω₂₃, and (c) ω₁ ≈ ω₁₃, ω₂ ≈ ω₂₃, ω₃ ≈ ω₁₂.

Download Full Size | PDF

A rough estimate of the nonlinear susceptibility can be obtained by considering the conduction band profile at a specified design bias where the nonlinearity is assumed to be governed by a single subband triplet which is in near-resonance with the MIR pump beams [1, 5]. In Figs. 3(b) and 3(c), the two cases of near-resonant DFG are illustrated. In Fig. 3(b), ω₁ ≈ ω₁₃, ω₂ ≈ ω₁₂, and ω₃ ≈ ω₂₃ is assumed. Keeping only the corresponding quasi-resonant terms in Eq. (3), we obtain [1, 5]

Analogously, for the case illustrated in Fig. 3(c) with ω₁ ≈ ω₁₃, ω₂ ≈ ω₂₃, and ω₃ ≈ ω₁₂, Eq. (3) can be simplified to

The subband triplets 9-10-3’ and 10-2’-3’ yielding the largest contributions to χ⁽²⁾ at 39kV/cm, see Fig. 2(a), correspond to the configurations illustrated in Figs. 3(b) and 3(c), respectively. In this context, it should be pointed out that the triplets 9-10-3’ and 10-2’-3’ contribute to the nonlinearity constructively although Eq. (5) features an additional minus sign, since z_3′10z_10,9z_3′9 > 0 and z_3′2 z_2′10z_3′10 < 0. Eqs. (4) and (5) can however not directly be applied since the resonances are detuned away from the MIR pump frequencies, as discussed above. Assuming MIR pump beams in resonance with the triplet 9-10-3’ and taking the values for the transition linewidths and subband occupations from Fig. 3(a), Eq. (4) yields |χ⁽²⁾| = 35nm/V which gives a reasonable estimate for the total nonlinear susceptibility.

5. Frequency dependence of the nonlinear susceptibility

Recently, widely tunable THz DFG has been demonstrated [4, 5, 9], with a maximum tuning range of 1.2−5.9THz [9]. This approach takes advantage of the fact that a broad spectral range in the THz regime can be covered by tuning one of the MIR pumps over a few THz, corresponding only to a fraction (~10%) of the MIR center wavelength. The optical nonlinearity has to be sufficiently broadband in order to cover the complete tuning range. In this context, the sensitivity of the active region nonlinear susceptibility to a change in applied bias can be used as an indicator for its frequency dependence, since a bias change always involves a change in transition frequencies. Accordingly, the moderate bias dependence of |χ⁽²⁾| in Fig. 1(a) serves as an indication for a weak frequency dependence of the active region nonlinear susceptibility. Also the available experimental results suggest that the optical nonlinearity is sufficiently broadband to achieve wide frequency tunability in the THz regime [4, 5, 9, 38].

In the following, we theoretically investigate the frequency dependence of the nonlinear susceptibility, focusing on above mentioned broadly tunable DFG design which has achieved a record THz tuning range [5, 9]. The two active regions, designed for center wavelengths of 8.2 μm and 9.2 μm, respectively, both exhibit a giant optical nonlinearity to increase the conversion efficiency. In the investigated experimental setup, one wavelength is fixed to 10.3 μm, while the second MIR pump is tuned between 8.6 and 9.8 μm [5]. The simulation is again based on EMC including carrier-light coupling. To reduce the numerical load, the MIR pump powers are here set to the experimentally measured values at the corresponding wavelengths [5], as discussed in Section 2. In Fig. 4, the simulated nonlinear susceptibilities and THz absorption coefficients are shown as a function of applied bias and frequency detuning of the MIR pump beams for both the 8.2 μm and 9.2 μm active region design. The peak values for |χ⁽²⁾| in Figs. 4(a) and 4(b) are 34nm/V and 75nm/V, respectively, which is comparable to the results for the DFG structure discussed above, see Fig. 1(a). The nonlinear susceptibility changes only slightly over an extended bias and spectral range. Specifically, for suitably chosen biases, |χ⁽²⁾| is close to its maximum value over the complete relevant frequency range. The underlying reasons have already been discussed in Section 4: The closely spaced subbands in the bound-to-continuum active region, with transition linewidths comparable to their energy separation, provide a broadband optical nonlinearity, where several subbands contribute to |χ⁽²⁾|. By comparison, the spectral and bias dependence of the absorption coefficients, shown in Figs. 4(c) and 4(d), is much more pronounced. In this context, it is important to consider that the total THz loss also depends on the other waveguide layers as already discussed in Section 3, and especially on the chosen waveguide geometry. In particular, the THz absorption can be significantly reduced by employing a Cherenkov emission scheme as used for the widely tunable THz source investigated in this section [5, 9], rather than a facet outcoupling geometry.

 

Fig. 4 Nonlinear susceptibility for the (a) 8.2μm and (b) 9.2μm design and THz absorption coefficent for the (c) 8.2μm and (d) 9.2μm design as a function of the applied bias and frequency detuning.

Download Full Size | PDF

6. Conclusion

Based on EMC simulations including carrier-light coupling, we have investigated the giant nonlinear susceptibility in QCL-based THz DFG structures. The optical nonlinearity exhibits a low temperature sensitivity, thus the temperature degradation of the THz output power in these structures is largely due to the decrease of MIR pump powers. Furthermore, the nonlinear susceptibility is close to its maximum value over an extended bias and frequency range, enabling frequency tuning over a wide THz spectral range. The underlying reason is that the bound-to-continuum active region providing the optical nonlinearity features closely spaced subbands with transition linewidths comparable to their energy separation. Thus, typically several levels contribute to |χ⁽²⁾| rather than a single quasi-resonant subband triplet, reducing the bias and frequency dependence of the susceptibility and providing a broadband optical nonlinearity.