Orbital-free approach for large-scale electrostatic simulations of quantum nanoelectronics devices

The route to reliable quantum nanoelectronic devices hinges on precise control of the electrostatic environment. For this reason, accurate methods for electrostatic simulations are essential in the design process. The most widespread methods for this purpose are the Thomas-Fermi approximation, which provides quick approximate results, and the Schr\"odinger-Poisson method, which better takes into account quantum mechanical effects. The mentioned methods suffer from relevant shortcomings: the Thomas-Fermi method fails to take into account quantum confinement effects that are crucial in heterostructures, while the Schr\"odinger-Poisson method suffers severe scalability problems. This paper outlines the application of an orbital-free approach inspired by density functional theory. By introducing gradient terms in the kinetic energy functional, our proposed method incorporates corrections to the electronic density due to quantum confinement while it preserves the scalability of a theory that can be expressed as a functional minimization problem. This method offers a new approach to addressing large-scale electrostatic simulations of quantum nanoelectronic devices.


I. INTRODUCTION
The continuous increase of complexity in quantum nanoelectronic devices has created the demand for a precise description of the electrostatic behavior of these systems, as this influences all other phenomena in the devices [1][2][3].Solving the electrostatic problem requires the solution of Poisson's equation (PE) for the electric potential ϕ coupled to a functional n[ϕ] that describes the electron density as a function of the electrostatic field itself where ρ fx is the eventual fixed charge in the device and e is the electron charge.
In the context of hybrid superconductor-semiconductor quantum devices, two methods are currently widely used: the Schrödinger-Poisson (SP) and the Thomas-Fermi (TF) methods [2][3][4][5][6][7][8].The SP method takes the quantum mechanical system into account within the Hartree approximation and provides precise results, but is computationally demanding [3,8,9].The TF method, instead, uses the free-fermions gas model to approximate the charge distribution.In this way, it reasonably approximates the device's electrostatic configuration at a much lower computational cost.While the SP method provides more precise results, it is computationally much more demanding since it requires the diagonalization of the Hamiltonian.Memory and time requirements can be so high that it is practically impossible to use on realistic systems of relevant size [2,10].However, the TF method fails to properly consider quantum confinement effects as the energy functional is completely local.For this reason, a common practice is to first focus on the electrostatic problem and solve it with the TF method, and then use the calculated electrostatic potential to solve the quantum mechanical problem [1,2].We will refer to this approach, which is not strictly self-consistent, as Schrödinger-Thomas-Fermi (STF) method.
The TF method is exact for uniform systems and holds as an approximation as long as the electron density slowly varies in space.The validity of this assumption for a specific system can be checked through the Thomas-Fermi error: where k F (r) = [3π 2 n(r)] 1/3 is the Fermi wavelength [11].
The TF approximation is valid in the limit R TF 1 [12,13].This condition is never satisfied at interfaces with vacuum or insulators where the density goes abruptly to zero.Two instances of typical devices are shown in Fig. 1, showing the density calculated with the TF method and the value of the TF error.
The TF method is the simplest model in the class of orbital-free (OF) density functional theories, and some of its shortcomings are overcome in more advanced models.The defining feature of these models is that their mathematical form is exactly an energy functional of the density.This means that explicit diagonalization of the Köhn-Sham Hamiltonian is not required to determine the ground state density of the system.For this reason, OF density functional theory methods are computationally efficient and can be applied to systems with a large number of electrons.However, the accuracy of OF density functional theory is generally lower, and it is typically used for approximate calculations and to study qualitative trends.
Traditional application domains of OF theories include atomic physics and material science [14][15][16][17][18]. Similar theories are also known as semiclassical methods in nuclear matter systems that include finite nuclei and astrophysical systems [19][20][21][22].This paper introduces these methods in a different field: quantum nanoelectronic devices design.We adapt The first row shows sketches of a cross-section of a common 2DEG device (a) and a simple circular quantum dot (d).In both cases, the devices are built on top of a semiconductor stack (schematic given in Fig. 2) that provides the vertical confinement of electrons.The lateral confinement is controlled by Au gates (yellow), HfO 2 oxide (gray), and an Al wire (dark gray).b) and e) show the electron density from a TF simulation.e) is plotted in the middle of the InAs well.c) and f) show R TF as defined in Eq. ( 2).
ideas from OF density functional theory to go beyond the TF method, improve on some of its drawbacks, and minimize any additional computational cost. .The goal is thus to develop and assess a numerically cheap model in a variational form that includes corrections to the electronic density due to quantum confinement.Moreover, the application of these ideas to nanoelectronic devices can be far easier than in the usual application domains, like material science, since a precise estimation of excitation energies is not required.We expect such a method to better predict the charge density compared to the TF method.Estimating the electron number is especially crucial for floating metallic parts of devices where charging energy is a relevant parameter.Following the literature, we will refer to this orbital-free method as the extended Thomas Fermi (ETF) method [20].
A comparable theory built following similar consideration is the effective-potential potential method [23,24] later expanded in the density-gradient theory [25,26].Unlike our approach, the density gradient theory tries to introduce quantum-confinement corrections to the driftdiffusion equations commonly used in modeling traditional semiconductor devices operated at room temperature.Our proposal is much simpler and more tailored to quantum nanoelectronic devices operated at low-temperature with small currents in a condition very close to equilib-rium.
As a specific application domain of our method, we consider the task of optimizing gate shapes and providing an intuitive understanding of how the potential landscape is affected by variation in voltage biasing.These tasks do not require precise calculations but rather quick and responsive evaluation of a configuration.In this respect, our tool provides, in addition to a precise electric field, a model for the electron fluid which is closer to the Schrödinger-Poisson evaluations.
An inevitable feature of an orbital-free approach is the incapacity of modeling the quantization of charge.While in many applications this is a disadvantage, it can be a useful feature in certain situations.The two outputs of the method, ϕ and n, will depend continuously on the geometrical and material parameters of the device.This means that any derived quantity L[n, ϕ] will be continuous as well with respect to the parameters of the model.This allows for the possibility of developing new devices by using parametric modeling, where the function L can be an objective function describing some desired property of the device (e.g., the transmission probability of a quantum point contact), allowing for computation of the optimal geometrical and material parameters by numerical optimization.
For a more precise treatment of quantization, a hybrid Schrödinger-ETF method, similar to the STF method, can be used.Additionally, the ETF method not only provides a precise estimate of the electrostatic potential ϕ(r) but also of the numerical density n(r), which can be used to generate an effective Hamiltonian through dimensionality reduction..For example, one can be interested in projecting a 3D Hamiltonian on a 2D subspace to model a 2DEG device.This could be done by taking the weighted average of the parameters of the Hamiltonian, like, for example, the effective mass, Rashba field, electrostatic potential, or conduction band minimum.

