Many-particle Hamiltonian for open systems with full Coulomb interaction : Application to classical and quantum time-dependent simulations of nanoscale electron devices

A many-particle Hamiltonian for a set of particles with Coulomb interaction inside an open system is described without any perturbative or mean-field approximation. The boundary conditions of the Hamiltonian on the borders of the open system in the real three-dimensional 3D space representation are discussed in detail to include the Coulomb interaction between particles inside and outside of the open system. The manyparticle Hamiltonian provides the same electrostatic description obtained from the image-charge method, but it has the fundamental advantage that it can be directly implemented into realistic classical or quantum electron device simulators via a 3D Poisson solver. Classically, the solution of this many-particle Hamiltonian is obtained via a coupled system of Newton-type equations with a different electric field for each particle. The quantum-mechanical solution of this many-particle Hamiltonian is achieved using the quantum Bohm trajectory algorithm X. Oriols, Phys. Rev. Lett. 98, 066803 2007 . The computational viability of the manyparticle algorithms to build powerful nanoscale device simulators is explicitly demonstrated for a classical double-gate field-effect transistor and a quantum resonant tunneling diode. The numerical results are compared with those computed from time-dependent mean-field algorithms showing important quantitative differences.


I. INTRODUCTION
The exact computation of a system of interacting electrons is extremely complicated 1,2 because the motion of one electron depends on the positions of all others ͑i.e., electrons are correlated 3 ͒.Thus, the prediction of the collective behavior of many electrons is still a very active field of research in nanoelectronics, quantum chemistry, nanobiology, quantum computing, materials science, etc.Several theoretical approximations have been proposed to improve the treatment of electron-electron correlations.
In quantum systems in equilibrium, the time-independent mean-field approximation appears as a successful solution to treat a set of interacting electrons.It simplifies the exact many-particle potential by an average or mean potential 2 that transforms the many-body Schrödinger equation into a much more simple set of time-independent single-particle Schrödinger equations.The Hartree-Fock method 4,5 is a successful example of such approximation.However, by construction, the correlations among electrons can only be treated approximately.In principle, the density-functional theory 6,7 provides an exact path to deal with full electron correlations using single-particle potentials.However, since the exact form of the single-particle potentials [6][7][8][9] is unknown, an educated guess for these average single-particle potentials is needed in all practical algorithms.Therefore, again, the electron-electron correlations are treated approximately. 9,10he accurate treatment of the electron-electron correlations in electric circuits is even a more difficult issue [11][12][13][14][15][16][17][18][19][20][21][22][23][24][25][26] because we deal with nonequilibrium open systems [11][12][13] ͑where the system interchanges energy and particles with its environment͒.In fact, the Coulomb interaction among electrons is directly not considered in many quantum transport formalisms 11,12 under the assumption that the open system behaves as a Fermi liquid. 14The well-known Landauer approach 15,16 is a very successful example of the applicability of this assumption.Nevertheless, the Fermi-liquid paradigm has difficulties dealing with high-frequency, 11,17 low-dimensionality, 1,2 or Coulomb blockade regimes. 11,18On the other hand, the nonequilibrium Green's-function formalism ͑also referred to as the Keldysh formalism͒ provides an interesting path to solve the Schrödinger equation with the Coulomb interaction introduced perturbatively. 19Alternatively, under the assumption that the system behaves like a capacitor, one can use a simple linear relationship between the number of electrons and the electrostatic potential in a particular region to introduce partially Coulomb effects. 12,18he mean-field approximation appears again as a solution for electron transport.For example, an average single-particle time-independent potential profile can be computed, selfconsistently, from the set of wave-function solutions of a single-particle time-independent Schrödinger equation. 2,12his represents a zero-order approximation ͑sometimes called the Hartree 4 approximation͒ to the complex problem of electron-electron correlations.Additionally, remarkable efforts have been done by Büttiker and co-workers [20][21][22] to include Coulomb interaction in their scattering matrix approach by applying different many-body approximations to provide self-consistent electron-transport theories with overall charge neutrality and total current conservation.Finally, extensions of the equilibrium density-functional theory to deal with electron transport, by means of a time-independent formalism, 23 or with a powerful time-dependent version [24][25][26] can also be found in the literature.The exact exchangecorrelation functionals needed in both formalisms are unknown and they have to be approximated.Therefore, in all the descriptions of nonequilibrium quantum systems mentioned here, the electron-electron correlations are approximated to some extent.
For classical electron devices, the electrostatic interaction among electrons is commonly obtained from an explicit solution of the Poisson ͑Coulomb͒ equation.However, again, this does not provide an exact treatment of the classical electron-electron correlations but only an average estimation. 27,28It is well known that the solution of a classical many-particle system can always be written as a coupled system of single-particle Newton-type equations.However, a classical mean-field approximation is explicitly assumed in semiclassical transport simulators in order to deal with a unique average electrostatic potential for all electrons. 27A successful application of the classical mean-field approximation appears, for example, in the semiclassical Boltzmann equation that describes the time evolution of the electron distribution function in a one-electron phase space. 279][30][31][32] A powerful time-dependent technique to solve the Boltzmann equation is the semiconductor Monte Carlo method applied to electron devices. 27n this work, we are interested in revisiting the computation of an ensemble of Coulomb-interacting particles in an open system without any of the approximations mentioned in the previous paragraphs.With this goal, we have developed an exact many-particle Hamiltonian for Coulomb-interacting electrons in an open system in terms of the solutions of the Poisson equation.To our knowledge, the type of development of the many-particle Hamiltonian proposed in this paper has not been previously considered in the literature because, up to now, it was impossible to handle the computational burden associated with a direct solution of a many-particle Hamiltonian.Here, we present a classical and also a quantum solution of the many-particle Hamiltonian, both of which are applicable to realistic three-dimensional ͑3D͒ electron devices.The classical solution is obtained by solving a coupled system of Newton-type equations with a different electric field for each particle.The quantum solution of the many-particle Hamiltonian is obtained from the use of quantum trajectories. 33The merit of the quantum solution is certainly remarkable because, nowadays, the computational burden associated with the direct ͑i.e., without any approximation͒ solution of the many-particle wave function is only accessible for very few ͑2,3,...͒ electrons. 1,2Our quantum algorithm is able to treat electron dynamics without any ͑mean-field or perturbative͒ approximation in the description of the electrostatic interaction among a larger number ͑ϳ50͒ of electrons.In this paper, we present the classical and quantum algorithms together because they solve exactly the same many-particle Hamiltonian and both share many technical details ͑such a 3D Poisson solver to treat spatialdependent permittivity scenarios͒ in their implementation into realistic 3D electron devices.
After this introduction, the rest of the paper is divided as follows.In Sec.II, we write the many-particle Hamiltonian for an ensemble of electrons in an open system.We discuss the role of the boundary conditions on the borders of the open system to include the Coulomb interaction between particles inside and outside of the open system in Sec.III.In Sec.IV, we discuss the solution of the many-particle Hamil-tonian in classical scenarios.The path for the quantum solution is provided in Sec.V using quantum ͑Bohm͒ trajectories.In Sec.VI, we show the numerical results for the classical and quantum solutions of the many-particle Hamiltonian for nanoscale electron devices and we compare the results with time-dependent mean-field approximation.We conclude in Sec.VII.Appendixes A and B discuss the technical details of the image-charge method and mean-field approximation.

