We present a method to perform accurate and efficient simulations of photovoltaic quantum cascade detector (QCD) structures based on the ensemble Monte Carlo (EMC) approach. Since the photocurrent is typically orders of magnitude smaller than the pump current in a quantum cascade laser (QCL), a direct application of QCL simulation techniques is not sufficiently accurate. We demonstrate that by exploiting thermodynamic equilibrium relations for the electron populations and scattering rates, the EMC method can be adapted to yield reliable results for the essential QCD figures of merit, such as responsivity and specific detectivity. The modeling approach is validated against available experimental results for various mid-infrared and terahertz QCD designs and furthermore compared to simulations based on the non-equilibrium Green’s function method.

For the detection of light in the mid- and far-infrared regime, intersubband detectors, which use optical transitions between quantized levels in the conduction or valence band of a multi-quantum well system, are especially suited. They were first realized in the form of photoconductive quantum well infrared detectors (QWIPs).1 These are biased periodic quantum well (QW) structures where the absorbed photon causes an intersubband transition in the corresponding well, which results in a photocurrent. By breaking the inversion symmetry of the band, QWIPs can also be operated at zero bias based on the photovoltaic effect, with potential advantages, e.g., in terms of noise characteristics.2,3 Another strategy for the development of photovoltaic intersubband detectors is to start from a quantum cascade laser (QCL) design,4 which by default features an asymmetric band structure,5 and optimize it for photodetection.6,7 The resulting device is commonly referred to as a quantum cascade detector (QCD). Besides the above mentioned superior noise behavior associated with photovoltaic operation, especially at elevated temperatures, the QCD has several advantages over standard QWIPs, such as the increased design flexibility due to the more complex band structure and the compatibility with standard QCL fabrication technology. Furthermore, established QCL modeling techniques can be adapted to QCDs,8,9 while in standard QWIPs, the carrier transition between quantum-confined and continuum states is a difficult theoretical problem usually treated by adjustable parameters,10 which is unsuitable for systematic design optimizations. In Fig. 1(a), the conduction band profile and subbands of a basic exemplary QCD structure are schematically shown, where the arrows indicate possible scattering transitions between the subbands. Photo-excitation occurs from the ground level g to an absorption level a. The further levels form a staircase of energy states so that the photoexcited electron cascades down to the ground state of the next period through scattering. Figure 1(b) contains an equivalent circuit diagram for noise modeling, where each scattering path is represented by a noise equivalent resistance.

FIG. 1.

(a) Schematic illustration of a QCD in photovoltaic operation mode. The rectangle denotes a single QCD period. (b) Equivalent circuit diagram for noise modeling.

FIG. 1.

(a) Schematic illustration of a QCD in photovoltaic operation mode. The rectangle denotes a single QCD period. (b) Equivalent circuit diagram for noise modeling.

Close modal

For a systematic investigation and optimization of QCD designs, an advanced carrier transport simulation is required, which does not rely on empirical or fitting parameters as an input. A main challenge is that the currents in photovoltaic QCDs typically are significantly smaller than in QCLs, preventing a direct simulation of the current for the irradiated design with the standard QCL simulation approaches, such as the ensemble Monte Carlo (EMC) or the non-equilibrium Green’s function (NEGF) method. In addition, a direct application of these modeling approaches does not yield information about the noise current, which is an essential quantity for photodetectors. In contrast to the ample literature existing in the context of QCL modeling,11 not much work has been dedicated to advanced QCD simulations.8,9 Based on the EMC method, we discuss how to adapt and extend the QCL carrier transport simulation framework to obtain accurate results for the detector signal and noise. In detail, by exploiting thermodynamic equilibrium relations, we derive expressions for the responsivity and specific detectivity, which require only the absorption coefficient and intersubband scattering rates as input from the carrier transport simulations. In particular, the direct simulation of current densities, which tend to be very low in photovoltaic QCD operation and are thus prone to numerical inaccuracies, is here circumvented. Importantly, we avoid design specific assumptions and approximations in our responsivity and specific detectivity models to maximize the accuracy and versatility of our approach. Notably, this does not significantly increase the associated computational burden, given that the main computational cost arises from the carrier transport simulations. Furthermore, specific aspects are discussed with regard to EMC-based modeling of QCDs, and also NEGF simulations are presented.

This paper is organized as follows: Sec. II introduces the responsivity and specific detectivity as the main figures of merit for characterizing QCD performance and discusses their computation based on the simulated scattering rates. In Sec. III, an EMC-based modeling approach for photovoltaic QCD operation is introduced, taking advantage of the transport and statistical properties in thermodynamic equilibrium. In Sec. IV, the model is validated against experimental and NEGF results for both mid-infrared and terahertz QCD designs. The paper is concluded in Sec. V.

A frequently used figure of merit for characterizing photodetector performance is the responsivity Rp=Is/Ps, i.e., the generated photocurrent Is divided by the incident light power Ps. Although design and optimization of QCD responsivity could in principle be carried out assuming large incident powers, the typical photocurrent of QCDs under real operating conditions is too small to be extracted from direct simulations of the irradiated structure based on, e.g., EMC or NEGF. In addition, important properties such as the detectivity and dark current, requiring knowledge of the current vs voltage under dark conditions, are inaccessible for these models. This problem can be circumvented by rewriting the responsivity in terms of the structure’s power absorption efficiency ηabs and extraction efficiency pe. For a structure with Np periods, Rp is then given by8,12

(1)

Here, e, h=2π, and f have the usual meaning of elementary charge, Planck’s constant, and frequency, respectively. The optical intersubband transition has a non-vanishing dipole matrix element only in the growth direction of the structure [see Fig. 1(a)], and thus only interacts with the corresponding electric field component. For QCDs, the facet is often chosen to be in in-plane direction; thus, the detector has to be tilted to obtain an optical response. For this standard mesa geometry, we obtain

(2)

for p-polarized light, if interference effects are neglected.8,12–14 Here, θ is the propagation direction of the light inside the structure relative to the growth direction, and nP denotes the number of passes of the optical field through the active region. In particular, nP=1 for a single-pass geometry, and nP=2 for a double-pass geometry, as commonly realized by exploiting internal reflection at the top contact. Higher values of nP can, for example, be obtained by using a multipass waveguide geometry.12,15Tf is the facet transmittance and can potentially also include losses due to the substrate and contact layers. In order to optimize Tf, Fresnel reflections are often minimized by polishing the entrance facet at an angle of θ.12,15 Furthermore, α denotes the power absorption coefficient, and Lp is the length of a period so that NpLp corresponds to the total active region thickness. In Eq. (1), the only input required from the EMC simulation is the absorption coefficient α, which can be straightforwardly extracted,16,17 and pe, which can be calculated from the intersubband scattering rates, as discussed in Sec. II A.