II. ORBITAL FREE METHOD
In the following, we will consider only direct bandgap semiconductors.We will neglect the hole accumulation process and model only the electron fluid.However, it is straightforward to extend the method by including a hole fluid.For low densities, we can assume that the electrons in the semiconductor are close to the Γ point and can be reliably described by a parabolic band.In this case, a general orbital-free theory for the conduction band electrons can be expressed by an energy functional similar to the one for free electrons: where ϕ is the electrostatic potential while n is the numerical density of electrons [10,[13][14][15]27].The energy functional is split into several pieces: K[n] is the kinetic energy functional of the electrons, V [n] is the potential energy of the electron fluid, U [ϕ] is the electrostatic field energy, and finally, the last term is the coupling between the electrostatic field and the charge which comprises the free electrons charge distribution −en and the fixed charge of the system ρ fx .This can be due to, e.g., doping or fixed surface charge, and can be fundamental in the electrostatic modeling of semiconductor systems [28,29].The computational domain Ω can be split into three different kinds of regions: metals, insulators, and semiconductors such that Ω = Ω M ∪ Ω Sm ∪ Ω I .
The electrostatic field energy takes the standard form where ε(r) is the permittivity, while the potential is where E CBM (r) is the conduction band minimum (CBM) of the semiconductor that acts as local chemical potential.
The exchange energy can be included within the local density approximation through the V ex (n(r)) term, which we neglect in this work as including this would make the comparison with the SP method more complicated.The form of the kinetic energy functional is the most complicated part of the orbital-free theory.Finding precise kinetic energy functionals is the focus of much modern research [14][15][16][17][18].We opted for the simplest model that goes beyond the homogeneous electron gas case of the Thomas-Fermi method 5 is the Thomas-Fermi constant, λ vW is the von Weizsäcker parameter and m * is the effective mass.The kinetic energy functional incorporates a TF term and the so-called von Weizsäcker correction [10,18,20,27,30,31].This term captures the energy cost of rapid variation of the density in space and was originally introduced by von Weizsäcker to correct the TF method issues when applied to rapidly varying electron densities.
In the literature, there has been much discussion on the value of the vW coefficient, λ vW , which works as a weight of the gradient-dependent vW term [27].In the limit λ vW → 0, the TF method is recovered.In Ref. [32], the response function of a uniform system of independent fermions is investigated, and it is shown that λ vW = 1 is exact in the limit of short-wavelength perturbations, whereas λ vW = 1/9 is exact in the limit of long-wavelength perturbations.Other analysis pointed to the value of λ vW = 1/5 as the most adequate [33].
In this application, we decided to treat λ vW as a parameter of the model and empirically select a value in the [0, 1] interval that agrees with SP simulations in simple geometries.The complete form of the energy functional and some of its extensions are further discussed in Ref. [20].
It is worth noting that a functional theory that takes this form is not, strictly speaking, a density functional theory, as the electric potential is explicitly included and cannot be easily transformed into a pure density functional theory.This is caused by the fact that metal parts (like gates or floating metallic islands) are present in nanoelectronic devices.These modify the boundary condition of the electric potential equation such that the PE Green's function becomes a complicated object generally not expressible in an analytic form.Including the electric potential explicitly avoids this problem.A more detailed discussion is included in Appendix A.
Moreover, when treating a closed system like an atom or a molecule, the minimization problem is characterized by the constraint where N is the total number of electrons, that is fixed.The existence of a solution is not guaranteed when including an exchange and correlation potential [34].We do not have such a constraint in our case.The total number of electrons N is not fixed, but actually one of the results we want to determine.Moreover, it can also take fractional values.
To move from an optimization problem to a boundary value problem, we will use functional minimization.It is convenient to define the matter field ψ = √ n before proceeding.Note that ψ is a real field defined as the square root of the density.It cannot be interpreted as a wavefunction as it does not carry information about the phases of the electrons.Therefore this theory cannot describe long-range interference effects, but it just includes a local correction to the electronic density due to quantum confinement effects.The functional derivative δE δϕ = 0 returns the Poisson equation while δE δψ = 0 returns a partial differential equation that has the form of a nonlinear Schrödinger equation (NLSE).The system of coupled PDEs is then The insulators are usually large gap semiconductors that are depleted for low enough bias voltages.There-fore, we can exclude insulator regions from the NLSE domain and assign the boundary condition of ψ(∂Ω I ) = 0 to all semiconductor-insulator interfaces, where with the symbol ∂Ω i we denote the surface of the region Ω i .We assume all metallic regions are described by ideal metals, and therefore the potential is homogeneous in these domains.When a voltage bias V i is applied to region Ω Mi , a Dirichlet boundary condition can be used that takes the form ϕ(∂Ω Mi ) = δµ i + V i where δµ i is a materialspecific constant modeling the Fermi energy difference between the metal in regions Ω Mi and the reference one.
When Ω Mi is a floating island, the Neumann condition ∇ ⊥ ϕ(∂Ω Mi ) = 0 fixes the electric field to be normal to the surface.Precise calibration of the band-offsets δ i requires a careful comparison of numerical simulations and experimental results.Since this work is a proof of concept, the offsets δµ i for metals are given arbitrary but reasonable values, and we refer to other references for the details of such procedure [35].The choice of material properties constants, including band offsets between semiconductors used to evaluate the conduction band minimum E CBM , are discussed in the Supplemental Material.Metal-semiconductor interfaces are a vast topic that is still actively being researched [36], and we have not found a convincing solution to the problem of including these interfaces in our method.We will first focus on the case where these interfaces are not present and assess the validity of the model while we postpone the discussion to Sec.IV.

