Fast nonadiabatic dynamics of many-body quantum systems.

Modeling many-body quantum systems with strong interactions is one of the core challenges of modern physics. A range of methods has been developed to approach this task, each with its own idiosyncrasies, approximations, and realm of applicability. However, there remain many problems that are intractable for existing methods. In particular, many approaches face a huge computational barrier when modeling large numbers of coupled electrons and ions at finite temperature. Here, we address this shortfall with a new approach to modeling many-body quantum systems. On the basis of the Bohmian trajectory formalism, our new method treats the full particle dynamics with a considerable increase in computational speed. As a result, we are able to perform large-scale simulations of coupled electron-ion systems without using the adiabatic Born-Oppenheimer approximation.

INTRODUCTION

Let us consider a many-particle electron-ion system at finite temperature. In calculating the dynamics of both the electrons and ions, we must account for the fact that the ions evolve multiple orders of magnitude more slowly than the electrons, as a result of their much higher masses. If we are interested in the long-time ionic dynamics (for example, the ion mode structure), then we face a choice of how to deal with this time scale issue. We can either model the system on the time scale of the electrons—nonadiabatically—and incur a substantial computational cost (a cost that is prohibitive in most simulation schemes), or model the system on the time scale of the ions—adiabatically—by treating the electrons as a static, instantaneously adjusting background. The latter approach is far cheaper computationally but does not allow for a viable description of the interplay of ion and electron dynamics.

The method we propose here (see Materials and Methods for details) enables us to use the former (nonadiabatic) approach, retaining the dynamic coupling between electrons and ions by reducing the simulation’s computational demands. We achieve this by treating the system dynamics with linearized Bohmian trajectories. Having numerical properties similar to those of molecular dynamics for classical particles (, ), our approach permits calculations previously out of reach: Systems containing thousands of particles can be modeled for long (ionic) time periods, so that dynamic ion modes can be calculated without discounting electron dynamics.

RESULTS AND DISCUSSION

To demonstrate the strength of our new method, we apply it to warm dense matter (WDM). With densities comparable to solids and temperatures of a few electron volts, WDM combines the need for quantum simulations of degenerate electrons with the description of a strongly interacting ion component. These requirements make WDM an ideal testbed for quantum simulations (). Further, as the matter in the mantle and core of large planets is in a WDM state (, ), and experiments toward inertial confinement fusion exhibit WDM states transiently on the path to ignition (, ), simulations of WDM are of crucial importance in modern applications.

Key dynamic properties of the WDM state can be represented by the dynamic structure factor (DSF) (). This quantity also connects theory and experiment: Probabilities for diffraction and inelastic scattering are directly proportional to the DSF (–), allowing testing of WDM models. Here, we focus on the ion-ion DSF that is defined viawhere N is the total number of ions and ρ(k, t) is the spatial Fourier transform of the ion density. In the following, we assume the WDM system to be isotropic and spatially uniform. Accordingly, the structure factor depends only on the magnitude of the wave number, k = ∣k∣. While the main contribution to S(k, ω) is due to direct Coulomb interactions between the ions, the modifications due to screening are strongly affected by quantum behavior in the electron component.

State-of-the-art calculations of the structure factor are typically carried out with variants of density functional theory (DFT) (). DFT’s Kohn-Sham formulation () has been the basis for many fundamental physical insights, and it has been successfully applied to fields as diverse as quantum chemistry, condensed matter, and dense plasmas (–). Recent work, however, has shown that predictions from standard DFT simulations for the ion-ion DSF are problematic due to the use of the Born-Oppenheimer approximation (). By using a Langevin thermostat alongside DFT, it has been found that the dynamics of the electron-ion interaction may strongly change the mode structure—in particular, the strength of the diffusive mode. However, this approach requires a very simple, uniform frequency dependence for electron-ion collisions, which may not prove realistic in practice. It also contains an arbitrary parameter, the Langevin friction, and is of limited predictive power as a result.

We demonstrate here that our new method of Bohmian dynamics, which retains the dynamics of the electron-ion interaction, can overcome the shortcomings of previous approaches without introducing free parameters. The specific case that we consider is compressed liquid aluminum with a density of 5.2 g cm−3 and a temperature of 3.5 eV, which allows for direct comparison with previous results. Full details of the corresponding simulations and input parameters are given in the Supplementary Materials.

