Solver for the electronic V-representation problem of time-dependent density functional theory

TL;DR

This paper introduces a new solver for the scalar time-dependent Kohn-Sham potential in TDDFT, providing numerical results and insights into $V$-representability, with potential for improved quantum system simulations.

Contribution

The paper presents a novel solver for the time-dependent Kohn-Sham potential and offers a new characterization of $V$-representability for one-electron systems.

Findings

Successful numerical implementation for one- and two-electron systems

New insights into $V$-representability conditions

Potential pathways for improved quantum dynamics simulations

Abstract

One route to numerically propagating quantum systems is time-dependent density functional theory (TDDFT). The application of TDDFT to a particular system's time evolution is predicated on $V$-representability which we have analyzed in a previous publication. Here we describe a newly developed solver for the scalar time-dependent Kohn-Sham potential. We present and interpret the force-balance equation central to our numerical method, describe details of its implementation, and present illustrative numerical results for one- and two-electron systems. A new characterization of $V$-representability for one-electron systems is also included along with possible improvements and future directions.

Figures

Equations

$$n_t=\text{diag}\left({\text{Tr}}_{2\dots N}|{\Psi }_t⟩⟨{\Psi }_t|\right)$$

$$⟨x|s_m⟩=\frac{\text{mbox}sin\left[\frac{\pi \left(x-m\Delta x\right)}{\Delta x}\right]}{\pi \left(x-m\Delta x\right)}$$

$$⟨x=m^{\prime }\Delta x|s_m⟩=\frac{1}{\Delta x}{\delta }_{m,m^{\prime }}$$

$${\mathbf{T}}_{m{m}^{\prime }}$$

$${\mathbf{V}}_{m{m}^{\prime }}=\hat{V}\left(m\Delta x\right){\delta }_{m,m^{\prime }}=\text{diag}({\mathbf{v}}_m)$$

$$T_{\vec{m},{\vec{m}}^{\prime }}=\sum _i^D{T}_{m_i,m_i^{\prime }}\prod _{j\ne i}^D{\delta }_{m_j,m_j^{\prime }}$$

$$⟨x|{\varphi }_t^{(i)}⟩=\sum _m^{(i)}{c}_m⟨x|s_m⟩$$

$${\rho }_{jk}^{{\Phi }_t}=⟨{\Phi }_t|{\hat{a}}_j^{†}{\hat{a}}_k|{\Phi }_t⟩=\sum _i^N{}^{(i)}{\mathbf{c}}_j^{(i)}{\mathbf{c}}_k^{†}$$

$$q_{\vec{m}}=2\Re [\sum _i^N{}^{(i)}{\mathbf{c}}_{\vec{m}}^{\prime }{}^{(i)}{\bar{\mathbf{c}}}_{\vec{m}}^{\prime }-{}^{(i)}{\mathbf{c}}_{\vec{m}}{}^{(i)}{\bar{\mathbf{c}}}_{\vec{m}}^{\prime \prime }]$$

$$\begin{array}[]{l}\,{}^{(i)}{{{\bf c}}}^{\prime}={{{\bf T}}}\;^{(i)}{{{\bf c}}}\\ \,{}^{(i)}{{{\bf c}}}^{\prime\prime}={{{\bf T}}}\;^{(i)}{{{\bf c}}}^{\prime}\end{array}$$

$$i[\hat{V},{\partial }_t{\hat{n}}_j]=\sum _k(v_j-v_k)T_{kj}({\hat{a}}_j^{†}{\hat{a}}_k+{\hat{a}}_k^{†}{\hat{a}}_j).$$

$$i⟨{\Phi }_t|[\hat{V},{\partial }_t{\hat{n}}_r]|{\Phi }_t⟩=\sum _s{\left(-2\mathcal{M}(\Re [{\mathbf{T}}^{\Phi }])\right)}_{rs}{v}_s$$

$$\frac{\partial {\dot{n}}_i}{\partial {\theta }_j}$$

$$J_{ij}=\frac{\partial {\dot{n}}_i}{\partial {\theta }_j}$$

$$\mathbf{\vartheta }=\mathbf{\theta }-(\Delta \mathbf{\theta }-\Delta {\mathbf{\theta }}_p)$$

$$\Delta \mathbf{\theta }={\mathbf{J}}^{-1}(\dot{\mathbf{n}}-{\dot{\mathbf{n}}}^{aim})=-{\mathbf{K}}^{-1}(\dot{\mathbf{n}}-{\dot{\mathbf{n}}}^{aim})$$

$$\frac{\partial {\dot{n}}_j}{\partial {\theta }_k^{(i)}}=-K_{jk}^{(i)},$$

$${\mathbf{K}}^{{\Phi }_t}\Delta {\mathbf{\theta }}^{{\Phi }_t}=\dot{\mathbf{n}}-{\dot{\mathbf{n}}}^{aim},$$

$${\mathbf{\vartheta }}^{(i)}={\mathbf{\theta }}^{(i)}-(\Delta {\mathbf{\theta }}^{(i)}-\Delta {\mathbf{\theta }}_p^{(i)}).$$

$$\mathbf{K}=-2\Re [\mathbf{T}\circ {\mathbf{\rho }}^{\Phi }-\text{diag}\left(\sum _k{\left(\mathbf{T}\circ {\mathbf{\rho }}^{\Phi }\right)}_{jk}\right)]$$

$$-2\Re [\sum _{i=1}^N\text{diag}({}^{(i)}\mathbf{c})\mathbf{T}\text{diag}({}^{(i)}\bar{\mathbf{c}})]$$

$$\mathbf{c}{\mathbf{c}}^{†}=\mathbf{r}{\mathbf{r}}^T+i\mathbf{i}{\mathbf{r}}^T+\mathbf{i}{\mathbf{i}}^T-i\mathbf{r}{\mathbf{i}}^T$$

$$\sum _{i=1}^N-2\left(\text{diag}({}^{(i)}\mathbf{r})\mathbf{T}\text{diag}({}^{(i)}\mathbf{r})+\text{diag}({}^{(i)}\mathbf{i})\mathbf{T}\text{diag}({}^{(i)}\mathbf{i})\right)$$

$$\sum _j{\left(\mathbf{A}\circ \mathbf{B}\right)}_{jk}={\left(A{B}^T\right)}_{jj}$$

