Adaptive projected variational quantum dynamics

We propose an adaptive quantum algorithm to prepare accurate variational time evolved wave functions. The method is based on the projected Variational Quantum Dynamics (pVQD) algorithm, that performs a global optimization with linear scaling in the number of variational parameters. Instead of fixing a variational ansatz at the beginning of the simulation, the circuit is grown systematically during the time evolution. Moreover, the adaptive step does not require auxiliary qubits and the gate search can be performed in parallel on different quantum devices. We apply the new algorithm, named Adaptive pVQD, to the simulation of driven spin models and fermionic systems, where it shows an advantage when compared to both Trotterized circuits and non-adaptive variational methods. Finally, we use the shallower circuits prepared using the Adaptive pVQD algorithm to obtain more accurate measurements of physical properties of quantum systems on hardware.


I. INTRODUCTION
Simulation of static and dynamic properties of quantum systems is a notoriously hard task for classical computers.While analytical solutions are available only for specific cases, the amount of time and computing resources required in general by exact numerical methods is exponential in the system size, making the calculations quickly unfeasible.While several approximated many-body numerical techniques have been proposed [1][2][3][4], the accurate description of important physical and chemical phenomena is a very active research problem [5][6][7][8].
In recent years, quantum computers have seen significant developments [9][10][11], opening potential opportunities for scientific discoveries.Hardware capabilities continue to advance steadily, and we can already create and manipulate complex many-body quantum systems [12][13][14][15][16][17].However, large-scale fault-tolerant quantum computers remain far in the future, and contemporary devices show limitations in connectivity, size, and coherence times.
Accounting for these constraints, Variational Quantum Algorithms (VQAs) have emerged as the leading strategy to take advantage of near-term quantum devices [18][19][20][21].In this class of algorithms, the solution of a given problem (e.g.finding the ground state of a physical system) is encoded in a quantum circuit that depends on some parameters optimized with the aid of a classical device.VQAs have not only been proposed for quantum simulations but also for a variety of different applications, such as machine learning [22,23], combinatorial optimization [24,25], quantum error correction [26,27] and compila- * david.linteau@epfl.chtion [28][29][30].Variational schemes have also been introduced in quantum dynamics [31][32][33][34][35][36][37][38][39], as a more efficient alternative to Trotterization [40][41][42][43][44].The accuracy of a variational quantum simulation is then tied to the ability of a parameterized circuit to describe time-evolved wave functions.Even if the initial wave function is welldescribed by the chosen parameterized circuit, the complexity of the time-evolved wave functions varies with time and the chosen circuit may fail to describe them.The choice of the parameterized circuit is therefore crucial and it remains an open problem in variational simulations of quantum dynamics.
Adaptive schemes have been proposed in the context of variational ground state search [45][46][47][48] especially to avoid committing to a particular parameterized circuit.The key idea is to construct the parameterized circuit during optimization.By systematically appending specific quantum gates to the parameterized circuit, adaptive methods have been shown to surpass standard approaches in the number of operations required and in the accuracy of the final results.Moreover, adaptive methods provide flexible circuits suited for dynamics simulations [33,49].However, including an adaptive step for dynamics usually requires measurements of additional quantities, that might be difficult to perform, or auxiliary qubits.
In this work, we introduce an adaptive variational algorithm for real-time evolution based on the projected Variational Quantum Dynamics (pVQD) algorithm [36], denoted Adaptive pVQD.The method inherits all the properties of the original pVQD algorithm and integrates the adaptive modification of the parameterized circuit without requiring auxiliary qubits.The structure of this paper is as follows: in Section II we present the algorithm and describe how the adaptive routine is per-formed; in Section III we apply the method to study a time-dependent and a fermionic system, benchmarking the method against Trotter evolution and the original pVQD algorithm; Section IV concludes the paper with some considerations and outlooks on the proposed method.

