Spatial Entanglement of Fermions in One-Dimensional Quantum Dots

The time-dependent quantum Monte Carlo method for fermions is introduced and applied in the calculation of the entanglement of electrons in one-dimensional quantum dots with several spin-polarized and spin-compensated electron configurations. The rich statistics of wave functions provided by this method allow one to build reduced density matrices for each electron, and to quantify the spatial entanglement using measures such as quantum entropy by treating the electrons as identical or distinguishable particles. Our results indicate that the spatial entanglement in parallel-spin configurations is rather small, and is determined mostly by the spatial quantum nonlocality introduced by the ground state. By contrast, in the spin-compensated case, the outermost opposite-spin electrons interact like bosons, which prevails their entanglement, while the inner-shell electrons remain largely at their Hartree–Fock geometry. Our findings are in close correspondence with the numerically exact results, wherever such comparison is possible.


Introduction
During the past few decades there has been an increasing interest in developing new models and computational tools to address the fundamental and practical challenges related to quantum correlations and entanglement, in connection with their potential application in newly emerging quantum technologies [1]. The properties of composite systems of quantum particles are expected to play an important role in information processing, as well as in devices for manipulating systems of atoms and molecules. While various algebraic operator methods have been used to characterize entanglement in spin systems [2][3][4], an efficient approach to assess the spatial entanglement in many-body quantum systems together with its evolution over time is still lacking. It is well known that the correlated non-relativistic particle motion described by the time-dependent Schrödinger equation (SE) is tractable for only a limited number of cases. While solvable numerically for few particles in 1D and 2D, the direct numerical solution of the SE scales exponentially with system size, and is therefore beyond the capabilities of today's computers. That exponential time scaling is usually attributed to the nonlocal quantum effects that result from the dependence of the wave function Ψ(r 1 , . . . , r N , t) on the coordinates of all interacting particles. The standard approaches to ameliorate the workload are to reduce the many-body SE to a set of coupled single-body equations, where the most prominent are the mean-field approaches: the Hartree-Fock (HF) method [5] and density-functional theory (DFT) [6]. That reduction, however, occurs at the price of neglecting the detailed fluctuating forces between the electrons and replacing them with averages, thus totally ignoring the dynamic quantum correlations in the HF method, while DFT reduces the many-body problem to a single-body problem of non-interacting electrons moving in an effective exchange-correlation potential, which is generally unknown and suffers self-interaction issues. More accurate but more computationally expensive are the multi-configuration time-dependent Hartree-Fock method [7], and the full configuration interaction method [8]. Another approach that has gained much attention lately is the density matrix renormalization group method [9], which allows one to treat correlated 1D many-body problems with good accuracy; however, its application to higher dimensions and arbitrary potentials has been challenging thus far.
A different class of methods to tackle the quantum many-body problem includes the quantum Monte Carlo methods [10], which allow one to accurately calculate the electronic structures of atoms, molecules, nanostructures, and condensed-matter systems at a fully correlated level. For example, the diffusion quantum Monte Carlo (DMC) method uses random particles (walkers) whose evolution towards the ground state of the system involves a combination of diffusion and branching-which, however, prevent its use for real-timedependent processes where causality is of primary importance. Moreover, the artificial nature of the many-body wave function in configuration space used in the DMC method prevents the calculation of some important quantities other than the energy. The recently introduced time-dependent quantum Monte Carlo (TDQMC) method [11][12][13] employs concurrent ensembles of walkers and wave functions for each electron, where each wave function is associated with a separate walker (particle-wave dichotomy), and these evolve in physical spacetime where no initial guess for the many-body wave function is needed. Recently, we have applied the TDQMC method to analyze the ground state preparation for simple bosonic systems with several particles in one and two dimensions, with a good tradeoff between scaling and accuracy [14]. In this work we apply the TDQMC method to several interacting fermions in one-dimensional quantum dots, where the obtained reduced density matrices allow us to quantify the entanglement between the electrons considered to be identical or distinguishable particles. We consider some proof-of-principle aspects of the TDQMC method for fermions, rather than practical matters concerning quantum dots. Entanglement in two-electron quantum dots, atoms, and molecules has been studied elsewhere [15][16][17][18].

