Electrothermal Monte Carlo Simulation of a GaAs Resonant Tunneling Diode

Abstract

This paper deals with the electron transport and heat generation in a Resonant Tunneling Diode semiconductor device. A new electrothermal Monte Carlo method is introduced. The method couples a Monte Carlo solver of the Boltzmann–Wigner transport equation with a steady-state solution of the heat diffusion equation. This methodology provides an accurate microscopic description of the spatial distribution of self-heating and its effect on the detailed nonequilibrium carrier dynamics.

1. Introduction

Due to the continued miniaturization of integrated circuits and the current trend toward nanoscale electronics, power densities, heat generation, and chip temperatures will reach levels that will prevent the reliable operation of such circuits. In order to minimize self-heating effects, the development of accurate electrothermal simulators is required, which takes into account the coupling between electronic and lattice dynamics. In the active regions of such small devices, heat generation is a direct consequence of the nonequilibrium carrier transport. In high electric field regions, the electrons are accelerated and collide with the lattice in such a way that the emission of a large number of phonons contributes to heat transport in the device. In the framework of semiclassical charge transport, electrothermal simulators are based on drift–diffusion or hydrodynamic models [1,2], which are able to capture nonequilibrium transport effects. Alternatively, the direct simulation Monte Carlo (MC) can provide an accurate nonequilibrium charge transport simulation, which is free from the approximations made in the drift–diffusion or hydrodynamic model. Electrothermal Monte Carlo simulators have been developed during these years [3,4,5] but not in quantum regimes where the Boltzmann Transport Equation must be replaced by the Wigner Transport Equation (WTE). Since electron devices are quantum systems outside of thermodynamic equilibrium, scattering by phonons should be included in the WTE for a realistic simulation. Many proposals for the collision operator can be found in the literature [6,7], which provide an accurate description of the phenomena at the price of a high requirement of computational resources. Because of that, the use of such operators is restricted to very simple (idealized) systems. In this paper, the effects of scattering with phonons are taken into account via a semiclassical Boltzmann collision operator, which employs transition rates calculated using Fermi’s golden rule, obtaining the so called Boltzmann–Wigner transport equation (BWTE). Numerical solvers of the WTE can be based on finite-difference schemes [8,9,10,11,12,13], where scattering was restricted to the relaxation time approximation and the momentum space to one dimension. The Monte Carlo method allows for scattering processes to be included on a more detailed level, assuming a three-dimensional momentum–space. In this paper, we shall use the so-called Signed Particle Monte Carlo method (SPMC) [14,15] in which the effect of the Wigner potential is interpreted as a probabilistic generation of couples of positive and negative particles, where the quantum information is carried by their sign. The huge number of generated particles can be controlled by an annihilation process: two particles with an opposite sign entering a given phase space cell are canceled. Recently, this method has also been understood in terms of the Markov jump process theory [16], producing a class of new stochastic algorithms. Algorithms that belong to this class are a standard time-splitting algorithm and a new no-splitting algorithm that avoids errors due to time-discretization [17,18].

Taking advantage of previous Electrothermal Monte Carlo semiclassical models, in this paper, we shall study the heating effect in a Resonant Tunneling Diode (RTD), coupling the SPMC solver of the BWTE with a steady-state solution of the heat diffusion equation. To the author’s knowledge, this model is the first of its kind in terms of model accuracy. The paper is organized as follows. Details of the Boltzmann–Wigner transport equation are provided in Section 2, and in Section 3 we deal with the Signed Particle Monte Carlo method. In Section 3, we introduce the Resonant Tunneling diode structure and in Section 5 the Electrothermal Signed Particle Monte Carlo Method. Simulation results are shown in Section 6, and conclusions are drawn in Section 7.

2. The Boltzmann–Wigner Transport Equation

The BWTE writes [19]

$$\frac{\partial }{\partial t}{f}_w(t,x,k)+\frac{\hslash }{m^{*}}k·{\nabla }_x{f}_w(t,x,k)+\frac{e}{\hslash }{\nabla }_x\phi ·{\nabla }_k{f}_w(t,x,k)=\mathcal{Q}\left(f_w\right)+\mathcal{C}\left(f_w\right).$$

