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.
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 being and (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 order , with 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 .
2. Setting of the Optimal Design Problems
Consider a nanoscale semiconductor electronic device composed of layers occupying positions . The local potential energy at the th layer is denoted by , . For , the potential energy is denoted by and for it is . It is assumed that a single electron propagating from is incident at and that a voltage bias 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:where is the vector of local layer potentials in the device, and is the characteristic function of the interval , .
The transmission coefficient of the device is defined as the ratio of current density transmitted from the device at and the incident one at . As explained in detail in [1], may be expressed aswhere and , for values of the energy . The cases and may be treated in a similar way. Here is the effective mass of the electron, denotes the electron charge, is Planck’s constant, is the electron energy, and finally solves the following boundary-value problem for the Schrödinger equation: Here denotes the unit imaginary number and the amplitude of the transmitted wave at .
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 coefficient , which is defined for , and lower, , and upper, , bounds for the local layer potentials, with ,where is given by (3) with and , .
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 variablesto the vector of local layer potentials. Here represents a random event and thus is a small unknown error in manufacturing the local potential . Hence, the cost functional considered in problem (5) becomes a random variable given byIn order to obtain a design of the potential energy profile less sensitive with respect to fabrication of unknown fluctuations, the new cost functional is considered:with being a weighting parameter. Here and var denote the expectation and variance operators, respectively. Then, the robust optimization problem is formulated aswhere 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 , 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 allows us to prove its smoothness with respect to 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 Figure 1.
The WKB method proposes a solution of the Schrödinger equation in the form with , if for all , and for all .
In the context of quantum electronic devices, solutions to (4) admit a smooth representative in their classes (of regularity class ). Hence, continuity of and its first derivative at the interface point leads towhereBy writing each matrix in (11) as the product of two matrices as and solving (13) for and ,Proceeding in the same way at the point , one obtainswhere Substituting (15) into (14), we get Denoting by the product of matrices from to , the transmission coefficient (3) takes the form , where is the first entry of . More precisely, the following explicit formula for that transmission coefficient, in the case for all , is obtained:where is given by (12) andThe case for all is completely analogous. Denoting bythe transmission coefficient is given bywhere
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, ]). For the case , this condition can be expressed asand for ,In particular, (23) and (24) constrain the values of and for which the WKB method applies. Assume that : since , for all . Hence, by introducing the functioncondition (23) is satisfied whenever . Figure 2 displays the function and its associated transmission coefficient , as given by (18), for and .
(a)
(b)
Analogously in the case , by considering the functioncondition (24) holds for . Figure 3 plots the functions and the transmission coefficient given in (21) for two values of : namely, and , with the energy of an electron eV.
(a)
(b)
3.1.3. Case of Multiple Potential Barriers
Consider now the configuration described at the beginning of Section 2.1, which is composed of potential barriers with energy given by (1). Let us denote bythe potential energy of the th barrier and by
The linear relationship between the coefficients of the electron wave function at and the ones at is expressed aswhere the form of the transfer matrices depends on the relationship between and . For the cases in which WKB method may be applied, 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 is approximated by piecewise constant potential for which explicit expressions of are well known [2].
Finally, as in the case of a single potential barrier, the transmission coefficient 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 and , which have been considered in problems (5) and (9), respectively, are smooth (of regularity class ). In addition, the set of admissible designs is a compact set of . 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 , as defined in (7), is analytic with respect to the random variable . For this reason, we propose an adaptive, isotropic, sparse grid, stochastic collocation method to approximate both the cost functional 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 , we denote a complete probability space. is the set of outcomes, is the -algebra of events, and is a probability measure. , , are the images spaces of the sample space through the real-valued random variables considered in (6), and is the product space. Assuming that the distribution measure of is absolutely continuous with respect to the Lebesgue measure, there exists a joint probability density function for . Hence, is mapped onto , where is the -algebra of Borel sets on , and is the Lebesgue measure. Finally, the expectation and variance of take the formFollowing [13, 14], the isotropic sparse grid of sampling quadrature nodes is defined as follows. Starting from an integer (called the level), the index set is considered, with . The level determines the number of collocation points in the th stochastic direction, which for the case of Smolyak rule is given by Smolyak quadrature rule applied to a generic function giveswhere , with , is a quadrature rule in which the coordinates of the nodes are the same as those for the one-dimensional quadrature formula and its associated weights are the difference between those for the and levels.
It remains to analyze the question on how to properly choose the quadrature level . For this, we have the following:(1)A positive, large enough integer and a tolerance level are fixed.(2)The first- and second-order statistical moments of , as given by (30) and the first term in the right-hand side of (31), are approximated by using (34), with level . These two approximations, which are denoted by and , respectively, play the role of enriched or reference values for the exact values of the first two statistical moments of , which, obviously, cannot be explicitly computed.(3)Finally, the level is linearly increased from to until the stopping criterionis satisfied. Here and denote, respectively, approximations, by using (34) with level , of the first two statistical moments of . If (35) is not satisfied for the tolerance and the initial , then the reference level 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 ( nm) is considered, so that the total length is nm. The desired linear transmission coefficient is Quadratic and square root transmission coefficients may be treated analogously. The design is based on equally spaced bias voltages. Hence, , . The electron mass is , with Kg (which is appropriate for an electron in the conduction band of ), its energy is eV, and its charge is C. Planck’s constant is . The lower and upper bounds for the design variable are taken as eV and 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 in Section 3, and (b) the gradient is numerically computed by using finite differences. In both cases, the algorithm is initiated with eV, , and the stopping criterion, provided by the MATLAB routine fmincon, is fixed to .
First column in Table 1 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 the corresponding one obtained via the numerical gradient.
4.2. Deterministic Design versus Design under Uncertainty
Manufacturing uncertainties are modeled by random variables uniformly distributed in : that is, for all . As an illustration, the cases and are considered. As in the preceding example, the algorithm is initiated with eV, , and the stopping criterion is fixed to .
Case 1 (, which represents of error in manufacturing each one of the potentials ). The algorithm described in Section 3.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 , is satisfied for and , which corresponds to collocation nodes in the random domain Figure 4 displays the rates of convergence of the isotropic sparse grid algorithm.
Table 2 shows the results, after convergence of the algorithm, for the expectation and the standard deviation of the cost functional given by (7). As expected, the optimal design in mean () 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.
Case 2 (, which corresponds to of noise). The same procedure as in the preceding case has been applied. In this case, (35) is satisfied, with , for , and , which corresponds to collocation nodes.
The same qualitative results as in the preceding case are observed in Table 3. 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 of noise, the reduction in the mean value of for the mean optimal value (), in comparison with the deterministic optimal design, is of the order of . For the case of of noise, this reduction is of the order of . The decreasing if the standard deviation of the computed designs is also much more significant in the current case than in the case of 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 class ) 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 Section 4.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 about over brute-force, l-based, approaches. Moreover, the robustness of the design is manifest even under weight values of 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].
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
Ociel Morales was supported by CONACyT (Consejo Nacional de Ciencia y Tecnología, Mexico) (Grant no. 726714), under program Movilidad en el Extranjero (291062). Francisco Periago was supported by Ministerio de Economía y Competitividad (Spain) (Projects DPI2016-77538-R and MTM2017-83740-P) and Fundación Séneca (Agencia de Ciencia y Tecnología de la Región de Murcia (Spain)) (19274/PI/14). José A. Vallejo was supported by a CONACyT Project CB-179115.