II. MANY-PARTICLE HAMILTONIAN IN OPEN SYSTEM
In this section, we develop an exact expression for the many-particle Hamiltonian that describes a set of electrons in an open system.Throughout this paper, we will assume that the magnetic field is negligible and that the particle velocity is small enough to assume a nonrelativistic behavior.In addition, in order to provide a discussion valid for either classical or quantum systems, we will assume spinless particles.Let us clarify that the exchange interaction is always present in a system of identical particles ͑electrons͒, but it will not be mentioned in this section because it does not affect explicitly the expression of the ͑first-quantization͒ many-particle Hamiltonians discussed here.The exchange interaction is introduced into the symmetry ͑when electron positions are in-terchanged͒ of the many-body wave function.We will briefly revisit this issue in Sec.V when dealing with the quantum solution.

A. Many-particle Hamiltonian for a closed system
First, we start our discussion with a set of M particles in a closed system.The many-particle Hamiltonian contains kinetic plus Coulomb energies,

͑1͒
The factor 1 2 that appears in the second term of the right-hand side is due to the fact that each two-particle interaction is counted twice ͓i.e., V͑r ជ k , r ជ j ͒ is identical to V͑r ជ j , r ជ k ͔͒.The con- dition j k takes into account the obvious restriction that a particle cannot interact with itself.The kinetic energy K͑p ជ k ͒ that appears in Eq. ͑1͒ is defined, for a classical system, as while for a quantum system

͑2b͒
Let us notice that the position and momentum in Hamiltonian ͑1͒ can be either classical variables, r ជ k and p ជ k in Eq. ͑2a͒, or quantum ͑real-space representation͒ operators, r ជ k and −iបٌ r ជ k in Eq. ͑2b͒.In particular, it is important to emphasize that when we refer to r ជ k as the electron position, we are not referring to a fixed position but a variable vector.On the contrary, when we are interested in specifying a fixed electron position, we will write r ជ k ͓t͔.The parameter m k is the particle mass that, in Sec.VI, will be understood as the particle effective mass.Identically, we define the Coulomb potential in Eq. ͑1͒ as where q j is the particle charge and is the permittivity.Although we are always thinking about electrons in semiconductors, the development of this section is valid for arbitrary particles with different masses and charges.
A complete electronic circuit ͑including the devices, the wires, and the batteries͒ behaves as a closed system with a large ͑M → ϱ͒ number of electrons.However, since we can only deal with a finite number of electrons, we restrict our system to a small part of the circuit, for example, the channel of a transistor.Thus, we need to develop the Hamiltonian that describes the dynamics of a subensemble of the whole set of M particles in an open system inside a limited volume ⍀ ͑see Fig. 1͒.

B. Many-particle Hamiltonian for an open system
We divide the previous ensemble of M particles into a subensemble ͕1,2,3, ... ,N͑t͖͒ of particles whose positions are inside the volume ⍀ and a second subensemble ͕N͑t͒ +1, ... , M͖ which are outside it ͑see Ref. 34͒.We assume that the number of particles inside, N͑t͒, is a time-dependent function that provides an explicit time dependence in the many-particle ͑open-system͒ Hamiltonian.As drawn in Fig. 1, we assume a parallelepiped where the six rectangular surfaces S = ͕S 1 , S 2 , ... ,S 6 ͖ are the boundaries of ⍀.Along this paper, we use r ជ l as the "boundary" vector representing an arbitrary position on the surfaces S l .
Since we are only interested in the dynamics of the first N͑t͒ particles, the kinetic energy and the Coulomb interaction between the particles of the second subensemble do not appear in the new Hamiltonian of the open system.Nevertheless, the Coulomb interaction between particles of the first and second subensembles must explicitly appear.Thus, the many-particle Hamiltonian for the first N͑t͒ particles can be written as

͑4͒
Let us notice also that the factor 1 2 disappears in the last term of Eq. ͑4͒ because there is no double counting of interactions between electrons inside and outside ⍀.For convenience, we rewrite Eq. ͑4͒ as

͑5͒
Although up to this point we have discussed the manyparticle Hamiltonian in terms of the Coulomb force, this approach is inconvenient to deal with solid-state scenarios with a spatial-dependent permittivity. 35,36For this reason, we rewrite our many-particle Hamiltonian in terms of the more generic Poisson ͑or Laplace͒ equation, which can be applied to systems with ͑or without͒ a spatial-dependent permittivity ͓by substituting → ͑r ជ͒ in the Poisson equation͔.Each term V͑r ជ k , r ជ j ͒ that appears in Eq. ͑5͒ can be explic- itly obtained from a Poisson ͑or Laplace͒ equation inside the volume ⍀.Then, using the superposition property of the Poisson ͑or Laplace͒ equations, we can rewrite Eq. ͑5͒ as H͑r ជ 1 , ... ,r ជ N͑t͒ ,p ជ 1 , ... ,p ជ N͑t͒ ,t͒

͑6͒
where the term W k ͑r ជ 1 , ... ,r ជ k , ... ,r ជ N͑t͒ ͒ is a particular solution of the following Poisson equation: The term k ͑r ជ 1 , ... ,r ជ N͑t͒ ͒ in Eq. ͑7͒ depends on the position of the first N͑t͒ electrons, 34 but Eq. ͑8͒ is independent of the position of the external particles because they only affect the boundary conditions of Eq. ͑7͒.Let us notice that there are still terms, V͑r ជ k , r ជ j ͒, in Eq. ͑6͒ that are not computed from Poisson equation but from Eq. ͑3͒.However, we will show that these terms V͑r ជ k , r ជ j ͒ have no role in the classical ͑i.e., Sec.IV͒ or quan- tum ͑i.e., Sec.V͒ solutions of Eq. ͑6͒.
By construction, comparing Eqs.͑5͒ and ͑6͒, the term W k ͑r ជ 1 , ... ,r ជ N͑t͒ , t͒ can be rewritten as The dependence of W k ͑r ជ 1 , ... ,r ជ N͑t͒ , t͒ on the positions of the external particles is explicitly written in the last sum in Eq. ͑9͒, while in Eq. ͑7͒ this dependence is hidden in the boundary conditions of W k ͑r ជ 1 , ... ,r ជ k , ... ,r ជ N͑t͒ ͒ on the surface S = ͕S 1 , S 2 , ... ,S 6 ͖.In fact, the boundary conditions are a delicate issue that we will discuss in Sec.III.Finally, we want to remark that this discussion is valid for either classical or quantum Hamiltonians because the expression ͑9͒ ͓or its equivalent definition in ͑7͒ and ͑8͔͒ of W k ͑r ជ 1 , ... ,r ជ N͑t͒ , t͒ at r ជ 1 , ... ,r ជ k , ... ,r ជ N͑t͒ is identical for a classical or a quantum system.