III. CALIBRATION AND BENCHMARK
Before considering the ETF theory, we analyze the limits of the TF method (λ vW = 0) by simulating two simple but relevant geometries for 2DEG devices: a nanowire and a circular dot.For the 2DEG, we choose a semiconductor stack similar to the ones used in many modern experiments, e.g., Refs.[37][38][39].Fig. 1(a) display the schematic of a cross-section of a 2DEG nanowire.Two Au gates (yellow) serve the purpose of depleting the areas next to the Al wire (dark gray).The two Au gates and the Al wire are separated by HfO 2 dielectric (gray).Fig. 1(b) shows the electron density in the system, simulated using a top gate voltage of −3 V with respect to the grounded Al wire, while in Fig. 1(c) we plot the TF error as defined in Eq. ( 2).From this simulation, it is clear that the TF error does not satisfy R TF 1, and thus the electron density is varying too quickly in space to justify the use of the Thomas-Fermi approximation.This behavior is consistent when trying different top gate voltages.In Figs.1(d), (e), and (f), we simulate a quantum dot on the same semiconductor stack.We apply a voltage of −0.25 V to the outer gate and 0 V to the inner gate.In Fig. 1 e) and f), we show, respectively, the electron density in a plane located in the middle of the InAs layer and the TF error that approaches 1 at the boundary to the depleted region.
After having assessed the need for a more elaborated kinetic energy functional than the Thomas-Fermi one, we now consider the ETF method.We start with the problem of determining the optimal value of λ vW for the use case of electrostatic simulations of nanoelectronic devices.We considered a 2D translational invariant MOS device formed by a semiconductor heterostructure quantum well with an insulator layer and a metallic top-gate as shown in Fig. 2(a).We study the electrostatic problem with the TF and ETF method with various λ vW and compare the results with a simulation done with the self-consistent SP method.The density per unit area is shown in Fig. 2(b).
To assess quantitatively the optimal value for the von Weizsäcker parameter, we introduce two metrics: the absolute difference of the density per unit area where N = n dx is the total number of electrons per unit area, and the quantity that takes into account the difference in the shape of the density profile.The results are shown in Fig. 2(c).In general, we find that low values of λ vW in the interval [0.05, 0.2] provide the best agreement with SP results.We decided to elect λ vW = 1/9 as the standard parameter because of the theoretical works backing this choice [27,32].We simulated a quantum dot shown in Fig. 3 with the TF, ETF, and STF methods to evaluate the accuracy of the method in a more realistic situation.Here we fix the outer gate to −0.35 V and consider the number of electrons in the dot as a function of the inner gate voltage.The results can be seen in Fig. 3(c).Notice that the device under consideration does not show a dot-like behavior as the charge increases almost linearly with the gate voltage, as expected for a 2D system.This suggests that the lateral confinement induced by the gate system is not able to strongly confine the electrons.
Of the three methods, the TF method generally predicts the largest number of electrons, while the STF method predicts the lowest.The ETF results are intermediate between the two.Even though the ETF and STF methods do not overlap for all voltages, the two methods seem to have great compliance in the moderate filling regime.
One important difference between the TF and ETF methods is that the TF method predicts a steep jump in the differential capacitance C(V ) = dQ dV as the voltage increases.This happens as electrons at a particular voltage start accumulating in a previously classically forbidden region, while the ETF method, by allowing exponentially suppressed tails in the classically forbidden layers, prevents this from happening.This makes the ETF method less prone than the TF method to show unphysical behavior in semiconductor heterostructures.Thus the ETF FIG. 2. Calibration of λ vW .(a) shows a schematic of the simulated semiconductor stack inspired by the one used in [39] (not to scale).An Au gate (yellow) is separated from the four semiconducting layers by a layer of oxide (gray).(b) shows the results of TF, ETF (with various λ vW ), and SP simulations.(c) shows a comparison of different values of λ vW , using the two metrics defined in Eqs. ( 9) and (10).For all simulations the built-in bias at the interface is 0.1 V. method would, for instance, also be more appropriate to calculate capacitance compared to the TF method.In the setup considered here, we can thus conclude that the ETF method is superior to the TF method for simulating the number of electrons and calculating elements of the capacitance matrix.