For many applications involving photodetectors, the signal to noise ratio rather than the absolute signal strength is relevant. A suitable parameter is the specific detectivity D=Rp(AdΔf)1/2/IN, which gives the signal (per unit incident power) to noise ratio, normalized to the measurement electrical bandwidth Δf and the detector area Ad. Here, IN is the root mean square noise current. For photovoltaic (i.e., unbiased) detectors where Johnson noise dominates at elevated temperatures, we obtain IN=(4kBTΔf/Rd)1/2 where Rd is the detector resistance, T denotes the temperature and kB is the Boltzmann constant. Thus, we can write

(3)

Since Rd can be calculated from the intersubband scattering rates as discussed in Sec. II B, only these and the calculated absorption spectrum in thermodynamic equilibrium, i.e., for zero bias and no incident light, are required from EMC for our detector simulations, which greatly facilitates the modeling.

In the following, we discuss how to calculate the pe and Rd directly from the intersubband scattering rates delivered by EMC.

The extraction efficiency pe gives the probability that an electron which is photo-excited from the ground level g to an absorption level a in Fig. 1(a) is transported by scattering to the next period, rather than recombining to the ground state. In the pioneering work of Harrer et al.,8pe is determined by monitoring the scattering-driven evolution of the electron distribution in a closed three-period structure if initially an electron is placed in the excited state, in order to determine if it relaxes to the ground state of the same or an adjacent period. Here, the evolution has to be stopped at a carefully chosen point in time since the electron distribution otherwise inevitably evolves to thermodynamic equilibrium, which does not contain information about pe. For this reason, we have developed an alternative method, providing a robust and compact approach to determine pe from the EMC scattering rates. We assume that the QCD structure is periodic with N subbands per period, i.e., the transition rates from subband i to j fulfill rij=ri+N,j+N. Formulating the rate equation for a reference period and taking into account transitions to and from the left- and the right-neighboring period, we obtain in the stationary case11,18

(4)

Here, the transition rate matrix Q contains the elements qii=ji(ri,j+N+rij+ri+N,j) for i=1,,(N1), qNj=1 for j=1,,N, and qij=rj+N,i+rji+rj,i+N otherwise. Furthermore, vector p contains the subband occupation probabilities pi, and the vector elements of b are given by bi=δiN, where δiN denotes the Kronecker delta. For incident light, we obtain a change in rate δrga=rp from the ground level g to the absorption level a due to photon-induced transitions, yielding a modified matrix Q+δQ where the non-zero elements of δQ are given by δqag=rp, δqgg=rp. The associated stimulated emission process is considered by setting δqga=Γgarp, δqaa=Γgarp with Γij=gi/gj, where gi denotes the degeneracy of a level i. For intersubband transitions, Γij=1 holds in the case of parallel subbands, as discussed in detail in Sec. III A. The change in stationary level occupations δp resulting from stimulated photon transitions is found by evaluating (Q+δQ)(p+δp)=b. Re-inserting Eq. (4) and neglecting the higher order term δQδp, which is adequate for moderate intensities, i.e., if the detector operates in the linear response regime, we obtain

(5)

which yields the population changes δp in the irradiated structure. Here δQp=(Γgapapg)rpv, with vg=1, va=1, and vi=0 otherwise. The occupation changes can also be written in an explicit form by applying Cramer’s rule to Eq. (5), which yields

(6)

where Qi is the matrix obtained by replacing the ith column of Q by vector v, and || denotes the matrix determinant. The extraction efficiency is then given by the ratio of current change through the period to change in current from ground to upper level,

(7a)
(7b)

where Eq. (7b) has been obtained by inserting Eq. (6). As can be seen from Eq. (7b), pe is independent of the chosen value for rp, as is to be expected in the linear response regime assumed above. Thus, when evaluating Eq. (5) with a linear equation system solver rather than using Cramer’s rule, rp can be set to a numerically convenient value. We note that for photovoltaic QCDs, the stimulated emission contribution is negligible, i.e., pa0 in Eqs. (6) and (7b), if the energetic spacing Eag=EaEg of the optical transition fulfills EagkBT. This is fulfilled for mid-infrared detectors also at elevated temperatures, while terahertz QCD functionality has anyway only been observed under cryogenic conditions.19 

Often a QCD has two or more relevant upper absorption levels a.15 Then, in the linear response regime assumed above, the occupation change δpi in Eq. (7a) can be obtained by summing over the individual contributions of the absorption transitions. Here, the proportionate photon absorption rate carp must be used in Eqs. (5) or (6) for each transition, where ca gives the relative absorption contribution of level a at the wavelength of interest, i.e., aca=1. This is equivalent to calculating the effective pe from a weighted average,

(8)

where the pe,a for each of the absorption levels a has been separately computed, e.g., from Eq. (7b). Apart from Γga, pe,a only depends on the scattering rates rij, which also determine the associated stationary level occupations pi through Eq. (4). Hence, since typically pa0 or Γga1 for QCDs as mentioned above, pe,a in Eq. (8) is wavelength independent. However, this does not apply to pe if there are two or more relevant absorption levels, since the absorption coefficient of each transition, and thus also ca, is wavelength dependent.

For noise modeling, the QCD structure can be represented by a resistance network as illustrated in Fig. 1(b), where the subbands become the nodes, and the intersubband transitions are represented by equivalent resistances that generate the same noise as caused by the statistical fluctuations of the carriers.20 For unbiased QCD structures, the Johnson noise dominates at elevated temperatures,20 and the conductance per unit area of the resistances is given by

(9)

Here, ns is the sheet doping density per period. The total detector resistance Rd in Eq. (3) can then be calculated by solving Kirchhoff’s equations for the equivalent resistance network. Assuming a network of conductances and setting σii=0, we can formulate Kirchhoff’s law for any node i as

(10)

