Self-excited current oscillations in a resonant tunneling diode described by a model based on the Caldeira–Leggett Hamiltonian

Quantum coherence and its destruction by coupling to a dissipative environment play an important role in the transport phenomena of a particle moving in a potential [1–3]. Well known examples include electron transfer in molecular and biological systems [4, 5], many chemical reactions [6–8], SQUID rings [9, 10], quantum ratchets [11, 12], nonlinear optical processes [13–17] and tunneling processes in device systems [18, 19]. Such systems are commonly modeled as one-dimensional or two-dimensional potential systems coupled to heat bath degrees of freedom, which drive the systems toward the thermal equilibrium state. The heat bath degrees of freedom are then reduced using such methods as the projection operator method or the path integral method, for example. Many equations of motion have been derived for the purpose of understanding the quantum aspects of dissipative dynamics [20–30].

Because a complete picture of quantum dissipative dynamics must treat phenomena that can only be described in real time, a great deal of effort has been dedicated to the problem of numerically integrating these equations of motion in real time [25–34]. Although these equations are analogous to the classical kinetic equations, which have proved to be useful for classical transport problems, such equations cannot be derived in a quantum mechanical framework without significant approximations and/or assumptions. For example, the quantum Boltzmann equation is based on the assumption that the effects of collisions between electrons can be described by the rates determined from Fermi's golden rule, and hence it is regarded as a semi-classical equation [20–22]. Similarly, the quantum Fokker–Planck equation can be derived from the Caldeira–Leggett Hamiltonian under a Markovian approximation, but in order for this to be possible, the heat bath must be at a sufficiently high temperature, in which case most of the important quantum dynamical effects play a minor role [23, 24]. Treatments of these kinds are therefore not sufficient to construct fully quantum mechanical descriptions of broad validity.

To circumvent this problem, we present a quantum mechanical approach, which is valid for arbitrary temperatures. This treatment employs the reduced hierarchy equations of motion (HEOM), and it can be used in application to systems for which fully quantum mechanical description is necessary [35, 36]. In particular, the reduced HEOM approach can be used to numerically treat non-Markovian system–bath coupling in a non-perturbative manner [37–46]. Here, we apply this approach to study the dynamics of a resonant tunneling diode (RTD) described by the Caldeira–Leggett Hamiltonian.

The RTD system that we consider is modeled by a double barrier structure with an electrostatic potential representing a region consisting of an undoped layer positioned between two doped layers (see figure 1). The double barrier structure constitutes a single quantum well with discretized energy states for electrons which are called resonant levels. The number of resonant levels depends on the height of the barriers. When a bias is applied to the RTD system, as long as the energy of electrons which flow in the RTD is lower than the energy of the resonant level, the most of electrons are reflected by the barrier because of the small transmission coefficient. When the energy of the electrons matches the resonant energy, electrons can go through the barriers efficiently due to resonant tunneling, and the current acquires the maximum value. On the other hand, the current decreases after the energy of the electrons exceeds the resonant energy. As a result, RTD systems exhibit characteristic negative differential resistance (NDR) in the current–voltage (I–V ) relation [47]. Until now, RTDs have been mostly used as high-frequency oscillators device using NDR characteristics because the tunneling is the fastest charge-transport mechanism in semiconductors [48]. RTDs are presently the highest-frequency active semiconductor devices in existence [49–51] and oscillation frequencies above 1 THz have recently been realized [52–54].

Zoom In Zoom Out Reset image size

From a fundamental physics point-of-view, the RTD system provides a simple and convenient 'context' for studying and testing various methods of analysis for nanoscale quantum devices [18, 19]. Frensley discovered NDR in the I–V curve through a numerical computation treating a quantum Liouville equation in the Wigner representation that adopted open boundary condition and ignored phonon-scattering processes [55–57]. Kluksdahl et al [58] incorporated dissipative and self-consistent effects, employing the Poisson–Boltzmann equation by adopting a relaxation time approximation, and succeeded in modeling the experimentally observed hysteresis behavior of the I–V curve. Jensen and Buot developed a numerical scheme to treat systems of the same kind and found evidence that the current oscillation and plateau-like behavior arise from intrinsic bistability [59–64].

When plateau-like behavior and hysteresis of the I–V curve in the NDR region, which were thought to arise from the feedback of the electrostatic field, were experimentally observed [65, 66], Sollner claimed that they result merely from resonance with the external circuit [67]. Although theoretical calculations have provided evidence of intrinsic bistability and self-excited current oscillations in the NDR region, such phenomena have not been justified by experimental means. In addition, because there exists no well-established methodology that can be applied rigorously to this type of model and includes the effect of dissipation, which is the origin of Joule heat, previous theoretical results have not been well justified. The HEOM approach is ideal to clarify a role of bistability in the NDR region and for detailed analysis of the RTD system.

This paper is organized as follows. In section 2, we introduce the reduced HEOM applicable to the resonant tunneling problem. We then present the computational details for the numerical simulations in section 3. Numerical results for the I–V curves and current oscillations are presented in section 4. Section 5 is devoted to concluding remarks.

We consider the following Caldeira–Leggett Hamiltonian [1], which describes the dynamics of an electron subjected to a thermal environment:

$$\begin{equation} \skew3\hat{H} = \frac{\skew3\hat{p}^2}{2m} + U(\skew3\hat{q};t) + \sum_j \left[ \frac{\skew3\hat{p}_j^{2}}{2m_j} + \frac{m_j\omega_j^2}{2} \left( \skew3\hat{x}_j - \frac{a_j V(\skew3\hat{q})}{m_j\omega_j} \right)^2 \right]. \end{equation} \tag{ 1 }$$

Here, m, $\hat {p}$ and $\hat {q}$ are the mass, momentum and position variables of the electron, and $m_j, \hat {p}_j, \hat {x}_j$ and ω_j are the mass, momentum, position and frequency variables of the jth phonon oscillator mode. In equation (1), the electron–phonon interaction is given by

$$\begin{equation} \skew3\hat{H}_{\rm I} = - V(\skew3\hat{q}) \sum_j a_j \skew3\hat{x}_j . \end{equation} \tag{ 2 }$$

