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.
I. INTRODUCTION
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 to an absorption level . 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.
(a) Schematic illustration of a QCD in photovoltaic operation mode. The rectangle denotes a single QCD period. (b) Equivalent circuit diagram for noise modeling.
(a) Schematic illustration of a QCD in photovoltaic operation mode. The rectangle denotes a single QCD period. (b) Equivalent circuit diagram for noise modeling.
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.
II. FIGURES OF MERIT
A frequently used figure of merit for characterizing photodetector performance is the responsivity , i.e., the generated photocurrent divided by the incident light power . 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 and extraction efficiency . For a structure with periods, is then given by8,12
Here, , , and 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
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 denotes the number of passes of the optical field through the active region. In particular, for a single-pass geometry, and for a double-pass geometry, as commonly realized by exploiting internal reflection at the top contact. Higher values of can, for example, be obtained by using a multipass waveguide geometry.12,15 is the facet transmittance and can potentially also include losses due to the substrate and contact layers. In order to optimize , Fresnel reflections are often minimized by polishing the entrance facet at an angle of .12,15 Furthermore, denotes the power absorption coefficient, and is the length of a period so that 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 , 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 , which gives the signal (per unit incident power) to noise ratio, normalized to the measurement electrical bandwidth and the detector area . Here, is the root mean square noise current. For photovoltaic (i.e., unbiased) detectors where Johnson noise dominates at elevated temperatures, we obtain where is the detector resistance, denotes the temperature and is the Boltzmann constant. Thus, we can write
Since 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 and directly from the intersubband scattering rates delivered by EMC.
A. Extraction efficiency
The extraction efficiency gives the probability that an electron which is photo-excited from the ground level to an absorption level 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.,8 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 . For this reason, we have developed an alternative method, providing a robust and compact approach to determine from the EMC scattering rates. We assume that the QCD structure is periodic with subbands per period, i.e., the transition rates from subband to fulfill . 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
Here, the transition rate matrix contains the elements for , for , and otherwise. Furthermore, vector contains the subband occupation probabilities , and the vector elements of are given by , where denotes the Kronecker delta. For incident light, we obtain a change in rate from the ground level to the absorption level due to photon-induced transitions, yielding a modified matrix where the non-zero elements of are given by , . The associated stimulated emission process is considered by setting , with , where denotes the degeneracy of a level . For intersubband transitions, holds in the case of parallel subbands, as discussed in detail in Sec. III A. The change in stationary level occupations resulting from stimulated photon transitions is found by evaluating . Re-inserting Eq. (4) and neglecting the higher order term , which is adequate for moderate intensities, i.e., if the detector operates in the linear response regime, we obtain
which yields the population changes in the irradiated structure. Here , with , , and otherwise. The occupation changes can also be written in an explicit form by applying Cramer’s rule to Eq. (5), which yields
where is the matrix obtained by replacing the th column of by vector , 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,
where Eq. (7b) has been obtained by inserting Eq. (6). As can be seen from Eq. (7b), is independent of the chosen value for , 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, can be set to a numerically convenient value. We note that for photovoltaic QCDs, the stimulated emission contribution is negligible, i.e., in Eqs. (6) and (7b), if the energetic spacing of the optical transition fulfills . 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 .15 Then, in the linear response regime assumed above, the occupation change in Eq. (7a) can be obtained by summing over the individual contributions of the absorption transitions. Here, the proportionate photon absorption rate must be used in Eqs. (5) or (6) for each transition, where gives the relative absorption contribution of level at the wavelength of interest, i.e., . This is equivalent to calculating the effective from a weighted average,
where the for each of the absorption levels has been separately computed, e.g., from Eq. (7b). Apart from , only depends on the scattering rates , which also determine the associated stationary level occupations through Eq. (4). Hence, since typically or for QCDs as mentioned above, in Eq. (8) is wavelength independent. However, this does not apply to if there are two or more relevant absorption levels, since the absorption coefficient of each transition, and thus also , is wavelength dependent.
B. Resistance
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
Here, is the sheet doping density per period. The total detector resistance in Eq. (3) can then be calculated by solving Kirchhoff’s equations for the equivalent resistance network. Assuming a network of conductances and setting , we can formulate Kirchhoff’s law for any node as
As illustrated in Fig. 1(b), the QCD is represented by a periodic network with . 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 , where 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
Here, the conductance matrix elements are given by
and the elements of 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 , the linear equations represented by Eq. (11) are not independent, which reflects the fact that there is no unique solution for since the zero potential can be arbitrarily chosen. Thus, we set and eliminate the th equation from Eq. (11), i.e., we now have , and the order of and length of and are reduced to . As for Eq. (6), Cramer’s rule can be applied to Eq. (11) to write the in explicit form, which yields
where is the matrix obtained by replacing the th column of by vector .
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,
where Eq. (13b) has been obtained by inserting Eq. (12). From Eq. (13b), we see that only depends on the , which can with Eqs. (4) and (9) be computed directly from the EMC intersubband scattering rates. In particular, is independent of the chosen value for , which approaches zero for QCDs in the photovoltaic mode. The detector resistance in Eq. (3) is given by
where is the in-plane cross section area of the quantum wells, which coincides with the detector area in Eq. (3) if the facet is in in-plane direction.
III. MODELING OF QCDS IN THERMODYNAMIC EQUILIBRIUM
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 and associated subband occupation probabilities . 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 is formed by electron states , where the quasi-free electron motion in the QW planes is described by the two-dimensional wavevector with the associated kinetic energy , and is the quantized energy in growth direction. The occupation probability of each state is given by the Fermi–Dirac distribution function
Quantum cascade devices typically utilize energy states in the conduction band valley, where the dispersion relation can for moderate kinetic electron energies be approximated by . Here, nonparabolicity effects are to first order considered by assigning a different effective mass to each subband.11 This model is schematically illustrated in Fig. 2. The electron sheet density in a subband is then obtained from Eq. (14) by multiplying with the density of states and integrating over ,11 which yields
Equation (15) only depends on , , , and the chemical potential , which can be computed from the sheet doping density per period via .11 The thermal occupation probabilities are then obtained as
Schematic illustration of the in-plane dispersion relation for ground and absorption subband and of the electron distribution function . The density of states in each subband is represented by the density of horizontal lines. The arrow indicates an optical transition.
Schematic illustration of the in-plane dispersion relation for ground and absorption subband and of the electron distribution function . The density of states in each subband is represented by the density of horizontal lines. The arrow indicates an optical transition.
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 and ,10,23,24
A. Correction factor
As mentioned above, the correction factor in Eqs. (6) and (7b) for optical transitions between two subbands and is in general not exactly unity. The total number of stimulated emission (, ) and absorption (, ) events per unit time scales as , where we use as integration variable the kinetic energy in the ground subband and thus have .24–26 Here, with , and also the dephasing rate generally depends on . can then be expressed as
For elevated temperatures and moderate doping densities, the Pauli blocking factor in Eq. (18) can be neglected, and given by Eq. (14) can be approximated by a Maxwell–Boltzmann distribution .25 In the case of parallel subbands, i.e., , Eq. (18) then yields . This is typically a good approximation for terahertz QCDs, where nonparabolicity effects tend to be negligible due to the small energetic spacing 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 in Fig. 2. accounts for two factors: stimulated photon emission from subband to 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 decreases more slowly as a function of in the upper than in the lower subband,25 as also becomes apparent from Fig. 2. Thus, the bigger , i.e., the smaller , the greater the ratio of stimulated photon emission to absorption. However, since for mid-infrared QCDs, and thus even at elevated temperatures, the term can be neglected in Eqs. (6) and (7b).
B. Ensemble Monte Carlo simulations
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 .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 of the high-lying subbands are very small, which leads to difficulties in EMC since the and corresponding scattering rates 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 and . 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 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 is imposed on the occupation probabilities, yielding after normalization modified values for all subbands. Here, 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 of the redistributed electron is randomly chosen in EMC by evaluating Eq. (14) with the inverse transform sampling method, which yields
Here, is a random variable, uniformly distributed between 0 and 1.
The (single-electron) scattering rates are evaluated in the usual way by counting the tabulated transition events from subband to over a given time interval and normalizing the result to the number of electrons in the initial level . In this way, the increased number of scattering events from high-lying subbands with artificially enhanced is automatically compensated for, yielding the correct values for . In an analogous way, the 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 are obtained since also the randomly selected scattering partner from the electron ensemble obeys the modified probability distribution . We avoid this issue by assigning a virtual scattering partner with (unmodified) thermal distribution 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 to avoid double counting does then not apply30,31).
Due to remaining fluctuations of the sampling process and other numerical limitations, the extracted rates 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 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 for , and thus tends to be determined more accurately in EMC due to the lower statistical error. Hence, for the calculation of and , the rate from the higher-energy subband to the lower-energy subband is taken from EMC, while the back-scattering rate is computed from Eq. (17) to avoid violations of the detailed balance principle resulting from numerical inaccuracies.
IV. RESULTS AND DISCUSSION
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
A. Mid-infrared quantum cascade detectors
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 and . As explained in Sec. II A, the weighted average for the extraction efficiency of the transitions and 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 . More specifically, the quantum well structure is engineered to optimize detectivity, which involves a trade-off between close to unity extraction efficiency and high detector resistance . 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.
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 () underlined is . 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.
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 () underlined is . 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.
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 W, in comparison to a measured value of 1.69 mA W. For lower temperatures, the simulated and measured responsivities show considerably better agreement. The simulated peak responsivity increases to 3.1 mA W for 100 K, as compared to a measured value of 3.8 mA W. 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.
Simulated responsivity of the QCD structure N1021 as a function of wavenumber in the temperature range of 100–300 K.
Simulated responsivity of the QCD structure N1021 as a function of wavenumber in the temperature range of 100–300 K.
To validate our detector model, we evaluate the simulation results for the extraction efficiency , power absorption efficiency and device resistance , and compare them with the experimentally measured values. Results are illustrated in Fig. 5 for 100–300 K.
Comparison of simulated (x marks) and experimental (squares) results for the detector structure N1021. (a) Extraction efficiency (solid lines) and peak absorption (dashed lines) as a function of temperature . (b) Detector resistance-area product as a function of the inverse temperature . The lines represent the Arrhenius plots with calculated activation energies .
Comparison of simulated (x marks) and experimental (squares) results for the detector structure N1021. (a) Extraction efficiency (solid lines) and peak absorption (dashed lines) as a function of temperature . (b) Detector resistance-area product as a function of the inverse temperature . The lines represent the Arrhenius plots with calculated activation energies .
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 smaller than in experiment. Due to broadening effects and reduced population of the ground level with rising temperatures, the absorption efficiency and thus the responsivity get reduced. For decreasing temperature, a shift of the probability density maximum of state to the right thinner QW and of state to the left active QW occurs due to the change in band bending associated with space charge effects. The peak absorption of the transition thus increases significantly and becomes the dominating absorption transition. The responsivity depends furthermore on the extraction efficiency , which is also displayed in Fig. 5(a). Over the plotted range, the simulated decreases linearly with increasing temperature, whereas the experimental values saturate to 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 to the ground level . 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 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 is presented in Fig. 5(b). Experimentally, the device resistance 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 meV is in good agreement with the experimentally obtained value of 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 and detector resistance-area product for the temperature range of 100–300 K are presented in Figs. 6(a) and 6(b), respectively. The simulated peak responsivity decreases from 5.3 mA W at 100 K to 1.0 mA W at 300 K similar to the experimental data. At 300 K, the values and are obtained from the simulated scattering rates, which is in good agreement with the experimentally obtained values of and at room temperature. As illustrated in Fig. 6(b), the device resistance shows the Arrhenius plot behavior. The high-temperature activation energy of meV calculated from the simulated detector resistance values compares well with the experimentally obtained value of 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
Comparison of simulated (x marks) and experimental (squares) results for the predefined figures of merit of the detector structure N1022. (a) Peak responsivity as a function of temperature . (b) Detector resistance-area product as a function of the inverse temperature . The lines represent the Arrhenius plots with calculated activation energies .
Comparison of simulated (x marks) and experimental (squares) results for the predefined figures of merit of the detector structure N1022. (a) Peak responsivity as a function of temperature . (b) Detector resistance-area product as a function of the inverse temperature . The lines represent the Arrhenius plots with calculated activation energies .
B. Terahertz quantum cascade detector
Furthermore, our approach has been applied to a GaAs/AlGaAs based terahertz QCD operating at a wavelength of around .19 Here, we assume again a 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 W at 109 cm, which compares well with the experimental value of 8.6 mA W at 119 cm. The simulated extraction efficiency of is in reasonable agreement with the experimental estimate of , which was obtained under idealized assumptions.19 For the two relevant transitions and , 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.
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 () underlined is .
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 () underlined is .
Responsivity spectrum 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.
Responsivity spectrum 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.
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/cm ) 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 cm) 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 W at 109 cm and 14.6 mA W at 129 cm from the NEGF simulations are in good agreement with the EMC results of 7.2 mA W at 110 cm and 6.7 mA W at 129 cm. 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 41 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 cm and 53.5% at 129 cm, 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.
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.
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.
From the experimentally measured dark current at ,19 we estimate a zero-bias resistance-area product at 10 K under the assumption that this bias is already in the linear regime. The simulated resistance-area products of (EMC) and (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 extracted from the measured current–voltage characteristics and obtained with EMC and NEGF is displayed logarithmically as a function of the inverse temperature .
Detector resistance-area product as a function of the inverse temperature , as obtained from EMC and NEGF simulations and extracted from experimental dark current measurements.
Detector resistance-area product as a function of the inverse temperature , as obtained from EMC and NEGF simulations and extracted from experimental dark current measurements.
A simulated specific detectivity is calculated from the EMC scattering rates at 10 K. By substituting the experimentally measured responsivity and resistance-area product into Eq. (3), we obtain a specific detectivity 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 () of the THz QCD. With the measured nW, a specific detectivity of 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.
V. CONCLUSION
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
ACKNOWLEDGMENTS
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).
DATA AVAILABILITY
The data that support the findings of this study are available from the corresponding author upon reasonable request.
AUTHOR DECLARATIONS
Conflict of Interest
The authors have no conflicts to disclose.