IV. SEMICONDUCTOR-METAL BOUNDARY CONDITION
As is briefly described in the method presentation, the treatment of the interfaces between metals and semiconductors is an open problem [36].Moreover, the electron density shows a strong dependence on the thickness of thin metallic films that afflicts heterostructures when placed in contact with clean metallic films [2].Usually, the electrostatic problem of the charge distribution in the metal is not taken into account, and metallic parts are assumed neutral by assigning a Dirichlet boundary condition at the surface [1,2].We will consider the problem starting with a complete treatment of the metal-semiconductor interface and discuss the issues of this solution.Next, we will search for an approximate solution able to reproduce the important physical behavior of the electrostatics in the semiconductor.

A. Complete treatment
A complete treatment of the Sm-M interface can be formulated assuming that ψ 2 represents the conduction band electron density in both the semiconductor and the metal.We denote by E F,M the bulk Fermi energy of the metal defined as We define E W as the difference between the Fermi energy of the metal and the conduction band minimum of the semiconductor such that the local CBM takes the form where 1 Ωi is the indicator function of region Ω i .
In addition, since in the metal the Fermi energy lies in the conduction band, the equilibrium bulk electron density needs to be compensated by a positive background that we assumed homogeneous and equal to , (12) such that a homogeneous metallic system is neutral at equilibrium with an electron density equal to the one predicted by the Thomas-Fermi model.However, there are issues with this model when it is applied to a metalsemiconductor system.Since ψ 2 is the electron density in both materials, there will be a huge jump in its value at the interface, violating the assumption of slowly varying electron density.