The validity and accuracy of our implementation of Bohmian dynamics are strongly supported by the excellent reproduction of static ion-ion correlations from DFT simulations. Figure 1 illustrates the static ion-ion structure factor obtained with the Bohmian trajectory technique. This quantity is the Fourier transform of the pair distribution function and, thus, represents the degree of correlations present in the system (). The comparisons with orbital-free DFT and the computationally more intensive Kohn-Sham DFT both yield agreement within the statistical error of the simulations. This match was achieved by a single parameter fit defining λ (see Materials and Methods). The different simulation techniques predict almost the same thermodynamics as shown by the small pressure difference.

Fig. 1

Static ion-ion structure factors for aluminum.

The static structure factor is defined as S(k) = ∫ S(k, ω)dω. The main graph compares our results from Bohmian dynamics with data obtained by density functional theory molecular dynamics (DFT-MD) with orbital-free DFT (OFDFT) () for a density of 5.2 g cm−3 and a temperature of 3.5 eV. The lower insets compare our results to data from full Kohn-Sham DFT (KS-DFT) simulations at solid density and two different temperatures. The excellent agreement of the methods is also demonstrated by the very small differences in pressure as quantified by the parameter R: These values give the difference in ionic pressure between the methods normalized to the difference of the DFT pressures and the pressure of an ideal gas, that is, R = (PBohm − PDFT)/(PDFT − P0).

Figure 2A shows calculations of the fully frequency-dependent DSF. One can clearly notice the appearance of side peaks in the DSF that correspond to ion acoustic waves. Their dispersion for smaller wave numbers, and the corresponding sound speed, is very sensitive to the interactions present in the system. Thus, they reflect the screening of ions by electrons as well as dynamic electron-ion collisions. For larger wave numbers k, these modes cease to exist because of increased damping. The data also exhibit a diffusive mode around ω = 0, although it is not as prominent as predicted in ().

Fig. 2

Results for the dynamic ion structure for aluminum at 3.5 eV and 5.2 g cm−3.

(A) The frequency-resolved DSF from the Bohmian dynamics. (B) Comparison of the dispersion relation of the ion acoustic modes from our Bohmian approach with the data from the Langevin model of ().

The dispersion relation of the ion acoustic modes is displayed in Fig. 2B, which also contains results from the Langevin approach. The latter approach requires ad hoc friction parameters that were chosen to cover the range between the classical and quantum Born limits [see ()]. The strength of the Bohmian approach lies in the fact that it does not require a free parameter, thereby allowing access to a self-consistent prediction of the sound speed. This comparison may also be used to assess the quality of the friction parameter applied in the Langevin approach. For the case considered, one finds that neither the classical nor the weak coupling Born limit is applicable—a finding that is typical of the WDM regime.

Figure 2A demonstrate the strength of the Bohmian approach in modeling quantum systems with strong interactions and nonlinear ion dynamics. For static and thermodynamic properties, we obtain results in very close agreement with DFT simulations. In addition, while the standard implementation of DFT involves the Born-Oppenheimer approximation, our Bohm approach can treat electrons and ions nonadiabatically, retaining the full coupling of the electron and ion dynamics. As a result, we can investigate the changes of the ion modes due to dynamic electron-ion correlations that are inaccessible to standard DFT. In contrast to a Langevin model, we have no free parameters and can thus predict the strength of the electron drag to the ion motion. Simulations based on time-dependent DFT () represent another way to avoid the Born-Oppenheimer approximation. However, this method is numerically extremely expensive, drastically limiting particle numbers and simulation times; at present, this limitation precludes results for the ion modes as presented here.

The principal advantage of our approach is its relative numerical speed, which allows for the modeling of quantum systems with large numbers of particles. For comparison, the recent time-dependent DFT simulation of () models a system of 128 electrons for approximately 0.001 as per central processing unit (CPU) core and second. The comparative Bohmian dynamics system models eight times as many electrons for approximately 20 as per CPU-core and second. This vast difference in computation time enables our method to access a new class of correlated quantum systems. These calculations are not only relevant for WDM but also address core problems in chemical and biological systems (e.g., protein folding), as well as radiation damage of materials (–).

