Mixed Quantum-Classical Dynamics for Near Term Quantum Computers

Mixed quantum-classical dynamics is a set of methods often used to understand systems too complex to treat fully quantum mechanically. Many techniques exist for full quantum mechanical evolution on quantum computers, but mixed quantum-classical dynamics are less explored. We present a modular algorithm for general mixed quantum-classical dynamics where the quantum subsystem is coupled with the classical subsystem. We test it on a modified Shin-Metiu model in the first quantization through Ehrenfest propagation. We find that the Time-Dependent Variational Time Propagation algorithm performs well for short-time evolutions and retains qualitative results for longer-time evolutions.


I. INTRODUCTION
Quantum computers have found great success in electronic structure theory through the variational quantum eigensolver [1] and subsequent algorithms known as variational quantum algorithms (VQAs) [2].Adding additional electrons to a system greatly increases its complexity, this is something we hope quantum computers could handle.If one were to also consider the full nuclear dynamics, the problem becomes unmanageable much faster, potentially even for quantum computers [3].One way to reconcile this is to partition the system into interacting quantum and classical parts.This is the realm of mixed quantum-classical (MQC) approaches, which are a widely used set of tools for understanding chemical systems [4,5].In quantum computing, this area is less researched than the electronic structure problem, but it is actively being explored [3,6,7].In this work we propose and explore a noisy intermediate-scale quantum (NISQ) friendly algorithm that can be used to study MQC dynamics.
Using quantum computers alongside classical computers is the backbone of VQAs, but splitting a system into sections treated separately by each machine is not new.A DFT embedding scheme with a quantum computer expansion of the active space [8] and has been found to outperform certain types of state-of-the-art approximate techniques such as CASSCF [9] in finding ground state energies.Furthermore, ground state dynamics, geometry relaxation, and force measurements for MD applications have been explored with success in [7].In [3], dynamics are explored in both first and second quantization, but using a time-independent Hamiltonian.This is also the case for various other time propagation techniques [10][11][12]; we draw inspiration from and use p-VQD [11] in this work.
Our contribution is the presentation of a general algorithmic structure to tackle non-adiabatic molecular dynamics (NAMD) by offloading the QM part to a quantum computer and evolving the classical system by the Ehrenfest method.Observables from the quantum mechanical (QM) subsystem are measured and used to update the classical system, which in turn will update the time-dependent Hamiltonian that is used to evolve the QM state in turn.This is all done in first quantization, which saves the algorithm from needing to measure nonadiabatic couplings, as these are treated directly within the wave function and its evolution in this setting.The algorithm is demonstrated in the Shin-Metiu model [13], which is often used to test various non-adiabatic techniques [14][15][16][17].We modify this to be a NAMD-like problem by partitioning the system into a classical nucleus and quantum electron.The major contribution is the study of how the interaction between observable measurement and system updates play out as well as introducing a scheme that is suited to begin exploring other time-dependent phenomena on NISQ machines.
The theoretical advantage of using quantum computers is that they have access to an exponentially growing computational space for each additional qubit in the system [18].Current machines have access to hundreds of qubits, which would ideally allow them to already outperform current supercomputers.This is not the case due to noise coming from interactions with the environment and imperfect gate implementations.As such, NISQ algorithms [19] have to contend with limits on the number of imperfect operations that can be made.But even if this were not the case and full quantum dynamics could be simulated, we probably will always want to tackle a problem bigger than current machines can handle, so these kinds of approximations will always be used.
It should be noted that existing supercomputers by far outperform existing quantum computers in handling large quantum-chemistry problems, and applications to chemical problems will have to be deferred until there is a provable quantum advantage.Approximations like limiting the simulation to a selected active space [8] can extend the reach of quantum computers, but analogues of these ideas apply to classical computers as well.For now, quantum algorithms in chemistry mostly study systems that are comfortably computable on current classical hardware.