II. METHOD
Consider a physical system governed by a Hamiltonian H.For clarity of exposition, we focus on timeindependent Hamiltonians.However, this is not a requirement of the algorithm, as we explicitly show in Section III.To simulate the dynamics of quantum systems on a quantum computer, we have to prepare the time-evolved wave function |Ψ(t) = U (t)|ψ 0 , where |ψ 0 = U 0 |0 ⊗N is the initial state, N indicates the number of qubits representing the physical system and U (t) is the unitary time evolution operator.The Adaptive pVQD algorithm aims to approximate the state |Ψ(t) using parameterized states of the form where each real parameter θ i ∈ θ is associated to a Hermitian generator A i ∈ A. The parameterized state is therefore specified by the set of parameters and operators {θ, A}, and it can be implemented as a quantum circuit.
To simulate a physical model until a final time t f , we divide the evolution into small time intervals ∆t.We further assume that the parameterized state |ψ(θ) is a good approximation of the time-evolved wave function at time t.The wave function at time t + ∆t can thus be represented by U TS (∆t)|ψ(θ) , where U TS (∆t) is a Trotter-Suzuki decomposition of the time evolution operator U (∆t) [40,41].In this manuscript we use a first order decomposition, but higher orders can be considered.The choice of the optimal ∆t is problem dependent and will be discussed in Section III.We then approximate the evolution step t → t + ∆t using a new set of parameters θ → θ + dθ that maximizes the overlap between U TS (∆t)|ψ(θ) and |ψ(θ + dθ) .This can be achieved by minimizing, with respect to dθ, the infidelity where the fidelity can be measured on a quantum device [36].
At each time step, the initial parameters and operators {θ, A} are those obtained at the previous time step.Assuming that the set of operators A is sufficient to describe the state at time t + ∆t, we find the parameter shift dθ * that minimizes I(dθ, ∆t).Details about the minimization routine can be found in Appendix A. If the minimization routine is not successful, new gates built using generators from the operator pool are added to the parameterized circuit following the adaptive procedure described in Section II A. This adaptive procedure is repeated up until the convergence criteria are met.
The algorithm starts with the initial state |ψ 0 represented by an empty set of operators.As needed, new gates are added through the time evolution until the chosen final time t f .The complete procedure is illustrated in Fig. 1.We note that the original pVQD scheme [36] can be recovered by fixing the set of operators A through the entire simulation.

A. Adaptive step
When the parameterized circuit |ψ(θ) is not expressive enough to accurately describe the time step evolution by only shifting the variational parameters, we add new gates to it.This is referred to as the adaptive step of the algorithm.Given an operator pool, we determine the best gate to grow the quantum circuit.As first proposed in [45], we look for the operator whose gate maximizes the derivative of the cost function with respect to its parameter.This is achieved by iterating over all the operators in the pool, a step that can be performed in parallel even on different quantum devices.
For ground state methods, the cost function is the energy of the system, while the gradient is obtained by measuring the expectation value of the commutator between the trial operator and the Hamiltonian [45,50].We must ensure that is possible to apply a similar procedure when dynamics is considered.In the adaptive scheme proposed in [33], this step requires an additional measurement of the variance of the Hamiltonian with respect to the nonadaptive case.In our method, the gradient of the fidelity with respect to the shift dθ a of parameter θ a associated with a trial operator A a has the form where we define the projector P 0 = |ψ 0 ψ 0 | and the state |φ(θ, ∆t) = U † (θ)U TS (∆t) |ψ(θ) (see Appendix B for the full derivation).To ensure continuity of time evolution, we initially set θ a = 0. We note that measuring the derivative of the fidelity corresponds to measuring the Hermitian operator [P 0 , iA a ] with respect to the pVQD circuit U † (θ)U TS (∆t) |ψ(θ) modified by the addition of the gate e idθaAa .However, we evaluate the derivative using the parameter shift rule [51], as for the minimization routine (for more details, see Appendix A).This operator search is still parallelizable on multiple devices and does not require auxiliary qubits.
The adaptive step has been lately extended and optimized [46,48,50], with new protocols that greatly reduce the computational resources required with respect to the first proposal.In particular, we adopt the scheme presented in [48], which increases the depth of the parameterized circuit |ψ(θ) by 1 at every adaptive step.While the infidelity defined in Eq. ( 2) remains above a fixed threshold ε, additional adaptive steps are performed.For a detailed description, see Appendix C.

B. Operator pool
The choice of the operator pool is a key ingredient in the success and efficiency of adaptive variational algorithms.Having a complete pool of operators is exponentially complex in the size of the physical system, therefore, one has to make some restrictions in its selection.Many different strategies have been proposed, such as the creation of a minimally complete pool [46,52], the inclusion of symmetries directly in the operator pool [53], or the extension of a complete pool acting on a subsystem of the studied model [47].
In the study of the dynamics, we can refer to the Trotterization of the time evolution operator to select the pool.In particular, we consider local (L) and non-local (NL) operator pools, respectively, given by where X i , Y i and Z i are the Pauli gates acting on site i.Given that A L ⊂ A NL , we expect that A NL will generate more flexible parameterized states.However, not only the choice of A NL leads to a measurement overhead, but the non-local nature of this pool may add long-range controlled-NOT (CNOT) gates to the circuit, according to the device connectivity.In Section III, we report the comparison of the two pools in the study of a fermionic system.