While our initial implementation of Bohmian dynamics focuses on establishing dynamic correlations of systems in thermal equilibrium, generalization to nonequilibrium systems is also possible through dynamically updating the system potentials to account for local, time-dependent thermodynamic conditions. In particular, the electron-ion or electron-phonon energy exchange in two-temperature systems is amenable to this approach.

MATERIALS AND METHODS

To begin constructing our method, we consider Bohm’s formulation of quantum mechanics. We can imagine a classical N-body system as an “N-trajectory” moving through 3N-dimensional configuration space. Bohm’s theory treats an N-body quantum system as an ensemble of these classical N-trajectories, interacting through an additional N-body potential VB. This Bohm potential is a functional of the density of N-trajectories in configuration space, Φ, and a function of the spatial position x, that is, VB = VB(x∣Φ). Provided matching initial conditions, the Bohm ensemble of N-trajectories reproduces the dynamics of the probability density—as given by the Schrödinger equation—exactly () (see also the Supplementary Materials for more details).

In its exact form, Bohm’s formulation is as intractable as the N-body Schrödinger equation, requiring simulation of a huge number of N-body interacting classical systems. However, we can construct a fast computational method solving Bohm’s theory by introducing a thermally averaged, linearized Bohm potential. The exact (but inaccessible) calculation for a pure quantum state with many particles—based on the theory above—would require us to propagate an array of N-trajectories through time, at each step recalculating their density and, thereby, VB (see Fig. 3A). Reliable calculations of density in 3N-dimensional space would require a prohibitive number of N-trajectories, however, making this unfeasible. Here, we propose an alternative: As opposed to applying this theory to a pure quantum state, we consider propagating an array of thermally coupled N-trajectories in a similar manner, as a model of a system at finite temperature. Our core assumption is that the time evolution of a finite-temperature quantum system can be approximated by a similar procedure to that used for the pure state; we consider an array of N-trajectories, each coupled to a heat bath setting its temperature, evolving under a potential that is a functional of N-trajectory density in configuration space. This procedure can be seen as a trajectory-based analog of the linearization of the Bohm potential over states in quantum hydrodynamics ().

Fig. 3

Schematic of the applied linearization approximation.

(A) The time evolution of an N-trajectory in an exact Bohmian representation of a pure quantum state (top) and the Bohm potential, VB, that it experiences (bottom). VB isa functional of the density of N-trajectories in configuration space, Φ. At each time step, all of the N-trajectories in the ensemble must be updated, Φ calculated, and the updated Bohm potential determined. (B) The time evolution of an N-trajectory in our linearized Bohmian representation of a thermal state (top) and the Bohm potential that it experiences (bottom). We need only to track a single N-trajectory: Its coupling to a heat bath ensures its ergodicity, so that Φ becomes equal, in equilibrium, to its time-averaged density in configuration space. As a result, the N-trajectory evolves with a time-independent Bohm potential, generated self-consistently by its own time-integrated density. (C) The block panel summarizes our scheme to determine the Bohmian dynamics and gives the sources for the different potentials needed as input.

For simplicity, we focus initially on systems in thermal equilibrium. We allow our thermal N-trajectories—together modeling the probability density of our finite-temperature system—to evolve in time under a linearized mean Bohm potential of the underlying pure states. In this linear approximation, we replace the mean Bohm potential experienced, expressible as a sum over functionals of individual states, with a functional of a sum over individual states. In this way, we construct an effective configuration space potential that is an approximate average over the corresponding potentials of the exact N-body wave functions. This averaging scheme then acts as a direct estimate of the full thermal system.

In addition to moving focus from a pure state to a more practically important finite-temperature state, the key feature of our approximation is that it markedly reduces the computational expense required to simulate the system (as compared with the exact pure state case). Crucially, now that we are considering finite-temperature N-trajectories, we need only to track a single N-trajectory through time, rather than an impractical number of them. This follows from two observations:

1) Properties of a classical system with correlation-dependent potentials can be determined self-consistently. For an arbitrary system of N well-localized particles, the N-particle correlation function can be written as g(x) = P(x)/P0(x), where P denotes the joint positional probability distribution of the particles and P0 is the distribution for a noninteracting classical system with equal particle densities. Here, the variable x is the set of particle positions, x = {x1, x2, …, x}. We may construct interparticle potentials that are functionals of g: V(x) = V0(x) + V(x∣g), where V0 denotes pair interactions and external forces and V is a contribution that varies with g. The equilibrium properties of a system in this potential can be found self-consistently. Starting from an initial guess for V(), we can calculate g() through a Monte Carlo simulation or a similar method (). This value of the N-particle correlation function g then gives rise to a new approximation for the potential. Iterating this procedure allows for both g and V to be found.