$$-2\Re [-\text{diag}\left(\sum _k{\left(\mathbf{T}\circ {\mathbf{\rho }}^{\Phi }\right)}_{jk}\right)]=\text{diag}(2\Re [\sum _i^N{}^{(i)}\mathbf{D}])$$

$${}^{(i)}{D}_{jk}={}^{(i)}{\mathbf{c}}_j^{\prime }{}^{(i)}{\bar{\mathbf{c}}}_k$$

$$M_{\vec{m}{\vec{m}}^{\prime }}=\frac{1}{\left|K_{\vec{m}\vec{m}}\right|}{\delta }_{\vec{m}{\vec{m}}^{\prime }}$$

$${\mathbf{M}}^{1/2}\mathbf{K}{\mathbf{M}}^{1/2}({\mathbf{M}}^{-1/2}\mathbf{v})={\mathbf{M}}^{1/2}\mathbf{s}$$

$${{{\bf M^{\Phi}}}}=\left[\begin{array}[]{cccc}{{{\bf M^{(1)}}}}&0&0&0\\ 0&{{{\bf M^{(2)}}}}&0&0\\ 0&0&\ddots&0\\ 0&0&0&{{{\bf M^{(N)}}}}\end{array}\right]$$

$$M_{\vec{m}{\vec{m}}^{\prime }}^{(i)}=\frac{1}{\left|K_{\vec{m}\vec{m}}^{(i)}\right|}{\delta }_{\vec{m}{\vec{m}}^{\prime }}$$

$${\mathbf{M}}^{1/2}{\mathbf{K}}^{\mathbf{\Phi }}{{\mathbf{M}}^{\mathbf{\Phi }}}^{1/2}({{\mathbf{M}}^{\mathbf{\Phi }}}^{-1/2}\mathbf{\Delta }{\mathbf{\theta }}^{\mathbf{\Phi }})={\mathbf{M}}^{1/2}\left(\dot{\mathbf{n}}-{\dot{\mathbf{n}}}^{aim}\right).$$

Full text

††preprint: APS/123-QED

1. I Introduction

2. I.1 Lattice basis

3. I.2 The force-balance equation

4. II $V^{KS}$-representation solver

5. II.1 Preparing consistent initial KS states

6. II.1.1 One-electron

7. II.1.2 Multi-electron systems

8. II.2 Numerical inversion of the force-balance matrix

9. II.3 Iterative solver

10. II.3.1 Preconditioning

11. II.4 Time Propagator

12. II.5 Obtaining the partial derivative of the potential

13. III One-electron test

14. IV Two-electron test

15. V $V$-representability

16. VI Conclusions

${\scriptstyle\Delta}x$, Spacing between basis functions in position space. 2. 2.

$m$, Index for location of phase-space basis function in position space at $m{\scriptstyle\Delta}x$, not necessarily an integer but $m$ has integer spacing. 3. 3.

$\vec{m}=\left[m_{1},m_{2},...m_{D}\right]$, The vector of all one-dimensional indices for the multidimensional basis function with $D$ dimensions. 4. 4.

$L$, Total number of one dimensional phase space functions. 5. 5.

$N$, Size of basis. 6. 6.

$\hat{H}$, operators have the hat symbol. 7. 7.

${{{\bf H}}}$, Matrices will be capital letters in bold and non-italicized. 8. 8.

${{{\bf v}}}$, Vectors will be small letters in bold and non-italicized. 9. 9.

${{{\bf v}}}^{\dagger}$ is the conjugate transpose of ${{{\bf v}}}$ 10. 10.

$V^{(i)}$: the $(i)$ superscript refers to a specific particle.

I Introduction

Important classes of time-evolution algorithms widely employed by chemists and physicists are based on reduced descriptions of the wave function. These include methods focused on the two-body reduced density matrix and the electron density (the diagonal of the one-body reduced density matrix).

Since all interactions of non-relativistic Hamiltonians are between at most two electrons, the $N$-electron wave function, $\Psi$, contains more information than necessary. For this reason, the two-electron reduced density matrix (2RDM) contains enough information to characterize properties of non-relativistic quantum systems Mazziotti (2012). However to minimize the energy, one must characterize the set of 2RDMs corresponding to a valid $N$-electron wave function otherwise nonphysical 2RDMs allow unbounded solutions. The characterization of the 2RDMs that came from $N$-electron wave functions is known as the $N$-representability problem Coulson (1960). The $N$-representability problem was proven to be QMA-complete Lui et al. (2007) highlighting the theoretical difficulty of 2RDM methods. Nonetheless, there has been successful efforts to perform time evolution using 2RDM methods Lackner et al. (2015).

An even more concise description is afforded by the ground state one-electron probability density,

[TABLE]

which we will simply refer to as the density. The Hohenberg-Kohn theorems dictate that the ground state density, $n_{gs}$, is sufficient to characterize all properties of the quantum system Hohenberg and Kohn (1964). This provides the basis for density functional theory (DFT). While theoretically compelling, many functionals to efficiently compute properties from the density are unknown. Moreover, the universal functional necessary for evaluating the energy is unlikely to be determined numerically even to only polynomially accurate in the size of the system. Despite the numerous approximations to the universal functional, computational complexity arguments Schuch and Verstraete (2009); Whitfield et al. (2014a) showed that obtaining the numerically exact functional is intractable even with quantum computation. It may be possible to obtain a functional using a particular form of the Kohn-Sham orbitals.Hatcher et al. (2019) The corresponding time-dependent result Whitfield et al. (2014b) states, for sufficiently well-behaved systems, the potential can be computed efficiently with access to a quantum computer. Unlike the ground state result, the time-dependent complexity analysis relied on the Kohn-Sham (KS) construction lying at the heart of nearly all practical schemes for DFT. In this letter, we return to the analysis begun in our previous work Whitfield et al. (2014b) using a combination of theory and numerics.

The Kohn-Sham (KS) construction considers a system of non-interacting electrons that can be used to replace a given interacting system of electrons. An electronic system with with two-body interactions, $\hat{W}$ is governed by Hamiltonian $\hat{\bf T}+\hat{\bf V}+\hat{\bf W}$. The KS system is a new collection of non-interacting electrons with a new additional potential term $\hat{\bf V}^{KS}$. The new potential in the KS Hamiltonian, $\hat{\bf T}+\hat{\bf V}+\hat{\bf V}^{KS}$, is selected such that the one-electron probability density of the non-interacting system matches that of the given interacting system of interest.

