1 Introduction

Resonant tunnelling diodes (RTDs) (Chang et al. 1974) have been widely studied, and have been demonstrated in InGaAs/InAlAs lattices, matched to InP and GaAs/AlGaAs based systems (Mizuta and Tanoue 1995). In addition, the usage of InAs/AlSb and GaSb/AlSb systems has been recently demonstrated (Bennett et al. 2005), whilst GaN/AlGaN systems becomes a promising system (Grier et al. 2015). They provide a set of important information about the quantum mechanics of electron transport in semiconductor layered heterstructures. A short time ago, optoelectronic devices based on intersubband transitions were demonstrated in an Al-free III-V antimonide-based material system, where In\(_{0.53}\)Ga\(_{0.47}\)As was used as the well for electrons and GaAs\(_{0.51}\)Sb\(_{0.49}\) acted as the barriers (Detz et al. 2010). The understanding of the physical properties of such a system and developed experimental techniques opens the need for theoretical models and computational techniques for electron tunnelling transport in in InGaAs/GaAsSb RTDs which would be of great help in the characterization of these systems and, hence, in the optimisation of more complex heterostructure intersubband devices like quantum-cascade lasers or quantum well infrared photodetectors based on this material combination (Deutsch et al. 2010, 2012). Structures based on an Al-free material system would result in electrons with a lower effective mass in the barrier layers and, therefore, higher optical gain and improved device performance would be expected.

2 Theoretical overview

The behaviour of an electron’s wavefunction through the multilayer heterostructure is modelled by the envelope-function effective mass approximation which can be described by the common form of the Schrödinger equation

$$\begin{aligned} \left( -\frac{\hbar ^{2}}{2m_j}\frac{\partial ^{2}}{\partial z^{2}} + V_j(z)\right) \psi_j (z) = E\psi_j (z), \end{aligned}$$
(1)

where \(\hbar\) is the reduced Planck’s constant, mj is the electron effective mass, \(V_j(z)\) is the potential, and \(\psi_j(z)\) is the electron wavefunction in the \(j^{th}\) layer along the heterostructure.

The conduction band profile of the analysed InGaAs/GaAsSb double-barrier RTD under a range of applied voltages is shown in Fig. 1. If a voltage is applied across the device, the electric field F causes the conduction band to tilt. When a small voltage, U, is applied across the resonant structure (see Fig. 1b), the effective potential slopes slightly, resulting in a drop of resonant levels. However, eventually the voltage will cause the potential to slope to an extent where the Fermi energy \(E_f\) coincides with the first resonant level, \(E_1\). This means that electrons from the highly doped emitter will start to tunnel efficiently through the barrier (Fig. 1c). As the voltage continuously increases, the resonant level will drop below the conduction band minima \(E_c\), (Fig. 1d) and the tunnelling current will decrease as the tunnelling is suppressed, therefore negative differential resistivity (NDR) will be experienced. Eventually, with a further increase of the voltage, the RTD will return to a normal quasi-linear increase in tunnelling current once the next resonant level is approaching the value \(E_f\).

Fig. 1
figure 1

Left panel: The conduction band structure of a double barrier InGaAs/GaAsSb RTD under a range of applied voltages, U (electric fields F). The device consists of a quantum well sandwiched between two thin barriers. Either side of the barriers, the highly doped contact emitter/collector layers of semiconductors form a highly concentrated zone of electrons. \(E_c\) is conduction band minima and \(E_f\) is the Fermi energy. Right panel:Band structure of type-II \(\hbox {In}_{0.53}\hbox {Ga}_{0.47}\hbox {As/GaAs}_{0.51}\hbox {Sb}_{0.49}\) heterojunction. The parameters for the electronic structure are taken from Refs. Detz et al. (2010), Vurgaftman et al. (2001)

The solutions for electron wavefunctions in the analysed RTD structure from Fig. 1 are based on the solution of Eq. 1 in a more general structure which, after discretisation in a piece-wise constant manner, can be schematically presented as in Fig. 2. The solution is presented as a linear combination of travelling waves moving in opposite directions: one in the \(+\hat{z}\) direction with amplitude \(A^{+}_{j}\) and the other moving in the -\(\hat{z}\) direction with amplitude \(A^{-}_{j}\). These waves form a superposition everywhere except in the final zone N, due to the lack of a reflected wave to interfere with the transmitted component, since there are no further potential barriers