ON THE BORDERS OF THE OPEN SYSTEM
Since we want to deal with solutions of the Poisson equation ͓Eq.͑7͔͒, the boundary conditions for the N͑t͒ terms W k ͑r ជ 1 , ... ,r ជ k , ... ,r ជ N͑t͒ ͒ have to be specified on the border surface S = ͕S 1 , S 2 , ... ,S 6 ͖.Such boundary conditions will provide, somehow, information on the electrostatic effect that outside electrons ͓i.e., N͑t͒ +1, ... , M͔ have on the electrons inside ⍀.In order to provide a clear notation for discussing the boundary conditions of W k ͑r ជ 1 , ... ,r ជ k , ... ,r ជ N͑t͒ ͒, we distinguish between the "source" vectors r ជ 1 , ... ,r ជ k−1 , r ជ k+1 , r ជ N͑t͒ and the additional "observation" vector r ជ that runs over all space. 36In particular, the electrostatic potential that appears in Hamiltonian ͑6͒ is defined as the value of the potential Our goal is to find an educated guess for all the N͑t͒ terms W k ͑r ជ 1 , ... ,r ជ k−1 , r ជ , r ជ k+1 , r ជ N͑t͒ , t͒ at all observation points r ជ = r ជ l on all surfaces l = 1 , . . ., 6.For example, the information of such boundary conditions can come from the value of the total voltage ͑due to internal and external electrons͒ at position r ជ l and time t.We define the total voltage B͑r ជ 1 , ... ,r ជ N͑t͒ , ... ,r ជ M , r ជ l , t͒ as the electrostatic potential asso- ciated to an additional probe charge q M+1 situated on that boundary, r ជ l ϵ r ជ M+1 S l ͑see Fig. 2͒.The electrostatic poten- tial "seen" by this extra charge due to the presence of the rest of the particles is just

͑11͒
where the expected restriction j M + 1 is hidden in the limit of the sum.Once relationship ͑11͒ is established, we can easily define the boundary conditions of any of the N͑t͒ electrostatic potential In particular, from Eq. ͑9͒, we know that The discussion done here is valid for either classical or quantum systems ͑see Ref. 37͒.In the previous discussion we have assumed Dirichlet boundary conditions; however it is a straightforwardly procedure to develop the same argumentations for Neumann boundary conditions.
The reader can be surprised by the fact that the right-hand side of expression ͑12͒ tends to infinite, V͑r ជ l , r ជ k ͒ → ϱ, when r ជ k → r ជ l .However, when r ជ k → r ជ l , the extra particle at r ជ l ϵ r ជ M+1 S l will also provide an infinite value of the electro- static potential, i.e., B͑r ជ 1 , ... ,r ជ N͑t͒ , ... ,r ជ M , r ជ l , t͒ → ϱ, due to the presence of the k particle on the surface.Therefore, the first infinite, V͑r ជ l , r ជ k ͒ → ϱ, is canceled by the second infinite, B͑r ជ 1 , ... ,r ជ N͑t͒ , ... ,r ជ M , r ជ l , t͒ → ϱ.This discussion will be rel- evant in Sec.VI when we discuss the numerical implementation of these boundary conditions. .͑Color online͒ The electrostatic potential B͑r ជ 1 , . . .,r ជ N͑t͒ , . . .,r ជ M , r ជ l , t͒ measured on a surface S l at position r ជ l and time t ͑due to internal and external electrons͒ by an additional probe charge q M+1 situated on the boundary r ជ l ϵ r ជ M+1 S l .
Up to here, our argumentation might seem somehow tricky.We have defined the value of W k ͑r ជ 1 , ... ,r ជ , ... ,r ជ M ͉͒ r ជ=r ជ l on the volume boundaries without mentioning the position of the external particles but using the function B͑r ជ 1 , ... ,r ជ N͑t͒ , ... ,r ជ M , r ជ l , t͒ which is unknown.This strategy transforms the complexity of finding the boundary conditions of N͑t͒ electrostatic potential W k ͑r ជ 1 , ... ,r ជ , ... ,r ជ N͑t͒ ͒ into pro- viding an educated guess for a unique function B͑r ជ 1 , ... ,r ជ N͑t͒ , ... ,r ជ M , r ជ l , t͒.In our numerical results in Sec.VI, we will fix B͑r ជ 1 , ... ,r ជ N͑t͒ , ... ,r ជ M , r ជ l , t͒ based on standard arguments for electron devices.We will assume a uniform value of the voltage on the l surface independent of the external electrons, B͑r ជ 1 , ... ,r ជ N͑t͒ , ... ,r ជ M , r ជ l , t͒ Ϸ B l ͑r ជ 1 , ... ,r ជ N͑t͒ , t͒.Such value can be obtained taken into account the voltage fixed by the external battery and the requirement of charge neutrality at the contacts. 35,38inally, we want to enlighten the physical interpretation of the many-particle Hamiltonian ͑6͒ and boundary conditions of Eq. ͑12͒ developed here.To do this, we compare our approach with the image-charge method applied to electron transport.The image-charge method is a basic solving tool in electrostatics 36 that has been successfully applied, for example, in the calculation of the electric field in field-emission devices 39 or the barrier-reduction in the metal-semiconductor Schottky contacts. 40The name of the method originates from the replacement of certain "real" charges by a set of few "imaginary" charges that replicate the real boundary conditions at the surface ͑see Fig. 3͒.From the uniqueness theorem of electrostatics, 36 once the charge of the 1 , . . ., N͑t͒ particles inside a volume is fixed and the correct electrostatic potential ͑or electric field͒ is specified at the boundaries of that volume, B͑r ជ 1 , ... ,r ជ N͑t͒ , ... ,r ជ M , r ជ l , t͒, the solution of the Poisson equation inside the volume is unique and does not depend on whether the external charges are real or imaginary.Then, according to the image-charge method, the electrostatic potential seen by the k particle is equal to the electrostatic potential generated by the sum of the imaginary plus the real particles except the k particle.Thus, identically to our many-particle Hamiltonian, the image-charge method goes beyond the mean-field approximation ͑discussed in Appendix B͒ because each particle feels its own potential profile that excludes the Coulomb interaction with itself.In Appendix A, we show in detail that the boundary conditions in Eq. ͑12͒ are identical to the boundary conditions found with the image-charge method.
However, although the outcome of the image-charge method and our many-particle Hamiltonian are identical, the Hamiltonian presented in this paper has an unquestionable advantage over the image-charge method: the former can be directly implemented into practical 3D electron device simulators as we will see in Sec.VI while the latter cannot.For example, let us consider the numerical simulation of the transistor done in Sec.VI.The system consists in a large number ͑ϳ20͒ of electrons inside a volume ⍀ limited by surfaces ͕S 1 , S 2 , ... ,S 6 ͖ with Dirichlet and Neumann boundary conditions.Then, the exact application of the image-charge method faces up to the following serious difficulties.The computation of the imaginary charges in an arbitrary surface ͑different from the standard infinite plane whose imaginary charges are found in textbooks 41 ͒ is not at all obvious. 42Let us notice that each imaginary charge that provides the correct value of the Neumann boundary condition on S i does also affect the Neumann ͑or Dirichlet͒ boundary condition on all other surfaces ͕S 1 , S 2 , ... ,S 6 ͖.Finally, even after assuming that we would be able to find somehow the density distribution of imaginary charges that reproduces simultaneously the boundary conditions on all six surfaces, the practical application of this method in a time-dependent simulator with a 3D Poisson solver ͑to be able to deal with spatial-dependent permittivity scenarios͒ would require simulating much more particles ͑the real plus the imaginary͒ in a larger simulation box ͑to include the location of those imaginary particles outside of ⍀͒.In summary, the image-charge method is an excellent approach to obtain analytical expression for the many-body description of electron transport in simple systems ͑such as one electron crossing an infinite ideal 41 metallic surface͒, but it is not practically possible to implement it in simulators for actual 3D electron devices. 42On the contrary, as we will show in our numerical result in Sec.VI, the many-particle Hamiltonian ͑6͒ and the boundary conditions of Eq. ͑12͒ can be implemented in a extremely simple and transparent way in, either classical or quantum, realistic electron device simulators using a 3D Poisson solver for arbitrary surfaces.