Let us define the general $V$-representation problem of a KS system as the task of constructing a model system which has the same expectation values on selected observables as a target system. $V$-representability refers to the existence of solutions to this problem when different constraints are placed on the model system. The time-dependent subset of $V$-representation problems considers as input an initial state and the target trajectory of selected observables, and the task is to find the correct time-dependent fields for a specified control Hamiltonian. This general framework is not limited to electronic systems as illustrated by a study of this problem in the context of spin systems Tempel and Aspuru-Guzik (2012). Here, attention will focus on the electronic $V$-representation problem where the tasks is to construct a KS system governed by a time-dependent potential $V^{KS}(t)$ such that the KS density matches the density of a specified interacting many-electron system at all times. We refer to this as the $V^{KS}$-representation problem.

The $V^{KS}$-representation problem associated with the model KS system is the task of constructing the potential for a model system such that the non-interacting model attains a target density expectation value. This task is not always possible and questions of $V$-representability refer to the existence of solutions to the $V^{KS}$-representation problem when different constraints are placed on the KS or on the target density. Here, attention will primarily focus on the time-dependent $V^{KS}$-representation problem where the task is to construct a time-dependent potential $V^{KS}(t)$ such that the KS density matches a given target density at all times.

For fermionic simulations, Ref. van Leeuwen (1999) was first to give a constructive solution to time-dependent $V^{KS}$-representation problem. This was challenged in Ref. Baer (2008) where counter-examples were presented to this construction. These counter-examples were largely addressed by Refs. Li and Ullrich (2008); Farzanehpour and Tokatly (2012) through a detailed analysis of densities evolving on lattices. In separate work, implicit solutions using a fixed point mapping has been formulated directly in the continuum limit using various techniques Ruggenthaler and van Leeuwen (2011); Ruggenthaler et al. (2012); Ramsden and Godby (2012); Nielsen et al. (2013, 2018); Hodgson et al. (2013, 2016). There is also a PDE constrained method Jensen and Wasserman (2016). Here, we will present an explicit method based on the lattice algorithm analyzed in Whitfield et al. (2014b).

The paper begins with the description of the lattice basis used, followed by the description of the force-balance equation. We then turn towards the implementation details of the solver for the $V^{KS}$-representation problem. A novel method for preparing consistent initial KS states for one-electron systems is given. We give some numerical examples before discussing single-electron $V$-representability theorems. Finally, we draw conclusions and point to next steps as we close the article.

I.1 Lattice basis

The one dimensional basis functions $s_{m}\left(x\right)=\langle x|s_{m}\rangle$ that we use are defined in Ref Colbert and Miller (1992) as

[TABLE]

Here the functions are defined on domain $-\infty<x<\infty$, $m=0,\pm 1,\pm 2,...$. Each of the basis functions possess the property that

[TABLE]

If infinitely many basis functions are used, the kinetic energy operator matrix ${{{\bf T}}}$ is given exactly as

[TABLE]

which is equivalent to an infinite order finite difference approximation Colbert and Miller (1992). The potential energy matrix (${{{\bf V}}}$) in this basis is

[TABLE]

where $\hat{V}\left(x\right)$ is the potential operator for the given system and ${{{\bf v}}}$ is the vector of diagonal elements of ${{{\bf V}}}$.

In practice, the range of $m$ is restricted to a finite number which introduces errors to ${{{\bf T}}}$ and ${{{\bf V}}}$ but contribute minimally if $m$ extends over the range where the wavefunction is appreciable. The total number of one-dimensional grid points we use is given by $L$ and is always an odd integer such that $m=0,\pm 1,\pm 2,...,\pm\frac{L-1}{2}$.

In multiple dimensions, a product basis is used such that $s_{\vec{m}}\left(\vec{x}\right)=s_{m_{1}}\left(x_{1}\right)s_{m_{2}}\left(x_{2}\right)...s_{m_{D}}\left(x_{D}\right)$, where $\vec{m}=\left[m_{1},m_{2},...,m_{D}\right]$, $\vec{x}=x_{1},x_{2},...,x_{D}$. The total number of basis functions is $L^{D}$. The kinetic energy operator matrix is then given by

[TABLE]

The potential energy matrix remains diagonal with elements ${{{\bf V}}}_{\vec{m}\vec{m}}=\hat{V}\left(m_{1}{\scriptstyle\Delta}x,m_{2}{\scriptstyle\Delta}x,...m_{D}{\scriptstyle\Delta}x\right)={{{\bf v}}}_{\vec{m}}$. In practice, this product basis will only be viable up to around four or five dimensions. For more dimensions, it will be necessary to prune the indices $\vec{m}$. The scaling of matrix-vector products in this basis is $\mathcal{O}\left(DL^{D+1}\right)$ for the kinetic energy operator and $\mathcal{O}\left(L^{D}\right)$ for the potential (as it is diagonal)Larsson et al. (2016).

Before discussing the force-balance equation in the next section, we must introduce a few notations and definitions. First, we will denote the KS system’s wave functions as $\Phi_{t}$. The one-particle wave functions are defined via an expansion over the basis functions listed above as

[TABLE]

Corresponding to each basis function $s_{m}$, there exists a creation operator $\hat{a}_{m}$ with the property the $\hat{a}_{m}\hat{a}_{j}^{\dagger}+\hat{a}_{j}^{\dagger}\hat{a}_{m}=\delta_{mj}$ and $\hat{a}_{j}\hat{a}_{i}+\hat{a}_{i}\hat{a}_{j}=0$. These anti-commutation relations ensure that a wave function given by $\prod_{j}\left(\hat{a}_{j}^{\dagger}\right)^{k_{j}}|vac\rangle$ is antisymmetric when $|vac\rangle$ is the state with no particles occupied and $k_{j}\in\{0,1\}$.

I.2 The force-balance equation

The non-interacting $V^{KS}$-representation problem requires that the fictitious system’s wave function, $|\Phi_{t}\rangle$, evolve such that its density expectation value, $\langle\Phi_{t}|\hat{\bf n}|\Phi_{t}\rangle$, matches a target evolution, ${\bf n}^{aim}(t)$. Here, and throughout, $\hat{\bf n}$ is the operator-vector with components $\hat{n}_{j}=\hat{a}_{j}^{\dagger}\hat{a}_{j}$. The force-balance equation determines the required instantaneous potential to correctly construct the KS system. Note that forces enter at second order of evolution as anticipated by Newton’s law: $\bm{F}=m\bm{a}$.

