## Computation of many-particle quantum trajectories with exchange interaction: Application to the simulation of nanoelectronic devices

## A. Alarc´n, S. Yaro, X. Cartoix` and X. Oriols o a

Departament d'Enginyeria Electr`nica, Universitat Aut`noma de Barcelona, 08193, o o Bellaterra, SPAIN

E-mail:

xavier.oriols@uab.es

Abstract. Following Ref. [Oriols X 2007 Phys. Rev. Lett., 98 066803], an algorithm to deal with the exchange interaction in non-separable quantum systems is presented. The algorithm can be applied to fermions or bosons and, by construction, it exactly ensures that any observable is totally independent from the interchange of particles. It is based on the use of conditional Bohmian wave functions which are solutions of singleparticle pseudo-Schr¨dinger equations. o The exchange symmetry is directly defined by demanding symmetry properties of the quantum trajectories in the configuration space with a universal algorithm, rather than through a particular exchange-correlation functional introduced into the single-particle pseudo-Schr¨dinger equation. o It requires the computation of N 2 conditional wave functions to deal with N identical particles. For separable Hamiltonians, the algorithm reduces to the standard Slater determinant for fermions, or permanent for bosons. A numerical test for a two-particle system, where exact solutions for non-separable Hamiltonians are computationally accessible, is presented. The numerical viability of the algorithm for quantum electron transport (in a far-from equilibrium time-dependent open system) is demonstrated by computing the current and fluctuations in a nano-resistor, with exchange and Coulomb interactions among electrons.

## 1. Introduction

A system with N identical particles gives rise to a host of fascinating phenomena. Only those wave functions whose probability density remains unchanged under permutations of particles are a good description of such system. For separable Hamiltonians, these wave functions can be constructed from single-particle wave functions. However, for non-separable Hamiltonians, the computational burden associated with getting the N -particle wave function makes the exact solution inaccessible in most practical situations. This is known as the many-body problem [1].

There has been a constant effort among the scientific community to provide solutions to the many-body problem. The quantum Monte Carlo solutions of the Schr¨dinger o equation provide approximate solutions to exact many-particle Hamiltonians [2, 3]. The Hartree-Fock (HF) algorithm [4, 5] approximates the manyparticle wave function by a single Slater determinant of non-interacting single-particle wave functions. Although it is known that the Hartee-Fock wave function cannot approach the original many-particle wave function, it can provide useful information on the original ground state. Alternatively, density functional theory (DFT) shows that the charge density can be used to compute any observable without the explicit knowledge of the many-particle wave function [6, 7]. Practical computations within DFT make use of the Kohn-Sham theorem [8], which defines a system of N non-interacting single-particle wave functions that are able to provide a system of equations to find the exact charge density of the interacting system. However, the complexity of the many-body system is still present in the so called exchange-correlation functional, which is unknown and needs to be approximated. DFT has had a great success, mostly, in chemistry and material science [9], both, dealing with equilibrium systems. Similar ideas can also be used for non-equilibrium time-dependent scenarios, through the Runge-Gross theorem [10], leading to the time-dependent density functional theory (TDDFT). In contrast to the stationary-state DFT, where accurate exchange functionals exist, approximations to the time-dependent exchange-correlation functionals are still in their infancy. TDDFT has been reformulated in terms of the current density [11, 12] and extended into a stochastic time-dependent current density when the system is interacting with a bath [13].

The common strategy in all many-particle approximations is to obtain the observable result from mathematical entities defined in a real space, R 3 , (single-particle wave functions for HF and charge density for DFT) rather than from the many-particle wave function, whose support is defined in the configuration space R 3 N .

Bohmian mechanics [14, 15, 16, 17] is a consistent explanation of quantum phenomena based on the use of wave functions and trajectories. Apart from its ontological implications, Bohmian mechanics is nowadays used as a mathematical machinery that is able to reproduce the wave function evolution from fluid lines [18, 19, 20, 21, 22]. This is the point of view used in this work to study exchange interaction in many-particle systems [23]. In Bohmian mechanics one can naturally

find a single-particle wave function defined in R 3 , while still capturing many-particle features of the system. Such an entity is named conditional wave function [15], and it is built by substituting all degrees of freedom present in the many-particle wave function, except one, by its corresponding Bohmian trajectories. This substitution produces a single-particle wave function with a complicated time-dependence [18]. Recently, manyparticle Bohmian trajectories associated to the conditional wave function have been investigated by Oriols et al. [18, 24] and the idea of introducing exchange interaction into non-separable systems through conditional wave functions was briefly indicated in the seminal work of Ref. [18].

The purpose of this paper is to present an algorithm to introduce exchange interaction into non-separable systems through conditional wave functions following the idea of Ref. [18]. This paper includes physical discussions, technical details and numerical results, omitted in Ref. [18], that justifies the physical soundness of the proposal. The paper also includes the implementation of the exchange algorithm into a numerical simulator of quantum electron transport, justifying its numerical viability in practical systems. The paper is organized as follows. In section 2 we provide an introduction to many-particle wave functions and Bohmian mechanics. For such introduction, we will use many-particle wave functions for separable Hamiltonians. From a didactic point of view, these simple systems will be useful to discuss how the exchange interaction determines the behavior of the Bohmian trajectories of identical particles. In section 3, we will explain how to compute many-particle Bohmian trajectories for identical particles for non-separable Hamiltonians, without computing the manyparticle wave function. Ensemble results for the kinetic, classical and quantum potential energies will be discussed for systems with and without exchange interaction. Finally, in section 4, we show the numerical viability of the algorithm to include exchange and Coulomb interaction for electron transport simulators. In section 5 we present the conclusions and some additional discussions.

## 2. Many-particle trajectories from many-particle wave functions

In this section, we introduce many particle wave functions and Bohmian mechanics to explain general properties of Bohmian trajectories associated to identical particles. These discussions will be of great utility in the subsequent sections.

## 2.1. Summary of many-particle wave function

For non-relativistic open systems of N -particles, a general expression for a many-particle wave function, Φ ≡ Φ( /vector r , .., /vector 1 r N , t ), with or without exchange interaction, is:

$$\Phi & = C \sum _ { s _ { z j } } \Psi _ { s _ { z 1 }, \dots, s _ { z N } } ( \vec { r } _ { 1 },.., \vec { r } _ { N }, t ) \gamma ( s _ { z 1 },.., s _ { z N } ), & ( 1 ) \\ \intertext { h o r o c o n t e t h e n o i t i o n o f t h a t i o d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d }$$

where /vector r j represents the position of the j -th particle and s zj is the z -component of its spin, which can take the value s zj = /planckover2pi1 / 2 (or ↑ j ) for spin up and s zj = -/planckover2pi1 / 2 (or ↓ j )

for spin down. The normalization constant is C . The sum in (1) is over all possible combinations of spin [25].

In most discussions of this paper (except the numerical results discussed in section 4) we will assume that the quantum system is described by just one of the terms in (1). In particular, we will consider the term where all spins are parallel, e.g., s zj = ↑ j for j = 1 , ...N . In order to simplify our notation, the orbital part of this term will be written as Ψ ≡ Ψ( /vector r , .., /vector 1 r N , t ), without any reference to the spins because their interchange becomes irrelevant. Therefore, the (orbital) wave function is solution of the following many-particle Schr¨dinger equation: o

$$\mu _ { \mu } ^ { \mu } & \quad i \hbar { \frac { } { } } \partial t = \left ( \sum _ { k = 1 } ^ { N } - \frac { \hbar { \ } ^ { 2 } } { 2 m } \nabla _ { k } ^ { 2 } + U ( \vec { r } _ { 1 },.., \vec { r } _ { N }, t ) \right ) \Psi, \\ \text{where $m$ is the free electron mass and $U(\vec{r}_{1},\dots,\vec{r}_{N},t)$ is a non-separable potential. By}$$

where m is the free electron mass and U /vector ( r , . . . , /vector 1 r N , t ) is a non-separable potential. By construction, we know that the solution of (2) satisfies the following continuity equation:

$$\text{there } \vec { I } _ { \infty } \, = \vec { I } _ { \infty } ( \vec { r } _ { k } ) \vec { J } _ { \vec { r } _ { k } } = 0, \\ \text{there } \vec { I } _ { \infty } \, = \vec { I } _ { \infty } ( \vec { r } _ { k } ) \quad \vec { r } _ { \infty } \, t ) \, \text{is the environment values of the current non-habilitiv} \, \text{density}$$

where /vector J /vector r k ≡ /vector J /vector r k ( /vector r , ..., /vector 1 r N , t ) is the expectation values of the current probability density [26] and | Ψ | 2 the presence probability density. This last result will be relevant in section 2.2 when presenting Bohmian trajectories.

Two particles are said to be identical if there are no experiments that can detect differences between them. This restriction on observable results can be satisfied by imposing the following property into the wave function Ψ of identical particles:

$$\Psi (., \vec { r } _ { j },., \vec { r } _ { h },., t ) = e ^ { i \gamma } \Psi (., \vec { r } _ { h },., \vec { r } _ { j },., t ),$$

for any j and h indices. We consider γ = 0 ( mod 2 π ) for bosons (symmetry) and γ = π ( mod 2 π ) for fermions (antisymmetry).

We say that the system has exchange interaction when the wave function satisfies (4). For physical systems of identical particles, the many-particle potential in (2) remains invariant under the permutation of two positions, i.e. U .,/vector ( r , ..., /vector j r h , ., t ) = U .,/vector ( r h , ..., /vector r , ., t j ) for any j and h , and the symmetry or antisymmetry property of the wave function in (4) for time t holds for all instants. Next, as a simple example of the difference between systems with and without exchange interaction, we discuss on the total energy, which will be useful later in section 4.

2.1.1. Example: The effect of exchange interaction on total energy We consider a system of N particles in free space. For simplicity, we consider 1D particles where its position is defined in R . Then, the many-particle wave function Ψ( x , ..., x 1 N , 0) at t = 0 can be constructed from the following single-particle Gaussian wave packets:

$$\psi _ { j } ( x _ { j }, 0 ) = \frac { \exp \left ( i k _ { o j } x _ { j } \right ) } { \left ( \pi \sigma _ { x _ { j } } ^ { 2 } \right ) ^ { 1 / 4 } } \exp \left ( - \frac { \left ( x _ { j } - x _ { o j } \right ) ^ { 2 } } { 2 \sigma _ { x _ { j } } ^ { 2 } } \right ),$$

where σ xj is the spatial dispersion, x oj the central position, E oj = ( /planckover2pi1 k oj ) 2 / (2 m ) the central energy of each wave packets and k oj the central wave vector.

In particular, the N -particle wave function Ψ( x , ..., x 1 N , 0) with exchange interaction can be defined from:

$$\Psi ( x _ { 1 },.., x _ { N }, 0 ) = C \, \sum _ { n = 1 } ^ { N! } \prod _ { j = 1 } ^ { N } \psi _ { j } ( x _ { p ( n ) _ { j } }, 0 ) \, \text{sign} ( \vec { p } _ { n } ), \\ \text{where the sum is over all $N$}! \, \text{permutations $\vec{p}_{n}$} \ = \ \{ p ( n ) _ { 1 } \dots p ( n ) _ { N } \} \text{ and } C \text{ is a}$$

where the sum is over all N ! permutations /vector p n = { p n ( ) 1 , ..., p ( n ) N } and C is a normalization constant. For fermions, the sign( /vector p n ) = ± 1 means the sign of the permutations, i.e. (6) is the Slater determinant. Alternatively, we will consider sign( /vector p n ) = 1 for bosons, meaning that (6) has to be interpreted as the permanent.