IV. EXACT MANY-PARTICLE HAMILTONIAN FOR CLASSICAL OPEN SYSTEMS
In this section, we discuss the classical solution of the many-particle open-system Hamiltonian of expression ͑6͒.Interestingly, the results obtained here can partially be used for the quantum solution of the many-particle Hamiltonian developed in Sec.V.
The classical description of the particle dynamics subjected to the many-particle Hamiltonian ͑6͒ can be computed by using the well-known Hamilton's equations.In particular, we can obtain the ͑Newton-type͒ description of the classical trajectory r ជ i ͓t͔ in the real space through For the many-particle Hamiltonian studied in this work, expression ͑13b͒ gives the trivial result m • v ជ i ͓t͔ = p ជ i ͓t͔, while the evaluation of expression ͑13a͒ requires a detailed development.We know that the r ជ i gradient of the exact many-particle Hamiltonian ͑6͒ can be written as We define the multidimensional vector R ជ = ͑r ជ 1 , ... ,r ជ N͑t͒ ͒ to account, in a compact way, for the classical trajectories of N͑t͒ electrons, R ជ ͓t͔ = ͑r ជ 1 ͓t͔ , ... ,r ជ N͑t͒ ͓t͔͒.Substituting the defi- nition of W k ͑r ជ 1 , ... ,r ជ N͑t͒ , t͒ done in expression ͑9͒ into Eq.͑14͒, we find Note the elimination of the factor 1 2 in the last term of the right-hand part of Eq. ͑15͒ that accounts for those terms q k • V͑r ជ k , r ជ j ͒ in Eq. ͑14͒, where r ជ k r ជ i and r ជ j = r ជ i that are iden- tical to the term q i • V͑r ជ i , r ជ k ͒ in Eq. ͑15͒.For the same reason, we include a factor 2 on the first term of right hand of Eq. ͑15͒.From expressions ͑9͒ and ͑15͒, we realize that Only the term W i ͑r ជ 1 , ... ,r ជ N͑t͒ ͒ of whole Hamiltonian ͑6͒ be- comes relevant for a classical description of the i particle.In fact, since we only evaluate a r ជ i gradient, the rest of particle positions can be evaluated at their particular value at time t, i.e., r ជ k → r ជ k ͓t͔ for all k i.Therefore, we define the single- particle potential W ¯i͑r ជ i , t͒ from the many-particle potential as

͑17͒
We use a "hat" to differentiate the ͑time-dependent͒ singleparticle electrostatic potential from the many-particle poten-tial.Each i term of the single-particle electrostatic potential, W ¯i͑r ជ i , t͒, is a solution of one particular 3D-Poisson equation, where the single-particle charge density is defined as and the boundary conditions ͓Eq.͑12͔͒ are adapted here as where we have included the approximation B͑r ជ 1 , ... ,r ជ N͑t͒ , ... ,r ជ M , r ជ Ј, t͒ϷB l ͑r ជ 1 , ... ,r ជ N͑t͒ , t͒ mentioned in Sec.III.Let us remind that expressions ͑17͒-͑20͒ provide an exact treatment of the many-particle correlations in classical scenarios.They imply a coupled system of Newton equations.The N͑t͒ Newton equations are coupled by N͑t͒ Poisson equations.Therefore, as mentioned in Sec.I, the manyparticle Hamiltonian of Eq. ͑6͒ can be written exactly ͑without mean-field approximation͒ as a sum of singleparticle Hamiltonian for classical scenarios,

͑21͒
The term W ¯k͑r ជ k , t͒ in Eq. ͑21͒ means that each particle "sees" its own electrostatic potential ͑or electric field͒, which is different from that of others.

V. EXACT MANY-PARTICLE HAMILTONIAN FOR QUANTUM OPEN SYSTEMS
The many-particle open-system Hamiltonian developed in expression ͑6͒ is also valid for quantum systems.In this section, we will explain its practical quantum solution using the recent quantum ͑Bohm͒ trajectory formalism in Ref. 33.For convenience, we rewrite the many-particle Hamiltonian in Eq. ͑6͒ as

͑22͒
where we explicitly write the electron momentum as p ជ k =−iបٌ r ជ k in the kinetic energy ͓as mentioned in Eq. ͑2b͔͒.

͑24͒
The practical utility of expression ͑24͒ in understanding quantum scenarios can seem quite doubtful because its direct solution becomes computationally inaccessible for more than very few electrons. 1,2,43However, one of the authors has recently developed a transport formalism 33 in terms of Bohm trajectories that simplifies the complexity of evaluating Eq. ͑24͒.
Some introductory explanations about Bohm trajectories in single-particle and many-body scenarios can be found in Refs.44-47.Here, we go directly to the main result of Ref. 33 where it is shown that a many-particle electron Bohm trajectory r ជ a ͓t͔ computed from the many-particle wave func- tion, ⌽͑r ជ 1 , ... ,r ជ N͑t͒ , t͒, solution of the Eq.͑24͒ can be equivalently computed from the single-particle wavefunction ⌿ a ͑r ជ a , t͒ solution of the following single-particle Schrödinger equation:

͑25͒
where we have defined R ជ a ͓t͔ = ͕r ជ 1 ͓t͔ , r ជ a−1 ͓t͔ , r ជ a+1 ͓t͔ , r ជ N ͓t͔ , t͖ as a vector that contains all Bohm trajectories except r ជ a ͓t͔.The exact definition of the other potentials that appear in Eq. ͑25͒, G a ͑r ជ a , R ជ a ͓t͔ , t͒ and J a ͑r ជ a , R ជ a ͓t͔ , t͒, can be obtained from Ref. 33.The total many-particle electrostatic potential in Eq. ͑24͒ has been divided into two parts, From expressions ͑9͒ and ͑23͒, we realize that U a ͑r ជ a , R ជ a ͓t͔ , t͒ can be greatly simplified as The rest of the terms V͑r ជ j ͓t͔ , r ជ i ͓t͔͒ of expression ͑26͒ appear in U b ͑R ជ a ͓t͔ , t͒ and they are included in the potential G a ͑r ជ a , R ជ a ͓t͔ , t͒.However, this term U b ͑R ជ a ͓t͔ , t͒ has no role on the single-particle wave function ⌿ a ͑r ជ a , t͒ because it has no dependence on r ជ a and it only introduces an irrelevant global phase on ⌿ a ͑r ជ a , t͒.
Let us notice that, in the right-hand side of expression ͑27͒, we have used the same definition of the potential profile as in classical expression ͑17͒.The only difference here is that R ជ a ͓t͔ are not classical trajectories but quantum ͑Bohm͒ trajectories.Therefore, the computation of the potential profile W ¯a͑r ជ a , R ជ a ͓t͔ , t͒ that appears in the single-particle Schrödinger equation ͓Eq.͑25͔͒ just needs 3D Poisson equations ͓Eqs.͑18͒ and ͑19͔͒ with the boundary conditions ͓Eq.͑20͔͒.Interestingly, since the term W ¯a͑r ជ a , R ជ a ͓t͔ , t͒ is computed from a Poisson equation, our quantum-trajectory algorithm can also be directly extended to spatial-dependent permittivity systems.
Finally, according to the summary done in Sec.I of this work, we want to emphasize the similarities between the ͑open-system͒ Bohm-trajectory computational algorithm discussed in Ref. 33 and the ͑electron-transport version͒ density-functional theory mentioned in Sec.I.For the latter, the decomposition of the many-particle system into a set of coupled single-particle Schrödinger equations is exact and demonstrated by the Hohenberg-Kohn-Sham theorem. 6,7owever, from a practical point of view, the exact exchangecorrelation functionals that appear in the single-particle Schrödinger equations are unknown and, therefore, require some approximation.Identically, in the quantum ͑Bohm͒ trajectory algorithm, the use of single-particle Schrödinger equation ͓Eq.͑25͔͒ is exact to treat many-particle system as demonstrated in the theorem considered in Ref. 33.However, the exact values of the terms G a ͑r ជ a , R ជ a ͓t͔ , t͒ and J a ͑r ជ a , R ជ a ͓t͔ , t͒ that appear in Hamiltonian ͑25͒ are unknown ͑because they require the partial knowledge of the shape of the many-particle wave function͒.Thus, from a practical point of view, they need to be approximated by some educated guess as in the density-functional theory.
Finally, the exchange interaction among the ͑fermions͒ electrons can also be considered in the present quantum algorithm.A brief explanation of how the exchange interaction can be introduced in the present quantum ͑Bohm͒ trajectory algorithm is mentioned in Ref. 33, but the discussion of this issue is far from the goal of the present paper.