Here, $V(\hat {q})$ is any function of $\hat {q}$ and the quantities a_j are coefficients that depend on the nature of the electron–phonon coupling.

The heat bath can be characterized by the spectral distribution function, defined by

$$\begin{equation} J(\omega) \equiv \sum_j \frac{a_j^2}{2m_j\omega_j} \delta(\omega-\omega_j) \end{equation} \tag{ 3 }$$

and the inverse temperature, β ≡ 1/k_BT, where k_B is the Boltzman constant. We assume the Drude distribution, given by

$$\begin{equation} J(\omega) = \frac{m\zeta}{\pi} \frac{\gamma^2\omega}{\gamma^2 + \omega^2}, \end{equation} \tag{ 4 }$$

where the constant γ represents the width of the spectral distribution of the collective phonon modes and is the reciprocal of the correlation time of the noise induced by phonons. The parameter ζ is related to the electron–phonon coupling strength. For the collective heat bath coordinate $\hat {X} = \sum _j a_j \hat {x}_j$, the canonical and symmetrized correlation functions, respectively defined by $\Psi (t) \equiv \beta \langle \hat {X} ; \hat {X}(t) \rangle _{\mathrm {B}}$ and $C(t) \equiv \frac {1}{2} \langle \hat {X}(t) \hat {X}(0) + \hat {X}(0) \hat {X}(t) \rangle _{\mathrm {B}}$, where $\hat {X}(t)$ is the Heisenberg representation of $\hat {X}$, and 〈··· 〉_B represents the thermal average over the bath modes, are given by [35, 36]

$$\begin{equation} \Psi(t) = { m \zeta \gamma}\, {\rm e}^{ - \gamma \left| t \right|} \end{equation} \tag{ 5 }$$

and

$$\begin{equation} C(t) = c_0 \mathrm{e} ^{ - \gamma \left| t \right|} + \sum\limits_{k = 1}^\infty {c_k } \mathrm{e} ^{ - \nu _k \left| t \right|}. \end{equation} \tag{ 6 }$$

Here, ν_k ≡ 2πk/βℏ are the Matsubara frequencies, and we have

$$\begin{equation} c_0 = \frac{\hbar m\zeta \gamma ^2 }{2}\left[ {\frac{2}{{\beta \hbar \gamma }} + \sum_{k=1}^\infty{\frac{{4\beta \hbar \gamma }}{{(\beta \hbar \gamma )^2 - (2\pi\! k)^2 }}} } \right] \end{equation} \tag{ 7 }$$

and

$$\begin{equation} c_k = - \frac {\hbar m \zeta \gamma ^2 } {2} \frac{{ 8 \pi k}} {{(\beta \hbar \gamma )^2 - (2\pi\! k)^2 }}. \end{equation} \tag{ 8 }$$

The function C(t) is analogous to the classical correlation function of X(t) and corresponds to the correlation function of the bath-induced noise, whereas Ψ(t) corresponds to dissipation. The noise C(t) is related to Ψ(t) through the quantum version of the fluctuation–dissipation theorem, C[ω] = ℏω cot h(βℏω/2)/2Ψ[ω], which ensures that the system exists in the thermal equilibrium state for finite temperatures [68]. Note that in the high temperature limit, βℏγ ≪ 1, the noise correlation function reduces to $C(t) \propto {\mathrm { e}}^{ - \gamma \left | t \right |}$. This indicates that the heat bath oscillators interact with the system in the form of Gaussian–Markovian noise.

To derive the equation of motion for the electron, we use the reduced density operator of the system by taking the trace over the heat bath degrees of freedom:

$$\begin{equation} \skew3\hat{\rho}(t) = \mathrm{Tr}_{\mathrm{B}}\, {\skew3\hat{\rho}_{\mathrm{tot}}(t)}. \end{equation} \tag{ 9 }$$

In the path integral representation, the reduced density matrix elements are written

$$\begin{eqnarray} &&\fl \rho (q,\,q';\,t)= \int {D[q( \tau )]\int {D[q'( \tau )]\int {\mathrm{d}q_i } \int {\mathrm{d}q'_i } } } \rho (q_i ,q'_i )\rho_{\mathrm{CS}} (q,\,q_i,\,q',\,q'_i;\,t)\nonumber\\ &&\times {\rm e}^{\frac{\rm i}{\hbar }S_{\mathrm{A}} [q;\,t]}F[q,\,q';\,t]{\rm e}^{ - \frac{\rm i} {\hbar }S_{\mathrm{A}} [q';\,t]},\end{eqnarray} \tag{ 10 }$$

where S_A[q; t] is the action for the Hamiltonian of the system, $\hat {H}_{\mathrm {A}}=p^2/2m+U(q; t)$, expressed as

$$\begin{equation} S_{\mathrm{A}} [q;\,t] \equiv \int_{t_i}^{t} \mathrm{d}\tau \left[ \frac{1}{2} m \dot q^2(\tau) - U(q (\tau); \tau) \right], \end{equation} \tag{ 11 }$$

ρ(q_i,q'_i) is the initial state of the system at time t_i, F[q, q'; t] is the influence functional [2] and ρ_CS(q, q_i, q', q'_i; t) is the initial correlation function between the system and the heat bath [3]. The functional integrals for q(τ) and q'(τ) are carried out from q(t_i) = q_i to q(t) = q and from q'(t_i) = q_i' to q'(t) = q', respectively. In the HEOM approach, we can specify ρ_CS(q, q_i, q', q'_i; t) by non-zero hierarchy elements. To simplify the derivation of the HEOM, here we set ρ_CS(q, q_i, q', q'_i; t) = 1 and regard ρ(q_i, q'_i) as a temporal initial condition. Then, after deriving the HEOM, we take into account ρ_CS(q, q_i, q', q'_i; t) through implementation of a hierarchy of initial conditions that can be evaluated numerically [35]. The influence functional for the inverse temperature β is given by [2, 3]