On the other hand, the wave function for particles without exchange interaction can be written as:

$$\Phi ( x _ { 1 }, \dots, x _ { N }, 0 ) = \prod _ { j = 1 } ^ { N } \psi _ { j } ( x _ { j }, 0 ), \\ \text{which, by construction, is already well normalized to unity. The ensemble value of the}$$

which, by construction, is already well normalized to unity. The ensemble value of the kinetic energy of the j -th particle belonging to a system of particles without exchange interaction is computed as:

$$\langle T _ { j } \rangle & = \int \dots \int \Psi ^ { * } \hat { T } _ { j } \Psi \ d x _ { 1 } \dots d x _ { N }, \\ \theta \, \text{line-interior over operator is } \hat { T }. \, - \frac { h ^ { 2 } } { } \frac { \partial ^ { 2 } } { } \text{    Horator                                                                                                                                                                                                        
                                                                                                                                                                                                        } \text{                                                                                                                                                                                                        \text{                                                                                                                                                                                                    
   \text{                                                                                                                                                                                                  
     \text{                                                                                                                                                                                                
       \text{$$

The kinetic energy operator is ˆ T j = -/planckover2pi1 2 2 m∂x ∂ 2 2 j . Hereafter, unless specified, the spatial integrals are assumed to extend over the whole configuration space. The same expression (8) can be used for identical particles defined from the wave function in (6). Then, one can easily realize that 〈 T j 〉 = 〈 T h 〉 for any j and h indexes. As expected, one cannot discern between identical particles from the measurement of their kinetic energies.

We compute the behavior of the total kinetic energy 〈 T 〉 = 〈 T 1 〉 + 〈 T 2 〉 + 〈 T 3 〉 for three electrons (with parallel spins) with and without exchange interaction, as a function of the distance among the wave packets in the configuration space (see inset in figure 1). We define the normalized phase-space distance among the central positions and central wave vectors of two wave packets as [27]:

$$d ( 1, j ) ^ { 2 } & = \frac { ( k _ { o 1 } - k _ { o j } ) ^ { 2 } } { 2 \sigma _ { k } ^ { 2 } } + \frac { ( x _ { o 1 } - x _ { o j } ) ^ { 2 } } { 2 \sigma _ { x } ^ { 2 } } ; \ j = 2, 3, \\. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &.. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &. &.$$

where σ kj = 1 /σ xj is the wave vector dispersion. In figure 1, we plot, in a square (black) line, the mean value of the total kinetic energy of three electron (fermions) with exchange interaction, whose wave function is defined from (6). The result is repeated for different values of the distance d = d (1 , 2) = d (1 , 3) with the condition x o 1 -x o 2 = x o 3 -x o 1 and k o 2 -k o 1 = k o 3 -k o 1 seen in the inset of figure 1. Identically, we plot in up triangle (blue) line the total kinetic energy computed for three particles without exchange interaction, whose wave function is defined from (7). For large d , the values of the kinetic energy of the three electrons with and without exchange interaction are identical. For such large values of d , all electrons are placed far away from each other in the phase-space and

the exchange interaction has no effect. However, this is not true for small values of d . Then, the difference between the kinetic energy of electrons with or without exchange interaction increases as we place the three electrons closer inside the phase-space.

Figure 1: (Color online) Ensemble value of the kinetic energy for a 3-particle system with (square solid black line) and without (up triangle solid blue line) exchange interaction as function of their normalized phase-space distance d . The inset shows the positions of the central position and central wave vector of each wave packet in the phase-space, which are used to define the distance d among them.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0388v1-figure-efd36846-98e9-4ff9-a13b-78bfc3d59f3e.jpg)

The result plotted in figure 1 is just (the wave packet version of) the celebrated Pauli exclusion principle: identical fermions cannot be in the same quantum state. The discussion has been done with three particles, instead of two, because in Appendix B we generalize the present example to three electrons with exchange interaction and different spins orientations.

## 2.2. Summary of many-particle trajectories

In Bohmian mechanics [14, 15, 16, 17], each particle of the system is represented by a trajectory guided by a wave. The wave is the many-particle wave function discussed above, Ψ( /vector r , ..., /vector 1 r N , t ), with all its computational difficulties. Such wave function satisfies the continuity equation, written in (3), that relates current and probability presence densities. From such continuity equation, one can easily define a (Bohmian) velocity /vector v j ( /vector r , ..., /vector 1 r N , t ) at each position of the configuration space as:

$$\psi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi\Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phis \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Psi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi
\Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Psilon \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \ Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phie \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi(\Phi)\Phi)\Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phip \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \psi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phis\Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \phie \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi\Psi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phis
\Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi\Phis\Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phis\Phis\Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Ph$$

The (Bohmian) trajectory of the j -th particle, /vector r l j [ t ], in real space can be defined by time-integrating (10) as:

$$\overset { \dots } { \dots } \overset { \dots } { \dots } \circ \overset { \dots } { \dots } \circ \overset { \dots } { \dots } \overset { \dots } { \dots } \overset { \dots } { \dots } \overset { \dots } { \dots } \overset { \dots } { \dots } \overset { \dots } { \dots } \overset { \dots } { \dots } \overset { \dots } { \dots } \overset { \dots } { \dots } \overset { \dots } { \dots } \overset { \dots } { \dots } \overset { \dots } { \dots } \overset { \dots } { \dots } \overset { \dots } \dots \\ \vec { r } _ { j } ^ { \prime } l t ] & = \vec { r } _ { j } ^ { \prime } [ 0 ] + \int _ { 0 } ^ { t } \vec { v } _ { j } ^ { \prime } ( \vec { r } _ { 1 } ^ { \dots } [ t ^ { \prime } ], \dots, \vec { r } _ { N } ^ { \dots } [ t ^ { \prime } ], t ^ { \prime } ) d t ^ { \prime }.$$

Obviously, one has to select the initial position /vector r l j [0] to perfectly specify the trajectory. The super-index l = 1 , ..., M on the trajectory accounts for the M →∞ different initial positions that can be selected. We refer to /vector r l j [ t ] as the Bohmian trajectory in R 3 , while we will refer to { /vector r l 1 [ t , .., /vector ] r l N [ t ] } as a many-particle (or N -particle) Bohmian trajectory in R 3 N .

The relevant property of these Bohmian trajectories that makes them meaningful for quantum computations is the fact that, by construction, a proper ensemble of them (with different initial positions) does exactly reproduce the time-evolution of the manyparticle wave function, at any time. A proper ensemble means that the initial positions, { /vector r l 1 [0] ..../vector r l N [0] } , are selected according to the probability distribution | Ψ( /vector r , ..., /vector 1 r N , 0) | 2 . This last condition is called 'quantum equilibrium hypothesis' [14, 28].

We can now deduce an important property of these trajectories that will be very relevant later. Since Ψ( /vector r , ..., /vector 1 r N , t ) is a single-valued wave function, the Bohmian velocity computed from (10) in each point of the configuration space is unique. This means that if two trajectories coincide at some point of the configuration space, then, they will coincide forever (because their velocities become identical). This well-known result can be summarized in a simple sentence: two many-particle Bohmian trajectories (with different initial positions) do not cross in the configuration space, either for bosons, fermions or non-identical particles [29].

Equivalently, the presentation of such trajectories can be done by introducing the polar form of the many-particle wave function ψ /vector ( r , .., /vector 1 r N , t ) = R /vector ( r , .., /vector 1 r N , t ) e iS /vector ( r 1 ,..,/vector r N ,t ) / /planckover2pi1 into (2). The modulus R ≡ R /vector ( r , .., /vector 1 r N , t ) and the phase S ≡ S /vector ( r , .., /vector 1 r N , t ) are real functions. Then, one obtains again, from the imaginary part of (2), the continuity equation defined in (3) in polar form:

$$\cdot \quad \frac { \partial R ^ { 2 } } { \partial t } + \sum _ { j = 1 } ^ { N } \vec { \nabla } _ { \vec { r } _ { j } } \left ( R ^ { 2 } \frac { \vec { \nabla } _ { \vec { r } _ { j } } S } { m } \right ) = 0,$$

where we recognize the velocity of the j -th particle as:

$$\vec { v } _ { j } ( \vec { r } _ { 1 },.., \vec { r } _ { N }, t ) & = \frac { \vec { \nabla } _ { \vec { r } _ { j } } S ( \vec { r } _ { 1 },.., \vec { r } _ { N }, t ) } { m }. \quad \ \. \quad \. \quad \. \quad \. \quad \. \quad \. \quad \. \quad \. \quad \. \quad \.$$

By construction [15], the velocity definition in (13) is identical to that in (10). On the other hand, the real part of the Schr¨dinger equation leads to a many-particle version o of the quantum Hamilton-Jacobi equation:

$$\bar { \ e r e s t e r e } \frac { \partial S } { \partial t } + U + \sum _ { j = 1 } ^ { N } \left ( K _ { j } + Q _ { j } \right ) = 0, \\ \intertext { r e r e $ U \equiv U ( \vec { r } _ { 1 }, \dots \vec { r } _ { s s. t } ) $ i s $ the $ \text{notential in (2)} $ and we have defined the (local)} $ \text{Bohmian}$$

where U ≡ U /vector ( r , .., /vector 1 r N , t ) is the potential in (2) and we have defined the (local) Bohmian kinetic energy as:

$$K _ { j } \equiv K _ { j } ( \vec { r } _ { 1 },.., \vec { r } _ { N }, t ) = \frac { 1 } { 2 } m \vec { v } _ { j } ( \vec { r } _ { 1 },.., \vec { r } _ { N }, t ) ^ { 2 },$$

and the (local) quantum potential energy:

$$Q _ { j } \equiv Q _ { j } ( \vec { r } _ { 1 },.., \vec { r } _ { N }, t ) = - \frac { \hbar { ^ { 2 } } } { 2 m } \frac { \vec { \nabla } _ { \vec { r } _ { j } } ^ { 2 } R ( \vec { r } _ { 1 },.., \vec { r } _ { N }, t ) } { R ( \vec { r } _ { 1 },.., \vec { r } _ { N }, t ) }.$$

When dealing with Bohmian trajectories, the ensemble kinetic energy defined in (8) is divided into two parts, 〈 ˆ T j 〉 = 〈 ˆ K j 〉 + 〈 ˆ Q j 〉 . The first part:

related to the local (Bohmian) kinetic energy K j ≡ K x ,...x j ( 1 N , t ), and the second part:

$$\left \langle \hat { K } _ { j } \right \rangle = \int \dots \int R ^ { 2 } \, K _ { j } \, d x _ { 1 } \dots d x _ { N }, \\ \text{lated to the local (Bohmian) kinetic energy } K _ { j } \equiv K _ { j } ( x _ { 1 }, \dots x _ { N }, t ), \, \text{and the second part} \colon$$

$$\left \langle \hat { Q } _ { j } \right \rangle = \int \dots \int R ^ { 2 } \, Q _ { j } \, d x _ { 1 } \dots d x _ { N }, \\ \mathbb { m } \text{ the quantum potential energy } Q _ { j } \equiv Q _ { j } ( x _ { 1 }, \dots, x _ { N }, t ).$$

to the quantum potential energy Q j ≡ Q x ,..., x j ( 1 N , t ).

2.2.1. Properties of many-particle Bohmian trajectories with exchange interaction Now, we can list a series of important properties for those ensembles of Bohmian trajectories that represents identical particles, i.e., when exchange interaction is present. In order to simplify the notation, we define /vector X = { /vector r , .., /vector 1 r N } . Identically, we define the N -particle Bohmian trajectory at time t = 0 as /vector X l [0] = { /vector r l 1 [0] , .., /vector r l N [0] } . Another set of initial conditions will be refereed as /vector X f [0] = { ., /vector r f h [0] , ., /vector r f j [0] , . } when it contains the same initials positions as /vector X l [0], but the two initial positions, /vector r l j [0] and /vector r l h [0], are interchanged. Because of (4), the modulus of the many-particle wave function satisfies:

$$R ( \vec { X } ^ { l } [ 0 ], 0 ) = R ( \vec { X } ^ { f } [ 0 ], 0 ),$$

for any such type of two set of initials conditions l and f . Identically, the phase satisfies:

$$S ( \vec { X } ^ { l } [ 0 ], 0 ) = \gamma + S ( \vec { X } ^ { f } [ 0 ], 0 ),$$

where γ = 0 ( mod 2 π ) for bosons (symmetry) and γ = π ( mod 2 π ) for fermions (antisymmetry). As discussed for the wave function, the requirements in (19) and (20) are satisfied at any time t . The property of (20) togther with the definition of the velocity in (13) implies:

$$\vec { v } _ { j } ( \vec { X } ^ { l } [ t ], t ) = \vec { v } _ { h } ( \vec { X } ^ { f } [ t ], t ).$$

This condition on the Bohmian velocities, which is valid for either bosons or fermions, has two relevant consequences. First, let us compare the two sets of many-particle trajectory with different initial positions mentioned above: the l -set and the f -set. Their difference are only /vector r l j [0] = /vector r f h [0] and /vector r l h [0] = /vector r f j [0]. Then, we realize from (21) that all Bohmian trajectories with identical initial conditions will be equal independently of the initial conditions, except the two trajectories which have their initial positions interchanged. For these trajectories, we get /vector r l j [ t ] = /vector r f h [ t ] and /vector r l h [ t ] = /vector r f j [ t ].

The second consequence of (21) is valid for those many-particle trajectories that have, at least, two equal components, i.e. /vector r l j [0] = /vector r l h [0] ≡ /vector a . Because of this coincidence,

we have /vector X l [0] = { ., /vector , a .,/vector , . a } and also /vector X f [0] = { ., /vector , a .,/vector , . a } which in fact are the same. Then, the condition /vector v j ( /vector X t , t l [ ] ) = /vector v h ( /vector X t , t f [ ] ) can be written as:

$$\vec { v } _ { j } ( \vec { X } ^ { l } [ t ], 0 ) = \vec { v } _ { h } ( \vec { X } ^ { l } [ t ], 0 ) \equiv v _ { a },$$

because /vector X l [0] = /vector X f [0]. Then, the trajectory /vector r l j [ t ] at the subsequent time /vector r l j [0 + dt ] = /vector a + /vector v dt a is identical to the other trajectory /vector r l h [0 + dt ] = /vector a + /vector v dt a . This result, means /vector r l j [ t ] = /vector r l h [ t ] at any time.

Because of the previous property and the non-crossing property of Bohmian trajectories discussed before [29], we have an important corollary. We define 'diagonal' many-particle trajectories as those trajectories where at least two components, /vector r l j [ t ] = /vector r l h [ t ], are identical (the rest of components can be different). Since other Bohmian trajectories cannot cross such 'diagonal' trajectories, all Bohmian trajectories are restricted to remain in subspaces of the configuration space. According to Ref. [30], Bohmian mechanics for identical particles can be described in a 'reduced' space R 3 N /S N , with S N the permutation space of N -particles.

Finally, we want to mention that in Bohmian computations, even with the symmetrization postulate, trajectories of particles are obviously distinguishable. One labels the trajectory of particle 1 as /vector r l 1 [ t ] and that of particle 2 as /vector r l 2 [ t ]. We have shown that, by construction, the Bohmian trajectories have special symmetry requirements. Then, all results for particle 1 computed from an ensemble of these trajectories will be identical to those computed for particle 2. In simple words, for a system of identical particles, Bohmian trajectories are distinguishable, while observable results associated to different particles become indistinguishable.

2.2.2. Example: The effect of exchange interaction on Bohmian trajectories Let us discuss, with some numerical examples, the previous properties of Bohmian trajectories of identical particles. In all the numerical examples of this subsection, we consider two free particles propagating, each one, in 1D physical space. The single-particle wave packets that will be used to construct the many-particle wave function at the initial time t = 0 are defined from (5).

First, we consider two electrons with a wave function Ψ( x , x 1 2 , t ) computed from (7), without any symmetry. See the initial modulus of the 2-particle wave function in figure 2. In particular, we consider E o 1 = 0 12 eV, . x o 1 = +50 nm and σ x 1 = 25 nm for the first wave packet, and E o 2 = 0 08 eV, . x o 2 = -50 nm and σ x 2 = 25 nm for the second. In order to see the spatial interaction of the two particles, the momentum of the first particle is negative and that of the second positive. We consider a free electron mass for both electrons. Once we know Ψ( x , x 1 2 , t ), we compute the (two-particle) Bohmian trajectory from (11) with different initial positions. As seen in figure 3, they correspond to roughly parallel lines.

In figure 4, we see that the quantum potential of either the first or the second particles are nearly zero. The total quantum potential, as the sum of the two particles, is plotted in dashed square (red) line. The ensemble (Bohmian) kinetic energy remains

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0388v1-figure-6d981886-0c4b-4f95-bef3-d0f941ff3cce.jpg)

Figure 2: (Color online) Modulus of the wave function for two particles without exchange interaction in the 2D configuration space at t = 0 fs.

Figure 3: (Color online) Two-particles Bohmian trajectories with different initials conditions for particles without exchange interaction in a free space.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0388v1-figure-41461886-6a4c-4620-97f2-497b79264173.jpg)