B. Effective boundary treatment
Since a complete model of the semiconductor-metal interface is not available and the boundary conditions for the ETF method are unknown, we explored the possibility of finding an effective way to treat the semiconductormetal boundary.
A first approach consists in excluding the metal from the NLSE.In this way, ψ 2 models the electron density only in the semiconductor, and the problem reduces to finding an effective boundary condition by trial and error.We tested three options.If the interface is not transparent, setting the density at the interface to zero ψ(∂Ω M ) = 0 can be an acceptable boundary condition even in the metallic case.Alternatively, if we assume that the metal is not perturbed at all, we can impose that the density should continuously go to the unperturbed metal density at the interface √ n M .Finally, we tried a Neumann boundary condition by fixing the change in electron density to zero at the semiconductor-metal interface ∂ψ(∂Ω M ) = 0.
These approaches cannot be applied if the metal has to be included, for example, because it is a floating part.In this case, we found that a mixed TF and ETF approach can be used, where we solve for the electron density in both the semiconductor and the metal, but we allow λ vW to vary in space.The idea of promoting λ vW to an inhomogeneous field has already been considered in more traditional application fields of OF density functional theory [40,41].
Since metals have an extremely large electron density compared to that of the semiconductor, there will be an extremely large gradient of the electron density at the very interface.Since the vW-correction of the energy functional is proportional to the absolute value of this gradient squared, the vW-term will be extremely large here and thus essentially penalize the energy functional, trying to remove the abrupt change in electron density.However, this abrupt change in electron density at the Sm-M interface is what we would expect of the system and should thus not be removed.A solution to this could be to change λ vW through the stack to diminish the correction where we expect large gradients.We dubbed this method λ vW -sweep.In the metal, we expect an extremely large electron density (when compared to the semiconductor) that is only weakly perturbed by being in contact with a semiconductor.Changes are thus slow in space on the length scale of the Fermi wavelength, and the electron density in the metal can be effectively described by the TF method, i.e., assigning λ vW 0 in the metal while using a finite λ vW in the semiconductor.This circumvents the problem of assigning a boundary condition to the interface.In the bulk of the semiconductor, we use λ vW = 1/9, and in the aluminum we set λ vW to zero (for convergence reasons, we use a small but non-zero λ vW of 2 • 10 −8 ).
Simulation of the semiconducting stack in Fig. 2 (a), with a 5 nm layer of aluminum on top.The upper plot shows the electron density in the semiconducting stack, while the lower one shows the electron density in the metal.6 different ways of treating the semiconductor-metal interface are simulated: Three ETF simulations with different boundary conditions at the interface (ψ = 0, ψ = √ n Al , and ∂ψ = 0), one TF simulation, one SP simulation that is solved only in the semiconductor, and one ETF simulation where the value of λ vW is swept such that it is 1/9 in the semiconductor and 0 in the metal.
To test the method above, we considered a 2DEG system built with the semiconducting stack from Fig. 2 with a 5 nm layer of Al deposited on top as shown in Fig. 4. In the cases where the metal is included in the density part, we set the Fermi energy to E F,Al = 11.6 eV and neutralized the bulk charge with a positive background.With this value, the bulk electron density of Al is roughly n M = 1.8 × 10 29 m −3 .The electron density simulated for all the cases described can be seen in Fig. 4. For comparison, we included the result of a SP simulation.
In the top plot of Fig. 4, we see how the results for the ψ = 0 and ∂ψ = 0 boundary conditions are similar to the ones obtained with λ vW -sweep model, and how all of these methods seem to provide a good agreement with the SP method.The Dirichlet boundary condition ψ = √ n Al approach seems to give an erroneously small electron density in the well.This is caused by the extreme change in electron density that has to be neutralized entirely in the semiconductor.This causes an excess of negative charge near the interface that is left unbalanced.The λ vWsweep method makes it possible to avoid this problem since the electron density in the metal is now described by the TF method and, thus, no boundary condition is required at the semiconductor-metal interface.As can be seen, even if the depletion in the metal is localized in a tiny region at the interface, given the extremely high bulk density, the decrease causes a strong dipole at the interface that strongly changes the electrostatics restoring a result closer to the SP simulation.We can conclude that if the semiconductor layer closer to the metal is classically forbidden, both ψ(∂Ω M ) = 0 and ∂ψ(∂Ω M ) = 0 are acceptable effective boundary conditions for the orbitalfree theory.If the metal has to be explicitly included, a λ vW -sweep can be used for this purpose.We ruled out the ψ(∂Ω M ) = √ n M option as it does not provide comparable results.

