Quantum algorithms for grid-based variational time evolution

The simulation of quantum dynamics calls for quantum algorithms working in first quantized grid encodings. Here, we propose a variational quantum algorithm for performing quantum dynamics in first quantization. In addition to the usual reduction in circuit depth conferred by variational approaches, this algorithm also enjoys several advantages compared to previously proposed ones. For instance, variational approaches suffer from the need for a large number of measurements. However, the grid encoding of first quantized Hamiltonians only requires measuring in position and momentum bases, irrespective of the system size. Their combination with variational approaches is therefore particularly attractive. Moreover, heuristic variational forms can be employed to overcome the limitation of the hard decomposition of Trotterized first quantized Hamiltonians into quantum gates. We apply this quantum algorithm to the dynamics of several systems in one and two dimensions. Our simulations exhibit the previously observed numerical instabilities of variational time propagation approaches. We show how they can be significantly attenuated through subspace diagonalization at a cost of an additional $\mathcal{O}(MN^2)$ 2-qubit gates where $M$ is the number of dimensions and $N^M$ is the total number of grid points.


Introduction
Simulating quantum dynamics is of foremost importance to understand a multitude of chemical processes.Despite an impressive progress in the development of computational methods [1,2,3,4], accurate calculations remain restricted to molecules with less than a few tens of collective degrees of freedom [5,6,7,4].An alternative path towards the efficient modeling of quantum dynamics is to switch towards a new computational paradigm.In particular, quantum computers have the potential to simulate quantum systems in polynomial time and memory [8,9,10].Following Feynman's thesis, Wiesner [11] and Zalka [12] first designed a framework to simulate molecular quantum dynamics on a digital quantum computer with a grid encoding of real space.This framework was then made concrete and applied to the simulation of several small nuclear quantum systems in rectangular or harmonic potentials [13,14,15,16,17].More recently, it was also extended to the simulation of the non-adiabatic dynamics of nuclear wavepackets [18,19].All these approaches share the same circuit representation of the time evolution operator, obtained from Trotter approximation of the latter [20].However, in general, this leads to very deep quantum circuits, which greatly exceed the capacities of present quantum hardware due to their limited coherence times.This is particularly true when the Hamiltonian is given in first quantization as, in this case, an efficient encod-ing of even a single Trotter step into a quantum circuit is cumbersome [18,21,22].As a consequence, the type of implementable potentials is restricted.For instance, all the works cited above only present dynamics under potentials that can be defined with a polynomial function of the position.It is worth emphasising here that active research is undergoing to implement Coulombic potentials [23,24,25,22].
To address the issue of circuit depth, variational time evolution (VTE) quantum algorithms were proposed.For a comprehensive overview see for instance Refs.[26,27].Relying on an iterative exchange of information between a classical and a quantum computer, these algorithms allow to work with shallower quantum circuits of constant depth in time.In particular, Li and Benjamin [28] first showed how to use a variational principle to simulate the real time dynamics of quantum systems and applied it to a quantum Ising model.
The goal of the present work is to extend this approach to the simulation of wavepacket quantum dynamics in a grid-based encoding for general Hamiltonians given in first quantization.The impact of the resulting algorithm lies beyond the ability of performing simulations on noisy near-term quantum hardware as it enables grid-based quantum dynamics in the first place, which can be a highly non-trivial task with a Trotter approximation.We highlight the emergence of strong numerical instabilities in this context and present a local diagonalization scheme for making the variational time evolution algorithm stable and efficient.We expect this novel approach to be easily extendable and beneficial to any problem instance similarly discretized, as is the case for example in nonlinear problems [29], quantum field theory [30], electron-phonon systems [17], or quantum risk analysis [21].

Wavepacket dynamics in position space
In this work, we study the wavepacket dynamics governed by the time-dependent Schrödinger equation (TDSE) where x is a spacial dimension mapped here to a one or two-dimensional grid.The wavefunction is normalized, |ψ(x, t)| 2 dx = 1, so that the |ψ(x, t)| 2 becomes the probability probability for finding the system at position x of the grid at time t.The time variable is considered as a continuous parameter.The Hamiltonian in Eq. ( 1) is given by where the first term is the kinetic energy of the system and the second term describes a static external potential.In this framework, the timeevolution operator takes the form For an initial wavepacket ψ(x, t 0 ) (mostly not an eigenstate of the Hamiltonian) the time evolution is given by ψ(x, t) = e −iH(t−t 0 )/ℏ ψ(x, t 0 ) which can be evaluated numerically by means of standard integration schemes (such as the Trotter decomposition).Note that the kinetic operator (first term in Eq. ( 2)) becomes diagonal in the momentum representation, i.e., after applying the Fourier transformation: x → p = −iℏ∂/∂x; it is therefore a common practice to define integrators of the time-evolution evolution operator, which leverage the advantages of both position and momentum representations of the system Hamiltonian.
Alternatively, one can also express the time evolution through a variational principle.We will delve deeper into this concept in the upcoming section.

