eFF, a method to simulate large scale excited electron dynamics

Ground versus excited state wavefunctions

eFF computes the energy of a collection of point charge nuclei and spherical Gaussian electrons ( $ \psi = \exp((\mathbf{r} - \mathbf{x_{i}})^{2}/s_{i}^{2})$) as a function of the nuclear and electron positions plus the electron sizes. To obtain ground states and adiabatic dynamics, the electron positions and sizes are varied as a parametric function of the nuclear positions to minimize the overall energy of the system.

To obtain excited state dynamics, electrons are represented as wave packets:

$\displaystyle \Psi\left(\mathbf{r}, t\right) \propto \prod_{j} \,\exp\left[-\le... ...(\mathbf{r}-\mathbf{x})^{2}\right]\cdot\exp[i \mathbf{p_{x}} \cdot \mathbf{r}].$    
which substituted into the time-dependent Schrodinger equation and assuming a locally harmonic potential produces the semiclassical equations of motion [2] (Figure 4):
    $\displaystyle d\mathbf{{p}_{x}}/dt = -\nabla_{x}V, \ \mathbf{p_{x}} = m_{\mathrm{elec}} \dot{\mathbf{x}}$  
    $\displaystyle d{p}_{s}/dt = -\partial V/\partial s, \ p_{s} = (3 m_{\mathrm{elec}}/4) \dot{s}$  
Thus the average position of a wave packet obeys classical dynamics (Ehrenfest's theorem) with the addition that the size of the wave packet obeys classical dynamics as well. The radial effective mass factor $ 3 m_{\mathrm{elec}} / 4$ depends on the dimension of the wave packet, and becomes $ 2 m_{\mathrm{elec}}/4$ for a 2D Gaussian and $ m_{\mathrm{elec}}/4$ for a 1D Gaussian. Nuclei move via Ehrenfest dynamics [17,18].

Figure 4: Gaussian wave packet traveling in a one-dimensional harmonic potential; in this case, the semiclassical dynamics is the exact dynamics.

Figure 5 shows a comparison of exact versus semiclassical Gaussian wave packet dynamics propagated in the one-dimensional double well potential $ V(x) = \frac{1}{20} x^{4} - \frac{1}{2} x^{2}$. Over a short time interval, the two models match well. However, over longer periods the exact wave packet interferes with itself and eventually delocalizes over both wells, while the Gaussian wave packet bounces back and forth, with no signs of damping.

Figure 5: Exact versus semiclassical Gaussian wave packet dynamics in a one-dimensional double well potential, showing (a) the electron density and (b) the mean position $ x(t)$ and size $ s(t)$ of the wave packets as a function of time.

In systems where electrons are well-localized, semiclassical dynamics should be valid. However, quantum interference and tunneling effects will not be captured by the approximation.

From the equations of motion, it follows that if we define a total kinetic energy $ T$ as

$\displaystyle T = \sum_{i} \frac{1}{2} m_{\mathrm{elec}} \dot{x}_{i}^{2} + \frac{1}{2} \cdot \frac{3}{4} m_{\mathrm{elec}} \cdot \dot{s}_{i}^{2}$ (1)
then the total energy $ T + V(x, s)$ is a constant of motion. The kinetic energy of motion $ T$ is not to be confused with the electronic kinetic energy $ 3/(2 s^{2})$ which appears in the potential energy $ V(\mathbf{x}, s)$.

Energy expression

In eFF, the total energy is the sum of a Hartree product kinetic energy, a Hartree product electrostatic energy, and an antisymmetrization (Pauli) correction:

$\displaystyle E = E_{\mathrm{ke}} + E_{\mathrm{nuc \cdot nuc}} + E_{\mathrm{nuc \cdot elec}} + E_{\mathrm{elec \cdot elec}} + E_{\mathrm{Pauli}} $
which can be broken down further as follows:
$\displaystyle E_{\mathrm{ke}}$ $\displaystyle =$ $\displaystyle \sum_{i} \frac{3}{2} \ \frac{1}{s_{i}^{2}}$  
$\displaystyle E_{\mathrm{nuc \cdot nuc}}$ $\displaystyle =$ $\displaystyle \sum_{i<j} \frac{Z_{i} Z_{j}}{R_{ij}}$  
$\displaystyle E_{\mathrm{nuc(i)\cdot elec(j)}}$ $\displaystyle =$ $\displaystyle -\sum_{i,j} \frac{Z_{i}}{R_{ij}}\ \mathrm{Erf}\left(\frac{\sqrt{2} \, R_{ij}}{s_{j}}\right)$  
$\displaystyle E_{\mathrm{elec(i) \cdot elec(j)}}$ $\displaystyle =$ $\displaystyle \sum_{i<j} \frac{1}{r_{ij}}\ \mathrm{Erf}\left( \frac{\sqrt{2} \, r_{ij}}{\sqrt{s_{i}^{2} + s_{j}^{2}}}\right)$  
$\displaystyle E_{\mathrm{Pauli}}$ $\displaystyle =$ $\displaystyle \sum_{\sigma_{i} = \sigma_{j}} E(\uparrow \uparrow)_{ij} + \sum_{\sigma_{i} \neq \sigma_{j}} E(\uparrow \downarrow)_{ij}$  
where $ E(\uparrow \uparrow)$ and $ E(\uparrow \downarrow)$ are the Pauli potential functions:
$\displaystyle E(\uparrow \uparrow)_{ij}$ $\displaystyle =$ $\displaystyle \left(\frac{S_{ij}^{2}}{1-S_{ij}^{2}} + (1 - \rho) \frac{S_{ij}^{2}}{1 + S_{ij}^{2}}\right) \Delta T_{ij}$  
$\displaystyle E(\uparrow \downarrow)_{ij}$ $\displaystyle =$ $\displaystyle \frac{\rho S_{ij}^{2}}{1 + S_{ij}^{2}} \Delta T_{ij}$  
and where $ \Delta T_{ij}$ is a measure of the change in the kinetic energy of the electrons upon antisymmetrization, and $ S_{ij}$ is the overlap between the two wave packets, which are both dependent on the scaled sizes $ \bar{s}_{i}$ and $ \bar{s}_{j}$ of the electrons and the scaled distance separating them $ \bar{r}_{ij}$:
    $\displaystyle \Delta T_{ij} = \frac{3}{2}\left(\frac{1}{\bar{s}_{1}^{2}}+\frac{... ...2}+\bar{s}_{2}^{2})-2 \bar{r}_{12}^{2})}{(\bar{s}_{1}^{2}+\bar{s}_{2}^{2})^{2}}$  
    $\displaystyle S_{ij} = \left(\frac{2}{\bar{s}_{i}/\bar{s}_{j}+\bar{s}_{j}/\bar{s}_{i}}\right)^{3/2} \exp(-\bar{r}_{ij}^{2}/(\bar{s}_{i}^{2}+\bar{s}_{j}^{2}))$  