Methods
The TDQMC method transforms the standard Hartree-Fock (HF) equations into a set of coupled stochastic equations capable of describing the correlated particle motion. That transformation is based on the physical assumption that the modulus square of the single-body wave function in coordinate (physical) space may be thought of as an envelope (or kernel density estimation) of the distribution of a finite number of particles (walkers). In this way, for each electron in an atom a large set of single-body wave functions that reside in physical spacetime is created, where each wave function responds to the multicore potential due to both the nucleus and the walkers of the rest of the electrons. The crucial point in this picture is that it allows for each walker for a given electron to interact with the walkers of any other electron through weighted Coulomb potential, thus naturally incorporating the spatial quantum nonlocality. Then, from the evolution of the walker distributions one can evaluate quantum observables without resorting to the many-body wave function. Formally, we start from the HF equation for the i th from a total of N electrons, within a single-determinant ansatz [19]: where V en (r i ) is the electron-nuclear potential, and the HF electron-electron potential reads: where: is the Hartree potential, and V X ee (r i , t) is the exchange potential: where ϕ i (r, t) satisfy the orthonormality property ϕ i (r, t)ϕ * j (r, t)dr = δ i,j , and the indices s i , s j denote the spins of the corresponding electrons. The inequality j = i in the sums of Equations (3) and (4) stresses the fact that even though the self-interaction between the electrons is naturally canceled in the HF approximation, it is also not present in the Hartree approximation, where there is no exchange potential [19]. It is known that the wave functions ϕ i (r, t) of Equation (1) variationally minimize the system energy: where the Hartree and exchange energies read: and respectively. It is known that the Hartree-Fock approximation does not account for the dynamic electron-electron correlations beyond those due to the exchange interaction. In order to correct for this in the TDQMC methodology, we replace the HF wave function for each electron ϕ i (r, t) with a family of slightly different wave functions ϕ i (r, t) → ϕ k i (r, t) ; k = 1, . . . ,M [11][12][13], which allows us to further lower the system energy below the HF level. This is accomplished by applying a stochastic windowing to the distribution ϕ j (r j , t) 2 in the Hartree potential V H ee (r i , t) of Equation (3), by using a "window" function K r j , r k j (t), σ j,i centered at a certain walker's trajectory r k j (t), which samples the distribution given by ϕ j (r j , t) 2 . The parameters σ j,i determine the widths of those "windows", such that the product ϕ j (r j , t) 2 K r j , r k j (t), σ j,i is different for each separate trajectory r k j (t). In this way, for each electron, Equations (1)-(4) are transformed into a set of M Hartree-Fock-like equations for the different replicas ϕ k i (r i , t) of the initial HF wave function ϕ i (r i , t), each one attached to a separate trajectory r k j (t) (particle-wave dichotomy [13]): where ϕ k i (r, t)ϕ k * j (r, t)dr = δ i,j ; i = 1, . . . ,N, k = 1, . . . ,M, and where: is the effective electron-electron interaction potential represented as a Monte Carlo (MC) convolution that incorporates the spatial quantum nonlocality by allowing each walker for a given electron to interact with a group of walkers of any other electron. In fact, it can be seen from Equation (9) that the effective potential "seen" by the kth wave function for the ith electron involves interactions with a number of walkers that belong to the jth electron and lie within the nonlocal length σ j,i around r j (t) [12][13][14]. For the Gaussian kernel we have: which determines the weighting factor in Equation (9) to be: As seen from Equations (10) and (11), the limit σ j,i → ∞ where K r j , r k j (t), σ j,i → 1 recovers the Hartree-Fock approximation-as opposed to the local interaction, where It is clear, therefore, that σ j,i may serve as variational parameters to minimize the system energy between these two limiting cases.
The connection between the trajectories r k i (t) and the wave functions ϕ k i (r i , t) is given by the walker's velocities (de Broglie-Bohm equation, e.g., in [20] for real-time propagation, and: for imaginary-time propagation, where the drift velocity reads: and η(τ) is a Markovian stochastic process (see also the appendix in [21]). The striking similarity between the drift velocity of Equation (14) and the de Broglie-Bohm Equation (12) comes from the fact that both equations describe drift-diffusion processes in imaginary and in real time, respectively. It is seen that although the individual wave functions guide the corresponding walkers through Equation (12), the TDQMC method solves coupled single-body Hartree-Fock-like equations (e.g., Equation (8)) instead of using quantum potentials as in Bohmian mechanics [20]. Following the particle-wave dichotomy described above, the system energy can be calculated conveniently using both particle trajectories and wave functions: During the preparation of the ground state of the quantum system, the initial Monte Carlo ensembles of walkers and wave functions propagate in imaginary time (τ) toward a steady state, in accordance with Equations (8)- (14), for different values of the nonlocality parameters σ j,i , until the energy of Equation (15) displays a minimum. Another alternative to account for the exchange effects is to apply a short-range screening to the interaction potential in order to modify the repulsion between the same-spin electrons [22].
Considering the ensemble of waves ϕ k i (r i , t) delivered by the TDQMC as random variables, one can build a reduced density matrix for the ith electron, which may serve as the variance-covariance matrix in the Hilbert space that carries important statistical information [13,23]: For example, the density matrix of Equation (16) allows one to easily calculate the spatial entanglement of a given electron state, which may serve also as a good measure for the overall accuracy of the calculation (e.g., [24]). Without entering the ongoing debate on entanglement witnesses, for opposite-spin electrons-where there are no exchange terms in HF and TDQMC equations-we employ the linear quantum entropy for distinguishable (non-identical) particles as a conventional measure for the spatial entanglement [25]: while for N same-spin electrons the component of the entanglement that reflects the trivial minimum correlation due to the anti-symmetrization of the wave function can be eliminated [26][27][28], yielding: This definition ensures that for the wave functions used in HF approximation (or in general for any Slater rank 1 many-body state [29,30]) the linear entropy should vanish.
For indistinguishable (identical) particles, the spatial part of the 2N-body wave function can be represented in the simplest case as a product of normalized spin-up and spin-down Slater determinants [10]: Thus, the entanglement between-for example-spin-up only electrons can be estimated using the reduced density matrix, averaged over the configurations provided by the TDQMC algorithm: where D k↑ i are spin-up Slater determinants composed by the individual wave functions.