VTE in grid encodings
The variational approach to quantum dynamics aims to approximate the solution of the TDSE on a low-dimensional submanifold of the full Hilbert space.The trial wavefunction defined on this manifold |ψ(θ(t))⟩ is parameterized by a set of n p time-dependent parameters, θ ≡ θ(t) = {θ 1 (t), ..., θ np (t)}.For a given Hamiltonian H, the McLachlan variational principle [26,31,32,33,34] leads to the following equations of motion with Further details on these equations are available in the supplementary information.The calculation of the F kj and V k matrix and vector elements are classically hampered by the high dimension of |ψ(θ)⟩ which ensures accurate simulations.Instead, they can be efficiently measured on a trial wavefunction encoded in the state of a qubit register.In this case, the trial wavefunction is defined as where |ϕ⟩ is a reference state and U(θ) is a unitary operator (e.g., the quantum circuit) depending on real parameters θ ≡ θ(t).The measurement of analytic gradients on quantum computers was discussed in Ref. [35].If the parameters θ are chosen as the rotation angles of single qubit gates, there always exists a simple circuit In the supplementary information we recall the canonical approach to computing F kj and V k when the Hamiltonian is written as a weighted sum of Pauli tensor strings, H = p c p P p .
Here, instead, we will focus on Hamiltonians expressed in first quantization, H = p 2 /2m + V(r), where p is the momentum, m the mass, and V(r) the potential given as a function of the position r.H describes an M -dimensional systems.The time-evolution is directly performed in momentum and position space, discretized on a grid.The N points, per dimension, of the grid are encoded in the basis states of N q = log 2 (N ) qubits.The total number of qubits is then M N q .In this case, no explicit transformation to a basis representation is required, saving the numerical effort associated with such computations, which would scale as the square of the basis set size [36].Note that the encoding (real space discretization on a grid) employed here is different from the discrete variable representation (DVR) of Ref. [37].We argue in favor of the present grid encoding as the DVR representation leads to an exponentially growing number of Hamiltonian terms.In our case, the expectation value of the Hamiltonian is simply obtained from two sets of measurements, one in the momentum basis and one in the position basis.In fact, this constitutes an important advantage compared to previous implementations of the same VTE quantum algorithm, e. g., in second quantization or in the DVR representation.Indeed, the high number of measurements required to perform the VTE becomes quickly impractical when increasing the dimensionality of the problem [38,39,40].This advantage appears in the calculation of the right-hand side of Eq.S3, which requires measuring The first term is straightforward and reduces to measuring the expectation values ⟨ψ(θ)| p 2 |ψ(θ)⟩ and ⟨ψ(θ)| V(r) |ψ(θ)⟩.In practice, this is obtained by repeatedly preparing and measuring the state |ψ(θ)⟩ in the computational (position) basis.For each outcome of the binary representation of r, V(r) is classically computed and ⟨ψ(θ)| V(r) |ψ(θ)⟩ is obtained by averaging over all realizations of V(r).The same is done for the kinetic term by measuring in the momentum basis, i. e., applying a quantum Fourier transform (QFT) right before the measurement (see Ref. [18] and supplementary information).
The second term in V k can be rewritten as The scalar ℜ ⟨ϕ| W † (θ)V(r)U(θ) |ϕ⟩ is then obtained from the calculation of E[(−1) s V(r)] with r and s the measurement outcomes of the quantum register and ancilla qubit, respectively, obtained from the following quantum circuit .
The proof is given in the supplementary information.The same is done for the kinetic term by measuring in the momentum basis, i. e., applying a QFT right before the measurement.

The wavefunction representation and quantum advantage
The quantum advantage of the approach presented herein relies on the possibility to accurately approximate the wavefunction with a quantum computer, using a number of variational parameters, n p , that is much smaller than the size of the full Hilbert space.The parameters can then be propagated efficiently according to Eq. (S3), as long as n p scales favorably (i.e., polynomially) with the system size.
When working in first quantization, the translation of a Trotter operator T into a quantum circuit can be very costly in terms of the gate count.Indeed, an efficient decomposition of for an arbitrary potential energy function V(r) into a set of quantum gates is not guaranteed to exist.In fact, the quantum circuit for encoding a general function requires either exponentially many gates or ancilla qubits to perform quantum arithmetic.In the latter case, the targeted function is generally approximated and the number of ancilla qubits scales exponentially with the inverse of the desired accuracy [41,21,22,18].
On one hand, this speaks in favor of the use of variational approaches for time-propagation over Trotter-based algorithms since the former does not require the implementation of an accurate (and costly) gate decomposition of the operator T .On the other hand, finding appropriate, physically motivated, variational Ansätze is still a challenge, and heuristic, hardware efficient approaches are therefore preferred.Hence, as a first demonstration in low dimensions, we will limit ourselves to heuristic variational forms, which can be implemented on quantum hardware.These can be systematically improved by repeating the same quantum circuit units with independent parameters, i. e., by increasing the circuit depth d.
The different heuristic circuits employed in this work are detailed in the supplementary information.

