## Abstract

We present an extended ensemble Monte Carlo approach, allowing for the self-consistent modeling of terahertz difference frequency generation in quantum cascade lasers. Our simulations are validated against available experimental data for a current room temperature design. Tera-hertz output powers in the mW range are predicted for ideal light extraction.

© 2013 Optical Society of America

## 1. Introduction

Quantum cascade laser (QCL) sources based on difference frequency generation (DFG) constitute a very promising approach for the generation of coherent terahertz (THz) radiation [1]. The DFG process is here implemented by dual-wavelength mid-infrared (MIR) QCLs pumping a giant optical nonlinearity, which is commonly directly integrated into the QCL gain medium [2–5]. With such intracavity DFG sources, room temperature operation has been demonstrated [2, 4, 5], whereas in conventional THz QCLs operating temperatures are still restricted to below 200 K [6]. On the other hand, while output powers of several 100mW have been obtained with THz QCLs [7], DFG structures still suffer from low THz powers in the *μ*W range [2, 4, 5]. For a further improvement of these sources, a detailed understanding of the MIR generation, nonlinear frequency conversion and THz outcoupling is essential. In this context, self-consistent state of the art carrier transport simulations which do not rely on free fitting parameters or require experimental input are especially valuable.

Here we employ the ensemble Monte Carlo (EMC) method, which has proven very helpful for the design, analysis and optimization of MIR [8–10] and THz [11–16] QCLs. For the investigation of DFG-based THz sources, an inclusion of the optical dynamics is crucial, which is typically neglected in EMC and other advanced carrier transport simulation methods. Recently we have developed an extended EMC approach, accounting for the carrier transport and the optical cavity field on an equal footing [10, 17, 18]. Here we adapt our tool to the simulation of DFG-based THz sources, self-consistently including the nonlinear frequency conversion process. We employ our approach to analyze a recent DFG structure which achieved record THz output powers at room temperature [4]. The simulation yields good agreement with available experimental data. Specifically, our approach allows us to identify the main factors limiting the emitted THz power, and to explore the theoretical limit for optimum outcoupling. From our calculations, room temperature output powers in the mW range are found to be accomplishable.

## 2. Modeling of difference frequency generation

As mentioned above, in DFG structures two MIR QCLs at frequencies *ω*_{1} and *ω*_{2} are used to pump a medium with giant optical nonlinearity, generating THz radiation at a frequency *ω* = *ω*_{1} − *ω*_{2}. The effective second order susceptibility in the nonlinear medium is given by [19]

*m*

^{*},

*ε*

_{0}and

*ħ*denote the effective mass, vacuum permittivity and reduced Planck constant, respectively. Furthermore,

*ω*= (

_{mn}*E*−

_{m}*E*) /

_{n}*ħ*and

*d*are the resonance frequency and the dipole matrix element of the optical transition between the subbands

_{mn}*i*=

*m*,

*n*with corresponding eigenenergies

*E*. The optical linewidth

_{i}*γ*(

_{mn}*ε*) and the occupation probability

*f*(

_{i}*ε*) depend on the kinetic electron energy

*ε*.

To evaluate Eq. (1), we can proceed in analogy to EMC-based optical gain calculations [20]: The occupation probability is directly extracted from the simulation, and the linewidth is evaluated based on lifetime broadening, yielding *γ _{mn}* (

*ε*) = [

*γ*(

_{m}*ε*) +

*γ*(

_{n}*ε*)]/2 with the intersubband outscattering rates

*γ*(

_{i}*ε*) [20]. For periodic structures such as QCLs, it is sufficient to consider the states

*ℓ*of a single central period, along with all available states

*m*,

*n*(including those in neighboring periods);

*L*

_{P}is then the length of a QCL period. Nonparabolicity effects are included in our Schrödinger-Poisson solver [21], and are also considered in the EMC simulation by assigning individual effective masses ${m}_{\ell}^{*}$ to each subband

*ℓ*[10]. We implement the different ${m}_{\ell}^{*}$ into Eq. (1) by taking ${n}_{E,\ell}^{2D}={m}_{\ell}^{*}/\left(\pi {\hslash}^{2}\right)$, and using modified definitions

*ω*(

_{mn}*ε*) = (

*E*+

_{m}*ε*−

_{m}*E*−

_{n}*ε*) /

_{n}*ħ*,

*γ*(

_{mn}*ε*) = [

*γ*(

_{m}*ε*) +

_{m}*γ*(

_{n}*ε*)] / 2 with ${\epsilon}_{i}=\epsilon {m}_{\ell}^{*}/{m}_{i}^{*}$.