As illustrated in Fig. 1(b), the QCD is represented by a periodic network with σij=σi+N,j+N. In the following, we exploit the periodicity of the network, which greatly reduces the numerical cost. To this end, we assume that the voltages of equivalent nodes in adjacent periods are related by ui+N=ui+up, where up is the voltage drop per period. This applies if the structure contains many periods so that the boundary conditions at the end contacts can be neglected. Formulating Eq. (10) for a reference period and taking into account transitions to and from the left- and the right-neighboring period, we then arrive at

(11)

Here, the conductance matrix elements are given by

and the elements of d are

A model equivalent to Eqs. (10) and (11) was previously used to assess the QCD noise due to longitudinal-optical (LO) phonon scattering.21 Taking i,j=1,,N, the N linear equations represented by Eq. (11) are not independent, which reflects the fact that there is no unique solution for u since the zero potential can be arbitrarily chosen. Thus, we set uN=0 and eliminate the Nth equation from Eq. (11), i.e., we now have i,j=1,,(N1), and the order of G and length of u and d are reduced to N1. As for Eq. (6), Cramer’s rule can be applied to Eq. (11) to write the ui in explicit form, which yields

(12)

where Gi is the matrix obtained by replacing the ith column of G by vector d.

The conductance per unit area of a single QCD period can be calculated from the ratio of net current density flowing between adjacent periods to the voltage drop per period,

(13a)
(13b)

where Eq. (13b) has been obtained by inserting Eq. (12). From Eq. (13b), we see that σp only depends on the σij, which can with Eqs. (4) and (9) be computed directly from the EMC intersubband scattering rates. In particular, σp is independent of the chosen value for up, which approaches zero for QCDs in the photovoltaic mode. The detector resistance in Eq. (3) is given by

where A is the in-plane cross section area of the quantum wells, which coincides with the detector area Ad in Eq. (3) if the facet is in in-plane direction.

Photovoltaic QCDs are not biased and are thus in thermodynamic equilibrium if there is no incident light. As discussed in Sec. II, the figures of merit in the linear response regime can then be obtained from the corresponding scattering rates rij and associated subband occupation probabilities pi. These probabilities, as well as the electron distribution within each subband, are also governed by thermodynamic equilibrium relations, which can directly be exploited for QCD modeling.21,22 In a QW structure, each subband i is formed by electron states |i,k, where the quasi-free electron motion in the QW planes is described by the two-dimensional wavevector k with the associated kinetic energy εi(k), and Ei is the quantized energy in growth direction. The occupation probability of each state is given by the Fermi–Dirac distribution function

(14)

Quantum cascade devices typically utilize energy states in the conduction band Γ valley, where the dispersion relation can for moderate kinetic electron energies εi be approximated by εi=2|k|2/(2mi). Here, nonparabolicity effects are to first order considered by assigning a different effective mass mi to each subband.11 This model is schematically illustrated in Fig. 2. The electron sheet density in a subband i is then obtained from Eq. (14) by multiplying fi with the density of states niE=mi/(π2) and integrating over εi,11 which yields

(15)

Equation (15) only depends on Ei, mi, T, and the chemical potential μ, which can be computed from the sheet doping density per period ns via ns=i=1Nns,i.11 The thermal occupation probabilities are then obtained as

(16)
FIG. 2.

Schematic illustration of the in-plane dispersion relation for ground and absorption subband and of the electron distribution function f. The density of states in each subband is represented by the density of horizontal lines. The arrow indicates an optical transition.

FIG. 2.

Schematic illustration of the in-plane dispersion relation for ground and absorption subband and of the electron distribution function f. The density of states in each subband is represented by the density of horizontal lines. The arrow indicates an optical transition.

Close modal

A further important aspect of device modeling at thermodynamic equilibrium is that the principle of detailed balance holds for the scattering rates between each pair of states and thus also between each pair of subbands i and j,10,23,24

(17)

As mentioned above, the correction factor Γga in Eqs. (6) and (7b) for optical transitions between two subbands g and a is in general not exactly unity. The total number of stimulated emission (i=a, j=g) and absorption (i=g, j=a) events per unit time scales as 0γag(γag2+Ω2)1fi(εi)[1fj(εj)]dε, where we use as integration variable the kinetic energy ε=εg in the ground subband and thus have εa=εmg/ma.24–26 Here, Ω=ω[Eag+(mg/ma1)ε]/ with Eag=EaEg, and also the dephasing rate γag generally depends on ε. Γga can then be expressed as

(18)

For elevated temperatures and moderate doping densities, the Pauli blocking factor 1fj in Eq. (18) can be neglected, and fi given by Eq. (14) can be approximated by a Maxwell–Boltzmann distribution fi(pi/mi)exp[εi/(kBT)].25 In the case of parallel subbands, i.e., mg=ma, Eq. (18) then yields Γga=1. This is typically a good approximation for terahertz QCDs, where nonparabolicity effects tend to be negligible due to the small energetic spacing Eag of the optical transition. On the other hand, for mid-infrared QCDs, nonparabolicity must be considered to obtain realistic modeling results.27 Since the effective mass increases with the subband energy, we assume ma>mg in Fig. 2. Γga accounts for two factors: stimulated photon emission from subband a to g sees a lower density of final electron states than the associated absorption transition, as illustrated in Fig. 2. On the other hand, the occupation probability fi[2|k|2/(2mi)] decreases more slowly as a function of |k| in the upper than in the lower subband,25 as also becomes apparent from Fig. 2. Thus, the bigger |k|, i.e., the smaller ω, the greater the ratio of stimulated photon emission to absorption. However, since for mid-infrared QCDs, ωkBT and thus papg even at elevated temperatures, the term Γgapa can be neglected in Eqs. (6) and (7b).

The EMC method is widely used for the simulation of quantum cascade structures.11,28 In EMC, scattering is self-consistently evaluated based on Fermi’s golden rule.11 Our code includes all relevant scattering mechanisms, such as the interaction of electrons with other electrons as well as LO and acoustic phonons, and inelastic scattering due to interface roughness, alloy disorder and impurities. Thus, the simulation does not rely on empirical or fitting parameters but only uses the design specifications such as the layer widths and material compositions along with the corresponding material parameters. Also, nonparabolicity is considered through energy-dependent effective masses in the Schrödinger–Poisson solver, which results in modified subband eigenenergies and wavefunctions, as well as in different subband effective masses mi.11 