$$\begin{eqnarray} &&\fl F[q,\,q';\,t] ={\rm exp}\left\{ \left( { - \frac{\rm i}{\hbar }} \right)^2 \int_{t_i }^t {\mathrm{d}\tau } V^{\times}(q,\,q';\,\tau)\left[\int_{t_i }^{\tau } {\mathrm{d}\tau'} \,\frac{\mathrm{i}\hbar}{2} \Psi (\tau - \tau')\frac{\partial}{\partial\tau'}V^{\circ}(q,\,q';\,\tau')\right. \right.\nonumber\\ &&\left. \left. + \int_{t_i }^{\tau } {\mathrm{d}\tau' } C(\tau - \tau' ) V^{\times}(q,\,q';\,\tau') \vphantom{ \left( { - \frac{\rm i}{\hbar }} \right)^2}\right] \right\}, \end{eqnarray} \tag{ 12 }$$

where $V^{\times }(q,\,q';\,\tau ) \equiv {V({q(\tau )}) - V( {q'(\tau )})}$ and $V^{\circ }(q,\,q';\,\tau ) \equiv {V({q(\tau )}) + V ( {q'(\tau )})}$. If we choose K so as to satisfy ν_K = 2πK/(βℏ) ≫ ω_c, where ω_c is the characteristic frequency of the system such as the frequency of self-excited current oscillations, the factor e^−ν_k|t| in equation (6) can be replaced with Dirac's delta function, using the approximation ν_k e^−ν_k|t| ≃ δ(t) (for k ⩾ K + 1). Therefore, C(t) can be expressed as

$$\begin{equation} C(t) = c_0 \mathrm{e} ^{ - \gamma |t|} + \sum\limits_{k = 1}^{K} {c_k } \mathrm{e} ^{ - \nu _k |t|} + \delta(t) \sum\limits_{k = K+1}^{\infty} \frac{c_k}{\nu_k} . \end{equation} \tag{ 13 }$$

By choosing 2πK ≫ βℏω_c, the above expression allows us to evaluate C(t) for finite K with negligible error at the desired temperature 1/β.

The reduced HEOM can be obtained by considering the time derivative of the reduced density matrix with the kernel given in equations (5) and (13). The HEOM have been used to study chemical reactions [37, 38, 69, 70], linear and nonlinear spectroscopy [39–45, 71–73], exciton transfer [74–79], electron transfer [80–83], quantum dots [73, 84], quantum ratchet [46] and quantum information [85–87]. A variety of numerical techniques have been developed for the HEOM approach in order to accelerate numerical calculations [88–95]. The accuracy of the HEOM approach has been justified for a Brownian oscillator system [41–45], a displaced Brownian oscillators system [40] and a spin-boson system [39, 71, 86] via linear and nonlinear response functions by comparing the analytical solutions of the response functions [15–17]. The validity of the HEOM are also confirmed with other numerically methods such as the iterative quasi-adiabatic propagator path-integral scheme and a time-convolution less master equation in the relevant crossover regime from weak to strong system–bath coupling [96–98].

The HEOM are ideal for studying quantum transport systems, in conjunction with the Wigner representation, characterized by the Wigner distribution function

$$\begin{equation} W(p,q;t) \equiv \int_{-\infty}^\infty \mathrm{d}r\,\mathrm{e}^{-\frac{\mathrm{i} pr}{\hbar}} \rho\left(q+\frac{r}{2},q-\frac{r}{2};t \right), \end{equation} \tag{ 14 }$$

because they allow us to treat continuous systems utilizing open boundary conditions and periodic boundary conditions [57, 99]. Although the Wigner distribution function is not positive definite, it is the quantum analogue of the classical distribution function in the phase space [55–64, 99–101]. Its classical limit can be computed readily. This is helpful, because knowing the classical limit allows us to identify the purely quantum mechanical effects [37, 38, 45, 46, 64].

While we can handle any form of V (q) [41–45], here we consider the linear–linear system–bath coupling case defined by V (q) = q. In the Wigner representation, the equations of motion are expressed in hierarchical form as follows [35, 45]:

