# Feynman’s clock, a new variational principle, and parallel-in-time quantum dynamics

###### Abstract

We introduce a new discrete-time variational principle inspired by the quantum clock originally proposed by Feynman, and use it to write down quantum evolution as a ground state eigenvalue problem. The construction allows one to apply ground state quantum many-body theory to quantum dynamics, extending the reach of many highly developed tools from this fertile research area. Moreover this formalism naturally leads to an algorithm to parallelize quantum simulation over time. We draw an explicit connection between previously known time-dependent variational principles and the new time embedded variational principle presented. Sample calculations are presented applying the idea to a Hydrogen molecule and the spin degrees of freedom of a model inorganic compound demonstrating the parallel speedup of our method as well as its flexibility in applying ground-state methodologies. Finally, we take advantage of the unique perspective of the new variational principle to examine the error of basis approximations in quantum dynamics.

### .1 Introduction

Feynman proposed a revolutionary solution to the problem of quantum simulation: use quantum computers to simulate quantum systems. While this strategy is powerful and elegant, universal quantum computers may not be available for some time, and in fact, accurate quantum simulations may be required for their eventual construction. In this work, we will use the clock Hamiltonian originally introduced by FeynmanFeynman (1982, 1985) for the purposes of quantum computation, to re-write the quantum dynamics problem as a ground state eigenvalue problem. We then generalize this approach, and obtain a novel variational principle for the dynamics of a quantum system and show how it allows for a natural formulation of a parallel-in-time quantum dynamics algorithm. Variational principles play a central role in the development and study of quantum dynamicsLippmann and Schwinger (1950); Heller (1976); Kerman and Koonin (1976); Jackiw and Kerman (1979); Balian and Vénéroni (1981); Deumens et al. (1994); Poulsen (2011); Haegeman et al. (2011), and the variational principle presented here extends the arsenal of available tools by allowing one to directly apply efficient approximations from the ground state quantum many-body problem to study dynamics.

Following trends in modern computing hardware, the simulations of quantum dynamics on classical hardware must be able to make effective use of parallel processing. We will show below that the perspective of the new variational principle leads naturally to a time-parallelizable formulation of quantum dynamics. Previous approaches for recasting quantum dynamics as a time-independent problem include Floquet theory for periodic potentialsMilfeld and Wyatt (1983); Autler and Townes (1955); Dion and Hirschfelder (2007) and more generally the formalism of Peskin and MoiseyevPeskin and Moiseyev (1993). However, the approach proposed in this manuscript differs considerably from these previous approaches. We derive our result from a different variational principle, and in our embedding the dynamics of the problem are encoded directly in its solution, as opposed to requiring the construction of another propagator. Considerable work has now been done in the migration of knowledge from classical computing to quantum computing and quantum informationNielsen and Chuang (2000); Aspuru-Guzik et al. (2005); Kassal et al. (2008); Wang et al. (2008); Kassal and Aspuru-Guzik (2009). In this paper, we propose a novel use of an idea from quantum computation for the simulation of quantum dynamics.

The paper is organized as follows. We will first review the Feynman clock: a mapping stemming from the field of quantum computation that can be employed for converting problems in quantum evolution into ground-state problems in a larger Hilbert space. We then generalize the Feynman clock into a time-embedded discrete variational principle (TEDVP) which offers additional insight to quantum time-dynamics in a way that is complementary to existing differential variational principles Dirac (1930); Frenkel (1934); McLachlan (1964). We then apply the configuration interaction method Slater (1929); Boys (1950) from quantum chemistry to solve for approximate dynamics of a model spin system demonstrating convergence of accuracy of our proposed approach with level of the truncation. We demonstrate how this construction naturally leads to an algorithm that takes advantage of parallel processing in time, and show that it performs favorably against existing algorithms for this problem. Finally we discuss metrics inspired by our approach that can be used to quantitatively understand the errors resulting from truncating the Hilbert space of many-body quantum dynamics.

### .2 Physical dynamics as a sequence of quantum gates

Consider a quantum system described by a time-dependent wavefunction . The dynamics of this system are determined by a Hermitian Hamiltonian according to the time-dependent Schrödinger equation in atomic units,

(1) |

A formal solution to the above equation can be generally written:

(2) |

Where is the well known time-ordering operator and is a unitary operator that evolves the system from a time to a time . These operators also obey a group composition property, such that if then

(3) |

From the unitarity of these operators, it is of course also true that where indicates the adjoint. Thus far, we have treated time as a continuous variable. However, when one considers numerical calculations on a wavefunction, it is typically necessary to discretize time.

We discretize time by keeping an ancillary quantum system, which can occupy states with integer indices ranging from to where is the number of discrete time steps under consideration. This quantum system has orthonormal states such that

(4) |

This definition allows one to encode the entire evolution of a physical system by entangling the physical wavefunction with this auxiliary quantum system representing time, known as the “time register”. We define this compound state to be the history state, given by