VI. NUMERICAL RESULTS FOR THE MANY-PARTICLE HAMILTONIAN IN CLASSICAL AND QUANTUM ELECTRON DEVICES
In Secs.I-V of this paper, a many-particle Hamiltonian has been developed for an arbitrary ensemble of particles with Coulomb interaction among them.In this section, our numerical examples will deal with electrons in solid-state semiconductors.Therefore, first of all, we have to specify under what approximations the theoretical many-particle Hamiltonian developed in the first part can be used to describe solid-state semiconductors.Only the dynamics of the free electrons will be studied in our numerical results.The interaction with the rest of the charges ͑associated to core electrons and ions͒ will be considered as average polarization charges via a spatial-dependent permittivity. 36We do also assume an effective mass 48 for the free electrons that accounts for their interaction with the periodic atomic structure under the standard Born-Oppenheimer approximation 49 ͑that neglects the interaction of valence electrons with other kind of particle such as phonons͒.These are reasonable approximations in most electron-transport models of ballistic devices. 12,27e solve the many-particle ͑open-system͒ Hamiltonian from expression ͑6͒ to compute the current-voltage characteristic for classical and quantum electron devices.We use the classical algorithm discussed in Sec.IV for the simulation of a double-gate field-effect transistor 50 ͑DG-FET͒, while we use the quantum algorithm discussed in Sec.V for a resonant tunneling diode 47 ͑RTD͒.As mentioned above, no phonon, impurity, or roughness scattering mechanism is included in the simulations and only the full Coulomb interaction is considered explicitly.In order to emphasize the importance of our treatment of the electron-electron correlations in such nanoscale devices, we will compare these current-voltage characteristics with the results obtained with a time-dependent "mean-field" approach that will be discussed in Appendix B. We refer to "many-electron" results to describe the simulation done with either classical ͑Sec.IV͒ or quantum ͑Sec.V͒ algorithms that requires solving N͑t͒ Poisson equations with N͑t͒ different boundary con-ditions and charge densities ͓expressions ͑17͒-͑20͔͒ at each time step.Alternatively, we refer to the time-dependent ''mean-field'' results when a single Poisson equation ͓expressions ͑B1͒-͑B4͔͒ is solved for all electrons at each time step of the simulation.
For all simulations ͑quantum, classical, mean field, or many electron͒, the same electron injection model is used.We use an injection model applicable to systems with arbitrary electron confinement, which is a time-dependent version of the Landauer boundary conditions, valid for degenerate and nondegenerate systems.We inject electrons according to the Fermi-Dirac statistics defined by a Fermi level ͑an electrochemical potential͒ deep inside the contacts. 38The applied bias provides a difference between the values of the Fermi level at each injecting surface.Our injection model, coupled to the boundary conditions of the Poisson equation, also ensures charge neutrality at the contacts. 35ll simulations use a 3D finite-difference Poisson solver scheme.The whole volume ⍀ of the active region drawn in Fig. 1  with the educated guess that the component of the electric field normal to the surface is zero there, E ជ l ͑r ជ l , t͒ • n ជ l = 0, where n ជ l is a unit vector normal to the mentioned surfaces and pointing outward.The continuity of the displacement vector normal to surfaces justifies the assumption E ជ l ͑r ជ l , t͒ • n ជ l =0 at the boundaries when the relative permittivity inside ⍀ is much higher than the corresponding value outside.On the other hand, in the surfaces S 1 and S 4 of Fig. 1 we use the Dirichlet boundary conditions discussed in Sec.III, B l ͑r ជ 1 , ... ,r ជ N͑t͒ , t͒, with final expression ͑20͒.
Finally, a technical remark about the application of expression ͑20͒ in the classical or quantum many-electron simulations is mandatory.Strictly speaking, our assumption that the potential at one particular surface is position independent, B͑r ជ 1 , ... ,r ជ N͑t͒ , ... ,r ជ M , r ជ l , t͒ϷB l ͑r ជ 1 , ... ,r ជ N͑t͒ , t͒, is not completely accurate because we known from the discussion in Sec.III that the original function B͑r ជ 1 , ... ,r ជ N͑t͒ , ... ,r ជ M , r ជ l , t͒ in expression ͑11͒ has to repro- duce, somehow, the atomistic charge distribution on the surface. 41n particular, one can expect B͑r ជ 1 , ... ,r ជ N͑t͒ , ... ,r ជ M , r ជ l , t͒ → ϱ when the electron is close to the border, r ជ k → r ជ l .However, due to our ignorance about the atomistic description of the contact interface, 41 we apply the boundary conditions ͓Eq.͑20͔͒ assuming that the distance between r ជ k and r ជ l is always greater than 1 nm ͑this value is interpreted as a measure of range of the atomistic pseudopotential 2 in the spatial-dependent permittivity scenarios discussed here͒.