III. RESULTS
We apply the Adaptive pVQD method to the study of the 1D Heisenberg XYZ model with an external driving field and the 2D Fermi-Hubbard model.Both have non-trivial dynamics and open the pVQD method to the study of time-dependent and fermionic systems.In both cases, open boundary conditions were imposed.

A. Driven Heisenberg model
Given an open chain of L spins, the driven Heisenberg XYZ Hamiltonian can be written as: where J x , J y and J z are coupling parameters and D(t) is the time-dependent driving term.Many different driving terms can be applied to the system.Among those we choose where ω is the driving frequency.
First, we investigate the performance of the Adaptive pVQD algorithm with a local pool on a perfect simulator and compare to Trotterized circuits and the original implementation of pVQD.We consider J x = 1, J y = 0.8, J z = 0.6, an antiferromagnetic initial state |ψ 0 = |0101 and a final evolution time t f = 2.In the clas-sic version of the pVQD algorithm, we have to choose an ansatz for the time evolved wave function.We consider a circuit equivalent to a Trotter step where all the rotations are defined by variational parameters.The Trotter step circuit implementation for this model is shown in Appendix E. Both the Trotter and the pVQD full circuits are then obtained repeating this structure n TS times.In particular, we fix n TS = 10 for the Trotter circuit and n TS = 3 for the pVQD ansatz.
After running the algorithms, we compare the different circuits obtained and use them to measure expectation values of single-and two-spin observables.The results are shown in Fig. 2. The Trotter circuit lags behind variational methods both in terms of accuracy and resource required.The pVQD method instead achieves accurate results until t = 1.0,where the associated circuit becomes shallower than the one of Adaptive pVQD.This phenomenon suggests that in that time step the fixed representation power is the main source of error in the variational calculations.
In order to show the flexibility of the Adaptive pVQD, we implement a naive modification of the pVQD algorithm, that we indicate as pVQD with block extensions.In this case, a new step of the Trotterized variational ansatz is added to the circuit once the optimization procedure does not reach the desired accuracy.While this approach does improve the performance of the pVQD algorithm, we remark that it is not general, as it depends on the ansatz structure we have chosen.Furthermore, we can see from the bottom panel of Fig. 2 that the Adaptive pVQD method always produces shallower circuits, with resources tailored to the needs of the specific time step.
Then, we extend the study to systems with different sizes.
To this end, we define the integrated exact infidelity with respect to the exact wave function |Ψ(t) computed on a classical device.We again fix a final evolution time t f = 2 and evaluate ∆ ex I (t f ) for each method for systems of L ∈ [3,11] spins.In particular, we consider a Trotter circuit with a fixed depth of n TS = 10 and one with fixed Trotter step size dt = J x t/n TS = 0.05, the same we use in the Trotter step of the pVQD algorithm.The results are shown in Fig. 3, together with the circuit depth at the end of the time evolution.
We note that the depth of the Adaptive pVQD circuits increases with the system size and converges to the Trotter circuit with fixed depth, while having a lower integrated exact infidelity.We highlight that Fig. 3  Final depth steps are required to evolve the system to t f = 2, resulting in circuits almost one order of magnitude deeper than any other.We performed multiple pVQD simulations with different variational ansätze equivalent to n TS = 1, 2, 3, 8 Trotter steps.We note that the integrated exact infidelities of pVQD with n TS = 1, 2, 3 have a steep transition when the number of gates becomes smaller than the adaptive circuit.This phenomenon suggests that the ansatz limitation is the main source of error in the variational calculations, while the adaptive circuit is able to increase effectively its representation power.
On the other hand, the standard pVQD calculation with n TS = 8 never undergoes this transition.While the integrated exact infidelity is always lower than the adaptive approach, we have to note that the entire time evolution is performed with a deeper circuit.Finally, we note a plateau in the depth of the circuit required by the adaptive algorithm when L > 8.This is similar to what observed in [33], where the system size at which the number of gates required saturates depends on the evolution time.
The adaptive method is able to produce circuits that are orders of magnitude shallower than Trotterization while keeping the accuracy comparable to it.Those circuit can be used to improve the measurement of observables at long times on current quantum devices, which are otherwise limited by the depth of the Trotterization.For this reason, we first run the Adaptive pVQD algorithm on the simulator and use the resulting sets of variational parameters to prepare quantum circuit on the hardware for a system of L = 4 spins.In Fig. 4, we compare observables measured both on those variational wave functions and on Trotterized circuits with a fixed Trotter step size of dt = 0.2.
In this experiment, the final Trotter circuit has 180 CNOTs.This circuit is beyond what is currently accessible on quantum devices, settling the expectation value of the correlator close to 0 for J x t > 0.8.On the other hand, the Adaptive pVQD parameterized circuit |ψ(θ) has 28 CNOTs at the end of the evolution.This improvement in the number of gates is crucial for the application of error mitigation techniques, especially at longer times.
In particular, zero noise extrapolation (ZNE [31,54]) was applied both on the noisy simulations and hardware experiments.We choose a quadratic fit on values obtained with noise scaling factors [1,2,3].Moreover, when running our algorithm on hardware, we dynamically decouple the idle qubits from the active ones using the standard procedure available in Qiskit [55].We expect that more advanced noise mitigation techniques, such as the one presented in [56], will improve the results on the Trotter circuit.However, this is also true for the variational circuit prepared by the Adaptive pVQD.