$$\begin{aligned} \psi _{j}(x) = A^{+}_{j}e^{ik_{j}z} + A^{-}_{j}e^{-ik_{j}z} \end{aligned}$$
(2)

where j is the layer index and \(k_j = \frac{\sqrt{m_j(E-V_{j}(z))}}{\hbar }\) is the wave vector. (See Fig. 2)

Fig. 2
figure 2

A demonstration of the indexing and formation of a multi-layer structure with piece-wise constant potential steps: the wavefunctions, electron effective masses and potentials in each layer. (Color figure online)

The transfer matrix method involves approximating a given potential barrier as a series of smaller piece-wise constant step layers. In the realistic model, the electron effective mass mismatch between layers must also be considered. Thus, the common Ben-Daniel Duke boundary conditions of continuities of the wavefunction and the ratio between wavefunction derivative and effective mass (Ando and Itoh 1987) are used to account for the possible effective mass mismatch at \(z= z_{j+1}\) (see red section of Fig. 2). Therefore, the transfer matrix, \(\mathbf {D}_{j}\), can be formed which describes how the potential barrier height affects the transmission of the wavefunction (Mizuta and Tanoue 1995; Ando and Itoh 1987)

$$\begin{aligned} \left( \begin{array} { c } { A _ { j } ^ { + } } \\ { A _ { j } ^ { - } } \end{array} \right) = \frac{ 1 }{ 2 } \left( \begin{array} { c c } { \left( 1 + \frac{ m _ { j } k _ { j + 1 } }{ m _ { j + 1 } k _ { j } } \right) e ^ { i \left( k _ { j + 1 } - k _ { j } \right) z _ { j + 1 } } } &{} { \left( 1 - \frac{ m _ { j } k _ { j + 1 } }{ m _ { j + 1 } k _ { j } } \right) e ^ { - i \left( k _ { j + 1 } + k _ { j } \right) z _ { j + 1 } } } \\ { \left( 1 - \frac{ m _ { j } k _ { j + 1 } }{ m _ { j + 1 } k _ { j } } \right) e ^ { i \left( k _ { j + 1 } + k _ { j } \right) z _ { j + 1 } } } &{} { \left( 1 + \frac{ m _ { j } k _ { j + 1 } }{ m _ { j + 1 } k _ { j } } \right) e ^ { - i \left( k _ { j + 1 } - k _ { j } \right) z _ { j + 1 } } } \end{array} \right) \left( \begin{array} { c } { A _ { j + 1 } ^ { + } } \\ { A _ { j + 1 } ^ { - } } \end{array} \right) \end{aligned}$$
(3)

To produce a numerical result for the value of the wavefunction at a position z along the system after N potential steps, the following closed form can be used

$$\begin{aligned} \left( \begin{array} { c } { A _ { 1 } ^ { + } } \\ { A _ { 1 } ^ { - } } \end{array} \right) = \left( \prod _ { j = 1 } ^ { N - 1 } \mathbf { D } _ { j } \right) \left( \begin{array} { c } { A _ { N } ^ { + } } \\ { A _ { N } ^ { - } } \end{array} \right) = \mathbf { Q } \left( \begin{array} { c } { A _ { N } ^ { + } } \\ { A _ { N } ^ { - } } \end{array} \right) . \end{aligned}$$
(4)

It is therefore appropriate to determine the transmission coefficient, T(E), of the wavefunction after passing through a series of potential barriers. This can be calculated using the elements of the matrix \(\mathbf {Q}\) in Eq. 4. An electron is introduced from the left hand side, and there is no electron introduced from the right, meaning that we can use \(A^{+}_{1} = 1\) and \(A^{-}_{N} = 0\) which will result in

$$\begin{aligned} T(E) = \left| \frac{A^{+}_{N}}{A^{+}_{1}}\right| ^{2} = \frac{1}{| Q_{11}|^{2}} \end{aligned}$$
(5)