As mentioned above, a special challenge of photovoltaic QCD simulations is that the currents are much smaller than in QCLs, impeding a direct simulation of the carrier transport in the irradiated structure due to the limited numerical accuracy. This problem can be circumvented by computing the QCD characteristics from the simulated absorption coefficient and scattering rates rather than from the output current, as described in Sec. II. Here, it must be taken into account that the thermal occupation probabilities pi of the high-lying subbands are very small, which leads to difficulties in EMC since the pi and corresponding scattering rates rij cannot be accurately evaluated due to insufficient statistics in the stochastic sampling process. For low temperatures, it even happens regularly that not a single electron enters the highest-lying subbands within reasonable simulation times due to their extremely low occupation probabilities. However, the computation of the detector figures of merit in Sec. II requires highly accurate values of pi and rij. This issue can, for example, be addressed by employing a weighted EMC scheme, where certain scattering rates are artificially enhanced in order to reduce the statistical fluctuations.24 Here, we choose a different approach, taking advantage of the thermodynamic equilibrium in the system. In particular, we compute the pi required in Eqs. (6)–(9) directly from Eq. (16), rather than extracting the occupations from the EMC simulation. In fact, we also fix the electron distribution in EMC to a Fermi–Dirac statistics, albeit with artificially increased upper-level occupations to enhance the accuracy of the obtained scattering rates. Thus, after a scattering process has occurred and been tabulated in EMC, the electron involved is randomly assigned to a subband according to the probability distribution given by Eq. (16), where, however, a lower bound of pmin=0.002 is imposed on the occupation probabilities, yielding after normalization modified values pi for all subbands. Here, pmin is chosen so that sufficient sampling statistics is obtained also for scattering from the high-lying states but occupation dependent contributions to scattering, such as Pauli exclusion or screening, are not significantly altered. Subsequently, the kinetic energy εi of the redistributed electron is randomly chosen in EMC by evaluating Eq. (14) with the inverse transform sampling method, which yields

(19)

Here, ξ is a random variable, uniformly distributed between 0 and 1.

The (single-electron) scattering rates rij are evaluated in the usual way by counting the tabulated transition events from subband i to j over a given time interval and normalizing the result to the number of electrons in the initial level i. In this way, the increased number of scattering events from high-lying subbands i with artificially enhanced pi is automatically compensated for, yielding the correct values for rij. In an analogous way, the k resolved scattering rates, required for calculating the absorption spectrum,16,17 are obtained. Special care has to be taken for electron–electron scattering: in its typical implementation as a two-body process,11,29 wrong values for rij are obtained since also the randomly selected scattering partner from the electron ensemble obeys the modified probability distribution pi. We avoid this issue by assigning a virtual scattering partner with (unmodified) thermal distribution pi and tabulating only the scattering process of the primary electron but not its virtual scattering partner (the correction of the scattering rate by a factor 1/2 to avoid double counting does then not apply30,31).

Due to remaining fluctuations of the sampling process and other numerical limitations, the extracted rates rij do not exactly fulfill the detailed balance principle, Eq. (17), which is not a shortcoming specific to EMC but constitutes a more general issue in numerical modeling.23 Consequently, Eq. (4) does not strictly hold if the exact thermal occupation probabilities pi computed from Eq. (16) are combined with the EMC scattering rates, which also affects the accuracy of Eqs. (6) and (7) derived from Eq. (4). Thus, we strictly enforce observance of the detailed balance principle in our simulations. From Eqs. (16) and (17), it follows that rij>rji for Ei>Ej, and thus rij tends to be determined more accurately in EMC due to the lower statistical error. Hence, for the calculation of pe and D, the rate rij from the higher-energy subband i to the lower-energy subband j is taken from EMC, while the back-scattering rate rji is computed from Eq. (17) to avoid violations of the detailed balance principle resulting from numerical inaccuracies.

In the following, we present the simulation results of QCDs in the mid-infrared and terahertz regime. The two simulated devices detecting in the mid-infrared region have their absorption maxima at wavelengths of 4.7 and 7.5 μm, respectively. Both devices are based on the InGaAs/InAlAs material system, lattice matched to the InP substrate.12,15,32,33 The terahertz structure consists of a GaAs/AlGaAs based superlattice structure on a semi-insulating GaAs substrate and exhibits a detection wavelength of 84 μm.19 

The investigated devices in the mid-infrared range have the sample identifier names N1021 (7.6 μm) and N1022 (4.7 μm). The active region of both designs consists of 30 periods. In this configuration, a 45° wedge geometry mesa-structure with a double-pass waveguide is considered for the simulation. The facet transmittance is assumed to be 70%. Further information about these designs, along with experimental data used for comparison purposes in this work, can be found in the literature, where also details on the measurements are given.12,15,32,33

The strongly delocalized Bloch-type wavefunction solutions of the Schrödinger equation for perfectly periodic potentials are not adequate for describing carrier transport in realistic QCD structures since any amount of disorder due to, e.g., growth irregularities and impurities, leads to the formation of localized states.34,35 In addition, scattering-induced dephasing suppresses multiple-tunneling processes and thus also contributes to a localization of states.36,37 Hence, we restrict the spatial simulation window in our Schrödinger–Poisson solver such that wavefunction solutions extending over multiple periods are avoided. The resulting conduction band profile and energy eigenstates are shown in Fig. 3 for the QCD N1021. The lattice constant of InGaAs at room temperature is 0.587 nm, and the thinnest well of the structure has a width of 2.8 nm. Interface roughness is included in our carrier transport simulations using typical values of 0.1 nm for the average root mean square roughness height and 10 nm for the in-plane correlation length. Absorption takes place in the active well of each period between the ground level g and the two closely spaced absorption levels a1 and a2. As explained in Sec. II A, the weighted average for the extraction efficiency pe of the transitions ga1 and ga2 is calculated using Eq. (8). Both absorption levels extend over two wells and give rise to resonant tunneling transport through the thick barrier between the active well and the subsequent well. Transport through the cascade is mediated by LO phonon assisted scattering and to a lesser extent by interface roughness scattering. Efficient extraction is achieved by arranging the energy spacing between two consecutive states close to the LO phonon energy of 32 meV in In0.53Ga0.47As. More specifically, the quantum well structure is engineered to optimize detectivity, which involves a trade-off between close to unity extraction efficiency pe and high detector resistance Rd. For the design N1021, the extractor states are localized in their respective wells. Since in this case, the overlap between adjacent states mainly results from the wavefunction leakage into the thick barrier regions, the corresponding scattering rates tend to be considerably lower than in QCLs, which typically feature thinner barriers and more delocalized wavefunctions.

