An adaptive multilevel Monte Carlo algorithm for the stochastic drift-diffusion-Poisson system

TL;DR

This paper introduces an adaptive multilevel Monte Carlo method for efficiently solving the stochastic drift-diffusion-Poisson system, combining goal-oriented mesh refinement and stochastic sampling to improve accuracy and computational efficiency.

Contribution

The paper develops a novel adaptive multilevel Monte Carlo algorithm with error estimation techniques for the stochastic drift-diffusion-Poisson system, enhancing convergence and efficiency over uniform refinement.

Findings

Achieves linear convergence of the H^1 error.

Improves error control through efficient error indicator estimation.

Demonstrates advantages over uniform mesh refinement in numerical examples.

Abstract

We present an adaptive multilevel Monte Carlo algorithm for solving the stochastic drift-diffusion-Poisson system with non-zero recombination rate. The a-posteriori error is estimated to enable goal-oriented adaptive mesh refinement for the spatial dimensions, while the a-priori error is estimated to guarantee \red{linear} convergence of the $H^1$ error. In the adaptive mesh refinement, efficient estimation of the error indicator gives rise to better error control. For the stochastic dimensions, we use the multilevel Monte Carlo method to solve this system of stochastic partial differential equations. Finally, the advantage of the technique developed here compared to uniform mesh refinement is discussed using a realistic numerical example.

Tables

Table 1. Table 1: Degrees of freedom 𝒩 𝒫 ℓ subscript 𝒩 subscript 𝒫 ℓ \mathcal{N}_{\mathcal{P}_{\ell}} for different levels comparing uniform MLMC-FE and adaptive MLMC-FE methods.

$\ell$	0	1	2	3	4	5	6
uniform	1 685	4 749	11 626	27 084	60 569	131 819	280 264
adaptive	1 685	2 839	5 531	10 876	22 128	45 885	98 123

Table 2. Table 2: Degrees of freedom 𝒩 𝒟 ​ 𝒟 ℓ subscript 𝒩 𝒟 subscript 𝒟 ℓ \mathcal{N}_{\mathcal{DD}_{\ell}} for different levels comparing uniform MLMC-FE and adaptive MLMC-FE methods.

$\ell$	0	1	2	3	4	5	6
uniform	755	2 141	5 293	12 372	27 743	60 509	128 833
adaptive	755	898	1 640	3 166	6 614	13 553	29 305

Table 3. Table 3: The optimal number of samples for the uniform MLMC-FE method.

$\epsilon$	$M_0$	$M_1$	$M_2$	$M_3$	$M_4$	$M_5$
0.080	17	5	–	–	–	–
0.040	90	26	8	–	–	–
0.020	418	121	35	11	–	–
0.010	1 671	481	137	42	–	–
0.007	3 759	1 080	308	94	30	–
0.005	7 366	2 117	604	183	58	–
0.002	5 3764	17 113	4 163	1 012	573	57

Table 4. Table 4: The optimal number of samples for the adaptive MLMC-FE method.

$\epsilon$	$M_0$	$M_1$	$M_2$	$M_3$	$M_4$	$M_5$	$M_6$
0.080	18	6	3	–	–	–	–
0.040	78	26	10	4	–	–	–
0.020	334	111	44	14	5	–	–
0.010	1 369	472	186	41	20	–	–
0.007	2 856	949	372	119	42	14	–
0.005	5 597	1 859	729	233	81	27	–
0.002	36 068	11 979	4 692	1 502	520	170	54

Equations

$$-\nabla \cdot (A(x,\omega )\nabla V(x,\omega ))=\rho (x,\omega ),$$

$$\rho =q(p-n+C).$$

$$n(x,\omega )$$

$$p(x,\omega )$$

$$D={\text{D}}_{\text{Si}}\cup {\text{D}}_{\text{ox}}.$$

$$V{|}_{\partial D_D}=V_D\text{and}\mathbf{n}\cdot \nabla V{|}_{\partial D_N}=0$$

$$\partial V(x)\cdot \mathbf{n}$$

$$J_n(x)\cdot \mathbf{n}$$

$$J_p(x)\cdot \mathbf{n}$$

$$-\nabla \cdot (A(x,\omega )\nabla V(x,\omega ))$$

$$\nabla \cdot J_n(x,\omega )$$

$$-\nabla \cdot J_p(x,\omega )$$

$$J_n(x,\omega )$$

$$J_p(x,\omega )$$

$$R_{\text{SRH}}(n,p):=\frac{np-n_i^2}{{\tau }_n(p+n_i)+{\tau }_p(n+n_i)},$$

$$-\nabla \cdot (A(x)\nabla V(x,\omega ))$$

$$U_T{n}_i\nabla \cdot ({\mu }_n{e}^{V/U_T}\nabla u(x,\omega ))$$

$$U_T{n}_i\nabla \cdot ({\mu }_p{e}^{-V/U_T}\nabla v(x,\omega ))$$

$$R_{\text{SRH}}(x,\omega )=n_i\frac{u(x,\omega )v(x,\omega )-1}{{\tau }_p(e^{V/U_T}u(x,\omega )+1)+{\tau }_n(e^{-V/U_T}v(x,\omega )+1)}.$$

$$u(x,\omega ){|}_{\partial D_{D,\mathrm{Si}}}=u_D(x)\text{and}v(x,\omega ){|}_{\partial D_{D,\mathrm{Si}}}=v_D(x).$$

$$V(0+,y,\omega )-V(0-,y,\omega )$$

$$A(0+){\partial }_{x}V(0+,y,\omega )-A(0-){\partial }_{x}V(0-,y,\omega )$$

$$-\nabla \cdot (A(x)\nabla V(x,\omega ))$$

$$U_T\nabla \cdot \left({\mu }_n{e}^{\frac{V\left(x,\omega \right)}{U_T}}\nabla u\left(x,\omega \right)\right)$$