A. Classical simulation of two-electron system: Many electron versus mean field
In this section we will explain the origin of the important differences that will appear later between the mean-field and many-electron algorithms using a simple semiclassical twoelectron system.We consider one electron ͑labeled as one electron͒ injected from the source surface, S 4 , at an arbitrary position.A second electron is injected, arbitrarily, from the drain surface, S 1 .A battery provides an external voltage equal to zero at the drain and source surfaces.A 3D cubic system with a volume of ⍀ = ͑20 nm͒ 3 is considered as the active device region.We consider silicon parameters for the numerical simulation.Within the mean-field approximation only the potential profile W ¯mean ͑r ជ , t͒ is calculated for the two- electron system using expressions ͑B1͒-͑B4͒.Then, we realize from Fig. 4 that each electron can be reflected by an artificial alteration of the potential profile related to its own charge.In Figs. 5 and 6 we have plotted the energy potential profile seen by the one electron, W ¯1͑r ជ 1 , t͒, and by the two electron, W ¯2͑r ជ 2 , t͒, using the many-electron algorithm de- scribed by expressions ͑17͒-͑20͒.Electrons are not affected by their own charge.We clearly see that, within the meanfield approximation, electrons can be unable to overcome the large potential barrier that appears at their own position ͑due to their own charge͒.In addition, the simple results confirm that the mean-field error is equal to expression ͑B7͒, i.e., the error of the mean-field potential profile at each position of the active region is Error k ͑r ជ , t͒ = V͑r ជ , r ជ k ͓t͔͒.
Finally, a discussion about the role of the spatial mesh used for the numerical solution of the Poisson equation is relevant.For an electron device with a length of hundreds of nanometers, we need a mesh of the 3D active region with spatial step DX ϳ DY ϳ DZ Ͼ 10 nm to deal with not more than 1000 nodes in the numerical solution of the Poisson equation.This computational limitation in the resolution of the potential is present when solving either the mean-field or the many-electron algorithm.With such spatial resolution, the short-range interaction is missing in both procedures because two electrons inside the same spatial cell will not repel each other.In addition, the error between both procedures, Error k ͑r ជ , t͒ = V͑r ជ , r ជ͓t͔ k ͒, is reduced because the numerical Coulomb potential profile is smoothed due to the low resolution ͓i.e., the diameter of the region where V͑r ជ , r ជ͓t͔ k ͒ has a strong influence is shorter than the cell dimensions͔.Therefore, we obtain roughly identical results with both schemes.In the subplots of Fig. 7, the same electron trajectory is presented for different mesh resolutions.As can be seen, for the best mesh resolution ͑DX = DY = DZ =2 nm͒, the differences between both treatments are maximized due to the important spurious autoreflection effect found in the mean-field trajectory.On the other hand, as the resolution of our mesh is reduced, differences between both treatments disappear, giving roughly equal trajectories for cell dimensions above 5 nm.
In summary, when the spatial cells are large, the meanfield and the many-electron schemes correctly model the long-range Coulomb interaction, but both neglect the shortrange component.On the contrary, with smaller spatial steps DX ϳ DY ϳ DZ Ͻ 5 nm, the many-electron resolution takes into account long-and short-range Coulomb interactions correctly, whereas the description of the short-range component within the mean-field approximation is completely incorrect ͑i.e., electrons are repelled by themselves͒.In other words, when DX , DY , DZ → 0 the mesh error in our many-electron algorithm reduces to zero, while the error in the mean-field approach tends to infinite, Error k ͑r ជ , t͒ → ϱ ͑see a schematic summary of the explanation of this discussion in Fig. 8͒.Finally, it is important to remark that electron trajectories in Fig. 7 are computed using the classical scheme of Sec.IV, but the electrostatic potential profiles are computed from a 3D Poisson solver that is identical for the classical or quantum algorithms.Therefore, the conclusions drawn here for the classical algorithm can be directly extrapolated to our quantum algorithm.In the classical algorithm, the wrong potential profile of Fig. 4 affects the electric field ͓Eqs.͑13a͒ and ͑13b͔͒ that modifies the electron dynamics.Identically, the wrong mean-field potential in expression ͑25͒ will affect the solution of the Schrödinger equation that will modify Bohm trajectories.

B. Classical simulation of a double-gate field-effect transistor
Now, we use the classical solution of the many-particle Hamiltonian ͑6͒ to provide a full simulation 51 for the DG-FET depicted in Fig. 9. Electron transport in the x direction ͑from source to drain͒ takes place along a silicon ͑100͒ orientation channel at room temperature.In particular, the electron mass is taken according to the six equivalent ellipsoidal constant energy valleys of the silicon band structure. 27,50The effective masses of the ellipsoid are m l = 0.9163m 0 and m t = 0.1905m 0 , with m 0 as the free-electron mass.For details on the particular effective mass taken by the electrons in each direction and valley see Ref. 38.The dimensions of the channel of devices L z and L y are both small enough, so that the active region becomes an effective one-dimensional ͑1D͒ system ͑a quantum wire͒ and the energy of an electron in one particular valley is represents the minimum energy of the first subband, whose value is E 1D = 0.182 eV for L z = 2 nm and L y = 5 nm.The energies of the next lowest subbands ͑E 1D q = 0.418 eV or E 1D q = 0.489 eV͒ are assumed high enough to keep a single band simulation.Therefore, we use a 3D Poisson solver for electrostatic, but a 1D algorithm to describe the velocity of each electron in the x direction.Due to the lateral electron confinement, the velocities in the y and z directions are zero. 52We solve N͑t͒ 1D Newton equation ͓Eq.͑13͔͒ coupled by N͑t͒ 3D Poisson equation ͓Eq.͑18͔͒.
We obtain the transistor current-voltage characteristics by computing the time evolution of many interacting electrons inside the 1D DG-FET.The classical many-electron algorithm is compared with the classical mean-field one.The details of the simulation are described in Table I.A total number of cells, Nx • Ny • Nz, on the order of 1000 and a number of electrons, N͑t͒, about 20-50, implies a simulation time on the order of 3-4 h for each bias point, 53 while it takes 20-30 min within the mean-field approximation.In Figs. 10 and 11 average values for current and number of particles are presented for different gate and drain voltages.When the FET remains under subthreshold region, the results are quite similar for both methods.However, interestingly, the average current of the FET system in the subthreshold region predicted by the many-electron algorithm is slightly larger than the result obtained by the mean-field approximation.In other words, the mean-field results remain in the subthreshold region, while the many-electron results show a DG-FET channel partially opened.In any case, the most important differences occur for higher gate voltages.In order to understand the results, we have to remind that the DG-FET tends to behave as a capacitor where the charge inside the channel is controlled by the gate voltage.In addition, the charge at the contacts is controlled by the injection process that achieves local charge neutrality there.Therefore, the number of electrons inside the channel tends to be identical within both methods.However, the average current that is sensible to electron dynamics is higher with the many-electron method than with mean-field approximation because fewer electrons are reflected in the former ͑i.e., there are no electrons reflected by its own charge͒.For the highest gate voltages, equal results for the mean current are obtained for both methods.

