Computational complexity of time-dependent density functional theory

Time-dependent density functional theory (TDDFT) is rapidly emerging as a premier method for solving dynamical many-body problems in physics and chemistry. The mathematical foundations of TDDFT are established through the formal existence of a fictitious non-interacting system (known as the Kohn-Sham system), which can reproduce the one-electron reduced probability density of the actual system. We build upon these works and show that on the interior of the domain of existence, the Kohn-Sham system can be efficiently obtained given the time-dependent density. Since a quantum computer can efficiently produce such time-dependent densities, we present a polynomial time quantum algorithm to generate the time-dependent Kohn-Sham potential with controllable error bounds. As a consequence, in contrast to the known intractability result for ground state density functional theory (DFT), the computation of the necessary time-dependent potentials given the initial state is in the complexity class described by bounded error quantum computation in polynomial time (BQP).

Despite the many successes achieved so far, the major challenge of time-dependent density functional theory (TDDFT) is to find good approximations to the Kohn-Sham potential, V KS , for a non-interacting system.This is a notoriously difficult problem and leads to failures of TDDFT in situations involving charge-transfer excitations [1], conical intersections [2] or photoionization [3].Naturally, this raises the following question: what is the complexity of generating of the necessary potentials?We answer this question and show that access to a universal quantum computer is sufficient.
The present work, in addition to contributing to on-going research about the foundations of TDDFT, is the latest application of quantum computational complexity theory to a growing list of problems in the physics and chemistry community [4].Our result emphasizes that the foundations of TDDFT are not devoid of computational considerations, even theoretically.Further, our work highlights the utility of reasoning using hypothetical quantum computers to classify the computational complexity of problems.The practical implications are that, within the interior of the domain of existence, it is efficient to compute the necessary potentials using a computer with access to an oracle capable of polynomial-time quantum computation.
Quantum computers are devices which use quantum systems themselves to store and process data.On the one hand, one of the selling points of quantum computation is to have efficient algorithms for calculations in quantum chemistry and quantum physics [5,6,7].On the other hand, in the worst case, quantum computers are not expected to solve all NP (non-deterministic polynomial time) problems efficiently [8].Therefore, it is an on-going investigation into when a quantum computer would be more useful than a classical computer.Our current result points towards evidence of computational differences between quantum computers and classical computers.In this way, we provide additional insights to one of the driving questions of information and communication processing in the past decades concerning practical application areas of quantum computing.
Our findings are in contrast to a previous result by Schuch and Verstraete [9], which showed that, in the worst-case, polynomial approximation to the universal functional of ground state density functional theory (DFT) is likely to be impossible even with a quantum computer.
Remarkably, this discrepancy between the computational difficulty of TDDFT and ground state DFT is often reversed in practice where for common place systems encountered by physicists and chemists, TDDFT calculations are often more challenging than DFT calculations.Therefore, our findings provide more reasons why quantum computers should be built.
The practical utility of our results can be understood in multiple ways.First, we have demonstrated a new theoretical understanding of TDDFT highlighting its relative simplicity as compared to ground state DFT computations.Second, we have introduced a V -representability parameter, which similar to the condition number of a matrix, diverges as the Kohn-Sham formalism becomes less applicable.Finally, for analysis purposes, it is often useful to know what the exact Kohn-Sham potential looks like in order to compare and contrast approximations to the exchange-correlation functionals.However, this has been limited to small dimensional or model systems and our results show that, with a quantum computer, one could perform such exploratory studies for larger systems.