$$U_T\nabla \cdot \left({\mu }_p{e}^{-\frac{V\left(x,\omega \right)}{U_T}}\nabla v\left(x,\omega \right)\right)$$

$$V(0+,y,\omega )-V(0-,y,\omega )$$

$$A(0+){\partial }_{x}V(0+,y,\omega )-A(0-){\partial }_{x}V(0-,y,\omega )$$

$$V(x,\omega )$$

$$\mathbf{n}\cdot \nabla V(x,\omega )$$

$$u(x,\omega )=u_D(x),v(x,\omega )$$

$$\mathbf{n}\cdot \nabla u(x,\omega )=0,\mathbf{n}\cdot \nabla v(x,\omega )$$

$$\underline{C}:=\underset{x\in D}{\inf }{C}_{\mathrm{dop}}(\cdot ,\omega )\le C_{\mathrm{dop}}(\cdot ,\omega )\le \underset{x\in D}{\sup }{C}_{\mathrm{dop}}(\cdot ,\omega )=:\overline{C}\forall \omega \in \Omega .$$

$$\frac{1}{K}\le u_D(x),v_D(x)\le K\forall x\in \partial D_{\mathrm{Si},D}.$$

$$0<{\mu }_n^{-}$$

$$0<{\mu }_p^{-}$$

$${\left(u,v\right)}_D:={\int }_{D}u(x)v(x)dx\forall u,v\in L^2(D)$$

$${\Vert u\Vert }_D={\left({\int }_{D}u(x{)}^{2}dx\right)}^{\frac{1}{2}}\forall u\in L^2(D).$$

$${\left(u,v\right)}_{k,D}:=\sum _{0\le |m|\le k}\left(\underset{D}{\int }{D}^{m}u(x)D^{m}v(x)dx\right)$$

$${\Vert u\Vert }_{H^k(D)}:={\left(\sum _{|\alpha |\le k}{\Vert D^{\alpha }u\Vert }_{L^2\left(D\right)}^2\right)}^{\frac{1}{2}}$$

Full text

AMS subject classifications: 60H15, 35R60, 65M50, 82D37

1 Introduction

The stochastic drift-diffusion-Poisson (DDP) system is a general model for charge transport in random environments. A leading example is the field-effect transistor (FET), where the stochastic coefficients can describe process variations, noise, and fluctuations in devices as diverse as transistors and sensors. Process variations, noise, and fluctuations are significantly important especially in devices scaled into the deca-nanometer regime, as random effects become more important in smaller devices. Among the many sources of noise, random-dopant fluctuations (RDF) [1, 2, 3] are one of the most important. Random-dopant fluctuations stem from the fact that the doping process in the semiconductors leads to a random number and random position of dopants. Therefore, each impurity atom influences the charge transport and the mobilities. A schematic diagram is shown in Figure 1.

In this paper we analyze a standard adaptive finite element method of the form SOLVE $\rightarrow$ ESTIMATE $\rightarrow$ MARK $\rightarrow$ REFINE [4] to discretize the stochastic drift-diffusion-Poisson. In order to do the MARK step we need to compute the a posteriori error estimate for each element.

A-priori error estimates yield knowledge about convergence and stability of the solvers and information on the asymptotic behavior of errors for different mesh sizes [5]. A-posteriori error estimates make it possible to control the mesh on the entire computational domain by using adaptive algorithms, i.e., by focusing computational effort on the parts of the domain which contribute most to the total error [6]. In adaptive mesh refinement, a-posteriori error estimators are used to indicate where the error is particularly high, and then more mesh elements are placed in those locations. Here we estimate the local error for a coupled system of equations. The error estimate indicates which elements should be refined or coarsened simultaneously for the Poisson equation and the drift-diffusion equations.

The mentioned stochastic problem is computationally very expensive; in order to obtain an acceptable error, thousands of simulations are necessary. In stochastic PDE, the Monte Carlo (MC) method is one of the most popular and straightforward numerical techniques. However, the main drawback of the MC method is its well-known convergence rate. The multilevel Monte Carlo finite-element (MLMC-FE) method [7, 8, 9] is an efficient numerical alternative. In [10], we introduced an optimal MLMC-FE method to model the random effects in a charge-transport model. In the optimal MLMC-FE method, in each level, we determine mesh sizes $h_{\ell}$ and the number of samples $M_{\ell}$ to minimize the computational costs such that the total error (i.e., statistical and discretization error) is less than a prescribed error tolerance. In [11], the efficiency of the method for three-dimensional simulations of various nanoscale devices was investigated in detail. Convergence can be improved by using a randomized rank-1 lattice rule [12, 13].

The first analysis of a finite-element method, a one-dimensional one, for solving the (deterministic) DDP system can be found in [14]. An extension of the analysis to the two-dimensional problem was presented in [15]. In [16], fixed points of finite-element discretizations were used to approximate the solutions of the steady-state drift-diffusion system, and the convergence rate in the energy norm was estimated. In [17], the optimal convergence rate and its stability were shown. However, in all these publications, the recombination rate is zero.

In [18], an adaptive stochastic Galerkin finite-element method for linear elliptic boundary-value problems was presented. The idea of using an adaptive MLMC method for weak approximations of solutions of stochastic differential equations was explained in [19]. In [20], an adaptive MLMC algorithm was introduced for PDEs with stochastic data.

