# Robust Optimal Design of Quantum Electronic Devices

Authors: Ociel Morales; Francisco Periago; José A. Vallejo
Journal: Mathematical Problems in Engineering
 (2018)

---

## Abstract

We consider the optimal design of a sequence of quantum barriers, in order to manufacture an electronic device at the nanoscale such that the dependence of its transmission coefficient on the bias voltage is linear. The technique presented here is easily adaptable to other response characteristics. There are two distinguishing features of our approach. First, the transmission coefficient is determined using a semiclassical approximation, so we can explicitly compute the gradient of the objective function. Second, in contrast with earlier treatments, manufacturing uncertainties are incorporated in the model through random variables; the optimal design problem is formulated in a probabilistic setting and then solved using a stochastic collocation method. As a measure of robustness, a weighted sum of the expectation and the variance of a least-squares performance metric is considered. Several simulations illustrate the proposed technique, which shows an improvement in accuracy over 69% with respect to brute-force, Monte-Carlo-based methods.


---

## Body

## 1. Introduction

Nanoelectronic devices operate with extremely low intensity currents. Under these circumstances, it is desirable to have at our disposal mechanisms to produce and control electronic currents with high precision. Electronic beams are relatively easy to produce, but their filtering to obtain nanocurrents with specified properties is much more difficult. A widely used approach consists in directing the beam on a sequence of quantum barriers with an externally adjustable bias voltage applied throughout the device, obtaining a nanocurrent by quantum tunneling. One expects to be able to control the response of the device in such a way that intensity depends, say, linearly on the applied bias. From an engineering perspective, this setting naturally leads to an optimal design problem: what must be the width and height of the layers composing the barriers, supposedly fixed in number, in order to achieve this linear response? (Of course, the problem is quite general, admitting a more complex relation between the external voltage and the current, but here we deal with the linear case just for simplicity.)There should be no need to stress the importance of the solution to this problem from a practical point of view, but it must be noticed right from the start that a closed-form, analytic solution is impossible to obtain in most cases. The use of numerical computations at some stage is unavoidable, and this leads to the question of which method to use in order to obtain a good approximation to the solution. In [1] the nonconstant potential energy profile is approximated by piecewise constant potentials. Then, the propagation matrix method [2, 3] is applied to compute the transmission coefficient and, finally, the gradient of a least-squares-type objective function (which is required by the numerical solution method) is computed using the adjoint method. It is important to point out that different discretization processes, which are used to approximate objective functions and/or its gradients, may lead to quite different results. Moreover, it has been observed in some optimization problems [4] that first approximating a cost functional and then computing the gradient of the approximated one in general differ from approximating the gradient of the exact cost functional. That is, the schemes “first discretize, then optimize” and “first optimize, then discretize” do not commute in general. Also, as it will be showed later on in this paper that optimizing for the same cost functional via its exact gradient gives different solutions compared to using an approximate one.Another issue, which cannot be obviated in a realistic mathematical model, is the presence of uncertainties. There are several sources of uncertainty in the problem under consideration, one of the most important regarding the influence on the computed optimal design being the manufacturing uncertainties. Due to the smallness of the currents involved, and the narrow width of the quantum barriers needed, methods such as MBE or CVD are used to grow thin layers (in many cases, monolayers) of some material to build the barriers, two of the preferred ones beingMoS2 and GaAs (see [5, 6], for example). These methods allow the growth of even monolayers, but the difficulties inherent in the manufacturing process at a semicommercial scale lead almost inevitably to inaccuracies that ultimately lead to a potential configuration that may be different from the numerically computed, optimal one [7]. For these reasons, the problem of computing an optimal quantum profile which, in addition, is robust against those uncertainties is an important one. If there is some statistical information about the uncertainties, then the machinery of probability theory gives a framework in which we can include uncertainties (by using random variables and/or random fields) and model objective functions (by means of expectation and variance operators, among other choices). In [8], this approach has been used for the case in which the cost functional only includes the averaging of a least-squares performance metric and by using the standard Monte-Carlo method for its numerical resolution.The present work addresses the problem of the optimal design of a quantum potential profile (modeling a nanoelectronic device) in order to obtain a transmission coefficient linearly depending on an externally applied bias voltage, in the presence of manufacturing uncertainties. The transmission coefficient is explicitly computed by using a semiclassical approximation based on the WKB method. As a consequence, an explicit formula for the gradient of the cost functional is obtained. The existence of many local minima for the problem under consideration has been already reported in [1, 8]. Although at a first glance the use of a search method for a global minimum (such as a genetic algorithm) could be reasonable, the dimension of the parameter space, depending on the number of layers, increases very quickly, and the optimization techniques for finding global minima become extremely complex due to the nonconvex character of the problem. Thus, a truly random search technique does not seem adequate, so in this paper we restrict ourselves to a gradient-based minimization algorithm.A weighted sum of expectation and variance of a random least-squares performance metric is considered as a measure of robustness. The inclusion of the second-order statistical moment in the cost functional amounts to a reduction the dispersion of the random transmission coefficient and hence an increase in the robustness of the optimal design. Since the resulting integrand in the cost functional is smooth with respect to a random parameter, a sparse grid stochastic collocation method is proposed for the numerical approximation of the involved integrals in the random domain. This method preserves the parallelizable character of Monte-Carlo sampling. However, in contrast to Monte-Carlo (which is computationally very expensive, of orderOM-1/2, with M the number of random sampling points), the stochastic collocation method shows an exponential convergence with respect to the number of sampling points. Several simulations illustrate the proposed approach, which shows itself to be an improvement in accuracy over brute-force, Monte-Carlo-based ones of about 69.4%.


## 2. Setting of the Optimal Design Problems