FIG. 3.

Calculated conduction band profile and probability densities of the investigated mid-infrared QCD structure N1021 detecting at 7.6 μm.32,33 Starting from the leftmost well, the layer sequence (in nm) of one period with barriers in boldface and n-doped layers (3×1017cm3) underlined is 8.2_/6.0/2.8/5.7/3.4/5.5/4.2/5.1/5.5/5.8. Space charge effects are included by self-consistently solving the Schrödinger and Poisson equation. The electron densities in the states are modeled accounting for thermal equilibrium under zero external bias and no incident light for a lattice temperature of 300 K.

FIG. 3.

Calculated conduction band profile and probability densities of the investigated mid-infrared QCD structure N1021 detecting at 7.6 μm.32,33 Starting from the leftmost well, the layer sequence (in nm) of one period with barriers in boldface and n-doped layers (3×1017cm3) underlined is 8.2_/6.0/2.8/5.7/3.4/5.5/4.2/5.1/5.5/5.8. Space charge effects are included by self-consistently solving the Schrödinger and Poisson equation. The electron densities in the states are modeled accounting for thermal equilibrium under zero external bias and no incident light for a lattice temperature of 300 K.

Close modal

The simulated responsivity spectrum for the detector structure N1021 is depicted in Fig. 4. At 300 K, we obtain a peak responsivity of 0.61 mA W1, in comparison to a measured value of 1.69 mA W1. For lower temperatures, the simulated and measured responsivities show considerably better agreement. The simulated peak responsivity increases to 3.1 mA W1 for 100 K, as compared to a measured value of 3.8 mA W1. At a temperature of 300 K, the simulated peak detection wavelength is 7.8 μm in comparison to 7.62 μm in the experiment. This slight blueshift of the experimental detection energy was also recognized in the design process15 and is attributed to process uncertainties in the fabrication, which can lead to thickness variations of the active QW. The observed redshift with increasing operation temperature, which is attributed to nonparabolicity and band filling, is replicated in the simulations by the inclusion of energy-dependent effective subband masses.

FIG. 4.

Simulated responsivity Rp of the QCD structure N1021 as a function of wavenumber in the temperature range of 100–300 K.

FIG. 4.

Simulated responsivity Rp of the QCD structure N1021 as a function of wavenumber in the temperature range of 100–300 K.

Close modal

To validate our detector model, we evaluate the simulation results for the extraction efficiency pe, power absorption efficiency ηabs and device resistance Rd, and compare them with the experimentally measured values. Results are illustrated in Fig. 5 for 100–300 K.

FIG. 5.

Comparison of simulated (x marks) and experimental (squares) results for the detector structure N1021. (a) Extraction efficiency pe (solid lines) and peak absorption ηabs,p (dashed lines) as a function of temperature T. (b) Detector resistance-area product AdRd as a function of the inverse temperature T1. The lines represent the Arrhenius plots with calculated activation energies Eact.

FIG. 5.

Comparison of simulated (x marks) and experimental (squares) results for the detector structure N1021. (a) Extraction efficiency pe (solid lines) and peak absorption ηabs,p (dashed lines) as a function of temperature T. (b) Detector resistance-area product AdRd as a function of the inverse temperature T1. The lines represent the Arrhenius plots with calculated activation energies Eact.

Close modal

The power absorption efficiency in Fig. 5(a) decreases linearly from 20% at 100 K to 11% at 300 K for the simulated results and is 3% smaller than in experiment. Due to broadening effects and reduced population of the ground level g with rising temperatures, the absorption efficiency and thus the responsivity get reduced. For decreasing temperature, a shift of the probability density maximum of state a1 to the right thinner QW and of state a2 to the left active QW occurs due to the change in band bending associated with space charge effects. The peak absorption of the transition ga2 thus increases significantly and becomes the dominating absorption transition. The responsivity depends furthermore on the extraction efficiency pe, which is also displayed in Fig. 5(a). Over the plotted range, the simulated pe decreases linearly with increasing temperature, whereas the experimental values saturate to pe,300K=5.9% for temperatures above 200 K. Evaluation of the simulated scattering rates shows that LO phonon emission is the dominant mechanism for relaxation from the two absorption levels a1,a2 to the ground level g. Extraction to the cascade is dominated by interface roughness and LO phonon scattering. Furthermore, the calculation of the individual extraction efficiencies in one period exhibits that backscattering in the extraction cascade affects significantly the behavior of pe over temperature. Backscattering by LO phonon absorption degrades due to the thermal reduction of the phonon occupation for low temperatures. By reducing temperature, backscattering in the extraction cascade gets suppressed, which results in an efficient charge extraction to the ground state of the adjacent period.11 Last, the device resistance Rd is presented in Fig. 5(b). Experimentally, the device resistance Rd was determined from dark current measurements at zero bias. In the temperature range of 100–300 K, the exponential behavior of the detector resistance is replicated. The calculated activation energy of Eact=135 meV is in good agreement with the experimentally obtained value of Eact=129 meV. The close correspondence of the experimentally obtained resistance values with the simulated ones confirms the validity of the resistance model presented in Sec. II B for the investigated temperature range.

The second mid-infrared detector N1022 exhibits a similar layer composition as N1021. Here, the thickness of the active well is reduced from 82 Å (N1021) to 51 Å, and the cascade is extended by three extractor levels. The structure modifications result in a higher detecting energy of 268 meV. The peak responsivity Rp and detector resistance-area product AdRd for the temperature range of 100–300 K are presented in Figs. 6(a) and 6(b), respectively. The simulated peak responsivity Rp decreases from 5.3 mA W1 at 100 K to 1.0 mA W1 at 300 K similar to the experimental data. At 300 K, the values pe=17.2% and ηabs=4.63% are obtained from the simulated scattering rates, which is in good agreement with the experimentally obtained values of pe=12.5% and ηabs=9% at room temperature. As illustrated in Fig. 6(b), the device resistance shows the Arrhenius plot behavior. The high-temperature activation energy of Eact=237 meV calculated from the simulated detector resistance values compares well with the experimentally obtained value of Eact=210 meV. Below 150 K, the experimental values for the resistance-area product deviate from the Arrhenius plot slope and exhibit a flattening trend, which cannot be replicated in the simulation. This behavior is ascribed in the literature to the presence of a parasitic parallel resistance, generated, e.g., by defects or surface currents.14 