In [10], we showed the effectiveness of using an MLMC method as an alternative to the MC method in the solving the stochastic DDP system. There the meshes were refined uniformly. The novelty of the present work is computing an accurate local-error estimator based on goal oriented error estimation for the coupled system of equations. In order to control the error, only a part of the computational geometry which has a large error must be refined and therefore, the method is more computationally effective. Randomness in the model problem considered here stems from the random position of dopants in a transistor, which affect the error and hence the refinement process. Since MLMC is a variance-reduction method, the faster decay in the variance of the MLMC method leads to a reduction of the statistical error and therefore the total computational cost. Here the effect of adaptive mesh refinement on the reduction of the variance is also taken into account. Finally, the developed numerical technique can be used for other stochastic problems as well, e.g., [21, 22, 23, 24, 25, 26].

The remainder of this manuscript is as follows. In Section 2, we give the system of model equations with stochastic coefficients and explain its boundary consitions. In Section 3, we formulate a finite element method for the space discretization of the equation. In Section 4, we present the error estimates in the finite-element space, namely an a-priori error estimate and an a-posteriori error estimate. In Section 5, we introduce the MLMC finite-element method for the DDP system and define an optimization problem to minimize the total computational cost. Then, we present numerical results for a transistor and quantify the random-dopants effect in Section 6. The adaptive MLMC-FE method is used to approximate the expected value of the solution of the system of equations with random coefficients. Also, the method is compared with MLMC-FE method with uniform mesh refinement. Finally, conclusions are drawn in Section 7.

2 The Stochastic Model Problem

The stochastic Poisson equation is used generally for the electrostatic potential

[TABLE]

on the bounded Lipschitz domain $D\subset\mathbb{R}^{2}$. In the equation, $V$ indicates the electrostatic potential, $A$ denotes the dielectric constant (permittivity), and $\rho$ is the charge concentration. In (1), $x\in D$ is the spatial variable, $\omega=(\omega_{1},\omega_{2},\ldots,\omega_{n})$ is an $n-$dimensional random variable defined on the complete probability space $(\Omega,\mathbb{A},\mathbb{P})$ equipped with $\mathbb{A}\subset 2^{\Omega}$ as the $\sigma-$algebra of events, $\mathbb{P}\colon\mathbb{A}\rightarrow[0,1]$ as a probability measure, the sample space $\Omega$. The randomness arises from the random distribution of dopant atoms (uniformly distributed) in source and drain areas (shown in Figure 1). Also, the coefficient $A\colon D\times\Omega\rightarrow\mathbb{R}^{2\times 2}$ is a matrix of real-valued functions, and $\rho$ is a real-valued scalar function. In a semiconductor, the charge concentration is derived by the free electron and hole densities (i.e., $n$ and $p$) and the doping concentration $C$; the total charge concentration is therefore

[TABLE]

We change the concentrations $n$ and $p$ to the so-called Slotboom variables $u$ and $v$, which are given by

[TABLE]

Here, $n_{i}$ is the intrinsic carrier density of the semiconductor (in the numerical examples, a value of $1.5\times 10^{10}\,\text{cm}^{-3}$ is used for silicon) and $U_{T}$ indicates the thermal voltage, which is at room temperature is about $26\,\mathrm{mV}$.

A schematic diagram of a sample computational geometry is shown in Figure 1. The domain is partitioned into two subdomains, i.e.,

[TABLE]

The first subdomain ${\text{D}_{\text{Si}}}$ consists of silicon, i.e., the channel and the source and drain areas, where the drift-diffusion-Poisson system models the charge transport. The gate contact is separated from the channel by an insulating silicon dioxide layer $\text{D}_{\text{ox}}$. In $\text{D}_{\text{ox}}$, there is no charge transport and therefore we only have the Poisson equation.

The boundary $\partial D$ of the domain $D$ is separated into $\partial D_{D}$ and $\partial D_{N}$, which denote the surfaces where the Dirichlet and Neumann boundary conditions

[TABLE]

hold, where $\bm{n}$ is the outward pointing unit normal on the boundary $\partial D_{N}$. Dirichlet boundary conditions $V|_{\partial D_{D}}$ are applied for the potential at the source, drain, and gate contacts, i.e., $V=V_{\text{S}}$, $V=V_{\text{D}}$, and $V=V_{\text{G}}$.

Neumann boundary conditions are applied on all other boundaries. On the Neumann parts $\partial D_{N}$ of the boundary the currents and the electric field are assumed to vanish in the normal direction to the surface. This yields the three Neumann boundary conditions

[TABLE]

The stochastic DDP system

[TABLE]

is employed to model self-consistent charge transport, where $J_{n}$ and $J_{p}$ are the current densities, $\mu_{n}$ and $\mu_{p}$ are the mobilities, and $D_{n}$ and $D_{p}$ are diffusion coefficients, which can be calculated by the Einstein relations $D_{p}=U_{T}\mu_{p}$ and $D_{n}=U_{T}\mu_{n}$.

Furthermore, $R$ is Shockley-Read-Hall recombination rate, which is defined by

[TABLE]

where $\tau_{n}$ and $\tau_{p}$ are the lifetimes of the free carriers (absolutely positive). For the purposes of the present work, any other recombination rate can be used as long as it satisfies modest assumptions [10].

Using the Slotboom variables $u$ and $v$ defined in (2), the DDP system (4) takes the form

[TABLE]

with the Shockley-Read-Hall recombination rate

[TABLE]

The Dirichlet boundary conditions for the Slotboom variables are

[TABLE]

Zero Neumann boundary conditions hold on the Neumann part $\partial D_{N}$ of the boundary. The interface conditions

[TABLE]

can be used to model the presence of a layer of charge carriers at the surface of a FET after homogenization [27]. Here $\Gamma$ is the interface or surface between ${\text{D}_{\text{Si}}}$ and $\text{D}_{\text{ox}}$, and the notation $0+$ and $0-$ denotes the limits from both sides of the interface $\Gamma$ located at $x=0$. The directions of $y$ are along the interface. The model has been used to model charge transport [28] in modern nanoscale devices, see, e.g., [29, 30, 31, 32]

Finally, for all $\omega\in\Omega$, we can write the following boundary-value problem (BVP)