Consider a nanoscale semiconductor electronic device composed ofN layers occupying positions x0=0<x1<⋯<xN=L. The local potential energy at the ith layer is denoted by Ui, i=1,2,…,N. For x<x0, the potential energy is denoted by U0 and for x>xN it is UN+1. It is assumed that a single electron propagating from -∞ is incident at x0 and that a voltage bias Vbias is applied across the device. A linear approximation of the underlying Poisson’s equation [1, 9] leads to the following expression for the resulting potential energy profile:(1)Vx=Vx,U,Vbias=U0,-∞<x<x0∑j=1NUjXjx-Vbiasx-x0L,x0≤x≤xNUN+1-Vbias,xN<x<+∞,where U=U1,…,UN is the vector of local layer potentials in the device, and (2)Xjx=1,xj-1≤x<xj0,otherwise,is the characteristic function of the interval [xj-1,xj[, 1≤j≤N.The transmission coefficient of the deviceT=TVbias,U is defined as the ratio of current density transmitted from the device at x=xN and the incident one at x=x0. As explained in detail in [1], T may be expressed as(3)TVbias,U=κN+1κ0ψxN,where κ0=2meE-U0/ℏ and κN+1=2meE-UN+1+Vbias/ℏ, for values of the energy E>maxU0,UN+1-Vbias. The cases E≤U0 and E≤UN+1-Vbias may be treated in a similar way. Here m is the effective mass of the electron, e denotes the electron charge, ℏ is Planck’s constant, E is the electron energy, and finally ψx solves the following boundary-value problem for the Schrödinger equation:(4)-ℏ22md2ψxdx2+Vxψx=Eψx,x0<x<xNiκ0ψx0+dψdxx0=2iκ0A0eiκ0x0,iκN+1ψxN-dψdxxN=0. Here i denotes the unit imaginary number and A0 the amplitude of the transmitted wave at x0.


### 2.1. Deterministic Optimal Design

The (deterministic) optimal design problem considered in this paper is formulated as the following nonlinear data-fitting problem: given a desired transmission coefficientT0Vbias, which is defined for Vmin≤Vbias≤Vmax, and lower, UL, and upper, UR, bounds for the local layer potentials, with 0≤UL<UR<∞,(5)MinimizeinU=U1,…,UN:JU=∑i=1MT0Vi-TVi,U2subjecttoUL≤Uj≤UR,j=1,…,N,where TVi,U is given by (3) with Vbias=Vi and Vmin≤Vi≤Vmax, i=1,…,M.


### 2.2. Optimal Design under Manufacturing Uncertainties

As indicated in the Introduction, it is very convenient to analyze the robustness of optimal designs with respect to manufacturing uncertainties. These may be modeled by adding a vector of random variables(6)Xω=X1ω,…,XNωto the vector U of local layer potentials. Here ω represents a random event and thus Xjω is a small unknown error in manufacturing the local potential Uj. Hence, the cost functional J considered in problem (5) becomes a random variable given by(7)JU,ω=∑i=1MT0Vi-TVi,U+Xω2.In order to obtain a design of the potential energy profile U less sensitive with respect to fabrication of unknown fluctuations, the new cost functional is considered:(8)JαU=EJU,·+αvarJU,·,with α≥0 being a weighting parameter. Here E and var denote the expectation and variance operators, respectively. Then, the robust optimization problem is formulated as(9)MinimizeinU=U1,…,UN:JαUsubjecttoUL≤Uj≤UR,j=1,…,N,where JαU is given by (8).


## 2.1. Deterministic Optimal Design

The (deterministic) optimal design problem considered in this paper is formulated as the following nonlinear data-fitting problem: given a desired transmission coefficientT0Vbias, which is defined for Vmin≤Vbias≤Vmax, and lower, UL, and upper, UR, bounds for the local layer potentials, with 0≤UL<UR<∞,(5)MinimizeinU=U1,…,UN:JU=∑i=1MT0Vi-TVi,U2subjecttoUL≤Uj≤UR,j=1,…,N,where TVi,U is given by (3) with Vbias=Vi and Vmin≤Vi≤Vmax, i=1,…,M.


## 2.2. Optimal Design under Manufacturing Uncertainties

As indicated in the Introduction, it is very convenient to analyze the robustness of optimal designs with respect to manufacturing uncertainties. These may be modeled by adding a vector of random variables(6)Xω=X1ω,…,XNωto the vector U of local layer potentials. Here ω represents a random event and thus Xjω is a small unknown error in manufacturing the local potential Uj. Hence, the cost functional J considered in problem (5) becomes a random variable given by(7)JU,ω=∑i=1MT0Vi-TVi,U+Xω2.In order to obtain a design of the potential energy profile U less sensitive with respect to fabrication of unknown fluctuations, the new cost functional is considered:(8)JαU=EJU,·+αvarJU,·,with α≥0 being a weighting parameter. Here E and var denote the expectation and variance operators, respectively. Then, the robust optimization problem is formulated as(9)MinimizeinU=U1,…,UN:JαUsubjecttoUL≤Uj≤UR,j=1,…,N,where JαU is given by (8).


## 3. Solving the Optimal Design Problems

The numerical resolution of the optimal design problems stated in the preceding section requires the computation of the transmission coefficient (3) and, therefore, the resolution of the boundary-value problem (4). This problem may be numerically approximated by standard numerical methods such as finite differences or finite elements. Another approach is proposed in [1] where, after approximating the potential Vx, as given by (1), by piecewise constant potentials, problem (4) is transformed into a two-dimensional linear nonautonomous difference equation. Here we propose a different approach based on the WKB method [10]. From the point of view of optimization, WKB method is very appealing since, within its range of validity, it provides an explicit form for the solution to (4). From this, explicit expressions for the gradients of the cost functionals considered in problems (5) and (9) are derived. In addition, having an explicit expression for JU,ω allows us to prove its smoothness with respect to U and ω. From this, both existence of solutions to (5) and (9) and designing a computationally very efficient numerical resolution method will be derived in this section.We begin by explicitly computing the transmission coefficient (3) and then describe the numerical resolution methods for problems (5) and (9).


### 3.1. Explicit Computation of Transmission Coefficient

#### 3.1.1. Case of a Single Potential Barrier

For the sake of clarity, consider first the case of a single potential barrier as illustrated in Figure1.Figure 1
Sketch of the one-dimensional potential energy barrierV(x)=U-Vbias(x-x1/x2-x1), x1≤x≤x2. An electron of mass m, charge e, and energy E, incident from left, has wave vector kj in region j.The WKB method proposes a solution of the Schrödinger equation in the form(10)ψx=ψ1x=A1ei/ħκ1x+B1e-i/ħκ1x,x<x1ψ2x=A2κ2xei/ħ∫x1xκ2sds+B2κ2xe-i/ħ∫x1xκ2sds,x1≤x≤x2ψ3x=A3ei/ħκ3x+B3e-i/ħκ3x,x>x2, with κ1=2meE, κ2(x)=2me(E-V(x)) if E>Vx for all x1≤x≤x2, and κ3(x)=2me(E+Vbias) for all x>x2.In the context of quantum electronic devices, solutions to (4) admit a smooth representative in their L2 classes (of regularity class C1). Hence, continuity of ψ and its first derivative at the interface point x1 leads to(11)ei/ħκ1x1e-i/ħκ1x1iħκ1ei/ħκ1x1-iħκ1e-i/ħκ1x1A1B1=1κ2x11κ2x1-C2κ22x1κ2x1+iκ2x1κ2x1-C2κ22x1κ2x1-iκ2x1κ2x1A2B2,where(12)C=CVbias=meVbiasx2-x1.By writing each 2×2 matrix in (11) as the product of two matrices as (13)11iħκ1-iħκ1︸Kκ1ei/ħκ1x100e-i/ħκ1x1︸Eκ1,x1A1B1=11-C2κ2x1+iκ2x1ħ-C2κ2x1-iκ2x1ħ︸Kκ2,x11κ2x1001κ2x1︸Eκ2,x1A2B2and solving (13) for A1 and B1,(14)A1B1=E-1κ1,x1K-1κ1Kκ2,x1Eκ2,x1A2B2.Proceeding in the same way at the point x2, one obtains(15)A2B2=E-1κ2,x2K-1κ2,x2Kκ3Eκ3,x2A3B3,where (16)Kκj,x=11-C2κj2x+iħκjx-C2κj2x-iħκjx,Eκj,x=ei/ħ∫xj-1xκjsdsκjx00e-i/ħ∫xj-1xκjsdsκjx. Substituting (15) into (14), we get(17)A1B1=E-1κ1,x1K-1κ1Kκ2,x1Eκ2,x1E-1κ2,x2K-1κ2,x2Kκ3Eκ3,x2A3B3. Denoting by M the product of matrices from E-1(κ1,x1) to E(κ3,x2), the transmission coefficient (3) takes the form T=|1/M11|2, where M11 is the first entry of M. More precisely, the following explicit formula for that transmission coefficient, in the case E>Vx for all x1≤x≤x2, is obtained:(18)TVbias,U=κ2x1κ2x2ħCcos⁡I/ħ4κ1κ22x1-ħκ2x1Ccos⁡I/ħ4κ1κ23x2-κ2x1sin⁡I/ħ2κ1-ħ2C2sin⁡I/ħ8κ1κ22x1κ23x2-κ3sin⁡I/ħ2κ2x22+cos⁡I/ħ2+κ3κ2x1cos⁡I/ħ2κ1κ2x2-ħCsin⁡I/ħ4κ23x2+ħCκ3sin⁡I/ħ4κ1κ22x1κ2x22-1,where C is given by (12) and(19)I=IVbias,U=22mex2-x13mVbiasE-U+Vbias3/2-E-U3/2.The case E<Vx for all x1≤x≤x2 is completely analogous. Denoting by(20)κ¯2x=2meVx-E,x1≤x≤x2,the transmission coefficient is given by(21)TVbias,U=κ¯2x1κ¯2x2κ¯2x1sinh⁡I¯/ħ2κ1-ħCcosh⁡I¯/ħ4κ1κ¯22x1+ħCκ¯2x1cosh⁡I¯/ħ4κ1κ¯23x2-ħ2C2sinh⁡I¯/ħ8κ1κ¯22x1κ¯23x2-κ3sinh⁡I¯/ħ2κ¯2x22+cosh⁡I¯/ħ2+κ¯2x1κ3cosh⁡I¯/ħ2κ1κ¯2x2+ħCsinh⁡I¯/ħ4κ¯23x2-ħCκ3sinh⁡I¯/ħ4κ1κ¯2x12κ¯2x22-1,where(22)I¯=I¯Vbias,U=22mex2-x13mVbiasU-E3/2-U-E-Vbias3/2.


#### 3.1.2. Range of Validity for WKB Method

The condition for the validity of WKB method and therefore for formulas (18) and (21) is that the change in the potential energy over the decay length be smaller than the magnitude of the kinetic energy (see, for instance, [9, p.483]). For the case E>Vx, this condition can be expressed as(23)ħdV/dxκ2x<Vx-E∀x1≤x≤x2,and for E<V,(24)ħdV/dxκ¯2x<Vx-E∀x1≤x≤x2.In particular, (23) and (24) constrain the values of U and Vbias for which the WKB method applies. Assume that U>Vbias: since Vx≥U-Vbias, κ2x≥2meE-U for all x1≤x≤x2. Hence, by introducing the function(25)FE>VVbias,U=E-U-ħVbiasx2-x12meE-U,condition (23) is satisfied whenever FE>VVbias,U>0. Figure 2 displays the function FE>VVbias,U and its associated transmission coefficient TVbias,U, as given by (18), for U=0.48 and U=0.55.Figure 2
(a) Picture ofFE>VVbias,U for U=0.48 eV (dashed line) and U=0.55 eV (continuous line). (b) Transmission coefficient TVbias,U for U=0.48 eV (dashed line) and U=0.55 eV (continuous line). Range of validity for WKB method is highlighted. In both pictures, effective mass is m=0.07m0, where m0 is the bare electron mass, electron’s energy is E=0.7 eV, and the thickness barrier is 1 nm.
(a) (b)Analogously in the caseE<V, by considering the function(26)FE<VVbias,U=E-U+Vbias-ħVbiasx2-x12meE-U+Vbias,condition (24) holds for FE<VVbias,U>0. Figure 3 plots the functions FE<VVbias,U and the transmission coefficient TVbias,U given in (21) for two values of U: namely, U=0.45 and U=0.55, with the energy of an electron E=0.0.26 eV.Figure 3
(a) Picture ofFE<VVbias,U for U=0.48 eV (dashed line) and U=0.55 eV (continuous line). (b) Transmission coefficient TVbias,U for U=0.48 eV (dashed line) and U=0.55 eV (continuous line). Range of validity for WKB method is highlighted. In both pictures, effective mass is m=0.07m0, where m0 is the bare electron mass, electron’s energy is E=0.026 eV, and the thickness barrier is 1 nm.
(a) (b)


#### 3.1.3. Case of Multiple Potential Barriers

Consider now the configuration described at the beginning of Section2.1, which is composed of N potential barriers with energy given by (1). Let us denote by(27)Vjx=Uj-Vbiasx-x0L,xj-1≤x<xj,the potential energy of the jth barrier and by(28)κjx=2meE-U0,-∞<x<x0,j=02meE-Vjx,xj-1≤x≤xj,j=1,2,…,N2meE-UN+1+Vbias,xN<x<+∞,j=N+1.The linear relationship between the coefficientsA0,B0 of the electron wave function at x0 and the ones AN,BN at xN is expressed as(29)A0B0=E-1κ0,x0K-1k0·∏j=1NMj·KκN+1ExN,κN+1ANBN,where the form of the transfer matrices Mj depends on the relationship between E and Vj. For the cases in which WKB method may be applied, Mj has been computed in the preceding section. For the spatial regions for which WKB method does not apply, appropriate connection formulas must be used. For instance, as in [3, 8], the linear potential Vjx is approximated by piecewise constant potential for which explicit expressions of Mj are well known [2].Finally, as in the case of a single potential barrier, the transmission coefficientT=TVbias,U is obtained from (29).


### 3.2. Numerical Resolution Method

Before describing the numerical methods proposed in this paper for solving (5) and (9), let us briefly consider the question concerning the existence of solution for such problems.From the explicit expressions (18) and (21), and taking into consideration the computations in Section 3.1.3, it is clear that the functionals JU and JαU, which have been considered in problems (5) and (9), respectively, are smooth (of regularity class C∞). In addition, the set of admissible designs Uad=U=U1,…,UN∈RN:UL≤Uj≤UL,1≤j≤N is a compact set of RN. As a consequence, the following existence result holds.Theorem 1.
Problems (5) and (9) have, at least, one solution.The nonlinear mathematical programming problem (5) is standard and may be solved by several methods, typically by a gradient-based method. In this paper, a subspace trust region method, which is based on the interior-reflective Newton method, as proposed in [11, 12], is used. This algorithm is implemented in the MATLAB constrained optimization routine fmincon.More challenging is the robust optimal design problem (9). The brute-force sampling Monte-Carlo is the most commonly method used to solve this kind of problems. However, for smooth (with respect to the random parameter) functions, sparse grid, stochastic collocation methods are able to keep the same accuracy as Monte-Carlo and, in addition, are computationally much more efficient [13]. In our case, again from (18) and (21) it is not hard to show that JU,ω, as defined in (7), is analytic with respect to the random variable Xω. For this reason, we propose an adaptive, isotropic, sparse grid, stochastic collocation method to approximate both the cost functional JαU and its gradient. Then, the robust optimal design problem (9) can be solved as in the deterministic case.In order to explain the method for approximating integrals in the random domain used in this paper, let us first introduce some notation. ByΩ,F,P, we denote a complete probability space. Ω is the set of outcomes, F is the σ-algebra of events, and P:F→0,1 is a probability measure. Γj=XjΩ, 1≤j≤N, are the images spaces of the sample space Ω through the real-valued random variables Xj considered in (6), and Γ=∏j=1NΓj is the product space. Assuming that the distribution measure of XΩ is absolutely continuous with respect to the Lebesgue measure, there exists a joint probability density function ρ:Γ→R+ for X=X1,…,XN. Hence, Ω,F,P is mapped onto Γ,B,ρzdz, where B is the σ-algebra of Borel sets on Γ, and dz is the Lebesgue measure. Finally, the expectation and variance of JU,ω take the form(30)EJU,·=∫ΓJU,zρzdz,(31)varJU,·=∫ΓJ2U,zρzdz-∫ΓJU,zρzdz2.Following [13, 14], the isotropic sparse grid of sampling quadrature nodes is defined as follows. Starting from an integer l (called the level), the index set (32)Il,N=i=i1,…,iN∈N+N:∑n=1Nin-1≤lis considered, with N+=1,2,3,…. The level l determines the number of collocation points Rin in the nth stochastic direction, which for the case of Smolyak rule is given by(33)Rin=1,forin=12in-1+1,forin>1. Smolyak quadrature rule applied to a generic function F:Γ→R gives(34)∫ΓFzρzdz≈∑i∈Il,NΔi1⊗⋯⊗ΔiNF=∑r1=1Ri1⋯∑rN=1RiNFz1r1,…,zNrNw1r1⋯wNrN,where Δin=Qin-Qin-1, with Q0=0, is a quadrature rule in which the coordinates znrn of the nodes are the same as those for the one-dimensional quadrature formula Qin and its associated weights wnrn are the difference between those for the in and in-1 levels.It remains to analyze the question on how to properly choose the quadrature levell. For this, we have the following:(1)
A positive, large enough integerl¯ and a tolerance level 0<ε≪1 are fixed.(2)
The first- and second-order statistical moments ofJU,z, as given by (30) and the first term in the right-hand side of (31), are approximated by using (34), with level l¯. These two approximations, which are denoted by M1,l¯JU and M2,l¯JU, respectively, play the role ofenriched or reference values for the exact values of the first two statistical moments of JU,z, which, obviously, cannot be explicitly computed.(3)
Finally, the levell is linearly increased from l=1 to lopt<l¯ until the stopping criterion(35)max⁡M1,l¯JU-M1,lJUM1,l¯JU,M2,l¯JU-M2,lJUM2,l¯JU≤ε,is satisfied. Here M1,lJU and M2,lJU denote, respectively, approximations, by using (34) with level l, of the first two statistical moments of JU,z. If (35) is not satisfied for the tolerance ε and the initial l¯, then the reference level l¯ is increased.For more details on this adaptive algorithm, including its convergence, we refer the reader to [13] and the references therein.


## 3.1. Explicit Computation of Transmission Coefficient

### 3.1.1. Case of a Single Potential Barrier

For the sake of clarity, consider first the case of a single potential barrier as illustrated in Figure1.Figure 1
Sketch of the one-dimensional potential energy barrierV(x)=U-Vbias(x-x1/x2-x1), x1≤x≤x2. An electron of mass m, charge e, and energy E, incident from left, has wave vector kj in region j.The WKB method proposes a solution of the Schrödinger equation in the form(10)ψx=ψ1x=A1ei/ħκ1x+B1e-i/ħκ1x,x<x1ψ2x=A2κ2xei/ħ∫x1xκ2sds+B2κ2xe-i/ħ∫x1xκ2sds,x1≤x≤x2ψ3x=A3ei/ħκ3x+B3e-i/ħκ3x,x>x2, with κ1=2meE, κ2(x)=2me(E-V(x)) if E>Vx for all x1≤x≤x2, and κ3(x)=2me(E+Vbias) for all x>x2.In the context of quantum electronic devices, solutions to (4) admit a smooth representative in their L2 classes (of regularity class C1). Hence, continuity of ψ and its first derivative at the interface point x1 leads to(11)ei/ħκ1x1e-i/ħκ1x1iħκ1ei/ħκ1x1-iħκ1e-i/ħκ1x1A1B1=1κ2x11κ2x1-C2κ22x1κ2x1+iκ2x1κ2x1-C2κ22x1κ2x1-iκ2x1κ2x1A2B2,where(12)C=CVbias=meVbiasx2-x1.By writing each 2×2 matrix in (11) as the product of two matrices as (13)11iħκ1-iħκ1︸Kκ1ei/ħκ1x100e-i/ħκ1x1︸Eκ1,x1A1B1=11-C2κ2x1+iκ2x1ħ-C2κ2x1-iκ2x1ħ︸Kκ2,x11κ2x1001κ2x1︸Eκ2,x1A2B2and solving (13) for A1 and B1,(14)A1B1=E-1κ1,x1K-1κ1Kκ2,x1Eκ2,x1A2B2.Proceeding in the same way at the point x2, one obtains(15)A2B2=E-1κ2,x2K-1κ2,x2Kκ3Eκ3,x2A3B3,where (16)Kκj,x=11-C2κj2x+iħκjx-C2κj2x-iħκjx,Eκj,x=ei/ħ∫xj-1xκjsdsκjx00e-i/ħ∫xj-1xκjsdsκjx. Substituting (15) into (14), we get(17)A1B1=E-1κ1,x1K-1κ1Kκ2,x1Eκ2,x1E-1κ2,x2K-1κ2,x2Kκ3Eκ3,x2A3B3. Denoting by M the product of matrices from E-1(κ1,x1) to E(κ3,x2), the transmission coefficient (3) takes the form T=|1/M11|2, where M11 is the first entry of M. More precisely, the following explicit formula for that transmission coefficient, in the case E>Vx for all x1≤x≤x2, is obtained:(18)TVbias,U=κ2x1κ2x2ħCcos⁡I/ħ4κ1κ22x1-ħκ2x1Ccos⁡I/ħ4κ1κ23x2-κ2x1sin⁡I/ħ2κ1-ħ2C2sin⁡I/ħ8κ1κ22x1κ23x2-κ3sin⁡I/ħ2κ2x22+cos⁡I/ħ2+κ3κ2x1cos⁡I/ħ2κ1κ2x2-ħCsin⁡I/ħ4κ23x2+ħCκ3sin⁡I/ħ4κ1κ22x1κ2x22-1,where C is given by (12) and(19)I=IVbias,U=22mex2-x13mVbiasE-U+Vbias3/2-E-U3/2.The case E<Vx for all x1≤x≤x2 is completely analogous. Denoting by(20)κ¯2x=2meVx-E,x1≤x≤x2,the transmission coefficient is given by(21)TVbias,U=κ¯2x1κ¯2x2κ¯2x1sinh⁡I¯/ħ2κ1-ħCcosh⁡I¯/ħ4κ1κ¯22x1+ħCκ¯2x1cosh⁡I¯/ħ4κ1κ¯23x2-ħ2C2sinh⁡I¯/ħ8κ1κ¯22x1κ¯23x2-κ3sinh⁡I¯/ħ2κ¯2x22+cosh⁡I¯/ħ2+κ¯2x1κ3cosh⁡I¯/ħ2κ1κ¯2x2+ħCsinh⁡I¯/ħ4κ¯23x2-ħCκ3sinh⁡I¯/ħ4κ1κ¯2x12κ¯2x22-1,where(22)I¯=I¯Vbias,U=22mex2-x13mVbiasU-E3/2-U-E-Vbias3/2.


### 3.1.2. Range of Validity for WKB Method

The condition for the validity of WKB method and therefore for formulas (18) and (21) is that the change in the potential energy over the decay length be smaller than the magnitude of the kinetic energy (see, for instance, [9, p.483]). For the case E>Vx, this condition can be expressed as(23)ħdV/dxκ2x<Vx-E∀x1≤x≤x2,and for E<V,(24)ħdV/dxκ¯2x<Vx-E∀x1≤x≤x2.In particular, (23) and (24) constrain the values of U and Vbias for which the WKB method applies. Assume that U>Vbias: since Vx≥U-Vbias, κ2x≥2meE-U for all x1≤x≤x2. Hence, by introducing the function(25)FE>VVbias,U=E-U-ħVbiasx2-x12meE-U,condition (23) is satisfied whenever FE>VVbias,U>0. Figure 2 displays the function FE>VVbias,U and its associated transmission coefficient TVbias,U, as given by (18), for U=0.48 and U=0.55.Figure 2
(a) Picture ofFE>VVbias,U for U=0.48 eV (dashed line) and U=0.55 eV (continuous line). (b) Transmission coefficient TVbias,U for U=0.48 eV (dashed line) and U=0.55 eV (continuous line). Range of validity for WKB method is highlighted. In both pictures, effective mass is m=0.07m0, where m0 is the bare electron mass, electron’s energy is E=0.7 eV, and the thickness barrier is 1 nm.
(a) (b)Analogously in the caseE<V, by considering the function(26)FE<VVbias,U=E-U+Vbias-ħVbiasx2-x12meE-U+Vbias,condition (24) holds for FE<VVbias,U>0. Figure 3 plots the functions FE<VVbias,U and the transmission coefficient TVbias,U given in (21) for two values of U: namely, U=0.45 and U=0.55, with the energy of an electron E=0.0.26 eV.Figure 3
(a) Picture ofFE<VVbias,U for U=0.48 eV (dashed line) and U=0.55 eV (continuous line). (b) Transmission coefficient TVbias,U for U=0.48 eV (dashed line) and U=0.55 eV (continuous line). Range of validity for WKB method is highlighted. In both pictures, effective mass is m=0.07m0, where m0 is the bare electron mass, electron’s energy is E=0.026 eV, and the thickness barrier is 1 nm.
(a) (b)


### 3.1.3. Case of Multiple Potential Barriers

Consider now the configuration described at the beginning of Section2.1, which is composed of N potential barriers with energy given by (1). Let us denote by(27)Vjx=Uj-Vbiasx-x0L,xj-1≤x<xj,the potential energy of the jth barrier and by(28)κjx=2meE-U0,-∞<x<x0,j=02meE-Vjx,xj-1≤x≤xj,j=1,2,…,N2meE-UN+1+Vbias,xN<x<+∞,j=N+1.The linear relationship between the coefficientsA0,B0 of the electron wave function at x0 and the ones AN,BN at xN is expressed as(29)A0B0=E-1κ0,x0K-1k0·∏j=1NMj·KκN+1ExN,κN+1ANBN,where the form of the transfer matrices Mj depends on the relationship between E and Vj. For the cases in which WKB method may be applied, Mj has been computed in the preceding section. For the spatial regions for which WKB method does not apply, appropriate connection formulas must be used. For instance, as in [3, 8], the linear potential Vjx is approximated by piecewise constant potential for which explicit expressions of Mj are well known [2].Finally, as in the case of a single potential barrier, the transmission coefficientT=TVbias,U is obtained from (29).


## 3.1.1. Case of a Single Potential Barrier

For the sake of clarity, consider first the case of a single potential barrier as illustrated in Figure1.Figure 1
Sketch of the one-dimensional potential energy barrierV(x)=U-Vbias(x-x1/x2-x1), x1≤x≤x2. An electron of mass m, charge e, and energy E, incident from left, has wave vector kj in region j.The WKB method proposes a solution of the Schrödinger equation in the form(10)ψx=ψ1x=A1ei/ħκ1x+B1e-i/ħκ1x,x<x1ψ2x=A2κ2xei/ħ∫x1xκ2sds+B2κ2xe-i/ħ∫x1xκ2sds,x1≤x≤x2ψ3x=A3ei/ħκ3x+B3e-i/ħκ3x,x>x2, with κ1=2meE, κ2(x)=2me(E-V(x)) if E>Vx for all x1≤x≤x2, and κ3(x)=2me(E+Vbias) for all x>x2.In the context of quantum electronic devices, solutions to (4) admit a smooth representative in their L2 classes (of regularity class C1). Hence, continuity of ψ and its first derivative at the interface point x1 leads to(11)ei/ħκ1x1e-i/ħκ1x1iħκ1ei/ħκ1x1-iħκ1e-i/ħκ1x1A1B1=1κ2x11κ2x1-C2κ22x1κ2x1+iκ2x1κ2x1-C2κ22x1κ2x1-iκ2x1κ2x1A2B2,where(12)C=CVbias=meVbiasx2-x1.By writing each 2×2 matrix in (11) as the product of two matrices as (13)11iħκ1-iħκ1︸Kκ1ei/ħκ1x100e-i/ħκ1x1︸Eκ1,x1A1B1=11-C2κ2x1+iκ2x1ħ-C2κ2x1-iκ2x1ħ︸Kκ2,x11κ2x1001κ2x1︸Eκ2,x1A2B2and solving (13) for A1 and B1,(14)A1B1=E-1κ1,x1K-1κ1Kκ2,x1Eκ2,x1A2B2.Proceeding in the same way at the point x2, one obtains(15)A2B2=E-1κ2,x2K-1κ2,x2Kκ3Eκ3,x2A3B3,where (16)Kκj,x=11-C2κj2x+iħκjx-C2κj2x-iħκjx,Eκj,x=ei/ħ∫xj-1xκjsdsκjx00e-i/ħ∫xj-1xκjsdsκjx. Substituting (15) into (14), we get(17)A1B1=E-1κ1,x1K-1κ1Kκ2,x1Eκ2,x1E-1κ2,x2K-1κ2,x2Kκ3Eκ3,x2A3B3. Denoting by M the product of matrices from E-1(κ1,x1) to E(κ3,x2), the transmission coefficient (3) takes the form T=|1/M11|2, where M11 is the first entry of M. More precisely, the following explicit formula for that transmission coefficient, in the case E>Vx for all x1≤x≤x2, is obtained:(18)TVbias,U=κ2x1κ2x2ħCcos⁡I/ħ4κ1κ22x1-ħκ2x1Ccos⁡I/ħ4κ1κ23x2-κ2x1sin⁡I/ħ2κ1-ħ2C2sin⁡I/ħ8κ1κ22x1κ23x2-κ3sin⁡I/ħ2κ2x22+cos⁡I/ħ2+κ3κ2x1cos⁡I/ħ2κ1κ2x2-ħCsin⁡I/ħ4κ23x2+ħCκ3sin⁡I/ħ4κ1κ22x1κ2x22-1,where C is given by (12) and(19)I=IVbias,U=22mex2-x13mVbiasE-U+Vbias3/2-E-U3/2.The case E<Vx for all x1≤x≤x2 is completely analogous. Denoting by(20)κ¯2x=2meVx-E,x1≤x≤x2,the transmission coefficient is given by(21)TVbias,U=κ¯2x1κ¯2x2κ¯2x1sinh⁡I¯/ħ2κ1-ħCcosh⁡I¯/ħ4κ1κ¯22x1+ħCκ¯2x1cosh⁡I¯/ħ4κ1κ¯23x2-ħ2C2sinh⁡I¯/ħ8κ1κ¯22x1κ¯23x2-κ3sinh⁡I¯/ħ2κ¯2x22+cosh⁡I¯/ħ2+κ¯2x1κ3cosh⁡I¯/ħ2κ1κ¯2x2+ħCsinh⁡I¯/ħ4κ¯23x2-ħCκ3sinh⁡I¯/ħ4κ1κ¯2x12κ¯2x22-1,where(22)I¯=I¯Vbias,U=22mex2-x13mVbiasU-E3/2-U-E-Vbias3/2.


## 3.1.2. Range of Validity for WKB Method

The condition for the validity of WKB method and therefore for formulas (18) and (21) is that the change in the potential energy over the decay length be smaller than the magnitude of the kinetic energy (see, for instance, [9, p.483]). For the case E>Vx, this condition can be expressed as(23)ħdV/dxκ2x<Vx-E∀x1≤x≤x2,and for E<V,(24)ħdV/dxκ¯2x<Vx-E∀x1≤x≤x2.In particular, (23) and (24) constrain the values of U and Vbias for which the WKB method applies. Assume that U>Vbias: since Vx≥U-Vbias, κ2x≥2meE-U for all x1≤x≤x2. Hence, by introducing the function(25)FE>VVbias,U=E-U-ħVbiasx2-x12meE-U,condition (23) is satisfied whenever FE>VVbias,U>0. Figure 2 displays the function FE>VVbias,U and its associated transmission coefficient TVbias,U, as given by (18), for U=0.48 and U=0.55.Figure 2
(a) Picture ofFE>VVbias,U for U=0.48 eV (dashed line) and U=0.55 eV (continuous line). (b) Transmission coefficient TVbias,U for U=0.48 eV (dashed line) and U=0.55 eV (continuous line). Range of validity for WKB method is highlighted. In both pictures, effective mass is m=0.07m0, where m0 is the bare electron mass, electron’s energy is E=0.7 eV, and the thickness barrier is 1 nm.
(a) (b)Analogously in the caseE<V, by considering the function(26)FE<VVbias,U=E-U+Vbias-ħVbiasx2-x12meE-U+Vbias,condition (24) holds for FE<VVbias,U>0. Figure 3 plots the functions FE<VVbias,U and the transmission coefficient TVbias,U given in (21) for two values of U: namely, U=0.45 and U=0.55, with the energy of an electron E=0.0.26 eV.Figure 3
(a) Picture ofFE<VVbias,U for U=0.48 eV (dashed line) and U=0.55 eV (continuous line). (b) Transmission coefficient TVbias,U for U=0.48 eV (dashed line) and U=0.55 eV (continuous line). Range of validity for WKB method is highlighted. In both pictures, effective mass is m=0.07m0, where m0 is the bare electron mass, electron’s energy is E=0.026 eV, and the thickness barrier is 1 nm.
(a) (b)


## 3.1.3. Case of Multiple Potential Barriers

Consider now the configuration described at the beginning of Section2.1, which is composed of N potential barriers with energy given by (1). Let us denote by(27)Vjx=Uj-Vbiasx-x0L,xj-1≤x<xj,the potential energy of the jth barrier and by(28)κjx=2meE-U0,-∞<x<x0,j=02meE-Vjx,xj-1≤x≤xj,j=1,2,…,N2meE-UN+1+Vbias,xN<x<+∞,j=N+1.The linear relationship between the coefficientsA0,B0 of the electron wave function at x0 and the ones AN,BN at xN is expressed as(29)A0B0=E-1κ0,x0K-1k0·∏j=1NMj·KκN+1ExN,κN+1ANBN,where the form of the transfer matrices Mj depends on the relationship between E and Vj. For the cases in which WKB method may be applied, Mj has been computed in the preceding section. For the spatial regions for which WKB method does not apply, appropriate connection formulas must be used. For instance, as in [3, 8], the linear potential Vjx is approximated by piecewise constant potential for which explicit expressions of Mj are well known [2].Finally, as in the case of a single potential barrier, the transmission coefficientT=TVbias,U is obtained from (29).


## 3.2. Numerical Resolution Method

Before describing the numerical methods proposed in this paper for solving (5) and (9), let us briefly consider the question concerning the existence of solution for such problems.From the explicit expressions (18) and (21), and taking into consideration the computations in Section 3.1.3, it is clear that the functionals JU and JαU, which have been considered in problems (5) and (9), respectively, are smooth (of regularity class C∞). In addition, the set of admissible designs Uad=U=U1,…,UN∈RN:UL≤Uj≤UL,1≤j≤N is a compact set of RN. As a consequence, the following existence result holds.Theorem 1.
Problems (5) and (9) have, at least, one solution.The nonlinear mathematical programming problem (5) is standard and may be solved by several methods, typically by a gradient-based method. In this paper, a subspace trust region method, which is based on the interior-reflective Newton method, as proposed in [11, 12], is used. This algorithm is implemented in the MATLAB constrained optimization routine fmincon.More challenging is the robust optimal design problem (9). The brute-force sampling Monte-Carlo is the most commonly method used to solve this kind of problems. However, for smooth (with respect to the random parameter) functions, sparse grid, stochastic collocation methods are able to keep the same accuracy as Monte-Carlo and, in addition, are computationally much more efficient [13]. In our case, again from (18) and (21) it is not hard to show that JU,ω, as defined in (7), is analytic with respect to the random variable Xω. For this reason, we propose an adaptive, isotropic, sparse grid, stochastic collocation method to approximate both the cost functional JαU and its gradient. Then, the robust optimal design problem (9) can be solved as in the deterministic case.In order to explain the method for approximating integrals in the random domain used in this paper, let us first introduce some notation. ByΩ,F,P, we denote a complete probability space. Ω is the set of outcomes, F is the σ-algebra of events, and P:F→0,1 is a probability measure. Γj=XjΩ, 1≤j≤N, are the images spaces of the sample space Ω through the real-valued random variables Xj considered in (6), and Γ=∏j=1NΓj is the product space. Assuming that the distribution measure of XΩ is absolutely continuous with respect to the Lebesgue measure, there exists a joint probability density function ρ:Γ→R+ for X=X1,…,XN. Hence, Ω,F,P is mapped onto Γ,B,ρzdz, where B is the σ-algebra of Borel sets on Γ, and dz is the Lebesgue measure. Finally, the expectation and variance of JU,ω take the form(30)EJU,·=∫ΓJU,zρzdz,(31)varJU,·=∫ΓJ2U,zρzdz-∫ΓJU,zρzdz2.Following [13, 14], the isotropic sparse grid of sampling quadrature nodes is defined as follows. Starting from an integer l (called the level), the index set (32)Il,N=i=i1,…,iN∈N+N:∑n=1Nin-1≤lis considered, with N+=1,2,3,…. The level l determines the number of collocation points Rin in the nth stochastic direction, which for the case of Smolyak rule is given by(33)Rin=1,forin=12in-1+1,forin>1. Smolyak quadrature rule applied to a generic function F:Γ→R gives(34)∫ΓFzρzdz≈∑i∈Il,NΔi1⊗⋯⊗ΔiNF=∑r1=1Ri1⋯∑rN=1RiNFz1r1,…,zNrNw1r1⋯wNrN,where Δin=Qin-Qin-1, with Q0=0, is a quadrature rule in which the coordinates znrn of the nodes are the same as those for the one-dimensional quadrature formula Qin and its associated weights wnrn are the difference between those for the in and in-1 levels.It remains to analyze the question on how to properly choose the quadrature levell. For this, we have the following:(1)
A positive, large enough integerl¯ and a tolerance level 0<ε≪1 are fixed.(2)
The first- and second-order statistical moments ofJU,z, as given by (30) and the first term in the right-hand side of (31), are approximated by using (34), with level l¯. These two approximations, which are denoted by M1,l¯JU and M2,l¯JU, respectively, play the role ofenriched or reference values for the exact values of the first two statistical moments of JU,z, which, obviously, cannot be explicitly computed.(3)
Finally, the levell is linearly increased from l=1 to lopt<l¯ until the stopping criterion(35)max⁡M1,l¯JU-M1,lJUM1,l¯JU,M2,l¯JU-M2,lJUM2,l¯JU≤ε,is satisfied. Here M1,lJU and M2,lJU denote, respectively, approximations, by using (34) with level l, of the first two statistical moments of JU,z. If (35) is not satisfied for the tolerance ε and the initial l¯, then the reference level l¯ is increased.For more details on this adaptive algorithm, including its convergence, we refer the reader to [13] and the references therein.


## 4. Numerical Simulations

In this section, numerical results for problems (5) and (9) are presented and discussed. In all experiments, a device consisting of four layers of the same thickness (1 nm) is considered, so that the total length is L=4 nm. The desired linear transmission coefficient is(36)T0Vbias=0.00002Vbias+0.0000099,0≤Vbias≤0.25. Quadratic and square root transmission coefficients may be treated analogously. The design is based on 10 equally spaced bias voltages. Hence, Vi=i0.25/10, 1≤i≤10. The electron mass is m=0.07×m0, with m0=9.10939×10-31 Kg (which is appropriate for an electron in the conduction band of AlξGa1-ξAs), its energy is E=0.026 eV, and its charge is e=1.602×10-19 C. Planck’s constant is ℏ=1.05457×10-34J·s. The lower and upper bounds for the design variable are taken as UL=0.7 eV and UH=1.7 eV, respectively.The goal of this section is twofold. On the one hand, it is aimed at analyzing the differences that may occur when using exact gradients or numerical gradients in the optimization algorithm. On the other hand, we want to analyze the influence of manufacturing uncertainties on the computed designs. We deal with these issues in the following subsections.


### 4.1. Exact versus Numerical Gradient

In this experiment, the deterministic problem (5) is solved, by using the MATLAB routine fmincon, in the cases where (a) the exact gradient is provided, as computed from the explicit expression for the cost functional JU in Section 3, and (b) the gradient is numerically computed by using finite differences. In both cases, the algorithm is initiated with Uj0=0.7 eV, 1≤j≤4, and the stopping criterion, provided by the MATLAB routine fmincon, is fixed to 10-15.First column in Table1 displays results for the values of the cost functional after convergence of the algorithm. The optimal energy potentials are showed in the remaining columns. Inspection of Table 1 reveals (as expected) that performance increases by using the exact gradient. Precisely, the value of the objective function obtained by using the optimal design as computed with the exact gradient improves in about 69.4% the corresponding one obtained via the numerical gradient.Table 1
Results for (5) using exact (first row) and numerical gradient (second row).
Gradient J U U 1 U 2 U 3 U 4 Exact 1.43 × 1 0 - 12 0.70 1.31 1.54 0.70 Numerical 4.68 × 1 0 - 12 1.03 0.96 0.96 1.02


### 4.2. Deterministic Design versus Design under Uncertainty

Manufacturing uncertainties are modeled by random variables uniformly distributed in-a,a: that is, Xj=U-a,a for all j=1,2,3,4. As an illustration, the cases a=0.05 and a=0.2 are considered. As in the preceding example, the algorithm is initiated with Uj0=0.7 eV, 1≤j≤4, and the stopping criterion is fixed to 10-15.Case 1 (a=0.05, which represents 5% of error in manufacturing each one of the potentials Uj).
The algorithm described in Section3.2 has been implemented by using the Sparse Grids MATLAB kit 15.8 (see https://csqi.epfl.ch and [13]). The stopping criterion (35), with ε=10-7, is satisfied for l¯=20 and lopt=15, which corresponds to 895 collocation nodes in the random domain Γ=-0.05,0.054. Figure 4 displays the rates of convergence of the isotropic sparse grid algorithm.Figure 4
Casea=0.05. Rates of convergence of the isotropic sparse grid algorithm. In the vertical axis, log10error represents the two terms for relative error considered in (35). Precisely, the continuous line corresponds to relative error for the first-order moment of JU,⋅, and the dashed-dotted line is its second moment. The horizontal dashed line represents the prescribed accuracy level in the stopping criterion (35).Table2 shows the results, after convergence of the algorithm, for the expectation and the standard deviation of the cost functional JU,ω given by (7). As expected, the optimal design in mean (α=0) provides a solution which reduces the impact of manufacturing errors in comparison with the deterministic approach (see first and second rows in the first column of Table 2). It is also observed that increasing the value of the weighing parameter α reduces the dispersion of the computed designs. These results are more significant when the level of uncertainty increases, as pointed out in Table 3.Table 2
Casea=0.05. Mean(first column) and standard deviation(second column) of the functional (7) for the optimal deterministic design(first row), the optimal design in mean(second row), and the optimal design for α=1012(third row).
Design E J U , ⋅ s t d J U , ⋅ U 1 U 2 U 3 U 4 Deterministic 1.737 × 1 0 - 11 2.143 × 1 0 - 11 0.70 1.31 1.54 0.70 α=0 1.664 × 1 0 - 11 1.920 × 1 0 - 11 0.95 0.90 1.60 0.70 α=1012 1.908 × 1 0 - 11 1.845 × 1 0 - 11 1.14 0.83 1.05 0.97Table 3
Casea=0.2. Mean(first column) and standard deviation(second column) of the cost functional (7) for the optimal deterministic design(first row), the optimal design in mean(second row), and the optimal design for α=1010  (third row). Columns 3 to 6 show the corresponding optimal potentials.
Design E J U , ⋅ s t d J U , ⋅ U 1 U 2 U 3 U 4 Deterministic 3.48 × 1 0 - 10 7.17 × 1 0 - 10 0.70 1.31 1.54 0.70 α=0 2.12 × 1 0 - 10 2.66 × 1 0 - 10 1.19 0.82 1.21 0.92 α=1010 2.35 × 1 0 - 10 2.15 × 1 0 - 10 1.20 0.88 1.07 1.05Case 2 (a=0.2, which corresponds to 20% of noise).
The same procedure as in the preceding case has been applied. In this case, (35) is satisfied, with ε=10-2, for l¯=20, and l=16, which corresponds to 1212 collocation nodes.The same qualitative results as in the preceding case are observed in Table3. However, since in this case the level of manufacturing uncertainties is higher, the differences of corresponding solutions are much more significant than in the preceding case. Indeed, for 20% of noise, the reduction in the mean value of JU,· for the mean optimal value (α=0), in comparison with the deterministic optimal design, is of the order of 39%. For the case of 5% of noise, this reduction is of the order of 4.2%. The decreasing if the standard deviation of the computed designs is also much more significant in the current case than in the case of 5% of manufacturing noise.


## 4.1. Exact versus Numerical Gradient

In this experiment, the deterministic problem (5) is solved, by using the MATLAB routine fmincon, in the cases where (a) the exact gradient is provided, as computed from the explicit expression for the cost functional JU in Section 3, and (b) the gradient is numerically computed by using finite differences. In both cases, the algorithm is initiated with Uj0=0.7 eV, 1≤j≤4, and the stopping criterion, provided by the MATLAB routine fmincon, is fixed to 10-15.First column in Table1 displays results for the values of the cost functional after convergence of the algorithm. The optimal energy potentials are showed in the remaining columns. Inspection of Table 1 reveals (as expected) that performance increases by using the exact gradient. Precisely, the value of the objective function obtained by using the optimal design as computed with the exact gradient improves in about 69.4% the corresponding one obtained via the numerical gradient.Table 1
Results for (5) using exact (first row) and numerical gradient (second row).
Gradient J U U 1 U 2 U 3 U 4 Exact 1.43 × 1 0 - 12 0.70 1.31 1.54 0.70 Numerical 4.68 × 1 0 - 12 1.03 0.96 0.96 1.02


## 4.2. Deterministic Design versus Design under Uncertainty

Manufacturing uncertainties are modeled by random variables uniformly distributed in-a,a: that is, Xj=U-a,a for all j=1,2,3,4. As an illustration, the cases a=0.05 and a=0.2 are considered. As in the preceding example, the algorithm is initiated with Uj0=0.7 eV, 1≤j≤4, and the stopping criterion is fixed to 10-15.Case 1 (a=0.05, which represents 5% of error in manufacturing each one of the potentials Uj).
The algorithm described in Section3.2 has been implemented by using the Sparse Grids MATLAB kit 15.8 (see https://csqi.epfl.ch and [13]). The stopping criterion (35), with ε=10-7, is satisfied for l¯=20 and lopt=15, which corresponds to 895 collocation nodes in the random domain Γ=-0.05,0.054. Figure 4 displays the rates of convergence of the isotropic sparse grid algorithm.Figure 4
Casea=0.05. Rates of convergence of the isotropic sparse grid algorithm. In the vertical axis, log10error represents the two terms for relative error considered in (35). Precisely, the continuous line corresponds to relative error for the first-order moment of JU,⋅, and the dashed-dotted line is its second moment. The horizontal dashed line represents the prescribed accuracy level in the stopping criterion (35).Table2 shows the results, after convergence of the algorithm, for the expectation and the standard deviation of the cost functional JU,ω given by (7). As expected, the optimal design in mean (α=0) provides a solution which reduces the impact of manufacturing errors in comparison with the deterministic approach (see first and second rows in the first column of Table 2). It is also observed that increasing the value of the weighing parameter α reduces the dispersion of the computed designs. These results are more significant when the level of uncertainty increases, as pointed out in Table 3.Table 2
Casea=0.05. Mean(first column) and standard deviation(second column) of the functional (7) for the optimal deterministic design(first row), the optimal design in mean(second row), and the optimal design for α=1012(third row).
Design E J U , ⋅ s t d J U , ⋅ U 1 U 2 U 3 U 4 Deterministic 1.737 × 1 0 - 11 2.143 × 1 0 - 11 0.70 1.31 1.54 0.70 α=0 1.664 × 1 0 - 11 1.920 × 1 0 - 11 0.95 0.90 1.60 0.70 α=1012 1.908 × 1 0 - 11 1.845 × 1 0 - 11 1.14 0.83 1.05 0.97Table 3
Casea=0.2. Mean(first column) and standard deviation(second column) of the cost functional (7) for the optimal deterministic design(first row), the optimal design in mean(second row), and the optimal design for α=1010  (third row). Columns 3 to 6 show the corresponding optimal potentials.
Design E J U , ⋅ s t d J U , ⋅ U 1 U 2 U 3 U 4 Deterministic 3.48 × 1 0 - 10 7.17 × 1 0 - 10 0.70 1.31 1.54 0.70 α=0 2.12 × 1 0 - 10 2.66 × 1 0 - 10 1.19 0.82 1.21 0.92 α=1010 2.35 × 1 0 - 10 2.15 × 1 0 - 10 1.20 0.88 1.07 1.05Case 2 (a=0.2, which corresponds to 20% of noise).
The same procedure as in the preceding case has been applied. In this case, (35) is satisfied, with ε=10-2, for l¯=20, and l=16, which corresponds to 1212 collocation nodes.The same qualitative results as in the preceding case are observed in Table3. However, since in this case the level of manufacturing uncertainties is higher, the differences of corresponding solutions are much more significant than in the preceding case. Indeed, for 20% of noise, the reduction in the mean value of JU,· for the mean optimal value (α=0), in comparison with the deterministic optimal design, is of the order of 39%. For the case of 5% of noise, this reduction is of the order of 4.2%. The decreasing if the standard deviation of the computed designs is also much more significant in the current case than in the case of 5% of manufacturing noise.


## 5. Conclusions

We have addressed the problem of determining optimal designs of nanoelectronic devices whose physical behavior is governed by the time-independent Schrödinger equation. Two situations are considered: a deterministic version of the problem, and the more realistic case where manufacturing uncertainties are accounted for. In both cases, the corresponding optimal design problems are formulated as the minimization of a least-squares performance metric. In the stochastic case, the variance of the random metric is incorporated in the cost functional as a measure of robustness. An explicit expression for that metric is computed using a semiclassical approximation. Having at our disposal an analytic expression for the cost functional has the following advantages:(a)
At the theoretical level, it is easily deduced that the cost functional is smooth (of regularity classC∞) with respect to the design variables. As a consequence, the existence of solutions for both optimization problems is proved.(b)
In the deterministic case, an explicit expression for the gradient of the cost functional is obtained. Numerical simulations in Section4.1 show the relevance of this issue, at least in the specific problem considered in this work, when gradient-based minimization algorithms are used as the numerical resolution method.(c)
In the stochastic optimization problem, it is also deduced that integrands, which appear in the considered cost functional, are analytic with respect to the random vector parameter. Thus, stochastic collocation methods (which possess an exponential rate convergence with respect to the number of sampling points and are, computationally, much more efficient than the classical brute-force Monte-Carlo method) are preferred.It is important to emphasize that, by using the WKB method, explicit expressions for the descent directions can be computed even for those cases for which the dependence of the current on the external voltage is more complex, containing quadratic or squared root expressions, for instance. In the same vein, the explicit expression for the integrals of the momentum, appearing in the phase of the WKB solution, can be computed by a CAS (Computer Algebra System) quite efficiently for a general class of potentials.The results obtained show an improvement of the accuracy in the linear response characteristic of about69.4% over brute-force, l-based, approaches. Moreover, the robustness of the design is manifest even under weight values of α=1012 in the variance (physically, as seen in in Table 2, this would amount to pass from a design using Germanium, with a band gap of ≃0.7 eV, to one using Silicon, whose band gap is ≃1.1 eV, so that extreme case is still physically feasible).As noticed along the paper, the semiclassical approach used in this work has the drawback that constrains the values of the applied bias voltage and those of the local energy potentials. Accordingly, the whole range of possible energies and potential profiles may be covered by using the approach proposed in this work in combination with appropriate connection formulas [1].



---

*Source: 3095257-2018-04-05.xml*