(5) |

where subscripts are used to emphasize when we are considering a time-independent state of a system at time . That is, we define . From these definitions, it is immediately clear from above that the wavefunction at any time can be recovered by projection with the time register, such that

(6) |

Additionally, we discretize the action of our unitary operators, such that and we embed this operator into the composite system-time Hilbert space as . While the utility of this discretization has not yet been made apparent, we will now use this discretization to transform the quantum dynamics problem into a ground state eigenvalue problem.

### .3 Feynman’s clock

In the gate modelNielsen and Chuang (2000); Farhi et al. (2000) of quantum computation, one begins with an initial quantum state state, applies a sequence of unitary operators , known as quantum gates. By making a measurement on the final state , one determines the result of the computation, or equivalently the result of applying the sequence of unitary operators . The map from the sequence of unitary operators in the gate model, to a Hamiltonian that has the clock state as its lowest eigenvector is given by a construction called the Clock HamiltonianFeynman (1985). Since its initial inception, much work has been done on the specific form of the clock, making it amenable to implementation on quantum computersKitaev et al. (2002). However, for the purposes of our discussion that pertain to implementation on a classical computer, the following simple construction suffices

(7) |

where is a penalty term which can be used to specify the state of the physical system at any time. Typically, we use this to enforce the initial state, such that if the known state at time is given by , then

(8) |

One may verify by action of on the history state defined above, , that the history state is an eigenvector of this operator, with eigenvalue .

### .4 A discrete-time variational principle

We now introduce a new time-embedded discrete variational principle (TEDVP) and show how the above eigenvalue problem emerges as the special case of choosing a linear variational space. Suppose that one knows the exact form of the evolution operators and wavefunctions at times . By the properties of unitary evolution it is clear that the following holds:

(9) |

From this, we can see that if the exact construction of is known for all , but the wavefunctions are only approximately known(but still normalized), it is true that

(10) |

where equality holds in the case that the wavefunction becomes exact at each discrete time point. Reintroducing the ancillary time-register, we may equivalently say that all valid time evolutions embedded into the composite system-time space as minimize the quantity

(11) |

where (script font for operators denotes they act in the composite system-time space) is the operator given by

(12) |

It is then clear from the usual ground state variational principle, that the expectation value of the operator

(13) |

is only minimized for exact evolutions of the physical system. This leads us immediately to a time-dependent variational principle for the discrete representation of a wavefunction given by:

(14) |

It is interesting to note, that this is a true variational principle in the sense that an exact quantum evolution is found at a minimum, rather than a stationary point as in some variational principlesMcLachlan (1964). This is related to the connection between this variational principle and the McLachlan variational principle that is explored in the next section. However, to the authors knowledge, this connection has never been explicitly made until now.

To complete the connection to the clock mapping given above, we first note that this operator is Hermitian by construction and choose a linear variational space that spans the entire physical domain. To constrain the solution to have unit norm, we introduce the Lagrange multiplier and minimize the auxiliary functional given by

(15) |

It is clear that this problem is equivalent to the exact eigenvalue problem on with eigenvalue . The optimal trajectory is given by the ground state eigenvector of the operator . From this construction, we see that the clock mapping originally proposed by Feynman is easily recovered as the optimal variation of the TEDVP in a complete linear basis, under the constraint of unit norm. Note that the inclusion of as a penalty term to break the degeneracy of the ground state is only a convenient construction for the eigenvalue problem. In the general TEDVP, one need not include this penalty explicitly, as degenerate allowed paths are excluded, as in other time-dependent variational principles, by fixing the initial state.

We note, as in the case of the time-independent variational principle and differential formulations of the time-dependent variational principle, the most compact solutions may be derived from variational spaces that have non-linear parameterizations. Key examples of this in chemistry include Hartree-Fock, coupled cluster theoryCoester and Kümmel (1960); Cizek (1966); Bartlett and Musial (2007), and multi-configurational time-dependent HartreeBeck et al. (2000). It is here that the generality of this new variational principle allows one to explore new solutions to the dynamics of the path without the restriction of writing the problem as an eigenvalue problem, as in the original clock construction of Feynman.

### .5 Connection to previous time-dependent variational principles

In the limit of an exact solution, it must be true that all valid time-dependent variational principles are satisfied. For that reason, it is important to draw the formal connection between our variational principle and previously known variational principles.

Considering only two adjacent times and , the operator is given by:

(16) |

Now suppose that the separation of physical times between discrete step and , denoted is small, and the physical system has an associated Hamiltonian given by , such that

(17) |

By inserting this propagator into equation 14, rewriting the result in terms of , and dropping terms of order we have

(18) |

After defining the difference operator such that , we can factorize the above expression into

(19) |

In the limit that , these difference operators become derivatives. Defining , this is precisely the McLachlan variational principleMcLachlan (1964).

(20) |

