Quantum Monte Carlo simulation of resonant tunneling diodes based on the Wigner distribution function formalism

A tool for the simulation of resonant tunneling diodes (cid:126) RTDs (cid:33) has been developed. This is based on the solution of the quantum Liouville equation in the active region of the device and the Boltzman transport equation in the regions adjacent to the contacts by means of a Monte Carlo algorithm. By accurately coupling both approaches to current transport, we have developed a quantum simulation tool that allows the use of simulation domains much larger and realistic than those previously considered, without a signiﬁcant increase in computational burden. The main characteristics expected for the considered devices are clearly obtained, thus supporting the validity of our tool for the simulation of RTDs. © 1998 American Institute of Physics


Quantum Monte Carlo simulation of resonant tunneling diodes based on the Wigner distribution function formalism
J. García-García, F. Martín, a) X. Oriols, and J. Suñ é Departament d'Enginyeria Electrònica, Universitat Autònoma de Barcelona 08193 Bellaterra (Barcelona), Spain ͑Received 27 July 1998; accepted for publication 14 October 1998͒ A tool for the simulation of resonant tunneling diodes ͑RTDs͒ has been developed. This is based on the solution of the quantum Liouville equation in the active region of the device and the Boltzman transport equation in the regions adjacent to the contacts by means of a Monte Carlo algorithm. By accurately coupling both approaches to current transport, we have developed a quantum simulation tool that allows the use of simulation domains much larger and realistic than those previously considered, without a significant increase in computational burden. The main characteristics expected for the considered devices are clearly obtained, thus supporting the validity of our tool for the simulation of RTDs. © 1998 American Institute of Physics. ͓S0003-6951͑98͒00550-6͔ The scaling down of vertical and lateral dimensions of electron devices up to the nanometer scale, associated with recent progress in heteroepitaxial growth and nanopattern fabrication techniques, has opened the door to design strategies based on quantum-mechanical phenomena ͑such as tunneling or quantum interference effects͒. Among these devices, resonant tunneling diodes ͑RTD's͒ have received much attention due to their potential applications as microwave sources and to high-speed electronics. 1,2 To predict and explain device behavior, as well as to aid in RTD-based circuit design, the development of efficient and reliable simulation tools is of great importance. In this regard, the nonequilibrium Green function theory due to Lake et al. 3 and the Wigner distribution function ͑WDF͒ formulation of quantum mechanics [4][5][6][7] have been demonstrated to be a powerful basis for the numerical simulation of vertical transport quantum devices, in which dissipative effects and self-consistency have to be considered. For instance, the intrinsic bistability and current oscillations that are experimentally found in resonant tunneling diodes biased in the negative conductance region have been reproduced in the framework of the WDF formalism. 4 However, the computational burden associated with the iterative solution of the Liouville and Poisson equations has imposed severe limitations on the spatial domains, these being too small for reliable simulations of quantum electron transport in RTD's. 8 The need to extend the simulation domain is not related to the fact that quantum effects are present at the boundaries ͑since scattering events tend to destroy phase coherence͒, but to the fact that the usual open system boundary conditions 9 are not applied at boundaries where the potential profile is flat. To overcome this computer time-related limitation, one possibility is to extend the simulation domain by solving the discretized Liouville equation in a region of the device containing the double barrier ͓quantum window ͑QW͔͒, and the Boltzmann transport equation ͑BTE͒ using a Monte Carlo ͑MC͒ solver in the external regions, where a classical behavior for current transport is expected. With this in mind, we have developed a onedimensional simulation tool applicable to any vertical transport quantum device. In what follows, we will discuss how the classical Monte Carlo algorithm and the Liouville solver have been coupled, and the main results, derived on the basis of RTD's. Although the idea of such a coupling was previously advocated by Salvino and Buot, 10 it was never subsequently developed.
Self-consistency has been implemented by iteratively solving the Poisson's equation in the whole device and the Liouville/Boltzmann equations in the quantum/classical regions. In the quantum window, the Liouville equation ͑including dissipative effects by means of a relaxation time term 11 ͒ is solved by discretizing the phase space, so that the WDF is obtained from a matrix equation at each time step, provided the boundary conditions are known. 12 These are inferred from the Monte Carlo distribution of carriers projected to the component of the wave vector in the direction of motion in the adjacent cells to the QW. A transformation of the Monte Carlo distribution to the WDF formalism is needed. For this purpose, the classical regions are discretized in k space with identical intervals as those considered in the QW. In this way, the distribution of carriers f ki can be obtained. From it, the boundary condition for the WDF solver is found to be where ⌬x, ⌬k, and A are the cell width, the k-space interval, and the cross-sectional area of the device, respectively. As previously indicated, the BTE is solved by means of a Monte Carlo algorithm in the classical regions. At each iteration step, charge carriers have to be injected from the reservoirs ͑external contacts to the device͒ and from the QW. Injection from the external contacts has been implemented according to an Ohmic contact model, 13 where charge neutrality in the adjacent cells to the boundaries dictates carrier injection. Obviously, this model cannot be applied for carrier injection from the QW since, in general, charge neutrality does not hold at the boundaries between the QW and the classical regions. On the other hand, injection from the QW to the classical region has to be dependent on the behavior of carriers in the neighboring to the classical regions. We have a͒ used the portion of the WDF at the QW interfaces with wave-vector pointing to the classical regions, weighted by carrier momentum, as the injecting distribution. According to it, the component of the momentum of injected carriers in the direction of motion has been randomly determined, while a Maxwellian distribution has been used to obtain the transverse components. The number of carriers injected at each time step, has been obtained from where q is the electron charge, ⌬t the time step, and J 1/2 is a current density due to carriers entering the classical region from the QW. This current is obtained from the WDF at the boundary according to the expression provided by the WDF formalism, 9 but limiting the integral to the portion of the k axis corresponding to carriers flowing to the classical region. Our tool has been applied to the simulation of RTD's with the aim to demonstrate the advantages of the technique over previous simulators based on the WDF, and to point out the potential applications of it to the simulation of quantum devices requiring extensive simulation domains. In this regard, we have simulated the I -V characteristic ͑Fig. 1͒ of a RTD with a deliberately small nominal doping level, since under these conditions the electric field is expected to extend above the boundaries of the QW, and therefore, the device is an appropriate test structure to analyze the behavior of our simulation tool. At the moment, scattering in the QW has been purposely switched off to magnify the quantum effects. Figure 1 has been obtained by the application of successive voltage steps of 0.02 V increments to the structure. The evolution of the current up to steady state that results after the sudden voltage changes ͑averaged from the beginning of the switch for noise filtering͒ shows that a stable current is approximately achieved after 500 iterations. Also in Fig. 1, the static I -V curve that results without the coupling to the Monte Carlo algorithm, i.e., by considering only the quantum Liouville equation applied to the QW, is shown for comparison. The differences in the peak-to-valley ratio and resonant current reveal the limitations of previous simulation tools, which are based on the solution of the Liouville equation in small integration boxes ͑due to computational bur-den͒, and are not able to reproduce the more reliable results obtained by extending the simulation domains according to the technique proposed in this work. Besides these improvements in the reliability of the simulation results, the fundamental point to emphasize is that no significant CPU time is added by coupling the Liouville solver to the classical Monte Carlo algorithm. This is so because, for typical Monte Carlo parameters, it is the WDF that dictates the rate at which the simulations proceed. As an example, for the conditions considered in Fig. 1, the time needed by the Liouville solver to give a single WDF is higher than the time required by the Monte Carlo scheme to move the particles during an iterative step by a factor of 10. This is very important since improvement in accuracy can be achieved without a penalty in efficiency.
In Fig. 1͑b͒, the steady-state potential profile and electron concentration obtained at resonance ͑0.17 V͒ show that the main qualitative features expected for the devices are reproduced. In particular, the flatband at the extremes of the devices reveals charge neutrality at the external contacts, contrary to what occurs if smaller domains are considered. In view of Fig. 1͑b͒, charge accumulation in the well results at resonance, in accordance with quantum-mechanical consid- erations. Note that the electron density peak in the emitter is slightly separated from the edge of the emitter barrier, such as one expects for a quantum system. It is also very important to mention the smooth transition of the charge density and potential profile at the QW boundaries, which means that the coupling model is well behaved. The momentum distribution of carriers at resonance ͓Fig. 1͑c͔͒ shows that the tunneling ridge ͑already found in previous WDF simulations 14 ͒ is also reproduced, although it progressively vanishes due to thermalization of carriers in the collector.
We have also simulated a RTD with more realistic device parameters and accounted for scattering in the QW through a momentum relaxation time reasonable for operation at room temperature. Specifically, a nominal doping level of 2ϫ10 18 cm Ϫ3 and a well width of 6 nm have been considered. The price to pay at these high doping levels is the use of a large number of MC particles per cell to start up the simulations and a much lower time step, which means that longer simulation times are needed. These simulation conditions are required for the elimination of noise-related instabilities, which are expected to be magnified at high doping levels, since the charge ''carried'' by each MC particle increases. In spite of this, the potential profile has to be treated with some caution to avoid divergence of the current. In this regard, the actualization of the potential profile has been done by averaging the potential profile of the previous iteration with the new one, giving a higher weighing factor to the former. In such a way, divergences are avoided, although time evolution towards steady state can be slightly affected by this artifact. The static I -V curve, obtained by the evolution of the average current up to steady state after the application of successive 0.04 V step increments, is depicted in Fig. 2͑a͒. In the inset of Fig. 2͑a͒, the evolution of the current is shown without time averaging. In spite of the presence of noise, an oscillatory behavior can be clearly observed in the negative conductance region. In our opinion, this behavior ͑found previously by other authors 4 ͒ is not due to instabilities of the tool, but to the continuous filling and depletion of the electron density in the well driven by self-consistency. This conclusion is supported by the fact that out of the negative conductance region, instabilities are not found. The frequency of oscillations ͑approximately 15 THz͒ has to be taken as orientative since the free evolution of the system at high doping levels is limited due to the averaging of the potential profile. Work is in progress to relax this bond in order to use the simulator not only to obtain the dc I -V characteristic of RTDs, but also to make accurate predictions of the time evolution of these devices. Finally, in Fig. 2͑b͒, the electron density and potential profile obtained in the valley are represented. The higher doping level gives a more pronounced transition of the electrostatic potential, and charge depletion in the well results, as is expected, under the considered applied bias.
In conclusion, we have presented a tool for the simulation of RTDs, also extensible to other vertical transport quantum devices. On the basis of the results obtained on RTDs we have demonstrated the validity of our technique, since the main characteristics expected for these devices have been perfectly reproduced. We have pointed out the advantages that our simulator offers over previous tools based on the WDF, which are related to the consideration of larger simulation domains without paying a penalty in computational efficiency. In view of the obtained results, our proposal is a good candidate for the simulation of vertical transport quantum devices oriented to device design optimization.
This work has been supported by the Dirección General de Investigación Científica y Técnica ͑DGICYT͒ under project Contract No. PB94-0720.