A. Time-Dependent Hamiltonian Variational
Quantum Propagation The time-dependent variational quantum propagation (TDVQP) algorithm builds on the circuit compression idea of "projected variational quantum dynamics" (p-VQD) [11] by allowing the Hamiltonian to be timedependent.For many large problems of interest to theoretical chemistry, especially in MD, it is impossible to fully simulate the system of interest quantum mechanically.As such, the system is subdivided into classical and quantum components.The evolution of both systems occurs in locked steps, with the classical system defining the Hamiltonian for the quantum evolution, and the quantum system then feeding back into the classical system in the way of some observable, usually the energy gradient (force).One mustn't limit themselves to molecular or even physical systems, as this algorithm would work with any set of observables that can be used to update the classical system of interest.
The algorithm begins with a parameterized circuit initialized to some desired state.This is denoted as |ψ 0 ⟩, which will have been generated according to a Hamiltonian based on an initial vector of classical parameters ⃗ q 0 of the classical coordinates, which we denote H(⃗ q 0 ) := H 0 .This is done by choosing some sufficiently expressive parameterized circuit ansatz Ĉ( ⃗ θ) which takes the quantum computer's initial state, denoted |0⟩, to |ψ 0 ⟩.This can be done using a VQA to find a chosen state with respect to H 0 , which returns the circuit parameters 0 }.These observables are used to evolve the classical state of the system, generating a new vector of classical parameters ⃗ q 1 .These can then be used to generate H 1 and { Ô(s) 1 }.Now, one evolves the state from |ψ 0 ⟩ to |ψ 1 ⟩ by applying the time evolution operator to the state, The physical implementation of the time evolution can either be the Trotterized form of the operator or some other approximate time evolution.In this work, we simply use Ĥ0 to evolve the state with a first-order trotter expansion.The Hamiltonian used for the time evolution could be of higher order.For example, ( Ĥ0 + Ĥ1 )/2 can be used with no extra cost in this scheme, but higherorder integrators will require an additional evolution and observable measurement for each timestep beyond H 1 .
In return, one gets higher-order symplectic integration.Now a single step of the p-VQD algorithm is applied which generates the new circuit parameters θ 1 such that |ψ 1 ⟩ ≈ Ĉ(θ 1 ) |0⟩ to some desired threshold.This process is repeated until the desired timestep is reached.The entire process is more precisely described in Algorithm 1, and a depiction of the quantum circuit can be seen in Fig. 1.The overall number of circuit evaluations is linear with respect to the number of iterations, circuit parameters, timesteps and Pauli terms of the observables.This is treated in greater depth in S-I A.
Algorithm 1: Time-Dependent Variational Quantum Propagation (TDVQP) 14 Function TDVQP(Hgen(⃗ q ), C( ⃗ θ), ⃗ q0): i }, ⃗ qi in arrays θ, O, q θ, O, q θ, O, q TDVQP should be thought of as a meta-algorithm that has replaceable components.The most directly replaceable part is the choice of ansatz Ĉ( ⃗ θ), which at the moment is generally a heuristic choice for most problems in NISQ devices.More advanced ansatze such as the family of adaptive ansatze, which changes the ansatz throughout the evolution would work, but could not use the previous step's θ parameters as effectively.The very costly time evolution is currently a Trotterized form of the time evolution operator, as in this work, and in [11,12].This can be replaced by a plethora of more NISQ-friendly time evolutions as is done in [20,21] if the form of the Hamiltonian allows this.The limit is the no-fast forwarding theorem [22][23][24], which states that you cannot achieve a time evolution of time t in a sublinear gate count for a general Hamiltonian, but for shorter time evolutions, limited sizes and specific cases of Hamiltonians, includ- 1. Sketch of the TDVQP process with slices showing the state after each gate at the initial condition and final condition.The initial guess for the parameter vector θ ′ i is θi and its final value is denoted θi+1.Then observables Ô can be measured on |ψi+1⟩ to update the Hamiltonian or the classical state of the system.
ing our sparse Hamiltonian using short time evolutions, this likely is not the case [21].
In the classical evolution, the choice of integrator and the actual Hamiltonian used in the time evolution will depend on the type of problem and desired accuracy.Integrators like the Velocity Verlet algorithm require no additional resources.TDVQP becomes exact when Ĉ( ⃗ θ) can express the system perfectly for any configuration of classical parameters ⃗ q, given that the exact parameters θ can be found by optimization.This is only a statement of the best-case scenario.In reality, finding a good ansatz, VQA and shot-efficient optimizer is at the forefront of research in this area [2], and it is out of scope for this work.

Error propagation in TDVQP
The TDVQP algorithm inherits all of the errors of its constituent parts.This includes the chosen circuit compression algorithm, time evolution approximation, and in the classical propagator.Nonetheless, it is important to have an intuition of the potential pitfalls of the algorithm.This section illustrates the sources of error in the wavefunction and observables and their interaction velocity Verlet integrator.A more thorough derivation and explanation can be found in the supplementary materials, S-I.
When running the algorithm, any coherent error on the wavefunction representation in the quantum computer ψ can be represented as a superposition of the desired state |ψ⟩ and some combination of undesired orthogonal states |ϕ⟩, such that ψ = √ 1 − I 2 |ψ⟩ + I |ϕ⟩, where I is the infidelity.When we measure the expectation value of a Hermitian observable O on ψ we will get How this translates to the actual measured observable used here is completely system dependent.This will lead to an error in the observable, which in the case of the velocity Verlet integrator with a force error F i ϵ ∝ I 2 in 1 dimension will give a new position Ri of shifted from the expected true position R i .This is linear in the error of the force and quadratic with respect to the timestep ∆t.The following time evolution Hamiltonian and observable operator will be based on this position with an error, which is again, system dependent.The effect is illustrated by the equations Even in the one-dimensional model used in this work, this effect is not analytically computable, but it is small if the timesteps are sufficiently small.This then enters the velocity ( Ṙ) update as which is linear in the error and timestep.Assuming a constant error over all time of F ϵ , that is to say, that the force deviates from the correct one by a constant offsetthis is equivalent to having an additional linear term on the potential.This has the overall effect on the position at iteration i of This expression is quadratic in i and quadratic in timestep.The effect on the fidelity of the TDVQP wavefunction compared to an exact propagation is nontrivial, but numerical examples are provided in the supplementary materials S-I.
The other main source of error inherent to the p-VQD algorithm is that the optimizer never finds a perfect representation of the time-evolved wavefunction, but rather an approximation that meets some fidelity threshold T < 1.If this threshold is met exactly at each p-VQD step, assuming all observable measurements are unaffected, then the decrease in the fidelity is modelled by When the algorithm is run under limited quantum resources and thus subject to finite sampling noise, both the observable and p-VQD step fidelity measurements will have some Gaussian distribution, which will feed into the errors above on a simulation by simulation basis.The effect of this has been analysed numerically in the case of our modified Shin-Metiu Model.

B. The Shin Metiu Model
The Shin-Metiu model is a numerically exactly solvable minimal model which captures essential nonadiabatic effects [13].It is often used as a benchmark system for new techniques and is used to study the effects of different environments as has been done for polaritonic dynamics, coupling to cavities, and the effect of electromagnetic fields [14][15][16]25].It is simple to change its parameters for it to exhibit adiabatic to strongly non-adiabatic dynamics.
In its simplest and original conception, the model shown in Fig. 2 consists of two stationary ions separated by a distance of L, specifically located at L 2 and −L 2 .These enclose a mobile ion p of mass M at distance R from the origin and an electron e − at distance r.The modified Coulomb potential is parameterized by the constants R l , R r and R f , as shown in Eq. 8.This is done to avoid singularities and make the system numerically simpler to simulate.
FIG. 2. Illustration of the Shin-Metiu model with fixed ions at −L 2 , L 2 as stationary boundaries, the mobile ion p of mass M at distance R from the origin and the electron e − at distance r from the origin.R l , Rr and R f are constants for the regularized Coulomb potential in eq. 8.
The full Hamiltonian of the system is with the electronic part being The equation uses atomic units, setting e = Z = ℏ = 1, we also take m = 1 and M = 1836 in the simulation.The constants R l , R r and R f , as shown in Figure 2 are chosen to create specific adiabatic surfaces with transitions we would like to observe, as in Figure 3.
We use the values R f = 5.0, R l = 4.0 and R r = 3.2, which resulting in avoided crossing around R = −1.9 when the distance between the ions is L = 19.These parameters were chosen to be similar to those used in several studies of the model [14,17].The shape of the Born-Oppenheimer potential energy surfaces (BOPES) can be seen in Fig.