FIG. 6.

Comparison of simulated (x marks) and experimental (squares) results for the predefined figures of merit of the detector structure N1022. (a) Peak responsivity Rp as a function of temperature T. (b) Detector resistance-area product AdRd as a function of the inverse temperature T1. The lines represent the Arrhenius plots with calculated activation energies Eact.

FIG. 6.

Comparison of simulated (x marks) and experimental (squares) results for the predefined figures of merit of the detector structure N1022. (a) Peak responsivity Rp as a function of temperature T. (b) Detector resistance-area product AdRd as a function of the inverse temperature T1. The lines represent the Arrhenius plots with calculated activation energies Eact.

Close modal

Furthermore, our approach has been applied to a GaAs/AlGaAs based terahertz QCD operating at a wavelength of around 84μm.19 Here, we assume again a 45° double-pass waveguide mesa configuration with a transmittance at the air/GaAs interface of 70%.

In Fig. 7, the conduction band diagram and energy eigenstates obtained with a Schrödinger–Poisson solver for a lattice temperature of 10 K are shown. Here, we choose a conduction band offset of 135 meV for the calculation of the wavefunctions. The thinnest barrier has a width of 1.5 nm, which is approximately three lattice constants (0.565 nm for GaAs). Figure 8 displays the simulated and measured spectral responsivity at a temperature of 10 K, showing good agreement. Our simulation yields a peak responsivity of 10.6 mA W1 at 109 cm1, which compares well with the experimental value of 8.6 mA W1 at 119 cm1. The simulated extraction efficiency of pe=0.36 is in reasonable agreement with the experimental estimate of pe0.5, which was obtained under idealized assumptions.19 For the two relevant transitions ga1 and ga2, absorption efficiencies of 2.4% and 3.8% are obtained from the simulation, which agrees reasonably well with the estimate of 3% for both transitions, inferred from the calculated dipole matrix elements.19 However, the experimentally measured absorption values were below 2%, which was in part attributed to trapping of electrons by impurities.19 In addition to experimental uncertainties, certain effects may be relevant which are not considered in the simulation, in part because their quantitative influence is not well known. For example, the measured photocurrent and, therefore, the derived responsivity depends on the extraction from the device contacts, which are not included in the model.

FIG. 7.

Calculated conduction band profile and probability densities of the investigated terahertz QCD structure detecting at 84 μm for a lattice temperature of 10 K. Starting from the leftmost well, the layer sequence (in nm) of one period with barriers in boldface and n-doped layers (6×1015cm3) underlined is 18.0/3.4/21.0/1.5/10.0/4.5/11.0/4.0/13.5/3.8/14.5_/3.6.

FIG. 7.

Calculated conduction band profile and probability densities of the investigated terahertz QCD structure detecting at 84 μm for a lattice temperature of 10 K. Starting from the leftmost well, the layer sequence (in nm) of one period with barriers in boldface and n-doped layers (6×1015cm3) underlined is 18.0/3.4/21.0/1.5/10.0/4.5/11.0/4.0/13.5/3.8/14.5_/3.6.

Close modal
FIG. 8.

Responsivity spectrum Rp as a function of wavenumber for the 3.5 THz QCD structure at a temperature of 10 K, as obtained from the simulation and extracted from experimental photocurrent measurements. Furthermore, a comparison between EMC and NEGF simulation results is shown for 50 K.

FIG. 8.

Responsivity spectrum Rp as a function of wavenumber for the 3.5 THz QCD structure at a temperature of 10 K, as obtained from the simulation and extracted from experimental photocurrent measurements. Furthermore, a comparison between EMC and NEGF simulation results is shown for 50 K.

Close modal

Furthermore, we have compared the EMC results to simulation results obtained by a non-equilibrium Green’s function (NEGF) approach,38–40 where the responsivity is calculated by simulating the photocurrent under irradiation. This model has been developed for simulating QCLs under operating conditions, where small (several A/cm2 ) uncertainties in the current density due to the limited number of periods and basis states considered are acceptable. However, the dark currents and photocurrents in QCDs are on a similar order of magnitude as this uncertainty. Therefore, we have subtracted from the photocurrent a constant offset defined as the photocurrent at a high frequency (160 cm1) where the simulated absorption is close to zero. The obtained responsivity spectra of both the EMC and NEGF simulations for a lattice temperature of 50 K are illustrated in Fig. 8. The two responsivity peaks of 10.4 mA W1 at 109 cm1 and 14.6 mA W1 at 129 cm1 from the NEGF simulations are in good agreement with the EMC results of 7.2 mA W1 at 110 cm1 and 6.7 mA W1 at 129 cm1. Notably, both the EMC and NEGF spectra are red-shifted with respect to the measured one. The depolarization shift, not considered in our simulations, is estimated to be 0.1meV41 and is thus too small to account for the discrepancy. Another explanation might be that similar fabrication uncertainties as described for the mid-infrared designs lead to slightly higher measured peak absorption energies. In general, the peak responsivity decreases with increasing temperature due to the thermal activation of electrons and the resulting occupation of higher lying states. According to the EMC simulations, this results for 10–50 K in a peak responsivity reduction of 32.1% at 109 cm1 and 53.5% at 129 cm1, respectively.

Figure 9 shows the simulated current–voltage characteristics of the terahertz QCD without illumination. The simulated results of the NEGF and EMC models are compared to the experimentally measured values for the three temperatures 50, 100, and 150 K. As mentioned above, the currents in QCDs are significantly smaller than in QCLs, which makes the dark current simulation of such terahertz QCD devices with NEGF and EMC at small biases rather unpractical. In fact, the NEGF simulations only converged for certain bias points and only showed robust results at 50 K when electron–electron interaction was included.42 As illustrated in Fig. 9, the experimentally measured dark currents for applied biases below 0.2 V cannot be reproduced by both simulation approaches. This problem becomes even more manifest for temperatures of 50 K and below, resulting in even smaller dark current densities. For higher biases, both simulation models can replicate the measured current–voltage characteristics.

FIG. 9.

The measured current–voltage characteristics of the THz detector design (dots) is compared to results obtained from EMC (squares) and NEGF including electron–electron scattering (crosses) simulations for 50, 100, and 150 K. Fits to the experimental results (solid lines) as well as to the simulation results for EMC (dashed lines) and NEGF (dotted lines) are included using a polynomial function of second degree.

