Here, a dichotomy of particles and waves is employed in a quantum Monte Carlo calculation of interacting electrons. Through the creation and propagation of concurrent stochastic ensembles of walkers in physical space and in Hilbert space, one can correctly predict the ground state and the real-time evolution of a single electron interacting with a larger quantum system. It is shown that such walker ensembles can be constructed straightforwardly through a stochastic sampling (windowing) applied to the mean-field approximation. Our calculations reveal that the ground state and the real-time evolution of the probability distributions and the decoherence due to the Coulomb interaction in the presence of strong ultrashort laser pulses can be accounted for correctly by calculating the density matrix of the electron, without referencing the quantum many-body state of the whole system.
© 2017 Optical Society of America
The recent experimental progress of laser technologies has triggered substantial research in light–matter interactions on ultrafast time scales with theoretical efforts directed toward a better understanding of the intimate mechanisms that control the motions of electrons in atoms, molecules, and condensed matter . A significant part of that research focuses on ab initio methods, where in the ideal case, the direct (numerical) solution of the time-dependent Schrödinger equation would provide not only the ground state but also the real-time evolution of the many-body wave function, which next would allow all quantities of interest to be calculated. Unfortunately, the wave function is a physical object only for an isolated particle, while in the many-body case, it is defined in configuration space, which introduces both formal and computational difficulties in the quantum mechanics framework in that the solution of the many-body quantum problem by a classical computer requires effort that scales exponentially with the system size. Therefore, substantial progress in both methodology and numerical algorithms would be needed to deliver accurate and numerically efficient methods for real-time evolution of many-body quantum systems.
While obvious for the many-body Schrödinger equation, the exponential scaling may also concern the particle-based methods (e.g., quantum Monte Carlo) and, in general, any method pursuing an exact solution [2,3]. Therefore, it is expected that much effort has been devoted to reformulate quantum mechanics in terms of physical space (one-body) wave functions, which are easier to calculate under different conditions as required by the experiment. The most prominent of those are the Hartree–Fock and the density functional approximations. The reduction of the many-body quantum problem to lower dimensionality comes, however, at the price of losing essential quantum properties such as the quantum correlations due to the particle interactions, which are simply neglected in Hartree–Fock and introduced semi-empirically in the density functional. This leads, for example, to inaccurate estimations of the ground state and/or binding energies of some atoms and molecules. One class of improved accuracy many-body methods involves multiple configurations (e.g., in Ref. ), which can be considered to be a generalization of Hartree–Fock, where the quantum correlations can be accounted for. However, these methods are significantly less numerically efficient for realistic systems, which may involve multiple electrons. The situation is somewhat improved with the help of quantum Monte Carlo (QMC) methods (e.g., the diffusion QMC, [5,6]), where variationally optimized but essentially artificial many-body trial functions of Slater–Jastrow-type are used to guide the diffusing Monte Carlo walkers to their most energetically favorable positions in the presence of interaction potentials. In this way, very accurate estimates for the ground and excited state energies have been obtained, but the real-time evolution still remains beyond reach. Besides that, for convergence of the results for fermions, additional constraints have to be imposed close to the nodes of the ground state wave function (fixed node approximation), which are aimed to deal with the infamous “sign problem.” Somewhat complementary to the diffusion QMC is the auxiliary field QMC (AFQMC), where the random walk occurs in the state (or wave-function) space where an imaginary-time propagator projects the trial state to the ground state . At the heart of AFQMC lies the continual Fourier transform (Hubbard–Stratonovich transformation), which recasts the system of interacting particles to a system of noninteracting particles moving in stochastic fields. However, for repulsive interactions (e.g., for those between each two electrons), the stochastic potentials turn out to be complex quantities, which is counterintuitive for the ground state and causes a “phase problem” analogous to the “sign problem” above. The essentially stochastic nature of AFQMC also prevents its use for real-time propagation. In an attempt to address in a numerically efficient way some of the existing problems in the treatment of interacting quantum particles for both imaginary- and real-time propagation, a new quantum Monte Carlo method was introduced recently [8,9]. The method is called time-dependent quantum Monte Carlo (TDQMC) and it involves concurrently ensembles of particles and wave functions, such that these evolve self-consistently in physical spacetime. The symmetric way in which particles and waves are treated in TDQMC is called here particle–wave dichotomy. The method is well suited for calculations of both ground state properties and real-time evolution of quantum averages and coherences for many-body quantum systems.
This paper provides a systematic derivation of the working equations of the TDQMC method and presents numerical results for the ground state and the time evolution of a simple quantum system. The results for a one-dimensional helium atom in an external field are compared with those from the exact numerical solution of the many-body Schrödinger equation.
2. UNLOCKING THE QUANTUM CORRELATIONS FROM THE MEAN-FIELD APPROXIMATION
A. Time-Dependent Quantum Monte Carlo
We start here with the many-body Schrödinger equation for an -electron atom in an external field:1) is a sum of electron–nuclear, electron–electron, and external potentials:
From here, we use as a starting point the mean-field approximation, which can be derived either variationally or by substituting a simple representation of the many-body wave function expressed as a product of one-body wave functions1) and integrating over the spatial variables (except ) 4). The fact that we assign a single wave function to the th electron implies that we presume that the electron be in a pure state. It is clear, however, that going beyond the mean-field approximation, the interaction of each electron with the rest of the electrons should, in principle, lead to a mixture of states due to the correlated electron motion, which remains hidden in the mean field because the Hartree potential depends on the entire distribution for the th electron in Eq. (5). One simple intuitive way to unlock the electron–electron correlations from the mean-field approximation would be to introduce randomness in the local environment experienced by each electron. Accordingly, it is assumed in TDQMC that portions of the th electron cloud may experience the Coulomb potential due to random parts of the th electron cloud and vice versa. One way to accomplish this is to introduce a stochastic windowing of in Eq. (5) by a kernel function with a width , centered at a point in space, such that the product is different for each separate . For example, for Gaussian kernel we have
Then, the windowing kernel transforms the probability density under the integrals in Eq. (5) into a mixture of densities, according to the rule7) defines the “windowing,” where is the “window” in space centered at and weighted by . Then, substituting Eq. (7) into Eq. (5) introduces a new (effective) potential: 4), into a set of equations for the different replicas of the original one-body wave function :
The next key assumption in TDQMC is that, in fact, the trajectories may be represented by a set of Monte Carlo walkers in physical space, each of which samples its own probability density given by , where represents a concurrent walker in Hilbert space. This allows us to transform the integrals in Eq. (8) into sums according to the Monte Carlo integration rule, thus allowing us to express the effective potential as a Monte Carlo convolution:
Finally, the TDQMC equations for the wave functions, Eqs. (9)–(11), need to be complemented by equations that connect the motion for the Monte Carlo walkers in physical space with those in Hilbert space, where the most natural choice is the de Broglie–Bohm guiding equations (see, e.g., ):12), the walkers in Hilbert space will be called, henceforth, “guide waves.” In this way, Eqs. (9) through (12) present a method to evolve concurrent ensembles of particles and guide waves in a self-consistent manner in space and time.
Since the windowing parameter in Eq. (10) determines the amount of walkers belonging to the th electron, which contribute to the electron–electron potential experienced by the th electron, the parameter is called the nonlocal quantum correlation length in TDQMC. Clearly, in the limiting case where , the set of TDQMC equations for the guide waves [Eq. (9)] is transformed into a set of linear Schrödinger equations with a classical pairwise potential , which basically overestimates the electron–electron repulsion. The opposite case, where , returns us to the mean-field approximation, which underestimates the electron–electron repulsion. While the results are not expected to be sensitive to the exact shape of the kernel in Eq. (10), for, e.g., bell-shaped kernels, it is physically justified to specify the parameter to be proportional to some global statistical parameter that characterizes the coordinate uncertainty for each separate quantum particle, such as the kernel density estimation bandwidth [or simply the standard deviation ] of the corresponding Monte Carlo sample, and to be inversely proportional to the average distance between the walkers for the different electrons. In case of localized distributions for most of the time (e.g., in an atom at ground state), one can write5] but with less variational parameters. For practical implementation, it should be noted that a Monte Carlo sampling occurs during the preparation of the ground state of the system through a Markovian diffusion process, consistent with the kinetic part of Eq. (9) (see also ): 12) and is a vector random variable with zero mean and time-dependent variance that tends to zero toward steady state. Clearly, for imaginary-time propagation, the velocity in Eq. (14) vanishes from Eq. (12). Note that the random motion of the walkers gives rise to randomness of the kernel windowing function of Eq. (10), which may thus be considered as a stochastic auxiliary field to participate in the one-body Schrödinger equation for the guide waves [Eq. (9)]. However, unlike in AFQMC, the resulting potential in Eq. (10) stays real for all times and it is stabilized toward the end of the ground state propagation stage, which is of crucial importance for the stability of the whole method. Also, an important sampling that involves branching (birth and death of couples of Monte Carlo particles and waves) was found to further improve the convergence and the accuracy of the TDQMC calculation during the preparation stage, similar to in other quantum Monte Carlo methods. Note that in case of branching, an additional term should appear in the RHS of Eq. (14) to reflect the “quantum force” . Once steady state is established, the walkers are next guided by the corresponding guide waves [see Eq. (12)] for real-time propagation. This contrasts other approaches that use stochastic wave functions, e.g., in the quantum jump method [11,12], where a significant stochastic component is introduced also during the real-time propagation of the quantum system, which limits the propagation time to the so-called “useful time” .
B. Statistics of Particles and Waves
We have seen so far that the quantum particle–wave dichotomy can be reformulated for interacting particles where instead of stating that each physical particle is guided by a guide wave we can say that each physical particle is represented by a set of fictitious particles (walkers) whose distribution in physical space reproduces its probability density, and each of those walkers is guided by its own wave function (guide wave), as drawn schematically in Fig. 1. Then, as indicated in Fig. 1, in the mean-field approximation, all walkers are guided by the same wave function and the electron is in a pure state. In other words, together with an ensemble of walkers in physical space, for each physical particle there exists a joined statistical ensemble of guide waves, which differ from each other if that particle is in a mixed state. This, on the other hand, warrants the introduction of appropriate density matrices and correlation functions based on these guide waves. The definition of the density matrix in our case is simplified by the fact that the classical probability for occurrence of a given wave function (guide wave) from the mixture is imposed by the distribution of the Monte Carlo walkers in space and therefore we may construct the density matrix for each separate electron directly, e.g., by using a simple mean over the guide waves :12], Chap. 3). It is clear, therefore, that the one-body density matrix of Eq. (15) may provide an instrumental approach to access probability densities and coherences for parts of a complex quantum system without having to access the many-body quantum state. For example, averages such as the dipole moment or the atomic ionization of an electron exposed to an external electric field can be calculated using either the walkers in physical space or those in Hilbert space (the guide waves) provided by the TDQMC methodology.
Besides for the calculation of densities (distributions), the walker trajectories can be conveniently used to calculate various many-body effects without having to calculate multidimensional integrals. For example, starting from the two-electron density15] 3), gives after a simple integration (for spherical symmetry)
3. NUMERICAL RESULTS
Before we present the results of the numerical implementation of TDQMC for simple systems, a few remarks are needed. As already mentioned, Eq. (14) accomplishes the diffusion of the walkers in physical space where each walker samples the modulus square of its own guide wave and where the walkers velocity tend to zero at the ground state. It is important that the time dependence of the stochastic variable is not too abrupt because the walkers would not manage to thermalize to the correct ground state distribution. It was found that a linear dependence is well suited, where is a random vector with zero mean and unit variance and is the moment where steady state is established. Also, besides the diffusion, most of the quantum Monte Carlo methods implement branching of walkers, which reflects the influence of the potential term and is based on the estimate of the local energy that rules whether a given move is acceptable . Although the branching improves the convergence of the walk, we found that at the variational stage, the TDQMC algorithm operates better without using it for more precise evaluation of the parameter of Eq. (13) through minimization of the ground state energy. Once the value of is determined, the branching can be turned on in the standard way , except that here we have birth and death of both particles and guide waves based on the local energy given by
In order to compare the TDQMC results with the numerically exact solution of the many-body Schrödinger equation, here we employ a simplified model for a one-dimensional helium atom with a soft core potential, which has proven to be one of the most strongly correlated systems known (e.g., ; atomic units are used henceforth). We have
The parameters and determine the strength of the electron–nuclear and electron–electron interaction, respectively. For helium in the -state (antiparallel spins), we have and in Eqs. (21) and (22). Then, the set of TDQMC equations [Eqs. (9)–(11)] reads as10)]
We start here with the calculation of the ground state of the atom by minimizing the total energy by varying the parameters ; , 2. Since the two electrons are in the -state, we have . The size of the grid for this calculation is 30 a.u., and the imaginary time step is 0.1 a.u., where 200 time steps are sufficient to reach a stable ground state. Figure 2 depicts the minimization curve, which indicates that the optimal value of is indeed close to unity, as expected for such a simple quantum system where the spatial uncertainty determines the range of spatial quantum nonlocality. The ground state energy we found is ., which is close to the exact numerical value of . Next, Fig. 3(a) shows a snapshot of the walker distribution at the end of the ground state calculation for , while Fig. 3(b) shows a distribution smoothed by a kernel density estimation (red line) compared with the exact solution obtained from Eqs. (20)–(22) (blue line). It is seen from Fig. 3(b) that the two distributions match quite satisfactorily and that the TDQMC calculation correctly reflects the electron repulsion along the axis [the main diagonal in Fig. 3(b)]. The total number of walkers and guide waves for this calculation is 64,800, distributed among 216 compute processes, where the ground state calculation takes less than a minute.
In order to further test the accuracy of TDQMC for real-time propagation, we carried out a calculation of the evolution of the one-body density matrix of a single electron [see Eq. (15)] for ultrafast electron scattering where the nuclear potential is suddenly set to zero after the ground state preparation has ended; see Fig. 4. We have4 also the numerically exact result (blue line) obtained from the direct integration of the two-body time-dependent Schrödinger equation [Eq. (20)] and next use the standard definition of the reduced (single-electron) density matrix expressed through the resulted wave functions:
It is seen from Figs. 4(a) and 4(b) that the TDQMC result matches almost perfectly with the exact one, which indicates that not only the distributions [Fig. 3(b)] but also the coherences are calculated with very good accuracy by the TDQMC algorithm, without referring explicitly to the many-body quantum state. In this calculation, the space grid spans 50 a.u. and a total of 43,000 walkers and guide waves take part while the execution time is somewhat longer during the real-time propagation due to the calculation of two-dimensional quantities, such as in Eq. (26), which are not an essential part of the TDQMC algorithm.
The last numerical result here deals with the ionization and the degree of quantum coherence of one-dimensional helium exposed to a strong few-cycle laser pulse, using the diagonal and the off-diagonal elements of the density matrix, Eq. (26), respectively. The external potential in dipole approximation reads as , where is the incoming pulse with amplitude 0.16 a.u., carrier frequency 0.1 a.u., and duration of two periods of the carrier. Physically, during each half-cycle the electric field pushes the electron cloud away from the core, which corresponds to the steps in the ionization curves in Fig. 5(a), while at the same time the electron loses coherence due to its scattering from the core potential and the electron–electron interaction [Fig. 5(b)]. Here the degree of coherence is quantified as a sum of the elements along the main antidiagonal of the density matrix. It is seen from Figs. 5(a) and 5(b) that both the ionization and the degree of coherence stay in close correspondence with the exact result (blue lines) calculated from Eqs. (20)–(22). This indicates that, despite the huge number of essentially stochastic wave functions involved in the TDQMC calculation, the elements of the resulting one-body density matrix of Eq. (15) are calculated with a very good accuracy.
It is shown in this paper that concurrent ensembles of quantum trajectories in physical and Hilbert spaces can be introduced within the time-dependent quantum Monte Carlo method such that the predicted spatial distributions of walkers and one-body density matrices match well with the exact solution of the many-body Schrödinger equation. Since thousands of replicas of the guide waves are generated and propagated under random effective potentials, the fully correlated picture is recovered in TDQMC, without the need to introduce multiple configurations in advance as in other methods. The calculated density matrix correctly predicts the decoherence resulting from electron scattering and due to the ionization of an electron in a two-electron atom by a powerful external field. This implies that by creating sufficiently rich statistics of walkers in physical space and in Hilbert space for each quantum species, one can facilitate building an efficient approach to solving many-body time-dependent quantum problems. Our calculations reveal that although TDQMC is essentially a stochastic method, it can perform well for both imaginary-time and real-time propagations as compared to other (imaginary-time) methods using stochastic auxiliary fields. The stochasticity of the electron–electron potential is introduced through the random motion of the Monte Carlo walkers in TDQMC, which makes it milder as compared to in auxiliary field quantum Monte Carlo, which greatly improves the stability of the time propagation, also allowing real-time evolution. The TDQMC algorithm offers an affordable time scaling that is almost linear with the number of electrons, which is decisive for larger systems. This is owing to the intrinsic parallelism where the computational load is distributed almost independently between the compute processes where the only communication needed is to calculate the nonlocal quantum correlations. Since the mass of the nucleus is much bigger than the electron mass, here we have considered the former to be a classical particle and thus we have neglected the electron–nuclear quantum correlations which, however, can be accounted for readily by including the nuclear wave functions in the product of Eq. (3) .
Bulgarian National Science Fund (BNSF) (FNI T02/10).
1. T. Brabec, ed., Strong Field Laser Physics (Springer, 2008).
2. M. Troyer and U. Wiese, “Computational complexity and fundamental limitations to fermionic quantum Monte Carlo simulations,” Phys. Rev. Lett. 94, 170201 (2005). [CrossRef]
3. A. Montina, “Exponential complexity and ontological theories of quantum mechanics,” Phys. Rev. A 77, 022104 (2008). [CrossRef]
4. H. D. Meyer, U. Manthe, and L. S. Cederbaum, “The multi-configurational time-dependent Hartree approach,” Chem. Phys. Lett. 165, 73–78 (1990). [CrossRef]
5. B. Hammond, W. Lester, and P. Reynolds, Monte Carlo Methods in Ab Initio Quantum Chemistry (World Scientific, 1994).
6. B. M. Austin, D. Y. Zubarev, and W. A. Lester, “Quantum Monte Carlo and related approaches,” Chem. Rev. 112, 263–288 (2012). [CrossRef]
7. G. Sugiyama and S. E. Koonin, “Auxiliary field Monte-Carlo for quantum many-body ground states,” Ann. Phys. 168, 1–26 (1986). [CrossRef]
8. I. P. Christov, “Correlated non-perturbative electron dynamics with quantum trajectories,” Opt. Express 14, 6906–6911 (2006). [CrossRef]
9. I. P. Christov, “Dynamic correlations with time-dependent quantum Monte Carlo,” J. Chem. Phys. 128, 244106 (2008). [CrossRef]
10. P. R. Holland, The Quantum Theory of Motion (Cambridge University, 1993).
11. K. Molmer, Y. Castin, and J. Dalibard, “Monte Carlo wave-function method in quantum optics,” J. Opt. Soc. Am. 10, 524–538 (1993). [CrossRef]
12. H.-P. Breuer and F. Petruccione, The Theory of Open Quantum Systems (Oxford University, 2002).
13. P. Deuar and P. Drummond, “First-principles quantum dynamics in interacting Bose gases: I. The positive P representation,” J. Phys. A 39, 1163–1168 (2006). [CrossRef]
14. I. P. Christov, “Double-slit interference with charged particles: density matrices and decoherence from time-dependent quantum Monte Carlo,” Phys. Scr. 91, 015402 (2015). [CrossRef]
15. R. P. Feynman, Statistical Mechanics (Benjamin, 1972).
16. R. Grobe and J. H. Eberly, “Photoelectron spectra for two-electron system in a strong laser field,” Phys. Rev. Lett. 68, 2905–2908 (1992). [CrossRef]
17. I. P. Christov, “Molecular dynamics with time dependent quantum Monte Carlo,” J. Chem. Phys. 129, 214107 (2008). [CrossRef]