[TABLE]

In order to state the main result, the coefficients and boundary conditions in the boundary-value problem (7) have to satisfy the following assumptions.

Assumptions 1**.**

The bounded computational domain $D\subset\mathbb{R}^{2}$ has a $C^{2}$ Dirichlet boundary $\partial D_{D}$, the Neumann boundary $\partial D_{N}$ consists of $C^{2}$ segments, and the Lebesgue measure of the Dirichlet boundary $\partial D_{D}$ which is nonzero. The $C^{2}$ manifold $\Gamma\subset D$ separates the domain $D$ into two nonempty regions $D^{+}$ and $D^{-}$; therefore, $\mathop{meas}(\Gamma\cap\partial D)=0$ and $\Gamma\cap\partial D\subset\partial D_{N}$ hold. 2. 2.

*The coefficient $A(x,\cdot)$ is assumed to be a strongly measurable mapping from $\Omega$ into $L^{\infty}(\Omega)$. Also, the coefficient $A$ is *symmetric, satisfying the ellipticity condition, and bounded w.r.t. position $x\in D$ and elementary events $\omega\in\Omega$. Moreover, {\color[rgb]{0,0,0}A}|_{D^{+}\times\Omega}\in C^{1}(D^{+}\times\Omega,\mathbb{R}^{2\times 2}) and {\color[rgb]{0,0,0}A}|_{D^{-}\times\Omega}\in C^{1}(D^{-}\times\Omega,\mathbb{R}^{2\times 2}). 3. 3.

The doping concentration $C_{\mathrm{dop}}(x,\omega)$ is bounded above and below with the bounds

[TABLE] 4. 4.

There is a constant $K\geq 1$ which satisfies

[TABLE] 5. 5.

The electron and hole mobilities are uniformly bounded functions of $x\in D$ and $\omega\in\Omega$. Therefore, $\forall x\in D$ as well as $\forall\omega\in\Omega$ we have

[TABLE]

*where $\mu_{p},\mu_{n}\in C^{1}({\text{D}_{\text{Si}}}\times\Omega,\mathbb{R}^{2\times 2})$.

Moreover, the inclusions $V_{D}\in H^{1/2}(\partial D)\cap L^{\infty}(\Gamma)$ and $u_{D},v_{D}\in H^{1/2}(\partial{\text{D}_{\text{Si}}})$ hold.*

Remark 1. For all $\omega\in\Omega$, the matrix $A(\cdot,\omega)$ is uniformly elliptic, i.e., there is {\color[rgb]{0,0,0}\hat{A}(\cdot,\omega)\in C^{1}(D,\mathbb{R}^{2\times 2})} such that $A(\cdot,\omega)=\hat{A}(\cdot,\omega)*\hat{A}(\cdot,\omega)$.

3 Formulation of the finite element method

This section is devoted to present a finite element method for solving (7) and separated into two subsections. In the first one, we obtain a variational formulation and provide some results regarding the existence and uniqueness of the weak solutions. In the second subsection, the finite element discretization based on the mentioned variational formulation is presented.

3.1 The random variable dependent weak formulation

Assume $D\subset\mathbb{R}^{2}$ satisfies Assumption 1 . To present the variational formulation of Equation (7), we consider the Hilbert space $L^{2}\left(D\right)$ equipped with the inner product

[TABLE]

and hence the norm

[TABLE]

Also, we denote the Hilbert space equipped with the inner product

[TABLE]

and the norm

[TABLE]

by $H^{k}(D)$.

Moreover, we define

[TABLE]

The random variable dependent weak formulation of (7a)–(7d), i.e., the Poisson equation on $D$ and the drift-diffusion equations for electrons and holes on ${\text{D}_{\text{Si}}}$ for a fixed random variable $\omega\in\Omega$ can be written as

[TABLE]

We mention the following lemma to state existence and uniqueness of the continuous solution.

Lemma 1**.**

[10, 27]** Considering Assumptions 1, for $V_{D}\in H^{1/2}(\partial D)$ and $u_{D},v_{D}\in H^{1/2}(\partial D_{Si})$, there exists a unique random variable-dependent weak solution

[TABLE]

for Eq. (8)

3.2 The random variable dependent finite element formulation

Assume $D$ is triangulated by a regular (in the sense of Ciarlet), locally quasi-uniform mesh $\mathcal{T}_{\ell}$ with mesh width $h_{\ell}$. Also, we assume that the mesh is $\gamma$-shape regular in the sense that $\text{diam}(T)\leq\gamma|T|^{1/2}$ for all $T\in\mathcal{T}_{\ell}$. Let $\mathcal{P}^{1}(T)$ be the space of first degree polynomials on $T\in\mathcal{T}_{\ell}$. Then, we define

[TABLE]

Hence, based on the weak formulation (8), we arrive at the random variable dependent Galerkin discretization of finding $(V_{\ell}(\cdot,\omega),u_{\ell}(\cdot,\omega),v_{\ell}(\cdot,\omega))\in(X_{1,\ell},X_{2,\ell},X_{3,\ell})$ such that

[TABLE]

The existence and uniquness of a FEM formulation of (10) is shown in [33].

4 A priori and a posteriori error estimation

In this section, we derive both a-priori and a-posteriori error estimates. The main feature of an a-priori error estimate is providing knowledge about the asymptotic behavior of the discretization error. Considering an a-posteriori error estimator, an adaptive mesh-refinement process consists of the following steps:

Define an initial mesh $\mathcal{T}_{\ell}$. Let $\ell=0$. 2. 2.

Solve the discrete problem (10) on $\mathcal{T}_{0}$. 3. 3.

For each element $T\in\mathcal{T}_{\ell}$, obtain the computed error estimate. 4. 4.