Time-dependent Kohn-Sham systems
To introduce TDDFT and its Kohn-Sham formalism, it is instructive to view the Schrödinger equation as a map [10] { V (t), Ψ(t 0 )} → {n(t), Ψ(t)}. ( The inputs to the map are an initial state of N electrons, Ψ(t = t 0 ), and a Hamiltonian, Ĥ(t) = T + Ŵ + V (t) that contains a kinetic-energy term, T , a two-body interaction term such as the Coulomb potential, Ŵ , and a scalar time-dependent potential, V (t).The outputs of the map are the state at later time, Ψ(t) and the one-particle probability density normalized to N (referred to as the density), TDDFT is predicated on the use of the time-dependent density as the fundamental variable and all observables and properties are functionals of the density.The crux of the theoretical foundations of TDDFT is an inverse map which has as inputs the density at all times and the initial state.It outputs the potential and the wave function at later times t, This mapping exists via the Runge-Gross theorem [11] which shows that, apart from a gauge degree of freedom represented by spatially homogeneous variations, the potential is bijectively related to the density.However, the problem of timedependent simulation has not been simplified; the dimension of the Hilbert space scales exponentially with the number of electrons due to the two-body interaction Ŵ .As a result, the time-dependent Schrödinger equation quickly becomes intractable to solve with controlled precision on a classical computer.Practical computational approaches to TDDFT rely on constructing the noninteracting time-dependent Kohn-Sham potential.If at time t the density of a system described by potential and wave function, { V (t), Ψ(t)}, is n Ψ(t) , then the noninteracting Kohn-Sham system ( Ŵ = 0) reproduces the same density but using a different potential, V KS .The key difficulty of TDDFT is obtaining the timedependent Kohn-Sham potential.
Typically, the Kohn-Sham potential is broken into three parts: The first potential is the external potential given in the problem specification and the second is the Hartree potential The third is the exchange-correlation potential and requires an approximation to be specified wherein lies the difficulty of the Kohn-Sham scheme.In this article, we discuss how difficult approximating the full potential is but we make note that only the exchangecorrelation is unknown.While we discuss the computation of the full Kohn-Sham potential from a given external potential and initial density, we will not construct an explicit functional for the exchange-correlation potential.
The route to obtaining the Kohn-Sham potentials we focus on is the evaluation of the map, Here, the wave function of the Kohn-Sham system, Φ(t) = A[φ 1 (t)φ 2 (t)...φ N (t)], is an anti-symmetric combination of single particle wave functions, φ i (t), such that for all times t, the Kohn-Sham density, n KS (t) = n Φ(t) = N i=1 |φ i (t)| 2 , matches the interacting density n Ψ(t) .If such a map exists, we call the system V -representable while implicitly referring to non-interacting V KS -representablity.
As the map in Eq. ( 4) is foundational for TDDFT implementations based on the Kohn-Sham system, there are many articles [12,13,14,15,16,17] examining the existence of such a map.Instead of attempting to merely prove the existence of the Kohn-Sham potential, we will explore the limits on the efficient computation of this map and go beyond the scope of the previous works by addressing questions from the vantage of computational complexity.
The first approach to the Kohn-Sham inverse map found in Eq. ( 4), was due to van Leeuwen [12] who constructed a Taylor expansion in t of the Kohn-Sham potential to prove its existence.The construction relied on the continuity equation, −∇ • ĵ = ∂ t n, and the Heisenberg equation of motion for the density operator to derive the local force balance equation at a given time t: where Q = i[ T , ∂ t n] is the momentum-stress tensor.In the past few years, several results have appeared extending van Leeuwen's construction [13,14,15,16,17] to avoid technical problems (related to convergence and analyticity requirements).Here previous rigorous results by Farzanehpour and Tokatly [17] on lattice TDDFT are directly applicable to our quantum computational setting.

The discrete force balance equation
We summarize the details of the discretized local force-balance equation from [17].More detailed derivations are found in [17] and as well as a more general derivation we provide in Appendix A. Consider a system discretized on a lattice of M points forming a Fock space.In second quantization, the creation âi and annihilation â † j operators for arbitrary sites i and j must satisfy âi âj = −â j âi and âi â † j = δ ij − â † j âi .We define a discretized onebody operator as Â = M n M m A mn â † m ân and designate A as the coefficient matrix of the operator.The matrix elements are A mn = m| Â|n where |m and |n are the single electron sites corresponding to operators âm and ân .Similar notation and definitions hold for the two-body operators.
The Hamiltonian, the density at site j, and the continuity equation are then given respectively by For the density of the Kohn-Sham system, n KS (t) = n Φ(t) , to match the density of the interacting system, n(t) = n Ψ(t) , the discretized local force balance equation [17] must be satisfied, Here Γij = â † i âj + â † j âi is twice the real part of the one-body reduced density operator.A complete derivation of this equation is found the Appendix A. The vector S aim is defined as . The force balance coefficient matrix, K = K Φ(t) , is defined through Eq. ( 10) and Eq.(11).Since the target density enters only through the second derivative appearing in S aim , the initial state Φ(t 0 ) must reproduce the initial density, n Ψ(t0) , and the initial time-derivative of the density, The system is non-interacting V -representable so long as K is invertible on the domain of spatial inhomogeneous potentials.Moreover, the Kohn-Sham potential is unique [17].
Hence, the domain of To ensure efficiency, we must further restrict attention to the interior of this domain where K is sufficiently well-condition with respect to matrix inversion.The cost of the algorithm grows exponentially as one approaches this boundary but can in some cases be mitigated by increasing the number of lattice points.

Quantum algorithm for the Kohn-Sham potentials
We consider an algorithm to compute the density with error in the 1-norm to be efficient when the temporal computational cost grows no more than polynomially in 1/ , polynomially in (max 0<s<t H(s) )t, polynomially in M , the number of sites, and polynomially in, N , the number of electrons.We will describe such an algorithm within the interior of the domain of V -representability.
To ensure that the algorithm is efficient, we must assume that the local kinetic energy and the local potential energy are both bounded by constant E L and that there is a fixed number, κ such that Note that, as we work in the Fock space, this condition does not preclude Coulombic interactions with nuclei so long as the site orbitals have finite spatial extent.
We will show that as long as E L ≤ √ log N , the algorithm remains efficient for fixed κ.As is typical in numerical matrix analysis [18,19], the inversion of a matrix become extremely sensitive to errors as the condition number, C = K K −1 , grows.The Lipschitz constant of the Kohn-Sham potential must also scale polynomially with the number of electrons.
The Lipschitz constant of the Kohn-Sham system could be different than that of the interacting system [20,10] and understanding of the relationship between these timescales requires a better understanding of the initial state Φ(t 0 ) dependence.What can be done, in practice, is to begin with an estimate of the maximum Lipschitz Force2balance+equa/on+ Kohn2Sham+equa/on+ constant and if any two consecutive Kohn-Sham potentials violate this bound, restart with a larger Lipschitz constant.Our efficient algorithm for computing the time-dependent potential, is depicted in Figure 1.There are two stages.The first stage involves a quantum computer and its inputs are the initial many-body state Ψ(t 0 ) and the external potential V (t) on a given interval [t 0 , t 1 ].The quantum computer then evolves the initial state with the given external potential and obtains the time-evolved wave function at a series of discrete time-steps.The detailed analysis of the expectation estimation algorithm found in Ref. [21] is used to bound errors in the measurement of the density and to estimate its second time derivative.In order to rigorously bound the error term, we assume that the fourth time derivative of the density is bounded by a constant, c 4 .
The total cost of both stages of the algorithm is dominated by the cost of obtaining the wave function as this is the only step that depends directly on the number of electrons.Fortunately, quantum computers can perform time-dependent simulation efficiently [22,23,24].The cost depends on the requested error in the wave function, δ ψ , and depends on the length of time propagated when time is measured relative to the norm of the Hamiltonian being simulated.The essential idea is to leverage the evolution of a controllable system (the quantum computer) with an imposed (simulation) Hamiltonian [6].It should be highlighted that obtaining the density through experimental spectroscopic means is equivalent to the quantum computation provided the necessary criteria for efficiency and accuracy are satisfied.
The second stage involves only a classical computer, with the inputs being a consistent initial Kohn-Sham state Φ(t 0 ) and the interacting ∂ 2 t n Ψ(t) on the given interval [t 0 , t 1 ].The output is the Kohn-Sham potential at sufficiently many time steps to ensure the target accuracy is achieved.The classical algorithm performs matrix inversion of a M by M matrix.The cost for the matrix inversion is O(M 3 ) regardless of the other problem parameters (such as the number of electrons).
In our analysis detailed in the next section, we only consider errors from the quantum and classical aspects of our algorithm and we avoided some unnecessary complications by omitting detailed analysis of the classical problem of propagating the non-interacting Kohn-Sham system.Kohn-Sham propagation in the classical computer is well studied and can be done efficiently using various methods [25].Further, we have also assumed that errors in the measured data are large enough that issues of machine precision do not enter.Thus, we have ignored the device dependent issue of machine precision in our analysis and refer to standard treatments [18,19] for the proper handling of this issue.

Overview of error bounds
We demonstrate that our algorithm has the desired scaling by bounding the final error in the density.We follow an explicit-type marching process to obtain the solution at time q∆t from the solution at (q − 1)∆t.The full technique is elaborated in the next section.
As the classical matrix inversion algorithm at each time step is independent of the number of electrons and the quantum algorithm requires poly(N, t 1 − t 0 , δ −1 ψ , −1 ) per time step (recall that δ ψ is the allowed error in the wave function due to the quantum simulation algorithm), we can utilize error analysis for matrix inversion and an explicit marching process to get a final estimate of the classical and quantum costs for the desired precision The parameter r is the number of repetitions of the quantum measurement required to obtain a suitably large confidence interval.We define the V -representability parameter as R = κE 2 L and if R is bounded by a constant, then the algorithm is efficient.The intractability of the algorithm with growing R indicates the breakdown of V -representability.Despite the exponential dependence of the algorithm on the representability parameter, the domain of V -representability is known to encompass all time-analytic Kohn-Sham potentials in the continuum limit [13,14,15,16].Examining the exponential dependence, it is clear that increases in κ can be offset by decreases in the local energy.

Description of techniques used to bound cost
Before diving into the details, let us give an overview of our techniques and what is to follow.In the first subsection, we look at the error in the wave function at time t.In each time step, the error is bounded from the errors in the previous steps.This leads to a recursion relation which we solve to get a bound for the total error at any time step.This error is propagated forward because we must solve KV = S = Q + ∂ 2 t n for V based on the data from the previous time step.The error in ∂ 2 t n is due to the finite precision of the quantum computation and is independent of previous times.In the second subsection, the error in the density is then derived followed by a cost analysis in the final subsection.
We rescale time by factor c such t 1 − t 0 = 1 to get the final time step z = 1/∆t.This rescaling is possible because there is no preferred units of time.That said the rescaling of time cannot be done indefinitely for two reasons.First, the Lipschitz constant of both the real and the KS system must be rescaled by same factor of c.Since the cost of the algorithm depends on the Lipschitz constant, increasingly long times will require more resources.Second, the quantum simulation algorithm does have an intrinsic time scale set by the norm of the H and its time derivatives [22,23,24].Rescaling time by c increases the norm of H by the same factor; consequently, the difficulty of the quantum simulation is invariant to trivial rescaling of the dynamics.
It is important to get estimates which do not directly depend on the number of sites.To do this, we assume that the lattice is locally connected under the hopping term such that there are at most d elements per row of T (since T is symmetric, it is also d-col-sparse).This is equivalent to a bound for the local kinetic energy.
Throughout, we work with the matrix representations of the operators and the states.The L p vector norms [18] with p = 1, 2, and ∞ are defined by |x| p = ( |x i | p ) 1/p .The induced matrix norms are defined by A p = max |x|p=1 |Ax| p .Induced norms are important because they are compatible with the vector norm such that |M x| p = M p |x| p .The vector 1-norm is appropriate for probability distributions and the vector 2-norm is appropriate for wave functions.The matrix 2-norm is also called the spectral norm and is equal to the maximum absolute value of an eigenvalue.For a diagonal matrix, D, the matrix 2-norm is the vector ∞-norm of diag(D).Note that |x| p ≥ |x| p for p < p .Important, non-trivial characterizations of the infinity norms are |x| ∞ = max i |x i | and A ∞ = max i j |A ij |.

Error in the wave function via recursion relations
We bound the error of the evolution operator from time k∆t to (k − 1)∆t, denoted ∆U (k, k − 1) 2 , in terms of the previous time step in order to obtain a recursion relation.We first bound the errors in the potential due to the time discretization and then those due to the computation errors using Lemma 1 found in Appendix B. The computation errors will depend on the error at the previous time step which will lead to the recursion relation sought after.
To bound the error in ∆U 2 we must bound the error in the potential Here, {V (t k )} is the discretized potential with time step |t j − t j+1 | = ∆t.The error due to temporal discretization can be controlled assuming a Lipschitz constant L for the potential such that for all t and t , The computational error Now we need to bound the errors in |∆Q| ∞ and ∆K ∞ in terms of the error The error bound for |∆Q| ∞ is obtained as ≤ max i pq The product d max |T ij | is the maximum local kinetic energy and is, by assumption, bounded by E L .Similarly, = max We convert from errors in the real part of the 1-RDM to errors in the wave function via The inequality Eq. ( 21) follows because the maximum eigenvalue of a † i a j ψ for all ψ is bounded by 1 and Γ ij = 2 real a † i a j ψ .Taking the maximum over all i, j we have Here δ Φ k−1 bounds the error in the two-norm |∆Φ| 2 at time step k − 1.
Putting together Eq. ( 15), Eq. ( 17), Eq. ( 19), and Eq. ( 22) gives To obtain the desired recursion relation, we note that at time step k the error can be bounded via (24) obtained using an expansion similar to the one found in Eq. (20).Utilizing Lemma 1 (see Appendix B) and bound Eq. ( 23), we arrive at To obtain a recursion relation we let the LHS of Eq. ( 25) define the new upper bound at time step k.

Recursion relations of the form f
. Thus, we have for the bound at time step k Now consider the final time step at z = 1/∆t, and e x ≥ (xz −1 + 1) z for z < ∞, We applied Lemma 3 from Appendix B to obtain the last line.This bound is similar to the Euler formula for the global error but arises from the iterative dependence of the potential on the previous error; not from any approximate solution to an ordinary differential equation.
To ensure that the cost is polynomial in M and N for fixed κ, we must insist that E L ≤ √ log N .Consider the exponential factor and assume that E L > 1.Then exp(16κE 2 L ) ≤ exp(16κ log N ) = N 16κ is a polynomial for fixed κ.

Error bound on the density
To finish the derivation, we utilize our bound for the wave function at the final time to get a bound on the error of the density at the final time.This will translate into conditions for the number of steps needed and the precision required for the density.The error in the density is bounded by the error in the wave function through the following, Now consider the i-th element, n i = a † i a i , and the Cauchy- Finally, from the definition of the 1-norm, For final error in the 1-norm of the density, we allow error /2 due to the time step error and /2 error due to the density measurement.Following Eq. ( 29) and Eq. ( 30), we have for the number of time steps, The bound for the measurement precision also follows as,

Cost analysis
To obtain the cost for the quantum simulation and the subsequent measurement, we leverage detailed analysis of the expectation estimation algorithm [21].To measure the density at time t ∈ [t 0 , t 1 ], a quantum simulation [22,23,24] of ψ(t 0 ) → ψ(t) is performed at cost q ≤ poly(N, t 1 − t 0 , δ −1 ψ ) following an assumption that H(t) is simulatable on a quantum computer which is usually the case for physical systems.In order to simplify the analysis, we assume that δ ψ is such that δ n + δ ψ ≈ δ n is a reasonable approximation.Given the recent algorithm for logarithmically small errors [24], this assumption is reasonable.
The expectation estimation algorithm (EEA) was analyzed in [21].The algorithm EEA(ψ, A, δ, c) measures ψ|A|ψ with precision δ and confidence c such that Prob(ã− δ ≤ ψ|A|ψ ≤ ã + δ) > c , that is, the probability that the measured value ã is within δ of ψ|A|ψ is bounded from below by c.The idea is to use an approximate Taylor expansion: The confidence interval is improved by repeating the protocol r = | log(1−c)| times.If the spectrum of A is bounded by 1, then the algorithm requires on the order O(r/δ 3/2 ) copies of ψ and O(r/δ 3/2 ) uses of exp(−iAs) with s = √ 3δ/2.To perform the measurement of the density, we assume that the wave function is represented in first quantization [6] such that the necessary evolution operator is: exp(−in j s) = N k exp(−i|j j| (k) t).Here each Hamiltonian |j j| (k) acts on site j of the kth electron simulation grid.Hence, each operation is local with disjoint support.Since there are N M sites, this can be done efficiently.Comparing the costs, we will assume that the generation of the state dominates the cost.
Combining these facts, we arrive at the conclusion that the cost to measure the density to within δ n precision is Pairing this with Eq. ( 31) and Eq. ( 32), we have an estimate for the number of quantum operations The classical computational algorithm is an [M × M ] matrix inversion at each time step costing

Quantum computation and the computational complexity of TDDFT
Since the cost of both the quantum and classical algorithms scale as a polynomial of the input parameters, we can say that this is an efficient quantum algorithm for computing the time-dependent Kohn-Sham potential.Therefore, the computation of the Kohn-Sham potential is in the complexity class described by bounded error quantum computers running in polynomial time (BQP).This is the class of problems that can be solved efficiently on a quantum computer.Quantum computers have long been considered as a tool for simulating quantum physics [26,27,5,6,7].The applications of quantum simulation fall into two broad categories: (1) dynamics [28,29,30] and (2) ground state properties [31,32,33].The first problem is in the spirit of the original proposal by Feynman [26] and is the focus of the current work.
Unfortunately, unlike classical simulations, the final wave function of a quantum simulation cannot be readily extracted due to the exponentially large size of the simulated Hilbert space.The retrieval of the full state would require quantum state tomography, which in the worst case, requires an exponential number of copies of the state and would take an exponentially large amount of space to even store the data classically.If, instead, the simulation results can be encoded into a minimal set of information and the simulation algorithm can be efficiently executed on a quantum computer, then the problem is in the complexity class BQP.Extraction of the density [21] is the relevant example of such a quantity that can be obtained.Note that the density's time-evolution is dictated by wave function and hence the Schrödinger equation.
In summary, what we have proven is that computing the Kohn-Sham potential at bounded κE 2 L is in the complexity class BQP.To be precise, two technical comments are in order.First, we point out that we are really focused on promise problems since we require constraints on the inputs to be satisfied (i.e.κE 2 L <constant).Second, computing the map Eq. ( 4) is not a decision problem and cannot technically be in the complexity class BQP.However, we can define the map to b bits of precision by solving M log b accept-reject instances from the corresponding decision problem, which is in BQP.These concepts are further elaborated in [34,35,4].While the quantum computer would allow most dynamical quantities to be extracted without resorting to the Kohn-Sham formalism, we have attempted to understand the difficulty of generating the Kohn-Sham potential.We only consider a polynomial time quantum computer as a tool for reasoning about the complexity of computing Kohn-Sham potentials.In essence, the Kohn-Sham potentials are a compressed classically tractable encoding of the quantum dynamics that allows the quantum simulation to be performed in polynomial time on a classical computer.This may have implications for the question of whether a classical witness can be used in place of quantum witness in the quantum Merlin Arthur game [35] (i.e.QMA ?=QCMA).A second useful by-product of our result is the introduction of the Vrepresentability parameter which has general significance for practical computational settings.

Concluding remarks
In this article, we introduced a V -representability parameter and have rigorously demonstrated two fundamental results concerning the computational complexity of time dependent density functional theory with bounded representability parameter.First, we showed that with a quantum computer, one need only provide the initial state and external potential on the interval [t 0 , t 1 ] in order to generate the timedependent Kohn-Sham potentials.Second, we show that if one provides the density on the interval [t 0 , t 1 ], the Kohn-Sham potential can be obtained efficiently with a classical computer.
We point out that an alternative to our lattice approach may exist using tools from partial differential equations.Early results in this direction have been pioneered using an iterated map whose domain of convergence defines V -representability [15,16].The convergence properties of the map have been studied in several one-dimensional numerical examples [36,15,16].Analytical understanding of the rate of convergence to the fixed point would complement the present work with an alternate formulation directly in real space.
While this paper focuses on the simulation of quantum dynamics, the complexity of the ground state problem is interesting in its own right [34,35,4,9].In this context, ground state DFT was formally shown [9] to be difficult even with polynomial time quantum computation.Interestingly, in that work, the Levy-minimization procedure [37] was utilized for the interacting system to avoid discussing the non-interacting ground state Kohn-Sham system and its existence.We have worked within the Kohn-Sham picture, but it may be interesting to construct a functional approach directly.
Future research involves improving the scaling with the condition number or showing that our observed exponential dependence on the representability parameter is optimal.Our work can likely be extended to bosonic and spin systems [38] since we have relied minimally on the fermionic properties of electrons.Finally, preconditioning the matrix K can also help increase the domain of computationally feasible V -representability.
Our findings provide further illustration of how the fields of quantum computing and quantum information can contribute to our understanding of physical systems through the examination of quantum complexity theory.
So now consider the LHS as vector Ŝ with components Ŝj = ∂ 2 t nj − Qj .Similarly consider the potential V as a vector with components V i , then we can write Eq. (A.11) as Ŝ = KV .Examining Eq. (A.10), if V k = V k for all k, k then the RHS of Eq. (A.10) vanishes.Hence, K always has at least one vector in the null space, namely the spatially constant potential.
Farzanehpour and Tokatly [17] study the existence of a unique solution for the non-linear Schrödinger equation which follows from Eq. (A.11): (A.12) In the space where K has only one zero eigenvalue, the Picard-Lindelöf theorem [M.E. Lindelöf, C. R. Hebd.Sances Acad.Sci.116, 454 (1894)] guarantees the existence of a unique solution.
The Picard-Lindelöf theorem concerns the differential equation ∂ t y(t) = f (t, y(t)) with initial value y(t 0 ) on t ∈ [t 0 − ε, t 0 + ε].If f is bounded above by a constant and is continuous in t and Lipschitz continuous in y then, according to the theorem, for ε > 0, there exists a unique solution y(t) on [t 0 − ε, t 0 + ε].This solution can be extended until either y becomes unbounded or y is no longer a solution.The conditions of the theorem are satisfied because K(Φ) and Ŝ are quadratic in Φ, the RHS is Lipschitz continuous in Φ in the domain where K has only one zero eigenvalue, and the continuity of K and Ŝ in time follows immediately from the continuity of Φ.
A nice connection of Eq. (A.11) to master equations in probabilistic processes can be drawn.In Eq. (A.11), K has the form of a master equation for a probability distribution P , The key difference is that the entries of K are not strictly positive ( Φ(t)|â † i âj |Φ(t) can be positive or negative).Since K is Hermitian and its null space contains the uniform state, if all transition coefficients were positive, then K would satisfy detailed balance.
Appendix B. Lemmas Lemma 1.For two time-dependent Hamiltonians H(t) = H 0 + V (t) and H(t) = H 0 + Ṽ (t), the error in the evolution from t 0 to t 1 is bounded as To obtain the statement in Eq. (B.1), recall that for a diagonal matrix, the induced matrix 2-norm is the infinity norm of the corresponding vector of diagonal elements.Noting that V is diagonal gives V 2 = |V | ∞ to complete the proof.

Figure 1 .
Figure1.In part a, the quantum computer takes as inputs the initial state and the time-dependent Hamiltonian and outputs the density at sufficiently many times.The output allows the numerical computation of the second derivative of the density at each time step which is then utilized by the classical computer to solve the discrete force balance equation Eq.(11).A consistent initial state at time t = 0 must also be given which reproduces n(0) and ∂tn(0).Note that while the wave function is obtained from the quantum computation, it cannot be processed for use in the classical part of the computation.The classical algorithm uses the density to obtain the Kohn-Sham potential at each subsequent time step through an iterated marching process as depicted in part b.