(1)

$x\in {\mathbb{R}}^3$ and $\hslash k\in {\mathbb{R}}^3$ are the electron position and momentum, respectively, $m^{*}$ is the electron effective mass, and $\phi$ the slowly-varying potential satisfying the Poisson equation

$$\nabla ·\left[{\epsilon }_0{\epsilon }_r\nabla \phi \left(x\right)\right]=-e(N_D-N_A-n),$$

(2)

where e is the elementary charge, ${\epsilon }_0$ the absolute dielectric constant, ${\epsilon }_r$ the relative dielectric constant, $N_D$, $N_A$ are the donors and acceptors’ doping profiles, and n the particle density

$$n(t,x)=\int f_w(x,k,t)dk.$$

(3)

$\mathcal{C}\left(f_w\right)$ is the Boltzmann scattering operator which, in the not-degenerate case, is as follows [20]:

$$\begin{array}{c}\hfill \mathcal{C}\left(f_w\right)=\int \left[w_s(k^{\prime },k)f_w\left(k^{\prime }\right)-w_s(k,k^{\prime })f_w\left(k\right)\right]d{k}^{\prime },\end{array}$$

(4)

where $w_s(k,k^{\prime })$ is the scattering rate at which electrons suffer with phonons and impurities, given by the Fermi’s golden rule. The quantum evolution is taken into account by the term

$$\mathcal{Q}\left(f_w\right)={\int }_{{\mathbb{R}}^d}{V}_w(x,k-k^{\prime })f_w(t,x,k^{\prime })d{k}^{\prime },$$

(5)

where $V_w$ is the Wigner potential

$$V_w(x,k)=\frac{1}{i\hslash {\left(2\pi \right)}^d}{\int }_{{\mathbb{R}}^d}d{x}^{\prime }{e}^{-ik·x^{\prime }}\left[V\left(x+\frac{x^{\prime }}{2}\right)-V\left(x-\frac{x^{\prime }}{2}\right)\right],$$

(6)

and V is the rapidly-varying term of the potential energy.

3. The Signed Particle Monte Carlo Method

The quantum evolution term (5) can be interpreted like the Gain term of the collisional operator of the Boltzmann transport equation, in which the Loss term is missing. However, the Wigner potential (6) is not always positive and, for this reason, cannot be considered a scattering term. The main idea of the Signed Particle Monte Carlo method [14] consists of separating $V_w$ into a positive and negative part $V_w^{+},V_w^{-}$ such that

$$V_w=V_w^{+}-V_w^{-},V_w^{+},V_w^{-}\ge 0$$

(7)

Consequently, we can define an integrated scattering probability per unit time as

$$\gamma \left(x\right)=\int d{k}^{\prime }{V}_w^{+}(x,k-k^{\prime })=\int d{k}^{\prime }{V}_w^{-}(x,k-k^{\prime })$$

(8)

and rewrite the quantum evolution term as the difference between Gain and Loss terms, i.e.,

$$\begin{array}{ccc}\hfill \mathcal{Q}\left(f_w\right)& =& \int d{k}^{\prime }w(k^{\prime },k)f_w(t,x,k^{\prime })-\gamma \left(x\right)f_w(t,x,k)\hfill \end{array}$$

(9)

$$\begin{array}{ccc}\hfill w(k^{\prime },k)& =& V_w^{+}(x,k-k^{\prime })-V_w^{-}(x,k-k^{\prime })+\gamma \left(x\right)\delta (k-k^{\prime }).\hfill \end{array}$$

(10)

The interpretation of the scattering term $w(k^{\prime },k)$ is that a particle produces, in the same position, a couple of new particles with weight u and $-u$ according to a generation rate given by the function $\gamma \left(x\right)$. The momentum of the new particles is generated with probability $V_w^{+}(x,k)/\gamma \left(x\right)$. Since usually $\gamma$ is rapidly oscillating, an exponential growth of particle numbers is expected and, in order to control the particle number, a cancellation procedure is mandatory.

This procedure has been understood using the theory of the piecewise deterministic Markov processes [16], where the state space is