The transfer matrix method is now used to solve the double barrier \(\hbox {In}_{0.53}\hbox {Ga}_{0.47}\hbox {As/GaAs}_{0.51}\hbox {Sb}_{0.49}\) heterostructure (de Sousa et al. 2011, 2010). Structural parameters are given in Fig. 1 (right panel) with band nonparabolicity effects (de Sousa et al. 2010) included, using an approach of energy dependent effective mass (Nelson et al. 1987) following the calculated nonparabolicity parameters in Ref. Demic et al. (2016).

The next step in the simulation was to apply a range of electric fields across the structure (terminal voltage), and record the current density, calculated using the Esaki-Tsu formula (Chang et al. 1974). This simplified approach is based upon a coherent picture of electron tunnelling, assuming that electron transport is not affected to any phase-coherence breaking scattering effects (Chang et al. 1974; Mizuta and Tanoue 1995). Despite its simplicity this transport model has been widely used to characterise both conventional GaAs-based (Chang et al. 1974; Mizuta and Tanoue 1995; Ando and Itoh 1987) and new material RTDs (for example GaN-based Sakr et al. 2011).

$$\begin{aligned} J = \frac{ e k _ { B } T }{ 2 \pi ^{2}\hbar ^ { 3 } } \int _ { E_c } ^ { E _ { m a x } } m ^ { * } T ( E ) \ln \left[ \frac{ 1 + \exp \left[ \left( E _ { f } - E \right) / \left( k _ { B } T \right) \right] }{ 1 + \exp \left[ \left( E _ { f } - E - e U \right) / \left( k _ { B } T \right) \right] } \right] d E, \end{aligned}$$
(6)

where for the reference level \(E_c=0\) is used, T is the crystal lattice temperature, T(E) is the transmission coefficient as mentioned before, \(E_{f}\) is the Fermi energy in the highly doped emitter (see Fig. 1), U is the potential drop across the structure (such that \(U=F \times\) length of resonant structure) and \(m^{*}\) is the non-parabolic effective mass in InGaAs. The Fermi energy \(E_f\) is calculated using Fermi-Dirac statistics in a highly doped emitter/collector assuming that all donors \(N_D\) are ionised i.e. that electron concentration in the emitter/collector is \(n=N_D\) (Blakemore 1982).

If the cross-sectional area of the diode is known, the tunnelling current can be calculated from the current density, and qualitatively compared to measurements of a realistic device. Naturally, this is qualitative due to the fact that the magnitude of the current obtained from this method does not take into account other contributing factors to the total current, such as the scattering current and the thermionic current, thus only the trends and negative differential resistivity behaviour can be identified and compared with experiment.

Furthermore, staying within the coherent tunnelling approach, we can also express the electron density in the RTD structure from Fig. 1 following Ref. Mizuta and Tanoue (1995) as

$$\begin{aligned} n(z)= \frac{ k _ { B } T }{ 2^{3/2} \pi ^{2}\hbar ^ { 3 } } \bigg ( \int _ { E_c } ^ { E _ { m a x } } m ^ { * } E^{-1/2} |\psi (z,E)|^{2} \ln \left[ 1+\exp \left( \frac{ E _ { f } - E}{k _ { B } T} \right) \right] dE \end{aligned}$$
(7)
$$\begin{aligned} + \int _ { E_c-eU } ^ { E _ { m a x } } m ^ { * } E^{-1/2} |\psi (z,E)|^{2} \ln \left[ 1+\exp \left( \frac{ E _ { f } -E-eU}{k _ { B } T} \right) \right] dE \bigg ) \end{aligned}$$
(8)

where, again, the reference level \(E_c=0\) is used.

3 Results and discussion

In the first set of calculations we analysed an InGaAs/GaAsSb double barrier RTD structure similar to one in experimental papers by de Sousa et al. (2011, 2010), which was grown lattice-matched on InP substrates. The structure has a 13 nm InGaAs quantum well sandwiched between GaAsSb barriers each of 9 nm thickness . The double barrier structure was placed between injector and collector InGaAs layers. The doping of the InGaAs injector/collector was increased in steps from \(10^{16}\) cm\(^{-3}\) to \(10^{18}\) cm\(^{-3}\) by three layers having a thickness of 100 nm each. To simplify numerical analysis, the doping of the InGaAs emitter/collector was set to be \(N_D=10^{17}\) cm\(^{-3}\) in our calculations, corresponding to an average doping of emitter/collector in the mentioned experimental studies. Thickness of each layer in digitised potential was 0.1 nm, giving the initial number of transfer matrix steps of \(N=\) 6310, which is then significantly reduced, assuming that potential was constant in the injector/collector, to \(N=\) 510.

