A new method to calculate ultrafast electron dynamics in molecules based on matrix product states

. We propose a new method to describe electron dynamics in molecules on the scale of femtoseconds. It is based on factorizing the electronic wave function into a matrix product state and using this factorization to solve the time dependent Schrödinger equation.


Introduction
Electronic correlations are the main source for numerous interesting effects in solid-state physics and quantum chemistry. With experimental facilities like the Linear Coherent Light Source [1] in Stanford, US and the European X-Ray Free-Electron Laser [2] in Hamburg, Germany we are now able to observe effects caused by these correlations on the atto-to femtosecond timescale. Efficient methods to describe ultrafast dynamics of correlated electron systems will be therefore essential for the future correspondence of theory and experiment.
Here, we present benchmark calculations comparing a time dependent method based on matrix product states to calculations including all configurations. As a first step we focus on electron dynamics, however, moving nuclei will be included in a later work.

Matrix Product States
The electronic wave function of a molecule may be written in orbital basis as where σ i ∈ {·, ↑, ↓, ↑↓} are the occupation numbers of orbital i, c σ 1 ···σ L (t) are the expansion coefficients and |σ 1 · · · σ N are the many-body basis states representing the electronic configurations. The orbital basis may be molecular orbitals, atomic orbitals, or any other set of orthogonal orbitals. This form of the electronic state includes the full configurational space (FCI). The exponentially growing number of possible configurations makes working with |ψ(t) so challenging, as it requires too many coefficients to store and handle for medium to large sets of orbitals. Today there are many approximations to reduce the number of configurations included in Eq. 1, for example, the configurational interaction [3] and complete active space [4] methods. An approach that just recently got adapted to quantum chemical problems and which is * e-mail: lars-hendrik.frahm@physnet.uni-hamburg.de * * e-mail: daniela.pfannkuche@physnet.uni-hamburg.de  extremely successful in calculating ground states is a factorization of the coefficient tensor c σ 1 ···σ L (t) into a matrix product state [5] (MPS) with matrices A σ i (t) getting optimized to include the most import configurations only. In most physical cases, using matrices of some limited dimension D approaches the FCI state A σ 1 (t) · · · A σ L (t) ≈ c σ 1 ···σ L (t) efficiently. We start from a small dimension (for example a product state with D = 1) and increase until convergence.

Benchmark
We want to compare results using the FCI state from Eq. 1 to our MPS approach from Eq. 2 with different matrix dimensions. To solve the time-dependent Schrödinger Equation, we apply a Krylov space method [6] with a Krylov space size of 8 and a time evolution step ∆t = 3as. With an equivalent implementation for MPS and FCI states, any observed disagreements between the MPS and the FCI result are exclusively caused by the limited dimension of the matrices in the MPS.

Hydrogen Chain
First benchmark is a well studied example in quantum chemistry density matrix renormalization group [7,8], the one dimensional H 10 molecule with bond distance r = 1.8a 0 . We use molecular orbitals derived from a Hartree-Fock + Pipek-Mezey localization calculation using the Molpro software package [9] on a STO-6G level. The initial state is the Hartree-Fock ground state, meaning 10 electrons occupying the lowest 5 molecular orbitals. We start by comparing the occupation number of the HOMO-3 orbital as an example. The orbital is doubly occupied at t 0 and becomes depopulated by about 5% within the first few atto seconds. In Figure 1 although a complete factorization of the FCI coefficient tensor in Eq. 1 results in significantly larger matrices with D FCI = 462 (employing particle number and total spin symmetry as implemented in CheMPS2 [8]). We see, only a fraction of configurations participate in the dynamics during the first femtosecond and our MPS approach is able to find it. As a more profound comparison between the MPS and the FCI state we also compare the two electron reduced density matrix Γ i jkl (t) = στ ψ(t)|ĉ † iσĉ † jτ c kτ c lσ |ψ(t) . As a measure, we compare the relative error between MPS and FCI matrices with respect to the Frobenius matrix norm. We see in Figure 1 (dashed line, right axis) that the error remains smaller than 5% for D = 100 during the first femtosecond.

Doubly Ionized Water
As a second example we study the electron dynamics in the water molecule. The initial state is Hartree-Fock orbitals derived from the 6-31G basis set, however with two electrons removed from the 1s orbital at the oxygen. Due to the large charge hole at the oxygen, we expect more complex dynamics compared to the previous example. We use conservative Krylov parameters with a space size of 4 and a time evolution step ∆t = 0.3as. Figure 2 shows the dynamics of the occupation number of the HOMO-1 orbital calculated using FCI and MPS with a maximum dimension of D = 100. The MPS approach has all possible configurations included at D FCI = 316.
The HOMO-1 orbital is almost depopulated by 50%, where one electron tries to fill the charge hole at the oxygen atom during the first few atto seconds. Both approaches show this behavior with the MPS approach fitting to the FCI result almost perfectly. The depopulation, as well as the frequency of the following oscillation is well reproduced using MPS.