Simulation of THz Oscillations in Semiconductor Devices Based on Balance Equations

Instabilities of electron plasma waves in high-mobility semiconductor devices have recently attracted a lot of attention as a possible candidate for closing the THz gap. Conventional moments-based transport models usually neglect time derivatives in the constitutive equations for vectorial quantities, resulting in parabolic systems of partial differential equations (PDE). To describe plasma waves however, such time derivatives need to be included, resulting in hyperbolic rather than parabolic systems of PDEs; thus the fundamental nature of these equation systems is changed completely. Additional nonlinear terms render the existing numerical stabilization methods for semiconductor simulation practically useless. On the other hand there are plenty of numerical methods for hyperbolic systems of PDEs in the form of conservation laws. Standard numerical schemes for conservation laws, however, are often either incapable of correctly handling the large source terms present in semiconductor devices due to built-in electric fields, or rely heavily on variable transformations which are specific to the equation system at hand (e.g. the shallow water equations), and can not be generalized easily to different equations. In this paper we develop a novel well-balanced numerical scheme for hyperbolic systems of PDEs with source terms and apply it to a simple yet non-linear electron transport model.


Introduction
In recent years instabilities of electron plasma waves in high-electron-mobility transistors (HEMT) have attracted a lot of attention as a possible solution for closing the THz gap [4,7,8]. By solving the full Boltzmann Transport Equation (BTE) the behavior of such devices can be simulated in a very accurate way [13]. As transient simulations of the BTE however are computationally very expensive, especially for two or three dimensional devices, simpler transport models are needed. By using balance equations derived from the BTE together with appropriate closure relations, models of arbitrary complexity can be specified [1,11]. The simplest of these models is the well-known Drift-Diffusion (DD) model, where only two moments of the distribution function of the BTE are considered. However, the standard DD model is not able to capture plasma wave effects, because the time derivative in the current constitutive equation is usually neglected [16]. Furthermore, two moments are not enough since heating effects of the electron gas can not be described. The closure relation is also too simplistic and inconsistent with a non-zero current density [20].
By using more balance equations as well as more advanced closure relations, better models (e.g. the full hydrodynamic model) can be developed, at the expense of higher complexity. Due to non-linear terms such models are also quite difficult to stabilize numerically, and the well-known Scharfetter-Gummel scheme [27] used for the standard DD model can not be applied anymore.
Due to the strong built-in electric fields present in semiconductor devices at the transitions between differently doped regions, the balance equations considered also contain strong and in some regions even dominating source terms, where the electron density will change by many orders of magnitude on a small length scale. Standard numerical schemes from computational fluid dynamics (CFD) are not able to capture these strong source terms in a stable way, and often lead to unphysical negative densities.
In this paper we develop a novel numerical scheme specifically for the purpose of simulating electron plasma oscillations by solving balance equations. By splitting the equations into a stationary and a dynamic part, well-balancing can be achieved even for enormous source terms, where conventional numerical schemes fail. To that end we use a one-dimensional modified hyperbolic DD model (sometimes referred to as the isothermal hydrodynamic model) including non-linearity in the flux term due to a convective derivative as a testbed for the scheme, while maintaining a simple extensibility to more advanced models including further balance equations.