The transmission T(E) is given in Fig. 3 for a set of different voltage drops across the barrier-well-barrier resonant tunnelling structure. The application of an external voltage across the original two barriers causes them to slope in the negative E direction as z increases, as can be seen in Fig. 1, alongside the effects of changing the potential drop on the calculated transmission coefficient, and the location of its peaks. As expected, the position of the T(E) resonances moves towards lower energies when the electric field (applied voltage) is increased. This is important to understand the results of an increasing bias/electric field on the tunnelling current density and that when the transmission coefficient’s first peak moves past the band edge in the emitter (see Fig. 1d), the current begins to decline from the peak.

Fig. 3
figure 3

Transmission coefficient for different energies in the two barrier system for three different values of voltage drop U (applied electric filed F) across resonant barrier-well-barrier InGaAs/GaAsSb structure. As an illustration the position of Fermi level in emitter doped with \(N _D=10^{17}~\rm{cm}^{-3}\) at \(T=4.2\) K is shown. (Color figure online)

In Fig. 4 current density-voltage characteristics for three different temperatures, and the Fermi energies that correspond to emitter doping of \(N_D=10^{17}\) cm\(^{-3}\), are shown. To overcome potential numerical issue with integration of a function with the very narrow and sharp peaks we used a highly accurate global adaptive quadrature method (Shampine 2008). Calculated current density at cryogenic temperature of 4.2 K shows a good qualitative agreement with measurements in Ref. de Sousa et al. (2010) owing to the fact that here, only the resonant tunnelling current is considered in our calculation. In the present model, the voltage U represents a potential drop across the barrier-well-barrier structure only, which, in reality, is only a proportion of the total voltage applied to the whole RTD structure, as significant potential decrease can be seen in the emitter and collector regions (Ohnishi et al. 1986). This is a reason that in our calculation given in Fig.4 the peak current density at \(T=4.2\) K at around 50 mV is identified, whilst the measured one was around 150 mV(de Sousa et al. 2010). Furthermore, calculated peak current of \(1 \times 10^{-4}\) A, assuming the cross section of the RTD is \(100 \mu\)m\(\times 100 \mu\)m, is somehow lower than measured \(2.8 \times 10^{-4}\) A in de Sousa et al. (2010), again because our model takes into account only the resonant tunnelling component of the total current. As expected, an onset of current density at 4.2 K at around 30 mV corresponds to an alignment of the first resonant level \(E_{1}\) with the Fermi energy in the emitter, shown also in Fig.3. Qualitatively similar results have been obtained when emitter/collector doping density was varied between \(10^{16}\) cm\(^{-3}\) and \(10^{18}\) cm\(^{-3}\) and peak current density changed between 0.05 A/cm\(^2\) and 10 A/cm\(^2\) at \(T=4.2\) K. We should note that an inclusion of the Poisson equation in the model, thus forming a complex self-consistent Schrödinger-Poisson solver (Mizuta and Tanoue 1995), would further improve agreement with the experimental findings in de Sousa et al. (2010) but with a significant CPU-time cost.

Fig. 4
figure 4

The tunnelling current density against the voltage drop across the double barrier InGaAs/GaAsSb structure with two T(E) peaks crossing the Fermi \(E_f\) and conduction band minima \(E_c\) levels as the barriers slope due to increased voltage. Calculation at the different temperatures are shown. (Color figure online)

In order to illustrate further electron tunnelling processes, we have calculated electron concentration distribution along the same RTD structure at \(T=4.2\) K, and the results are shown in Fig. 5. For zero, or very low, voltages when the energy of the first resonant level is still far from the Fermi energy, (see also Fig. 1) the tunnelling electron transport is very weak, and there are no electrons accumulated in the centre of the quantum well (dotted blue line). However, with an increase of the voltage and as the lowest resonant level draws closer to the Fermi energy, reaching an onset voltage at \(U=30\) mV, the electron accumulation in the well goes up (red dot-dashed line). A further increase of the voltage results in a significant accumulation of the electrons in the well, approaching its highest values at the voltage \(U=50\) mV that corresponds to the peak of the current. With further increase of the voltage, the electron concentration in the well starts to decrease (green dashed line) due to a less effective tunnelling process once level \(E_1\) reaches conduction band minima \(E_c\) (Fig. 1).

