1. INTRODUCTION

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 [1]. 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. [4]), 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 [7]. 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 $N$-electron atom in an external field:

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 functions

Then, the windowing kernel transforms the probability density under the integrals in Eq. (5) into a mixture of densities, according to the rule

The next key assumption in TDQMC is that, in fact, the trajectories ${\mathbf{r}}_j^k(t)$ may be represented by a set of $M$ Monte Carlo walkers in physical space, each of which samples its own probability density given by ${|{\varphi }_j^k({\mathbf{r}}_j,t)|}^2$, where ${\varphi }_j^k({\mathbf{r}}_j,t)$ 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., [10]):

Since the windowing parameter ${\sigma }_j^k$ in Eq. (10) determines the amount of walkers belonging to the $j$th electron, which contribute to the electron–electron potential experienced by the $i$th electron, the parameter ${\sigma }_j^k$ is called the nonlocal quantum correlation length in TDQMC. Clearly, in the limiting case where ${\sigma }_j^k\to 0$, 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 $V_{e-e}^{\text{eff}}({\mathbf{r}}_i,t)=V_{e-e}[{\mathbf{r}}_i,{\mathbf{r}}_j^k(t)]$, which basically overestimates the electron–electron repulsion. The opposite case, where ${\sigma }_j^k\to \infty$, 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 ${\sigma }_j^k$ 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 ${\sigma }_j(t)$] 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 write

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 [14]:

 

Fig. 1. Schematic of the transition from the mean-field approximation where all Monte Carlo walkers are guided by the same wave function to the TDQMC representation where each walker is guided by its own guide wave.

Download Full Size | PDF

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 density

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 $\mathbf{\eta }(t)$ 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 $\mathbf{\eta }(t)=(\tau -t){\mathbf{\eta }}_0$ is well suited, where ${\mathbf{\eta }}_0$ is a random vector with zero mean and unit variance and $\tau$ 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 [5]. 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 ${\alpha }_j$ of Eq. (13) through minimization of the ground state energy. Once the value of ${\alpha }_j$ is determined, the branching can be turned on in the standard way [5], 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., [16]; atomic units are used henceforth). We have

The parameters $a$ and $b$ determine the strength of the electron–nuclear and electron–electron interaction, respectively. For helium in the $s$-state (antiparallel spins), we have $a=2$ and $b=1$ in Eqs. (21) and (22). Then, the set of TDQMC equations [Eqs. (9)–(11)] reads as

We start here with the calculation of the ground state of the atom by minimizing the total energy $E=\sum _k^M{E}_{\text{loc}}^k$ by varying the parameters ${\alpha }_j$; $j=1$, 2. Since the two electrons are in the $s$-state, we have ${\alpha }_1={\alpha }_2=\alpha$. 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 $\alpha$ 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 $-2.2391\text{a.u.}$., which is close to the exact numerical value of $-2.2393\text{a.u.}$. Next, Fig. 3(a) shows a snapshot of the walker distribution at the end of the ground state calculation for $\alpha =0.9$, 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 $x_1=x_2$ 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.

 

Fig. 2. Energy of the ground state of a 1D helium atom as function of the parameter $\alpha$ of Eq. (13), with corresponding error bars (for 20 consecutive runs). The rightmost point shows the transition to the mean-field ground state energy.

Download Full Size | PDF

 

Fig. 3. Snapshot of the walker distribution for (a) the ground state and (b) the corresponding smoothed distribution (red line). The blue line in (b) depicts the exact ground state.

Download Full Size | PDF

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 $V_{e-n}$ is suddenly set to zero after the ground state preparation has ended; see Fig. 4. We have

 

Fig. 4. Contour maps of (a) the moduli of the density matrices for the ground state and (b) after free diffraction of the interacting electron waves. Exact result—blue contours; TDQMC result—red contours.

Download Full Size | PDF

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 $V_{\text{ext}}(x_i,t)=-x_{i}E(t)$, where $E(t)$ 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.

 

Fig. 5. (a) Survival probability for a single electron in a one-dimensional helium atom exposed to a powerful ultrashort laser pulse, and (b) the quantum coherence calculated as sum of the antidiagonal elements of the density matrix in coordinate representation. Red lines—from TDQMC; blue lines—exact result.

Download Full Size | PDF

4. CONCLUSIONS

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) [17].