Figure 4: (Color online) Time evolution of the total and individual (ensemble average) energies of two-electron system without exchange interaction in free space.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0388v1-figure-10855e8e-db2a-4b9b-8c07-f82cd566bce7.jpg)

equal to its initial value, 0 12 eV for the first wave packet in solid circle (blue) line and . 0 08 eV for the second in solid plus (blue) line. . The total energy in solid (green) line remains constant and equal to 0 2 eV, i.e., the sum of kinetic energies. . These simple Bohmian trajectories for non-identical particles move roughly like classical particles.

Next, we consider a wave function Ψ( x , x , t 1 2 ) of two identical electrons computed from the Slater determinant of (6) for N = 2. We use the same two initial gaussian wave packets discussed above. In figure 5(a), we plot the (symmetric) modulus of the manyparticle wave functions at t = 267 8 fs. . In particular, we get Ψ( a, a, t ) = 0 at any point { a, a } of the diagonal. In figure 6(a), we plot a set of Bohmian trajectories. The initials positions { x l 1 [0] , x l 2 [0] } are selected symmetrically with respect to the 'diagonal'. First, we observe that x l 1 [ t ] = x f 2 [ t ] and x l 2 [ t ] = x f 1 [ t ] when x l 1 [0] = x f 2 [0] and x l 2 [0] = x f 1 [0]. As discussed in section 2.2.1, the Bohmian trajectories corresponding to interchanged initial positions become symmetrical with respect to the diagonal points of the configuration space. Second, we observe that the Bohmian trajectories do not cross the diagonal.

In figure 7(a), we plot the energies of this two-particle fermion system. The total energy of the identical particles is equal to that of the particles without exchange interaction discussed in figure 4. The reason, as explained in section 2.1.1, is because the momentum of the wave packets are very different. One momentum is positive and the other negative and no Pauli effect is observed in the energy. However, since Bohmian trajectories are 'reflected' at the diagonal, their (bohmian) velocity becomes zero at that time. Then, the ensemble average of K x ,x ,t j ( 1 2 ) in (17) is almost zero, while the ensemble average of Q x ,x , t j ( 1 2 ) grows to keep the total energy constant. The same result can be argued by noting that the quantum potential in (16) depends on the curvature of the modulus, which becomes large at that points. In addition, in contrast to the two particles without exchange interaction discussed in figure 4, the (Bohmian kinetic plus quantum) energies of the first particle are identical to those of the second particle. The observable results of the energy of the individual particles are indistinguishable, while we can perfectly distinguish the trajectories in figure 5(a).

/negationslash

Finally, in figure 5(b), figure 6(b) and figure 7(b), we plot the same result as in the previous figures but considering two identical bosons. We use exactly the same parameters for the wave packets discussed in the previous figures. The only difference is that the initial wave function is computed from (6) when the sign( /vector p n ) is substituted by 1. Let us notice again the symmetric property of the modulus of the wave function in the configuration space. Although we have Ψ( a, a, t ) = 0 at the diagonal points { a, a } , we see in figure 6(b), that Bohmian trajectories do not cross that diagonal. This is an expected result because our discussions on the properties of Bohmian trajectories in section 2.2.1 do not depend on the bosonic or fermionic nature of particles. There is a Bohmian trajectory located along the diagonal points of the configuration space (not plotted) that does not allow to be crossed by other trajectories. The initials positions { x l 1 [0] , x l 2 [0] } are selected symmetrically with respect to the 'diagonal' and identical to the ones used for the fermions. Again, trajectories are symmetric under the exchange of initial positions. In figure 7(b), we plot the energies of the two-particle bosonic system.

The numerical results of the (Bohmian) kinetic energy for bosons (0.019 eV) are slightly lower than fermions (0.022 eV) when the wave packet is close to the diagonal of the configuration space. The reason is because there are more bosonic Bohmian trajectories that arrive closer to the diagonal in figure 5(b) and figure 6(b) than the fermionic ones in figure 5(a) and figure 6(a).

## 3. Many-particle trajectories without many-particle wave functions

As commented in the introduction, many attempts have been developed in the literature to provide accurate solutions to the many-body problem. Here, we briefly review one of this approximations presented by one of the authors in Ref. [18]. Then, we explain how the exchange interaction can be included in the mentioned approximation.

## 3.1. The conditional wave function

The main idea behind the many-body approximation mentioned in Ref. [18] is the fact that the computation of the Bohmian velocity for the /vector r a [ t ] trajectory from (10) only requires the spatial derivatives of Ψ( /vector r , .., /vector 1 r , .., /vector a r N , t ) on the /vector r a directions, and not on the rest of degrees of freedom. Thus, in principle, the trajectory /vector r a [ t ] can be equivalently computed from the many-body wave function Ψ( /vector r , ..., /vector 1 r N , t ) or from the following conditional wave function:

$$\Psi _ { a } ( \vec { r } _ { a }, t ) = \Psi ( \vec { r } _ { a }, \vec { X } _ { a } [ t ], t ),$$

where /vector X t a [ ] = { /vector r 1 [ t , /vector ] r a -1 [ t , /vector ] r a +1 [ t , /vector ] r N [ t ] } is a vector that contains all Bohmian trajectories except /vector r a [ t ]. We also use /vector X a = { /vector r , .../vector 1 r a -1 , /vector r a +1 , .., /vector r N } when referring to all the degrees of freedom except /vector r a . When not relevant, we avoid the superindex l in the Bohmian trajectory that specifies the initial positions of the trajectory. Certainly, the conditional wave function in (23) is defined in a much smaller configuration space, R 3 , than the many-body wave function. Thus, in principle, the conditional wave function needs much less computational effort than the explicit many-particle wave function. Following Ref. [18], the single-particle wave function Ψ ( a /vector r , t a ), that we will use to compute /vector r a [ t ], can be obtained as a solution of the single-particle Schr¨dinger equation: o