Simulations in position and momentum spaces
In the following, we study the performance of the quantum algorithm for concrete applications in first quantization.We seek to perform this study without introducing any quantum hardware noise bias.Hence, the simulations are obtained in a perfect classical emulator of a quantum computer.The expectation values are computed from matrix-vector multiplications.For this reason, we approximate the wavefunction derivatives with forward finite differences of stepsize 10 −8 (which should not be confused with the adaptive time-step used for the time propagation).We show in the supplementary information that this step-size leads to converged results.The equations of motion (Eq.(S3)) are solved with a least-squares approach as implemented in NumPy [42] with a cutoff ratio for small singular values, or reconditioning number, arbitrarily set to 10 −6 .To solve the ordinary differen-tial equations, we employ a state-of-the-art adaptive Runge-Kutta solver of order 5(4) available in SciPy-routines [43].We first test the VTE approach in one-dimensional systems defined on a grid as described in Ref. [18] and in the supplementary information.The length of the box is L = 14.We study three different problems: a freely moving particle, a particle in a harmonic potential (harmonic oscillator), and a particle colliding with an Eckart barrier.The Hamiltonians for these three systems are given by We work in atomic units, taking m = 1, c 1 = 1, c 2 = 13, and c 3 = 3/2.In all cases, the state is initialized to a Gaussian wavepacket, with A the normalization factor, B the width of the wavepacket, and p 0 and x 0 the initial momentum and position, respectively.The evolution is carried out for a total time t tot = 1.5.Snapshots of the exact time evolution (from exact exponentiation of the Hamiltonian matrix) for each of the three systems are shown in Fig. 1 (left panels) for times t = 0, 0.45, 0.91, and 1.5.
The space is discretized with N q = 6 qubits (corresponding to 64 grid points).In this case, the full Hilbert space can be described by n full p = 2(2 6 − 1) = 126 real parameters.It is clear that a quantum advantage can only be achieved with a much smaller number of variational parameters, satisfying n p ≪ n full p .We perform VTE simulations for each of the three aforementioned systems with the variational form vf1 (see supplementary information) for different depths (or circuit repetitions, see supplementary information).
The initial parameters θ 0 are found by maximizing the fidelity where |0⟩ is the vacuum state.Note that this procedure, which is sub-optimal in general, is implemented here for convenience.An extensive discussion on state initialization goes beyond the scope of this work and is left for future investigations.The initial conditions are (x 0 , p 0 ) = (−3.5, 5) for the free particle and the Eckart barrier, and (x 0 , p 0 ) = (−3.5, 2) for the harmonic oscillator.In all cases, the width of the initial wavepacket is set to B = 1/ √ 2. The fidelities as a function of the simulation time, are shown in Fig. 1 (right panels).These results demonstrate that, in general, the number of variational parameters required to maintain a fidelity above 95% throughout the entire simulation time always approaches n full p , particularly so when tackling hard problems such as the scattering off an Eckart barrier, reducing the possibility for quantum advantage.The evolution of all parameters is also given in the supplementary information and shows sharp changes in their trajectories for each of the three systems.To our understanding, supported by the detailed study given in the supplementary information, the above observations can be rationalized as follows.The chosen heuristic variational forms have enough flexibility to accurately and efficiently, i. e., with few variational parameters, approximate the targeted wavefunctions at all times of their dynamics.Therefore, the loss of accuracy observed throughout the different simulations is not due to the variational forms, but is an inherent effect of the VTE.When the number of variational parameters is insufficient, the dynamics strongly depends on the numerical setup of the simulation (grid mesh, reconditioning number, initial parameters, etc.).These observations are in agreement with Ref. [37] where strong numerical instabilities were also put forward.Surprisingly, the correct dynamics are always recovered when increasing the number of variational parameters.In this case, the algorithm is stable.In addition, we observe that the more complex the dynamics gets, the more parameters are needed to obtained a robust and stable dynamics.Here, the term complex refers to systems with a large support in position and frequency spaces in their canonical representation, which lack of any particular symmetry reduction property.For completion, we also show in the supplementary information the robustness of the variational approach to the introduction of hardware noise.
We repeated our simulations in the momentum representation by defining a variational Ansatz of the form The results are given in Fig. 1.As expected, the evolution of the free particle can now be performed accurately with very few parameters, since the introduced momentum basis diagonalizes the Hamiltonian.Interestingly however, the evolution is also smoother in the case of the harmonic oscillator and the Eckart barrier.Note that we do not expect the harmonic oscillator to be symmetric in position and momentum representations, since the position and momentum grids are different and the corresponding potentials come with different prefactors.In the supplementary information, we show that this momentum space representation improves the dynamics compared to the position space in all cases considered in this work.Results obtained by combining both the position and the momentum representation in a common Ansatz are also given in the supplementary information but do not show substantial improvements.
In order to search for efficient heuristic Ansätze, we systematically change the characteristics of the variational circuit, namely the single-qubit rotations, the type of entangling gates, and the connectivity of the entangling block (see supplementary information for details).We apply all these different variational forms to the simulation of the wavepacket interacting with the Eckart barrier.The results presented in the supplementary information show that, in general, the accuracy is improved by working in the momentum representation and increasing the circuit depth.However, no clear trend for the design of more efficient variational forms could be identified.

Simulations in local diagonal space
The introduction of unitary transformations that achieve a (at least partial) diagonalization of the Hamiltonian seems essential to make the quantum VTE efficient in the grid representation.The reasoning behind this approach is that a significant part of the dynamics can be associated with the eigenstate evolution in uncoupled subspaces, such that only the remaining, often smooth evolution needs to be treated in a variational setting.
We elaborate on this concept in the supplementary information.
In the supplementary information, we show some preliminary one-dimensional tests demonstrating that an approximate diagonalization of the system Hamiltonian can already give a sufficiently accurate and stable time evolution.Of more interest is the application of this approach to multi-dimensional systems.Following our previous observations, we expect that, for a multidimensional configuration space, an improvement can already be achieved by diagonalizing each dimension independently, i. e., by fixing all but one coordinate of the Hamiltonian.This leads to a quantum circuit of the form We call the space obtained from this basis transformation the local diagonal (LD) space.In the next paragraph, we discuss the scaling of this approach.
As stated earlier, we employ N grid points, corresponding to N q = log 2 (N ) qubits, for each of the M dimensions.Hence, the total number of qubits is M N q .The total size of the problem is N M , prohibitively large for direct classical simulations in high dimensions.Efficient algorithms should scale to low polynomial order in N and M .Our approach first requires the diagonalization of M matrices of size (N, N ).This step generally scales as O(N 3 ).The resulting unitary is then translated into a quantum circuit.The classical scaling of this operation is O(N q N 3 ) [44], leading to O(N 2 ) CNOT gates [45].Finally, the VTE is performed with a variational form comprising n p parameters.The partial diagonalization step ensures the scaling of n p to be of low polynomial order, i. e., O (M N q ) P where P is small enough such that (M N q ) P ≪ N M (see the supplementary information for a more comprehensive description of the LD basis and its effect on n p ).At each time step of the VTE, the equations of motion are reconstructed with O(n 2 p ) measurements and solved classically at cost O n 3 p .Recent approaches such as the one of Ref. [39] focus on reducing the scaling in n p for obtaining the M matrix and could be extended to the present case.Finally, we obtain an efficient hybrid quantumclassical time evolution algorithm with polynomial scaling in N and M .As a proof of concept, we test this approach on the evolution of a two-dimensional system on a Mexican Hat potential energy surface (PES).The Hamiltonian reads with r 2 = x 2 + y 2 , c 4 = 0.1, c 5 = 1, and again m = 1.The space is discretized with 8 qubits (4 per dimension), corresponding to a total of 256 grid points.Note that the Hilbert space can be fully represented by 510 real parameters.The wavepacket is evolved in the region of space x ∈ [−5, 5] and y ∈ [−5, 5].It is initialized at (x 0 , y 0 ) = (−3.0,0) with no initial momentum and width B x = B y = 1/ √ 2. The momentum increases in subsequent time steps as the wavepacket slides down the brim of the PES.A graphical representation of the system is given in Fig. 2a together with snapshots of the exact evolution for times t = 0, 0.91, 1.8, and 3.0.
We simulate the quantum dynamics with VTE up to t tot = 3.0, employing the same settings as before, and using both the position and the LD Ansatz.The one-dimensional Hamiltonians H x = H(x, 0) and H y = H(0, y) are diagonalized to get the quantum circuits D x and D y .For both the evolutions in position and LD space, the parameterized part of the quantum circuit U(θ) corresponds to the variational form vf1 with depth 20 and depth 25 (336 and 416 parameters, respectively).The results are displayed in Fig. 2b  and c.Significant improvements are obtained in LD space.Note that these results were obtained from a choice of one-dimensional Hamiltonians that are very simple to diagonalize.They could be improved even further by exploiting symmetries of the system or by diagonalizing low-dimensional mean-field Hamiltonians.