The mixing parameter $ \rho = -0.2$ and the size and distance scaling parameters $ \bar{r}_{ij} = r_{ij} \cdot 1.125$ and $ \bar{s}_{i} = s_{i} \cdot 0.9$ were chosen to reproduce the ground state geometries of small molecules such as $ \mathrm{CH_{4}}$, $ \mathrm{C_{2}H_{6}}$, LiH, and $ \mathrm{B_{2}H_{6}}$. They are the same for all electrons, and for all systems studied.

Electrostatics, Heisenberg uncertainty, and the Pauli principle

The eFF potential energy includes electron kinetic energies and electrostatic interactions between the individual electron densities, computed assuming the electrons act as independent particles.

The electron kinetic energies are a manifestation of the Heisenberg uncertainty principle, and vary as the inverse square of each electron's size. The finite size of a hydrogen atom, for example, is a balance between the electrostatic attraction of the electron to the nuclei, and the kinetic energy of the electron, which prevents the electron from collapsing onto the nucleus (Figure 6).

Figure 6: The balance of an electron's kinetic energy, governed by its size, and its potential energy, which comes from its electrostatic attraction to the proton, determines the size of a hydrogen atom.

A Pauli exclusion term is needed because electrons are not independent particles -- they are indistinguishable fermions, and the overall antisymmetry of the wavefunction upon interchange of any two particles gives rise to a Pauli repulsion between electrons.

The Pauli expression approximates the energy difference between the true antisymmetric wavefunction and the corresponding Hartree product wavefunction (Figure 7). It is derived by considering the difference between the antisymmetric and symmetric combinations of valence bond states, and assuming that kinetic energy differences predominate; analogous approximations have been reported elsewhere [11,12,6,9,4].

Figure 7: Pauli repulsion between two electrons with size $ s$ = 1 bohr, as a function of their separation $ r$. The Pauli repulsion between two same spin electrons increases more steeply at short distances than the electrostatic repulsion.

Dynamic electron mass

In eFF, the parameter $ m_{\mathrm{elec}}$ corresponds to a dynamic electron mass and governs the overall time scale of electron motions. It does not affect the electronic kinetic energy which appears in $ V(x, s)$, and so has no effect on the electron sizes and bond lengths in ground state molecules.

For electrons in harmonic potentials, the dynamic mass equals the true electron mass. However, in general this need not be the case, as in semiconductors, where the effective electron mass is $ \hbar^{2} (\partial^{2}\epsilon(k)/\partial k^{2})^{-1}$, which may be much larger than the true electron mass; or in metals, where quasiparticles have an effective mass determined by many-body theory. Currently, we take $ m_{\mathrm{elec}}$ to be an adjustable parameter.

In our studies on the Auger process in diamondoids, we found that core hole lifetimes were consistent with those obtained from experimental linewidth measurements only when $ m_{\mathrm{elec}} \sim m_{\mathrm{proton}}$. Aside from an overall time-scaling, changing $ m_{\mathrm{elec}}$ did not change the basic events observed: core hole filling, ejection of a secondary electron, and excitation of neighboring electrons. The partitioning and distribution of electron energies were also unchanged.

It is also convenient to use a higher $ m_{\mathrm{elec}}$ when simulating adiabatic dynamics, after Car-Parrinello dynamics, which often use an effective mass of $ \sim 500$ a.u. In the case of the warm dense hydrogen studies [1], we found that changing $ m_{\mathrm{elec}}$ did not affect the thermodynamic quantities averaged over time, and so used $ m_{\mathrm{elec}} = m_{\mathrm{proton}}$.

There are cases where it is essential to use a lower electron mass. For example, an ion passing through a solid excites electrons in the region it passes through, which slows the ion down. Here a rigorous separation of nuclear and electron time scales is needed to obtain the correct stopping power.