$$z_j\left(t\right)=(u_j\left(t\right),x_j\left(t\right),k_j\left(t\right)),t\ge 0,j=1,...,N\left(t\right),$$

(11)

and $u_j\in \{-1,+1\}$ is the weight. The time evolution of the particle system (11) is assigned by a deterministic motion according to the flow

$$F(t,z)=(u,x+v\left(k\right)t,k),v=\frac{\hslash }{m}k$$

(12)

and a jump kernel $Q\left(z_j\left(t\right)\right)$. The random waiting time $\tau$ until the next jump satisfies

$$\mathbb{P}(\tau \ge t)=exp\left(-{\int }_0^{t}Q\left(F(s,z)\right)ds\right).$$

(13)

For numerical purposes, we introduce a majorant ${\widehat{V}}_w$ such that

$$|V_w(x,k)|\le {\widehat{V}}_w(x,k)\forall x,k\in {\mathbb{R}}^3.$$

(14)

If the j-th particle generates two new particles with

$$z_1^{{}^{\prime }}=\left(u_j\mathrm{sign}{V}_w(x_j,k),x_j,k_j+k\right),z_2^{{}^{\prime }}=\left(u_j\mathrm{sign}{V}_w(x_j,k),x_j,k_j-k\right)$$

(15)

the jump kernel takes the form [17]

$$Q\left(z_j\right)=\frac{1}{2}\int {\widehat{V}}_w(x_j,k)dk$$

(16)

and Equation (13) writes

$$\mathbb{P}(\tau \ge t)=exp\left(-{\int }_0^t\widehat{\gamma }(x_j+v\left(k_j\right)s)ds\right)$$

(17)

where

$$\widehat{\gamma }\left(x\right)=\frac{1}{2}\int {\widehat{V}}_w(x,k)dk$$

(18)

represents the generation probability. It is possible to prove that functionals of the solution of the Wigner equation are expressed in terms of the particle system using the representation [16]

$$\int \int \varphi (x,k)f(t,x,k)dkdx=\frac{1}{N_{ini}}\mathbb{E}\left(\sum _{j=1}^{N\left(t\right)}{u}_j\left(t\right)\varphi (x_j\left(t\right),k_j\left(t\right))\right)$$

(19)

where $\varphi$ is an appropriate test function, and $N_{ini}=N\left(0\right)$ is the initial particle number. In order to to separate the transport and the jump processes, usually a splitting time step $\mathsf{\Delta }t$ is used at the expense of a discretization error. This can be avoided by using a no-splitting algorithm recently introduced in [17,18]. By introducing a majorant for the generation process (8) and one for the total phonon scattering rate

$${\mathsf{\Gamma }}_s\ge max\lambda \left(k\right),\lambda \left(k\right)=\sum _{\alpha }\int w_{\alpha }(k,k^{\prime })d{k}^{\prime }$$

(20)

The total majorant is

$$\mathsf{\Gamma }={\mathsf{\Gamma }}_s+\widehat{\gamma }$$

(21)

and Equation (17), for all particles, now is as follows:

$$\mathbb{P}(\tau \ge t)=exp\left(-\sum _{j=1}^N{\int }_0^t\left[{\mathsf{\Gamma }}_s+\widehat{\gamma }(x_j+v\left(k_j\right)s)\right]ds\right).$$

(22)

In the case in which $\widehat{\gamma }$ does not depend on the position, we have

$$\mathbb{P}(\tau \ge t)=exp\left(-\mathsf{\Gamma }Nt\right)\to \tau =-\frac{1}{\mathsf{\Gamma }N}logr$$

(23)

where $r\in U[0,1]$, and $\tau$ is completely determined. With respect to the splitting case, now the transport and the generation process can not be separated, and the results shall not be affected by any discretization error.

4. The Resonant Tunneling Diode