Conclusions
We introduced a quantum VTE algorithm to perform nuclear wavepacket dynamics on a grid in first quantization.
This variational algorithm manifests a crucial advantage compared to Trotter-like quantum approaches to this problem class, namely the fact that it does not require the direct implementation of the time-evolution operator (exponentiating the Hamiltonian) in the qubit register.Furthermore, we stressed the advantage of our method in relation to the need of sampling expectation values in only two bases representations (position and momentum), irrespective of the system size.We studied the performance of the quantum algorithm in classical emulations for several one-and two-dimensional systems.In general, we observed strong numerical instabilities when performing the dynamics with an efficient number of variational parameters.However, we could demonstrate that the accuracy of the quantum algorithm can be improved by expressing the variational quantum circuit in a problem specific basis.This basis is obtained by diagonalizing each dimension of the system independently, without introducing significant computational overhead.Moreover, adaptive approaches for variational time evolution [31] might help to further reduce ansatz cost in situations where our diagonalization procedure is not applicable.However, selecting a suitable operator pool for constructing the adaptive ansatz could be challenging, since these approaches typically rely on the Trotterization of the time evolution operator.As mentioned previously, this is not generally efficient in first quantization.
We discussed the overall cost of the proposed approach, which shows an effective polynomial scaling in both the number of grid points per dimension, N , and the number of system dimensions, M (for a total of N M grid points).It is worth noting that, treating higher dimensional grids, the qubit connectivity of the respective quantum hardware will have significant impact on implementation details of our approach and the ordering of qubits will have to be taken into consideration.Going further, this approach can be easily extended to the simultaneous treatment of electronic and nuclear degrees of freedom, opening up new opportunities for the simulation of non-adiabatic dynamics beyond the Born-Oppenheimer approximation.When dealing with identical particles, the system wavefunction needs to be adapted accordingly: antisymmetrized for fermions and symmetrized for bosons upon particle exchange.The generation of quantum circuits fulfilling invariance under (anti-)symmetrization is discussed in detail in the literature [24].

A The grid encoding
In this work, real space is discretized on a grid which points are mapped to the basis states, |j⟩, of the quantum register.We illustrate the grid encoding of this work in a one dimensional example.The length of the grid is L.There are N points and N q = log 2 (N ) qubits.The initial point in position space is x 0 = −L/2.Each subsequent point is defined as x j = x 0 + L N −1 j. j is an integer which binary representation is given by |j⟩.In the reciprocal (momentum) space the points are then defined as p j = p 0 + 2π L j with p 0 = − N π L .This shift allows to account for nagtive values of the momentum.This choice implies the use of a centered Quantum Fourier Transform (cQFT) operator to implement the switch from the position to the momentum space.This cQFT can simply be implemented by adding a X gate on the last qubit right before and after the QFT (for cQFT, i. e. cQFT ≡ (I ⊗ I ⊗ ... ⊗ X)QFT(I ⊗ I ⊗ ... ⊗ X)) and QFT −1 (for cQFT −1 ) such that they respectively undergo a cyclic permutation.The equation defining this permutation is given below for QFT as an illustration:

B The time-dependent variational principle in first and second quantization
The variational approach to quantum dynamics aims to approximate the solution of the TDSE on a low-dimensional submanifold of the full Hilbert space.The trial wavefunction defined on this manifold |ψ(θ(t))⟩ is parameterized by a set of n p time-dependent parameters, θ(t) = {θ 1 (t), ..., θ np (t)}.A time-dependent variational principle (TDVP) defines the optimal evolution of the parameters within the submanifold of the full Hilbert space.There exist different formulations of the TDVP.For Kähler manifolds, i. e., when the tangent space is a complex subspace [34], all the different formulations lead to the same equations of motion for the variational parameters.However, this is not the case for unitary parameterizations of the type where |ϕ⟩ is a reference state and U(θ(t)) is a unitary operator (e.g., the quantum circuit) depending on real parameters θ(t) [34].In this case, the equations of motion differ and hold distinct properties such as the conservation of the norm or the energy.Recent works [26,31,32,33,34] promoted the use of the equations of motion derived from the McLachlan variational principle due to of their higher numerical stability [34].For a given Hamiltonian H, these equations, when accounting for a global phase mismatch (see Ref. [26] for a thorough derivation), read with For a regular parametrization, Eq. (S4) represents a quantum metric (the Fubini-Study metric) in parameter space.In conjuction with Eq. (S4), Eq. (S5) shows that Eq. (S3) results from the orthogonal projection of the exact time-derivative of |ψ⟩ on the variational manifold.Note that the variational manifold is entirely general, and any parametrized form of the wave function can be used In this work, however, we will focus on Hamiltonians expressed in first quantization, where p is the momentum, m the mass, and V(r) the potential given as a function of the position r.
H describes an M -dimensional systems.The time-evolution is directly performed in momentum and position space, discretized on a grid.The N points, per dimension, of the grid are encoded in the basis states of N q = log 2 (N ) qubits.The total number of qubits is then M N q .In this case, no explicit transformation to a basis representation is required, saving the numerical effort associated with the calculation of the integrals, which would scale as the square of the basis set size [36].
The expectation value of the Hamiltonian is now simply obtained from two sets of measurements, one in the momentum basis and one in the position basis.This gives a clear advantage in the calculation of the right-hand side of Eq.S3, which requires measuring ⟨ψ(θ)| H |ψ(θ)⟩ and ℑ(⟨∂ The first one is straightforward and reduces to measuring the expectation values ⟨ψ(θ)| p 2 |ψ(θ)⟩ and ⟨ψ(θ)| V(r) |ψ(θ)⟩.In practice, this is obtained by repeatedly preparing and measuring the state |ψ(θ)⟩ in the computational (position) basis.For each outcome of the binary representation of r, V(r) is classically computed and ⟨ψ(θ)| V(r) |ψ(θ)⟩ is obtained by averaging over all realizations of V(r).The same is done for the kinetic term by measuring in the momentum basis, i. e., applying a QFT right before the measurement (see Ref. [18]).
Let's now consider the calculation of the components V k , where In what follows, we simplify the notations of the quantum circuits as U = U(θ) and W = W(θ).
The scalar ℜ ⟨ϕ| W † HU |ϕ⟩ is obtained in a similar way as in the second quantization case, namely with the circuit of Eq. (S7) but by measuring both the ancilla qubit (in the computational basis) and the variational state register.This is given by the following quantum circuit, We prove this by first considering the potential part of the Hamiltonian.In this case, ⟨ϕ| W † VU |ϕ⟩ is calculated from E[(−1) s V(j)] where s is the measurement outcome of the ancilla qubit, while j is that of the register.Indeed, the probability of measuring states s and j from the above circuit is With this, we finally obtain The same can be shown for the kinetic part of the Hamiltonian, provided the register is measured in the momentum basis, i. e., by introducing a QFT before the measurement.

C Heuristic variational forms
The different variational forms employed in this work are summarized in Tab. 1.They all comprise alternating layers of single qubit rotations and entangling blocks.They differ in the choice of single qubit rotation gates, entangling gates, and coupling map, i. e., in the geometrical way the qubits are coupled to each other.In the linear coupling map, each qubit is coupled to its two nearest neighbors only.The circular map is similar but adds a coupling between the first and last qubits.Finally, in the full coupling map, each qubit is coupled to all other qubits.For the sake of clarity, Fig. S1 gives concrete examples of the resulting circuits for each variational form for four qubits.

Name Single qubit gates Entangling gates Coupling map
# params.for Nq qubits and depth Table 1: Definitions of the different variational forms employed throughout this work.The blue colour highlights a parameterized gate.

D Comparing variational forms
One of the main advantages of the variational approach described in this work for first quantized Hamiltonians is the possibility to perform the dynamics without explicitly including the Hamiltonian in the quantum circuit.Instead, heuristic variational forms are exploited.There exist many ways the heuristic variational forms can be defined.Here, we test several ones for the dynamics of the nuclear wavepacket colliding with an Eckart barrier (6 qubits) and compare the resulting fidelities over the simulation time.All variational forms are defined in section C of this supplementary information and we explicitly give the depths and corresponding number of parameters in Tab. 2. From the results shown in Fig. S2 one can see that, in general, an increase in accuracy is obtained by enlarging the depth and by going to momentum space.However, no clear trends can be identified as to how to construct the variational form.Nonetheless, variational forms 1 and 5 (with linear connectivity of controlled-Z entangling gates), seem to generally perform the best.

E Time evolution of the variational parameters
In this section, we show the dynamics of the parameters obtained with the different VTEs to supplement the fidelity results shown in the main text.In Fig. S3, we report the evolution in position space, in Fig. S4 those in momentum space, and, finally in Fig. S5, the evolutions in diagonal space.The first two cases (position and momentum) are obtained with variational form vf1, while in the last case (diagonal), we use variational form vf2. For all systems, in position space, we observe sharp changes in the dynamics.In general, the parameter values also diverge with time.The same observations can be made in momentum space for the harmonic oscillator and the Eckart barrier.
Whenever the circuit diagonalizes the Hamiltonian as in Fig. S4(a) and Fig. S5(a) and (b), the parameter evolutions become smooth or even trivial.This is unsurprising as, in this case, there is no transfer of amplitudes but only evolution of the phases.Note that in Fig. S5(b), when cut > 0, the partial diagonalization of the Hamiltonian suffices to ensure a smoother evolution of the parameters.

F Study on the accuracy of the time evolution in one dimension
To better understand the origin of the errors occurring during the VTE, we perform additional numerical simulations.Unless explicitly stated otherwise, the results of this section are obtained with variational forms in position space.We start with the Eckart barrier case.The first step is to identify whether the variational form is flexible enough to represent the state at each time step.For this, we optimize the parameters obtained from VTE in 6 qubits with depth 5 and variational form vf1 to maximize F(t) at each time step (see definition in main text).The results displayed in the top panel of Fig. S6 show that we can indeed get a better fidelity for times t > 0.8.This suggests that the variational form is good enough to represent, with good accuracy, the exact wavefunction through the entire dynamics.
There exist several degenerate sets of variational parameters which represent the same wavefunction.In other words, the initial state can be prepared with a given accuracy from different sets of parameters using the same Ansatz .In the bottom panel of Fig. S6, we show how the accuracy of the dynamics can be affected by the choice these initial parameters.For the three different trials displayed in Fig. S6, we find those initial parameters by maximizing F(t = 0).Each one of these three optimization processes starts with random guesses and converges to above 99% fidelity to different parameters solutions.Fig. S6 clearly shows a difference in the fidelities of the dynamics for the three trials.However, we observe that choosing the best set of initial parameters is non-trivial and cannot be made by simply looking at the initial fidelity or early time observables such as the local-in-time error [48].
As a second step, we aim to study the effect of the different numerical parameters on the accuracy.We choose to work with the simplest system: a free particle.The space is discretized with 5 qubits.The initial conditions are (x 0 , p 0 ) = (0, 5).The width of the initial wavepacket is set to B = 1/ √ 2. The parametrized circuit corresponds to variational form vf2 with depth 3.
We first run the VTE for different values of the finite difference step size, ϵ.The results are displayed in Fig. S7, showing identical fidelities in the relevant part of the evolution when the accuracy is above 95%.
We then fix ϵ back to its original value of 10 −8 and change the reconditioning number, rc, the ratio for cutting off small singular values in the least-squares algorithm.The results displayed in Fig. S8(a) show that the fidelity increases with smaller rc.This implies numerical errors coming from instabilities in the inversion of the matrix F when solving Eq.S3 (main text).
We then study the performance of the ordinary differential equation solver.We employ the same Runge-Kutta 5(4) solver, but now fix the maximum time step to 10 −4 , and compare to an explicit Runge-Kutta solver of order 8 (DOP853) [43].As shown in Fig. S9, the results of these two simulations are identical in the first part of the dynamics (before the accuracy drops below 90%) and differ afterwards.
Finally, we vary the initial width B of the wavepacket.Interestingly, we observe that the results are much improved when increasing the initial width (see Fig. S10).The latter influences the overall spread of the wavepacket during the evolution.The difference between initial and final width is 0.841, 0.148 and 0.028 when B = 1/ √ 2, B = 2/ √ 2 and B = 3/ √ 2, respectively.These results suggest that the spread of the wavepacket is difficult to capture with the grid-based VTE.
To validate those observations, we run the VTE of a wavepacket oscillating in a harmonic potential without spread.The space is discretized in 6 qubits as in the simulations of the main text.We also keep the same initial conditions (x 0 , p 0 ) = (−3.5, 2), and the same variational form, vf1 depth 5 in position space.This time, however, the initial width is taken to be B = B gs = 0.6, where B gs is the width of Hamiltonian's ground state.Fig. S11 confirms our previous observations: the results in the case of a non-changing width are much improved compared to the results of the main text in which the wavepacket's width changes over time.This leads to the conclusion that the VTE performs reasonably  well (with few parameters) for simple dynamics.
We then study the effect of increasing the precision in the grid mesh on the dynamics of the harmonic oscillator.The initial conditions are the ones of the main text: B = 1/ √ 2 and (x 0 , p 0 ) = (−3.5, 2).The space is discretized with 6, 7 and 8 qubits corresponding to 64, 128 and 256 grid points, respectively.The variational form is always taken to be vf1.The results are shown in Fig. S12.We employ both the position and momentum representations of the wavefunction (indicated with x and p in Fig. S12, respectively).The different depths and corresponding number of parameters are also indicated in the legend.The fidelities shown in Fig. S12 are computed from the variationally time evolved wavefunctions with respect to the exactly evolved ones discretized on the same grid.In other words, both the reference and the variational wavefunctions are expressed in the same number of qubits.In all cases, we see that the correct dynamics are recovered in the limit of the number of parameters, n p , approaching the size of the Hilbert space.In the momentum space, the number of parameters needed to maintain an accuracy of F > 95% over the whole time range remains low for N q = 6 and N q = 7.It also shows a good scaling behavior when going from 6 to 7 qubits.Indeed, if 60 parameters are necessary in the 6-qubit case, this number only raises to 84 when we double the number of grid points (7 qubits).However, we observe a drop in accuracy when going to 8 qubits.In this case, 416 parameters (for a full Hilbert space represented with 510 real parameters) were not enough to get accurate dynamics both in position and momentum spaces.These results show the influence of the grid mesh on the accuracy of the VTE.The observed non-monotonic behavior pinpoints the strong correlation existing between the different numerical factors, such as N q and rc.It is important to note that those numerical instabilities are always corrected for when the number of parameters is high enough.
We repeat the simulation of the same harmonic oscillator system for different values of rc, the cutoff  ratio of small singular value for the inversion of the matrix F. Those dynamics are obtained with 6, 7, and 8 qubits.The results are shown in Fig. S13.We employ variational form vf1 in position (x) and momentum (p) space, and for different depths (d) as indicated in the legend of Fig. S13.As opposed to the results of Fig. S8, in this case, decreasing the value of rc does not improve the results but even worsens them slightly.On the other hand, we observe improved results when increasing rc with thinner a grid mesh (Fig. S13(c)).This shows again that the numerical effects are correlated and system dependent.
In conclusion, we highlight the following point: • Heuristic variational forms have the flexibility to accurately and efficiently, i. e., with few variational parameters, describe the targeted wavefunctions at all times of their dynamics.The loss of accuracy observed throughout the different simulations is an inherent effect of the method.
• When the number of variational parameters is insufficient, the dynamics strongly depend on the numerical setup of the simulation (grid mesh size, reconditioning number, initial parameters, etc).
• The correct dynamics are always recovered by increasing the number of variational parameters.4) solver used throughout this work but with a maximal time step fixed at 10 −4 .On the other hand, the dashed line was obtained with an explicit Runge-Kutta method of order 8 as implemented in SciPy [43].These results are obtained for the case of the free particle with 5 qubits and variational form vf2, depth 3.
In this case, the algorithm is stable.
• The more complex the dynamics is, the larger is the number of needed parameters.The term complex relates to the size of the energy eigenspace involved in the dynamics (number of eigenvectors of H spanning the subspace in which the dynamics is defined).In practice, few symmetries and a large range of positions/frequencies involved are characteristics of complex dynamics when working with the canonical first quantized representation of the Hamiltonian.
• In all cases presented here, the expression of the wavefunction in momentum space improves the results.