The force-balance equation is easily derived from the second derivative of the density following the Heisenberg equation Whitfield et al. (2014b); Farzanehpour and Tokatly (2012). If we aim for a target evolution, then we should have that $\partial_{t}^{2}{\bf n}^{aim}$ (from the target interacting system) is equal to $i\langle\Phi_{t}|[\hat{\bf H},\partial_{t}\hat{\bf n}]|\Phi_{t}\rangle$ (from the KS system). Expanding this commutator, we have two terms, $i[\hat{\bf T},\partial_{t}\hat{\bf n}]$ and $i[\hat{\bf V},\partial_{t}\hat{\bf n}]$ which we will physically interpret as well as giving some guidance on numerical implementation.

We first discuss the acceleration which the forces must cause. The free acceleration, also called the momentum-stress tensor van Leeuwen (1999), $\hat{q}_{x}=i[\mathbf{\hat{T}},\partial_{t}\hat{n}_{x}]=-[\hat{\bf T},[\hat{\bf T},\hat{n}_{x}]]$ is independent of the potential operator. To evaluate the free acceleration, we must evaluate the expression $\hat{\bf q}=2\Re[\textrm{diag}\left(\hat{\bf T}\bm{\hat{\rho}}\hat{\bf T}+\bm{\hat{\rho}}\hat{\bf T}^{2}\right)]$ evaluated against state $\Phi_{t}$. This can be evaluated efficiently using matrix-vector products. The expectation value of the operator $\rho$ is

[TABLE]

where $N$ is the number of non-interacting electrons in the Kohn-Sham system. Evaluating $\bf q$ for each particle before obtaining the real part is

[TABLE]

where

[TABLE]

and $\bar{{{{\bf c}}}}$ signifies the complex conjugate of ${{{\bf c}}}$. Calculating $\bf q$ has a total scaling cost of $\mathcal{O}\left(2NDL^{D+1}\right)$.

The forced acceleration is given by the difference between the free acceleration and the target acceleration ${\bf s}=\partial_{t}^{2}{\bf n}^{aim}-{\bf q}$. The forced acceleration then determines the forces required from the potential.

The forces enter through the term $i[\hat{V},\partial_{t}\hat{\bf n}]$ which can be recast into two useful forms; one illustrating a connection to forces and the other geared towards determining the potential. The first form we examine is

[TABLE]

This form gives a nice analogy to the real space forces as $F(x)=-\nabla V(x)$. We note that $i\langle\Phi_{t}|[\hat{V},\partial_{t}\hat{n}_{j}]|\Phi_{t}\rangle=2\sum_{k}({v}_{j}-{v}_{k})\;\Re[{T}_{kj}^{\Phi_{t}}]$ has the form of a generalized discrete gradient. Here, the real part of ${T}^{\Phi_{t}}_{jk}={T}_{jk}\langle\Phi_{t}|\hat{a}_{j}^{{\dagger}}\hat{a}_{k}|\Phi_{t}\rangle$ includes both the influence of the probability mass at each site as well as the underlying spatial metric.

The second form is more applicable to numerical simulation: $i\langle\Phi_{t}|[\hat{V},\partial_{t}\hat{n}_{j}]|\Phi_{t}\rangle$ $=\sum_{r}\langle\Phi_{t}|\left(i[\hat{n}_{r},\partial_{t}\hat{n}_{j}]\right)|\Phi_{t}\rangle{v}_{r}$. To introduce a simpler expression, it will be advantageous to define $\mathcal{M}(\bm{A})_{ij}={A}_{ij}-\delta_{ij}\left(\sum_{k}{A}_{jk}\right)$. Then, the force-balance equation can be expressed as:

[TABLE]

For consistency with our publication Whitfield et al. (2014b), we define ${{{\bf K}}}=-2\mathcal{M}\left(\Re[{{{\bf T}}}^{\Phi}]\right)$ as the force-balance operator. For symmetric matrices, ${{{\bf A}}}={{{\bf A}}}^{T}$, $\mathcal{M}({{{\bf A}}})$ will have the constant vector in its null space. The gauge freedom physically stems from the irrelevance of the zero of energy. Since we are concerned with time-dependent quantum mechanics, the constant potential only imprints an unobservable global phase on the wave function.

II $V^{KS}$-representation solver

Elsewhere Whitfield et al. (2014b), we presented and analyzed an explicit solution for the time-dependent potential necessary for TDDFT provided with the density time-trace of a $V$-representable system. The algorithm was found to scale polynomially in all input parameters except for the $V$-representability parameter Whitfield et al. (2014b) which diverges exponentially when the interacting system no longer has a corresponding KS system.

The algorithm requires, as inputs, the complete time-trace of the density and a consistent initial state $\Phi$ that reproduces the initial density and the initial time derivative of the density. For numerical implementation, a representation of the kinetic operator in the lattice basis is also needed.

II.1 Preparing consistent initial KS states

II.1.1 One-electron

Unlike previous numerical studies Li and Ullrich (2008); Baer (2008); Verdozzi (2008); Ruggenthaler et al. (2012); Schmitteckert et al. (2013); Nielsen et al. (2013); Fuks and Maitra (2014), we consider preparation of initial KS states with non-zero initial momentum. Suppose that the initial KS system should yield an initial time derivative given by $\dot{{{{\bf n}}}}^{aim}$. Using the Heisenberg equation of motion it can be shown Whitfield et al. (2014b), that for a one-electron wave function in a lattice basis, $\psi$, the density derivative is given by $\dot{n}_{j}=-i\sum_{k}{T}_{ki}(\langle\hat{a}_{i}^{\dagger}\hat{a}_{k}\rangle_{\psi}-\langle\hat{a}_{k}^{\dagger}\hat{a}_{i}\rangle_{\psi})$. Given that we consider a single electron wave function, the 1-RDM elements can be written as $\langle\hat{a}_{i}^{\dagger}\hat{a}_{j}\rangle_{\psi}=\sqrt{n_{i}n_{j}}e^{i(\theta_{j}-\theta_{i})}$.