Model
We consider a double-gate metal-oxide-semiconductor field-effect transistor (MOSFET) with a quasi two-dimensional channel homogeneous perpendicular to the transport direction and ideal ohmic source and drain contacts shown in Fig. 1. We solve the two-dimensional Poisson equation (PE) together with a one-dimensional DD model including an additional time derivative and a convective term in the constitutive current equation. Holes can be neglected in this case. The PE is given by [28] where is the electrostatic permittivity, ϕ the quasistationary potential, e the elementary charge, N D the donor concentration and n is the electron density. At the contacts Dirichlet where V S , V D , V G are the respective source, drain and gate voltages and V T = k B T e is the thermal voltage, with the Boltzmann constant k B and the temperature T . φ MS is a constant work function difference, and n int is the intrinsic carrier density.
The DD model consists of the continuity and constitutive equation [6] ∂ t n + ∂ x j = 0 ( 3 ) where j is the electron current density, c = k B T m * is the thermal velocity with the effective mass m * . The temperature T is fixed to the lattice temperature. Furthermore, E = −∂ x ϕ is the longitudinal electric field in the channel and τ is the macroscopic relaxation time [1], which, for simplicity, is also assumed to be constant. Note that the inclusion of the convective term ∂ x j 2 n leads to a non-linear equation, which complicates the numerical solution process immensely.
At the contacts we use a Dirichlet condition for the electron density (ideal ohmic contacts) [28] n ct = n int exp arsinh N D 2n int .
The terminal currents including displacement currents are calculated by the Ramo-Shockley theorem [17,26,29]. It should be noted that this model can only serve as a rough approximation and is not sophisticated enough to model transport in nanoscale devices [23], especially in the case of high mobilities, that are necessary for plasma waves (quasi-ballistic transport). However, we chose this model because it is simple due to its inclusion of only two moments of the distribution function, while it still contains some non-linearity and can describe plasma instabilities [7]. This alleviates the development and validation of a suitable numerical solver, which can then be generalized to more advanced models as e.g. the full hydrodynamic model or other models with even higher complexity.
A critical limitation is the assumed constant thermal velocity, leading to the formation of supersonic regions in the device at sufficiently large applied voltages, when the drift velocity v = j n becomes larger than c. The transition between a sub-and a supersonic region as well as the opposite case, which leads to stationary shock-waves [21], is hard to capture in a numerical scheme. We will therefore limit our analysis to the subsonic case |v| < c near equilibrium by restricting the allowed drain-source voltage to rather small values, which leads to smooth solutions [6]. This also means that the maximum voltage we can apply to savely stay in the subsonic regime is limited by the temperature through c.
We write eq. (3) and eq. (4) in conservation form [30] Equations of this form are frequently encountered in the field of computational fluid dynamics, for example the shallow-water equations. Due to the presence of a source term SU on the right-hand side, special care must be taken to achieve so called well-balancedness [12,18]; meaning that stationary solutions, where the flux derivative is exactly canceled out by the source term, should also be exactly preserved in the discrete system. This is especially difficult in our case, because due to the large differences in the doping concentration the density can vary over many orders of magnitude in a small region, resulting in large built-in fields. A simple implementation of the finite volume scheme [10] will then result in large errors in these regions, because the flux derivative and the source term are discretized in different ways (difference of two fluxes compared to numerical integral of the source term over the cell), which can lead to numerical instability.
Since we want to perform transient simulations of plasma oscillations around a steady state, we do not need to achieve well-balancedness for all possible stationary solutions, but only for that particular steady state. To that end we split U and S as well as ϕ into a stationary and a dynamic part where U s (x) is the solution of the stationary equation which needs to be solved prior to a transient simulation. The stationary and dynamic source matrices are given by The stationary and dynamic potentials both have to satisfy the corresponding version of eq. (1): In the case of the extended DD model eq. (8) leads to an ordinary differential equation for n s and a constant current density j s Inserting eq. (8) into eq. (6) results in where the purely stationary contribution to the source term was replaced by a flux derivative. This modified equation can now be discretized more easily in a well-balanced manner, because the remaining source term is small near the stationary solution. Note that the last term on the right-hand side of eq. (13) is not neglected, i.e. the dynamical part is not assumed to be small compared to the stationary part. Even though the source term in eq. (6) is linear, a similar procedure can be used if a more sophisticated model of the form is considered, where a non-linear source term T (U, ϕ) is included. The stationary equation is then given by Subtracting eq. (15) from eq. (14) yields corresponding to eq. (13).

Discretization
For simplicity we use an equidistant grid in x and y direction with the coordinates for the direct and for the adjoint grid points. Note that our scheme is not fundamentally limited to the equidistant case and could easily be extended to non-equidistant grids. The PE is discretized on the 2D grid in the usual way (vertex centered finite volumes [28]), resulting in discrete values for the potential defined on the direct grid nodes.