C. Quantum simulation of a resonant tunneling diode
In this section, we will provide a numerical example of the solution of the quantum many-particle Hamiltonian ͑6͒ for an ensemble of electrons in a RTD of Fig. 12.We again compare our many-electron method with the mean-field approximation.We consider a RTD composed of two highly doped drain-source GaAs regions, two AlGaAs barriers, and a quantum well ͑see Table II͒.We assume a constant effective mass m = 0.067m 0 , with m 0 as the electron free mass along the whole structure.Transport takes place from emitter to collector in the x direction.The lateral dimensions are small enough to consider electron confinement in y and z directions. 52The energy of an electron in one particular valley is represents the minimum energy of the first subband, whose value is E 1D q = 0.137 eV for L z = 9 nm and L y = 9 nm.The energies of the next lowest subbands are inaccessible to electrons ͑E 1D q = 0.551 eV or E 1D q = 1.239 eV͒.Again, room temperature is assumed.
The practical quantum algorithm for the RTD implies solving numerically N͑t͒ time-dependent single-particle 1D Schrödinger equation ͓Eq.͑25͔͒ for the transport direction x.Due to the confinement in the lateral directions, we assume that the Bohm velocity in y and z directions is negligible. 52ince expression ͑25͒ deals with time-dependent potential profiles, its solution must be computed with a numerical finite-difference scheme method ͑see the numerical algorithm presented in Appendix A of Ref. 54͒.In particular, as discussed in Ref. 33, we assume the zero-order Taylor ap-  proximations G a ͑r ជ a , R ជ a ͓t͔ , t͒ϷG a ͑r ជ a ͓t͔ , R ជ a ͓t͔ , t͒ and J a ͑r ជ a , R ជ a ͓t͔ , t͒ϷJ a ͑r ជ a ͓t͔ , R ជ a ͓t͔ , t͒ in expression ͑25͒.We do also emphasize that the term U a ͑r ជ a , R ជ a ͓t͔ , t͒ that appears in Eq. ͑25͒ contains the full ͑long and short ranges͒ Coulomb interaction with the particular boundary conditions developed in Sec.III.All Schrödinger equations are coupled by N͑t͒ 3D Poisson equations with N͑t͒ different boundary conditions and charge densities ͓expressions ͑17͒-͑20͔͒.The total number of cells, Nx • Ny • Nz, on the order of 1000, and the number of electrons, N͑t͒, about 10-20, implies a computational time on the order of 2-3 days for each bias point, 53 while it takes 10 h within the mean-field approximation.The calculations of the mean-field approximations and our manyelectron approach are identical except in the computation of the potential profile.In the former a unique potential profile is computed, while in the later there is one potential profile for each electron.Finally, let us notice that we consider resonant tunneling, electron confinement, and Coulomb interaction in our quantum solution of the many-particle Hamiltonian, but we do not include the exchange interaction among electrons.The algorithm to include such interaction in our quantum ͑Bohm͒ trajectory proposal is presented in Ref. 33.In Fig. 13 average values for current are presented for different biases with the many-electron and mean-field algorithms.We compute the average current at each bias point using a detailed version of the Ramo-Shockley theorem 54 in surface S4 ͑emitter͒ and surface S1 ͑collector͒.As expected, identical results are obtained from both surfaces showing the numerical accuracy of our simulator.When the RTD remains far from the resonant voltage, the results are quite similar for both methods, but the many-electron approach provides a bit higher current because it avoids the self-reflected electrons in the contact that are found in the mean-field approach, as mentioned previously in Figs.4-6.On the contrary, in the resonant region, the correct consideration of the electronelectron interaction is very relevant because the quantum transport is very sensible to the quantum well electrostatics.Now, the potential profile determines the shape of the quantum well and, therefore, the resonant energies ͓dashed horizontal line in insets ͑a͒ and ͑b͒ of Fig. 13͔ of the electrons.When a "mean-field" electron tries to traverse the "empty" double barrier structure, it "feels" a perturbation in the quantum well due to its own charge implying an increase in the resonant energy and the possibility of being finally reflected by its own charge.In other words, the "mean-field" electron can be Coulomb blockaded by itself.Our many-electron al-gorithm is free from this pathological behavior.This important difference explains the spurious reduction in the current with the mean-field method at resonance.It also explains the movement of the position of the resonant voltage ͑i.e., the voltage at the maximum current͒ as schematically explained in the insets of Fig. 13.The inset of Fig. 13͑a͒ shows how the electron that traverses an empty quantum well feels its own repulsion when the mean-field approximation is used ͑increasing the resonant energy in dotted line͒, while in inset of Fig. 13͑b͒, the electron with a many-electron simulation is free from this pathological effect.
In summary, the relation between the shape of the potential profile and the behavior of the electron can be much more complex in the quantum regime than in its classical counterpart ͑where the spatial derivative of the potential profile directly defines the electron acceleration͒.Figure 13 shows the importance of providing the exact many-particle Coulomb description of electrons in the current-voltage characteristic of a RTD.The time-dependent mean-field approach used here provides spurious effects on the correlations of electrons that are evident even in the dc behavior of RTD simulations.To be fair, let us notice that the time-dependent mean-field approach ͑the same as in Ref. 47͒ does improve the treatment of the Coulomb correlations when compared with the standard Fermi-liquid approaches 13,14 because, in spite of providing a pathological "autointeraction" with itself, it captures the Coulomb correlation between one electron and the others. 47In any case, the many-electron approach developed here, with the exact description of the electron Coulomb interactions, is greatly preferred.

VII. CONCLUSIONS
The prediction of the collective behavior of many electrons is a very active field of research and several theoretical FIG. 13. ͑Color online͒ Average current through surfaces S1 and S4 for the RTD of Fig. 12 as a function of bias using the manyelectron ͑solid symbols͒ and mean-field ͑open symbols͒ algorithms ͑lines are a visual help͒.Nonuniform voltages steps are used to focus on the resonant region.Insets show schematically the effect of an electron crossing an empty well on its own electrostatic potential using ͑a͒ the mean-field or ͑b͒ the many-electron approaches.
approximations have been proposed to improve the treatment of electron-electron correlations.In this work, an exact many-particle Hamiltonian for N͑t͒ electrons inside an open system is developed, without any mean-field or perturbative approximation.The many-particle Hamiltonian ͑6͒ is built from a sum of N͑t͒ electrostatic potentials W k ͑r ជ 1 , ... ,r ជ k , ... ,r ជ N͑t͒ ͒ solutions of N͑t͒ Poisson equation ͓Eq.͑7͔͒ in a 3D volume.We use the Poisson equation to define W k ͑r ជ 1 , ... ,r ជ k , ... ,r ជ N͑t͒ ͒ instead of the Coulomb law because the former is valid for scenarios with ͑or without͒ a spatial-dependent permittivity.In particular, it is shown that the boundary conditions ͓Eq.͑12͔͒ are different for each term W k ͑r ជ 1 , ... ,r ជ k , ... ,r ជ N͑t͒ ͒.It is shown that these particle- dependent boundary conditions ͓Eq.͑12͔͒ of the electrostatic potentials provide the same electron dynamics than the image-charge method applied to electron transport.However, our many-particle approach has the fundamental advantage that it can be directly implemented into 3D realistic ͑classical or quantum͒ electron device simulators, while the imagecharge method is an excellent analytical approach applicable only to very simple systems ͑such as one electron crossing an ideal 41 infinite metallic surface͒.
A classical solution of the many-particle Hamiltonian is presented for a DG-FET.The results are compared with a time-dependent mean-field approach described in Appendix B. Within the mean-field approximation only one potential profile W ¯mean ͑r ជ , t͒ is calculated for all electrons.Then, each electron can be reflected by an artificial alteration of the potential profile due to its own charge.On the contrary, in the many-electron algorithm described here, electrons are not affected by their own charge.The average current and the number of particles are computed for the DG-FET showing that the differences between the mean-field approximation and the exact many-electron approach become important when small geometries ͑that imply stronger electrostatic in-teraction͒ are involved.
A quantum solution of the many-particle Schrödinger equation with the exact many-particle ͑open-system͒ Hamiltonian developed here is presented in terms of the quantum ͑Bohm͒ trajectory algorithm mentioned in Ref. 33.The relevant point of the quantum-trajectory formalism is that Bohm trajectories can be computed without the full knowledge of the many-particle wave function, ⌽͑r ជ 1 , ... ,r ជ N͑t͒ , t͒, but with the knowledge of the single-particle wave function, ⌿ a ͑r ជ a , t͒.It is emphasized that the approach presented in Ref. 33 has similarities with the density-functional theory.In both, the decomposition of many-particle system into a coupled set of single-particle Schrödinger equations is exact, but both need an approximation for the single-particle potentials that appear in their equations ͓i.e., the exchangecorrelation functionals in the latter and the terms G a ͑r ជ a , R ជ a ͓t͔ , t͒ and J a ͑r ជ a , R ជ a ͓t͔ , t͒ in the former͔.We do also emphasize that the electrostatic term U a ͑r ជ a , R ជ a ͓t͔ , t͒ that appears in the time-dependent Schrödinger equation with timedependent potentials, expression ͑25͒, contains the full ͑long and short ranges͒ Coulomb interaction with the particular boundary conditions developed in Sec.III.Numerical results are presented for a RTD and compared with time-dependent mean-field approach developed in Appendix B. The many-electron approach developed here is greatly preferred because it avoids the "self-interaction" found in the timedependent mean-field approach discussed in Appendix B.
Finally, since either the classical or the quantum manyelectron solutions of the Hamiltonian are time-dependent Coulomb-interacting algorithms, apart from the average ͑dc͒ current shown in this work, both many-particle approaches are a really valuable simulation tools to obtain reliable information on the high-frequency and ͑dc and ac͒ noise performances of the state-of-the-art nanoscale devices.Future work will follow this direction.with our algorithm and both have identical boundary conditions, the uniqueness electrostatic theorem guarantees that the potential distributions obtained by both algorithms are identical not only in the boundaries but in any point inside the volume ⍀.Therefore, we obtain the identity W ¯i image ͑r ជ , t͒ = W ¯i͑r ជ , t͒; ∀r ជ ⍀.An explanation in terms of the electric field, rather than the electrostatic potential, follows identical steps.
Finally, as mentioned in the paper, let us remind that the demonstration of the identity is quite simple, but the relevant point to compare both methods is that finding the imaginary charges, ͚ j Ј =N͑t͒+1 MЈ q j Ј / 4͉r ជ l − r ជ j Ј ͉, which fulfill expression ͑A1͒ in 3D realistic scenarios is not at all trivial. 41,42On the other hand, the ability of our many-particle Hamiltonian to be included into 3D realistic devices is explicitly demonstrated in Sec.VI.