To assign the phases, roughly speaking, we must solve the equation: $\Delta\bm{\theta}=\frac{\partial\bm{\theta}}{\partial\dot{{{{\bf n}}}}}\Delta\dot{{{{\bf n}}}}$. This is the content of Newton’s method. We will describe the modifications needed to handle the gauge degree of freedom after deriving the Jacobian, $J_{ij}=\frac{\partial\dot{n}_{i}}{\partial\theta_{j}}$.

The Jacobian is also given by (minus) the force-balance matrix:

[TABLE]

Before applying the Newton method, we must account for the gauge corresponding to the global phase of the wave function. Other manifestations of this gauge degree of freedom are 1) the one-dimensional null space of the Jacobian and 2) the constraint that $\sum\dot{n}_{j}=0$.

We can fix the phase of one of the $M$ wave function components in order to fix the gauge. Suppose the fixed phase is the component with the largest $n_{m}$ at position $p$ (which we assume is the best described lattice point), then we only update $\vec{\theta}_{g}$ on the $M-1$ remaining components. We have found this to be superior to just choosing an arbitrary component to fix the phase.

Putting it all together, the Newton rule for updating the phase vector $\vec{\theta}$ to a new assignment $\bm{\vartheta}$ is

[TABLE]

where component $\Delta\bm{\theta}_{p}$ is subtracted so that the $p$th component is not updated. We have defined

[TABLE]

Note that dropping the gauge component to get $\dot{n}_{g}$ and $\dot{n}^{aim}_{g}$ loses no information due to the constraint that $\sum\dot{n}_{j}=0$.

Numerical results with a straightforward implementation of Newton’s method works quite well provided that the wave function’s initial momentum is somewhat close to the target momentum.

II.1.2 Multi-electron systems

For $N$ non-interacting one-electron wave functions in a lattice basis with wavefunctions $\phi^{(i)}$, the total density derivative is given by $\dot{n}_{j}=-i\sum_{i}^{N}\sum_{k}{T}_{kj}(\langle a_{i}^{\dagger}a_{k}\rangle_{\phi^{(i)}}-\langle a_{k}^{\dagger}a_{j}\rangle_{\phi^{(i)}})$. The 1-RDM elements can be written as $\langle a_{j}^{{\dagger}}a_{k}\rangle_{\Phi}=\sum_{p}^{N}\sqrt{{}^{(p)}n_{j}\,^{(p)}n_{k}}e^{{i\mkern 1.0mu}(\theta_{k}^{(p)}-\theta_{j}^{(p)})}$.

Following the section above, the Jacobian is given as

[TABLE]

where ${{{\bf K}}}^{(i)}$ denotes the force-balance matrix that only includes the i-th one-electron wavefunction of the $N$ non-interacting wavefunctions. To obtain the next approximation of phases, we need to solve the under-determined least squares problem

[TABLE]

where ${{{\bf K}}}^{\Phi_{t}}=\left[{{{\bf K}}}^{(1)},{{{\bf K}}}^{(2)},...,{{{\bf K}}}^{(N)}\right]$ and $\Delta\bm{\theta}^{\Phi_{t}}=\left[\Delta\bm{\theta}^{(1)},\Delta\bm{\theta}^{(2)},...,\Delta\bm{\theta}^{(N)}\right]^{T}$. In this paper, the solution is obtained iteratively using a least squares method as implemented by LABEL:SOL or a Morse-Penrose pseudo inverse. In two and three dimensions, the iterative method will likely be preferred. The Newton steps taken for each particle are then

[TABLE]

Here $\Delta\bm{\theta}_{p}$ is the zero vector with the $p$th entry equal to $\Delta\theta_{p}$. Note that the gauge needs to be fixed individually for all $N$ electrons.

II.2 Numerical inversion of the force-balance matrix

The inversion of the force-balance matrix has a few problems. It always includes a null space of the constant vector and can be poorly conditioned when parts of the density are very small. If the KS system does not have a single potential corresponding to a given momentum-stress tensor, then other singular vectors can be introduced.

In order to exploit the structure of the force-balance matrix, and the corresponding $\hat{\rho}$, we use iterative Krylov-space methods to solve for ${{{\bf K}}}{{{\bf b}}}={{{\bf x}}}$. The 1RDM $\bm{\rho}$ is a rank-$N$ matrix with matrix elements given in Eq. (10). The force-balance matrix is

[TABLE]

where $\textrm{diag}({{{\bf v}}})$ is a diagonal matrix with elements $v_{j}$ as the $j$th entry and $\circ$ denoting the Hadamard product.

Using a property of Hadamard products of matrices, the first term of Eq. (21) can be evaluated as

[TABLE]

To obtain the real part only, we note that

[TABLE]

where ${{{\bf r}}}$ and ${{{\bf i}}}$ are real vectors such that ${{{\bf c}}}={{{\bf r}}}+{i\mkern 1.0mu}{{{\bf i}}}$. This means that the complex rank-1 matrix ${{{\bf c}}}{{{\bf c}}}^{\dagger}$ is equivalent to the sum of two real rank-1 matrices and two imaginary rank-1 matrices. By retaining only the real matrices, Eq. (22) can be rewritten as

[TABLE]