Results
As an example here we calculate the ground state of a quantum dot with parabolic core potential V en (r i ) = ω 2 r 2 i /2 and with soft-core electron-electron Coulomb repulsion [31]: where r ≡ r i − r j . Within the formalism of Section 2 the degree of spatial correlation and, hence, the spatial entanglement, is controlled in TDQMC by the quantum nonlocal length σ j,i where, for bound electrons, higher σ j,i lead to lower correlation (entanglement) between the ith and the jth electron, and vice versa. Since the spatial extent of the electron cloud for the jth electron is determined by the standard deviation s j of the corresponding MC ensemble, the nonlocal length σ j,i is expected to be close to s j : where α j,i may now serve as the variational parameters to minimizing the energy. Since for parallel-spin electrons the eigenstates are orthogonal to one another, their overlap is small and, hence, the dynamic correlation between such states is expected to be smaller compared to the correlation between opposite-spin electrons. Here, we consider 1D quantum dots in two basic configurations: one is a spin-polarized configuration where each energy level is filled with just one same-spin electron, as seen in Figure 1a, and the other is a spin-compensated configuration where each level is occupied by two opposite-spin electrons (Figure 1b). We start with the ground state and the three excited states of a total of four same-spin electron configurations (Figure 1a), where due to the orthogonality of the spatial wave functions the dynamic correlations between the same-spin electrons are expected to be rather small. The system of coupled non-linear TDQMC equations (Equation (8)) is solved iteratively using both spilt-step Fourier and Crank-Nicolson numerical schemes, which give close results. The final states obtained are statistical variants of the corresponding Hartree-Fock wave functions for the different energy levels, with certain distortions due to the different effective interaction potentials V k e f f (r i , t) in Equation (8). Starting from preliminary calculated Hartree-Fock wave functions, after 200 steps of imaginary-time propagation of Equations (8)- (11) and Equations (13) and (14), for ω = 1, and for different values of the variational parameter α j,i of Equation (22), we find the energy minima for one, two, three, and four occupied levels in succession. Figure 2a shows these energies (green line), which are in a very good correspondence with the numerically exact energies (blue line) obtained from the direct numerical solution of the Schrödinger equation for up to four electrons in one spatial dimension. Our calculations reveal that accuracy of three significant digits for the energy can be attained by varying α 1,i , while α 2,i , α 3,i , and α 4,i are set to infinity, which practically keeps the ground state (level 1) at its Hartree-Fock geometry. The optimal values of α 1,i for the ground level |1 are shown in Figure 2c for two, three, and four electrons, also showing that both α 1,i and the nonlocal length σ 1,i are almost independent of the number of electrons-except for one electron at the ground state where there is no e-e interaction-and α 1,i is set to zero. The degree of entanglement for the four configurations of Figure 1a is quantified by the linear quantum entropy of Equation (18), as shown in Figure 2b. There are two distinct cases: In the first, the electrons are considered identical (blue line), to be compared with the result from the exact numerical solution of the Schrödinger equation (green line). It can be seen that the linear entropy in this case remains almost constant, in close agreement with the exact numerical result. The second case, plotted with red dots in Figure 2b, depicts the linear entropy for the different electrons considered as distinguishable particles with the density matrix of Equation (17). It can be seen that the linear entropy for the distinguishable electrons increases due to the screening effect of the inner electrons, which causes larger fluctuations in the shape of the outer wave functions. For the filled-shell configuration of Figure 1b, each level contains two opposite-spin electrons, whose wave functions overlap in space almost completely, and therefore these electrons interact more like bosons [14]. It is therefore reasonable to calculate the only spatial entanglement that is due to same-shell states that is expected to significantly exceed the entanglement of the same-spin electrons at different levels. For the configuration of Figure 1b we have found that the major source of entanglement is the Coulomb repulsion between the two outermost electrons. For the ground state (level 1), where only two opposite-spin electrons are present in the vicinity of the core, the system energy exhibits a well-defined minimum as a function of the nonlocal length σ ↑ 1,1 = σ ↓ 1,1 , as seen in Figure 3a. When adding electrons at the higher levels, the corresponding wave functions acquire zeros, which is a source of larger fluctuations of the energy, as seen in Figure 3b-d, where the red curves represent the polynomial least squares fit for better visualization of the energy minimum.
Note that in the process of adding correlated electrons at the outer shells, the innershell electrons remain largely intact due to their stronger localization (confinement) to the core. Therefore, to a good approximation, the inner-shell wave functions need not be recalculated, and these may remain at their self-consistent Hartree-Fock configurations. Figure 4a depicts the system energy (blue line), which almost perfectly matches the numerically exact energies (green line) obtained using the standard DMC method [10]. The linear entropy predicted by the TDQMC method decreases when adding new excited states to the electron configuration (blue line in Figure 4b), which contrasts with the case of bosonic quantum dots, where it increases for more electrons at the ground state [14]. This behavior is also confirmed by the numerically exact results (green line) for levels 1 and 2 (two and four identical electrons), and can be explained by the orthogonality of the wave functions in the fermionic calculation, which causes weaker interaction for the outer-shell electrons.