If the computed global error is sufficiently small, then stop. Otherwise, determine which elements have to be refined and obtain the next mesh $\ell+1$ and return to the second step (2).

For the sake of simplicity, from here we use $h$ instead of $h_{\ell}$ to denote the mesh width.

Theorem 2** (A-priori error estimate).**

Let $(V(\omega),u(\omega),v(\omega))\in\left(H^{2}(D\setminus\Gamma)\cap L^{\infty}(D\setminus\Gamma)\right)\times\left(H^{2}({\text{D}_{\text{Si}}})\cap L^{\infty}({\text{D}_{\text{Si}}})\right)^{2}$ be the solution of (8) and $(V_{\ell}(\omega),u_{\ell}(\omega),v_{\ell}(\omega))\in X_{1,\ell}\times X_{2,\ell}\times X_{3,\ell}$ be the solution of (10). Then, there exists a constant $c\in\mathbb{R}^{+}$ depending on the doping concentration as well as the shape regularity of the mesh such that the inequality

[TABLE]

holds for a fixed random variable $\omega\in\Omega$.

Proof.

For the sake of simplicity, the proof is done in four steps.

Step 1: In this step, we provide four equations containing the error terms. Let

[TABLE]

where $\Pi_{1,\ell}$ is the $L^{2}$ projector of $L^{2}(D)$ onto $X_{1,\ell}$ and $\Pi_{i,\ell}$ projects $L^{2}(D_{Si})$ onto $X_{i,\ell}$ for $i=2,3$. To simplify the notations, we drop the random variable $\omega$ from the arguments of all functions in the following.

Substituting $\varphi_{1}=e_{1,\ell}$ into (8a) and $\psi_{1}=e_{1,\ell}$ into (10a) as well as subtracting these two equations, we have the following error equation

[TABLE]

Substituting $\varphi_{1}=e_{1,\ell}$ into (8b) and $\psi_{1}=e_{1,\ell}$ into (10b) as well as subtracting these two equations results in

[TABLE]

Furthermore, substituting $\varphi_{2}=e_{2,\ell}$ into (8c) and $\psi_{2}=e_{2,\ell}$ into (10c) and subtracting these two equations leads to

[TABLE]

Again, substituting $\varphi_{3}=e_{3,\ell}$ into (8d) and $\psi_{3}=e_{3,\ell}$ into (10d) and subtracting these two equations yields

[TABLE]

Step 2: In this step, we simplify some terms of the error equations from the previous step. For this purpose, we consider the following notations

[TABLE]

Then, using the notations from (14)-(16), we can write

[TABLE]

and

[TABLE]

Also, It can be easily seen

[TABLE]

and

[TABLE]

Moreover, for the sake of simplicity we then define

[TABLE]

Substituting (4)–(19) into (4)–(13), together with the above notation we infer that

[TABLE]

Step 3: In this step, we strive to find an upper bound for the terms given in the right hand side of (4).

We denote the functional $M_{1}:X_{1,\ell}\times X_{1,\ell}\times X_{3,\ell}\rightarrow\mathbb{R}$ by

[TABLE]

Moreover, we define

[TABLE]

where $0<\theta<1$. Using the Taylor expansion for the functional $M_{1}$ over $\Pi_{1,\ell}V$ leads to

[TABLE]

Since the $F$ is differentiable, considering the directional derivative of $M_{1}$ in the direction of projection of $V$ leads to

[TABLE]

Here, we should note that $F$ is Lipschitz continuous (see Eq. (3.13) of [33]) and the constant $c(\theta_{1})$ is independent of the mesh width. Considering (22) and (4), we conclude that

[TABLE]

Therefore, Eq. (4) and the triangle inequality give rise to

[TABLE]

In the next step, we define the functional

[TABLE]

Similarly, we define the variable

[TABLE]

where $0<\theta<1$. Taylor expansion of $M_{2}$ over $\Pi_{2,\ell}u$ shows that

[TABLE]

Again, since the function $f$ is differentiable, the directional derivative of $M_{2}$ is given by

[TABLE]

where again the constant $c(\theta_{2})$ is independent of $h$. Therefore, we have

[TABLE]

and applying the triangle inequality to (25) results in

[TABLE]

Now using the triangle inequality and applying the the approximation and the stability properties of $\Pi_{1,\ell}$, $\Pi_{2,\ell}$, and $\Pi_{3,\ell}$, we have

[TABLE]

as well as

[TABLE]

and

[TABLE]

Moreover, using the approximation and stability properties of $\Pi_{i,\ell}u$, i=1,2,3, we can write

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Also, we can get the similar results for $\left|\left(e^{\frac{\Pi_{1,\ell}V}{U_{T}}}\Pi_{2,\ell}u-e^{\frac{V}{U_{T}}}u,e_{1,\ell}\right)_{{\text{D}_{\text{Si}}}}\right|$, $\left|\left(e^{\frac{\Pi_{1,\ell}V}{U_{T}}}\nabla\Pi_{2,\ell}u-e^{\frac{V}{U_{T}}}\nabla u,\nabla e_{2,\ell}\right)_{{\text{D}_{\text{Si}}}}\right|$ as well as $\left|\left(e^{\frac{\Pi_{1,\ell}V}{U_{T}}}\Pi_{2,\ell}u-e^{\frac{V_{\ell}}{U_{T}}}u_{\ell},e_{1,\ell}\right)_{{\text{D}_{\text{Si}}}}\right|$.

Step 4: Combination of the above arguments gives us

[TABLE]

Finally, using the triangle inequality yields the desired results. ∎

In the next step, we use a residual based a-posteriori error estimation technique to estimate the local error $\eta_{T}$ on each finite element $T\in\mathcal{T}_{\ell}$. The error indicator will serve as the foundation for a refinement strategy in order to control and minimize the errors in the Poisson equation (7a–7b) and in the continuity equations (7c–7d).