$$\text{unpow } r _ { a [ \nu ] }, \, \text{can use a non-zero} \, \text{on any} \, \text{equation}. \\ i \hbar { \frac { \partial \Psi _ { a } ( \vec { r } _ { a }, t ) } { \partial t } } & = \left ( - \frac { \hbar { ^ { 2 } } } { 2 m } \nabla _ { \vec { r } _ { a } } ^ { 2 } + U _ { a } ( \vec { r } _ { a }, \vec { X } _ { a } [ t ], t ) + G _ { a } ( \vec { r } _ { a }, \vec { X } _ { a } [ t ], t ) \\ + i J _ { a } ( \vec { r } _ { a }, \vec { X } _ { a } [ t ], t ) \right ) \Psi _ { a } ( \vec { r } _ { a }, t ). \\ \text{The exact definition of the terms } U _ { a } ( \vec { r } _ { a }, \vec { X } _ { a }, t ), \, G _ { a } ( \vec { r } _ { a }, \vec { X } _ { a }, t ) \text{ and } J _ { a } ( \vec { r } _ { a }, \vec { X } _ { a }, t ) \text{ can be}$$

The exact definition of the terms U a ( /vector r , X , t a /vector a ), G a ( /vector r , X a /vector a , t ) and J a ( /vector r , X a /vector a , t ) can be found in Ref. [18].

In brief, we have been able to decompose an irresolvable N -particle Schr¨dinger o equation into a set of N -single-particle Schr¨dinger o equation with time-dependent potentials [18]. At this point we realize that the extraordinary numerical simplification comes at the prize that there are terms in (24) which are unknown and need pertinent

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0388v1-figure-ef38d358-ac8e-4e5e-af2d-3fe484660800.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0388v1-figure-91239600-d3fe-4baf-bf4d-c06090361731.jpg)

Figure 5: (Color online) Modulus of the wave function for two identical particles in the 2D configuration space. (a) two fermions at t = 267 8 fs and (b) two bosons at . t = 178 5 fs. .

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0388v1-figure-50073589-d36b-4cd6-8cac-deb430c39857.jpg)

Figure 6: (Color online) Two-particle Bohmian trajectories with different initial conditions in free space. (a) Fermions and (b) bosons. The inset is a zoom of the diagonal non-crossing properties of Bohmian trajectories.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0388v1-figure-c3d0b03a-859b-498c-853a-f1e9626190dd.jpg)

Figure 7: (Color online) Time evolution of the total and individual (ensemble average) energies of the two-particle system in a free space. (a) Fermions and (b) bosons.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0388v1-figure-207b8ac2-1c48-43e6-b4d9-e5b2e49442a9.jpg)

approximations, G a ( /vector r , X a /vector a [ t , t ] ) and J a ( /vector r , X a /vector a [ t , t ] ). This is a similar situation to that in DFT discussed in the introduction.

3.1.1. Test for non-separable harmonic potentials without exchange interaction Next, in order to clarify the use of conditional (Bohmian) trajectories discussed above, we applied it to a simple system of two electrons without exchange interaction under a non-separable Hamiltonian. We consider two 1D particles so that the configuration space is R 2 . We use the non-separable potential energy:

$$U ( x _ { 1 }, x _ { 2 } ) = c ( x _ { 1 } - x _ { 2 } ) ^ { 2 },$$

where the factor c will allow us to modify arbitrarily the strength of the non-separable interaction. In particular, we will use c = 10 12 eV/m 2 . The many-body wave function Ψ( x , x , t 1 2 ) can be solved exactly from (2) with N = 2. Once the exact 2D wave function Ψ( x , x 1 2 , t ) is known, we can compute the exact 2D Bohmian trajectories straightforwardly from (10).

In figure 8, we have plotted the ensemble results of the (Bohmian) kinetic energy, (17), the quantum potential energy, (18), for the two electrons. We compute the results directly from the 2D exact wave function solution of (2). We emphasize that there is an interchange of kinetic energies between the first and second particles (see their kinetic energy in the first and second oscillations). This effect clearly manifests that the Hamiltonian of that quantum system is non-separable.

Alternatively, we can compute the trajectories used to compute figure 8 without knowing the many-particle wave function, but computing the conditional wave function Ψ ( a x , t a ) solution of (24) with the proper approximation for terms G a and J a . Here, we consider a zero order Taylor expansion around x a [ t ] for the unknown potentials terms G a and J a . In other words, we consider them as purely time-dependent potential terms, G a ( x , x a b [ t , t ] ) ≈ G ′′ a ( x a [ t , t ] ) and J a ( x , x a b [ t , t ] ) ≈ J ′′ a ( x a [ t , t ] ). This is the simplest approximation. Then, we know that the (complex) purely time-dependent terms G ′′ a ( x a [ t , X ] /vector a [ t , t ] ) and J ′′ a ( x a [ t , X ] /vector a [ t , t ] ) in the Hamiltonian of (24) only introduce a (complex) purely time-dependent phase. Then, we can write Ψ ( a x , t a ) as:

$$\Psi _ { a } ( x _ { a }, t ) = \tilde { \psi } _ { a } ( x _ { a }, t ) \exp ( z _ { a } ( t ) ),$$

where the term z a ( ) t is the (complex) purely time-dependent term that has no effect on the Bohmian trajectory x a [ t ], because this phase has no spatial dependence. Then, under the previous approximation, (24) can be simplified into the following equation for the computation of ˜ ( ψ a x , t a ):

$$i \hbar { \frac { \partial \tilde { \psi } _ { a } } { ( x _ { a }, t ) } } & = \left ( - \frac { \hbar { ^ { 2 } } } { 2 m } \frac { \partial ^ { 2 } } { \partial x _ { a } ^ { 2 } } + U _ { a } ( x _ { a }, x _ { b } [ t ] ) \right ) \tilde { \psi } _ { a } ( x _ { a }, t ), \\ \intertext { a r a b a t i o n i t i o n a r e i o n } \ e n \ h a \ I I I ( r _ { a }, r _ { b } [ t ] ) \, = \, \ e ( r _ { a }, \ e _ { b } [ t ] \vee 2 \ f o r \ a \ - \ 1 \ e n d }$$

Here, the potential energies can be U 1 ( x , x 1 2 [ t ]) = c x ( 1 -x 2 [ t ]) 2 for a = 1 and U 2 ( x , x 2 1 [ t ]) = c x ( 1 [ t ] -x 2 ) 2 for a = 2. The initial wave functions are ψ 1 ( x , 1 0) and ψ 2 ( x , 2 0) defined, both, form (5). In particular, we consider E o 1 = 0 06 eV, . x o 1 = 50 nm and σ x 1 = 25 nm for the first wave packet, and E o 2 = 0 04 eV, . x o 2 = -50 nm

and σ x 2 = 25 nm for the second. In general, we need N -conditional wave functions to compute one N -particle Bohmian trajectory. If we change the initials positions, we need new N -conditional wave functions.

Figure 8: (Color online) Time evolution of the total and individual (ensemble averaged) energies of two-electron system without exchange interaction under a non-separable potential.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0388v1-figure-6db39f7c-a91a-440d-ba65-278a6af377ae.jpg)

In figure 9, we have plotted the same information than in figure 8 with our singleparticle 1D approximation algorithm explained in section 3.1. For this particular scenario, our simplest approximation for the unknown terms works perfectly and the agreement between 2D exact results and our 1D approximation is excellent. In general, potentials with small spatial variations are better adapted to the simplest 1 D approximation of the term G a and J a used in this work. We emphasize that the kinetic energy of the first and second particles are clearly distinguishable. We have compute the ensemble energies in order to justify that the algorithm is accurate not only for an arbitrarily selected set of Bohmian trajectory, but for most of them. In particular, the ensemble results are computed from 160000 two-particle Bohmian trajectories.

Figure 9: (Color online) Time evolution of individual (ensemble averaged) Bohmian kinetic energies of identical two-electron system without exchange interaction under non-separable potential computed from 2 D exact and 1 D approximate solutions.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0388v1-figure-0bb84ccc-44dc-4e5c-9d75-5ed09dc0a21e.jpg)

3.2. Algorithm to include exchange interaction in many-particle Bohmian trajectories

Since (24) is valid for system with or without exchange interaction, one could look for approximations to the terms G a and J a different from the simplest one mentioned above in order to incorporate the exchange interaction directly into (24). However, from a computational point of view, such approximations seems quite difficult to implement. For example, in a system of fermions, we have seen that the wave function Ψ ( a /vector r , t a ) becomes zero everywhere /vector r a is equal to the position of another trajectory. At these positions, because of their dependence on the inverse of the modulus, we would obtain G a → ±∞ and J a → ±∞ . These infinities are difficult to treat numerically. In this subsection we present a different strategy that will be able to capture the exchange interaction avoiding the previous difficulties. The algorithm can be explained in four steps:

1.The first step is developing an expression for Ψ( /vector r , .../vector 1 r N , t ) as a sum of wave functions. Each one of these wave function without symmetry. For example, let us define Ψ ns ( /vector r , .../vector 1 r N , t ) as a many-particle wave function without any (bosonic or fermionic) symmetry. Then, we can construct a global wave function with exchange interaction using a sum of the term Ψ ns ( /vector r , .../vector 1 r N , t ) with all possible permutations of the positions:

$$\Psi & = C \sum _ { n = 1 } ^ { N! } \Psi _ { n s } ( \vec { r } _ { p ( n ) _ { 1 } },.., \vec { r } _ { p ( n ) _ { N } }, t ) \, \text{sign} \left ( \vec { p } ( n ) \right ), \\ \text{et us emphasize that each term } \Psi _ { n s } ( \vec { r } _ { n ( n ) },., \vec { r } _ { n ( n ) _ { n } },.., \vec { r } _ { n ( n ) _ { x } }, t ) \text{ is also a solution of}$$

Let us emphasize that each term Ψ ns ( /vector r p n ( ) 1 , /vector r p n ( ) 2 , ..., /vector r p n ( ) N , t ) is also a solution of a many-particle Schr¨dinger o equation for non-separable Hamiltonian without special exchange symmetry requirements. Finally, the conditional wave function Ψ ( a /vector r , t a ) extracted from (28) can be written as:

$$\Psi _ { a } ( \vec { r } _ { a }, t ) = C \sum _ { n = 1 } ^ { N! } \Psi _ { n s } ( \vec { r } _ { p ( n ) _ { 1 } } [ t ],., \vec { r } _ { p ( n ) _ { e } },., \vec { r } _ { p ( n ) _ { N } } [ t ], t ) \times \text{sign} \left ( \vec { p } ( n ) \right ). \quad \text{(29)} \\ \text{Ne have substituted all positions} \, \text{hv the corresponding trajectory event} \, \text{the degree of}$$

We have substituted all positions by the corresponding trajectory except the degree of freedom /vector r p n ( ) e = /vector r a .

2.- The second step is solving each wave function Ψ ns ( /vector r p n ( ) 1 [ t , ..., /vector ] r p n ( ) e , ..., /vector r p n ( ) N [ t , t ] ) present in (29) as a solution of (24). Since Ψ ns has no exchange interaction, we can look for a solution similar to the one mentioned in the example in section 3.1.1. Then, we can write Ψ ns ≡ Ψ ns ( /vector r p n ( ) 1 [ t , ..., /vector ] r p n ( ) e , ..., /vector r p n ( ) N [ t , t ] ) as:

$$\Psi _ { n s } = \tilde { \psi } _ { p ( n ) _ { e }, a } ( \vec { r } _ { a }, t ) \exp ( z _ { p ( n ) _ { e }, a } ( t ) ),$$

where the term z p n ( ) e ,a ( ) t is the (complex) purely time-dependent term related to G ′′ a ( /vector r a [ t , X ] /vector a [ t , t ] ) + J ′′ a ( /vector r a [ t , X ] /vector a [ t , t ] ). The subindex a in ˜ ψ p n ( ) e ,a ( /vector r , t a ) specifies which are the potential U a ( /vector r , X a /vector a [ t , t ] ) used when solving (27). The other subindex p n ( ) e identifies the initial wave function, as explained in next step.

3.- The third step is finding the initial wave function ˜ ψ p n ( ) e ,a ( /vector r , a 0). When dealing with quantum transport, we can assume that the initial many-particle wave function is

located far from the active region (deep inside the reservoirs in a free space region) where it can be written as a Slater determinant (or permanent) as in (6) only during the time t = 0. Then, we can easily realize that the initial state defining ˜ ψ p n ( ) e ,a ( /vector r , a 0) ≡ ψ /vector, l ( r 0) is the particular wave packet of the ones defined in (5) which accomplishes p n ( ) e = a (see Refs. [15] and [31]).

4.The fourth step, once we know all wave functions ˜ ψ p n ( ) e ,a ( /vector r , t a ), is to compute the many-particle wave function Ψ ( a /vector r , t a ) = Ψ( /vector r 1 [ t , .../vector ] r a -1 [ t , /vector ] r , /vector a r a +1 [ t , .. ] ) in (29) as:

$$\Psi _ { a } ( \vec { r } _ { a }, \vec { X } _ { a } [ t ], t ) & = C \sum _ { n = 1 } ^ { N! } \tilde { \psi } _ { \vec { p } ( n ), a } ( \vec { r } _ { a }, t ) \times \exp \left ( z _ { \vec { p } ( n ), a } ( t ) \right ) \, \text{sign} \left ( \vec { p } ( n ) \right ). \quad \text{$31$} \\ \text{We have to specify the values of the unknown hashes $z_{1}$} \ \ ( t \, \cdots \ \text{$We will $fv$}$$

We have to specify the values of the unknown phases z p n ( ) e ,a ( ). t We will fix these phases trying to satisfy the symmetry requirements of the Bohmian trajectories discussed in section 2.2.1. In particular, we will demand that the observable results associated to different particles are indistinguishable. The following phases z /vector p n ,a ( ) ( ) t accomplish the previous symmetry condition:

/negationslash

In the Appendix A we show that this condition is enough to ensure that ensemble results of different particles are identical.

$$\mu \quad & \quad \exp \left ( z _ { \vec { p } ( n ), a } ( t ) \right ) = \prod _ { k = 1, k \neq a } ^ { N } \tilde { \psi } _ { p ( n ) _ { e }, a } ( \vec { r } _ { k } [ t ], t ). \\ \intertext { the Appendix A we show that this condition is enough to ensure that ensemble results }$$

These are the four necessary steps needed to compute an N -particle Bohmian trajectory with exchange interaction for non-separable Hamiltonians. Let us discuss the number of conditional wave functions that we need for each N -particle Bohmian trajectory. We realize that we have N possible initial wave functions ψ /vector l ( r , a 0) in (31). Since the potential U a ( /vector r , X a /vector a [ t , t ] ) is invariant under the exchange of trajectories different than /vector r a [ t ], there are only N different potentials needed. Then, when computing the N ! functions ˜ ψ p n ( ) e ,a ( /vector r , t a ) present from (31), we realize that there are many repeated solutions. Therefore, there are N × N different wave functions ˜ ψ l,a ( /vector r , t a ) that we have to solve in order to compute (31). The N × N correspond to e = 1 , ..., N different potentials and a = 1 , ...., N different initial wave packets.

In addition, it is important to notice what is the result of our algorithm when the non-separability of the Hamiltonian becomes negligible but the exchange interaction is still present. Then, we directly recover the Slater determinant (or permanent) defined in (6). Finally, we want to emphasize that the algorithm for the inclusion of the exchange interaction is universal in the sense that exactly the same 4 steps have to be followed for any system.

3.2.1. Test for non-separable harmonic potentials with exchange interaction In order to clarify the explanation of the exchange algorithm, we applied it to the same system discussed in section 3.1.1 but with the exchange interaction included. In figure 10, we have plotted the ensemble results of the (Bohmian) kinetic energy, (17), the quantum

potential energy, (18), computed directly from the two-particle wave function 2D exact solution of (2) for two identical electrons (with parallel spins). In particular, we consider fermions with the some potential and initial wave packets that we discuss in section 3.1.1. Now, the energies of particle 1 and 2 become indistinguishable. In addition, we realize that the fact that Bohmian trajectories cannot cross the diagonal of the configuration space, implies a decrease/increase of the (Bohmian) kinetic/quantum energy when the wave function crosses the diagonal. This is the same effect discussed previously in section 2.2.2.

Figure 10: (Color online) Time evolution of the total and individual (ensemble averaged) energies of identical two-electron system with exchange interaction under non-separable potential.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0388v1-figure-856940bd-5036-4608-9783-064287d3894c.jpg)

Finally, we compute the same results as in figure 10 with the 1D approximated conditional wave functions discussed in section 3.2. Now, we have to compute 4 different functions, ˜ ψ l,a ( x , a 0). The wave function ˜ ψ 1 1 , ( x , 1 0) has initial wave function ψ 1 ( x , 1 0) defined from (5) and the potential U 1 ( x , x 1 2 [ t ]) = c x ( 1 -x 2 [ t ]) 2 . The wave function ˜ ψ 1 2 , ( x , t 2 ) has the same initial state ψ 1 ( x , 2 0) but different potential energy U 2 ( x , x 2 1 [ t ]) = c x ( 1 [ t ] -x 2 ) 2 . Finally, ˜ ψ /vector 2 1 , ( x , t 1 ) has initial state ψ 2 ( x , 1 0) and potential U 2 ( x , x 2 1 [ t ]) = c x ( 1 [ t ] -x 2 ) 2 , while ˜ ψ /vector 2 2 , ( x , t 2 ) has the same initial state ψ 2 ( x , 2 0) and U 1 ( x , x 1 2 [ t ]) = c x ( 1 -x 2 [ t ]) 2 . The final wave functions Ψ ( l 1 x , t 1 ) and Ψ ( l 2 x , t 2 ) for the computation of the Bohmian trajectory, x l 1 [ t ] and x l 2 [ t ] are, respectively:

$$\Psi _ { 1 } ^ { l } ( x _ { 1 }, t ) & = C \left ( \tilde { \psi } _ { 1, 1 } ^ { l } ( x _ { 1 }, t ) \tilde { \psi } _ { 2, 2 } ^ { l } ( x _ { 2 } [ t ], t ) \ - \tilde { \psi } _ { 2, 1 } ^ { l } ( x _ { 1 }, t ) \tilde { \psi } _ { 1, 2 } ^ { l } ( x _ { 2 } [ t ], t ) \right ), \\ \Psi _ { 2 } ^ { l } ( x _ { 2 }, t ) & = C \left ( \tilde { \psi } _ { 1. 1 } ^ { l } ( x _ { 1 } [ t ], t ) \tilde { \psi } _ { 2. 2 } ^ { l } ( x _ { 2 }, t ) - \tilde { \psi } _ { 2. 1 } ^ { l } ( x _ { 1 } [ t ], t ) \tilde { \psi } _ { 1. 2 } ^ { l } ( x _ { 2 }, t ) \right ).$$

We emphasize that we require four single-particle wave functions for each 2-particle trajectory { x l 1 [ t , x ] l 2 [ t ] } . As discussed previously, the algorithm with exchange scales as N 2 . In figure 11, we have plotted the information about the energies for the same electrons discussed in section 3.1.1. The agreement between the exact 2D results and the approximate 1D ones for identical particles is acceptable. Let us emphasize that the excellent agreement in figure 9 and the results of figure 11, both, have been computed

$$\Psi _ { 2 } ^ { l } ( x _ { 2 }, t ) = C \, \overset { \bullet } { \left ( \tilde { \psi } _ { 1, 1 } ^ { l } ( x _ { 1 } [ t ], t ) \tilde { \psi } _ { 2, 2 } ^ { l } ( x _ { 2 }, t ) - \tilde { \psi } _ { 2, 1 } ^ { l } ( x _ { 1 } [ t ], t ) \tilde { \psi } _ { 1, 2 } ^ { l } ( x _ { 2 }, t ) \right ) } ^ { \prime }. \quad \text{(34)} \\ \text{We emphasize that we require four single-particle wave functions for each $2$-particle}$$

with the mentioned approximations on the terms G a and J a in (24). However, the algorithm with exchange interaction is more sensible to the approximations because we have to deal with Ψ ( l a x , t a ) that are very close to zero at the diagonal of the configuration space. A small deviation (due to the approximate G a and J a ) in the value of the modulus close to zero becomes an amplified deviation in the velocity, as seen in (10), which is inversely proportional to the modulus. This difficulty is not present in the results of figure 9 because, there, trajectories are not forced to be closer to regions where the modulus is zero. This difficulty is also manifested in the quantum potential. See the importance of the quantum potential in the 2 D exact result with exchange interaction at 1100 fs in figure 10 (red lines), while the quantum potential is negligible in the 2 D exact results without exchange interaction, as seen in figure 8 (red lines).

All (ensemble) results presented in this work can be explained in terms of individual trajectories (each trajectory with different initial conditions). Thus, the error in the ensemble results in figure 11 is due to errors in some individual trajectories, not all (mainly those trajectories starting outside of the center of the wave packet and arriving at regions where the wave function is almost zero). The total number of trajectories is 2 × 160000 (computed from 4 × 160000 conditional wave functions). A wrong trajectory will never be converted into a correct one at a later time. On the contrary, one can expect that a correct trajectory at some particular time can become a wrong one at a later time due to the approximations in G a and J a . In any case, a better approximation of the unknown terms G a and J a in (24) for systems without exchange will improve the accuracy of the algorithm with exchange. An interesting path to improve the approximations in the unknown terms G a and J a of the equation of the conditional wave function can be obtained by following Ref. [35]. It is showed there that the terms G a and J a in (24) can be computed, in principle, from a (infinite) set of coupled differential equations. A practical implementation will certainly require cutting the infinite set somewhere.

Figure 11: (Color online) Time evolution of the individual (ensemble averaged) Bohmian kinetic energies of two-electron system with exchange interaction under non-separable potential computed from 2 D exact and 1 D approximate solutions.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0388v1-figure-b558acfc-e881-4558-bbdf-da0966643231.jpg)

The most relevant feature is that the ensemble results for the first and second particle become indistinguishable with our 1D approximated conditional wave functions (see circle, plus, cross and square symbols in figure 11), although we perfectly distinguish Bohmian trajectories labeled as { x l 1 [ t , x ] l 2 [ t ] } .

## 4. Numerical results for electron transport simulation