_{n}For the optical propagation in the resonator along the *z* direction, we write the electric field as *E _{i}* (

*x*,

*y*,

*z*) =

*A*(

_{i}*z*)

*F*(

_{i}*x*,

*y*)exp(i

*k*− i

_{i}z*ω*) + c.c. [22], with (∬|

_{i}t*F*|

_{i}^{2}d

*x*d

*y*)

^{1/2}= 1.

*A*and

_{i}*F*denote the complex amplitude and transverse mode distribution of the field component with frequency

_{i}*ω*and wavenumber

_{i}*k*, and c.c. refers to the complex conjugate. The optical power at frequency

_{i}*ω*is then given by

_{i}*P*= 2

_{i}*cε*

_{0}

*n*|

_{i}*A*|

_{i}^{2}, with the effective refractive index of the corresponding mode

*n*and the vacuum speed of light

_{i}*c*. The propagation of the THz field component generated by the DFG process is in the slowly varying amplitude approximation described by [23]

*=*

_{k}*k*

_{3}−

*k*

_{1}+

*k*

_{2}is the phase mismatch,

*a*denotes the THz resonator loss coefficient at frequency

*ω*, and

*χ*

^{(2)}is the peak value of the nonlinear susceptibility

*χ*

^{(2)}(

*x*,

*y*) in the QCL waveguide. Furthermore, the asterisk denotes the complex conjugate. The overlap integral

*f*

_{321}is defined as ${f}_{321}={\left({\chi}^{\left(2\right)}\right)}^{-1}\iint {\chi}^{\left(2\right)}\left(x,y\right){F}_{3}^{*}{F}_{2}^{*}{F}_{1}\text{d}x\text{d}y$[1, 22].

For solving Eq. (2), *P*_{1} and *P*_{2} are commonly taken to be constant over the resonator length *L*, which holds for moderate MIR outcoupling at the facets. Furthermore, the depletion of the MIR pump beams due to the DFG process is negligible. The calculation is greatly simplified for *L* ≫ *a*^{−1}. With *a* ∼ 100cm^{−1} and *L* in the mm range [1,2,4], this condition is typically fulfilled for DFG structures. Introducing the effective interaction area *S*_{eff} = |*f*_{321}|^{−2} and the coherence length
${L}_{\text{coh}}={\left({\mathrm{\Delta}}_{k}^{2}+{a}^{2}/4\right)}^{-1/2}$, the THz power at the facet *P*_{f} = *P*_{3} (*L*) is then obtained as [1]

*∂*

_{z}P_{3}|

_{a}= −

*aP*

_{3}in the resulting equation corresponds to the THz power lost in the waveguide. Integration over the propagation length yields the total absorbed power

*P*

_{a}. For

*L*≫

*a*

^{−1}, we can approximate

*P*

_{3}=

*P*

_{f}and thus obtain

*P*

_{a}= −2

*LaP*

_{f}, where 2

*L*is the cavity roundtrip length and

*P*

_{f}is given by Eq. (3). The total generated power is then

The large losses due to the strong THz absorption in the waveguide can be reduced by out-coupling the generated THz radiation along the resonator surface rather than at the facet. This can for example be achieved by a grating etched into the laser ridges [3], or in a Cherenkov emission scheme [5]. For surface-emitting structures characterized by a THz outcoupling coefficient *κ*, the theoretical description can easily be adapted by replacing *a* with *a* + *κ* in Eq. (2). In analogy to Eq. (4), we obtain for *L* ≫ (*a* + *κ*)^{−1} the generated power

*P*

_{f}is given by Eq. (3) with a modified coherence length ${L}_{\text{coh}}={\left[{\mathrm{\Delta}}_{k}^{2}+{\left(a+\kappa \right)}^{2}/4\right]}^{-1/2}$. In Eq. (5), the first contribution 2

*LaP*

_{f}corresponds to the power absorbed in the waveguide, and the second contribution is the power outcoupled at the surface

## 3. Analysis of a terahertz source based on difference frequency generation

In the following, we present a self-consistent analysis of a room temperature DFG source, producing a THz output power of 8.5*μ*W in single mode operation [4]. The active region consists of a bound-to-continuum (BTC) section which also contains the giant optical nonlinearity for the DFG process (see Fig. 1(a)), and a double-resonant-phonon (DRP) section. A distributed feedback grating is integrated into the resonator to enforce single mode lasing of the BTC and DRP design at frequencies of 32 and 28THz, respectively, resulting in DFG at 4THz [4].

#### 3.1. Device simulation