Ehrenfest propagation of the model
To perform Ehrenfest propagation of the Shin-Metiu model, we split the system in two.The nucleus (p ) and the electron (e − ).The electron subsystem is treated as a quantum particle described by Hamiltonian 8, where H e is parameterized by the nuclear position (R).The nuclear subsystem is treated classically by tracking parameters of position (R) and velocity ( Ṙ).
For initial coordinates R 0 and Ṙ0 , we first prepare the electronic Hamiltonian H e (R 0 ), which is used to compute the initial state of the electron |ψ 0 ⟩.Thanks to the simplicity of the model, we can use exact diagonalization to compute the eigenvectors and choose any arbitrary superposition of eigenvectors as the initial state.
The nucleus is evolved using the velocity Verlet method [26] with the acceleration being computed from the Coulombic repulsion from the fixed ions and the force from the electronic state.The electronic state is evolved by unitary time evolution with the Hamiltonian at the nuclear position.We use a timestep ∆t, and the system state at timestep i is denoted by under scripts i, where the time is simply i • ∆t.We set our initial conditions at timestep 0, and for the i th step, we compute C. Numerical results The results shown in this section are the result of two types of potential cases.The first, which is referred to as 'single' is the evolution of a single set of initial conditions meant to represent the precision of this algorithm to exactly reproduce a quantum-classical system.Although this is not the intended use case of TDVQP, it is nonetheless the most instructive to determine its behaviour.The second, referred to as 'MD' is the molecular dynamics-like use case, we use a single TDVQP evolution per trajectory.The trajectories are picked from random pairs of normal thermal distributions around the same initial state as the 'single' simulations, with the specific values written in Section IV A. Their average behaviour is taken as the approximation to the true system evolution.For NISQ devices there are some potential cases for different approaches to MD such as those mentioned in [27].
FIG. 4. Mean relative TDVQP energy shown for "Single" initialization (blue) and for "MD" initialization (orange).The highlighted areas show the standard deviation of the distribution of 100 separate runs of 1000 timesteps of 0.5 a.u. at the infinite shot limit.The other lines show the ideal initial state energy evolution.The energy continually increases in the TDVQP as higher energy levels are increasingly populated through leakage.
An important gauge for the validity of simulations of closed systems is whether they conserve energy or not.We use a symplectic integrator in the classical system (velocity Verlet), and in the exact diagonalization case, we see energy conservation for up to 50,000 timesteps.As can be seen in Figure 4 the TDVQP algorithm does not conserve energy.This is because the populations are not preserved in the diagonal basis in the p-VQD step, as the optimization is limited to a finite number of iterations and the ansatz is system agnostic.The effect of this can be seen clearly in Figure 5, where it can be seen that the population in higher states increases much faster than in the ideal case.Although it is not shown, starting the exact evolution from the VQE state does begin with some population spread, but this does not change as the evolution progresses.The population plot is shown at the infinite shot limit for clarity, but the finite shot cases can be found in S-I B.
FIG. 5. Time evolution of state populations between the ideal simulation (solid) and 100 instances of TDVQP evolution (dashed) over 1000 timesteps of 0.5 a.u. at the infinite shot limit for the "single" simulation.The populated states are the same as those seen in Figure 3.The main graph is in logarithmic scale showing faint lines for higher energy levels populated by TDVQP, with the inset showing a linear scale of the two most populated levels.
As a consequence of the higher energy levels being increasingly populated as the evolution progresses, it is the case that the fidelity decreases gradually.This is indeed the case and can be seen in Figure 7.This general degradation of quality is not optimal and strategies could be employed in the optimization to mitigate this, such as measuring the energy and allowing the cost function to penalize when the system is not conserving energy.This would require measuring the expectation value of the system Hamiltonian which would increase the cost of this algorithm.
Despite the problem with energy conservation, using such an algorithm to measure an observable such as the force exerted on the nucleus by the electron (F el ) can still lead to reasonable results.Figure 6 shows the mean of the electron force measurements from TDVQP compared to the ideal measurements at different per-circuit shot counts.It is clear that the mean value slowly deviates from the ideal evolution in even the infinite shot limit, and that you require 10 5 shots per circuit to reach qualitatively relevant results at longer times.Efficiently estimating energy gradients is a huge undertaking, and this work does not implement some of the NISQ-friendly techniques that have been developed [28,29], but it is expected to be a problem even in the fault-tolerant regime [30].
We see in Figure 7 that the fidelity decays in all cases over time and that for long time evolutions, one requires more than 10 5 shots when not using any additional tech- FIG. 6. Electron force observable for the single initialization.The ideal simulation (black) and 100 instances of TDVQP evolutions (coloured) over 1000 timesteps of 0.5 a.u. at varying shot counts per circuit.We observe qualitative agreement of TDVQP with the ideal case at 10 5 shots and the infinite shot limit, with a shift downwards due to leakage to higher energy levels, which in this case biases the force in this direction.niques to better measure the force or better preserve the populations when not undergoing a transition.At lower shot counts the fidelity falls quickly, following eq.7 until the equal superposition is approached, which sets a higher floor than zero for the decay of the fidelity.The potential effects of other noise sources are described in more detail in S-I. Figure 8 zooms into the two more reasonable fidelity lines, those of the infinite and 10 jectories in an MD sense somewhat improves the simulation fidelity compared to the single trajectory case, and more quantum-tailored algorithms like [27] may improve this further.In the infinite shot case, the difference is minimal -but due to the larger variance of the MD simulations compared to the ideal trajectory, the performance tends to be minimally worse.The relationship between shots and fidelity is also illustrated in Figure 9, where one can more clearly see that the MD simulation slightly improves the simulation at longer time evolutions when using finite shots.However, this improvement is not massive.It also highlights the large jump in fidelity gained when using higher shot counts.The p-VQD result [11] on which we base our time evolution evolves its system for 40 iterations (20 a.u.here), where we see very high compression fidelities beyond 10 4 shots per circuit evaluation.
Finally, Figure 10 illustrates that in the simulation the maximum number of iterations (100) is quickly reached before the 200 th iteration at the infinite shot limit.The overall mean final infidelity is 10 −5 , although the fitted threshold of eq.7 for the overall algorithm is slightly lower at 3•10 5 .The infidelity is I = 1−F, where F is the fidelity.This implies there is an additional error, likely due to the drift of the exact simulation of the system from the TDVQP simulation.This is consistent with a 10 −6 to 10 −7 shift in the force as described in the numerical simulations in S-I, which is also roughly the difference in the force observable seen in Figure 6.
Overall the results show some interesting behaviour.The number of timesteps modelled in this work is very high, and this results in only qualitative agreement in the region of interest where there is a significant population transfer.Furthermore in this simple model, the energy levels are well separated and we begin in the ground  state.This results in a strong unidirectional contribution from population leakage to higher energy levels.In a more complex molecule, one would begin, for example, from a thermal ensemble of not only velocities but also states.In turn, the leakage would then result in deviations from both higher and lower energy levels and may be less detrimental to the ensemble average than here.