FIG. 9.

The measured current–voltage characteristics of the THz detector design (dots) is compared to results obtained from EMC (squares) and NEGF including electron–electron scattering (crosses) simulations for 50, 100, and 150 K. Fits to the experimental results (solid lines) as well as to the simulation results for EMC (dashed lines) and NEGF (dotted lines) are included using a polynomial function of second degree.

Close modal

From the experimentally measured dark current at 50mV,19 we estimate a zero-bias resistance-area product RdAd=300Ωcm2 at 10 K under the assumption that this bias is already in the linear regime. The simulated resistance-area products of RdAd=3Ωcm2 (EMC) and RdAd=1Ωcm2 (NEGF) are two orders of magnitude smaller. For temperatures below 20 K, the dark current is expected to be dominated by a temperature independent contribution based on direct tunneling processes.19 This effect is not included in the EMC approach but fully accounted for in the NEGF simulations, thus it does not explain the observed discrepancy between theory and experiment. In Fig. 10, the temperature dependent resistance-area product RdAd extracted from the measured current–voltage characteristics and obtained with EMC and NEGF is displayed logarithmically as a function of the inverse temperature T1.

FIG. 10.

Detector resistance-area product AdRd as a function of the inverse temperature T1, as obtained from EMC and NEGF simulations and extracted from experimental dark current measurements.

FIG. 10.

Detector resistance-area product AdRd as a function of the inverse temperature T1, as obtained from EMC and NEGF simulations and extracted from experimental dark current measurements.

Close modal

A simulated specific detectivity D=1.05×109cmHz/W is calculated from the EMC scattering rates at 10 K. By substituting the experimentally measured responsivity Rp=8.6mAW1 and resistance-area product into Eq. (3), we obtain a specific detectivity D=6.3×109cmHz/W at 10 K. Experimentally, the specific detectivity was measured using a QCL with an emitting wavelength of 87 μm43 to estimate the noise equivalent power (NEP) of the THz QCD. With the measured NEP=7 nW, a specific detectivity of 5×107cmHz/W at 10 K was obtained,19 which is two orders of magnitude smaller than the value obtained from the dark current measurement. As it was pointed out above, the dominant noise mechanism for temperatures below 20 K is not temperature dependent and thus the simulated detectivity based on our Kirchhoff’s resistance model, taking into account Johnson noise as main noise mechanism, cannot replicate the experimentally measured specific detectivities. The resulting discrepancy between the measured specific detectivity and the one deduced from the dark current measurements implies a second noise source appearing in the photocurrent measurements, e.g., additional blackbody radiation entering through the cryostat window of the measurement setup.

In conclusion, we have presented an EMC-based modeling approach for the simulation of photovoltaic QCD operation. Our method takes advantage of the transport properties and carrier occupations in thermodynamic equilibrium, thus increasing the simulation accuracy and reducing the required runtime. Exploiting the periodicity of the QCD structure, compact expressions are derived for the extraction efficiency and device resistance, enabling the computation of the responsivity and specific detectivity directly from the EMC scattering rates and the related absorption coefficient. All in all, comparison with experimental data yields good agreement for both mid-infrared and terahertz QCD designs. For the investigated terahertz structure, the EMC results are furthermore compared to NEGF simulations, yielding good agreement. Notably, for low temperatures, the simulated zero-bias resistance of this structure is much lower than the value extracted from experiment. Due to its accuracy, versatility, and relative numerical efficiency, the presented simulation approach is well suited for the systematic optimization of QCD structures and also allows investigating the robustness of a design against fabrication tolerances by stochastic sampling of points in the design parameter space.44 

The authors acknowledge financial support by the European Union’s Horizon 2020 research and innovation programme under Grant Agreement No. 820419—Qombs Project “Quantum simulation and entanglement engineering in quantum cascade laser frequency combs” (FET Flagship on Quantum Technologies).

The data that support the findings of this study are available from the corresponding author upon reasonable request.

The authors have no conflicts to disclose.