2) Our linearized Bohmian system is equivalent to a classical system with correlation-dependent potentials. Consider a number of coupled thermal N-trajectories representing our linearized Bohmian system. Assuming that the system is in a temperature regime in which the particle motion is ergodic, we find that each N-trajectory has the same time-integrated correlations, that is, each has the same g. In the limit of infinite N-trajectories in our ensemble, it follows that the configuration space density of N-trajectories Φ is exactly proportional to this common g. As a result, each N-trajectory moves in a common static potential (see Fig. 3B), and, as this static potential is a functional of configuration space density, Φ, it is equivalently just a functional of g.

When combining the results above, our linearization approximation becomes a simple mappingwhere m is the mass of particle i and λ is a linearization factor that still needs to be determined. As g is common to all N-trajectories, our approximation scheme allows us to consider just a single N-body classical system (Fig. 3B). The required simulation is thus amenable to (computationally cheap) classical molecular dynamics (). The classical particle trajectories simulated then approximate the statistics of the full quantum system.

Before the scheme above can be implemented, we must overcome a final fundamental hurdle: The full correlation function g appearing in Eq. 2 is too complicated to be modeled directly (similar to the full N-body wave function). Therefore, we seek an approximate closure for this object in terms of lower-order correlations, which can be calculated accurately. For this goal, we use the pair product approximation, whereby the N-body correlation is replaced by a product of pair correlations. Furthermore, we generalize the dependence on λ to a set of λ to accommodate different particle species in the pair interactions.

The use of the pair product may appear to restrict our method to weakly coupled systems; however, the corresponding closure enters only into the calculation of the Bohm potential functional, rather than as a global restriction on the correlations treated. The pair correlation functions themselves are calculated with the full N-body system at each step of the algorithm, and the hierarchy of the N-body correlation effects is implicitly taken into account. We expect the pair product closure to begin to break down only when the system properties deviate strongly from those of a simple liquid—in such cases, a higher-order correlation closure should be constructed.

We need an additional correction to our potentials to fulfill the spin statistics theorem, as the Schrödinger equation—and thus Bohmian mechanics—does not incorporate particle spin directly. In particular, this correction will generate a Fermi distribution for the electrons in thermal equilibrium. Similar to successful approaches applied in quantum hydrodynamics and classical map methods (, , ), we introduce an additional Pauli potential term. This term is constructed such that exchange effects are reproduced exactly for a reference electron gas system. We also use pseudopotentials, commonly applied in modern DFT calculations, to represent core electrons bound in deep shells of the ions by an effective ion potential seen by the valence electrons.

Last, it remains to set the linearization parameters λ to fully determine the system’s Hamiltonian. In this work, we take λ = 1 for the ion-ion and ion-electron terms. To determine the electron-electron parameter λee, we match static ion correlations—that is, pair distribution functions—obtained by DFT calculations. This matching can be carried out rigorously with a generalized form of inverse Monte Carlo (see the Supplementary Materials). In this way, we determine λee with only static information of the system. Subsequent dynamic simulations can then be carried out without any free parameters.

We can now implement the Bohmian dynamics method with a molecular dynamics simulation inclusive of the potentials discussed above. Within the microcanonical ensemble, we could simply integrate the equation of motion for the particles. However, we want to simulate systems with a given temperature which requires the coupling of the system to a heat bath. Here, we used a modified version of the Nosé-Hoover thermostat. Its standard form is popular in classical molecular dynamics studies, achieving reliable thermodynamic properties while having minimal impact on the dynamics (). We use this standard version to control the ion temperature. The electrons, however, should relax to a Fermi distribution, which is not possible with this classical form. We have, thus, produced a modified version of the Nosé-Hoover thermostat for the dynamic electrons. It creates an equilibrium distribution equal to that of a noninteracting electron gas of the same density (see the Supplementary Materials for details and derivation). With this, we ensure the ions interact with an electron subsystem with a corrected energy distribution.