We then conclude that in the limit of infinitesimal physical time for a single time step, the TEDVP is equivalent to the McLachlan variational principle. Under the reasonable assumptions that the variational spaces of the wavefunction and its time derivative are the same and that the parameters are complex analytic, then it is also equivalent to the Dirac-Frenkel and Lagrange variational principlesBroeckhove et al. (1988). Moreover, as supplementary material (SI1), we provide a connection that allows other variational principles to be written as eigenvalue problems, and further discuss the merits of the integrated formalism used here.

To conclude this section, we highlight one additional difference between practical uses of the TEDVP and other variational principles: The TEDVP is independent of the method used to construct the operator . In quantum information applications, this implies it is not required to know a set of generating Hamiltonians for quantum gates. Additionally, in numerical applications, one is not restricted by the choice of approximate propagator used. In cases where an analytic propagator is known for the chosen basis, it can be sampled explicitly.Suppose that the dimension of the physical system is given by and the number of timesteps of interest is given by . Assuming that the time register is ordered, the resulting eigenvalue problem is block tridiagonal with dimension (See Fig. 1). This structure has been described at least once before in the context of ground-state quantum computationMizel (2004), but to the authors knowledge, never in the context of conventional simulation of quantum systems.

## I Many-Body application of the TEDVP

There has been a recent rise in the interest of methods for simulating quantum spin dynamics in chemistry Kuprov et al. (2007); Hogben et al. (2010). To study the properties of the clock mapping when used to formulate approximate dynamics, we chose a simple model spin system inside an inorganic molecule Luban et al. (2002). Specifically, we examine the spin dynamics of the vanadium compound depicted in Fig. 2. By choosing the three unpaired electron spins to interact with one another by means of isotropic exchange as well as uniform static external magnetic field , and a time-dependent transverse field , this system can be modeled with a spin Hamiltonian

(21) |

Where is the Pauli operator on spin , , is the Bohr magneton, is the measured spectroscopic splitting factor. The couplings as well as are fitted through experimental determinations of magnetic susceptibilityLuban et al. (2002). The fact that they are not equal is reflective of the isosceles geometry of the vanadium centers. The parameters of this model, are given by: , , and .
We will allow to vary to study the properties of the clock mapping in the
solution of approximate quantum dynamics.

The quantum chemistry community has decades of experience in developing and employing methods for obtaining approximate solutions of high-dimensional, ground-state eigenvector problems. By utilizing the connection we have made from dynamics to ground state problem, we will now borrow and apply the most conceptually simple approximation from quantum chemistry: configuration interaction in the space of trajectoriesHabershon (2012), to our model problem to elucidate the properties of the clock mapping.

For the uncorrelated reference, we take the entire path of a mean-field spin evolution that is governed by the time-dependent Hartree equations, and write it as:

(22) |

where is the reference spin-down state for spin , is the reference spin-down state after rotation at time , is the whole product system at time , and is determined from the mean field Hamiltonian. The equations of motion that determine are derived in a manner analogous to the time-dependent Hartree equations, and if one writes the wavefunction in the physical space as

(23) |

Then the equations of motion are given by

(24) | |||

(25) |

Where is the mean field Hamiltonian for spin formed by contracting the Hamiltonian over all other spins , is the expectation value of the Hamiltonian at time , and is the number of spins in the system.

To generate configurations, we also introduce the transformed spin excitation operator , which is defined by

(26) | ||||

(27) | ||||

(28) |

In analogy to traditional configuration interaction, we will define different levels of excitation (e.g. singles, doubles, …) as the number of spin excitations at each time that will be included in the basis in which the problem is diagonalized. For example, the basis for the configuration interaction singles(CIS) problem is defined as

(29) |

Note that for is simply defined to be the identity operator on all spins so that the reference configuration is naturally included. Similarly, the basis for single and double excitations(CISD) is given by

(30) |

Higher levels of excitation follow naturally, and it is clear that when one reaches a level of excitation equivalent to the number of spins, this method may be regarded as full configuration interaction, or exact diagonalization in the space of discretized paths.

The choice of a time-dependent reference allows the reference configuration to be nearly exact when , independent of the nature of the time-dependent transverse field. This allows for the separation of the consequences of time-dependence from the effects of two versus one body spin interactions.

After a choice of orthonormal basis, the dynamics of the physical system are given by the ground state eigenvector of the projected eigenvalue problem

(31) |

where we explicitly define the projection operator onto the basis as so that the projected operator is given by

Using these constructions, we calculate the quantum dynamics of the sample system described above. For convenience, we rescale the Hamiltonian by the value of . To add arbitrary non-trivial time dependence to the system and mimic the interaction of the system with a transverse pulse, we take . The magnitude of was taken to be in order to model perturbative two body interactions in this Hamiltonian. To propagate the equations of motion and generate the propagators for the clock operator we use the enforced time-reversal symmetry exponential propagatorCastro et al. (2004) given by

(32) |