1.
B.
Levine
,
K.
Choi
,
C.
Bethea
,
J.
Walker
, and
R.
Malik
,
Appl. Phys. Lett.
50
,
1092
(
1987
).
2.
H.
Schneider
,
P.
Koidl
,
F.
Fuchs
,
B.
Dischler
,
K.
Schwarz
, and
J.
Ralston
,
Semicond. Sci. Technol.
6
,
C120
(
1991
).
3.
H.
Schneider
,
C.
Schönbein
,
M.
Walther
,
K.
Schwarz
,
J.
Fleissner
, and
P.
Koidl
,
Appl. Phys. Lett.
71
,
246
(
1997
).
4.
J.
Faist
,
F.
Capasso
,
D. L.
Sivco
,
C.
Sirtori
,
A. L.
Hutchinson
, and
A. Y.
Cho
,
Science
264
,
553
(
1994
).
5.
C.
Deutsch
,
H.
Detz
,
T.
Zederbauer
,
A. M.
Andrews
,
P.
Klang
,
T.
Kubis
,
G.
Klimeck
,
M. E.
Schuster
,
W.
Schrenk
,
G.
Strasser
, and
K.
Unterrainer
,
Opt. Express
21
,
7209
(
2013
).
6.
D.
Hofstetter
,
M.
Beck
, and
J.
Faist
,
Appl. Phys. Lett.
81
,
2683
(
2002
).
7.
L.
Gendron
,
M.
Carras
,
A.
Huynh
,
V.
Ortiz
,
C.
Koeniguer
, and
V.
Berger
,
Appl. Phys. Lett.
85
,
2824
(
2004
).
8.
A.
Harrer
,
B.
Schwarz
,
S.
Schuler
,
P.
Reininger
,
A.
Wirthmüller
,
H.
Detz
,
D.
MacFarland
,
T.
Zederbauer
,
A.
Andrews
,
M.
Rothermund
,
H.
Oppermann
,
W.
Schrenk
, and
G.
Strasser
,
Opt. Express
24
,
17041
(
2016
).
9.
O.
Baumgartner
,
Z.
Stanojevic
,
K.
Schnass
,
M.
Karner
, and
H.
Kosina
,
J. Comput. Electron.
12
,
701
(
2013
).
10.
C.
Koeniguer
,
G.
Dubois
,
A.
Gomez
, and
V.
Berger
,
Phys. Rev. B
74
,
235325
(
2006
).
11.
C.
Jirauschek
and
T.
Kubis
,
Appl. Phys. Rev.
1
,
011307
(
2014
).
12.
F.
Giorgetta
, “Design, fabrication, and testing of intersubband infrared photodetectors operating at wavelengths between 2 μm and 17 μm,” Ph.D. thesis (School Université de Neuchâtel, 2007).
13.
A.
Harrer
,
B.
Schwarz
,
R.
Gansch
,
P.
Reininger
,
H.
Detz
,
T.
Zederbauer
,
A. M.
Andrews
,
W.
Schrenk
, and
G.
Strasser
,
Appl. Phys. Lett.
105
,
171112
(
2014
).
14.
B.
Schwarz
, “Monolithic integration of mid-infrared photonics,” Ph.D. thesis (School Technische Universität Wien, 2015).
15.
F. R.
Giorgetta
,
E.
Baumann
,
M.
Graf
,
Q.
Yang
,
C.
Manz
,
K.
Kohler
,
H. E.
Beere
,
D. A.
Ritchie
,
E.
Linfield
,
A. G.
Davies
et al.,
IEEE J. Quantum Electron.
45
,
1039
(
2009
).
16.
C.
Jirauschek
and
P.
Lugli
,
J. Appl. Phys.
105
,
123102
(
2009
).
17.
C.
Jirauschek
,
J. Appl. Phys.
122
,
133105
(
2017
).
18.
C.
Jirauschek
,
IEEE Photonics J.
10
,
1
(
2018
).
19.
M.
Graf
,
G.
Scalari
,
D.
Hofstetter
,
J.
Faist
,
H.
Beere
,
E.
Linfield
,
D.
Ritchie
, and
G.
Davies
,
Appl. Phys. Lett.
84
,
475
(
2004
).
20.
A.
Delga
,
L.
Doyennette
,
M.
Carras
,
V.
Trinité
, and
P.
Bois
,
Appl. Phys. Lett.
102
,
163507
(
2013
).
21.
A.
Buffaz
,
A.
Gomez
,
M.
Carras
,
L.
Doyennette
, and
V.
Berger
,
Phys. Rev. B
81
,
075304
(
2010
).
22.
S.
Saha
and
J.
Kumar
,
Superlattices Microstruct.
98
,
70
(
2016
).
23.
W. R.
Frensley
,
Rev. Mod. Phys.
62
,
745
(
1990
).
24.
R. C.
Iotti
and
F.
Rossi
,
Rep. Prog. Phys.
68
,
2533
(
2005
).
25.
V. B.
Gorfinkel
,
S.
Luryi
, and
B.
Gelmont
,
IEEE J. Quantum Electron.
32
,
1995
(
1996
).
26.
C.
Jirauschek
,
Appl. Phys. Lett.
96
,
011103
(
2010
).
27.
O.
Baumgartner
,
Z.
Stanojević
, and
H.
Kosina
, in Book of Abstracts of the 16th International Workshop on Computational Electronics (IWCE) (Society for Micro- and Nanoelectronics, 2013), p. 86.
28.
P.
Borowik
,
J.-L.
Thobel
, and
L.
Adamowicz
,
Opt. Quant. Electron.
49
,
96
(
2017
).
29.
S. M.
Goodnick
and
P.
Lugli
,
Phys. Rev. B
37
,
2578
(
1988
).
30.
M.
Moško
and
A.
Mošková
,
Phys. Rev. B
44
,
10794
(
1991
).
31.
X.
Gao
,
D.
Botez
, and
I.
Knezevic
,
J. Appl. Phys.
101
,
063101
(
2007
).
32.
D.
Hofstetter
,
F. R.
Giorgetta
,
E.
Baumann
,
Q.
Yang
,
C.
Manz
, and
K.
Köhler
,
Appl. Phys. B
100
,
313
(
2010
).
33.
D.
Hofstetter
,
F. R.
Giorgetta
,
E.
Baumann
,
Q.
Yang
,
C.
Manz
, and
K.
Köhler
, in Proceedings of SPIE (International Society for Optics and Photonics, 2010), Vol. 7608, p. 76081N.
35.
J.
Palmier
,
H.
Le Person
,
C.
Minot
,
A.
Chomette
,
A.
Regreny
, and
D.
Calecki
,
Superlattices Microstruct.
1
,
67
(
1985
).
36.
H.
Callebaut
and
Q.
Hu
,
J. Appl. Phys.
98
,
104505
(
2005
).
37.
H.
Willenberg
,
G. H.
Döhler
, and
J.
Faist
,
Phys. Rev. B
67
,
085315
(
2003
).
38.
A.
Wacker
,
M.
Lindskog
, and
D. O.
Winge
,
IEEE J. Sel. Top. Quantum Electron.
19
,
1
(
2013
).
39.
M.
Lindskog
,
J.
Wolf
,
V.
Trinite
,
V.
Liverini
,
J.
Faist
,
G.
Maisons
,
M.
Carras
,
R.
Aidam
,
R.
Ostendorf
, and
A.
Wacker
,
Appl. Phys. Lett.
105
,
103106
(
2014
).
40.
M.
Franckié
,
L.
Bosco
,
M.
Beck
,
C.
Bonzon
,
E.
Mavrona
,
G.
Scalari
,
A.
Wacker
, and
J.
Faist
,
Appl. Phys. Lett.
112
,
021104
(
2018
).
41.
G.
Pegolotti
,
A.
Vasanelli
,
Y.
Todorov
, and
C.
Sirtori
,
Phys. Rev. B
90
,
035305
(
2014
).
42.
D. O.
Winge
,
M.
Franckié
,
C.
Verdozzi
,
A.
Wacker
, and
M. F.
Pereira
,
J. Phys. Conf. Ser.
696
,
012013
(
2016
).
43.
G.
Scalari
,
L.
Ajili
,
J.
Faist
,
H.
Beere
,
E.
Linfield
,
D.
Ritchie
, and
G.
Davies
,
Appl. Phys. Lett.
82
,
3165
(
2003
).
44.
J.
Popp
,
M.
Haider
,
M.
Franckié
,
J.
Faist
, and
C.
Jirauschek
,
Opt. Quant. Electron.
53
,
287
(
2021
).