G Variational forms mixing position and momentum spaces
In the main text, we discuss the improvement of the results when defining the Ansatz in the momentum space, i. e., by adding a QFT at the end of the circuit.Here, we show the results obtained by mixing momentum and position space in the variational form.More explicitly, the Ansatz is now composed of several parts; first, a part of given depth in position space, followed by an inverse QFT, then another piece of variational circuit with its own depth, and finally a QFT closing the circuit.We refer to the part enclosed by the QFTs as the momentum part.We vary the depths of the position and momentum parts.The results, shown in Fig. S14, do not highlight any improvement in the performance.

H Time evolution in presence of noise
The robustness of variational approaches for quantum dynamics to hardware noise has been previously observed.[38] To evaluate the performance of our algorithm in presence of noise we modify the equations of motion as shown in Ref. [26] and work with the system's density matrix ρ(t).The unitary evolution of mixed states defined by ρ(t) under H is governed by the von Neumann equation The equations of motion for the variational parameters then become where and  In practice, ρ(θ) is extracted at each time step from a noisy quantum simulation performed with the AerSimulator (in density_matrix mode) of Qiskit [49] with a given NoiseModel.The equations of motion are then solved classically with exact matrix/vector multiplications.This allows us to study the effect of particular types of hardware noise without adding the one of finite sampling on top.The fidelities, F(t), shown in this section are calculated as where ρ ex (t) = |ψ ex (t)⟩ ⟨ψ ex (t)| with |ψ ex (t)⟩ = e −iHt |ψ 0 ⟩, the exact state at time t.We select a system for which there is no numerical instabilities in the noiseless case, namely the free particle system discretized with 5 qubits with variational form vf2 depth 3 and initial width B = 3/ √ 2 (see Fig. S10).In all the following noisy simulations the state is initialized with the parameters obtained without noise, i. e., the same ones employed in the results shown in Section F, Fig. S10.
First, we apply a depolarizing channel implemented in Qiskit [49], which is defined as where N q is the number of qubits and λ is a parameter, 0 ≤ λ ≤ 4 Nq /(4 Nq − 1).We repeat the simulations for various values of λ and show the results in Fig. S15 Because the initial parameters are obtained from a noiseless optimization, the fidelity at t = 0 does not lie above 99% anymore when λ increases.On the other hand, the fidelity does not significantly drop during the time evolution, suggesting that our algorithm is robust to a certain amount of noise.Second, we replace the depolarizing channel with thermal relaxation noise as implemented in Qiskit [49].The operation time are set to 0 ns for the Rz gates, since they can be virtually applied by changing the phase of the subsequent pulses, and 36 ns and 500 ns for the Ry and CZ gates, respectively.The thermal relaxation time constant T 1 and the dephasing time constant T 2 are chosen to be equal and their value is varied from 50 µs to 500 µs (see Fig. S16).In this case, we again only observe a small drop in accuracy during the dynamics for all noise regimes.This strengthens the conjecture that variational time evolution approaches show good resilience to quantum hardware noise.