The dynamics of some physical observables are displayed (Fig. 5) for the reference configuration, single excitations, double excitations, and full configuration interaction. The physical observables have been calculated with normalization at each time step. It is seen (Fig. 3) that as in the case of ground state electronic structure the physical observables become more accurate both qualitatively and quantitatively as the configuration space expands, converging to the exact solution with full configuration interaction. Moreover in Fig. 4 we plot the contributions from the reference configuration, singles space, doubles space, and triples space and observe rapidly diminishing contributions. This suggests that the time-dependent path reference used here provides an good qualitative description of the system. As a result, perturbative and other dynamically correlated methods from quantum chemistry may also be amenable to the solution of this problem.

In principle, approximate dynamics derived from this variational principle are not norm conserving, as is seen in Fig. 4, however this actually offers an important insight into a major source of error in quantum dynamics simulations of many-body systems, which is the truncation of the basis set as described in the last section. Simulations based on conventional variational principles that propagate within an incomplete configuration space easily preserve norm; however the trajectories of probability which should have left the simulated space are necessarily in error.

## Ii Parallel-in-time quantum dynamics

Algorithms that divide a problem up in the time domain, as opposed to spatial domain, are known as parallel-in-time algorithms. Compared to their spatial counterparts, such as traditional domain decompositionSmith et al. (2004), these algorithms have received relatively little attention. This is perhaps due to the counterintuitive notion of solving for future times in parallel with the present. However as modern computational architectures continue to evolve towards many-core setups, exploiting all avenues of parallel speedup available will be an issue of increasing importance. The most popular parallel-in-time algorithm in common use today is likely the parareal algorithmLions et al. (2001); Baffico et al. (2002). The essential ingredients of parareal are the use of a coarse propagator that performs an approximate time evolution in serial, and a fine propagator that refines the answer and may be applied in parallel. The two propagations are combined with a predictor-corrector scheme. It has been shown to be successful with parabolic type equationsGander and Vandewalle (2007), such as the heat equation, but it has found limited success with wave-like equationsGander (2008), like the time-dependent Schrödinger equation. We will now show how our variational principle can be used to naturally formulate a parallel-in-time algorithm, and demonstrate its improved domain of convergence with respect to the parareal algorithm for Schrödinger evolution of Hydrogen.

Starting from the TEDVP, minimization under the constraint that the initial state is fixed yields a Hermitian positive-definite, block-tridiagonal linear equation of the form

(33) |

where (to be solved for) contains the full evolution of the system and specifies the initial condition such that

(34) |

The linear clock operator, , is similar to before, where now we distinguish between a clock built from a coarse operator, , from a clock built from a fine operator, , as

(35) |

The spectrum of this operator is positive-definite and admits distinct eigenvalues, each of which is -fold degenerate. The conjugate gradient algorithm can be used to solve for , converging at-worst in a number of steps which is equal to the number of distinct eigenvaluesHestenes and Stiefel (1952); Faber and Manteuffel (1984). Thus application of the conjugate gradient algorithm to this problem is able to converge requiring at most applications of the linear clock operator . This approach on its own yields no parallel speedup. However, the use of a well-chosen preconditioner can greatly accelerate the convergence of the conjugate gradient algorithmEisenstat (1981).

If one uses an approximate propagation performed in serial, , which is much cheaper to perform than the exact propagation, as a preconditioning step to the conjugate gradient solve, the algorithm can converge in far fewer steps than and a parallel-in-time speedup can be achieved. The problem being solved in this case for the clock construction is given by

(36) |

To clarify and compare with existing methods, we now introduce an example from chemistry. The nuclear quantum dynamics of the Hydrogen molecule in its ground electronic state can be modeled by the Hamiltonian

(37) |

where , , and atomic unitsMakri (1993). The initial state of our system is a gaussian wavepacket with a width corresponding to the ground state of the harmonic approximation to this potential, displaced Å from the equilibrium position. To avoid the storage associated with the propagator of this system and mimic the performance of our algorithm on a larger system, we use the symmetric matrix-free split-operator Fourier transform method(SOFT) to construct block propagatorsFeit et al. (1982). This method is unconditionally stable, accurate to third order in the time step , and may be written as

(38) |

Here, and corresponds to the application of the fast Fourier transform (FFT) and its inverse over the wavefunction grid. The use of the FFT allows each of the exponential operators to be applied trivially to their eigenbasis and as a result the application of the propagator has a cost dominated by the FFT that scales as , where is the number of grid-points being used to represent . For our algorithm, we define a fine propagator, , and a coarse propagator, from the SOFT construction, such that for a given number of sub-time steps .

(39) | ||||

(40) |

We take for our problem the clock constructed from the fine-propagator and use the solution of the problem built from the coarse propagator as our preconditioner. In all cases, only the matrix free version of the propagator is used, including in the explicit solution of the coarse propagation.