V. GENERAL CONSIDERATION ON COMPUTATIONAL COMPLEXITY
Before concluding, we briefly discuss the convergence properties of ETF methods, and in general, OF methods, in comparison to TF and SP.Estimating the space and time complexity of these methods can be complicated, yet by assuming that the electric potential is discretized using N e degrees of freedom and the charge distribution on N q degrees of freedom, we can make some general observations.
The TF method is represented using only the electric potential, leading to a space complexity of O(N e ).The method requires solving a relatively simple nonlinear PDE for the electric potential, which generally requires iterative methods to converge, making it difficult to estimate the convergence rate.
In contrast, the ETF method is represented as a functional of both the electrostatic field and the matter field, leading to a space complexity of O(N e +N q ).The coupled PDEs that need to be solved are more complex than the one in the TF method.In some cases, it is possible to decouple the two PDEs and solve them separately and iteratively (using a so-called segregated or partitioned approach in contrast to a fully coupled/monolithic one) in order to optimize the convergence speed.
We assume that only k electrons are considered for the SP method, leading to a space complexity of O(N e +kN q ).The calculation is dominated by the solution of the Schrödinger equation, and assuming the use of an iterative solver like Arnoldi iteration, the time complexity is bounded from below by O(kN 2 q ).The cost of diagonalization is needed for each step of the iteration, together with the solution of the Poisson equation, which generally has a negligible cost compared to diagonalization.The solution strategy is inherently segregated, and this can cause potential instability problems that need to be addressed with special care [3].It is important to note that in the case of large systems, not only the number of degrees of freedom increase, but, in general, it is also necessary to increase the number of states k considered.
In summary, it is clear that the computational cost of OF methods increases more slowly with the size of the system compared to the SP method, making it a more attractive option for large-scale simulations.Moreover, the expression in variational form is amenable to the application of many convergence optimization techniques commonly used for general PDEs that are not easily applicable to diagonalization problems like multigrid methods, preconditioning, and in particular, domain decomposition.