A standard Resonant Tunneling Diode structure [21] has been implemented, as shown in Figure 1. The barriers have depth $b=$ 3 nm, height $a=$ 0.3, and the quantum well dimension is $b_w$ = 5 nm, symmetric with respect to the mid-point $L/2$ (total length $L=150\mathrm{nm}$). The barrier structure is embedded in a 30 nm lightly doped region ($N_D={10}^{16}{\mathrm{cm}}^{-3}$) which is connected to 60 nm highly doped regions on either side ($N_D^{+}={10}^{18}{\mathrm{cm}}^{-3}$).

In this case, the Wigner potential (6) can be easily evaluated in addition to the majorant (18) (see [18] for the details). The device considered is made by Gallium Arsenide (GaAs) (with $m^{*}=0.067$), and polar optical phonons (POP) within a single $\mathsf{\Gamma }$ band [20] in the parabolic band approximation used are taken into account. The total scattering rate is written as follows [20]:

$$\lambda (k,T_L)={\lambda }^{-}(k,T_L)+{\lambda }^{+}(k,T_L)$$

(24)

where the first term represents POP absorption and the second one emission

$${\lambda }^{-}(k,T_L)=\frac{e^2{\omega }_p\left(\frac{1}{{ϵ}_{\infty }}-\frac{1}{{ϵ}_r}\right)}{2\pi {ϵ}_0\hslash \sqrt{\frac{2\epsilon \left(\mathsf{k}\right)}{m^{*}}}}{n}_0{sinh}^{-1}\sqrt{\frac{\epsilon \left(k\right)}{\hslash {\omega }_p}}$$

(25)

$${\lambda }^{+}(k,T_L)=\frac{e^2{\omega }_p\left(\frac{1}{{ϵ}_{\infty }}-\frac{1}{{ϵ}_r}\right)}{2\pi {ϵ}_0\hslash \sqrt{\frac{2\epsilon \left(k\right)}{m^{*}}}}(n_0+1){sinh}^{-1}\sqrt{\frac{\epsilon \left(\mathsf{k}\right)}{\hslash {\omega }_p}-1}$$

(26)

The term $n_0\left(T_L\right)$ is the phonon equilibrium distribution, i.e.,

$$n_0\left(T_L\right)=\frac{1}{exp\left(\frac{\hslash {\omega }_p}{k_B{T}_L}\right)-1}$$

(27)

$\hslash {\omega }_p$ is the polar optical phonon energy (0.03536 eV) and $T_L$ the lattice temperature. The initial lattice temperature is 300 K, and ohmic boundary conditions are used.

5. The Electrothermal Signed Particle Monte Carlo Method

An important issue that arises from the coupling of an MC electronic transport algorithm to any thermal model is the significant difference in the characteristic time scales of electronic and thermal transport. Electronic transients in GaAs systems are of the order of picoseconds, whereas thermal transients may be of the order of nanoseconds, microseconds, or even longer. Performing MC computations for the duration of thermal transients across the whole semiconductor die would not be feasible. Consequently, the method used in this paper extracts steady-state electrothermal device characteristics only. The electrothermal SPMC method of simulation is an iterative approach:

• The initial SPMC iteration is run at a room temperature of 300 K for a few ps, in order to reach a steady-state;
• As the steady state is reached, electronic parameters are sampled for typically 15 ps, in order to evaluate the heat generation rate $H\left(x\right)$;
• The lattice temperature $T_L\left(x\right)$ is obtained by solving the steady-state heat diffusion equation

  $${\nabla }_x\left(\kappa {\nabla }_x{T}_L\left(x\right)\right)+H\left(x\right)=0$$

  (28)
  $\kappa$ being the thermal conductivity in GaAs;
• The SPMC solver is rerun, in the next iteration, with the new lattice temperature $T_L\left(x\right)$. We observe that the scattering rates (25) and (26) depend on the lattice temperature;
• We repeat this procedure until convergence is reached.

This model does not account for temperature changes beyond the semiconductor die. Radiation losses are neglected, as their contribution at the small die surface areas is insignificant.

The mechanism through which Joule heating occurs is that of electron scattering with phonons, and consequently only a simulation approach which deliberately incorporates all such scattering events will capture the complete microscopic, detailed picture of lattice heating. The phonon emission and absorption events during a simulation run are tallied and full heat generation statistics can be collected. We wait until the steady state has been reached at time $t_0$. Then, we count our events in the observation points $t_i,i=0,...,N_{obs}$. We evaluate the heat generation rate in two ways:

• Counting the phonon number.
  We introduce the quantity [22]

  $$H^c(t_{i-1},t_i,x)=\frac{n(t_i,x)}{N_p(t_i,x)}\frac{\hslash {\omega }_p\left[C^{+}-C^{-}\right]}{dt},$$

  (29)

  where $C^{+}(t_{i-1},t_i,x)$, $C^{-}(t_{i-1},t_i,x)$ are the numbers of the phonon emitted and absorbed in the time interval $(t_{i-1},t_i)$ in the x-th grid point, $n(t,x)$ the charge density, and $N_p(t,x)$ the particle number at time t in the x-th cell. Then, the heat generation rate is

  $$〈H^c\left(x\right)〉=\frac{1}{N_{obs}}\sum _{i=1}^{N_{obs}}{H}^c(t_{i-1},t_i,x)$$

  (30)
• using the integrated probability scattering function.
  From the integrated probability scattering (25) and (26) we can define

  $$H^F(t_i,x)=\frac{n(t_i,x)}{N_{ini}}\sum _{j=1}^{N\left(t_i\right)}{u}_{j}G\left(\epsilon \left(k_j\right)\right),G\left(\epsilon \right)=\hslash {\omega }_p\left[{\lambda }^{+}\left(k\right)-{\lambda }^{-}\left(k\right)\right]$$

  (31)
  Then, the heat generation rate is

  $$〈H^F\left(x\right)〉=\frac{1}{N_{obs}}\sum _{i=1}^{N_{obs}}{H}^F(t_i,x).$$

  (32)

The heat generation is reduced to the usual calculation of functionals according to Equation (19). This estimator enjoys better approximation properties due to reduced statistical fluctuations [5].

6. Numerical Results

In order to have a significant lattice temperature increase with respect to the equilibrium temperature of 300 K, the applied bias voltage $V_b$ must be greater than $0.8$ V. In Figure 2, we plot the heat generation rate versus the position, evaluated by means of the counting estimator (30) and the integrated probability estimator (32), for $V_b=0.8$ V. From this figure, we can see that the maximum heat is produced inside the quantum well region, representing a so-called hot spot region.

In Figure 3, we plot the corresponding standard deviation, proving the variance reduction of the integrated probability estimator (32). In Figure 4, we plot a zoom of Figure 2 with the error bar, proving that the integrated probability estimator is always inside the tolerance band of the counting estimator. Figure 5 shows the density for the first two iterations, showing no appreciable variation.

Figure 6 shows the lattice temperature evaluated by means of the heat diffusion Equation (28) for the first two iterations, which are enough to reach the convergence. We observe that the lattice temperature is decreasing with the iteration number. To explain this behavior, one must consider the function $G\left(\epsilon \right)$ in Equation (31). This function represents the difference between the emitted and absorbed phonon probability; if this quantity is positive, more phonons are released into the lattice and in turn the temperature increases. We plot this quantity in Figure 7 showing that, for this particular kind of scattering mechanism, it decreases with the lattice temperature. In Figure 8, we plot the current versus the iteration number, proving that this quantity is constant. If we double the applied voltage to $V_b=1.6$ V, the increase of temperature is of a factor 5 as shown in Figure 9.

7. Conclusions

The Electrothermal Signed Particle Monte Carlo algorithm provides an accurate tool for studying heat generation and quantum effects in nanometric semiconductor devices, at the expense of huge computational effort. The coupling between the MC charge transport and the heat diffusion equation is given by a term called heat generation rate obtained, usually, by counting the number of phonons emitted/absorbed during the steady-state. Alternatively, a new estimator of the heat generation rate, based on the integrated scattering probability function (32), can be used, which enjoys reduced statistical fluctuations. Simulation results for a Resonant Tunneling Diode are shown, proving that the heat is produced almost entirely inside the quantum well and estimating the lattice temperature, which depends on the applied voltage. The localization of hot spot regions can be useful in the design of such devices, in order to optimize the heat removal.