Now, we use the general algorithm (with exchange interactions in non-separable Hamiltonians) to study quantum electron transport. The numerical implementation of the algorithm has been included into our simulator BITLLES (Bohmian Interacting Transport for non-equiLibrium eLEctronic Structures) [32]. It is a general, versatile and time-dependent 3 D electron transport simulator that allows the computation of DC, AC, transient and current and voltage fluctuations (noise) for nanoelectronic devices [15]. We compute the influence of the exchange interaction on the static and dynamic performance of a simple nano-resistor. We consider explicitly the Coulomb and exchange interaction among electrons inside the device active region. No other scattering mechanics is considered in the simulations. The details of the injection model for electrons (with a Binomial distribution) are explained in Appendix C.

## 4.1. Definition of the simulation system with arbitrary spins

We consider a set of free electrons moving under the influence of the Coulomb and exchange interaction inside a device active region, when an external bias is applied between source and drain (see figure 12). Electrons are injected into the device with an arbitrary spin. Therefore, in principle, we would have to consider several terms (each one with different spin distributions) in the many-particle wave function (1) to treat properly the exchange interaction. The consideration of many terms in (1) would imply an intractable computational burden, as mentioned in Appendix B. Alternatively, we can assume that the effect of the exchange interaction on the dynamics of an electron with spin up ↑ (down ↓ ) is due only to the other electrons with spin up ↑ (down ↓ ). Therefore, in this work we will assume that the many-particle wave function can be separated into a product of spin-up and spin-down many-particle wave functions. For the example, we will consider that the wave function with arbitrary spins Ψ ↑ 1 , ↓ 2 , ↓ 3 , ↑ 4 ( /vector r , /vector 1 r , /vector 2 r , /vector 3 r , t 4 ) γ ( ↑ 1 , ↓ 2 , ↓ 3 , ↑ 4 ), can be approximated as the product Ψ ( ↑ /vector r , /vector 1 r , t 4 ) γ ( ↑ 1 , ↑ 4 ) by Ψ ( ↓ /vector r , /vector 2 r , t 3 ) γ ( ↓ 2 , ↓ 3 ):

$$& \Psi _ { \uparrow 1, \downarrow _ { 2 }, \downarrow _ { 3 }, \uparrow _ { 4 } } ( \vec { r } _ { 1 }, \vec { r } _ { 2 }, \vec { r } _ { 3 }, \vec { r } _ { 4 }, t ) \gamma ( \uparrow _ { 1 }, \downarrow _ { 2 }, \downarrow _ { 3 }, \uparrow _ { 4 } ) \approx \\ & \Psi _ { \uparrow } ( \vec { r } _ { 1 }, \vec { r } _ { 4 }, t ) \gamma ( \uparrow _ { 1 }, \uparrow _ { 4 } ) \Psi _ { \downarrow } ( \vec { r } _ { 2 }, \vec { r } _ { 3 }, t ) \gamma ( \downarrow _ { 2 }, \downarrow _ { 3 } ).$$

From this approximation, we can apply our algorithm presented in section 3.2 to Ψ ↑ and Ψ independently. However, there are some terms in the left hand side of (35) that are ↓ not present in the right hand side (see Appendix B). The approximation (35) has already

been tested in [33]. In the mentioned reference, it is shown that this approximation is almost an exact result when the normalized distance d defined in (9) is larger than 1. For closer electrons, d smaller than 1, an error appears in the computation of the Bohmian velocity [33]. On the other hand, we explicitly consider the Coulomb interaction among all the electrons, regardless of their spin. Therefore, we can expect that electrons try to be spatially separated, i.e. large d , because of the Coulomb repulsion among electrons. This argument provides an additional justification on the validity of (35) when study quantum electron transport (with Coulomb and exchange interaction). Finally, let us notice that the number of particles varies as the simulation progress, i.e., we are dealing with a more complex many-particle wave function than that represented in (1).

In figure 12 we show a scheme of the nano-resistor that we will simulate with two GaAs source/drain doped contacts with a Fermi level of 0 15 eV above the conduction . band and a device active region of intrinsic GaAs with length L x = 30 nm. Transport takes place in the x direction with room temperature in all simulations. A single spherical band with m ∗ GaAs = 0 067 . m 0 (and m 0 the electron free mass) is considered. Because of the geometry L x /greatermuch L , L y z , with L y = L z = 9 nm, energy confinement takes place in the lateral directions. We only take into account the first energy of the subband of GaAs with a value of E 1 = 0 13 eV just above the bottom of the conduction band. .

Figure 12: (Color online) Scheme of a nano-resistor with N + AsGa source/drain doped contacts and intrinsic AsGa in the device active region. The device dimension are L x = 30 nm, L y = L z = 9 nm. A Fermi level of 0 15 . eV and room temperature are considered.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0388v1-figure-a3e28b92-659c-414a-93c2-d7212cc28227.jpg)

## 4.2. Computation of I-V characteristic

In figure 13, we present the time-averaged current 〈 I 〉 as a function of the external bias for four different scenarios: Coulomb and exchange interactions (CEI), without exchange or Coulomb interactions (WI), with Coulomb interaction alone (CI) and with exchange interaction alone (EI). It is clear that the differences observed in the different curves are mainly a direct consequence of the Coulomb interaction. Somehow, for our particular device, the Coulomb interaction screens the effect of the exchange interaction because most of the electrons (not all) are already repelled by the Coulomb interaction. We emphasize that we are using Coulomb interaction beyond mean field, with selfinteraction correction [24]. The presence of Coulomb interaction tends to reduce the current because there are less electrons in the channel. Electrons repel each other.

For each simulated electron, we know when it enters the simulation box, t l in , and when it leaves t l out . The superindex l indicates which Bohmian trajectory is associated to the electron. Identically, we know whether the electron enters (leaves) the simulation box from the source (S) or drain (D) contacts. Then, we compute:

$$d _ { A / B } & = \sum _ { l _ { A / B } } \int \Theta ( t - t ^ { l } _ { i n } ) \Theta ( t ^ { l } _ { o u t } - t ) d t, \\ \intertext { w h o r o \ t h o \ H e a v i i d o \ e t o \ f n c t i o \ i e \ \Theta ( t ) } \, - \, 1 \ f o r \ \text{activo time and } \Theta ( t ) \, - \, 0 \ f o r$$

where the Heaviside step function is Θ( ) t = 1 for positive times and Θ( ) t = 0 for negatives ones. The time integral is over the whole simulation time. The sum is over all l -trajectories that have entered through the contact A = { S, D } and leave through B = { S, D } . In figure 14(a), we plot d S/D /d and d D/S /d , and in figure 14(b) d S/S /d and d D/D /d , where we have defined d = d S/D + d D/S + d S/S + d D/D . Electrons crossing the active region, d S/D and d D/S , are responsible for the DC (i.e. zero frequency) behavior. On the other hand, electrons that do not cross the device active region, d S/S and d D/D , do not contribute to DC, but only to high-frequency dynamics.

At zero bias, as seen in figure 14(a), without interaction (WI), half of electrons are transmitted from source to drain and half from drain to source. No reflected electrons. This is not true for the rest of scenarios. The EI simulation provides reflected electrons because of the effect seen in section 2.2.2, that forbids electrons from occupying the same positions (i.e. the diagonals points of the configuration space) and some of the electrons are finally bounced. Let us emphasize that the mean number of electrons in the active region of the quantum wire of figure 12 can be very small. From the current in figure 13, we can compute the rate of transmitted electrons which is quite similar to their transit time. The CI (and CEI) simulation shows reflected because of the Coulomb repulsion among electrons. Finally, we focus on set of plots d D/D /d (from drain to drain) in figure 14(b). At a 0.05 V, the EI simulation has a larger value d D/D /d than the others. This result does not implies a larger number of reflected particles with EI simulation, but only that these EI reflected particles spent more time inside the active region than, for example, the CI and CEI ones. The reflection in the EI simulation occurs when the particles are really very close, see figure 6(a), while the reflection in the CI simulation occurs for particles with a larger spatial separation. The Coulomb interaction has a longer range than the particular (exchange interaction) effect seen in figure 6(a). Thus, these CI and CEI particles spent less time in the active region. This result justifies again why the Coulomb interaction, somehow, screens the possible effects of the exchange interaction.

## 4.3. Computation of the noise

From the time-dependent current I ( ) provided by the Monte Carlo t BITLLES simulator we can also study the noise characteristics of this simple nano-resistor in a very simple way [34]. The fluctuations of the current can be easily obtained from the autocorrelation function R τ ( ) = ∆ I ( )∆ t I ( t + ) of the current fluctuations ∆ τ I ( ) = t I ( ) t -〈 〉 I being I ( ) t the instantaneous current provided by the quantum Monte Carlo BITLLES simulator

Figure 13: (Color online) The average current as a function of the applied bias for the simulated system in four different situations: Coulomb and exchange interactions, without interactions, Coulomb interaction, and exchange interaction.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0388v1-figure-ba9ee5d5-9b41-4222-ba07-d3e942034f37.jpg)

Figure 14: (Color online) Normalized (mean) time spent by the electrons inside the simulating box as a function of the applied voltage for the four different scenarios discussed in figure 13. (a) from drain to source, d D/S /d , and from source to drain, d S/D /d . (b) from the drain that have been finally bounced, d D/D /d and from the source that have been finally bounced, d S/S /d

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0388v1-figure-9cd8d7d1-093d-451c-ac48-c24167f9724f.jpg)

and 〈 I 〉 the time-average DC current computed above. The Fourier transform of this autocorrelation function R τ ( ) is the noise power spectral density of S f ( ) = ∫ ∞ -∞ R τ e ( ) -i 2 πfτ dτ . Finally, the Fano factor, defined as the ratio γ = S (0) /S schottky with S Schottky ( f ) = 2 q 〈 I 〉 , can be computed.

In figure 15 we show the noise power spectral density as a function of the frequency

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0388v1-figure-4b7531e1-83a1-486d-9a7c-f60d45e5a68e.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0388v1-figure-7099da66-7efd-4087-bbcc-acf99aac9f49.jpg)

Figure 15: (Color online) Noise power spectral density as a function of the frequency for four applied bias. (a) V DS = 0.01 V (b) V DS = 0.02 V (c) V DS = 0.1 V (d) and V DS = 0.2 V.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0388v1-figure-9b964944-e004-42c9-9ba3-27736350a21c.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0388v1-figure-60405414-b5cc-4d47-8f38-6229d51f4bf5.jpg)

for four different bias. The Fano factor decreases as we increase the voltage. Since we are simulating at room temperature, for low bias tending to zero, the current tends to zero while the thermal noise is still present giving a Fano factor tending to infinite. For high bias, the fluctuations follows the Binomial distribution (C.2) that roughly tends to a Poissonian one for electrons with energy above the Fermi level (as is the case for the sample with confinement discussed here). Deviations from this behavior are due to correlations originated from Coulomb and exchange interactions. First, we observe in figure 15 that, similarly to DC, the the most important effect on the Fano factor is due to the Coulomb interaction (see CI and CEI in figure 15). As we see in figure 15(d), the Fano factor is lower than one with CEI or CI simulations. It is well-known that the Coulomb repulsion between electrons tends to space them more regularly rather than strictly at random, and to evidence a sub-Poissonian statistics [36, 37]. However, we see that the small deviations between the CEI and CI results are influenced by the exchange interaction present in the simulation box. When EI is larger than CI, the CEI is larger than CI, and vice-verse. Secondly, we observe that the effect of the exchange interaction in the noise is more important at low bias. This result has the same explanation that we explain in the previous figure 13 and figure 14. Thirdly, in figure 15(b) and (c), at intermediate bias, we see a peak of the noise power spectral density above 1THz. One can realize that this peak in the fluctuations appears when only WI or EI are considered. The origin of this peak is the increment of the number of reflected electrons because