VI. CONCLUSIONS
In this work, we investigated orbital-free methods for the solution of the electrostatic problem for nanoelectronic devices.We checked that the widespread Thomas-Fermi method, regarded as the simplest orbital-free method, is often applied outside its range of validity.The most notable case of the approximation breakdown occurs at the interfaces with classically forbidden regions, like insulators, where the density of electrons has abrupt jumps.
To achieve a more accurate description of the density profile in these cases, we considered the extended Thomas-Fermi (ETF) method that includes the von Weizsäcker correction of the kinetic energy functional.This correction can significantly increase the precision of the density profile near interfaces.
We addressed the question of the optimal value of the λ vW parameter by studying a simple 2DEG system and found that the theoretically motivated value of λ vW = 1/9 provides a good agreement also in the practical example cases considered.By applying the method to the simulation of realistic device geometries, we found that the ETF method provides density profiles closer to the ones calculated with the Schrödinger-Poisson method than the density profile provided by the TF method.SP methods can be computationally expensive as they require explicit diagonalization of the Hamiltonian, whereas the predictive power of the TF method is poor due to the perfect local behavior of the energy functional.Therefore, the ETF method represents a good compromise in terms of computational speed and predictive power.orbital-free methods are an often neglected alternative for electrostatic simulations of nanoelectronic devices, which are useful to handle large systems since the problem can be nicely expressed in variational form and implemented on any finite element solver.

VII. ACKNOWLEDGMENTS
W. S. and A. M. acknowledge Georg Winkler for his help and support through this work, and Andreas Pöschl, Alisa Danilenko, and Andrea Giuseppe Landella for useful discussion.The authors acknowledge Microsoft Quantum for support and computational resources.This work.was supported by the Danish National Research Foundation, the Danish Council for Independent Research |Natural Sciences.

Appendix A: General orbital-free functional theory
As mentioned in the main text, the theory proposed is not strictly a density functional theory as the electric potential appears as an argument of the energy functional, differently from what is common in electronic structure calculations.In principle, it is still possible to eliminate the electric potential from the formulation, but the presence of metallic gates makes the application of the procedure extremely complicated and cumbersome.
To illustrate the procedure, we can start with the generic orbital-free functional provided by Eq. ( 3) and proceed to find the minimizing condition for ϕ.We get This density functional theory formulation is exactly equivalent to the orbital-free theory that is the object of this paper.However, it relies on the fundamental solution G(r, r ) that is an extremely complicated object in any real case since in almost any case inhomogeneous permittivity and metallic gates with nontrivial shapes are present in the system.
For electronic structure calculations, this is not the case and indeed the PE equations have the simple fundamental solution With this simplification, we can write the total energy functional as a functional of ρ only as we have seen that the functional minimization with respect to ϕ directly gives Poisson's equation.This means that we can now minimize the functional with respect to ρ, and couple this with Poisson's equation to capture both the physical behavior from ϕ and ρ.Thus we have decoupled the problem and can now write the functional with respect to ρ only.This is exactly the form that we implicitly assume in the calculations above and, thus, the functional we use to derive the TF and ETF methods.This result is exact in free space, as it is under the assumption of ε being constant and under the assumption of no boundary conditions.However, this is not correct in any real device.

FIG. 1 .
FIG. 1. Electrostatic simulations of two example devices.The first row shows sketches of a cross-section of a common 2DEG device (a) and a simple circular quantum dot (d).In both cases, the devices are built on top of a semiconductor stack (schematic given in Fig.2) that provides the vertical confinement of electrons.The lateral confinement is controlled by Au gates (yellow), HfO 2 oxide (gray), and an Al wire (dark gray).b) and e) show the electron density from a TF simulation.e) is plotted in the middle of the InAs well.c) and f) show R TF as defined in Eq. (2).

FIG. 3 .
FIG. 3. Results from the dot-probe device simulation, using three different methods, i.e., TF, ETF, and STF.(a) shows the semiconducting stack used (not to scale), and (b) shows the geometry.(c) shows how electrons accumulate in the semiconducting stack as the inner gate voltage increases.