B. Fermi-Hubbard model
The Hamiltonian of the Fermi-Hubbard model on a L x × L y rectangular lattice is given by where c † iσ (c iσ ) is the creation (annihilation) fermionic operator of spin σ ∈ {↑, ↓} at site i, n iσ = c † iσ c iσ counts the number of fermions with spin σ at site i and ij denotes nearest neighbor sites on the lattice.The first term in the Hamiltonian accounts for the hopping between nearest neighbor lattice sites, while the second term describes the on-site interactions.
There are several ways to encode fermionic Hamiltonians into qubit operators [57][58][59][60][61][62][63].In this work, we consider the Jordan-Wigner mapping [57] to encode each fermionic mode into a qubit.Since every lattice site can host two modes (↑, ↓), N = 2L x L y qubits are required to simulate the Fermi-Hubbard model on a L x × L y grid.Before performing a fermionic encoding, we eliminate the spin index via c i↑ → c i and c i↓ → c i+N/2 (and analogously for the number operator n iσ ).We then map each fermionic operator into a spin operator: where σ ± = (X ± iY )/2.The local occupation number can then be identified with the local spin number according to n i ∈ {0, 1} → Z i ∈ {↑, ↓}.More details on the fermionic indexing convention and implementing a Trotter step can be found in Appendix E.
Given that the mapping requires an ordering of the fermionic modes, operators that are local in space might generate very long Pauli strings.For example, considering the snake-like pattern, vertical hopping terms generate strings of Pauli Z with sizes up to 2L x − 2. This represents a bottleneck in studying fermionic systems with dimensionality higher than 1 on current quantum devices.By restricting the operator pool, we investigate the possibility of describing time-evolved wave functions of the 2D Hubbard model using only local gates.We perform noiseless simulations of a 2 × 2 square lattice, comparing local and non-local operator pools.In particular, we measure the expectation values of a local density operator and a density correlator and count the number of CNOTs in the circuits.We use a fixed-depth Trotter simulation and a pVQD with block extension as a benchmark.The results are shown in Fig. 5.
We do not restrict ourselves to specific quantum hardware to keep the comparison as general as possible.Instead, we count the number of CNOTs in a circuit by transpiling it into an abstract device with all-to-all connectivity that is able to perform arbitrary single qubit rotations and CNOTs.The local and non-local pool variants show different behavior over time in the count of CNOTs.We note that the non-local variant always requires fewer CNOTs than its local counterpart.However, some CNOTs are long-range, and their implementation on an actual device can be challenging on hardware with fixed topology and limited connectivity.In contrast, the circuit structure produced by the local pool variant is already suited for current hardware implementation.More details about the Adaptive pVQD output circuits can be found in Appendix D.Moreover, the plot highlights another limitation of the naive pVQD with block extensions approach.Indeed, it always prepare more expensive circuits than the Adaptive pVQD with non local pool and in the end it has similar CNOT requirement to the local variant, while being restricted to use long range gates as required by the Trotter step.