[TABLE]

Theorem 3** (A-posteriori error estimate).**

For $\omega\in\Omega$ let $(V(\omega),u(\omega),v(\omega))\in\left(H^{1}(D\setminus\Gamma)\cap L^{\infty}(D\setminus\Gamma)\right)\times\left(H^{1}({\text{D}_{\text{Si}}})\cap L^{\infty}({\text{D}_{\text{Si}}})\right)^{2}$ be the solution of (8) and $(V_{\ell}(\omega),u_{\ell}(\omega),v_{\ell}(\omega))\in X_{1,\ell}\times X_{2,\ell}\times X_{3,\ell}$ be the solution of (10). Then, for sufficiently small $h$, there exist constants $c_{i}$, $i\in\{1,\ldots,6\}$, depending on the doping concentration as well as the shape regularity of the mesh such that

[TABLE]

holds for a fixed event $\omega\in\Omega$ where $\mathcal{T}^{\text{Si}}_{\ell}:=\mathcal{T}_{\ell}\cap{\text{D}_{\text{Si}}}$ and

[TABLE]

Here $\zeta$ denotes the area of an element in $\mathcal{T}_{\ell}$ or $\mathcal{T}^{\text{Si}}_{\ell}$, $\gamma$ denotes the boundary of the element, and the brackets $\left[\cdot\right]$ indicate the jump at the element boundary.

Proof.

In the following, the dependence of the solutions on the random variable $\omega$ is not indicated in order to simplify notation. We first define

[TABLE]

Using the test functions $\varphi_{1}\in X_{1}$, $\varphi_{2}\in X_{2}$, and $\varphi_{3}\in X_{3}$, the weak formulation (8) yields

[TABLE]

Let $I_{i,\ell}$ be the projection operator on $X_{i,\ell}$ defined in Theorem 4.8.7 of [34].

Next we substitute $\psi_{1}:=I_{1,\ell}{\varphi_{1}}$ into (10a) and (10b), $\psi_{2}:=I_{2,\ell}{\varphi_{2}}$ into (10c), and $\psi_{3}:=I_{3,\ell}{\varphi_{3}}$ into (10d), and sum up the equations. Then we subtract them from (33)–(33d), which leads to

[TABLE]

Explicit error estimation involves the direct computation of the interior element residuals and the jumps at the element boundaries to find an estimate for the error. Now using Green’s theorem and the Cauchy-Schwarz inequality on (34) leads to

[TABLE]

where

[TABLE]

It can be easily seen that

[TABLE]

Now substituting $\varphi_{i}={\varepsilon_{i,\ell}}$ for $i\in\{1,2,3\}$ into (4) and using (36a–36d), Eqs. (24), (27), [34, Corollary 4.8.15 on page 123], and Cauchy-Schwarz inequality yield

[TABLE]

where $\mathcal{T}^{\text{Si}}_{\ell}:=\mathcal{T}_{\ell}\cap{\text{D}_{\text{Si}}}$. Finally, applying an inverse estimate, for sufficiently small $h$ we get

[TABLE]

which concludes the proof. ∎

5 Multilevel Monte Carlo Finite-Element Method

In a Monte Carlo finite-element (FE) method several evaluations are combined to obtain an approximation of the solution of the model equation or equations. Considering the Bochner space $L^{2}(\Omega;X)$ for the mappings $\mathcal{Y}\colon\Omega\to X$, we can define

[TABLE]

The standard MC estimator $\operatorname{E_{MC}}$ for $\mathbb{E}[u_{\ell}]$ is the sample mean

[TABLE]

where $u_{\ell}^{(i)}=u_{\ell}(x,\omega^{(i)})$ is the $i$-th sample (independent random variable) of the solution $u$.

In a multilevel Monte Carlo (MLMC) method, instead of calculating the expected value $\mathbb{E}[u]$ by $\mathbb{E}[u_{\ell}]$ on a constant triangulation $\mathcal{T}_{\ell}$, the MLMC method approximates the expected value $\mathbb{E}[u]$ using several $\mathbb{E}[u_{\ell}]$, $\ell\in\{0,1,\ldots,L\}$, estimated on the nested family $\{\mathcal{T}_{{\ell}}\}_{{\ell}=0}^{\infty}$. In fact, to overcome the drawback of the MC method, the MLMC estimator avoids prohibitively many expensive evaluations of $\mathbb{E}[u_{\ell}]$ on the finest level $L$.

The FE approximation of the expected value of $u_{L}$ at level $L$ can be written as

[TABLE]

Therefore, for different numbers $M_{\ell}$ of samples at level $\ell\in\{0,\ldots,L\}$, we have

[TABLE]

The mean square error (MSE) is estimated by

[TABLE]

as shown for example in [10], where the variance is given by $\sigma[u]^{2}:=\|\mathbb{E}[u]-u\|_{L^{2}(\Omega;D)}^{2}$. The first and second terms of (42) are the statistical error, while the last term is the discretization error. In the next section, we study the effect of mesh refinement, i.e., uniform refinement and adaptive refinement (using the error indicator (4)) on both terms of errors.

Here, we strive to use the a priori and a posteriori error estimates obtained above for the stochastic drift-diffusion-Poisson system. To do so, we consider the convexity of the norms and employ Jensen’s inequality as well as the Poincaré inequality. Therefore, taking the expectation in Theorems 2 and 3 with respect to a random variable leads to the following corollaries.

Corollary 3.1**.**

The error bound

[TABLE]

holds for a priori error estimation.

Corollary 3.2**.**

The error bound

[TABLE]

holds for a posteriori error estimation.

Corollary 3.2 makes it possible to estimate the a posteriori error estimator for the expected value of the solutions. Marking strategies such as the Dörfler strategy [35] can now be used to drive mesh adaptivity.