This term can be applied to a vector in $\mathcal{O}\left(2N(L^{D}+DL^{D+1}+L^{D}\right)$ by applying the matrices sequentially.Brown and Whitfield (2018)

For the second term in Eq. (21), it is best to obtain the diagonal elements of the matrix once. This is done using the property of Hadamard productsStyan (1973)

[TABLE]

This means that the second term can be evaluated as

[TABLE]

where

[TABLE]

where as before, $\,{}^{(i)}{{{\bf c}}}^{\prime}={{{\bf T}}}\,^{(i)}{{{\bf c}}}$, and ${}^{(i)}\bar{{{{\bf c}}}}$ signifies the complex conjugate of ${}^{(i)}{{{\bf c}}}$.

In one dimension, there is no advantage to performing the matrix-vector products this way. However, in multiple dimensions, the scaling of the matrix-vector products is much better and the storage requirements for the matrix are also better with order $\mathcal{O}(L+(2N+1)L^{D})$ instead of $L^{2D}$.

II.3 Iterative solver

In order to solve for the potential, we use the MINRES-QLP algorithmChoi et al. (2011). This method is designed to solve singular Hermitian problems like that of ${{{\bf K}}}{{{\bf v}}}=\partial_{t}^{2}{\bf n}^{aim}-\bf{q}$.

For obtaining the initial state phases, we use the LSQR algorithmPaige and Saunders (1982). This algorithm requires the application of ${{{\bf K^{\Phi}}}}\bf x$ and ${{{\bf K^{\Phi}}}}^{T}\bf y$ to solve the matrix problem ${{{\bf K^{\Phi}}}}\bf x=y$. This poses no problem when exploiting the structure of the force-balance matrix. Applying ${{{\bf K^{\Phi}}}}$ involves applying each ${{{\bf K}}}^{(i)}$ to the corresponding segment of $\bf x$ while ${{{\bf K^{\Phi}}}}^{T}\bf y$ involves applying each ${{{\bf K}}}^{(i)}$ to ${{{\bf y}}}$ in succession.

II.3.1 Preconditioning

Without preconditioning, the force-balance matrix is poorly conditioned due to the large changes in the density on the grid. In our testing, it appears that using a diagonal preconditioner works quite well. This is defined as

[TABLE]

where the new problem is

[TABLE]

For the assigning of phases, the other half of the preconditioning matrices is

[TABLE]

where

[TABLE]

with the problem now,

[TABLE]

If the value of $K_{\vec{m},\vec{m}}$ is zero (or $K^{(i)}_{\vec{m},\vec{m}}$) drops below $10^{-15}\times\mbox{max}\left({{{\bf K}}}\right)$ (or $10^{-15}\times\mbox{max}\left({{{\bf K}}}^{(i)}\right)$), then $M_{\vec{m},\vec{m}}$ is set to zero. The density at lattice point $\vec{m}$ is then deemed too small to resolve a potential and is removed from consideration. As both MINRES-QLP and LSQR are for singular problems, this poses no difficulties. In the future, it would be more efficient to prune those points from the required matrix-vector products to improve efficiency Brown and Whitfield (2018); Larsson et al. (2016).

II.4 Time Propagator

The method we use to solve the TDSE is based on LABEL:vanDijk2014. The Schrödinger equation that is solved is defined as

[TABLE]

where ${{{\bf H_{0}}}}$ is the time-independent portion of the Hamiltonian (that we assume is given) and $\tilde{{{{\bf V}}}}\left(t\right)$ is the time-dependent potential (which is the goal) at time $t$. ${{{\bf H_{0}}}}$ only includes ${{{\bf T}}}$, but can also include the external one-electron potential if it is given as part of the problem instance and is time-independent.

This method propagates wavefunctions using,

[TABLE]

where

[TABLE]

We only retain the next highest Euler-MacLaurin expansion term in $\bm{f}\left(t:\bm{\psi}\left(t\right)\right)$ which for our purposes is defined as,

[TABLE]

with Eq. (33) defining $\frac{\partial\bm{\psi}\left(t\right)}{\partial t}$ and the method to obtain $\frac{\partial\tilde{{{{\bf V}}}}\left(t\right)}{\partial t}$ described later.

Eq. (34) is an implicit equation in that ${\bm{f}}\left(t:\bm{\psi}\left(t\right)\right)$ depends on $\bm{\psi}\left(t\right)$. Here, both ${\bm{f}}\left(t:\bm{\psi}\left(t\right)\right)$ and $\tilde{{{{\bf V}}}}\left(t\right)$ are defined by $\bm{\psi}\left(t\right)$. Therefore the algorithm to update the potential and the wavefunction is an iterative process where,

[TABLE]

where $\bm{\psi}^{(+)}\left(t\right)=e^{-{i\mkern 1.0mu}{{{\bf H_{0}}}}\Delta t}\bm{\psi}^{(-)}\left(t-\Delta t\right)$ is calculated once using the short iterative Arnoldi propagator. Initially, ${\tilde{{{{\bf V}}}}^{(0)}\left(t\right)}=\tilde{{{{\bf V}}}}(t-\Delta t)$ and $\bm{\psi}^{(0)}\left(t\right)=\bm{\psi}\left(t-\Delta t\right)$. We first update only ${\bm{f}}\left(t:\bm{\psi}^{(i)}\left(t\right)\right)$ while keeping ${{{{\bf\tilde{V}}}}^{(k)}(t)}$ unchanged and iterate Eq. (37) until $\left|\bm{\psi}^{(k+1)}-\bm{\psi}^{(k)}\right|$ is below machine precision. We then update ${{{\bf V}}}^{(t)}$ using Eq. (29) and update ${\bm{f}}\left(t:\bm{\psi}^{(i)}\left(t\right)\right)$ followed by the iteration of Eq. (37). This process is repeated until the difference in $\bm{\psi}\left(t\right)$ before and after updating ${{{\bf V}}}\left(t\right)$ is below some threshold. Note that the time-independent external potential will need to be subtracted from the result of Eq. (29) if it was given as part of the problem instance specification.

II.5 Obtaining the partial derivative of the potential

The definition of ${\bm{f}}\left(t:\bm{\psi}^{(i)}\left(t\right)\right)$ includes a partial derivative of $\tilde{{{{\bf V}}}}(t)$ with respect to time. One can introduce an error in Eq. (37) by assuming that $\frac{\partial\tilde{{{{\bf V}}}}(t)}{\partial t}=0$. However, we use a finite difference approximation to $\frac{\partial\tilde{{{{\bf V}}}}(t)}{\partial t}$ which improves the accuracy of the results.

We initially define the zero of energy by some value $E_{0}$, which is usually $E_{0}=0$. Each time a potential $\tilde{{{{\bf V}}}}(t)$ is obtained, we shift the result such that the total energy of the system remains at $E_{0}$. This total energy is defined as

[TABLE]

and results in the correctly shifted potential of $\tilde{{{{\bf V}}}}(t)\rightarrow\tilde{{{{\bf V}}}}(t)-(E-E_{0})/N$. This method appears to correctly position the potential for the systems examined. The method would need to be modified if the multi-electron system has a source term imparting energy on the system.

To obtain the derivative of the potential, we use the backward finite-difference approximation. Therefore, $\frac{\partial\tilde{{{{\bf V}}}}(t=0)}{\partial t}=0$. $\frac{\partial\tilde{{{{\bf V}}}}(t=\Delta t)}{\partial t}=(\tilde{{{{\bf V}}}}(t)-\tilde{{{{\bf V}}}}(t-\Delta t))/\Delta t$, and $\frac{\partial\tilde{{{{\bf V}}}}(t=j\Delta t)}{\partial t}=(\tilde{{{{\bf V}}}}(t)-4\tilde{{{{\bf V}}}}(t-\Delta t)+\tilde{{{{\bf V}}}}(t-2\Delta t))/(2\Delta t)$ for $j>1$. Higher order finite-difference approximations did not improve results.

The method and techniques outlined above generally behave well in regions where the density is appreciable. However, it is helpful for stability to assign the phases after every time step which requires that the first derivative of the density is provided. This should not be a problem as the second derivative of the density is required for potential inversion. The first derivative can also be obtained by measuring the full one-body reduced density matrix. The first derivative of the density is then given as $i\partial_{t}\hat{\bm{\rho}}=\left[\hat{\bm{T}},\hat{\bm{\rho}}\right]$.

One stability issue that can occur is very large changes in the inverted potential where the density is small due to numerical artifacts. This can be attenuated by setting a maximum allowed change in the potential for each lattice point at each time step.

III One-electron test

The one-electron test we examine is a sum of two coherent wavepackets in a Harmonic potential. The potential is constant at all times as

[TABLE]

where $\omega=1$. The wavefunction initial state is

[TABLE]

where $A$ is a normalization factor. This is the superposition of the two coherent states with momentum of $\pm 5\sqrt{\pi}/2$. The grid used has 115 points and covers a range of $x\in\left[-11,11\right]$ with 1600 time steps to obtain a total time of $T/2$ where $T=2\pi/\omega$. The potential inversion is performed using a pseudoinverse of Eq. (29) instead of MINRES-QLP. For one-dimension, there is no advantage to using MINRES-QLP.

The exact solution is obtained by propagating the initial wavepacket using the Short Iterative Arnoldi (SIA) propagator. Snapshots of the method’s performance are shown in FIG. 1. At $t=0$ the density error is zero by definition. The potential error is inversely related to the magnitude of the density, or the error is smallest where the density is largest. At $t=T/4$, the wave packet is split into two distinct regions. Once again, the potential error is small where the density is large and larger where the density is small. However, there is another feature where a constant shift error occurs on the right side of FIG . 1(d). This is due to the inversion method not being able to resolve the middle (low density) region enough to connect the two sides with the same potential offset. At $t=T/2$, this offset error contributes to some inaccuracy in the recovered potential. The error in the recovered potential is approximated bounded by the shift error of $\approx 10^{-4}$ from $t=T/4$. However the results are still reasonable. If the initial state has a larger momentum than $\pm 5\sqrt{\pi}/2$, then the density connecting the two wavepackets tends to zero at $t=T/4$ and larger shift errors can occur. It would then be necessary to separate the problem into two smaller systems with different energies. The splitting of the system could be defined with the use of the preconditioner.

This one electron test shows that our method recovers the potential accurately under a variety of circumstances. It succeeds for large changes in density at individual grid points and can “discover” the correct potential as the density explores the physical space. It also can resolve the potential for two higher density areas with a small density in between, and works for complicated densities with non-negligible momentum.

IV Two-electron test

The two electron test that we examine is that first described in LABEL:Hodgson2013. The initial state is the first triplet state of

[TABLE]

where $V_{d}(x)=\alpha x^{10}-\beta x^{4}$, with $\alpha=5\times 10^{-11}$ and $\beta=1.3\times 10^{-4}$. At $t\geq 0$, the system becomes

[TABLE]

where $V_{ext}(x_{i})=V_{d}(x_{i})-x_{i}/10$ which drives the electrons to the right. The first triplet state is calculated using ARPACK to obtain the lowest energy state. Anti-symmetry of the spatial wave function was enforced for the starting vector and after each matrix-vector product.

The initial state for KS system was obtained by obtaining the full 1RDM of the two-electron initial triplet state and performing a singular value decomposition (SVD). The SVD of a $L\times L$ matrix obtains the optimal truncated $R$-rank representation

[TABLE]

where ${{{\bf M}}}={{{\bf U\Sigma V^{*}}}}$, ${{{\bf u}}}_{i}$ is the $i$th column of ${{{\bf U}}}$, ${{{\bf v}}}_{i}$ is the $i$th column of ${{{\bf V}}}$ and $\sigma_{i}$ is the $i$th largest singular value. As the 1RDM is hermitian positive semi-definite, ${{{\bf v}}}_{i}={{{\bf u}}}_{i}$. These are also the eigenvectors (natural orbitals) of the 1RDM.

In the case of the 1RDM, each $\sigma_{i}$ is a number that represents the occupancy of the $i$th natural orbital, ${{{\bf u}}}_{i}$. For the interacting system examined at $t=0$ of Eq. (41), the first two singular values were degenerate with occupancy of $0.999999799989864$. This signifies that the system is very close to non-interacting initially (i.e. describable by a state with $\sigma_{i}$ taking only values zero or one).

To obtain initial KS orbitals, an iteration scheme is implemented. Initially, the orbitals of the KS system are taken as the first two singular vector ${{{\bf u}}}_{1}$ and ${{{\bf u}}}_{2}$. Element by element, if $u_{1i}<{\bf n}^{(1)}_{i}\rightarrow u_{2i}=\mbox{sgn}(u_{2i})\sqrt{{\bf n}^{(1)}_{i}-u_{1i}^{2}}$. Likewise, if $u_{2i}<{\bf n}^{(1)}_{i}\rightarrow u_{1i}=\mbox{sgn}(u_{1i})\sqrt{{\bf n}^{(1)}_{i}-u_{2i}^{2}}$. The vectors ${{{\bf u}}}_{1}$ and ${{{\bf u}}}_{2}$ are orthonormalized after each sweep through the grid points. After a number of iterations, the vectors converge such that densities ${\bf n}^{(1)}$ and ${\bf n}^{(KS)}$ exactly match and ${{{\bf u}}}_{1}$ and ${{{\bf u}}}_{2}$ are orthonormal. Using the natural orbitals of the 1RDM for an interacting system should produce reasonable starting vectors for any number of electrons along with the iterative method described here.

The two-electron system is propagated using the short iterative Arnoldi method using a Krylov subspace of $50$ with exchange anti-symmetry enforced throughout the calculation. The basis functions are on a $271\times 271$ grid in the range of $x\in\left[-13.5,13.5\right]$ with ${\scriptstyle\Delta}x=0.1$ and $dt=0.005$. The KS system is propagated using the method of section II.4 with the phases assigned at every step. The inverse of ${{{\bf K}}}$ was obtained using the preconditioner ${{{\bf M}}}$ and MINRES-QLP. The inversion of ${{{\bf K}}}$ is still difficult even when using the preconditioner for this system. There is a very small ($\approx 10^{-10}$) eigenvalue of ${{{\bf M}}}^{-1/2}{\bf K}{\bf M}^{-1/2}$. All the other eigenvalues are $>10^{-3}$ aside from the constant vector which is known to be in the null space.

The results of the potential inversion are shown in FIG. 2. The initial $V^{KS}-V_{ext}$ is very similar to that of LABEL:Hodgson2013. It is slightly different due to the fact that the driving potential is included in our calculation at $t=0$ and we do not use vector potentials. The times of $t=3.5$ and $t=5.3$ are included to compare with the previous studies of Refs. Hodgson et al., 2013; Jensen and Wasserman, 2016. LABEL:Hodgson2013 obtains different Kohn-Sham potential in the middle region but still had the same $1/r$ decay. The primary difference is due to the fact that they used vector potentials and we have not. LABEL:Jensen2016 could not recover a potential for this system.

Our method fails around $t=5.75$ as shown in FIG.2(d). We suspect this is due to the formation of a cusp that prevents the solution from being continued. Shockwave formation in partial differential equations are known to cause similar issues [ref]. The cusp in the density forms at $x\approx 6.8$ starting at about $t\approx 5.70$ and appears for about 0.2 time units before disappearing. LABEL:Hodgson2013 does not discuss this region in time so it is unclear if their method can circumvent this issues. Despite the failure to continue the solution through the formation of this sharp feature in the density, we can easily restart the solver at time $5.9$ without further issues.

Note that we suspect that a cusp has formed since oscillatory behavior did not present a problem our method as seen in the one-electron test where the density fluctuates over its support.

V $V$-representability

We now make a few comments on the existence of solutions to the $V^{KS}$-representation problem. A previous theorem Farzanehpour and Tokatly (2012) showed that ground states of many-body interacting systems are always $V$-representable in the neighborhood of the initial time. They showed that the matrix ${\bf K}$ has only one zero eigenvalue and is positive definite in the space of inhomogeneous potentials. It should be noted that the theorem does not characterize the $V$-representability parameter thus numerical stablility is not ensured.

Here we use simpler arguments to provide additional characterizations of the spectrum of ${\bf K}$ in the single-electron case. This coincides with the previous theorem for the ground state but generalizes to all eigenstates.

Theorem: Given non-degenerate $\bm{\psi}$ such that ${\bf H}\bm{\psi}=E_{k}\bm{\psi}$ for ${\bf H}={\bf T}+{\bf V}\in\mathbf{M}(\mathcal{R})$, $\mathbf{K}(\bm{\psi})$ has $k-1$ negative eigenvalues and $M-k+2$ positive eigenvalues.

Proof: Assuming that $\bm{\psi}$ is an eigenstate of ${\bf H}={\bf T}+{\bf V}$ with eigenvalue $\lambda$, then $\mathbf{H}^{(\lambda)}=\mathbf{T}-(\lambda\mathbf{1}-\mathbf{V})=\mathbf{T}-\mathbf{D}$ has $\bm{\psi}$ in its null space. Rearranging, $(\mathbf{T}+\mathbf{D})\bm{\psi}=0$ implies $\sum_{k}T_{jk}\psi_{k}=D_{j}\psi_{j}$. Before using the definition of the force-balance equation, it is important to note that eigenvectors of symmetric matrices are real. Hence, for a single-particle in the eigenstate $\bm{\psi}$:

[TABLE]

Since $\psi$ has full support, $S_{mj}=\delta_{mj}\psi_{j}$ is non-singular and the number of $(+/0/-)$ eigenvalues are the same for $\mathbf{H}^{(\lambda)}$ and $\mathbf{K}(\bm{\psi})$ by Sylvester’s theorem Horn and Johnson (2005). $\square$

According to numerical tests, if $\bm{\psi}$ does not have full support then $\mathbf{K}$ will have an additional vector in the null space. This is consistent with the theorems from Ref. Farzanehpour and Tokatly (2012). The interacting extension of the present theorem does not seem to hold although we found that the many-body non-interacting ground state gives rise to $\mathbf{K}$ with the same inertia as $\mathbf{H}^{(\lambda_{0})}$ consistent with the previous findings Farzanehpour and Tokatly (2012).

VI Conclusions

In this article, we have described our software to implement potential inversion for a given time trace of quantum simulation data. Given the time-dependent density and an appropriate initial state, our method solves for the scalar Kohn-Sham potential. In this article, we have discussed the implementation details, a heuristic for preparing initial multi-electrons states with non-zero momentum, and preliminary results for one- and two-electron examples.

Next steps for extending our implementation are performing comparisons and combinations with the implicit fix-point methods Ruggenthaler and van Leeuwen (2011); Ruggenthaler et al. (2012); Nielsen et al. (2013) and methods using Kohn-Sham vector potentials Hodgson et al. (2013, 2016). We would also like to test our initial state preparation methods for more multi-electron systems.

An interesting open question concerns the possible formation of shockwaves within our two-electron example. Our method fails to continue the solution past the formation of a sharp feature which we suspect may be a cusp. Further investigations are warranted using tools from the partial differential equations literature. In our present work, we are unable to observe non-analytic spatial behavior due to the lattice formulation. Thus, an approach directly in in the continuum limit formulations could help address analytic behavior of this particular model.

We are also interested in connecting this work to the on-going excitement in the area of quantum computing. Previously, we have shown that the $V$-representation problem can be solved efficiently using a quantum computer but it still remains an open question if all quantum computations remains efficiently simulatable with TDDFT given access to efficient solutions to the $V$-representation problem. Our next steps will also include using simulation data extracted from cloud-accessible quantum computers.

Acknowledgements

This work was primarily supported by the U.S. Department of Energy, Office of Science, Office of Advanced Scientific Computing Research, under the Quantum Computing Application Teams program (Award 1979657). JDW was also supported by Walter and Constance Burke award, the NSF (Grant 1820747) and additional funding from the DOE (Award A053685).