$$\begin{eqnarray} &&\fl \frac{\partial}{\partial t} W^{(n)}_{j_1, \ldots, j_K}(t)= - \left[ \skew3\hat{L}_{\mathrm{qm}} + \skew3\hat{\Xi}' + n\gamma + \sum_{k=1}^K j_k \nu_k \right]W^{(n)}_{j_1, \ldots, j_K}(t) + \hat{\Phi}\left[\vphantom{\sum_{k=1}^K W^{(n)}_{j_1, \ldots, (j_k+1), \ldots j_K}}W^{(n+1)}_{j_1, \ldots, j_K}(t) \right. \nonumber \\ &&\left.+ \sum_{k=1}^K W^{(n)}_{j_1, \ldots, (j_k+1), \ldots j_K}(t) \right]+ n\gamma\hat{\Theta}_0 W^{(n-1)}_{j_1, \ldots, j_K}(t) + \sum_{k=1}^K j_k \nu_k \hat{\Theta}_k W^{(n)}_{j_1, \ldots, (j_k-1), \ldots, j_K}(t)\nonumber\\ \end{eqnarray} \tag{ 15 }$$

for non-negative integers n,j₁,...,j_K, where we have chosen K such that ν_K ≫ ω_c. In equation (15), $-\hat {L}_{\mathrm {qm}}$ is the quantum Liouvillian in the Wigner representation, given by

$$\begin{equation} \fl-\skew3\hat{L}_{\mathrm{qm}} W(p,q) \equiv -\frac{p}{m} \frac{\partial}{\partial q}W(p,q)- \frac{1}{\hbar}\int_{-\infty}^\infty \frac{\mathrm{d}p'}{2\pi\hbar} U_{\mathrm{W}}(p-p',q; t) W(p',q) \end{equation} \tag{ 16 }$$

with

$$\begin{equation} U_{\mathrm{W}}(p,q;t) \equiv 2\int_0^\infty \mathrm{d}r \sin\left( \frac{pr}{\hbar} \right) \left[ U\left(q+\frac{r}{2}; t\right) - U\left(q-\frac{r}{2}; t \right) \right]. \end{equation} \tag{ 17 }$$

The other operators appearing in equation (15) are the bath-induced relaxation operators, defined as

$$\begin{equation} \hat{\Phi} \equiv \frac{\partial}{\partial p}, \end{equation} \tag{ 18 }$$

$$\begin{equation} \hat{\Theta}_0 \equiv \zeta \left[p + \frac{m\hbar\gamma}{2} \cot \left( \frac{\beta\hbar\gamma}{2}\right) \frac{\partial}{\partial p}\right], \end{equation} \tag{ 19 }$$

$$\begin{equation} \hat{\Theta}_k \equiv \frac{c_k}{\nu_k} \frac{\partial}{\partial p} \end{equation} \tag{ 20 }$$

and

$$\begin{equation} \skew3\hat{\Xi}' \equiv \left\{-\frac{m\zeta}{\beta} \left[1 - \frac{\beta\hbar\gamma}{2}\cot\left( \frac{\beta\hbar\gamma}{2} \right) \right]+ \sum_{k=1}^K \frac{c_k}{\nu_k} \right\} \frac{\partial^2 }{\partial p^2}. \end{equation} \tag{ 21 }$$

In the case that the quantity $N \equiv n + \sum _{k=1}^K j_k$ satisfies the relation $N \gg {\omega _{\mathrm {c}}}/{\min (\gamma ,1/\beta \hbar )}$, this infinite hierarchy can be truncated with negligible error at the desired temperature 1/β by the terminator [45]

$$\begin{equation} \frac{\partial}{\partial t} W^{(n)}_{j_1, \ldots, j_K}(t) = - (\skew3\hat{L}_{\mathrm{qm}} + \skew3\hat{\Xi}') W^{(n)}_{j_1, \ldots, j_K}(t). \end{equation} \tag{ 22 }$$

The validity of the above truncation scheme and its extension for efficient numerical calculations have been discussed for the spin-Boson system [88–90]. Note that only W⁽⁰⁾_0,...,0(p,q;t) ≡ W(p,q;t) has physical meaning, and the other elements $W^{(n)}_{j_1, \ldots , j_K}(p,q;t)$ with $(n;j_1, \ldots , j_{{K}}) \ne (0;0,\ldots ,0)$ are auxiliary operators introduced to avoid the explicit treatment of the inherent memory effects that arise in the time evolution of the reduced density matrix. If the noise correlation is very short (γ → ∞) and the temperature is high (i.e. the noise is white), the quantum Fokker–Planck equation can be derived in a form similar to that of Kramers equation [23–25]. In the present case, however, we cannot employ the white noise approximation, because quantum effects play a dominant role in the low temperature regime (βℏω_c ≫ 1) [39, 46].

In equation (17), $U(\hat {q};t)$ is the effective potential for the electron, which can be written [58–63]

$$\begin{equation} U(\skew3\hat{q};t) = U_{\mathrm{static}}(\skew3\hat{q}) + U_{\mathrm{self}}(\skew3\hat{q};t), \end{equation} \tag{ 23 }$$

where $U_{\mathrm {static}}(\hat {q})$ and $U_{\mathrm {self}}(\hat {q};t)$ are the static and self-consistent parts, respectively. As the static potential, we employ the double-barrier structure depicted in figure 1(b). The self-consistent part, U_self(q;t) = −eϕ(q;t), is calculated from the electron distribution at each time step in the integration of equation (15) using the Poisson equation

$$\begin{equation} - \frac{\partial}{\partial q} \left[\epsilon \phi(q;t) \right] = e \left[n^{+}(q) - P(q,t) \right], \end{equation} \tag{ 24 }$$

where is the dielectric constant, n⁺(q) is the doping density and $P(q;t) = \int _{-\infty }^{\infty }\mathrm {d}p\,W^{(0)}_{0, \ldots ,0}(p,q;t) /(2\pi \hbar )$ is the electron density calculated from the Wigner distribution. Coupling the HEOM to the Poisson equation, we obtain a fully self-consistent model of quantum electron transport. This allows us to examine charge redistribution effects.

The equations of motion given in (15) were numerically evaluated using finite mesh representations of the Wigner distribution functions. The spatial derivative of the kinetic term in the Liouville operator, −(p/m)∂W(p,q)/∂q, was approximated by using a third-order left-handed or right-handed difference scheme. Note that the first-order difference scheme is not sufficiently accurate for the present problem, because this scheme introduces false diffusion for the wavepacket dynamics, which suppresses the self-excited current oscillations. Depending on the sign of the momentum, the expressions are given by

$$\begin{equation} \fl \frac{\partial W(p_k,q_j)}{\partial q} =\frac{1}{6\Delta q} ( {2}W(p_k,q_{j+1}) + 3W(p_k,q_j) - 6W(p_k,q_{j-1}) + W(p_k,q_{j-2})) \end{equation} \tag{ 25 }$$

for p_k > 0, and

$$\begin{equation} \fl \frac{\partial W(p_k,q_j)}{\partial q} =\frac{1}{6\Delta q} (- W(p_k,q_{j+2}) + 6W(p_k,q_{j+1}) - 3W(p_k,q_j) {-} 2W(p_k,q_{j-1}) ) \end{equation} \tag{ 26 }$$

for p_k < 0, in order to treat continuous systems utilizing the inflow and outflow boundary conditions with use of the first-order difference scheme at q = L (p > 0) and q = 0 (p < 0) [55–57].

The inflow boundary conditions were set by stipulating that W⁽ⁿ⁾_j1,...,jK(p < 0,q = L) and W⁽ⁿ⁾_j1,...,jK(p > 0, q = 0) are given by the equilibrium distribution of a free particle calculated from the HEOM with periodic boundary conditions. Due to fluctuations and dissipation, the flow of a wavepacket reaches a steady state even when there exists a non-zero bias voltage. The validity of the boundary conditions was verified by considering several system sizes.

Other derivatives with respect to p were approximated using the fourth-order centered difference scheme given by

$$\begin{eqnarray} &&\fl \frac{\partial W(p_k,q_j)}{\partial p}= \frac{1}{12\Delta p}( -W(p_{k+2},q_j) + 8W(p_{k+1},q_j) - 8W(p_{k-1},q_j) + W(p_{k-2},q_j) ) \end{eqnarray} \tag{ 27 }$$

and

$$\begin{eqnarray} &&\fl \frac{\partial^2 W(p_k,q_j)}{\partial p^2} = \frac{1}{12\Delta p^2} (- W(p_{k+2},q_j) + 16W(p_{k+1},q_j)- 30W(p_k,q_j) \nonumber \\ &&+ 16W(p_{k-1},q_j) - W(p_{k-2},q_j) ) . \end{eqnarray} \tag{ 28 }$$

The mesh size for the position, Δq, and momentum, Δp/(2πℏ), are respectively 0.2825 nm and 3.540 nm⁻¹.

We set the parameters used in the HEOM as γ = 24.2 THz (γ⁻¹ = 4.13 fs), ζ = 72.5 GHz (ζ⁻¹ = 13.8 ps) and T = 300 K in order to create conditions close to those used in previous theoretical studies. In appendix B, we also report the results of calculations for the strong coupling case, with ζ = 120.8 GHz (ζ⁻¹ = 8.28 ps), to elucidate the role of dissipation. The depth of the hierarchy and the number of Matsubara frequencies were chosen as N∈{2 − 6} and K∈{1 − 3}, respectively. To model GaAs, the effective mass of the electron was assumed to be constant across the device and equal to 0.067m₀, where m₀ is the electron mass in vacuum. The dielectric constant in equation (24) was set as  = 12.85.

As the static double-barrier potential, which models the hetero-structure of GaAs sandwiched between two thin AlGaAs layers, we set the widths of quantum well (undoped GaAs), barrier (undoped AlGaAs) and spacer layer (undoped GaAs) to be 4.520, 2.825 and 2.825 nm, respectively. The height of the potential barriers was 0.27 eV. The conduction band edge consists of a single quantum well bounded by tunneling barriers (figure 1(b)). The widths of the contact regions (the yellow parts in figure 1(a), where GaAs is doped with a concentration of 2 × 10¹⁸ cm⁻³) were chosen as 16.950, 33.900 and 42.375 nm. Note that, to adapt the one dimensional model, we rescaled the concentration by multiplying the doping density by unit area.

In a previous study [102], we chose a smaller value of γ, γ = 12.1 THz (γ⁻¹ = 8.26 fs), and fixed the width of the contact region as 42.375 nm. Note that, since the effective system–bath coupling strength is estimated as ∝ζγ²ω_c/(γ² + ω²_c), where ω_c is the characteristic frequency of the system [38], the damping strength in the present case is slightly larger than the previous case. In that case, we found hysteresis, double plateau-like behavior, and self-excited current oscillation in the NDR region of the current–voltage curve. We found that while most of the current oscillations decay in time in the NDR region, there exists a non-transient oscillation characterized by a tornado-like rotation in the Wigner space in the upper plateau of the NDR region. In this paper, we explore the cause of such current oscillations by considering several values of the width of the contact regions in cases of both weak and strong coupling, characterized by different values of ζ.

We determined the characteristics of the current–voltage (I–V ) according to the following procedure. First, we integrated equation (15) at zero bias voltage without the self-consistent part of the effective potential under the inflow boundary conditions specified above. When we obtained the temporal steady state, the obtained distribution was then used as the initial distribution for the self-consistent calculation, and we then integrated equation (15) again, with the effective potential $U(\hat {q};t)$ evaluated iteratively using the Poisson equation given in (24). Under this procedure, when the distributions reached the genuine steady state $W^{(n)}_{j_1, \ldots , j_K}(p,q; t\rightarrow \infty )$, the current was calculated by $I(t)=\int \mathrm {d}p\ pW_{0,\ldots ,0}^{(0)}(p,q;t)/2\pi \hbar m$ and then the genuine state was used as the initial distributions for the next bias step. While the temporal steady states were obtained for the static potential, $U_{\mathrm {static}}(\hat {q})$, the genuine steady states were calculated from the effective potential, $U_{\mathrm {static}}(\hat {q}) + U_{\mathrm {self}}(\hat {q};t)$. Since the value of the effective potential depended on the hysteresis of a physical process and since we wanted to use a uniquely determined steady state, we chose the temporal steady state as a temporal initial state to have the genuine steady state.

Following the above steps, we increased the bias from 0.000 to 0.500 V, and then decreased it to 0.000 V with bias steps of 0.01 V in the normal region and 0.002 V in the NDR region. The corresponding sweeping rates were 5 × 10⁹ V s⁻¹ in the normal region and 5 × 10⁷ V s⁻¹ in the NDR region. We found that the profiles of the I–V curves did not change for slower sweeping rate than the present values, whereas the width of plateau observed in the NDR region often became smaller for a faster sweeping rate. At each step, we integrated the equation of motion until the system exhibited the steady current. However, in some cases, in the NDR region, steady current oscillation arose. In such cases, the value of current in figure 2(a) was evaluated as a time average after stable oscillations were realized (between 30 and 40 ps in most cases).

Zoom In Zoom Out Reset image size

In figure 2, we present (a) current–voltage (I–V ) relations and (b) the time evolution of the self-excited current oscillation in the weak coupling case (ζ = 72.5 GHz) for three values of the width of the contact regions (the yellow parts in figure 1(a)). In these plots, the sizes of the contact regions are (i) 16.950 nm, (ii) 33.900 nm and (iii) 42.375 nm, respectively, with a fixed doping concentration of GaAs. These graphs reveal NDR behavior, hysteresis and plateaus in the I–V curve. Moreover, self-excited current oscillation appears in some regions of the plateau. While Jensen and Buot [59] observed only a single plateau similar to that in figure 2(i-a) with current oscillation, our results in figures 2(ii-a) and (iii-a) exhibit a double plateau structure. While the experimental result shown by Goldman et al [65] is similar to figure 2(i-a), those shown by Slight et al [49], Asada et al [52] and Suzuki et al [53] are similar to figures 2(ii-a) and (iii-a). For the sake of comparison, we present a graph corresponding to figure 2 calculated from the Boltzmann equation in appendix C. In the case of a single plateau structure (figure 2(i-a)), the width of the plateau is large, while in the case of a double plateau structure (figures 2(ii-a) and (iii-a)), the width is small. The I–V curves obtained in the strong coupling case (ζ = 120.8 GHz) are presented in appendix B. In that case, NDR behavior, hysteresis and a single plateau in the I–V curve are observed, but steady current oscillation does not appear, even in the plateau.

As in the case considered in the previous paper [102], we find that most of the current oscillations decay in time in the NDR region, but there also exist non-decaying oscillations in some regions of the plateau, as seen in figures 2(i-b) to (iii-b). The Fourier components of each persistent oscillation are plotted in figures 3(i-a) to (iii'-a). We can classify these oscillations into two types, according to the current amplitude. The first type is observed in the single plateau (figure 2(i)) and the lower part of the double-plateau (figure 2(iii)) with large amplitude (red). The plateau in this case is located in the middle of the NDR region. The second type is observed in the upper part of the double-plateau (figures 2(ii) and (iii)) with small amplitude (green). The plateau in this case is located at a current approximately three-quarters of the peak current. The second type of oscillation contains two Fourier components (figures 3(ii-a) and (iii-a)), whereas the first type contains just a single component (figures 3(i-a) and (iii'-a)). We find that as the size of the contact regions increases, the frequency of each peak decreases.

Zoom In Zoom Out Reset image size

As mentioned by Kluksdahl et al [58] and Zhao et al [62], the quantum well formed on the emitter side of the effective potential plays an important role in the realization of hysteresis and plateau-like structure in the NDR region. The time averaged electron densities (black dashed curves) and the effective potentials (black solid curves) calculated from the Wigner distribution function in the case of increasing bias are plotted in figures 3(i-b) to (iii'-b). For reference, in appendix A, we present a graph corresponding to figure 3(iii'-b) depicting the situation in the case of decreasing bias. A basin-like potential on the emitter side (emitter basin) is observed in the case of increasing bias, while the emitter basin does not exist in the case of decreasing bias. In the case of increasing bias, when the bias exceeds the peak point of the I–V curve, the kinetic energy of the inflowing electron becomes larger than the eigenenergy of the resonant tunneling state. As a result, the reflection of the current from the emitter side of the barrier becomes large, and the electron distribution function becomes concentrated near the barrier, as depicted in figures 3(i-b) to (iii'-b). When the electron density increases, the effective potential decreases. As a result, the emitter basin appears. When the emitter basin becomes sufficiently deeper, there appear resonant tunneling states between the emitter basin and the double-barrier well. In such a case, we find current oscillation and a plateau of the I–V curve in the NDR region. In the case of decreasing bias, however, because there is no resonant tunneling state between the emitter basin and double-barrier well, the current is much smaller than in the case of the increasing bias. This difference causes the hysteresis behavior.

To elucidate this point more clearly, we solve the steady-state Schrödinger equation for the regions of emitter basin and the double-barrier well, to obtain approximate eigenstates and eigenenergies of an electron whose energy is lower than the continuous energy band on the emitter side. Since bound states are not formed in the collector side of the potential, we exclude this region from the calculations. It should be noted that we employ the time-averaged potential for the purpose of the graph, but the effective potential and the corresponding eigenstates actually vary in time, because they depend on the electron distribution function, and it varies in time. Thus, for example, the identifications of the first (red) and second (green) eigenstates and the third (blue) and fourth (orange) eigenstates change in time, often becoming degenerate and interchanging. In addition, we ignore the continuous band on the collector side and the influence of the heat bath when calculating the eigenstates and eigenenergies. For this reason, the calculated eigenenergies are not precise, but we find that they are sufficient for determining the cause of the current oscillation, because each resonant frequency is rather isolated when estimating the oscillation frequency.

We find that each oscillation peak in both the large and small oscillation cases can be attributed to transitions between eigenstates in the emitter basin, as depicted in the insets of figures 3(i-a) to (iii'-a). As shown in figures 3(i-b) and (iii'-b), the first and second eigenstates are the tunneling states in the large oscillation case, while the third eigenstate in figures 3(ii-b) and (iii-b) is the tunneling state in the small oscillation case. When we compare the profiles of each eigenstate and the electron density distributions, we find that the tunneling state and the higher energy eigenstate close to the tunneling state are populated in both cases. Thus, we conclude that the current oscillation results from transitions between these two states, with the frequency of the oscillation determined by the frequency of these transitions. Since both the first and second eigenstates change in time, the amplitude of current becomes large in the large oscillation case. When the structure of the emitter basin becomes stable with respect to change of the bias voltage, a plateau forms. As the size of the contact regions increases, the transition frequencies decrease, because the size of the emitter basin increases, while the depth of the basin does not change. This accounts for the peak shift seen in figures 3(i-a), (iii'-a) and 3(ii-a), (iii-a). For a small size of contact regions, the emitter basin becomes more stable with respect to change in the bias voltage and as the result the plateau becomes larger.

Both figures 2(ii) and (iii) exhibit double plateau-like features in the NDR region, but we find current oscillation only in the upper plateau in the case of figure 2(ii). This is because the dissipation in the large oscillation case of figure 2(ii) is not sufficiently strong to create significant population of the tunneling states in the case of a small basin, for which the resonant frequencies between the tunneling states and adjacent states are large. In our previous study [102], we considered the same condition as in figure 2(iii), with a value of γ half as large (= 12.1 THz). In that case, however, there was only upper plateau oscillation, while we observed the bistability and the double plateau-like feature. This is because the effective system–bath coupling strength in the previous case is weaker than in the present case [38]. Thus, in that case, the ground tunneling states are not populated from the conduction band through dissipation. Note, however, that if the system–bath coupling is too strong, dissipation suppresses the current oscillation, as explained in appendix B.

One important aspect of the present methodology is that it allows elucidation of the dynamical behavior of the system through the time evolution of the Wigner distribution function. We have been able to characterize the patterns of the time evolution of the Wigner distribution for two types of current oscillations, with large and small oscillation amplitude. Here, we describe these two types in detail, using as one reference, the large oscillation case in figure 2(iii-a) and the small oscillation case in figure 2(ii-a). In figure 4, we display snapshots of the Wigner distribution function for the case of large oscillation at the times marked on the red curve in figure 2(iii-b). As illustrated in figures 3(i-b) and (iii'-b), the characteristic feature of this type of oscillation is the large electron density near the emitter side of the barrier, which is observed as a distinct peak separated from the conduction band (at q =  18 nm in figure 3(i-b) and at q =  43 nm in figure 3(iii'-b)). As a result of this feature, the effective potential possesses a deep emitter basin next to the barrier. The profiles of the eigenstates depicted in figure 3(iii'-b) indicate that this peak consists of the first (red) to fourth (orange) eigenstates. We find that a small peak near the edge of the conduction state (at q =  35 nm in figure 3(iii'-b)), which arises from the third (blue) and fourth (orange) eigenstates in figure 3(iii'-b), also plays an important role in the current oscillation. In the Wigner distribution plotted in figure 4(a), these two peaks are denoted by A and B, respectively. In the situation depicted in figure 4(a), the current flows into the system from the emitter side of the boundary, and then it is scattered by the emitter side of the barrier almost elastically, because the kinetic energy of the current electron is much higher than that of the tunneling state. The scattered current is trapped by the emitter basin and rotates clockwise around the peaks A and B, but, due to dissipation, some of it flows into peak B while losing energy. In the eigenstate representation given in figure 3(iii'-b), this behavior corresponds to population transfer from the fourth (orange) to the third (blue) eigenstate. In figure 4(b), when the third (blue) state decays further to the first (red) and second (green) tunneling states, the height of A increases. As a result, peak A becomes higher, while peak B becomes lower. In figure 4(c), because peak A is related to the tunneling state, the outflow current becomes larger whenever peak A becomes higher. Throughout this process, the peaks A and B become higher and lower by turns, in a manner reminiscent of the piston in a two-piston engine. As a result of this motion, the current exhibits oscillation. This behavior is typical for this large oscillation case. Because there is no current in the case of figure 4(a), the oscillation amplitude is large compared to that in the small oscillation case. If the dissipation is too strong, however, the piston-like motion is suppressed, and there is only steady current, as described in appendix B.

Zoom In Zoom Out Reset image size

In figure 5, we display snapshots of the Wigner distribution for the case of small oscillation at the times marked in figure 2(ii-b). In figure 5(a), while current flows into the system from the emitter side of the boundary, that part of it with energy larger than that of the tunneling state is scattered by the emitter side of the barrier. The remaining current, i.e. that whose energy is closer to the energy of the tunneling state, flows to the collector side in the form of steady current, through tunneling. In figure 5(b), it is seen how the scattered electron moves back and forth in the emitter basin, while losing energy due to dissipation. As a result, the electron flows into peak C, exhibiting tornado-like motion. Figure 5(c) depicts shaking motion of the effective potential that periodically accelerates the electron distribution in peak C to the tunneling state. Through this effect, current flows to the collector side through the barrier. Due to synchronization with this shaking motion, the current is enhanced periodically. This tornado-like motion is typical for the small oscillation case. Because there is a large contribution from steady current, the oscillation amplitude here is smaller than in the large oscillation case.

Zoom In Zoom Out Reset image size

In summary, we investigated current oscillations in the plateau structures of the NDR region for three sizes of the contact regions with a model that includes damping, employing the Caldeira–Leggett Hamiltonian. We found two distinct types of current oscillations. The first type is observed in the single plateau and in the lower part of the double-plateau structure. It is characterized by a large oscillation amplitude and a single Fourier component. The other type is observed in the upper part of double-plateau structure. It is characterized by a small oscillation amplitude and two Fourier components. An emitter basin that forms on the emitter side of the effective potential plays a key role in creating the current oscillation. Eigenstate analysis indicates that the first type is caused by transitions between the ground tunneling state and the adjacent excited state in the emitter basin, while the second type is caused by transitions between the intermediate tunneling state and higher states. Because the transition frequencies are large in the case of narrow emitter basin, there is high frequency oscillation in the case of small contact regions for a fixed basin depth. In Wigner space, these two types of oscillation are characterized by the two types of motion: two-piston engine-like motion and tornado-like motion. Dissipation plays an important role in the realization of current oscillation. In order for the ground tunneling state to be populated in the case of large oscillation, there must be fairly strong dissipation, whose strength is determined by the system–bath coupling and the noise correlation time. If the dissipation is too strong, however, the current oscillation vanishes due to damping. The key to have non-trivial behaviors such as hysteresis, single/double plateaus and self-excited current oscillations is on the existence of the resonant tunneling states, the charge redistribution effects and the dissipation. The present results may be helpful to design nano-devices including a molecular junction system [107].

Although many efforts have been made to improve Wigner transportation theory [100, 101], the quantum Boltzmann equation [103, 104] and other formalisms [105, 106] to study quantum dissipative dynamics in nano-devices, there are still a number of limitations and many subtle problems on such formalisms, which deserve further attention. On the basis of the reduced HEOM approach, our investigation was carried out through highly accurate numerical calculations applied to a tunneling device system in a non-Markovian environment at finite temperature. We have provided evidence of intrinsic bistability and self-excited current oscillations in the NDR region rigorously. While the current oscillation experimentally observed in the RTD induced by the resonance with external circuit [52, 53], however, these effects are not accounted for in our approach like many other theories based on the Boltzmann approach. To investigate the relation between the experimentally observed current oscillations and the existence of intrinsic current oscillation, further investigation is necessary. Although the validity of the Caldeira–Leggett Hamiltonian in the description of electron motion is not yet well established, we believe that the present results provide insight into the role of quantum mechanical phenomena in the type of system studied here.

The present approach can be used to treat a strong system–bath coupling non-perturbatively. In addition, any time-dependent external field can be added while taking into account the system–bath quantum coherence through the hierarchy elements. Such features are ideal for studying SQUID rings [9, 10] and quantum ratchet systems [11, 12, 46].

YT is grateful for financial support received from the Humboldt Foundation and a Grant-in-Aid for Scientific Research (B2235006) from the Japan Society for the Promotion of Science. AS acknowledges the research fellowship of Kyoto University.

To understand the origin of hysteresis in the NDR region, we plot the effective potential, electron density and energy eigenstates of the emitter basin and double-barrier well in the case of decreasing bias in figure A.1(a). In this case, the emitter basin is so shallow that there is no tunneling state between the emitter basin and the double-barrier well. The existence of the peak near the barrier indicates that the second excited state is significantly populated. The current arises through the transition from the second excited state to the ground quantum state through dissipation. The Wigner distribution is plotted in figure A.1(b). Due to dissipation, the distribution is in a steady state. The peak near the barrier arises because the barrier impedes the flow. The electron density then leaks to the collector side without oscillation through tunneling in the form of steady current.

Zoom In Zoom Out Reset image size

To see the effect of dissipation, we determined the I–V characteristics for the case of a stronger system–bath coupling, ζ = 120.8 GHz (ζ⁻¹ = 8.28 ps) (figure B.1). The values of the other parameters are the same as in figure 2. In contrast to the case considered in figure 2, in the present case current oscillation is not observed even in the plateau region. Also, the upper plateau does not exist. This is because the current oscillation decays quickly through the damping. Because the heat bath is strongly coupled to the system, the eigenstate picture of the electron system itself, as depicted in figure 3(b), is of questionable validity. As a result, the plateau is lost.

Zoom In Zoom Out Reset image size

For the sake of comparison, we present the I–V characteristics calculated from the Boltzmann equation and the Poisson equation [58–63] for the same physical conditions as in figure 2. The Boltzmann equation commonly used in the RTD problem is expressed as [59, 60]

$$\begin{equation} \frac{\partial}{\partial t} W(p, q; t) = - \skew3\hat{L}_{\mathrm{qm}}W(p, q; t) + \left( \frac{\partial W(p, q; t)}{\partial t} \right)_{\mathrm{coll}}, \end{equation} \tag{ C.1 }$$

where $\hat {L}_{\mathrm {qm}}$ is the quantum Liouvillian defined by equations (16) and (17) and

$$\begin{equation} \left( \frac{\partial W(p, q; t)}{\partial t} \right)_{\mathrm{coll}}=-\frac{1}{\tau} \left( W(p, q; t)- \frac{P(q; t)}{P_{\mathrm{eq}}(q)} W_{\mathrm{eq}}(p, q) \right)\end{equation} \tag{ C.2 }$$

is the modified collision operator under the relaxation time approximation. Here, τ is the relaxation time, W_eq(p,q) is the equilibrium Wigner distribution function, $P(q; t)=\int \mathrm {d}p W(p,q;t) /(2\pi \hbar )$ is the density of the electron and P_eq(q) is that of the equilibrium distribution, respectively. Because the collision term is determined by the Wigner distribution at time t, W(p,q;t), and does not depend upon the previous history of distribution, this equation describes Markovian dynamics. Because the Boltzmann equation does not have a fluctuation term that is related to a dissipation term through the quantum version of the fluctuation–dissipation theorem, the equilibrium distribution is not an intrinsic state of the Boltzmann equation. Moreover, we have to determine what the equilibrium distribution W_eq(p,q) is in an ad hoc manner.

We solve the Boltzmann equation for the effective potential calculated from equation (23) and the Poisson equation (24) following the same procedure as in section 4 with the same set of system parameters as in figure 2. Here, the equilibrium distribution, W_eq(p,q), is obtained from the quantum Liouville equation for the effective potential with the bias voltage zero [60, 62]. The boundary conditions is given by [56, 101]

$$\begin{equation} W(p,q=0 \ \mathrm{or} \ L) = \frac{m}{\pi\!\hbar^2\beta} \ln \left[1 + \exp\left( -\frac{\beta p^2}{2m}\right) \right]. \end{equation} \tag{ C.3 }$$

In the Boltzmann equation approach, the time constant τ was estimated from other theory and Buot et al set it to be τ = 525 fs at T = 77 K [60]. To compare with the HEOM result, here we solve the Boltzmann equation at T = 300 K for various τ to find the case that exhibits similar I–V profiles as in figures 2(i)–(iii). Note that the difference between the Fermi–Dirac and Bose–Einstein distributions is minor at this temperature. The obtained results for τ = 200 fs are presented in figure C.1. In figures C.1(i)–(iii), while we observed hysteresis, we could not find any plateau-like behavior and current oscillation in the NDR region as was shown in the HEOM calculations in figures 2(i)–(iii). The present results are also different from the result obtained from the Boltzmann equation at T = 77 K for τ = 525 fs, in which a single plateau behavior and current oscillations in the NDR region were observed [59, 62].

Zoom In Zoom Out Reset image size

To analyze the difference between the Boltzmann and HEOM results, we display snapshots of the steady-state Wigner distribution near the maximum and minimum of the I–V curves in figures C.2(i) and (ii). We find that the Wigner distribution in the collector region is smooth in the HEOM case, while there are many small peaks disturbing the flow in the Boltzmann case. This difference arises because the equilibrium state in the Boltzmann approach is fixed even when the effective potential is changed from the original one due to the self-consistent calculations. As the result, the difference between the imposed equilibrium state and the true equilibrium state becomes large especially in the collector region, where the distribution is far from the assumed equilibrium distribution, W_eq(p,q), due to the scattering from the potential barriers.

Zoom In Zoom Out Reset image size

As indicated in this appendix, dynamics described by the Boltzmann equation approach is different from the HEOM approach. This difference arises because the thermal equilibrium state of the Boltzmann equation is introduced as an assumption, while the thermal equilibrium state of the HEOM approach is an intrinsic state of the equation that is determined through the balance between the fluctuation term and dissipation term. This difference becomes significant for a system that exhibits hysteresis, since the equilibrium state of the system depends upon the pathway of process. This difference may also be significant if the system is driven by a time-dependent external field, in which the equilibrium state is not well-defined.