I Local diagonal space
In this section we discuss the addition of a unitary circuit for rotating the quantum state in a basis which partially diagonalizes the Hamiltonian.We first classically diagonalize the harmonic oscillator and Eckart barrier Hamiltonians.We then map the unitary matrix, made up of the eigenvectors sorted by increasing order of energy, to a quantum circuit using the isometry decomposition of Ref. [45] as implemented in Qiskit [49].The resulting circuits, D, are then appended to the variational form as |0⟩ ⊗Nq / U (θ) D .
To obtain the results of Fig. S17, we add D to the variational form vf2.Those results show a high state fidelity over the entire simulation time.In the case of the Eckart potential, we also tested unitaries D ′ , which only perform a partial diagonalization of the Hamiltonian.More specifically, prior to diagonalization, we set to zero all matrix elements with absolute values below cutoff thresholds of cut = 0.1 and 1.0.We stress here that cut should not be mixed up with rc which is instead the reconditioning number used in the least-squares algorithm.The density of non-zero elements in the resulting Hamiltonian matrix is then 0.98, 0.27, and 0.14 for cutoffs 0, 0.1, and 1.0, respectively.The results exhibit high accuracy for all cutoff values as seen from Fig. S17b.In the following we give insight on the reduction of the number of necessary variational parameters when working in the local diagonal space.The aim of the local diagonalization is to find a unitary transformation of the form D = M i=1 D i which captures a great part of the dynamics, and turns the variational problem into a smooth dynamical problem.This can be achieved by choosing D such as where the |Φ 0 n ⟩ are the eigenstates, associated to eigenvalues E 0 n , of the simpler Hamiltonian H 0 related to the original problem's Hamiltonian via a perturbation Ṽ as Note that H 0 can be adapted along the dynamics e. g., could be defined as a time-dependent mean-field Hamiltonian.
For the sake of clarity, let us first consider a variational Ansatz defined in the basis of the eigenstates |Φ 0 n ⟩ and with the variational parameters, α n , simply being the amplitudes associated with those eigenvectors: Clearly, under the above assumption the parameters are expected to vary smoothly in time with the limit αn − → 0 for Ṽ − → 0.Moreover, since H 0 is expected to describe H reasonably well, the energies E 0 n are close to the E n (eigenvalues of H).Because of the conservation of energy, the dynamics is approximately confined to an invariant subspace A (defined by the initial state).In other words, the |Φ 0 n ⟩ lying far in energy from the eigenstates spanning A are very unlikely to enter the variational wavefunction at any point during the dynamics.The subspace A is only spanned by a limited number of eigenvectors of H 0 , A = span{|Φ n ⟩} n∈I A , for a subset I A of indices n.Hence, a considerable reduction of complexity is apparent, since one can decrease the number of parameters by neglecting the α n for which n / ∈ I A .To stay on the safe side, one can use a sector A of the Hilbert space that is larger than the one strictly needed to describe the initial state, and that (in the language of perturbation theory) can accommodate the most important virtual transitions.Even in this case, the bottom line remains unaltered: a reduction of complexity is possible since H 0 captures the essence of the dynamics, hence its eigenspaces are quasi-invariant under the action of H.Although our problem is not in the form of Eq.S22, since our θ k are not directly the amplitude associated to a given |Φ 0 k ⟩ but real parameters entering the variational Ansatz , |ψ(θ)⟩ = U(θ) |ϕ⟩, the same observation applies.Our variational parameters only need to describe subspace A. As a result, the number of parameters, n p , must be on the order 2 dim(A) to fulfill the requirement that for each |Φ⟩ ∈ A there exists a set of parameter θ such that There can be some flexibility in n p depending on the target accuracy for the state representation.In any case, n p now scales with dim(A) leading to a reduction of the overall complexity.