Stationary Problem
To discretize the stationary DD model eq. (8) we use the point-wise defined values U s located at the adjoint grid points as the solution variables of the discrete equation system. Starting from an initial guess we calculate time derivatives which are then brought to zero by an iterative process. At first the function U s (x) is reconstructed from the solution variables by solving the stationary equation eq. (8) as an initial value problem separately in each interval starting from its midpoint and going in positive and negative x direction. This computationally costly process, which can be realized by using a high-order stiff ODE solver, is necessary to enforce well-balancedness. The accuracy should be set to almost machine precision, in order to enable convergence for the whole system later. Another requirement for the ODE solver is the supply of the derivatives of the end state with respect to the initial state and the electric field, which is needed for the Newton iteration later. Here we used the fifth order RADAU IIA method described in [31].
The result is a piece-wise defined reconstruction function and The source matrix S s k+ 1 2 depends on the electric field which is calculated by finite-differences from ϕ s and is assumed to be constant in each interval [x k , x k+1 ]. Due to this fact the local solutionsŨ , which is important for the efficiency and accuracy of the ODE solver.
Next, the equation system including the time derivative is integrated over one interval resulting in where the values denoted by an asterisk are located at the direct grid points between two intervals. Since there are always two values (left and right) available at each grid point, the interface states have to be calculated by a suitable averaging process. For this purpose we use an exact Riemann solver (a simple derivation can be found e.g. in [22]) which takes the full upwind information into account. The computational effort of an exact solver is quite limited in this case of only two equations, however for more advanced models an approximate solver can be used as well. The resulting upwinded interface states are then shared between two adjacent intervals, ensuring charge conservation. At the contact grid points the Dirichlet condition for the electron density leads to a semi-Riemann problem on the contact surface, since the state directly inside of the device and the density directly on the contact surface are given. A Riemann-Solver for the DD model works by connecting two initial states by two waves, where in between the two waves a constant interface state forms. In the subsonic case the two waves will travel in opposite directions. By fixing the density for the interface state to n ct , we can then calculate the missing current density on the contact surface by connecting the state on the contact to the state directly inside of the device by the respective inflowing wave. For more advanced models, more information on the contact boundaries must be supplied and more inflowing waves considered, in general one boundary condition is needed for every two balance equations.
Using the stationary reconstruction, the source term in eq. (23) can be replaced by a flux difference, resulting in where Note that even though it may be tempting to use eq. (25) for the dynamic case as well, this is not reasonable, because the reconstruction process explicitly involves solving the stationary equation for each interval. Even if we ignore the high computational cost of the ODE solver, it is not necessarily an accurate reconstruction for the dynamic case. Self-Consistent coupling to the PE is achieved by extracting the densities at the direct grid points from the interface states and using them for the right-hand side of the PE. The stationary coupled PE-DD system is then written in residual form where the residuals f s U of the DD model are the time derivatives of eq. (25). x s ϕ and x s U denote the discrete solution variables. The equation system is then solved by a damped Newton iteration, where the Jacobian is a large, sparse matrix containing the derivatives of the residuals with respect to the solution variables. These derivatives are computed as part of the solution process.

Dynamic Problem
After solving the stationary system, its solution is fixed and we can move on to the dynamic system. As opposed to the stationary case, for the dynamic case we use the average values of U d in each vertex centered cell as the solution variables, because it leads to an easier coupling to the PE. We use a simple piece-wise constant reconstruction of the dynamic variables in each cell. As an alternative we tested the CWENOZ3 scheme [2,3,19] in combination with a local characteristic decomposition [25]. Due to its third-order accuracy in smooth regions, it should be preferable to the first-order piece-wise constant reconstruction. Unfortunately, however, the resulting scheme proved to be numerically unstable in the regions with strong source terms, since the local characteristic decomposition only takes the flux Jacobian into account.
Integrating eq. (13) over one cell yields where are averaged over half of a cell. For the piece-wise constant reconstruction the half cell averages are equal to the full cell averages. The variables denoted by an asterisk are as before solutions to Riemann problems, however this time located at the adjoint grid points. At the boundaries we use the values of U directly on the contact vertices as solution variables instead of cell averages. The reconstruction in the contact (half) cells is done by a local polynomial reconstruction satisfying leading to linear and quadratic reconstruction polynomials forñ d where ct and ct ± 1 denote the contact and the adjacent cell.
Since the density is fixed on the contacts, its time derivative is zero. The time derivative of the current density however is not and it is obtained by evaluating eq. (13) directly on the contact, without integration over a finite volume. The spatial derivatives in that equation are evaluated consistently with the stationary and dynamic reconstructions. The full dynamic system can then be written in the form of where D U is a diagonal matrix. This system can then be solved by standard time-integration schemes, like the BDF methods [5] or in the frequency domain by harmonic balance [9].