The carrier transport in the QCL structure is modeled based on a semiclassical EMC approach, specifically adapted to the simulation of MIR QCLs [10]. All relevant scattering mechanisms are considered [10,20,24], and the dynamics due to the optical cavity field is taken into account [10, 17, 18]. The carrier transport simulation is coupled to a Schrödinger-Poisson solver [21], providing the subband energies and wave functions. The simulations are self-consistent, only relying on the device specifications and well-established material parameters.

Based on the given layer sequences of the BTC and DRP design [2,4], the EMC simulations provide the MIR powers and the electric current. The waveguide length is *L* = 3mm and the width is *w*_{0} = 16*μ*m, tapered to *w*_{f} = 60*μ*m toward the front facet with a small taper angle of 1° [4], so that the mode can expand adiabatically (see Fig. 1(b)). This results in a contact area of *A* = 0.076mm^{2} and an average waveguide width of *w*̄ = *A/L* = 25*μ*m, used to calculate the current and MIR powers from the current density and optical intensities delivered by EMC. Furthermore, we assume a confinement factor of 0.4 in both gain sections [2,4], and a combined waveguide and mirror loss of 6.4cm^{−1} + 2.1cm^{−1} for both MIR modes. Here the common assumption is made that the optical power in each mode is equally distributed over the resonator length [17]. The 32THz mode generated in the BTC section influences the carrier dynamics in the DRP gain medium, and analogously for the 28THz mode. This is considered by performing alternate simulations of the BTC and DRP sections: For the DRP region, the intensity of the 32THz mode is fixed to the value extracted from the BTC region simulation at the same current density, and vice versa. The simulations are repeated until convergence is reached.

The DFG process is simulated using the converged MIR EMC simulation results, since the generated THz power is so weak that its influence on the carrier transport can be neglected. The second order susceptibility is calculated from Eq. (1). Based on waveguide simulations using the effective refractive index approximation [25], *S*_{eff}, *a*, and the facet transmittance *T*_{f}[26] are calculated for the fundamental THz mode. Here, the Drude model is employed for the bulk layers [27], while the complex permittivity of the gain medium is directly extracted from the EMC simulation [18]. With THz losses of *a* ≈ 150cm^{−1}[4], only radiation generated at ≲ 0.1mm from the facet contributes to the THz emission. Thus, *S*_{eff}, *a* and *T*_{f} are calculated for the waveguide width at the facet *w*_{f} = 60*μ*m. Furthermore, a phase mismatch of Δ* _{k}* = 7cm

^{−1}≪

*a*is assumed [4], corresponding to near-perfect phase matching. The THz power outcoupled through the facet can then be calculated using

*P*

_{out}=

*T*

_{f}

*P*

_{f}, with

*P*

_{f}given by Eq. (3).

#### 3.2. Results

To validate our modeling approach, we compare our simulation results at 300K to available experimental data for a recent room temperature DFG source [4]. For an optimum bias where the THz output power *P*_{out} reaches its maximum value, the EMC simulation yields a current of 11.7A, which agrees very well with the measured value of 12A [4]. Also the simulated MIR powers *P*_{1} = 1.55W (32THz) and *P*_{2} = 0.60W (28THz) agree reasonably well with the measured values of 1.2 and 0.8W [4]. The THz performance of the DFG source is characterized by *P*_{out} and the conversion efficiency *η* = *P*_{out}/ (*P*_{1}*P*_{2}) [4]. With *S*_{eff} = 5300*μ*m^{2}, *T*_{f} = 0.23, |*χ*^{(2)}| = 44nm/V and *a* = 150cm^{−1}, we obtain *P*_{out} = 12*μ*W and *η* = 13 *μ*W/W^{2}, as compared to 8.5*μ*W and 10*μ*W/W^{2} found in experiment [4]. Furthermore, the simulation confirms the experimentally estimated values for the nonlinear susceptibility (|*χ*^{(2)}|≈ 40nm/V) and the waveguide loss (*a* ≈ 150cm^{−1}) [4]. Thus, good agreement with experiment is obtained. Remaining deviations are attributed to uncertainties in the MIR waveguide and interface roughness parameter values [10], as well as the implementation of the MIR mirror outcoupling by a distributed loss coefficient as mentioned above. In this context, we note that the experimental MIR and THz powers seem to strongly depend on the exact design parameters [2, 4].