they cannot occupy identical positions. These bounced electrons do not affect the zero (low) frequency fluctuations, but they affect the high-frequency values. This increment of reflected particles can be seen in figure 14(b).

Finally, as discussed in the injection process in Appendix C, we want to mention that wave packets with identical central positions and wave vectors are injected with a temporal separation between them equal to t 0 that tries to avoids relevant exchange interaction effects among them inside the simulating box.

## 5. Discussion and conclusions

In this work we have presented an algorithm for introducing the exchange interaction into the many-particle quantum (Bohmian) trajectories with an universal protocol valid for any quantum system, with separable or non-separable Hamiltonians, for either fermions or bosons. In principle, one could thought in introducing the exchange interaction into the terms G a and J a present in (24) that define the conditional wave function, in a similar way as the exchange-correlation functional is introduced into the Hamiltonian of a single-particle pseudo-Schr¨dinger equation used by DFT. However, There is no clear o prescription on how to define directly the terms G a and J a with exchange interaction. Alternatively, in this work we follow a different path. We compute N × N conditional wave function solutions of (24) [with a very simple approximation for G a and J a , leading to (27)] without symmetry requirements. However, the global conditional wave function constructed as a combination of them, (31), can satisfy the exchange requirements, after a proper guess for some phases. We have shown that the phases defined from (32) satisfy, by construction, three very satisfactory properties. First, in the case of separable Hamiltonians, it directly leads to the standard Slater determinants for computing manyparticle wave functions (or the permanent for bosons). Second, for fermions with nonseparable Hamiltonians, it guarantees that the probability presence of the many-particle wave function in the diagonal points of the configuration space is zero. Third, in any scenario, the Bohmian trajectories satisfies the expected symmetry property when initial positions are interchanged. This last property ensures that observable results computed from Bohmian trajectories are indistinguishable, as seen in figure 11. Note that the selection of the phases in (32) is not unique. In the Appendix A, we have shown another possibility, but it does not satisfy the second property discussed above.

An improvement on the simple approximation used to compute G a and J a when constructing the conditional (Bohmian) wave functions without exchange interaction would also improve the accuracy of the algorithm with exchange presented here.

As a practical demonstration of the numerical viability of the algorithm discussed here for quantum transport (in a far-from equilibrium open system), the current and its fluctuations are computed for a nano-resistor, with exchange and Coulomb interactions. For this simple device, the effects of the exchange interaction are mainly screened by the Coulomb interaction. In any case, the main conclusions that one can extract from these numerical results, applied to this very simple device, is that the algorithm explained in

this work can be perfectly implemented for a number of electrons on the order of 20-30. This requires solving 30 × 30 ≈ 1000 single particle conditional wave functions, which can still be handled with normal computing facilities and simulation times on the order of few hours for each simulation bias.

## acknowledgement

The present work have greatly beneficed from discussions with W. Struyve and G. Albareda. This work was supported through Spanish project MICINN TEC2012-31330.

## Appendix A. Indistinguishable results computed from the approximate conditional Bohmian wave functions using (32)

The requirement for ensuring that ensemble results computed from Bohmian trajectories are indistinguishable is that the N -particle Bohmian trajectories are symmetric under the interchange of their initial positions: We consider two different N -particle Bohmian trajectories, whose initial positions are /vector X l [0] and /vector X f [0]. In particular, we consider /vector r l k [0] = /vector r f k [0] for all k except /vector r l j [0] = /vector r f h [0] and /vector r l h [0] = /vector r f j [0]. Then, the sufficient condition to ensure that observable results of different particles are indistinguishable is ensuring that /vector r l k [ t ] = /vector r f k [ t ] = /vector r k [ t ] for all k except /vector r l j [ t ] = /vector r f h [ t ] and /vector r l h [ t ] = /vector r f j [ t ]. In this appendix, we show that this last property is guaranteed by our proposal of using (32) for fixing the unknown phases z p n ( ) e ,a ( ). t

At the initial time, by construction (see step 3 in section 3.2), it is obvious that (32) provides Bohmian trajectories with the desired property. Let us demonstrate that this relation between trajectories at the initial time is also true at a later time. We assume that the Bohmian trajectories have the desired property at the time t ′ (for example, t ′ = 0). Then, we want to demonstrate that /vector r l j [ t ] = /vector r f h [ t ] and /vector r l h [ t ] = /vector r f j [ t ] at a later time t = t ′ + dt . In fact, we only have to demonstrate that the velocities of the particles j and h at t ′ are interchanged. According to (32), the conditional wave function used to compute the velocity /vector v l j ( /vector X,t ′ ) | /vector X = /vector X t l [ ′ ] with the initial positions /vector X t l [ ′ ] = { /vector r 1 [ t ′ ] , ..., /vector r α [ t ′ ] , ..., /vector r β [ t ... ] } is:

$$\Psi _ { j } ^ { l } ( \vec { r } _ { j }, t ^ { \prime } ) & = C \sum _ { n = 1 } ^ { N! } \tilde { \psi } _ { p ( n ) _ { e }, \alpha } ( \vec { r } _ { 1 } [ t ], t ^ { \prime } ),.., \tilde { \psi } _ { p ( n ) _ { e }, \alpha } ( \vec { r } _ { j }, t ^ { \prime } ) \\ &,.., \tilde { \psi } _ { p ( n ) _ { e }, \alpha } ( \vec { r } _ { \beta } [ t ], t ^ { \prime } ).., \tilde { \psi } _ { p ( n ) _ { e }, \alpha } ( \vec { r } _ { N } [ t ], t ^ { \prime } ) \times \text{sign} \left ( \vec { p } ( n ) \right ), \quad \quad \quad$$

where we have defined /vector r l j [ t ] = /vector r f h [ t ] ≡ /vector r α [ t ] and /vector r l h [ t ] = /vector r f j [ t ] ≡ /vector r β [ t ]. Identically, we have /vector v l h ( /vector X,t ′ ) | /vector X = /vector X t l [ ′ ] :

$$\Psi _ { h } ^ { t } ( \vec { r } _ { h }, t ^ { \prime } ) & = C \sum _ { n = 1 } ^ { N! } \tilde { \psi } _ { p ( n ) _ { e }, \beta } ( \vec { r } _ { 1 } [ t ], t ^ { \prime } ),.., \tilde { \psi } _ { p ( n ) _ { e }, \beta } ( \vec { r } _ { \alpha } [ t ], t ^ { \prime } ) \\ &,.., \tilde { \psi } _ { p ( n ) _ { e }, \beta } ( \vec { r } _ { h }, t ^ { \prime } ).., \tilde { \psi } _ { p ( n ) _ { e }, \beta } ( \vec { r } _ { N } [ t ], t ^ { \prime } ) \times \text{sign} \left ( \vec { p } ( n ) \right ).$$

On the other hand, for the many-particle Bohmian trajectory with initial conditions /vector X f [0] = { /vector r 1 [0] , ..., /vector r β [0] , ..., /vector r α [0] ... } , we have for /vector v f j ( /vector X,t ′ ) | /vector X = /vector X f [ t ′ ] :

$$\Psi _ { j } ^ { f } ( \vec { r } _ { j }, t ^ { \prime } ) & = C \sum _ { n = 1 } ^ { N! } \tilde { \psi } _ { p ( n ) _ { e }, \beta } ( \vec { r } _ { 1 } [ t ], t ^ { \prime } ),.., \tilde { \psi } _ { p ( n ) _ { e }, \beta } ( \vec { r } _ { j }, t ^ { \prime } ) \\ &,.., \tilde { \psi } _ { p ( n ) _ { e }, \beta } ( \vec { r } _ { \alpha } [ t ], t ^ { \prime } ).., \tilde { \psi } _ { p ( n ) _ { e }, \beta } ( \vec { r } _ { N } [ t ], t ^ { \prime } ) \times \text{sign} \left ( \vec { p } ( n ) \right ).$$

and for /vector v f h ( /vector X,t ′ ) | /vector X = /vector X f [ t ′ ] :

$$\Psi _ { h } ^ { f } ( \vec { r } _ { h }, t ^ { \prime } ) & = C \sum _ { n = 1 } ^ { N! } \tilde { \psi } _ { p ( n ) _ { e }, \alpha } ( \vec { r } _ { 1 } [ t ], t ^ { \prime } ),.., \tilde { \psi } _ { p ( n ) _ { e }, \alpha } ( \vec { r } _ { \beta } [ t ], t ^ { \prime } ) \\ &,.., \tilde { \psi } _ { p ( n ) _ { e }, \alpha } ( \vec { r } _ { h }, t ^ { \prime } ).., \tilde { \psi } _ { p ( n ) _ { e }, \alpha } ( \vec { r } _ { N } [ t ], t ^ { \prime } ) \times \text{sign} \left ( \vec { p } ( n ) \right ). \\ & \dots \dots$$

As expected, it becomes obvious that /vector v l j ( /vector r , t j ′ ) = /vector v f h ( /vector r h , t ′ ) and /vector v l h ( /vector r h , t ′ ) = /vector v f j ( /vector r , t j ′ ). Q.E.D.

/negationslash

/negationslash

The use of expression (32) does also provide the additional property that Ψ ( l j /vector r , t j ) becomes zero whenever /vector r j = /vector r k [ t ] for j = k . It can also be shown that exp ( z /vector p n ,a ( ) ( )) = t ∏ N k =1 ,k = a ˜ ψ p n ( ) e ,k ( /vector r k [ t , t ] ) does also provide the required symmetry condition for the trajectories, but it does not guarantee that the conditional wave function is zero for fermions at the coincident points. Thus, (32) is preferred. Let us emphasize that ensuring that observable results become indistinguishable is not enough to ensure that our algorithm provides the correct results. Apart from the requirement mentioned here, Bohmian trajectories have to satisfy other properties as for example their non-crossing property discussed in section 2.2.1.

## Appendix B. Exchange interaction for electrons with different spin

Even for systems without spin-orbit interaction and when we are not interested in the time evaluation of the spins, we cannot neglect the spin degrees of freedom of the electrons (with arbitrary spins) because the symmetry of the overall wave function depends on the exchange properties of the orbital part and spin component. To understand the complexity of computing the antisymmetry wave function with spins in different directions, we present an example for three electrons, one with spin up ( ↑ j ) and the others two with spin down ( ↓ j ). Then, the global antisymmetric wave function in (1), for this particular case, can be written as:

$$\wedge \cdot \cdot \cdot \cdot \cdot \cdot \cd$$

We define the orbital wave functions ψ l ( x l ) as the Gaussian wave packets of (5). Equation (B.1) has 3! terms, each one composed of the product of an orbital function by

a spin function. Next, we compute the total norm taking into account the 3!3! products of permutations. For this purpose we have to multiply the orbital parts and the spin parts separately. Due to orthogonality, the product of the spin part can be either 0 and 1. The final result is:

$$\begin{array} { l l$$

Let us notice that, in principle, we had to keep 3!3! = 6 2 = 36 terms. However, only these terms whose product of spin parts is 1 are present in (B.3). The evaluation of the product of the spin parts have to be done explicitly, term by term, with no possibility of simplification. For example, the product of γ ( ↑ 1 , ↓ 2 , ↓ 3 ) by γ ( ↑ 1 , ↓ 3 , ↓ 2 ) is 1, while the product of γ ( ↑ 1 , ↓ 2 , ↓ 3 ) by γ ( ↓ 2 , ↑ 1 , ↓ 3 ) is 0. However, if we increase the number of electrons, the practical computation of the previous expression is computationally inaccessible. Note that N = 8 gives 8! 2 = 40320 2 terms.

In (35) we provide a (computationally accessible) approximation to treat wave functions with spin of different orientations. For example, if we assume that there is no exchange interaction between spin up and spin down components, then:

$$& \bar { \Phi } ( x _ { 1 }, x _ { 2 }, x _$$

Now, the norm would be:

$$\text{row, the norm would be: } \\ & \left | \bar { \Phi }$$

where there are terms on (B.3) that are not present in (B.4). However, in Ref. [33] we show that the approximation (35) provides a quite reasonable approximation for the computation of the Bohmian velocity. Expression (B.3) keep the most relevant terms of the exchange interaction.

## Appendix C. Electron injection probability

Our Bohmian algorithm requires each electron to be described by a (conditional) wave function plus a Bohmian trajectory. Every time an electron with a particular initial wave function is selected to enter the device active region, an initial position for the Bohmian trajectory associated to this wave packet has to be randomly selected according to 'quantum equilibrium hypothesis' [14, 28] discussed in section 2.2. Accordingly, this initial position is more frequently found to be around the center of the wave packet than in the borders. Next, we explain how the wave packets are selected.

The (central) kinetic energy E o of the wave packet is related to the (central) wave vectors k o by E o = ( /planckover2pi1 k o ) 2 / (2 m ∗ ) with m ∗ the effective mass. For each contact, we select a flat potential region in an (non-physical) extension of the simulation box (for x &lt; 0 in the source and x &gt; L x in the drain of figure 12). Let us notice that in a flat potential region without interaction, a single-particle wave packets is exactly the normalized conditional (Bohmian) wave function discussed in this work. The initial Gaussian wave packet is defined in this flat potential region (deep inside the contact) following the analytical expression (5). In particular, the central position of all wave packets is selected x o = 100 nm far from the border of the active region (inside the contacts) and the spatial dispersion of the wave packet is σ x = 25 nm (i.e. the wave packet is somehow similar to a scattering state). The only two additional parameters that we still have to fix to fully define the wave packet are the central kinetic energy E o of the wave packet and the injecting time when the electron effectively enters the simulation box. The selection of the energy E o has to satisfy the Fermi-Dirac occupation function f ( E o ) that depends on the (quasi) Fermi energy and temperature. The selection of the time when the electron is injected is a a bit more complex.

Let us define t 0 as the minimum temporal separation between the injection of two wave packets whose central wave vectors and central positions fit into the following particular phase-space cell k o ∈ [ k , b k b + ∆ ) and k x c ∈ [ x , x b b + ∆ ), being x x b the left border of the simulation region. For a 1D system, the value of t 0 can be easily estimated. The number of electrons n 1 D in the particular phase space cell ∆ k · ∆ x is n 1 D = 2 · ∆ k · ∆ x/ (2 π ) where the factor 2 takes into account the spin degeneracy. These electrons have been injected into ∆ x during the time interval ∆ t defined as the time needed for electrons with velocity v x = ∆ x/ ∆ = t /planckover2pi1 k /m o to travel a distance ∆ x . Therefore, the minimum temporal separation, t 0 , between the injection of two electrons

into the previous cell is ∆ t divided by the maximum number n 1 D of electrons:

$$t _ { 0 } = \frac { \Delta t } { n _ { 1 D } } = \$$

It is very instructive to understand the minimum temporal separation t 0 in (C.1) as a consequence of the wave packet version of the Pauli principle discussed in section 2.1.1. The simultaneous injection of two electrons with similar central positions and central momentums would require such a huge amount of energy that its probability is almost zero (see figure 1). In other words, a subsequent electron with central position and central momentum equal to the preceding ones can only be injected after a time interval given by t 0 .

The injection of electrons (from the mentioned phase-space cell) at multiple times of t o depends finally on the statistics imposed by the Fermi-Dirac function mentioned above. During each attempt of injection at multiples of t o , we select a random number r , and the electron is effectively injected only if f ( E o ) &gt; r . The mathematical definition of the rate and randomness of the injection process are given by the following binomial probability P E , N , τ ( o τ ) (See Ref. [38]):

$$P ( E _ { o }, N _ { \tau }, \tau ) = \frac { M _ { \tau }! } { N _ { \tau }! \cdot ( M _ { \tau } - N _ { \tau } )! } f ( E _ { o } ) ^ { N _ { \tau } } \left ( 1 - f ( E _ { o } ) \right ) ^ { M _ { \tau } - N _ { \tau } }.$$

This expression defines the probability that N τ electrons (from the mentioned phasespace cell) are effectively injected into the active region during the time interval τ . The parameter M τ is the number of attempts of injecting electrons during this time interval τ , defined as a natural number that rounds the quotient τ/t o to the nearest natural number towards zero. The number of injected electrons can be N τ = 1 2 , , .... ≤ M τ . More details can be found in Ref. [38]. Finally, let us clarify that (C.2) does only specify the injecting probability. The transmission probability with a certain energy depends on the injecting probability and also on all complex (Coulomb and exchange) phenomena explained along the text.

## References

- [1] 'The general theory of quantum mechanics is now almost complete. The underlying physical laws necessary for the mathematical theory of a large part of physics and the whole of chemistry are thus completely known,f and the difficulty is only that the exact application of these laws leads to equations much too complicated to be soluble.' Dirac P A M 1929 Proceedings of the Royal Society of London A 123 714
- [2] Nightingale M P Umrigar C J 1999 Quantum Monte Carlo Methods in Physics and Chemistry (Springer)
- [3] Suzuki M 1993 Quantum Monte Carlo Methods in Condensed Matter Physics (World Scientific)
- [4] Hartree D R 1928 The Wave Mechanics of an Atom with a Non-Coulomb Central Field. Part II. Some Results and Discussion Mathematical Proceedings of the Cambridge Philosophical Society 24 (01) 111
- [5] Fock V 1930 N¨herungsmethode zur L¨sung des quanten mechanichanischen Mehrk¨rperproblems a o o Z. Phys. 61 126
- [6] Kohn W 1998 Reviews of Modern Physics 71(5) 1253

- [7] Hohenberg P and Kohn W 1964 Inhomogeneous electron gas Physical Review 136 B864
- [8] Kohn W and Sham L J 1965 Phys.Rev. A 140 1163
- [9] Soler J M Artacho E Gale J D Garcia A Junquera J Ordejon P and Sanchez-Portal D 2002 Journal of Physics: Condensed Matter 14 2745
- [10] Runge E and Gross E K U 1984 Physical Review Letters 52 12 997
- [11] Giuliani G F and Vignale G 2005 Quantum Theory of electron Liquid (Cambridge University Press)
- [12] Capelle K and Gross E K U 1997 Phys. Rev. Lett. 78(10) 1872
- [13] Di Ventra M and D'Agosta R 2007 Phys. Rev. Lett. 98 22 226403
- [14] Holland P R 1993 The Quantum Theory of Motion (Cambridge University Press Cambridge)
- [15] Oriols X and Mompart J 2011 Applied Bohmian Mechanics: From Nanoscale Systems to Cosmology (Pan Stanford Publishing Singapore)
- [16] Bell J S 1987 Speakable and Unspeakable in Quantum Mechanics (Cambridge University Press Cambridge)
- [17] D¨ urr D and Teufel S 2009 Bohmian Mechanics: The Physics and Mathematics of Quantum Theory (Spinger Germany)
- [18] Oriols X 2007 Phys. Rev. Let. 98 066803
- [19] Wyatt R E 2005 Quantum Dynamics with Trajectories: Introduction to Quantum Hydrodynamics (Springer New York)
- [20] Wyatt R E and Bittner E R 2000 Journal of Chemical Physics 113 22
- [21] Goldfarb Y Degani I and Tannor D J 2006 Journal of Chemical Physics 125 231103
- [22] Gindensperger E Meier C and Beswick J A 2002 Journal of Chemical Physics 116 8
- [23] In this work, the quantum trajectories are only used to reproduce ensemble results without any ontological implication. Then, the adjective Bohmian and hydrodynamic are fully interchangeable in this work. However, we use the first one only because the concept of conditional wave function, which is fundamental in this work, is mainly developed for Bohmian mechanics.
- [24] Albareda G Su˜´ e J and Oriols X 2009 n Phys. Rev. B 79 075315
- [25] We use the basis of individual spins instead of the basis of the total spin because we are interested in open systems, where particles (spins) are interchanged with reservoirs. Then, the requirement of constant total spin (implicit in the total spin basis) is not at all evident.
- [26] Cohen-Tannoudji C Diu B and Laloe F 1978 Quantum Mechanics (Vol I and II) (Wiley, John and Sons)
- [27] Oriols X 2004 Nanotechnology 15 S167
- [28] D¨ urr D et al. 2004 Journal of Statistical Physics 116 9595
- [29] Oriols X Mart´ ın F and Su˜´ J 1996 ne Physical Review A 54 2595
- [30] Brown H R Sj¨ oqvist E 1999 Physics Letters A 251(4) 229
- [31] For example, the permutation /vector p n ( ) = { 3 2 1 , , } on the conditional wave function used to compute the trajectory /vector r 1 [ t ] has the initial wave packet ψ 3 ( /vector r 1 , t ) because p n ( ) 3 = 1.
- [32] http://europe.uab.es/bitlles
- [33] Alarc´n A Cartoix` X and Oriols X 2010 o a Physica Status Solidi C 7 11
- [34] By construction, an ensemble of Bohmian trajectories reproduces the probability presence of the many-particle wave function at any time and also the correlations. The only approximation in our computations of the current is that the (continuous) time-dependent measurement of the current provides a negligible back-action on the simulated system [15].
- [35] Norsen T 2010 Foundations of Physics 40 1858
- [36] Nazarov Y V and Blanter Y M 2009 Quantum Transport: Introduction to Nanoscience (Cambridge University press,f Cambridge)
- [37] Bulashenko O M Mateos J Pardo D and Gonz´ alez T Reggiani L Rub´ ı J M 1998 Phys. Rev. B 57 1366.
- [38] Oriols X Fernandez-Diaz E Alvarez A Alarc´n A 2007 o Solid State Electronics 51 306
A. Alarcón S. Yaro X. Cartoixà X. Oriols
Computation of many-particle quantum trajectories with exchange interaction: Application to the simulation of nanoelectronic devices
Following Ref. [Oriols X 2007 Phys. Rev. Lett., 98 066803], an algorithm to
deal with the exchange interaction in non-separable quantum systems is
presented. The algorithm can be applied to fermions or bosons and, by
construction, it exactly ensures that any observable is totally independent
from the interchange of particles. It is based on the use of conditional
Bohmian wave functions which are solutions of single-particle
pseudo-Schr\"odinger equations. The exchange symmetry is directly defined by
demanding symmetry properties of the quantum trajectories in the configuration
space with a universal algorithm, rather than through a particular
exchange-correlation functional introduced into the single-particle
pseudo-Schr\"odinger equation. It requires the computation of N^2 conditional
wave functions to deal with N identical particles. For separable Hamiltonians,
the algorithm reduces to the standard Slater determinant for fermions, or
permanent for bosons. A numerical test for a two-particle system, where exact
solutions for non-separable Hamiltonians are computationally accessible, is
presented. The numerical viability of the algorithm for quantum electron
transport (in a far-from equilibrium time-dependent open system) is
demonstrated by computing the current and fluctuations in a nano-resistor, with
exchange and Coulomb interactions among electrons.