III. DISCUSSION
Time-dependent evolution is an exceptionally interesting problem that can be well explored through quantum computers.Many techniques can be used for full quantum systems [11,12,20,[31][32][33] which are suitable to both near term and fault-tolerant machines.
Algorithms that are suitable for MQC dynamics do require efficient and accurate full quantum dynamics, but the interplay between the classical and quantum systems brings a new spate of challenges.To exchange information, one must measure observables from the quantum system, which is expensive and destroys the state, requiring at minimum an efficient way to measure energy gradients, which is an area of active research [28,29,34].This is a disadvantage, but it also means that one is limited to short-time evolutions between measurements.This makes it possible that a single trotter step is accurate enough [35,36], which is beneficial to near-term devices.
Even though larger timesteps may be possible, the longer the time evolution, the longer the optimizer takes to find the time-evolved ansatz parameters.This is because the previous timestep parameters are no longer as close to the evolved ones.At the same time, most classical MQC methods do not update the parameters that govern the classical system's evolution at the small time intervals we use [4].It would be advantageous to use the largest possible classical timestep for a given integrator.To do this, one could do multiple compression steps with short-time Trotterizations using a constant Hamiltonian and only measuring the desired observables after the quantum system has evolved for the standard timestep of your classical problem, performing updates after this point.This will leverage the underlying compression algorithm to its fullest and reduce the overall number of measurements required.
The algorithm we present takes advantage of the above facts and is highly modular.Although the results are shown using an algorithm like p-VQD [11] with Trotterization of the operator, there is no reason that other efficient time evolution algorithms couldn't be used.This is especially true if the time evolution operator could be efficiently represented by techniques other than the Trotterization of the Hamiltonian.The update step used here measures the Pauli string decomposition of the dH dR matrix to compute forces, but other techniques exist in the fault-tolerant regime [30,34], as well as in the NISQ regime [28,29].The main constraint with TDVQP is the fact that throughout the time evolution, there are inevitable inaccuracies in optimization, due to the compression step not preserving the populations in the diagonal basis as would have been expected as shown in Figure 5.This has the direct consequence that energy is not conserved, even though in the ideal simulation this is the case as shown in Figure 4. Due to this accumulated error, fidelity falls consistently, and the effect is compounded when quantum resources are finite.
These problems might be tackled by either increasing the threshold of the compression step or by measuring the energy and penalizing the optimizer when energy is not conserved.Another option that may be possible is designing an ansatz with problem-specific constraints [37].Such an ansatz considers properties such as particle preservation within their structure, which may remove the need for expensive additional iteration steps or measurements.Furthermore, it may be possible to replace the p-VQD propagation with other compression methods [12].Since we are working in the first quantization representation, it may be difficult to find what properties to conserve in the wavefunction, but with that disadvantage, we gain an advantage in not needing to measure non-adiabatic couplings.
We have also found that the error mostly comes from the compression step or due to finite sampling effects more than from the coupling to the classical system due to the small classical timestep.Although it is always interesting to see how an algorithm behaves under noisy conditions, the performance of p-VQD under noise has been explored for full quantum dynamics in [12].This work focuses on the interplay between the scheme under the effect of a Hamiltonian which depends on the measured observables.
Overall we have introduced the TDVQP algorithm for MQC dynamics with the quantum subsystem computed on a quantum computer and have explored it on the Shin-Metiu model as an example of Ehrenfest dynamics in first quantization.However, it is not limited to this setting.It reproduces the expected observables and state evolution qualitatively.The algorithm is modular and refinements to it may be tackled in future research.Inaccuracies of the quantum computer can also be mitigated when computing ensemble averages of the classical properties.This work shows that MQC simulations may be practically feasible on noisy quantum computers if it is proven that variational quantum algorithms can have an advantage in chemical problems.

IV. METHOD A. Numerical simulations
To gauge the performance of the scheme we implement the Shin-Metiu model as described in Section II B. In the BOPES, we see an avoided crossing at around R = −1.9a.u.We initialize the system with the nucleus at an initial position of R = −2 a.u. and an initial velocity of v 0 = 1.14 • 10 −3 a.u., the average nuclear velocity from the Boltzmann distribution at 300K.The electronic system is initialized through the VQE with a random set of parameters and is allowed 300 iterations to approximate the ground state.The system is then evolved through the TDVQP algorithm with a timestep of ∆t = 0.5 a.u.Each quantum time evolution step attempts to reach a fidelity threshold of 1 − 10 −5 or up to 100 iterations of stochastic gradient descent [38].to find the optimal circuit parameters to approximate the previous time evolved state.Gradients were computed through the parameter shift rule [39].All simulations are done on 16 grid points that can be represented by 4 qubits.
We examine two different situations.First, keep the initial conditions constant but sample different VQE ground state approximations, which we call the "Single Initial Condition" case.In the second case, we examine the MD-type approach in depth, where we sample a normal distribution of initial conditions for the initial velocity of the nucleus and allow one TDVQP evolution per sample.The velocity distribution is sampled from the Boltzmann distribution, only keeping positive velocities so that the nuclei approach the avoided crossing.The results shown are 100 samples that are evolved for 1000 timesteps which bring the classical trajectory beyond the avoided crossing point.
Additional examples are provided for longer-time evolutions as well as for non-ground state evolution in Appendix S-IV.Excited states and superpositions are prepared by using the uncomputation step of the TDVQP, but instead of starting the state with a known circuit from the VQE, the simulator is simply initialized to a desired arbitrary state, and the optimizer attempts to uncompute it with the ansatz and then those parameters are used as the initial step in place of the VQE.Various techniques to prepare excited states exist [40,41], but are not the focus of this work.
We use two different metrics to establish the accuracy of the TDVQP algorithm: the so-called "Ideal" evolution begins at the desired state to numerical precision and is evolved by exact diagonalization.But, precise state preparation is another area of intense study [42].To better gauge the performance of the TDVQP in isolation, we also perform an "Exact" evolution, which uses the VQE-optimized initial state for evolution via exact diagonalization.This allows us to remove any bias from a poorly optimized ground state.
The VQE uses an ansatz of the form shown in Figure 11, which was heuristically chosen as it can achieve ground state infidelities of up to 10 −5 on this system with 4 layers.We use the same ansatz as in [11], but various ansatze can be used, and for first quantization problems, in particular, there are some examples of how several different heuristic ansatze perform in [43].The number of repetitions of the Trotterization layer is another important parameter, but as the decomposition of the Trotterized operator into native gates is deep, we limit ourselves to one.Although for full quantum dynamics, this would be very inaccurate for larger timesteps, the interaction with a classical system necessitates that we use short time steps, so that the Hamiltonian of the system is kept up to date with the classical state of the system.This means that a single trotter step is all that is needed, and the number of layers of the ansatz can compensate as shown in S-II.
The simulations were run on the Qiskit state vector simulator (version 0.28) [44] using the parameter-shift rule [39,45] to determine the analytic gradients required for gradient-descent based optimization.Numpy [46] was used for the exact numerical simulations, to prepare the Hamiltonian and to compute the velocity Verlet steps.

B. Grid-based mapping to qubits
We treat the Shin-Metiu system on the quantum computer in first quantization.We use a finite difference method on an equidistant grid.For low-dimensional problems, this is an appropriate approximation, but in general discrete variable representations (DVR) are a better choice for problems in higher dimensions.In quantum computing DVRs have been used to explore first quantization simulations in [32] using the Colbert and Miller DVR [47].Issues exist with using DVRs on quantum computers as they generally require a full matrix Hamiltonian which is costly to measure and implement on quantum computers, and alternatives have been proposed [43].
Quantum computing in first quantization has the advantage that n g grid points can be represented by N = log 2 (n g ) qubits.We choose n g = 2 N to maximize the use of the N available qubits.In the simplest finite differences method each position is an integer multiple L /ng.The grid point g is by the quantum state |g⟩ and mapped as where j k is the k th bit value of the binary representation of g.The potential operator V is diagonal in this representation, simply sampling the potential at each grid point.The kinetic energy Hamiltonian is not diagonal in the position representation, and although one could use the split operator method [48] to make it diagonal in momentum space would require a quantum Fourier transform implementation, which, as far as we know, cannot be effectively implemented on existing quantum devices.
We sidestep all of the issues by using the finite differences method, in which the one-dimensional potential and the Laplacian form of the kinetic energy can be written as The tridiagonal matrix that results has the same value on the off-diagonal terms, which allows them to be decomposed into fewer Pauli strings than a full matrix via an elegant recursive form that is described in S-III following from [49].This is important because the number of non-zero entries is related to the number of terms in the Pauli decomposition, which should be kept minimal to reduce the length of the Trotterized time evolution operator and observable measurements.The finite difference method does require high grid densities to be accurate (although this depends on how oscillatory the system in question is), but this requirement will likely be met by the doubling of grid points per additional qubit.
Although not implemented in this paper, some algorithms can make the time evolution of Hamiltonians with this form more efficient on near-term devices [50].Depending on the particular Hamiltonian one chooses to study in this way, different efficient algorithms exist to lessen the cost of the time evolution such as variational fast forwarding and qubitization [20,21].It is also possible to efficiently solve for the eigenstates of tridiagonal matrices on quantum computers [51].

C. Circuit compression
A key building block of the presented algorithm and a fundamental aspect of quantum circuit optimization is the concept of circuit compression [52].Any operation on a quantum computer must be a unitary operation Û , but the quantum computer only has a finite set of few-qubit gates.An arbitrary Û must be expressed, or compiled, into a set of native gates [18], this can always be done, but it is an NP-hard problem.If you allow yourself to implement an approximation of Û within some threshold, you may find Ũ which might have a shorter circuit length than even an optimal decomposition of Û .This latter definition is what is generally known as circuit compression, although the term is sometimes used to refer to more optimal perfect decompositions [53].
With an error-corrected quantum computer, one could use arbitrary circuit depths, but NISQ hardware benefits from using short circuits to minimize errors from occurring.But designing hardware efficient ansatze for VQAs is an unsolved problem [54].This means that for most purposes we use a heuristic ansatz Ĉ parameterized by some vector θ brings the initial computational state |0⟩ to a desired state |ψ⟩ via Ĉ(θ) |0⟩ = |ψ⟩.The defining property of unitary matrices, namely ensures that C(θ) † reverses the action of C(θ).If we add a unitary Û , then it may be possible to find some parameters θ ′ such that Û Ĉ(θ) |0⟩ ≈ Ĉ(θ ′ ) |0⟩, thus compressing the action of the unitary back into the same quantum circuit, at least approximately.which is what is exploited by [10][11][12].
Another approach is to approximate an initial state with such an ansatz and then perform a short-time evolution via Trotterization of the time-evolved operator as is done in [25,26,37].The adjoint of the ansatz is appended to the circuit and its parameters varied such that the machine state is 'uncomputed' to its initial state.If one must assume that the chosen circuit ansatz is expressive enough to capture the entire state's time evolution, then such an approach is guaranteed to work.One then simply needs this new set of parameters and the original ansatz to express the new timestep without the time evolution operator, hence compressing the circuit.This idea has been implemented almost concurrently by Lin et al. [10], and by Barison et al.'s "projected variational quantum dynamics" (p-VQD) algorithm [11].Subsequent works, for example, [12], have built on the circuit compression idea.
Supplemental Materials: Mixed Quantum-Classical Dynamics for Near Term Quantum Computers

S-I. PROPAGATION AND SOURCES OF ERRORS IN THE TDVQP ALGORITHM
The TDVQP algorithm inherits all of the errors of its constituent parts.This includes the chosen circuit compression algorithm, time evolution approximation, and in the classical propagator.Nonetheless, it is important to have an intuition of the potential pitfalls of the algorithm.This section will begin with a short derivation of the effect of the magnitude of the infidelity of the wave function on an observable.This is followed by an analysis of the propagation through the velocity Verlet integrator, and finally, numerical experiments comparing the effect of various potential errors on the Shin-Metiu model.
We can think of an error as a superposition of the desired state |ψ⟩ some combination of undesired orthogonal states |ϕ⟩ are in a superposition of ψ = √ 1 − I 2 |ψ⟩ + I |ϕ⟩, where I is the infidelity.When we measure a Hermitian observable O we will get The actual measured observable is completely system dependent, so the effect of its magnitude on the rest of the algorithm cannot be estimated at this stage.Nonetheless, we can get an idea of the effect of this on the integrator by assuming that this directly translates to a worst-case error in the force observable, so that for an error of some magnitude we can replace F e in eq. 9 with F e (R i , |ψ i ⟩) = F i + F iϵ and propagate the first timestep.
This shows us that we are linear in the force error and quadratic in the timestep.This now enters the generation of the new position-dependent Hamiltonian such that we have H el (R 1 + Fϵ M ∆t 2 ) which is used on the subsequent step.The first evolution begins at a known position and the second is already affected by the previous error as Sadly even in this simple 1-dimensional model, it is difficult to analytically determine the effect on the evolution of the subsequent wavefunction, so we assume this effect is negligible within one timestep.We can now compute the effect on the velocity update, which yields The update can always be separated into the contributions of the ideal integration and the integration of the error.This shows us that the velocity estimation error is linear in the error in the force and linear in time.From the second timestep onwards, the error will accumulate leading to behaviour like This expression can then be expressed as the ideal contributions and the contributions from the force error as And if the error in the force is constant the expression simplifies to which is quadratic in the timestep, linear in the magnitude of the error, and quadratic with respect to the number of time steps.The error is also likely to be proportional to I 2 , which will initially be small but may increase unexpectedly as the desired position and subsequent error-prone time evolution will increasingly diverge from the ideal time evolution.
The above section illustrates that there is an effect due to the inherent interplay between the observables and the classical propagation which increases through time.But The fidelity of the wavefunction is always affected by the set optimizer threshold for fidelity.Assuming this is always met within the maximum allowed number of iterations per timestep, with a threshold of T and assuming no other errors, we expect to see the fidelity with the number of iterations i to fall as The overall effect of the threshold error is already illustrated in Figure 7, but the effect of the observable deserves numerical simulations.We have attempted two different types of errors.The first is a constant additive Multiplicative force factor Effect of a multiplicative force measurement error on ideal propagation.Here the force is a constant addition of a magnitude (ϵ) dependent error such that Fi = Fi + Fi • ϵ.Different colours represent different error magnitudes.
offset ( F = F + F ϵ ) that simply adds a force of the stated magnitude at each timestep and is shown in Figure S1.
The second is a multiplicative factor F = F + F • ϵ, which is force-dependent and its effect is shown in Figure S2.
The effect of both types of errors is quantified by comparing the fidelity of the ideal evolution to the one in which this error is injected at each timestep, but otherwise, the quantum evolution and integrator are the same as in the ideal case.We see that the effect in both is a straight line in the log-log plot of 1 − fidelity over time, which hints at there being a power law as in eq.S13.The gradient of the lines is identical in all cases within one error group (additive or multiplicative), but it is larger than 2. This may be due to the change in the Hamiltonian as the positions diverge which is not taken into account in eq.S13.
Another potential error that should be disentangled from the rest is the effect of the trajectory on the TD-VQP.To do this, the nucleus was set to follow the path it would have followed on the exact simulation irrespective FIG.S3.Effect of using a parameterized trajectory on the p-VQD algorithm, showing the infinite shot limit in blue and 10 5 shots in orange, with highlights showing the standard deviation.Solid lines are unparameterized while dashed lines are.In the infinite shot limit, there is little difference between the two, but in the finite shot case, fidelity does not fall as quickly nearing the transition point, likely due to a difference in the speed of approach of the nucleus.
of the measured observable.In effect, an evolution under an external time-dependent Hamiltonian.Figure S3 shows the effect of the parameterization, which is insignificant at the infinite shot limit, but quite noticeable in the finite shot case.The difference is particularly pronounced near the transition point (which is beginning to be approached around 200 a.u.), where the speed of the nucleus has a large effect on the transition probability, as expected by the Landau-Zener formula.The lower measured force over the simulation as seen in Figure 6 means that the speed of approach differs enough for the transition to cause a divergence in fidelity of the two.

A. Resource cost of TDVQP
To determine the overall number of circuit evaluations required by the algorithm, the analysis is split into two parts.The first one-time cost is in finding the initial state circuit parameters.This is the same as in VQE and is O(N H ), where N H is the number of Pauli strings required to express the Hamiltonian.The propagation then consists of finding the maximal overlap between the timeevolved state and its approximation, which only requires a single circuit to evaluate for the maximum number of iterations N iter , and is thus O(N iter ).The underlying optimizer will require O((N n param ), where N param is the number of circuit parameters, and n is generally small.The observable measurement requires the most circuit evaluations and is O(N obs ), where N obs is the number of Pauli strings required to express the observable in question.The whole propagation is linear in the number of timesteps desired N steps .Thus, the circuit evaluation cost is O(N H + N steps (N iter N n param + N obs )).In our example, omitting the first step, we have N obs = 7 which grows either linearly with grouped observables or exponentially without, then N param = 70 and N iter = 100, which are set by choice.This gives an overall cost of around 7 • 10 3 circuits per timestep.

B. Finite shot effects on the population
Figure S4 shows the effect of finite sampling error on the state populations.When shot noise is taken into account, the system tends to move towards the equal superposition state.Interestingly, the population transfer between states 0 and 1 is enhanced, likely due to the faster decrease of the 0 state to other states.

S-II. ON TROTTERIZATION AND ANSATZ LAYERS
Although one can have a near infinite amount of variability in the heuristic form of the ansatz, which we have chosen to be the one shown in Figure 11.But even for a single ansatz it is important to find what the optimal depth is for a given problem, and in this algorithm, we also want to see the importance of the number of steps in the Trotter approximation. Figure S5 shows that in the modified Shin-Metiu model we use the timestep of 0.05 a.u. the order of the Trotterization has a small effect compared to the depth of the ansatz.This is not unexpected, as the timestep is very small.For the simulations, we have used a depth of 5 and a single Trotter step for the best compromise between depth and precision.Another important point is whether it is beneficial to use longer timesteps within the limits of the chosen Trotterization depth or to use conservatively short steps.To see the effect of this, we can look at the fidelity of the same simulation as the single case in the main text but with 500 steps of ∆t = 0.5 [a.u.] and 5000 steps of ∆t = 0.05 [a.u.].Both have a total time of 250 [a.u.], but as should be evident from Figure S6 the two approaches show very different behaviours.When looking at the quality over 'simulated time', choosing larger timesteps is obvious.But if you look at fidelity over the number of timesteps, the shorter timesteps do have an advantage.The reason for this is because the optimizer has a fidelity threshold of ϕ compared to its previous step; as a first approximation, one can assume this fidelity is reached exactly at each timestep T , then after T steps, we would have a fidelity of ϕ T compared to our ideal situation.
Realistically, the size of the timestep also has a large effect on the iterations the optimizer must take to converge.1, 000 steps of ∆t = 0.5 take longer than 10, 000 steps of ∆t = 0.05 given the same treshold.There is also such a thing as a 'golden' initial state, which has the property that many parameters in the ansatz are initialized so that they stay constant or vary smoothly throughout the evolution.Such initializations optimize faster and retain higher fidelities than initial parameters that have much more chaotic 'spiky' evolutions.Although this is hard to quantify, it is something that could be used to filter out badly behaving initializations early on in the evolution.

S-III. PAULI STRING REPRESENTATION OF A TRIDIAGONAL HERMITIAN MATRIX
The real-space Hamiltonian for the Shin Metiu model is tridiagonal.There is a recursive solution to expressing tridiagonal Hermitian matrices in the Pauli basis.The off-diagonal matrices can be expressed as: Where X and Y are the Pauli matrices, and S π is the permutation function that returns a unique combination π of the Pauli string.That is to say that if we have the string XY Y , we would get the sum XY Y +Y XY +Y Y X.This grows exponentially, but with a qubit-wise recursive largest first commutator, it grows as O(2 n /2 ) and if we look at the fully commuting largest first approach, then the number of terms grows as O(n).The diagonal matrix is then the 2 n term weighted linear combination of all possible n length Pauli strings comprising of I and Z matrices.
Using either qubit-wise commutation relations or Pauli-word wise commutation relations, we can drastically reduce the number of Pauli strings required to implement and measure tridiagonal Hamiltonians.We have analyzed systems of up to 10 qubits via Qiskit and grouped the Pauli string decomposition through the above methods.These groupings can be used to reduce the number of measurements that have to be made, but the results are only presented here for completeness with no comment on how these measurements will be done in practice.Figure S7 shows the number of grouped terms using different existing techniques.These are not necessarily realizable on current machines.The number of grouped terms in a Hermitian tridiagonal matrix using different grouping strategies.If no strategy is used (solid blue) then the number of terms grows exponentially, which is also the case when using qubit-wise commuting groups (dashed orange).If one uses word-wise grouping then the number of terms grows linearly (dotted green).

S-IV. ADDITIONAL EXAMPLES A. Multiple transitions
The simulation parameters for these examples are synthetic, with a ∆t = 0.05 and an initial velocity of v = 0.2, which allows us to see if the algorithm can deal with more complex dynamics within 700 timesteps.All other parameters are kept as in the main text.Figure S8 shows the dynamics of the populations with 4 population crossings that are well described between states 0,1 and 2.
What is particularly interesting is that around the end of the simulation (t > 17) we see that there is a fall in occupation of the second state and an increase in the first state.This is nicely matched by what we can see in the energy in Figure S9.We see that now that the state is highly excited, population loss to lower energy states causes a drop in the energy rather than the increase we have generally seen.The plot shows the energy behaviour of the TDVQP algorithm (orange solid line) as it progresses over 700 timesteps of ∆t = 0.05 starting in the ground state, but with a fast-moving proton (with an initial velocity of v = 0.2).The ideal simulation is energy-conserving.The rise in TDVQP energy is due to leakage to higher energy levels, which then transitions to leakage in lower energy levels as the higher energy level is populated.
This fall in energy is not accompanied by an increase in fidelity, and as Figure S10 shows, the fidelity keeps decreasing at a steady rate, although likely that at very long times it would begin to oscillate at around 0.5.

B. Arbitrary state evolution
To prepare an arbitrary state we use the fact that we expect to be able to find some unitary U that can operate on a state |ψ⟩ such that U |ψ⟩ = |0⟩.If we make this circuit have the form of our desired ansatz and only allow ourselves to vary the parameters θ.We decide on some threshold or number of iterations and optimize the expression min θ (1 − | ⟨ψ| U (θ) |0⟩ | 2 ) to some threshold.This is done by running the circuit in Figure S11, where we simply initialize the quantum computer with the excited state in some way.For our simulation, we simply set the starting state to be our desired state, and we show an example for initializing in the first excited state in Figure S12 and a superposition of the two lowest BOPES in Figure S13.Figures S11 and S12 both show that the states of interest are initially well represented by the algorithm.As the evolution continues the evolution remains qualitatively similar, but degrades, especially when populations approach the 'noise floor' of the algorithm around 10 −3 where the higher energy levels are populated.Figure S14 shows a very sharp decrease in fidelity in the first timesteps and quite a large standard deviation compared to ground state results shown in Figure S10.This may either be due to not being able to initialize the ansatz as well in non-ground state settings with our approach here.Forces are still well followed qualitatively as in Figure S15.

FIG. 3 .
FIG. 3. Potential energy surfaces of the Shin-Metiu Model for the coefficients R f = 5.0, R l = 4.0 and Rr = 3.2 showing the avoided crossing around position R = −2 with L = 19 .

FIG. 10 .
FIG.10.Compression infidelity and number of iterations per step for the single TDVQP setup at the infinite shot limit, showing the mean number of iterations in red with scale to the left.The mean initial infidelity (dashed, blue) and final infidelity (solid, blue) with the scale on the right, where infidelity (I) is I = 1 − Fidelity.The number of iterations taken by the optimizer quickly reaches the limit of 100, while the final and initial infidelities gradually fall with the number of timesteps into the simulation.

FIG. 11 .
FIG. 11.One layer of the ansatz used for the VQE and TDVQP for three qubits.If multiple layers are used then the above circuit is repeated.If more qubits are used then the vertical motif is continued.In both cases, more parameters can be added as needed and ⃗ θ refers to the list of all parameters.The ansatz features x rotations and the parameterized ZZ rotation.

1 • 10 − 8 FIG. S1 .
FIG. S1.Effect of an additive force measurement error on ideal propagation.The force is a constant addition of Fϵ at each step such that Fi = Fi + Fϵ.Different colours represent different error factors.

FIG. S4 .
FIG. S4.Populations with finite sampling effect for the MD initialization simulation, showing different colours for the various shot counts, 0 being the infinite shot limit.Dashed lines represent the excited state while solid lines represent the ground state. 0 FIG. S7.The number of grouped terms in a Hermitian tridiagonal matrix using different grouping strategies.If no strategy is used (solid blue) then the number of terms grows exponentially, which is also the case when using qubit-wise commuting groups (dashed orange).If one uses word-wise grouping then the number of terms grows linearly (dotted green).
FIG. S8.State population evolution for multiple transitions.The simulation of the Shin-Metiu model is initialized in the ground state and evolved for 700 steps of ∆t = 0.05 with a fast-moving proton (v = 0.2).This induces 4 transitions.The lines show the mean BOPES state populations (solid for ideal and dashed for TDVQP) with faint lines representing higher energy levels of the TDVQP.
FIG. S9.Change in energy for the multiple transition simulation.The plot shows the energy behaviour of the TDVQP algorithm (orange solid line) as it progresses over 700 timesteps of ∆t = 0.05 starting in the ground state, but with a fast-moving proton (with an initial velocity of v = 0.2).The ideal simulation is energy-conserving.The rise in TDVQP energy is due to leakage to higher energy levels, which then transitions to leakage in lower energy levels as the higher energy level is populated.

2 Fidelity
FIG. S10.Fidelity for the multiple transition simulation.The fidelity (solid blue line) of the TDVQP wavefunctions compared to the exact evolution of the initial state for 700 timesteps of ∆t = 0.05.It can be seen that it decreases over time more quickly than the ground state evolution.
FIG. S11.Computing the ansatz parameters for arbitrary state preparation.The quantum simulator is initialized in the desired state, and the ansatz U parameters θ are varied until the final state is close to |0⟩.
FIG. S12.State population evolution of the first excited state.The simulation of the Shin-Metiu model initialized in the first excited state and evolved for 700 steps of ∆t = 0.05.Initial conditions are the same as in Section IV A and the lines show the mean BOPES state populations.The faint lines represent higher energy levels of the TDVQP.
FIG. S13.State population evolution of the superposition state.The simulation of the Shin-Metiu model initialized in an equal superposition of the 0 and 1 states showing the mean BOPES state populations and evolved for 700 steps of ∆t = 0.05.Initial conditions are the same as in Section IV A and the lines show the mean BOPES state populations.The faint lines represent higher energy levels of the TDVQP.
FIG. S15.Mean force measurements for an equal superposition of the lowest BOPES states (solid, blue) and first excited state (orange, dashed) compared to their respective ideal simulations (black).
Mean fidelity with respect to the depth of the ansatz and the number of trotter steps in the evolution for 10 samples of 200 timesteps.The lowest depth is 3 (blue), while the next step is of depth 4 (orange), and the deepest is 5 (green).The different Trotter approximations are various dashed or solid lines, but since the time evolution is short, this has little effect compared to the depth of the ansatz.