Figure 2 contains temperature dependent simulation results for the investigated DFG source. In Fig. 2(a), the MIR powers *P*_{1} and *P*_{2} are displayed for optimum bias where *P*_{out} is maximum. The corresponding outcoupled THz power *P*_{out} and conversion efficiency *η* are shown in Fig. 2(b). In agreement with experiment [2], *P*_{1}, *P*_{2} and *P*_{out} strongly degrade with increasing temperature, while *η* is insensitive over a wide temperature range.

As shown in Section 2, the generated THz power can be outcoupled more efficiently by using a surface outcoupling scheme, furthermore offering improved beam collimation [3]. The corresponding surface-emitted power *P*_{s} in Eq. (6) is maximized for near-perfect phase matching Δ* _{k}* ≪

*a*and a coupling coefficient

*κ*=

*a*. In the following, we assume ideal surface emission from our investigated structure in order to estimate the maximum extractable THz power. The corresponding surface grating is here considered to have a negligible effect on the MIR laser modes [3]. Thus, the EMC simulation results remain unaffected. For better comparability, we furthermore assume that also the THz mode profile and loss are the same as for the facet emitting device. Since the outcoupling now occurs across the whole surface area

*A*rather than at the facet, the average waveguide width

*w*̄ =

*A/L*instead of

*w*

_{f}is used for the waveguide simulations, yielding a reduced

*S*

_{eff}= 2200

*μ*m

^{2}while

*a*does not change significantly. In Fig. 2(c), the emitted THz power

*P*

_{s}and conversion efficiency

*η*

_{s}=

*P*

_{s}/ (

*P*

_{1}

*P*

_{2}) are shown. The outcoupled power at 300K is now

*P*

_{s}= 3.1mW, suggesting that room temperature output powers in the mW range should be feasible with an improved outcoupling scheme. From Eqs. (5) and (6) it follows that for

*κ*=

*a*, 50% of the generated THz power is emitted. By contrast, for the case of facet emission shown in Fig. 2(b), the outcoupling ratio is only

*P*

_{out}/

*P*

_{g}≈ 0.1%. On a practical note, experimental structures with grating-based outcoupling have up to now only yielded moderate improvement [3], partly because of the difficulties in implementing a strong grating with

*κ*≈

*a*[5]. An alternative surface emission scheme based on the Cherenkov effect has recently shown promising results [5].

## 4. Conclusion

An extended EMC approach for the self-consistent modeling of terahertz difference frequency generation in QCLs is developed and applied to a current experimental room temperature design. Good agreement with available experimental data is obtained. For optimized light extraction via surface outcoupling, room temperature output powers in the mW range are predicted.

## Acknowledgments

This work was funded by the Nanosystems Initiative Munich and the Emmy Noether program of the DFG ( JI115/1-1). We thank M. Belkin for fruitful discussions.

## References and links

**1. **M. Belkin, F. Capasso, A. Belyanin, D. Sivco, A. Cho, D. Oakley, C. Vineis, and G. Turner, “Terahertz quantum-cascade-laser source based on intracavity difference-frequency generation,” Nat. Photonics **1**, 288–292 (2007) [CrossRef] .

**2. **M. Belkin, F. Capasso, F. Xie, A. Belyanin, M. Fischer, A. Wittmann, and J. Faist, “Room temperature terahertz quantum cascade laser source based on intracavity difference-frequency generation,” Appl. Phys. Lett. **92**, 201101 (2008) [CrossRef] .

**3. **C. Pflügl, M. Belkin, Q. Wang, M. Geiser, A. Belyanin, M. Fischer, A. Wittmann, J. Faist, and F. Capasso, “Surface-emitting terahertz quantum cascade laser source based on intracavity difference-frequency generation,” Appl. Phys. Lett. **93**, 161110 (2008) [CrossRef] .

**4. **Q. Y. Lu, N. Bandyopadhyay, S. Slivken, Y. Bai, and M. Razeghi, “Room temperature single-mode terahertz sources based on intracavity difference-frequency generation in quantum cascade lasers,” Appl. Phys. Lett. **99**, 131106 (2011) [CrossRef] .

**5. **K. Vijayraghavan, R. W. Adams, A. Vizbaras, M. Jang, C. Grasse, G. Boehm, M. C. Amann, and M. A. Belkin, “Terahertz sources based on Čerenkov difference-frequency generation in quantum cascade lasers,” Appl. Phys. Lett. **100**, 251104 (2012) [CrossRef] .

**6. **S. Fathololoumi, E. Dupont, C. Chan, Z. Wasilewski, S. Laframboise, D. Ban, A. Mátyás, C. Jirauschek, Q. Hu, and H. Liu, “Terahertz quantum cascade lasers operating up to ∼200 K with optimized oscillator strength and improved injection tunneling,” Opt. Express **20**, 3866–3876 (2012) [CrossRef] [PubMed] .