From the construction of the coarse and fine propagators, with processors, up to communication time, the cost of applying the fine clock in parallel and solving the coarse clock in serial require roughly the same amount of computational time. This is a good approximation in the usual case where the application of the propagators is far more expensive than the communication required. From this, assuming the use of processors, we define an estimated parallel-in-time speed-up for the computational procedure given by

(41) |

where is the number of applications of the fine-propagator matrix performed in parallel and is the number of serial linear solves using the coarse-propagator matrix used in preconditioned conjugate gradient. The factor of 2 originates from the requirement of backwards evolution not present in a standard propagation. The cost of communication overhead as well as the inner-products in the CG algorithm are neglected for this approximate study, assuming they are negligible with respect to the cost of applying the propagator.

The equivalent parallel speedup for the parareal algorithm is given by

(42) |

where and are now defined for the corresponding parareal operators which are functionally identical to the clock operators without backward time evolution, and thus it lacks the same factor of 2.

As is stated above, in the solution of the linear clock without preconditioning, it is possible to converge the problem in at most steps, independent of both the choice of physical timestep and the size of the physical system by using a conjugate gradient method. However, with the addition of the preconditioner, the choice of timestep and total time simulated can have an effect on the total preconditioned system. This is because as the coarse (approximate) propagation deviates more severely from the exact solution, the preconditioning of the problem can become quite poor.

This problem exhibits itself in a more extreme way for the parareal algorithm, as the predictor-corrector scheme may start to diverge for very poorly preconditioned system. This has been seen before in the study of hyperbolic equationsGander (2008), and remains true in this case for the evolution of the Schrödinger equation. The construction derived from the clock is more robust and is able to achieve parallel-in-time speedup for significantly longer times. This marks an improvement over the current standard for parallel-in-time evolution of the Schrödinger equation.

To give a concrete example, consider the case where we divide the evolution of the nuclear dynamics of hydrogen into pieces containing evolution blocks, each of which is constructed from fine evolutions for a time as is described above. The dynamics over which we simulate are depicted in Fig. 8. We average the estimated parallel speedup for 10 time blocks (which we define as the whole time in one construction of the clock) forward and the results are for the speedup are given in Fig. 7. In this example we see that for small (or small total evolution times), the reduced overhead of having no backwards evolution yields an advantage for the parareal algorithm. However as the increases, the parareal algorithm is less robust to errors in the coarse propagation and performance begins to degrade rapidly. In these cases, our clock construction demonstrates a clear advantage. It is a topic of current research how to facilitate even longer evolutions in the clock construction.

## Iii Norm Loss as a measure of truncation error

Conservation of the norm of a wavefunction is often considered a critical property for approximate quantum dynamics, as it is a property of the exact propagator resulting from time-reversal symmetry. However, if norm conservation is enforced in a propagation which does not span the complete Hilbert-space, one must account for the components of the wavefunction that would have evolved into the space not included in the computation. It’s not immediately clear how population density which attempts to leave the selected subspace should be folded back into the system without being able to simulate the exact dynamics. This problem is sometimes glossed over with the use of exponential propagators that are guaranteed to produce norm-conserving dynamics on a projected Hamiltonian. Some more sophisticated approaches adjust the configuration space in an attempt to mitigate the errorWestermann and Manthe (2012).

This discrepancy is at the center of the difference between the approximate dynamics derived from the discrete variational principle here and the approximate dynamics derived from the McLachlan variational principle such as the multi-configurational time-dependent Hartree method. Mathematically, this results from the non-commutativity of the exponential map and projection operator defined above. That is for a Hermitian operator . In an approximate method derived from the McLachlan or any of the other differential time-dependent variational principles, the projection is performed on the Hamiltonian. As the projection of any Hermitian operator yields another Hermitian operator, the dynamics generated from the projection are guaranteed to be unitary if a sufficiently accurate exponential propagator is used. Conversely, projection of a unitary operator, as prescribed by the TEDVP, does not always yield a unitary operator. Thus for an approximate basis, one expects norm conservation to be violated, where the degree of violation is related to the severity of the configuration space truncation error. This leads us to define a metric of truncation error which we term the instantaneous norm loss. We define this as

(43) |

where is always assumed to be normalized, which in practice means that a renormalization is used after each time step here. However, as we proved above, in the limit of a short time step, with dynamics generated by a Hamiltonian, the TEDVP must contain essentially the same content as the McLachlan variational principle. For this reason, we propose an additional metric which is given by

(44) |

Where is the physical Hamiltonian. This is motivated by appealing to the McLachlan variational principle and substituting from the exact Schrödinger equation that where is the full (non-projected) Hamiltonian. By defining , we recognize this as a pertubation theory estimate of the error resulting from the configuration basis truncation at a given point in time.