IV. CONCLUSIONS
We presented an adaptive version of pVQD, called Adaptive pVQD, to simulate the real-time evolution of quantum systems.This algorithm importantly circumvents the need to choose a fixed ansatz from the beginning of the time evolution.The parameterized quantum circuits are grown adaptively to be both problems and hardwaretailored.This is obtained with a measurement overhead required to determine the best gate among those included in the operator pool.However, the gate search can be operated in parallel and, in our scheme, does not involve circuits with auxiliary qubits.This makes the Adaptive pVQD algorithm more hardware-efficient than standard methods, as exemplified in this work with the driven Heisenberg model on the IBM quantum hardware.Finally, we have simulated the dynamics of the 2D Hubbard model with only local gates, using the adaptive procedure to mitigate one of the bottlenecks that current quantum devices face in studying fermionic systems.Given the ease of introduction to the standard pVQD algorithm and its benefits, we believe that the adaptive procedure described here can be of great use in the simulation of dynamics both for current and future quantum devices.

DATA AVAILABILITY
The code used to run the simulations is open source and can be found at [64].It was written in Python using Qiskit [55].Exact classical simulations were performed using Qutip [65].
U (θ)|ψ 0 , we want to add the gate e −idθaAa to it, defining the new state |ψ(θ + dθ) = U (θ) e −idθaAa |ψ 0 .To obtain the gradient of the fidelity with respect to this added parameter dθ a , it is convenient to first rewrite the fidelity given in Eq. which precisely corresponds to Eq. ( 4).

Figure 1 .
Figure1.Flowchart of the time evolution of the Adaptive pVQD algorithm.Starting with a parameter-free circuit, we discretize the time evolution into multiple time steps.At each time step we optimize the parameters to approximate the real time evolution of the quantum system.If the optimization does not converge to the required accuracy, or the ansatz does not contain any parameter, then rotations {R A * i } based on the generators {A * i } are appended to the circuit according to the adaptive step procedure described in Section II A. The algorithm stops once the final time t f is reached.

Figure 2 .
Figure 2. Dynamics of the driven Heisenberg XYZ model studied with the Adaptive pVQD algorithm with local pool (L), compared to standard Trotter evolution, pVQD and pVQD with block extensions.The plot shows the results for an open chain of L = 8 spins with Jx = 1, Jy = 0.8 and Jz = 0.6.The top and middle panel show the measurements of a single spin observable and a correlator, respectively.The bottom panel shows the number of CNOTs in the circuit describing the time-evolved wave function.The simulation started in the antiferromagnetic state |ψ0 = |01010101 , and the infidelity threshold was set to ε = 10 −4 for all the variational methods.

Figure 3 .
Figure 3. Adaptive pVQD algorithm with local pool compared to standard Trotter evolution and pVQD for the driven Heisenberg XYZ model.We employ the same setup indicated in Fig. 2 for multiple systems of size L ∈ [3, 11].The top panel shows the integrated exact infidelity of pVQD and Trotterization over an entire time evolution with final time t f = 2 as a function of the system size.The bottom panel shows the circuit depth at the end of the time

1 Figure 4 .
Figure 4. Observables measured with the IBM Manila device for the driven Heisenberg XYZ model on an open chain with 4 sites, Jx = 1, Jy = 0.8, Jz = 0.6 and an antiferromagnetic initial state |ψ0 = |0101 .The Trotter simulation is performed with a fixed Trotter step size of dt = 0.2.The Adaptive pVQD circuits |ψ(θ) were obtained with a noiseless simulation that used a local operator pool.The shaded areas correspond to 50 noisy simulations using the noise model of IBM Manila.Each data point and error bar correspond to the mean and the standard deviation, respectively, of 50 experiments performed on hardware.Zero noise extrapolation was applied to both noisy simulations and hardware experiments.Idle qubits were also dynamically decoupled from the active ones.

Figure 5 .
Figure 5. Adaptive pVQD schemes for the Fermi-Hubbard model on a 2×2 open square lattice (8 qubits) with U/J = 0.8.Local (L) and non-local (NL) operator pools are used to perform noiseless simulations and the results are compared to a Trotter evolution with nTS = 5 Trotter steps and pVQD with block extensions.The system starts in the half-filled antiferromagnetic state |ψ0 = |n 0↑ n 1↑ n 2↑ • • • = |10100101 .We fixed the infidelity threshold to ε = 10 −4 .The top and middle panels show the expectation values of an on-site number density operator and a number density correlator over time.The bottom panel shows the number of CNOTs in the circuit describing the time-evolved wave function.