Results
We consider the device shown in Fig. 1, with a length of L = 100nm, where the gated length is  Driver circuit for resistive mixing. [15,24] respectively (see Fig. 2). Between each region of constant doping an exponential transition is assumed. We use an equidistant tensor-product grid with N x = 201 by N y = 21 grid points.
In Fig. 2 we tested our scheme by performing stationary simulations with two different applied drain voltages. The gate voltage was set to zero V G = 0V resulting in strong builtin fields to test the numerical stability of the scheme. For equilibrium conditions, where V S = V D = 0V, the electron density n eq roughly follows the doping concentration. In addition to this we calculated the density for V DS = 40mV by solving the stationary equation system eq. (27) directly as well as solving the dynamic system eq. (35) for vanishing time derivatives, where the fixed stationary solution was set to the equilibrium solution. The resulting density profiles are almost identical, confirming the validity of our approach of splitting the variables into stationary and dynamic parts.
For further analysis of our scheme we use the same driver circuit as in [15,24], shown in Fig. 3. Due to the high admittance at THz frequencies of the coupling capacitance between the gate and drain terminal, we make the same simplification as in [14,15] and apply the AC voltage to the source terminal instead, while keeping the gate and drain voltages constant (AC-wise grounded).  We start by calculating the small-signal drain-source admittance parameter Y DS , which is the drain current response to an AC voltage at source, where the device equations have been linearized around equilibrium with V DS = 0. To show the influence of the reconstruction procedure we compare the third-order CWENOZ3 scheme and the simple first-order piecewise constant reconstruction in Fig. 4. It is clearly visible, that the piecewise constant reconstruction introduces strong numerical damping at higher frequencies compared to the CWENOZ3 reconstruction. The stability problems we encountered for the high order reconstruction can not be seen in the small-signal regime, where the dynamic variables U d are treated as infinitesimal quantities, leading to a perfectly smooth solution. The non-linear weights introduced in the CWENO procedure here simply reduce to the linear weights of an optimal third-order parabolic reconstruction. By increasing the number of grid points the numerical damping of  Going further, we performed actual transient simulations of the device by using the secondorder BDF2-method [5]. The time-step size was chosen automatically by an adaptive timestep algorithm, which does not require additional information about the equation system, e.g. characteristic time scales (black-box approach). In Fig. 5    non-linear behavior of the device becomes visible in the asymmetric deformation of the drain current signal, which can not be captured by a simple small-signal analysis.
In Fig. 6 we show the internal variables of the periodic steady-state resolved over one period for a frequency of f = 1THz. In Fig. 6a the current density on the cell interfaces, namely the current density component of the Riemann solutions, is shown. These current densities are also used in the calculation of the terminal currents. In the highly doped regions it is nearly constant in space, leading to a constant electron density in time. In the lowly doped channel, the propagation of waves can be observed, leading to a phase shift of the current density at the source and drain terminals, which could already be seen in the drain-source admittance parameters. In Fig. 6b the total electron density at the cell centers normalized to the stationary density is shown. As expected, in the contact regions the density is essentially constant, and only in the channel oscillations occur.
We also calculated the current responsivity of the mixer given by the DC component of the current at the drain divided by the AC input power at the source terminal [15] averaged over one period T . In Fig. 7 the current responsivity at a DC bias of V DS = 0V for a frequency range of 0 to 3THz is shown. To highlight the effect the convective derivative has, we repeated the calculation without it, where both of these simulations were performed using the harmonic balance approach [9]. Even though in the small-signal regime the convective derivative has no effect at equilibrium, in this case the difference is quite obvious, highlighting the nonlinear nature of the responsivity effect. For small frequencies both models yield the same result, however for larger frequencies in the THz region inclusion of the convective term leads to much larger responsivities.

Conclusion
We developed a stable numerical method for the simulation of one dimensional balance equations with strong source terms. The equations were split into a stationary and a dynamic part, where well-balancedness was guaranteed for the stationary part. However, even for dynamic conditions different from the stationary solution, the scheme provides good results. We compared the results obtained from a low-order piecewise constant and a high order CWENOZ3 reconstruction and showed the damping effect the low-order reconstruction has at high fre-quencies. Furthermore, we demonstrated the capabilities of our scheme to perform nonlinear transient simulations which are able to capture important effects like the responsivity of a MOSFET mixer.
Funding Open Access funding enabled and organized by Projekt DEAL.
Availability of data and materials The datasets generated during and/or analysed during the current study are available from the corresponding author on reasonable request.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.