Conclusions
In conclusion, we have calculated the correlated ground state and excited states of 1D quantum dots with spin-polarized and spin-compensated electron configurations with up to 5 energy levels (10 electrons) within the time-dependent quantum Monte Carlo framework. Unlike the Hartree-Fock approximation, the TDQMC method accounts for the dynamic correlations between the constituents of the quantum system, and allows one to quantify the resultant spatial entanglement in a simple and efficient manner. By variationally minimizing the system energy with respect to the spatial quantum nonlocality, the optimal set of wave functions that describes each electron is found, which allows one further to calculate the reduced density matrices for the different electrons and, hence, to quantify the entanglement that they exhibit due to their mutual interactions. Using the linear quantum entropy as a measure for the entanglement, it was found that for a fully spin-polarized electron configuration the stochastic windowing applied to the ground state alone is sufficient to recover the entanglement of the excited states, in good agreement with the exact numerical result. For the spin-compensated electron system, the entanglement that is due to the interaction of the two outermost opposite-spin electrons is dominant, while the inner shells remain largely at their Hartree-Fock states. An essential advantage of this method is that it allows one to conceive quantum particles as identical as well as distinguishable objects. The theory presented here may find useful applications in treating quantum correlation effects in composite quantum systems such as molecules, clusters, and solid-state materials. Its accuracy could be further improved by using linear combinations of Slater determinants to better approach the Hilbert space of the quantum many-body problem.