To examine the quality of these metrics and to better understand the consequences of the non-commutativity of the exponential map and projection, we return to the sample spin system. In this case, we choose a basis for the propagation space based entirely on the initial state, and do not allow it to change dynamically in time as before. We perform simulations in the space of single excitations (S) from the initial state, double excitations (SD), and in the full Hilbert space (Ex). Dynamics from the TEDVP are generated by first building the exact propagator then projecting to the desired basis set while dynamics from the McLachlan variational principle (MVP) are modeled by projecting the Hamiltonian into the target basis set and exponentiating. A renormalization is used after each time step in the first case. Although one could perform the simulation with a timestep where timestep error is negligible, we remove this component of the calculation for this example by making the Hamiltonian time-independent. This allows direct analysis of the effect of timestep on non-commutativity deviations.

In Fig. 6 we show the dynamics of the Vanadium spin complex for the two approximate truncation levels (S and SD) with both methods and their associated error metrics( and ). Deviations in the qualitative features of the observable occur after even the first peak of the proposed metrics. In this particular simulation, the configuration interaction with singles and doubles spans all but one state in the full Hilbert space. The lack of this one state results in the large qualitative errors present, associated with the first and subsequent peaks present in these error metrics. The impact of later peaks is more difficult to discern, due to error in the wavefunctions, which accumulates as the propagation proceeds.

As predicted by the connection between the TEDVP and the McLachlan variational principle, while and are not identical for each case, in the short time limit they yield extremely similar information, which is highlighted in Fig. 6 displaying the resulting longer time dynamics for a time step of . In Fig. 9, however we explore the effects of a significantly larger time step and begin to discern the result of the non-commutativity discussed here. Recalling that because the Hamiltonian is time-independent in this case, the propagator used is numerically exact in both instances, so this effect is not a result of what would be traditionally called time step error, resulting from intrinsic errors of an integrator. Interestingly, it is observed that begins to decay to a nearly constant value. This occurs because the action of projection after exponentiation breaks the degeneracy of the spectrum of the unitary operator, resulting in eigenvalues with norms different than . As a result, repeated action and renormalization of the operator is analogous to a power method for finding the eigenvector associated with the largest eigenvalue. This effect is exacerbated by taking long time steps.

## Iv Conclusions

In this manuscript, we introduce a new variational principle for time-dependent dynamics inspired by the Feynman clock originally employed for quantum computation. Unlike other previously-proposed variational principles, the proposed TEDVP approach involves the solution of an eigenvalue problem for the entire time propagation. This perspective allows for readily employing many of the powerful truncation techniques from quantum chemistry and condensed-matter physics, that have been developed for the exact diagonalization problem. We show how this formulation naturally leads to a parallel-in-time algorithm and demonstrate its improved robustness with respect to existing methods. We introduce two novel error metrics for the TEDVP that allow one to characterize the basis approximations involved. The features of the method were demonstrated by simulating the dynamics of a Hydrogen molecule and a molecular effective spin Hamiltonian. Further research directions include the use of other approximate techniques for the time dynamics such as the use of perturbation theory Helgaker et al. (2002) or coupled-cluster approachesKvaal (2012), and further enhancement of parallel-in-time dynamics.

### iv.1 Acknowledgments

We would like to acknowledge David Tempel for his valuable comments on the manuscript. J.M. is supported by the DOE Computational Science Graduate Fellowship under grant number DE-FG02-97ER25308. A.A.-G. and J.P. thank the National Science Foundation for their support under award number CHE-1152291. A.A.-G. receives support from the HRL Laboratories, LLC under award number M1144-201167-DS and UCSD under grant number FA9550-12-1-0046. Research sponsored by United States Department of Defense. The views and conclusions contained in this document are those of the authors and should not be interpreted as representing the official policies, either expressly or implied, of the U.S. Government. We also thank the Corning Foundation, the Camille and Henry Dreyfus Foundation, and the Alfred P. Sloan Foundation for their continuing support of this research.

### iv.2 Appendix: On the general construction of eigenvalue problems from dynamics

We have presented one path for constructing eigenvalue problems from quantum dynamics problems so far, however it is instructive to illuminate precisely which part of our procedure allowed this. To do this, we will slightly generalize the Floquet-type Hamiltonians and demonstrate that the time embedding was the crucial feature that allows one to recast a dynamics problem as an eigenvalue problem. This is in addition to the choice to work in an integrated framework, which we will show simply allows for a convenient choice of basis.

Recall the definition of a Floquet-type Hamiltonian given by

(45) |

If one considers a finite time evolution for a length of time , this operator is Hermitian in the basis of Fourier functions states given by

(46) |

where is a time-independent state of the physical system when considering the generalized inner-product

(47) |

This is the the generalized Hilbert space first introduced by Sambe Sambe (1973) and generalized by Howland Howland (1974). Because this operator is Hermitian in this basis, by noting the similarity to the Lagrange Variational principle

(48) |