In order to estimate the computational errors, we continue with the degrees of freedom. It enables us to draw a fair comparison between adaptive MLMC-FE and uniform MLMC-FE methods. According to the error bound given in Corollary 3.1, we define the discretization error as

[TABLE]

Furthermore, at level $\ell$, we assume that

[TABLE]

where $\mathcal{N}_{\mathcal{P}_{\ell}}$ is the number of unknowns or the (degrees of freedom) for the Poisson equation and $\mathcal{N}_{\mathcal{D}_{\ell}}$ indicates the number of unknowns for the two continuity equations. The exponent $\alpha$ is the convergence rate of the error. For the statistical error, the following inequality

[TABLE]

is assumed to show the convergence of the statistical error (at level $\ell$). For $\ell=0$, the assumption

[TABLE]

is used as well.

Due to the computational challenge of solving a system of SPDEs, an effective computational strategy is crucial. We strive to determine the optimal number $M_{\ell}$ of samples which minimize the computational work when $\operatorname{MSE}\leq\varepsilon^{2}$. In other words, the optimal number of samples are defined such that the statistical error is less than $\varepsilon^{2}/2$. The optimal value of $L$ (the lowest possible number) determined in the sense that the discretization error ($\mathcal{E}_{L}$) is less than $\varepsilon/\sqrt{2}$. For this, the following optimization problem is solved

[TABLE]

where the optimization is over all $M_{\ell}>0$. Moreover, since the optimal numbers $M_{\ell}$ of samples at level $\ell$ are in general not integers, they are rounded up and replaced by $\lceil M_{\ell}\rceil$. The details of the optimal approach were given in Giles’s MLMC paper [12].

6 Numerical Example and Results

In this work, we have chosen a double-gate MOSFET (DG-MOSFET) as a realistic example to implement the multilevel adaptive method developed above and to investigate its behavior. In these semiconductor transistors, the width of the silicon channel is very small and two gate contacts are used in the both sides of the channel to control the channel efficiently [36, 37, 38]. Hence, the current can potentially be twice the current through a single-gate device, since inversion layers can exist at both gates. This device structure suppresses short-channel effects and leads to higher currents as compared to the usual MOSFET structure having only one gate. In this work we assumed that there is no charge fluctuation in the oxide layer (see 7b). However, due to charges in the dielectric subdomain and its boundaries (interface to the channel region), the right-hand side in Eq. 7b can also be assumed to be non-zero. In digital and non-digital applications, the sensitivity of nanowire conduction to trapping and de-trapping of charges at interface states can be considered as well.

The FET device (see Figure 1 for a schematic diagram) consists of two materials, namely silicon (${\text{D}_{\text{Si}}}$) in the channel and source and drain regions and silicon dioxide ($\text{D}_{\text{ox}}$) as the insulator. The purpose of the insulator is to suppress direct charge flow from the gate into the channel and vice versa. The permittivities of the materials are $A_{\text{Si}}=11.7A_{0}$ and $A_{\text{ox}}=3.9A_{0}$, where the vacuum permittivity (dielectric constant) is $8.85\cdot 10^{-12}\,\mathrm{Fm^{-1}}$. Moreover, the gates have a length of $30\,\mathrm{nm}$ and are separated from the silicon channel by a $2\,\mathrm{nm}$ thick oxide layer. The channel width is $W=15\,\mathrm{nm}$ and it is connected to the heavily n-type doped source and drain regions of length $L_{\text{SD}}:=10\,\mathrm{nm}$ in each region.

Regarding the boundary conditions of the model equations, we apply Dirichlet boundary conditions at the gates ($V_{\text{g}}=0.2\,\mathrm{V}$) and, the source-to-drain voltage is $V_{\text{SD}}=0.1\,\mathrm{V}$. Also, the thermal voltage is $U_{T}=0.026\,\mathrm{V}$. The contacts are illustrated in Figure 1. For the rest of the transistor, we apply zero Neumann boundary conditions. In the highly doped source and drain regions, the dopant atoms are randomly distributed and indicated in the figure by blue circles. Therefore, random-dopant effects are included due to the random position of the dopant atoms.

In the common models, the doping concentration is modeled as a macroscopic, deterministic quantity which averages out any microscopic non-uniformities due to the random placement and random number of dopants. For large devices with a large number of dopants, this continuum model is physically reasonable, since the electrostatic potential appears spatially homogeneous and is sufficiently well described by the averaged charge density. However, in nanoscale transistors, the randomly distributed dopant atoms lead to inevitable variations between the billions of transistors in an integrated circuit.

As mentioned already, the main source of device variation is the random motion of impurity atoms during the fabrication procedure of implantation and annealing. In order to model the stochastic coefficients in the model equations, each dopant is modeled as a Gaussian distribution such that the doping concentration at point $x$ is given by [39]

[TABLE]

where $x_{j}$ and $C_{j}$ are the position and the charge of the $j$-th dopant, respectively. In the source and drain, to determine the position of random dopant, two random points (according to the two-dimensional problem) are used to translate it. For instance, the random variable $\omega=(\frac{1}{2},\frac{1}{2})$ transforms the dopants to the center of a region ($x_{j}(\omega)$). Here, we assume that both regions have same equal number of dopants. Also, $\sigma:=0.35\,\mathrm{nm}$ corresponds to the extent of the electrostatic influence, and the results are not significantly sensitive to the value of $\sigma$. Finally, the source and drain regions contain n-type dopants corresponding to a continuous doping concentration of $1\cdot 10^{19}\,\mathrm{cm^{-3}}$ (heavily doped) and the doping concentration of the channel is $1\cdot 10^{16}\,\mathrm{cm^{-3}}$.

In the following, we strive to draw a comparison between adaptive and uniform MLMC-FE methods.