**7. **B. S. Williams, S. Kumar, Q. Hu, and J. L. Reno, “High-power terahertz quantum-cascade lasers,” Electron. Lett. **42**, 89–90 (2006) [CrossRef] .

**8. **R. C. Iotti and F. Rossi, “Carrier thermalization versus phonon-assisted relaxation in quantum-cascade lasers: A Monte Carlo approach,” Appl. Phys. Lett. **78**, 2902–2904 (2001) [CrossRef] .

**9. **X. Gao, D. Botez, and I. Knezevic, “X-valley leakage in GaAs-based midinfrared quantum cascade lasers: A Monte Carlo study,” J. Appl. Phys. **101**, 063101 (2007) [CrossRef] .

**10. **A. Mátyás, P. Lugli, and C. Jirauschek, “Photon-induced carrier transport in high efficiency midinfrared quantum cascade lasers,” J. Appl. Phys. **110**, 013108 (2011) [CrossRef] .

**11. **R. Köhler, R. C. Iotti, A. Tredicucci, and F. Rossi, “Design and simulation of terahertz quantum cascade lasers,” Appl. Phys. Lett. **79**, 3920–3922 (2001) [CrossRef] .

**12. **H. Callebaut, S. Kumar, B. S. Williams, Q. Hu, and J. L. Reno, “Analysis of transport properties of terahertz quantum cascade lasers,” Appl. Phys. Lett. **83**, 207–209 (2003) [CrossRef] .

**13. **O. Bonno, J.-L. Thobel, and F. Dessenne, “Modeling of electron-electron scattering in Monte Carlo simulation of quantum cascade lasers,” J. Appl. Phys. **97**, 043702 (2005) [CrossRef] .

**14. **J. T. Lü and J. C. Cao, “Coulomb scattering in the Monte Carlo simulation of terahertz quantum-cascade lasers,” Appl. Phys. Lett. **89**, 211115 (2006) [CrossRef] .

**15. **C. Jirauschek, G. Scarpa, P. Lugli, M. S. Vitiello, and G. Scamarcio, “Comparative analysis of resonant phonon THz quantum cascade lasers,” J. Appl. Phys. **101**, 086109 (2007) [CrossRef] .

**16. **A. Mátyás, M. Belkin, P. Lugli, and C. Jirauschek, “Temperature performance analysis of terahertz quantum cascade lasers: Vertical versus diagonal designs,” Appl. Phys. Lett. **96**, 201110 (2010).

**17. **C. Jirauschek, “Monte Carlo study of carrier-light coupling in terahertz quantum cascade lasers,” Appl. Phys. Lett. **96**, 011103 (2010) [CrossRef] .

**18. **C. Jirauschek, “Monte Carlo study of intrinsic linewidths in terahertz quantum cascade lasers,” Opt. Express **18**, 25922–25927 (2010) [CrossRef] [PubMed] .

**19. **Y. Shen, *The Principles of Nonlinear Optics* (Wiley-Interscience, 1984).

**20. **C. Jirauschek and P. Lugli, “Monte-Carlo-based spectral gain analysis for terahertz quantum cascade lasers,” J. Appl. Phys. **105**, 123102 (2009) [CrossRef] .

**21. **C. Jirauschek, “Accuracy of transfer matrix approaches for solving the effective mass Schrödinger equation,” IEEE J. Quantum Electron. **45**, 1059–1067 (2009) [CrossRef] .

**22. **G. Agrawal, *Nonlinear Fiber Optics* (Academic, 2001).

**23. **N. Bloembergen, *Nonlinear Optics* (World Scientific, 1996).

**24. **C. Jirauschek, A. Matyas, and P. Lugli, “Modeling bound-to-continuum terahertz quantum cascade lasers: The role of Coulomb interactions,” J. Appl. Phys. **107**, 013104 (2010) [CrossRef] .

**25. **K. Chiang, “Performance of the effective-index method for the analysis of dielectric waveguides,” Opt. Lett. **16**, 714–716 (1991) [CrossRef] [PubMed] .

**26. **J. Butler and J. Zoroofchi, “Radiation fields of GaAs-(AlGa)As injection lasers,” IEEE J. Quantum Electron. **10**, 809–815 (1974) [CrossRef] .

**27. **S. Kohen, B. S. Williams, and Q. Hu, “Electromagnetic modeling of terahertz quantum cascade laser waveguides and resonators,” J. Appl. Phys. **97**, 053106 (2005) [CrossRef] .