minimization of on this linear basis of Fourier states yields a Hermitian eigenvalue problem. Thus the time evolution can be reconstructed by solving the full time-independent eigenvalue problem in this basis, or by constructing a surrogate evolution operator as in the formalism of Peskin and Moiseyev Peskin and Moiseyev (1993). The use of Fourier basis states to express time dependence is natural given the form of the operator . That is, matrix elements of the derivative operator have a trivial analytic expression in this basis. This represents the same general idea we have been discussing here, which is to consider states in a system-time Hilbert space. However, as the solution of this variational problem will may yield a stationary point rather than a true minimum McLachlan (1964), ground state techniques are not appropriate for this particular operator. Moreover, this operator is not in general Hermitian when considering arbitrary basis functions of time. For example arbitrary choices of plane waves not corresponding to the traditional Fourier basis will yield a non-Hermitian operator. The operator , which has been used in the past for the construction of approximate dynamics Shalashilin and Burghardt (2008), is Hermitian, however it still suffers from a pathology that the optimal solution represents a stationary point rather than a minimum. However, the operator

(49) |

is always Hermitian and positive semi-definite by construction. Thus one can expand the system-time Hilbert space in any linear basis of time, and the optimal path in that space will be the ground state eigenvector of the operator utilizing the above generalized inner-product, assuming we have broken degeneracy by introducing the correct initial state. This can be viewed as an application of the McLachlan time dependent variational principle. From this, we see it is the expression systems in a system-time Hilbert space which allows for the eigenvalue problem construction. Moreover, we note that this is not limited to the use of Fourier time basis states, and that many more expressions of time dependence may be utilized to construct an eigenvalue problem within this framework.

One may ask then, what was the objective of working in the an integrated formalism, defined using unitary operators rather than differential operators. To see this, consider evaluating the system at a point in time with the Fourier construction above. One has to expand the system in all time basis functions and evaluate them at a time . When one usually considers time, however, they are thinking of a parametric construction where a number simply labels a specific state of a system. Embedding into the system-time Hilbert space with this idea would be most naturally expressed as delta-functions. However matrix elements of on this basis can be difficult to construct. In the ancillary time system framework used in the TEDVP, however, time is easily expressed as a discrete parametric variable. One might also consider the use of discretized derivatives in the operator . However, balancing the errors in numerical derivatives and the increasing difficultly of solving the problem with the number of simultaneously stored time steps can be practically difficult.

## References