Uniform and refined meshes corresponding to the device shown in Figure 1 are depicted in Figure 2. In adaptive refinement, we use the marking strategy introduced in [35]. Here, for each element $T\in\mathcal{T}_{\ell}$, the local refinement indicator $\eta_{T}$ satisfies

[TABLE]

where $\mathcal{M}$ is the set of elements marked for refinement and the associated error estimator is defined as

[TABLE]

In other words, we refine the smallest subset of elements whose corresponding error indicators in sum exceed the threshold $\theta\eta^{2}$.

The adaptive algorithm for the boundary-value problem is shown Algorithm 1. In the multilevel setting, the mesh is refined as long as $\mathcal{E}^{2}_{L}$ is greater than or equal to $\varepsilon^{2}/2$ and the number of samples are obtained according to the optimization problem (48). Also, the same number and positions of random variables are used on all levels $\ell\in\{1,\ldots,L\}$. In the numerical example, we set {\color[rgb]{0,0,0}\theta:=0.6} and the initial mesh $\mathcal{T}_{0}$ and its uniform refinement are depicted in Figure 2.

The adaptively refined meshes for $\ell\in\{1,\ldots,6\}$ for the coupled system of equations are shown in Figures 3, 4, and 5. As shown, most of the meshes have been refined due to the randomness in the source and drain areas. Similarly, the interface condition between the insulator and the channel $\Gamma$ gives rise to more refinements, where less mesh refinement occurred in the channel (green triangles). The corresponding degree of freedom for the Poisson ($\mathcal{N}_{\mathcal{P}_{\ell}}$) and drift-diffusion ($\mathcal{N}_{\mathcal{DD}_{\ell}}$) equations are summarized in Table 1 and Table 2, respectively. We compare the obtained degrees of freedom for adaptive and uniform refinement, where the initial mesh is the same in both cases.

In the next step, we compare the discretization error of uniform and adaptive refinement with 100 samples. Figure 6 indicates that the adaptive refinement reduces the discretization error, i.e., $\mathcal{E}_{\ell}$ (with respect to a reference solution), and obtains a convergence rate of $\alpha=1.46$. However, the uniform refinement leads to a smaller convergence rate of $\alpha=1.11$. Obviously, using adaptive refinement enables us to obtain a rate closer to the optimal convergence rate given in Corollary 3.1. Here, we also show the error estimator ($\eta$) for adaptive and uniform mesh-refinement to confirm its efficiency and reliability.

Moreover, we compare the statistical error of both multilevel methods. Figure 7 illustrates the decay of the variance for different degrees of freedom. The results show that similar to the discretization error, the variance in the adaptive approach is reduced faster ($\beta=2.27$) compared to uniform refinement ($\beta=1.73$). Again, the efficiency of the adaptive method is shown by the numerical results indicating error and error estimator for adaptive and uniform refinements. Also, $C_{0}=0.041$ is obtained as the variance of level $\ell=0$.

In order to estimate the optimal computational complexity, we solve the optimization problem (48). An interior point method can be used to solve the global optimization problem [10], where the results are the optimal number of samples. For different tolerances $\varepsilon$, the optimal values are summarized in Table 3 and Table 4 for uniform and adaptive refinements, respectively. In multilevel methods, most of the work is performed on coarse levels. The main reason is the reduction of variance on the finer grids.

The computational work $\sum_{\ell=0}^{L}M_{\ell}\mathcal{N}_{\ell}$ for both refinement methods are depicted in Figure 8. We here observe a significantly better efficiency of the adaptive model compared with the uniform approach. As depicted in the figure, the computational cost asymptotically behaves like $O(\varepsilon^{-2})$ for both multilevel techniques, which agrees with [40]. A more interesting computational result achieved regarding the CPU time. As Figure  8 shows since in the adaptive approach (compared to the uniform refinement) fewer degrees of freedom are needed, in different levels, lower computational time is used (to obtain same error tolerance). The difference between two computational times is more pronounced for lower prescribing errors that indicates the adaptive technique efficiency.

7 Conclusions

We have presented an adaptive MLMC-FE method for the numerical solution of the stochastic drift-diffusion-Poisson system. First, we proved an a-priori error estimate for the coupled system of equations with non-zero recombination rate. The error estimate points out how fast the error decreases as the mesh size decreases and can be considered as a useful measure of the efficiency of a given finite-element method. Also, using the stochastic numerical example, we estimated the convergence rate of the discretization error.

Secondly, a practically useful a-posteriori error indicator to bound the discretization error for the coupled system of equations was derived. From a computational point of view, the error estimator is inexpensive to estimate and guarantees the bounds on the error on all points of the geometry. The error indicator was used to design an adaptive refinement strategy to refine the mesh, where all coefficients in the system of equations can be random. In future works, the error analysis performed here can be implemented for extended systems of equations accounting for quantum effects at first order in the direction perpendicular to the gate (e.g., density gradient corrections to the charge term in the Poisson equation).

Regarding numerical examples, we implemented this adaptive MLMC-FE method to quantify noise and variations in nanoscale transistors as a real-world example. To this end, we defined a strategy to refine the meshes in the stochastic setting. The new technique was compared to the multilevel method with uniform refinement as a useful benchmark. Better convergence of the discretization error and better decay of variance were observed indicating the efficiency of the new approach. Finally, we employed an optimization problem to minimize the computational complexity. The optimal numbers of samples are obtained as the solution of the global optimization problem. The results indicate that in addition to a better control of error, a noticeable reduction of the computational work/time are achieved by the adaptive method.

8 Acknowledgments

The first and the last authors acknowledge support by FWF (Austrian Science Fund) START project no. Y660 PDE Models for Nanotechnology. The second author also acknowledges support by FWE project no. P28367-N35. The authors appreciate useful comments by Markus Melenk (TU Wien). The authors also acknowledge the helpful comments by three anonymous reviewers.