APPENDIX B: TIME-DEPENDENT MEAN-FIELD APPROXIMATION FOR THE MANY-PARTICLE HAMILTONIAN
As described in Sec.I, the mean-field approximation provides a single average potential for computing the dynamics of all the electrons.This average potential, which we label here by the suffix "mean" W ¯mean ͑r ជ , t͒, is still capable of pre- serving most of the collective effects of the Coulomb interaction.Here, we compare this approximation with our exact many-particle Hamiltonian.The term W ¯mean ͑r ជ , t͒ is computed by taking into account all charges inside the volume ⍀.However, since one particle cannot "feel" its own charge, in fact, W ¯mean ͑r ជ , t͒ can be interpreted as the electrostatic poten- tial seen for an additional probe charge whose position is r ជ, where the charge density is defined as ¯mean ͑r ជ,t͒ = ͚ j=1 N͑t͒ q j ␦͑r ជ − r ជ j ͓t͔͒,

͑B4͒
Let us notice that the time-dependent mean-field approximation discussed here can be applied to either classical or quantum systems.Both approaches share expressions ͑B2͒-͑B4͒ for the computation of the electrostatic potentials ͑change the classical trajectories by the quantum ones͒.We also want to remark the time-dependence of expression ͑B2͒.This is a common feature for classical ͑semiconductor Monte Carlo 51 ͒ simulations but less frequent for quantum meanfield approaches.Now, we estimate the error of our time-dependent meanfield approximation.First, we show that the mean-field potential can be written in terms of the potentials W ¯i͑r ជ i , t͒ men- tioned in Eq. ͑17͒.In particular, we can write the mean-field potential W ¯mean ͑r ជ , t͒ as Expression ͑B7͒ shows that Error k ͑r ជ , t͒ → ϱ when r ជ → r ជ k ͓t͔.The mean-field approximation implies that the po- tential "felt" by the k particle at r ជ → r ជ k ͓t͔ is its own potential profile.In fact, from a numerical point of view, the use of the mean-field approximation is not so bad.For example, classical simulators uses 3D meshes with cell sizes of few nanometers, DX Ϸ DY Ϸ DZ ӷ 10 nm.Then, the error of the meanfield approximation is smaller than the technical error ͑i.e., mesh error͒ due to the finite size of the cells.The long-range Coulomb interaction is well captured with the mean-field approximation, while this approximation is really bad strategy to capture the short-range Coulomb interaction [28][29][30][31][32] ͑see Fig. 8 in the paper͒.
Finally, let us remark another important point about the mean-field approximation.Looking at final expression ͑B7͒, rewritten here as W k ͑r ជ , t͒ = W ¯mean ͑r ជ , t͒ − V͑r ជ , r ជ k ͓t͔͒, it seems that W k ͑r ជ , t͒ can be computed from a unique mean-field so- lution of the Poisson equation W ¯mean ͑r ជ , t͒ when subtracting the appropriate two-particle potential V͑r ជ , r ជ k ͓t͔͒.In fact, this is the strategy used in several recent works [28][29][30][31][32] to improve the treatments of the short-range Coulomb interaction in electron device Monte Carlo simulators.However, this strategy is not as general as our procedure because these works need an analytical expression for the two-particle Coulomb interaction V͑r ជ , r ជ k ͓t͔͒.The analytical expression of V͑r ជ , r ជ k ͓t͔͒ written in expression ͑3͒ is only valid for scenarios with homogenous permittivity.On the contrary, our procedure with N͑t͒ electrostatic potentials computed from N͑t͒ different Poisson equations in a limited 3D volume ⍀ can be applied inside general scenario with ͑or without͒ spatialdependent permittivity.*Corresponding author: xavier.oriols@uab.es

FIG. 1 .
FIG. 1. ͑Color online͒ Schematic representation of the volume ⍀ = Lx • Ly • Lz and its limiting surface S = ͕S 1 , S 2 , . . .,S 6 ͖.There are N͑t͒ particles inside and M − N͑t͒ outside this volume.The vector r ជ l points an arbitrary position at the boundary surface S l .
FIG.2.͑Color online͒ The electrostatic potential B͑r ជ 1 , . . .,r ជ N͑t͒ , . . .,r ជ M , r ជ l , t͒ measured on a surface S l at position r ជ l and time t ͑due to internal and external electrons͒ by an additional probe charge q M+1 situated on the boundary r ជ l ϵ r ជ M+1 S l .
is divided into Nx • Ny • Nz cells.Each 3D cell has spatial dimensions of DX, DY, and DZ.Thus, the active region of our simulating device has a volume equal to ͑Nx • DX͒ • ͑Ny • DY͒ • ͑Nz • DZ͒ = L x • L y • L z .The boundary conditions of the Poisson equation on the six rectangular surfaces S = ͕S 1 , S 2 , ... ,S 6 ͖ of Fig. 1 are defined using either Dirichlet or Neumann arguments.In general, on the surfaces S 2 , S 3 , S 5 , and S 6 , Neumann boundary conditions are used

FIG. 4 .
FIG. 4. ͑Color online͒ Potential energy profile W ͑r ជ com- with Poisson solver using classical mean-field approximation on the plane X-Y of the active region ⍀ = ͑20 nm͒ 3 at z = 6 nm at 0.4 fs.The solid points are electron positions.

FIG. 10 .
FIG. 10. ͑Color online͒ Average current for the 1D DG-FET of Fig. 9 using the and mean-field algorithms.The open ellipses include results with the same gate voltage indicated on the

TABLE I .
Parameters for the DG-FET depicted in Fig.9.

TABLE II .
Parameters for the RTD depicted in Fig.12.