- Feynman (1982) R. Feynman, Int. J. Theor. Phys. 21, 467 (1982).
- Feynman (1985) R. P. Feynman, Opt. News 11, 11 (1985).
- Lippmann and Schwinger (1950) B. A. Lippmann and J. Schwinger, Phys. Rev. 79, 469 (1950).
- Heller (1976) E. J. Heller, J. Chem. Phys. 64, 63 (1976).
- Kerman and Koonin (1976) A. Kerman and S. Koonin, Annals of Physics 100, 332 (1976).
- Jackiw and Kerman (1979) R. Jackiw and A. Kerman, Phys. Lett. A 71, 158 (1979).
- Balian and Vénéroni (1981) R. Balian and M. Vénéroni, Phys. Rev. Lett. 47, 1353 (1981).
- Deumens et al. (1994) E. Deumens, A. Diz, R. Longo, and Y. Öhrn, Rev. Mod. Phys. 66, 917 (1994).
- Poulsen (2011) J. A. Poulsen, J. Chem. Phys. 134, 034118 (2011).
- Haegeman et al. (2011) J. Haegeman, J. I. Cirac, T. J. Osborne, I. Pižorn, H. Verschelde, and F. Verstraete, Phys. Rev. Lett. 107, 070601 (2011).
- Milfeld and Wyatt (1983) K. F. Milfeld and R. E. Wyatt, Phys. Rev. A 27, 72 (1983).
- Autler and Townes (1955) S. H. Autler and C. H. Townes, Phys. Rev. 100, 703 (1955).
- Dion and Hirschfelder (2007) D. R. Dion and J. O. Hirschfelder, “Time-dependent perturbation of a two-state quantum system by a sinusoidal field,” in Adv. Chem. Phys. (John Wiley and Sons, Inc., 2007) pp. 265–350.
- Peskin and Moiseyev (1993) U. Peskin and N. Moiseyev, J. Chem. Phys. 99, 4590 (1993).
- Nielsen and Chuang (2000) M. Nielsen and I. Chuang, Quantum Computation and Quantum Information, Cambridge Series on Information and the Natural Sciences (Cambridge University Press, 2000).
- Aspuru-Guzik et al. (2005) A. Aspuru-Guzik, A. D. Dutoi, P. J. Love, and M. Head-Gordon, Science 309, 1704 (2005).
- Kassal et al. (2008) I. Kassal, S. P. Jordan, P. J. Love, M. Mohseni, and A. Aspuru-Guzik, Proc. Natl. Acad. Sci. USA 105, 18681 (2008).
- Wang et al. (2008) H. Wang, S. Kais, A. Aspuru-Guzik, and M. R. Hoffmann, Phys. Chem. Chem. Phys. 10, 5388 (2008).
- Kassal and Aspuru-Guzik (2009) I. Kassal and A. Aspuru-Guzik, J. Chem. Phys. 131, 224102 (2009).
- Dirac (1930) P. A. M. Dirac, Math. Proc. Cambridge 26, 376 (1930).
- Frenkel (1934) J. Frenkel, Wave Mechanics (Claredon Press, Oxford, 1934).
- McLachlan (1964) A. McLachlan, Mol. Phys. 8, 39 (1964).
- Slater (1929) J. C. Slater, Phys. Rev. 34, 1293 (1929).
- Boys (1950) S. F. Boys, Proc. R. Soc. Lond. A Math. Phys. Sci. 200, 542 (1950).
- Farhi et al. (2000) E. Farhi, J. Goldstone, S. Gutmann, and M. Sipser, (2000), arXiv:quant-ph/0001106 .
- Kitaev et al. (2002) A. Kitaev, A. Shen, M. Vyalyi, and N. Vyalyi, Classical and Quantum Computation, Graduate Studies in Mathematics (American Mathematical Society, 2002).
- Coester and Kümmel (1960) F. Coester and H. Kümmel, Nucl. Phys. 17, 477 (1960).
- Cizek (1966) J. Cizek, J. Chem. Phys. 45, 4256 (1966).
- Bartlett and Musial (2007) R. J. Bartlett and M. Musial, Rev. Mod. Phys. 79, 291 (2007).
- Beck et al. (2000) M. Beck, A. Jäckle, G. Worth, and H.-D. Meyer, Phys. Rep. 324, 1 (2000).
- Broeckhove et al. (1988) J. Broeckhove, L. Lathouwers, E. Kesteloot, and P. V. Leuven, Chem. Phys. Lett. 149, 547 (1988).
- Mizel (2004) A. Mizel, Phys. Rev. A 70, 012304 (2004).
- Kuprov et al. (2007) I. Kuprov, N. Wagner-Rundell, and P. Hore, J. Magn. Reson. 189, 241 (2007).
- Hogben et al. (2010) H. J. Hogben, P. J. Hore, and I. Kuprov, J. Chem. Phys. 132, 174101 (2010).
- Luban et al. (2002) M. Luban, F. Borsa, S. Bud’ko, P. Canfield, S. Jun, J. K. Jung, P. Kögerler, D. Mentrup, A. Müller, R. Modler, D. Procissi, B. J. Suh, and M. Torikachvili, Phys. Rev. B 66, 054407 (2002).
- Habershon (2012) S. Habershon, J. Chem. Phys. 136, 054109 (2012).
- Castro et al. (2004) A. Castro, M. A. L. Marques, and A. Rubio, J. Chem. Phys. 121, 3425 (2004).
- Smith et al. (2004) B. Smith, P. Bjorstad, and W. Gropp, Domain decomposition (Cambridge University Press, 2004).
- Lions et al. (2001) J. Lions, Y. Maday, and G. Turinici, Comptes Rendus de l’Academie des Sciences Series I Mathematics 332, 661 (2001).
- Baffico et al. (2002) L. Baffico, S. Bernard, Y. Maday, G. Turinici, and G. Zérah, Phys. Rev. E 66, 057701 (2002).
- Gander and Vandewalle (2007) M. J. Gander and S. Vandewalle, SIAM J. Sci. Comp. 29, 556 (2007).
- Gander (2008) M. J. Gander, Boletín SEMA (2008).
- Hestenes and Stiefel (1952) M. R. Hestenes and E. Stiefel, “Methods of conjugate gradients for solving linear systems,” (1952).
- Faber and Manteuffel (1984) V. Faber and T. Manteuffel, SIAM J. Numer. Anal. 21, 352 (1984).
- Eisenstat (1981) S. C. Eisenstat, SIAM J. Sci. Stat. Comp. 2, 1 (1981).
- Makri (1993) N. Makri, J. Phys. Chem. 97, 2417 (1993).
- Feit et al. (1982) M. Feit, J. F. Jr., and A. Steiger, J. Comput. Phys. 47, 412 (1982).
- Westermann and Manthe (2012) T. Westermann and U. Manthe, J. Chem. Phys. 137, 22A509 (2012).
- Helgaker et al. (2002) T. Helgaker, P. Jorgensen, and J. Olsen, Molecular Electronic Structure Theory (Wiley, Sussex, 2002).
- Kvaal (2012) S. Kvaal, J. Chem. Phys. 136, 194109 (2012).
- Sambe (1973) H. Sambe, Phys. Rev. A 7, 2203 (1973).
- Howland (1974) J. S. Howland, Mathematische Annalen 207, 315 (1974).
- Shalashilin and Burghardt (2008) D. V. Shalashilin and I. Burghardt, J. Chem. Phys. 129, 084104 (2008).