Fig. 5
figure 5

The electron concentration distribution across the RTD structure at a lattice temperature of \(T=\)4.2 K. Calculated curves correspond to values of the voltage in characteristic points (see Fig. 4): unbiased structure (dotes blue line), current onset point (red dot-dashed line), current peak (black solid line) and current valley point (green dashed line). Doping density of emitter/collector layers was \(N_D=10^{17}\) cm\(^{-3}\). (Color figure online)

A reason why InGaAs/GaAsSb material system did not get as much attention as other III-V heterostructures until recently is perhaps down to the fact that crystal growth is more technologically challenging. Furthermore, layer thickness variation as well as interface roughness on the order of a fraction of a monolayer is another issue that arises (Detz et al. 2010). On the other hand, modern structures based on resonant tunnelling mechanisms like terahertz quantum-cascade lasers (THz QCLs) (Deutsch et al. 2010, 2012) often require very precise growth of the layer structures in order to provide efficient electron tunnelling and a selective injection transport process. To investigate briefly an influence of layer thickness fluctuation in the In\(_{0.53}\)Ga\(_{0.47}\)As/GaAs\(_{0.51}\)Sb\(_{0.49}\) resonant tunnelling structure, we have performed another set of the current-voltage characteristics simulations in a reference double-barrier structure with nominal barrier thickness of 3 nm and well thickness of 14.4 nm, a structure similar to the resonant injection layers in InGaAs/GaAsSb THz QCL (Deutsch et al. 2010, 2012). Doping density of the emitter side was \(N_D=10^{16}\) cm\(^{-3}\) with a lattice temperature of \(T=77\) K, giving a Fermi energy of \(-6\) meV with respect to the conduction band minima. Obtained results are shown in Fig. 6, notably the black solid line. Following that, a set of calculations with the barrier and well thickness fluctuation for one monolayer (\(\approx 0.6\) nm) (Detz et al. 2010) around nominal values was performed and results are also shown in Fig. 6. It has been noted that InGaAs quantum well thickness fluctuation has produced relatively mild changes in the magnitude of the current peak but with some impact on the position of the valley voltage (black dotted and dashed line). However, at the same time, monolayer fluctuation of GaAsSb barrier thickness, despite being Al-free, with a low electron effective mass, would produce a significant impact on the magnitude of the tunnelling current (blue and red solid line), being an important parameter in prospective InGaAs/GaAsSb THz QCL electron transport optimisation. It is important to recall that only qualitative analysis of the magnitude of the current and comparison of the trends and behaviour is possible here due to the fact that this model only takes into account the tunnelling current, whereas in real experiments, other contributing factors like electron scatterings and thermionic emission are also present.

Fig. 6
figure 6

The influence of barrier/well layer thickness fluctuation on the tunnelling current density as a function of the voltage drop across the double barrier InGaAs/GaAsSb structure, calculated at the temperature of \(T=77\) K. Doping density of the emitter side is \(N_D=10^{16}\)  cm\(^{-3}\). (Color figure online)

4 Conclusion

Experimental realisations of aluminium-free intersubband structures based on In\(_{0.53}\)Ga\(_{0.47}\)As/GaAs\(_{0.51}\)Sb\(_{0.49}\) system lattice matched to InP, opens a need for a simple theoretical model and analysis of the coherent tunnelling transport throughout the structure. We have modelled current-voltage characteristics at different lattice temperatures and analysed electron density distribution as a function of the voltage applied to the double-barrier resonant tunnelling structure. A good qualitative agreement with available experimental results are obtained. Calculations show that tunnelling current peak is very sensitive to small (monolayer) thickness variation, and, for example, could change almost by a factor of two when the thickness of both barriers vary from their target thickness of 3 nm by +/- one monolayer (2 monolayers total variation). This affects peak to valley ratio in a resonant tunnelling structure which may be useful information for optimising resonant tunnelling transport and injection efficiency in more complex structures like InGaAs/GaAsSb QCLs.