Figure 1 :
Figure 1: VTE dynamics for the three one-dimensional systems considered in this study.left panels Snapshots of the modulus square of the exact wavefunctions at times t = 0.00, 0.45, 0.91, and 1.5 (lightest to darkest curve, respectively) for a) the free particle, b) the harmonic oscillator, c) the Eckart barrier.right panels Fidelity, F as a function of time, t, of the VTE dynamics in position (x) and momentum (p) space for a) the free particle, b) the harmonic oscillator, c) the Eckart barrier.The results are obtained with 6 qubits and variational form vf1 at depth d (with corresponding number n p of variational parameters).

Figure 2 :
Figure 2: VTE dynamics of a two-dimensional system, representing a wavepacket evolving on a "mexican hat" potential.a) 3-dimensional representation of the system at the top and, at the bottom, snapshots of the modulus square of the exact wavefunction (darker color indicates larger amplitude) at times t = 0.00, 0.91, 1.8, and 3.0.Beside are the results of the VTE dynamics obtained in b) the position space and c) the LD space.The fidelity, F as a function of time, t, is obtained with 8 qubits and variational form vf1 at depth d (with the corresponding number n p of variational parameters).

Figure S1 :
Figure S1: Quantum circuits for each variational form with depth d on 4 qubits.The blue highlights a parameterized gate.

Table 2 :
Figure S2: Fidelity of the variationally time-evolved wavefunction with respect to the exact one as a function of time for top-left all variational forms with short depth in position space, top-right all variational forms with short depth in momentum space, bottom-left all variational forms with large depth in position space and bottom-left all variational forms with large depth in momentum space.All varioational forms are defined in Section C. The specific depths and corresponding number of parameters are given in Tab. 2 together.These results are obtained for the case of the Eckart barrier with 6 qubits.

Figure S3 :
Figure S3: Time evolution of the variational parameters for each one-dimensional system obtained with variational form vf1 in position space and at given depth, d (with n p variational parameters).

Figure S4 :
Figure S4: Time evolution of the variational parameters for each one-dimensional system obtained with variational form vf1 in momentum space and at given depth, d (with n p variational parameters).

Figure S5 :
Figure S5: Time evolution of the variational parameters for two of the one-dimensional systems obtained with variational form vf2 in diagonal space and at given depth, d (with n p variational parameters), and cutoff value cut.

Figure S6 :
Figure S6: Top: Fidelity of the variationally time-evolved wavefunction with respect to the exact one as a function of time (full line) and after being optimized at each time step to maximize the fidelity (dotted line).Bottom: Fidelity of the variationally time-evolved wavefunction with respect to the exact one as a function of time for three different sets of initial parameters.These results are obtained for the case of the Eckart barrier with 6 qubits and variational form vf1, depth 5.

Figure S7 :
Figure S7: Fidelity of the variationally time-evolved wavefunction with respect to the exact one as a function of time and for different values of the finite difference step, ϵ.These results are obtained for the case of the free particle with 5 qubits and variational form vf2, depth 3.

Figure S8 :
Figure S8: Fidelity of the variationally time-evolved wavefunction with respect to the exact one as a function of time.Different ratios, rc, for cutting off small singular values are employed.These results are obtained for the case of the free particle with 5 qubits and variational form vf2, depth 3.

Figure S9 :
Figure S9: Fidelity of the variationally time-evolved wavefunction with respect to the exact one as a function of time, obtained with different ordinary differential equation solvers.The full line corresponds to the Runge-Kutta 5(4) solver used throughout this work but with a maximal time step fixed at 10 −4 .On the other hand, the dashed line was obtained with an explicit Runge-Kutta method of order 8 as implemented in SciPy[43].These results are obtained for the case of the free particle with 5 qubits and variational form vf2, depth 3.

Figure S10 :
Figure S10: Fidelity of the variationally time-evolved wavefunction with respect to the exact one as a function of time and for different values of the width, B, of the initial wavepacket.These results are obtained for the case of the free particle with 5 qubits and variational form vf2, depth 3.

Figure S11 :
FigureS11: Top Snapshots of the exact (dashed lines) and the variationally time-evolved (full lines) modulus squared wavefunctions at times t = 0.00, 0.45, 0.91, and 1.5 (lightest to darkest curve, respectively) of the harmonic oscillator with a non-spreading width B = 0.6.Bottom Fidelity of the variationally time-evolved wavefunction with respect to the exact one as a function of time.These results are obtained for the case of the harmonic oscillator with width B = 0.6, 6 qubits, and variational form vf1 depth 5 in position space.

Figure S12 :
Figure S12: Fidelity of the variationally time-evolved wavefunction with respect to the exact one as a function of time.These results are obtained for the case of the harmonic oscillator with variational form vf1. The space is discretized with (a) 6 qubits, (b) 7 qubits, and (c) 8 qubits.

Figure S13 :
Figure S13: Fidelity of the variationally time-evolved wavefunction with respect to the exact one as a function of time.These results are obtained for the case of the harmonic oscillator with variational form vf1, and for different values rc for cutting off small singular values.The space is discretized with (a) 6 qubits, (b) 7 qubits, and (c) 8 qubits.

Figure S14 :
Figure S14: Fidelity of the variationally time-evolved wavefunction with respect to the exact one as a function of time.These results are obtained for the case the Eckart barrier with 6 qubits.The quantum circuit corresponds to variational form vf1 with alternating layers in position and momentum space of depths given in the legend.

Figure S15 :
Figure S15: Fidelity of the variationally time-evolved density matrix with respect to the exact one as a function of time and in presence of depolarizing noise.These results are obtained for the case of the free particle with 5 qubits and variational form vf2.

Figure S16 :
Figure S16: Fidelity of the variationally time-evolved density matrix with respect to the exact one as a function of time and in presence of thermal relaxation noise.time constants T 1,2 are given in µs.These results are obtained for the case of the free particle with 5 qubits and variational form vf2.

Figure S17 :
Figure S17: Fidelity, F as a function of time, t, of the VTE in LD space for (a) the harmonic oscillator, (b) the Eckart barrier.The results are obtained with 6 qubits and variational form vf2 at depth d (with corresponding number n p of variational parameters).In the case of the Eckart barrier, the fidelity results are given for different cutoffs, cut, of the Hamiltonian coefficients used to obtain the diagonalization unitary.