Carole Chatel
email: carole.chatel@iphc.cnrs.fr

Ludovic Mathieu

Mourad Aïche

Maria Diakaki

Olivier Bouland

Development of a small Time-Projection-Chamber for the quasiabsolute neutron flux measurement

Accurate actinides fission cross sections around 1 MeV are of primary importance for the safety of generation IV reactors. To have accurate measurements, the neutron flux must also be accurately estimated. This is usually done with respect to the 235 U(n,f) cross section. It is however possible to measure the neutron flux with respect to the 1 H(n,n)p cross section which is a primary standard, providing an independent and precise measurement. Typically, the usual proton recoil technique uses a silicon detector for neutrons of energy between 1 and 70 MeV. However, the high electron and gamma background due to neutron production under irradiation makes the use of this or any other detector not suitable for an accurate measurement below 1 MeV. To this end, the Gaseous Proton Recoil Telescope is developed and characterized. The goal is to provide quasi-absolute neutron flux measurements with an accuracy better than 3%. It consists of a double ionization chamber with a Micromegas segmented detection plane and the gaseous pressure can be adjusted to protonsand hence neutronenergy. The sensitivity to gamma and electrons background, the intrinsic efficiency as well as the resolution of this detector have been investigated.

Introduction

Discrepancies between different measurements around 1 MeV still exist for most of actinides fission cross sections, despite a general great effort from the nuclear community. The neutron fluxes being high in pressurized water reactors (generation II and III) and in fast neutron reactors (generation IV) from few hundreds keV to few MeV, it is of prior importance to reduce the uncertainties in this neutron energy range.

In most of the measurements, fission cross sections are measured with respect to secondary standards ( 235 U(n,f), 238 U(n,f), 237 Np(n,f)). Their accuracy ranges from 1.3% to 10% between 200 keV and 20 MeV [START_REF] Carlson | Evaluation of the neutron Data Standards[END_REF]. Moreover, in the energy range we are interested in, the usual reference is mostly 235 U(n,f). Thus, the constant use of the same standard creates a strong correlation between different measurements.

For this reason, it has become important to realize an independent and precise measurement using a different neutron flux measurement technique. The use of a primary standard instead of a secondary one, is thus suitable. The 1 H(n,n)p elastic scattering cross section is known with an accuracy ranging between 0.2% to 0.5% from 0 to 20 MeV [START_REF] Carlson | Evaluation of the neutron Data Standards[END_REF] and reproductible with ab initio calculations. This cross section is also structureless, creating no dependence between the measured structures and ones from the reference reaction used.

For a neutron energy range between 1 and 70 MeV, silicon junctions are perfectly adapted as proton recoil detectors [START_REF] Marini | [END_REF][3][4][5]. However, the desired measurements start around 200 keV (down to the fertile fission thresholds) and Si junction are here no more suitable. Indeed, they are sensitive to a high electron and gamma background under irradiation. As our setup is as close as possible to the neutron source to maximize the flux, many gamma and electrons produced by the source are detected in the recoil-proton detector [START_REF] Marini | [END_REF] and adapting the detector thickness is not a sufficient solution [3].

Even if other proton recoil detectors have been developed as plastic scintillator [6][7][8][START_REF] Kovash | Detecting sub-MeV neutrons in solid plastic scintillator with gamma-ray discrimination[END_REF][START_REF] Schmidt | [END_REF] or proportional counters [11][12][13], the inaccuracy of their efficiency makes these detectors not suitable for our goal (see [START_REF] Marini | [END_REF]14,15] for details).

As no detector seems to already exist to meet all our needs, the ACEN team of the LP2I decided to develop the Gaseous Proton Recoil Telescope (GPRT). This detector should not be sensitive to high electrons/gamma background and able to accurately measure neutron fluxes for neutron energies between 0.2 MeV and 2.5 MeV.

GPRT description

The developed detector is a double ionization chamber (E-E). The electric field is generated by the presence of two electrodes separated by 4 cm of gas. The amplification region is based on a Micromegas technology. It consists of the anode plane (12x4 cm²), segmented in 64 pads and an electrified mesh of the same dimensions, 125 µm above the segmented plane. A field cage is used to homogenize the electric field between the cathode and the mesh [16].

At the entrance of the detector stands a macor sample disk. Its use is to carry several H-rich foils upstream the detector. Then, a first collimator stands as the detector entrance, delimiting part of the detector geometry. A second collimator 2 cm farther delimitates the E chamber from the E one. A picture of the set-up is shown in Figure 1 and the segmented detection plane is zoomed in Figure 2.

Detection system specific characteristics

As the aim of this detector is to detect proton recoil, it is of primary importance that as less protons as possible come from another source than the H rich sample. For this reason, the structures have been chosen to be of Hfree materials: either of macor (collimators, sample disk, etc.), or metallic (cathode, field cage, etc.). For the same reason, the used gas is H-free: a mixture of N 2 and CO 2 . The exact composition will be discussed in 3.4.

The pressure of the gas can be adapted to the neutron energy so the proton stops inside the detector. However, the H-rich foil thickness must also be adapted to the neutron energy so the loss of energy of the recoil proton is minor in the foil. For this reason, the sample disk has been designed to carry up to 8 H-rich foils. A motor allows to change the sample without opening the detector chamber. Its axis, not arrowed on Figure 1, is protected by an insulating plastic cylinder. The electric field, by the polarization voltages, is then adapted to the chosen pressure: as high as possible to maximize the amplification gain but low enough to avoid electrical breakdowns. Breakdown voltages have been investigated as a function of the pressure and, to prevent micro-breakdowns, it has been chosen to work up to 20 V below the breakdown voltage.

To keep the proton recoil energy close to the neutron energy that created it, the proton scattering angle is limited thanks to the second collimator. The presence of this collimator allows a coincidence between the E chamber and the E chamber, making this detector a telescope.

Finally, once the recoil proton is created, it will first pass through the first collimator, that clearly defines the sample active area, then through the second collimator to be detected in the E chamber. The segmented in 64 pads detection plane will collect each signal, allowing the reconstruction of the proton's track. Hence, only those activating pads in the E and E chambers, passing through both collimators, are considered as correct. This allows to reject parasitic tracks.

Acquisition system

To get the proton tracks, the signals collected by the micromegas pads are treated by three successive cards: Zap, SAM and ZedBoard. The details of these cards are given in [14 (II.B) or 15]. Once treated by these cards, the signals are sent to the acquisition system composed of the Cobo software, developed by the CEA Saclay. Once saved, the events can be read by the CoboFrameViewer software, also developed by the CEA Saclay. As observed in the example of Figure 3, a frame is composed of 64 curves, each one corresponding to one pad response.

Detection system characterization

Several tests have been carried out to develop the detector and check its reliability.

Sensitivity under irradiation

The aim of developing this detector is to be insensitive to gamma and electrons background, unlike silicon detectors. This has been checked with a neutron beam produced by the AIFIRA accelerator of the LP2I. When no H radiator was placed upstream the GPRT, only few tracks were observed by the acquisition system. Indeed, when the threshold, given in the software parameters, was adapted to protons, only very few cosmic rays' events could be detected, none coming from gamma nor electrons. This proves the low sensitivity of the GPRT to background under irradiation. The results of this test have been published in [16].

Homogenization of the electric field

It has been checked that the electric field is stable and homogeneous under irradiation. This is possible thanks to the use of a field cage that can be adapted to the polarization voltages. Apart from minor edge effects, making the amplitude of the signals of the pads of the first and last columns dropping a bit, the field can be qualified as homogeneous. Under irradiation, the gain of the detector has proven to be constant. Indeed, issues of charge accumulation on insulating structures have been solved. The field is now stable under irradiation. The tests carried out and the results are presented in [16].

Tracks reconstruction and discrimination

As shown in Figure 3, each activated pad can be identified thanks to its curve. Hence, each track can be reconstructed with respect to the activated pads. An example of the track of a proton coming from the H-rich foil and passing through both collimators is shown in Figure 4. Thanks to the track reconstruction, it is possible to discriminate the desired proton tracks, coming from the H-rich foil and created by a direct neutron, from any other particle. Indeed, a particle not coming from the Hrich foil may not pass through both collimators, as observed for example in Figure 5 (left) the track of a cosmic ray detected by the GPRT. Moreover, the track of a recoil proton from the H-rich foil but created by a scattered neutron will not have the same length in the GPRT, as the neutron energy will be lower than that of a direct neutron, as observed in Figure 5 (right). Once the track of a wanted particle has been selected, it is possible to reconstruct the third dimension thanks to the electron drift time. Therefore, the reconstruction of the track in three dimensions makes the GPRT a Time Projection Chamber (TPC). However, to have access to an accurate third dimension, the chosen gas inside the chamber must be slow so the electrons drift slowly to the anode. Indeed, one can consider the particle passing through the GPRT instantaneously, given the speed of the particle, and thus ionizing the gas everywhere at the same moment. With this condition, the electrons getting first to the anode are the one created the closest from this electrode, allowing a three dimensions track reconstruction as schemed in Figure 6.

Choice of the gas

The gas needs to have good electric properties to allow the micromegas polarization without breakdowns and to enable an effective electrons multiplication. The chosen H-free gas is a mix gas composed of 70% of N 2 and 30% of CO 2 . First tests were performed with CF 4 . However, this gas is too fast inducing a poor time resolution which does not give access to the third dimension. Hence, the N 2 -CO 2 mixture has been chosen for its low electrons drift velocity and good insulation properties. However, tests are currently carried out to find the optimal N 2 /CO 2 proportions. Indeed, a higher ratio of CO 2 increases the gain (signal-to-background ratio increased) but a too high ratio of CO 2 creates breakdown problems due to a carbon deposit on the GPRT's elements. Investigations are still ongoing.

Tests with a proton micro-beam 4.1 Set up description

After the main characterizations and checks of the GPRT, the quantitative properties have been studied thanks to a proton micro-beam delivered by the AIFIRA facility. The rate of recoil protons during a cross section measurements would be around 1 proton/s. We thus need few protons/s entering the GPRT for our tests. To get such a beam, the initial beam was defocused and passed first through a micro collimator with a diameter of 20 µm standing at the entrance of a dedicated chamber (see Figure 7). An aluminum mylar window was used to separate the vacuum of the accelerator from the gas filled chamber. As both the first collimator and the window scatter protons, a second collimator of 2 cm upstream the detector was used. Downstream the GPRT, a silicon detector was used as reference. Finally, a set-up inside the chamber allows to move the GPRT on two directions so the beam can be sent closer or further from the anode or on the side.

Intrinsic efficiency

All the details of the intrinsic efficiency measurements can be found in [14,15]. The intrinsic efficiency has first been studied with an alpha source delivering 4.5 alphas/s inside the GPRT. The efficiency was thus (99.2 ± 3.7 sys ± 1.1 stat ) %, compatible with an intrinsic efficiency of 100%. This result has been checked with the proton micro-beam and confirmed for a rate of less than 3 protons/s. However, this efficiency decreases when the rate increases, down to 7% for 366 protons/s with a 50% efficiency for 38 protons/s. This problem has then been checked with an electronic bench test and it has been confirmed that it comes from the software or the hardware of the acquisition system. This issue is under investigation.

Energy and spatial resolution

The detection plane has quite large pads (0.5x0.5 cm² for the small ones) but it is still relevant to study the spatial resolution because the proton energy will be inferred via its range. It has been observed that the signal switched from one pad to its neighbor as the beam is shifted by 0.5 mm to the side. This effect is also observed when the primary electrons drift for 2 cm.

The energy resolution has also been studied thanks to the proton micro-beam. No quantitative results have been obtained but the energy resolution seems good enough for the requirements since the proton energy will not be reconstructed with the detected energy on each pad.

Conclusion

The GPRT is an almost-functional prototype that will be used to design the operational detector to measure neutron fluxes for fission cross section measurements. Its aim is to detect recoil protons for incident neutrons energy between 200 keV and 2.5 MeV.

As the GPRT allows the track reconstruction in three dimensions, this detector is a Time Projection Chamber.

The low sensitivity to gamma and electron background has been checked and validated under irradiation. Thanks to a field cage that can be adjusted, the electric field is uniform and stable under irradiation. Some minor edge effects are still present but they only impact the first and last columns of pads. The energy and spatial resolutions have been investigated and present no constraint for our use. The intrinsic efficiency has been studied and is 100% for low protons rates and decreases when the rate increases. An investigation is carried out to improve the acquisition system performances.

Finally, once the intrinsic efficiency issue is solved, a new GPRT will be designed with (i) a more robust behavior with respect to breakdowns issues since they prevent a good reliability of the GPRT, (ii) a better design geometry to reduce the uncertainties of the geometric efficiency and (iii) the number of pads may be increased to get a better track reconstruction.

Fig. 1 .Fig. 2 .

 12 Fig. 1. Picture of the Gaseous Proton-Recoil Telescope

Fig. 3 .

 3 Fig. 3. Example of a frame saved by the Cobo software. Each curve corresponds to one pad.

Fig. 4 .

 4 Fig. 4. Example of a track reconstructed based on the frame's curves.

Fig. 5 .

 5 Fig. 5. Example of a track reconstructed with respect to the frame's curves for a particle not coming from the H-rich foil (left) or coming from the H rich foil but from a scattered neutron (right).

Fig. 6 .

 6 Fig. 6. Reconstruction of the third dimension thanks to the electron drift velocity

Fig. 7 .

 7 Fig. 7. Scheme of the designed chamber designed for this experiment.

Acknowledgements: The authors would like to thank the AIFIRA team work, who helped for all the experiments, and the LP2I's electronic service, who is in help with all the acquisition's issues. This project has also received funding from the Euratom research and training program 2014-2018 under Grant No. 847552 (SANDA).


Julien Gacon

Jannes Nys

Riccardo Rossi

Stefan Woerner

Giuseppe Carleo

Variational Quantum Time Evolution without the Quantum Geometric Tensor

The real-and imaginary-time evolution of quantum states are powerful tools in physics, chemistry, and beyond, to investigate quantum dynamics, prepare ground states or calculate thermodynamic observables. On near-term devices, variational quantum time evolution is a promising candidate for these tasks, as the required circuit model can be tailored to trade off available device capabilities and approximation accuracy. However, even if the circuits can be reliably executed, variational quantum time evolution algorithms quickly become infeasible for relevant system sizes due to the calculation of the Quantum Geometric Tensor (QGT). In this work, we propose a solution to this scaling problem by leveraging a dual formulation that circumvents the explicit evaluation of the QGT. We demonstrate our algorithm for the time evolution of the Heisenberg Hamiltonian and show that it accurately reproduces the system dynamics at a fraction of the cost of standard variational quantum time evolution algorithms. As an application of quantum imaginary-time evolution, we calculate a thermodynamic observable, the energy per site, of the Heisenberg model.

I. INTRODUCTION

Quantum time evolution is a central task in physics. Real-time evolution provides detailed insight into properties of quantum mechanical systems, such as phase transitions [1][2][3] or thermalization [4,5]. Imaginary-time evolution is an important tool that enables the preparation of ground states or thermal states [6][7][8]. These can, in turn, be used for the calculation of thermodynamic observables [8,9]. In particular, combining real-and imaginary-time evolution would allow the direct calculation of dynamical correlation functions at thermal equilibrium.

The range of applications of imaginary-time evolution extends beyond the field of physics. Ground-state preparation with imaginary-time evolution for gapped, nondegenerate Hamiltonians is guaranteed to converge in the generic case of non-zero overlap between the ground state and the initial trial state. This makes it a promising candidate in settings where a good initial state can be constructed, e.g. in chemistry applications [10] or in classical optimization problems [11]. In quantum machine learning, the preparation of Gibbs states with imaginary-time evolution is a subroutine for quantum Boltzmann machines, which can, for example, be used in distribution learning or classification [12].

Since performing quantum time evolution generally requires representing the exponentially large wave function of a quantum system, quantum computers are a promising platform for developing efficient algorithms [13]. In fact, in 1996, the Trotter algorithm for real-time evolution was among the first proposed use cases for a quantum computer [14]. However, the complexity of the quantum circuits required for the Trotter algorithm depends on the Hamiltonian, and the circuit depth scales with the simulation time and accuracy [15,16]. This renders the algorithm currently unsuitable for general time evolution on near-term devices, which are characterized by limited qubit connectivity and coherence times. The imaginary- time counterpart of Trotter suffers from the same restriction [8].

Variational algorithms for quantum time evolution, on the other hand, allow to choose a parameterized circuit as an ansatz to approximate the wave function, that operates within the device's capabilities. Using a variational principle, variational quantum time evolution (VarQTE) maps the quantum state evolution to the evolution of parameters in the model [17], both for real -time evolution (VarQRTE) and imaginary-time evolution (VarQITE). The parameter update rules depend on the evaluation of the Quantum Geometric Tensor (QGT) and gradients of the current energy and state. For an ansatz with d variational parameters, the number of circuits required to evaluate the QGT and gradients scale as O(d 2 ) and O(d), respectively. While this does not pose a problem for small systems, the evaluation of the QGT quickly becomes a bottleneck once the system size, and therefore the number of variational parameters, increases.

Figure 1 shows a runtime estimation for VarQITE, assuming current superconducting processor specifications, see Appendix A for details on the derivation. For a few parameters the runtime is of the order of hours, but already for only 200 parameters the computation time around 1 week, which renders this algorithm currently impractical. With recent advances in processor sizes exceeding 100 qubits, such as the IBM Quantum Eagle or Osprey devices [START_REF]IBM Quantum[END_REF][START_REF]IBM Unveils 400 Qubit Plus Quantum Processor and Next Generation IBM Quantum System Two[END_REF], improving the resource requirements of quantum algorithms becomes crucial for finding practically relevant applications of quantum computers.

Recently, focus has shifted to optimization-based algorithms, which implement partial steps or approximations of the full Suzuki-Trotter step [START_REF] Barison | [END_REF][21][22][23][24]. In the case of realtime evolution, for example, the projected variational quantum dynamics (p-VQD) algorithm [START_REF] Barison | [END_REF] provides a scalable alternative to VarQRTE on near term-devices, if a single Trotter step can be efficiently implemented. However, the required quantum circuit gates in p-VQD reflect the couplings of the Hamiltonian. This means that, for Hamiltonians with long-distance interactions or numerous Pauli terms (e.g. in molecular dynamics), even a single step could involve global connections or deep circuits hindering the execution on near-term devices. Furthermore, the p-VQD algorithm is not directly applicable to imaginary-time evolution.

Other approaches concerned with real-time evolution are Variational Fast Forwarding (VFF) methods [25][26][START_REF] Gibbs | Long-time simulations with high fidelity on quantum hardware[END_REF][START_REF] Lim | [END_REF] and classical pre-processing approaches [29,30]. VFF methods rely on diagonalizing the Hamiltonian or the Trotterized time evolution operator with a variational ansatz. However, finding the diagonalizing unitary remains challenging in practice, which limits demonstrations to very few qubits. Classical pre-processing techniques, on the other hand, impose additional restrictions on the simulated system, such as translational invariance [30] or Hamiltonians with low entanglement [29]. Within such systems, these techniques scale to large systems, but they do not allow for general quantum time evolution.

Another line of work directly focuses on the preparation of thermal states by minimizing the free energy of a variational ansatz [31]. This approach, however, also does not implement general quantum time evolution.

In this paper, we propose a novel variational algorithm for quantum time evolution based on a dual optimization problem, which allows to replace the QGT by evaluating the overlap of the variational ansatz for different parameter values. This formulation applies equally real-and imaginary-time evolution and does not require additional qubits or connections than already present in the ansatz. We show that this new algorithm requires significantly fewer measurements and thereby drastically reduces the expected runtime compared to VarQTE. This is summarized in Fig. 1, where, under the same assumptions, our proposed method can reduce the expected runtimes from several weeks for VarQTE to only a few days. Following the naming conventions of VarQTE, we name the algorithm DualQTE with specifiers DualQITE for imaginarytime evolution and DualQRTE for real-time evolution.

The remainder of this paper is structured as follows. In Sec. II, we recap VarQTE based on variational principles, derive the proposed dual formulation, and discuss how to implement it on a quantum computer. Then, in Sec. III, we demonstrate our proposed algorithm for the imaginary-time evolution of the Heisenberg model and investigate the resource requirements. As a practical application, we use the quantum minimally entangled typical thermal states method (QMETTS) to calculate thermodynamic observables. Sec. IV demonstrates the dual formulation for real-time evolution, including the calculation of variational error bounds. Finally, Sec. V concludes the paper and gives an outlook on possible applications and further research directions.

II. DUAL FORMULATION OF VARIATIONAL TIME EVOLUTION

For a time-independent Hamiltonian H acting on n qubits, an initial quantum state |Ψ 0 ⟩ and an evolution time t, the real-time evolved quantum state is

|Ψ(t)⟩ = e -itH |Ψ 0 ⟩ . (1)
For an imaginary-time evolution, the time evolution operator is non-unitary, and the normalized state reads

|Ψ(t)⟩ = 1 ⟨Ψ 0 |e -2tH |Ψ 0 ⟩ e -tH |Ψ 0 ⟩ . (2)
Variational quantum time evolution maps the evolution of the quantum state |Ψ(t)⟩ to the evolution of parameters θ(t) ∈ R d of a parameterized quantum state |ϕ(θ(t))⟩. The parameters' dynamics can be derived with variational principles such as the Dirac-Frenkel, McLachlan, or time-dependent variational principle [17]. In McLachlan's formulation, the derivatives of the parameters are determined by the linear system of equations

g(θ(t)) θ(t) = b(θ(t)), (3)
where the matrix g = Re(G) ∈ R d×d is the real part of the QGT, and we call b ∈ R d the evolution gradient. The QGT is defined as

G ij (θ) = ⟨∂ i ϕ(θ)|∂ j ϕ(θ)⟩ -⟨∂ i ϕ(θ)|ϕ(θ)⟩ ⟨ϕ(θ)|∂ j ϕ(θ)⟩ , (4
) where we use the notation ∂ i := ∂/(∂θ i ) and do not explicitly state the time dependence of the parameters. The evolution gradient for VarQRTE is given by the expression

b R i (θ) = Im ⟨∂ i ϕ(θ)|H|ϕ(θ)⟩ -⟨∂ i ϕ(θ)|ϕ(θ)⟩ E(θ) , (5)
whereas a VarQITE evolution yields the following

b I i (θ) = -Re ⟨∂ i ϕ(θ)|H|ϕ(θ)⟩ = - ∂ i E(θ) 2 , (6)
with the energy E(θ) = ⟨ϕ(θ)|H|ϕ(θ)⟩. From hereon, we present general equations that apply to both real and imaginary-time evolution; thus, unless specified, we simply use b without a specific superscript.

Note that these equations are introduced for a timeindependent Hamiltonian, but they can also be applied to the time-dependent case H = H(t).

A. Dual formulation

Instead of solving the linear system defined in Eq. ( 3), we propose to solve the dual formulation of the problem [32,33] given by θ = argmin

θ θT g(θ) θ 2 -θT b(θ). (7)
The term ∥ θ∥ 2 g(θ) = θT g(θ) θ is the squared norm of the parameter derivative in the metric of the QGT. This quantity describes the magnitude of the derivative from an information geometric point of view and is derived from the Fubini-Study metric. For infinitesimal displacements δθ, we have

||δθ|| 2 g(θ) = δθ T g(θ)δθ = 1 -| ⟨ϕ(θ)|ϕ(θ + δθ)⟩ | 2 + O(∥δθ∥ 3 2 ), (8)
where ∥ • ∥ 2 is the ℓ 2 norm [33]. By writing θ = δθ/δτ , for some time perturbation δτ > 0, we can now reformulate the optimization in terms of the fidelity

F (θ, θ ′ ) = | ⟨ϕ(θ)|ϕ(θ ′ )⟩ | 2 as δθ ≈ argmin δθ 1 -F (θ, θ + δθ) 2(δτ ) 2 - δθ T b(θ) δτ = argmin δθ L(δθ) (δτ ) 2 , (9)
where we directly optimize for the parameter update δθ and we introduced the loss function

L(δθ) = 1 -F (θ, θ + δθ) 2 -δτ • δθ T b(θ). (10)
In practice, the optimization problem can be solved without the factor (δτ ) -2 , which decouples the shape of the locally quadratic infidelity term from the time perturbation and improves the numerical stability of the optimization.

Note that this dual formulation can alternatively be obtained from the derivation of quantum natural gradients [31,33], which is detailed in Appendix B. For an intuitive understanding of the relationship of the infidelity and QGT the effect of approximating ∥δθ∥ g(θ) ≈ 1-F (θ, θ+δθ) in an illustrative example is demonstrated in Appendix C.

Instead of computing the QGT at each timestep, which requires O(d 2 ) circuit evaluations, we now have to solve an optimization problem where the loss function requires only one fidelity evaluation. The required resources of DualQTE per timestep are therefore O(d) for the computation of the evolution gradient b, times the number of iterations in the optimization. Thus, we improve upon the direct QGT approach if the number of iterations scales better than O(d), which, as we show in the following sections, is the case for the examples we investigate in this work.

B. Evaluating the loss function

The evaluation of the loss function L, defined in Eq. (10), requires the calculation of the evolution gradient b and the fidelity of the ansatz |ϕ(θ)⟩ for two different parameter sets. For imaginary-time evolution, the evolution gradient can, for example, be evaluated with analytic gradient rules, such as the parameter-shift rule or a linear combination of unitaries (LCU), or with finite difference methods [34]. In the case of real-time evolution, however, we are restricted to an LCU approach, as this is the only method that allows the calculation of the imaginary part of gradients [12].

The fidelity F can, for example, be estimated using the swap test [35] and its variants [36], where the states are prepared in separate qubit registers followed by entangling gates across these registers, or with the Hadamard test, which adds only a single auxiliary qubit, but requires controlling the state-preparing unitary [37]. A more near-term-friendly option is the compute-uncompute method [38], which does not introduce additional global operations. If the states are given by |ϕ(θ)⟩ = U (θ) |0⟩ for a parameterized unitary U and two different parameter values θ and θ ′ , the fidelity can be calculated by preparing U † (θ)U (θ ′ ) |0⟩ and measuring the probability of obtaining |0⟩ on all qubits.

If the state |ϕ(θ)⟩ has n qubits and the preparing unitary U has depth m, the swap test variants require a circuit width of 2n with depth of m+O(1), whereas the evaluated circuits for the compute-uncompute method are of only width n, but of depth 2m. The Hadamard test for fidelities between the same circuit with different parameters can be evaluated by controlling the parameterized gates, resulting in depth of m and depth of n+1, plus the overhead of controlling the gates. For sparse device connectivities, this can be a challenge. To avoid increasing the circuit complexity, the overlap can also be estimated via randomized measurements of two independent state preparations [39]. However, this technique requires an exponential number of measurements.

Evaluating the QGT for VarQTE, however, suffers from similar issues. The QGT can be evaluated as the Hessian of the infidelity [40] using a parameter-shift or finite difference technique, which comes with the restrictions for fidelity evaluations as described above. Alternatively, Eq. ( 4) can be directly computed with an LCU approach, which adds two auxiliary qubits and two entangling gates [12]. This method is less demanding than e.g. a Hadamard test, but still comes with additional connectivity requirements. In practice, for both VarQTE and DualQTE a suitable combination of parameterized quantum state |ϕ(θ)⟩ and gradient method must be selected, such that the resulting circuits can be executed reliably.

Depending on the topology and coherence times of the available hardware and the structure and size of the unitary, either method for gradient and fidelity calculations can be advantageous. In this work, we focus on near-term friendly methods and use the parameter-shift rule for gradients (if possible) and the compute-uncompute method for the fidelity, as these do not require additional gate connections or an exponential number of measurements. Note that, for systems with a large number of qubits, this method might become unsuitable as it measures the global zero projector. Then, approaches using only local measurements, such as the Hadamard test, could be the better choice.

C. Solving for the update step

The infidelity-based loss function L is locally convex around δθ = 0, as its Hessian at this point is ∇∇ T L(0) = g/2, and g is positive semi-definite. To leverage this property, we use gradient descent as a local optimization routine, which also allows the use of analytic gradient formulas that have proven more stable in presence of shot noise. The gradient of L with respect to the parameter update δθ is

∇ δθ L(δθ) = - ∇ δθ F (θ, θ + δθ) 2 -δτ • b(θ).
The gradient of the fidelity can be evaluated with a parameter-shift rule

∂F ∂(δθ) i = F (θ, θ + δθ + e i s) -F (θ, θ + δθ -e i s) 2 sin(s)
, where e i is the i-th unit vector and s is the parameter shift, which can be chosen as, e.g., π/2 for single-qubit Pauli rotations [34].

At each timestep, the gradient descent update for the update step δθ is

δθ (k+1) = δθ (k) -η k ∇L δθ (k) ,
where η k > 0 is the learning rate at step k. This iteration is continued until a maximum number of iterations or a convergence criterion is met. An example of the latter is a minimum tolerance in the change of the cost function or the norm of the gradient.

An intuitive choice for the initial guess δθ (0) is the zero vector, which corresponds to no change in the parameters. However, a more efficient choice can be to introduce momentum by warm starting the optimization with the update step from the previous timestep. This heuristic is motivated by the idea that, especially for small timesteps, we do not expect the parameter derivatives θ to change significantly.

Methods that approximate the gradient, such as finite difference or SPSA [41], may face challenges in the optimization. For small timesteps, the fidelity is close to 1 and the noise in the readout, e.g. from finite sampling statistics or device noise, can easily mask changes in the cost function. Parameter-shift gradients suffer less from this problem, as they allow to evaluate the cost function over larger perturbations, and do not amplify the noise by dividing by a small constant.

The ideal choice of the time perturbation δτ is a tradeoff: the error in approximating the QGT scales with (δτ ) 3 , but a smaller perturbation amplifies any measurement noise in the loss function as the update step is obtained as δθ/δτ . Appendix C displays this trade-off for an illustrative example.

D. Trainability

Recently, there has been a lot of research showing that, in certain settings, the loss function gradients of variational algorithms decay to zero exponentially and cannot be evaluated efficiently, as they would require an exponential number of measurements. These so-called barren plateaus can be induced, for example, if the loss function requires measuring a global observable [42], if the quantum circuit preparing the parameterized state is too deep or generates too much entanglement [42][43][44], or if the measurements are too noisy [45].

Since variational quantum dynamics is driven by the evolution gradient defined Eqs. ( 5) and ( 6), it can be affected by a barren plateau and fail to track the true evolution of the quantum state. However, it is important to note that the gradients only vanish on average for a random initialization, whereas in time-evolution the initial quantum state is typically specifically chosen. Furthermore, Hamiltonians of physical systems are usually local, as they reflect the interactions of the quantum mechanical system, and exponentially vanishing gradients can be avoided by choosing a circuit depth scaling logarithmically in system size [42]. Alternatively, an application-specific ansatz with few variational parameters can help mitigate barren plateaus, such as circuits based on Hamiltonian evolutions [46,47].

In addition to the evolution gradient, the DualQTE loss function gradient ∇ δθ L depends on the gradient of the fidelity, which relies on measuring a global observable. This can be seen by writing the fidelity of two n-qubit states prepared by unitaries U (θ) and

U (θ ′ ) as | ⟨0|U † (θ)U (θ ′ )|0⟩ | 2 = ⟨λ|P 0 |λ⟩, where |λ⟩ = U † (θ ′ )U (θ) |0⟩ and P 0 = |0⟩ ⟨0|
⊗n is the global projector on the all-zero state. Thus, evaluating the fidelity gradient for two randomly selected parameter sets θ and θ ′ would exhibit barren plateaus at any circuit depth [42]. However, the optimization in DualQTE starts at zero perturbations, θ = θ ′ where the total state preparing unitary is the identity, |λ⟩ = U † (θ)U (θ) |0⟩ = I |0⟩, which is an initialization that is proven to not exhibit barren plateaus even for global cost functions [48]. Together with the fact that the DualQTE loss function is locally convex, the non-vanishing gradients at the initial point of the optimization is a strong motivation for the efficient trainability of DualQTE.

In Appendix E 4, we provide numerical evidence that for a local Hamiltonian and a logarithmic-depth circuit, neither the evolution gradient or the fidelity gradients decay exponentially with system size.

E. Sample complexity

The implementation of VarQTE on quantum hardware has several sources of errors: the model |ϕ(θ)⟩ could lack expressitivity to capture the dynamics, the time integration scheme introduces errors, the QGT and evolution gradient are subject to sampling error from a finite number of measurement, and each operation is affected by hardware noise. If we denote the ideal VarQTE parameters without sampling or hardware noise by θ(t) and the noisy parameters by θ(t), the error contributions can be split as

ε(t) = D B (ϕ( θ(t)), Ψ(t)) ≤ ε M (t) + ε S (t), (11)
where we measured the error in Bures distance

D B (ψ, ϕ) = 2(1 -| ⟨ψ|ϕ⟩ |), (12)
and distinguish in error due to lack of model expressitivity plus integration error ε M (t) = D B (ϕ(θ(t)), Ψ(t)) and error due to a noisy implementation of VarQTE ε S (t) = D B (ϕ( θ(t)), ϕ(θ(t)) [49].

Since the proposed DualQTE algorithm promises a reduction in measurement cost of VarQTE, but is not concerned with the ansatz selection or hardware noise, we here focus on investigating the scaling of the sampling error. The model error ε M can be bounded with a-posteriori errorbounds [50], which we also investigate for the real-time evolution case in Sec. IV.

Deriving a concrete bound in terms of system quantities such as the energy or the number of parameters requires an assumption on the circuit structure. Here we assume a circuit with only Pauli rotations R P (θ i ) where each parameter θ i is unique and does not have coefficients. Note that the bounds can be adjusted for different circuit structures and parameterizations. In addition, we assume a cutoff δ c > 0 on the smallest eigenvalue of g due to a regularization of the linear system.

Then, we can state the following upper bound on the number of samples required to achieve a sampling error of ε S ,

N ≤ O d 3 E 2 max ∆ 2 t δ 4 c ε 2 S , (13)
where E max is the maximal eigenvalue of the Hamiltonian.

In contrast to a similar approach described in Ref. [49], we state the upper bound in terms of the number of parameters in the model or the system's Hamiltonian, instead of the QGT and evolution gradients. Further, we are able to derive a tighter result by leveraging Latala's theorem from random matrix theory to upper bound the sampling error in g [51].

Since the DualQTE algorithm does not construct the QGT directly but only the evolution gradient, we expect a reduction of a factor d in the complexity, and an additional factor for the number of optimization steps K in each timestep. Indeed we can show that the upper bound for the number of samples is

N ≤ O d 2 K 2 ∆ 2 t δτ 2 ε 2 S 1 δτ + E max 2 . ( 14
)
The detailed derivation of both bounds is described in Appendix D. While it is possible to construct circuits where each component of these bounds are tight Sec. III shows that in practice the actual number of required samples scales less than this upper bound, which is further discussed in the Appendix.

III. IMAGINARY-TIME EVOLUTION

In this section, we show the results of DualQITE and investigate the circuit costs compared to VarQITE. As an application, we use our algorithm as a subroutine to prepare typical thermal states of the Heisenberg model, which are then used to calculate the energy per site as a thermodynamic observable. All circuits are constructed and simulated using Qiskit [START_REF]Qiskit: An open-source framework for quantum computing[END_REF].

A. Heisenberg model

We simulate the imaginary-time evolution of the Heisenberg model with nearest-neighbor interaction on a 12-qubit circle in a transverse field:

H = J ⟨ij⟩ (X i X j + Y i Y j + Z i Z j ) + g i Z i , (15)
with interaction strength J = 1/4 and transversal field strength g = -1. As a variational ansatz, we use a circuit with Pauli-Y and Pauli-Z single qubit rotation layers that alternate with pairwise CNOT entangling gates. 16)) for DualQITE and VarQITE, with mean and standard deviation of 5 experiments. The resources are measured in total number of measurements and are shown for evaluation of gradients (and the QGT, in the case of VarQITE) using parameter-shift rules (dashed) or a LCU approach (solid lines).

The circuit structure is shown in Appendix E 1 and we use r = 3 repetitions. The initial state for the evolution is the equal superposition of all qubits, |+⟩ ⊗n , which we prepare by setting the rotation angles of the last Pauli-Y layer to π/2 and the remaining angles to 0.

The optimization problems in DualQITE are solved using gradient descent with a fixed learning rate of η = 0.1 and time perturbation δτ = 0.01. The initial iteration performs 100 update steps and the subsequent, warmstarted iterations, only 10. These values are motivated by the simulations shown in the Appendix E 2 and are partially heuristic, as a termination criterion is challenging to define with access only to noisy loss function and gradient evaluations. The parameters are integrated with an explicit Euler scheme with timestep ∆ t = 0.01, i.e.

θ(t + ∆ t ) = θ(t) + ∆ t θ(t) = θ(t) + ∆ t δθ δτ .
Note that the integration timestep ∆ t , which determines the accuracy of the time integration, can be chosen differently from the time perturbation δτ , which affects the approximation error of the QGT metric with the infidelity.

We compare the performance to VarQITE with the same integration scheme and use an L-curve regularization [START_REF] Cultrera | [END_REF] for a stable solution of the linear system. Among all regularization techniques we attempted, such as adding a diagonal shift, truncating small or negative singular values or solving on a stable subsystem, the Lcurve regularization provided the most accurate and stable results.

In Fig. 2(a), we present the results for a varying number of shots along with the exact time evolution based on exact diagonalization. Already with as little as 100 mea-surements per circuit evaluation (shots) on the 12-qubit model, the dual time evolution is able to qualitatively follow the imaginary-time evolution and, up to time t ≈ 1, even outperform VarQITE with 1024 shots. Increasing the number of measurements of DualQITE to 1024 shots allows the dual method to closely track the exact solution towards the ground state, with a higher accuracy than VarQITE with 8192 shots.

B. Resource requirements

In the above experiment, DualQITE requires fewer circuit evaluations to achieve the same accuracy as Var-QITE. To investigate the total resource requirements, we perform both DualQITE and VarQITE with different resources and compute the achieved error. Since we are interested in following the imaginary-time dynamics as closely as possible at each timestep, we define the error as the average integrated Bures distance to the exact solution over the time evolution,

I B (T ) = 1 T T 0 D B (ϕ(θ(t)), ψ(t))dt. ( 16
)
The state fidelity is computed exactly, i.e., we compute the state vector of the model |ϕ⟩ at variational parameters θ(t), and take the inner product with the exact time-evolved state |Ψ(t)⟩.

The results for an integration time of T = 2 are shown in Fig. 2(b). We show the integrated Bures distance with respect to the total number of measurements recorded during the time evolution. In DualQITE, the resources can be split between using more optimization steps in each timestep or more shots to evaluate the gradients. The algorithm settings are detailed in Appendix E 3. The figure shows the resource counts for gradient calculations via the parameter-shift rule (PSR) and linear combination of unitaries (LCU). The LCU technique requires additional auxiliary qubits and additional non-local operations, but less overall circuits than PSR. For P Pauli terms in the Hamiltonian, the total number of required circuits C per timestep

C VarQITE LCU = d(d + 5) 2 + P d C VarQITE PSR = 2d(d + P + 1). ( 17
)
For DualQTE the number of circuits is

C DualQITE LCU = P d + Kd, (18)
and

C DualQITE PSR = 2C DualQITE LCU
, where K is the number of optimization steps per timestep. The total number of measurements N is obtained by multiplying the number of circuits with the number of shots.

We see that, on average, DualQITE requires about one order of magnitude fewer measurements to achieve the same accuracy as VarQITE. With an increasing number of parameters, we expect this difference to grow, since VarQITE scales as O(d 2 ) whereas our algorithm, with warm starting, only computes small corrections at each time step.

C. Sample complexity

In addition to the fixed-size model with 12 qubits, we investigate how the resource requirements scale with system size. We compare VarQITE and DualQITE for the Heisenberg model from Eq. ( 15) with varying number of spins n and the same circuit structure as before, but with an adjusted number of repetitions of r = ⌈log 2 (n)⌉ times, plus a final rotation layer. We then tune the settings of VarQITE and DualQITE to achieve a mean accuracy of I B ≤ 0.1 over 5 experiments and count the total number of required measurements N . This threshold corresponds to a per-timestep fidelity of 0.995.

The results are presented in Fig. 3, which show the improved scaling of DualQITE compared to VarQITE. For small system sizes and few parameters, the overhead of solving the optimization problem in DualQITE is larger than evaluating the QGT. But, as we increase the problem size, the quadratic scaling of VarQITE takes over and our algorithm becomes more efficient.

This experiment allows to validate the upper bound on the number of measurements of Sec. II E. As shown in Fig. 3, the model error ε M is negligible in comparison to the sampling error ε S and we approximately have ε S ≈ I B ≈ 0.1. The maximal energy of the Heisenberg model on a periodic chain scales with the number of spins, and can be bounded by

E max = O(n) ≤ O(n log(n)) = O(d).
Inserting these values in Eqs. ( 13) and ( 14) we expect the scaling to be upper bounded by O(d 5 ) for VarQTE and O(d 4 K) for DualQTE. In practice, we observe approximately a scaling of d 3.56 for VarQTE and d 2. 29 for DualQTE, which shows the expected improved scaling for our algorithm. The measured scaling also suggests that the bounds are not yet tight, which we discuss further in Appendix D.

D. Calculating thermodynamic observables

As an application of imaginary-time evolution, we calculate thermodynamic observables using the quantum minimally entangled thermal states algorithm (QMETTS) [8,54]. While the classical METTS algorithm has been specifically developed for Matrix Product State simulations, the thermal state preparation is still costly, and classical simulations fail if the system produces macroscopic entanglement during the imaginary time evolution (e.g. for low-temperature, 2D systems). Due to these restrictions, QMETTS is a promising application for quantum imaginary-time evolution algorithms.

For an observable A and inverse temperature β, the QMETTS algorithm generates samples {A m } m using a Markov chain whose average approximate the ensemble average:

⟨A⟩ ens = Tr(e -βH A) Tr(e -βH ) ≈ 1 M M m=1 A m .
The sampling process to obtain the sample A m is 1. Start from a product state |ϕ m (t = 0)⟩.

Evaluate the observable to obtain the sample

A m = ⟨ϕ m (β/2)|A|ϕ m (β/2)⟩ .
4. Measure |ϕ m (β/2)⟩ in some basis to obtain the next random product state |ϕ m+1 (0)⟩.

We investigate the Heisenberg model from Eq. ( 15) on a chain with n = 6 spins with parameters J = 1/4 and g = -1. As a thermodynamic observable we compute the energy per site, ⟨H⟩ /n. To reduce the auto-correlation length in the QMETTS Markov chain, and for faster convergence to the ensemble average, it is favorable to measure in different bases in each step. Since the Heisenberg Hamiltonian conserves the number of qubits in the |1⟩ state, avoiding the Z basis greatly reduces the standard deviation of the Markov chain. Thus, we here alternate between the X and Y basis for each sample.

As ansatz for DualQITE, we use problem-inspired circuits with pairwise CNOT couplings and r = 2 repetitions of rotation and entanglement layers, plus a final rotation layer, see Appendix E 1 for a circuit diagram. For evolutions of product states |±⟩ in the X basis, the rotation layers are single qubit R Y R Z gates, and for the states |±i⟩ in the Y basis, the layers implement R X R Z gates. The initial product states |ϕ m (0)⟩ are prepared by setting the parameters in the final layer rotation layer of the ansatz as follows,

|±⟩ → R Y ±π 2 R Z (0), |+i⟩ → R X π 2 R Z (π), |-i⟩ → R X π 2 R Z (0).
Each energy sample is evaluated with 1024 measurements per basis. The optimization problem in DualQITE is solved with a time perturbation δτ = 0.01 and gradient descent with a learning rate of η = 0.1 and 100 iterations in the first timestep, followed by 10 iterations in the following, warmstarted timesteps. We integrate with a fixed timestep of ∆ t = 0.01. Figure 4 shows the estimated energy per site, along with the standard deviation of the samples, for different inverse temperatures β. For the alternating X -Y basis, the Markov chain converges quickly and M = 25 samples suffice for an accurate estimate of the observable. For the imaginary-time evolution, we compare DualQITE with the same settings as in the previous sections to an exact evolution performed with matrix exponentials. It shows that using the dual method allows to reliably reproduce the mean and standard deviation of the Markov chain samples compared to the exact reference METTS.

IV. REAL-TIME EVOLUTION

The focus of this paper is on imaginary-time evolution as, to date, no other QGT-free time evolution algorithms exist in this setting. For real-time evolution, p-VQD has a similar structure as our algorithm and solves an optimization problem rather than evaluating the QGT. However, there are key differences to the dual algorithm applied to real-time evolution.

The p-VQD algorithm [START_REF] Barison | [END_REF] projects a single Suzuki-Trotter step onto the circuit model by solving the following optimization problem:

θ(t + ∆ t ) = argmax θ ′ ⟨ϕ(θ ′ )|e -iH∆t |ϕ(θ(t))⟩ 2 .
For Hamiltonians with many Pauli terms or long-range interactions, such as those arising in molecular dynamics, the single step might already lead to large circuits with non-local gates. While DualQRTE requires an LCU method to evaluate the imaginary part of the energy and state gradients, see Eq. ( 5), this only adds a single entangling gate, compared to a full Suzuki-Trotter step is required. Furthermore, our dual time evolution allows the evaluation of error bounds at almost no additional cost, which is not possible in p-VQD. Due to these differences, this section presents DualQRTE: the dual time evolution for real-time evolution.

A. Heisenberg model

We present the real-time evolution under the Heisenberg Hamiltonian of Eq. ( 15) on a linear chain with n = 4 spins with parameters J = 1/4, g = -1. As variational model, we use a circuit with alternating Pauli-X and Pauli-Y rotation layers, and Pauli-ZZ entangling gates that reflect the connectivity of the spins. The circuit structure is visualized in Appendix F and, in this experiment, all algorithms use r = 3 repetitions of the rotation as well as entangling gates. To prepare the initial state, |+⟩ ⊗4 , we set the parameters of the final Pauli-Y rotations to π/2 and the rest to 0.

During the evolution, we track the average magnetization in the X and Z direction,

⟨X⟩ = 1 n n i=1 ⟨X i ⟩ , ⟨Z⟩ = 1 n n i=1 ⟨Z i ⟩ .
Since this Heisenberg Hamiltonian preserves the qubit excitations, and the initial state is the equal superposition, the ⟨Z⟩ expectation value should remain 0 throughout the evolution. The results of the different time evolution algorithms for an integration time of T = 2 and timestep ∆ t = T /100 are presented in Fig. 5. Both DualQRTE and p-VQD accurately track the observables using only 200 shots per circuit. With the same resources, VarQRTE, on the other hand, has lower accuracy and we need to use 1024 shots per circuit to match the result of the optimization-based algorithms.

B. Error bounds

In variational real-time evolution, the model error in terms of Bures distance D B due to restriction to the variational manifold can be expressed as [50] εM :=

d k=1 θk |∂ k ϕ(θ)⟩ + iH |ϕ(θ)⟩ 2 2 = Var(H|ϕ(θ)) + θT g(θ) θ -2 θT b(θ), ( 19
)
where we set ℏ ≡ 1. Integrating this error rate provides an upper bound on the Bures distance, that is

D B (ϕ(θ(T )), Ψ(T )) ≤ T 0 εM (t)dt,
where |Ψ(t)⟩ is the exact time-evolved state and the timedependence of ε M is due to the time-dependence of the parameters θ = θ(t). 2 , this error is proportional to the loss function used in DualQRTE. By using the same expansion θ = δθ/δτ and using the infidelity to approximate the inner product with respect to the geometric tensor, we can rewrite the error as

εM = Var(H|ϕ(θ)) + 1 -F (θ, θ + δθ) (δτ ) 2 - 2δθ T b(θ) δτ + O(δτ ) = Var(H|ϕ(θ)) + 2L(δθ) (δτ ) 2 + O(δτ ).
Note that the error scales linearly in time perturbation δτ as the infidelity approximation has a cubic error term [33], which is divided by the square of the perturbation. If we, for example, use a forward Euler rule with timestep ∆ t to integrate the variational error, the integration error scales as O(∆ t + T δτ ). This highlights the importance of differentiating between the timestep ∆ t for the integration, and the time-perturbation δτ to approximate the derivative.

In Fig. 6, we show the error bounds along with the true error for the time evolution of the Heisenberg model. The bounds are computed for different timesteps ∆ t for VarQRTE, and for DualQRTE for a fixed time perturbations δτ = 10 -3 in exact simulations. Firstly, we can verify that the error bounds hold. Secondly, the larger the timestep relative to the time-perturbation, the more accurate the approximation of the dual time evolution, as the error O(∆ t + T δτ ) is dominated by the integration error.

V. CONCLUSION

In this paper, we present a novel algorithm for variational quantum time evolution that does not require the evaluation of the QGT, but instead solves a dual optimization problem in each timestep. The proposed dual time evolution algorithm, DualQTE, is particularly interesting for imaginary-time evolution, as there is currently no alternative variational algorithm able to circumvent the O(d 2 ) cost of VarQITE. For real-time evolution, p-VQD also offers an optimization-based approach by projecting a single Trotter step onto the variational form. In comparison, the dual time evolution has the advantage that no Suzuki-Trotter step has to be implemented, which could require deep circuits or non-local operations, depending on the Hamiltonian. Furthermore, our algorithm allows to evaluate variational error bounds [50], although how accurately they can be evaluated in the presence of shot noise remains an open question.

We demonstrated DualQTE for the imaginary-time evolution of a Heisenberg Hamiltonian on 12 qubits, and found that, in this setting, it requires about one order of magnitude less measurements to achieve the same accuracy as VarQITE. As a practical application of imaginary-time evolution, we calculated thermodynamic observables with the QMETTS algorithm and showed that the DualQITE is suitable to reproduce the sampling distributions. Finally, we applied our algorithm to an illustrative example for real-time evolution, where it produced comparable results to p-VQD for the same amount of resources, while both algorithms outperformed VarQRTE.

In the presented experiments, we used standard gradient descent algorithms with a fixed learning rate. We expect that the performance could be further improved by using more advanced optimization schemes, or methods that also take into account information from previous iterations. Another possible improvement would be a suitable termination criterion for noisy evaluations of the loss function. As for other optimization-based time evolution algorithms, such as p-VQD, it remains challenging to accurately measure the fidelity in the presence of hardware noise.

In conclusion, the proposed DualQTE is an efficient variational algorithm for quantum time evolution that does not suffer from the quadratic complexity of evaluating the QGT. This cost reduction enables scaling imaginary-time evolution to larger, practically relevant system sizes and allows the simulation and demonstration of a wide variety of important tasks such as Gibbs state preparation, mixed time evolution, or the evaluation of thermodynamic observables. Improving the resource requirements for near-term algorithms is an important step for scaling demonstrations to the full size of today's quantum computers and work towards practical applications.

VI. ACKNOWLEDGEMENTS

is chosen to be the Fubini-Study metric, which, as shown in Ref. [33], if locally approximated, yields the QGT:

d 2 (ϕ(θ), ϕ(θ + δθ)) = arccos 2 | ⟨ϕ(θ)|ϕ(θ + δθ)⟩ | = 1 -| ⟨ϕ(θ)|ϕ(θ + δθ)⟩ | 2 + O(||δθ|| 4 2 ) = ⟨δθ, g(θ)δθ)⟩ + O(||δθ|| 3 2 ).
The formulation in Eq. (B1) is then obtained by solving the minimization problem.

To circumvent the explicit evaluation of the QGT the natural gradient update can instead be calculated without the quadratic local approximation, and instead solve the optimization problem directly. If we use the infidelity as distance metric and replace the loss function gradient by the evolution gradient ∇ℓ(θ) = ∇E(θ)/2 = -b(θ), we obtain the same update rule as the main text

θ (t+1) = argmin θ -⟨b(θ (t) ), θ -θ (t) ⟩ + 1 -F (θ, θ (t) ) 2∆ t . Appendix C: Illustrative example
For an intuitive understanding of the approximations of the QGT norm, we investigate an illustrative example with the variational model |ϕ(θ)⟩ = R Z (θ)R Y (θ) |0⟩, the Hamiltonian H = Z and a timestep of δτ = 1/2. In Fig. 7(a) we compare the exact values of the loss function L for imaginary-time evolution around θ = π/4 obtained by using the metric ⟨δθ, g(θ)δθ⟩ or the infidelity 1 -F (θ, θ + δθ) as norm. At δθ = 0 the approximation is exact and in the vicinity the difference scales as (δθ) 3 , see also Fig. 7(b). Note, that the infidelity is periodic and bounded in [0, 1] but the linear term b T δθ is unbounded, which leads to the fact that the minimum of the infidelity-based loss function close to δθ = 0 is not the global minimum. This is well visible in Fig. 7(a) where the infidelity-based loss function achieves lower values for large δθ than the minimum of the QGT close to δθ = 0. Since we aim to find the same minimum as the QGT-based loss function using a local optimization routine, such as gradient descent, is crucial for the dual time evolution.

Impact of the time perturbation

The approximation error scales with the norm of δθ and, therefore, solving for the update step θ = δθ/δτ with a smaller time perturbation δτ should result in a smaller error in the update step. Remembering the definition of the loss function

L(δθ) = 1 -F (θ, θ + δθ) 2 -δτ • b T (θ)δθ,
we see that a smaller δτ moves the minimum closer to the minimum of the infidelity at δθ = 0, leading to a smaller approximation error. Since the fidelity is bounded but the linear part b T (θ)δθ is not, there is a maximum feasible range for the value of δτ . A necessary condition for the existence of the minimum is that the gradient of the loss function vanishes, ∇L = 0, which requires ∀i ∈ {1, . . . , d} :

1 2 ∂ ∂(δθ) i F (θ, θ + δθ) = -δτ • b i (θ). (C1)
For a circuit with unique parameters and only Pauli rotations gates, the gradient of the fidelity can be bounded via the parameter-shift rule to be in [-1/2, 1/2] (see also Appendix D). Thus, a necessary condition for the timestep perturbation is ∀i ∈ {1, . . . , d}

: δτ ∈ -1 4|b i (θ)| , 1 4|b i (θ)| , (C2)
which can be generalized to circuits with repeated parameters or other than Pauli gates. Note that this is only a necessary and not a sufficient condition for the existence of a minimum since, depending on the circuit structure, the fidelity gradient may not support the full range [-1/2, 1/2].

In Fig. 8(a) we visualize the impact of δτ on the loss landscape. For small time perturbations the QGT-based and dual loss landscapes almost coincide, but if δτ is chosen too large the dual loss function has no minimum. If the loss function can be evaluated exactly, choosing δτ as small as possible therefore minimizes the approximation error. In Fig. 8(b), we find the the error in the parameter derivative θ = δθ/δτ scales approximately as O(δτ ). In practice, however, the loss function is subject to measurement noise and errors in the solution δθ are amplified by 1/δτ . Hence, for a finite number of measurements there is a trade-off between QGT approximation error and controlling the noise amplification. Assuming a forward Euler integration, the Bures distance can be formulated in terms of the QGT as

ε S = 1 T T 0 2 1 -1 -∆ 2 t ∆ θT g(θ)∆ θ dt = 1 T T 0 ∆ t ∆ θg(θ)∆ θdt ≤ 1 T T 0 ∆ t ∥g(θ)∥ 2 ∥∆ θ∥ 2 dt ≤ ∆ t λ max ∥∆ θmax ∥ 2 , (D2)
where we introduced ∆ θ = θθ, λ max ≥ 0 is a bound on the largest eigenvalue of g for any parameter value θ and, similarly, ∥∆ θmax ∥ 2 an upper bound on the norm ∆ θ. In the first line, we dropped O(∆ 3 t ) error terms, in the second line we used a first order Taylor-expansion and in the last line we use the definition of the operator norm to bound the inner product of δθ in the metric of g.

VarQTE

In each VarQTE step we solve a linear system for the update step, where the measurements of the QGT and evolution gradient are subject to sampling error. We define the noisy quantities as g(θ) = g(θ)+∆g(θ) and b(θ) = b(θ)+∆b(θ) and, then, solve the noisy linear system

g(θ) θ = b(θ), (D3)
with the noisy update θ = θ + ∆ θ. To stabilize the linear system and ensure the QGT and it's estimate are invertible, we assume a regularization of g and g in form of a diagonal shift δ c . This shift is a trade-off of stability and bias, which is also discussed in Ref. [40], Appendix D. We can write the error in the update step using the difference of the noisy and exact linear system solutions, as

∥∆ θ∥ 2 = ∥(g + ∆g) -1 (b + ∆b) -g -1 b∥ 2 ≈ ∥(g -1 -g -1 ∆gg -1 )(b + ∆b) -g -1 b∥ 2 = ∥g -1 ∆b -g -1 ∆gg -1 b -g -1 ∆g -1 g -1 ∆b∥ 2 ≈ ∥g -1 ∆b -g -1 ∆g θ∥ 2 ≤ ∥g -1 ∥ 2 ∥∆b∥ 2 + ∥∆g∥ 2 ∥ θ∥ 2 , ( D4
)
where we dropped the explicit parameter dependence for legibility. In the second line we used the Neumann series to approximate (g + ∆g) -1 = g -1g -1 ∆gg -1 + O(∥∆g∥ 2 2 ∥g -1 ∥ 3 2 ) and dropped quadratic error terms on the fourth line. In the following we derive upper bounds on the maximal value of the individual contributions in the error bound, such that we finally obtain a bound ∥∆ θmax ∥ on the error in the update step.

a. Spectrum of g Each QGT entry can be computed as [40] g

ij (θ) = - 1 2 ∂ i ∂ j F (θ ′ , θ) θ ′ =θ . ( D5
)
For a circuit with unique, non-interacting parameter and only plain Pauli rotation gates R P (θ), we can use the parameter-shift rule [34] to write the entry as

g ij (θ) = - 1 2
F (++) ij -F (+-) ij -F (-+) ij + F (--) ij 4 , (D6)
where we abbreviated F (±±) ij = F (θ, θ ± e i π/2 ± e j π/2) for the ith unit vector e i (and jth unit vector e j ). Since the fidelity is in [0, 1] we can bound each entry by

- 1 4 ≤ g ij (θ) ≤ 1 4 , (D7)
for any value of θ. Gershgorin's circle theorem tells us that the maximal eigenvalue of g is bounded from above by the maximal sum over the columns or rows, which in this case is achieved by setting all elements of a column/row to 1/4. This gives the bound

λ max ≤ d i=1 1 4 = d 4 . (D8)
This bound can be generalized to circuits with coefficients or repeated parameters by applying the chain and product rules. For example, for a coefficient-free circuit where parameters can be repeated up to m times, the bound becomes md/4.

b. Norm of the update step The update step can be bounded as ∥ θ∥ 2 ≤ ∥g -1 ∥ 2 ∥b∥ 2 , where ∥g -1 ∥ 2 ≤ δ -1 c . The evolution gradient can be bounded using the parameter-shift rule, under the circuit structure assumptions as the previous section. Each element in the gradient is bounded by

|b i | = |E (+) i -E (-) i | 2 ≤ |E (+) i | + |E (-) i | 2 ≤ E max , (D9)
where

E (±) i = E(θ ± e i π/2
) and E max is the absolute maximum system energy. The norm over all elements is then ∥b∥ 2 ≤ √ dE max , leading to an overall bound of

∥ θ∥ 2 ≤ √ dE max δ c , (D10)
for any parameter value θ. c. Sampling errors Since the measurement noise is unbiased, the random variable ∆g = gg has zero mean with i.i.d. entries. This allows to apply Latala's theorem [51], which states that

E[∥∆g∥ 2 ] ≤ C   max i d j=1 E[(∆g) 2 ij ] + max j d i=1 E[(∆g) 2 ij ] + 4 d i,j=1 E[(∆g) 4 ij ]   , (D11)
for some constant C ∈ R.

Ref. [60] is concerned with the similar case of sampling the matrix [F (x i , x j )] d i,j=1 for a set of parameters {x i } d i=1 . There, the matrix entries are Bernoulli distributed with probability F (x i , x j ). Using QGT representation as Hessian and applying the parameter-shift rule, we can see that the entries of g ij are Poisson binomial distributed [61] with probabilities [F

(++) ij , 1 -F (+-) ij , 1 -F (-+) ij , F (--) ij ] over a shifted support [0, 1, 2, 3, 4] → [-2, -1, 0, 1, 2]
. Since the of this distribution are independent of the number of circuit parameters, it can be shown analogous to Ref. [60] that

E[|(∆g) ij | 2 ] = O 1 N and E[|(∆g) ij | 4 ] = O 1 N 2 , (D12)
which leads to a total bound of

E[∥∆g∥ 2 ] = O d N . ( D13
)
The bound on ∥∆b∥ 2 does not need to be tighter than ∥∆g∥ 2 ∥ θ∥ 2 , which is straightforward to achieve via the sampling error. Using the product rule we have

|∆b i | = | bi -b i | = | Ẽ(+) i - Ẽ(-) i -E (+) i + E (-) i | 2 ≤ | Ẽ(+) i -E (+) i | + | Ẽ(-) i -E (-) i | 2 = O   Var(E (+) i ) + Var(E (-) i ) 2 √ N   . (D14)
The variance of any state |ψ⟩ can be upper bounded by

Var(E) = ⟨ψ|H 2 |ψ⟩ -E 2 ≤ ⟨ψ|H 2 |ψ⟩ ≤ E 2 max . (D15)
Summing over all gradient elements we obtain

∥∆b max ∥ 2 = O √ dE max √ N . (D16)
d. Final bound Plugging the bounds in the previous paragraphs into Eq. (D4) and then into Eq. (D2), we obtain the final bound of

ε S ≤ O d 3/2 E max ∆ t δ 2 c √ N . (D17)
The same asymptotic bound can be derived by performing a moment expansion on the expectation E[ θθ].

As an example we investigate a simple product-state model, which allows to show the tightness of several of the above bounds. We look a the first timestep of the n-qubit Hamiltonian H = n i=1 Z i with an ansatz that consists of a single layer of Pauli-Y rotations, each with an individual parameter. The initial state is |+⟩ ⊗n which is prepared by setting each of the parameters to π/2. Each expectation value is computed with N = 1000 measurements and we use a regularization of δ c = 10 -2 . We then vary the number of qubits from n = 2 to 10 and measure the error term contributions over 10 averages, since ∆g and ∆b are random variables.

The QGT is measure of the correlation between the parameter derivatives and as there is not light-cone connecting any two parameterized gates in the product state ansatz, the QGT is a diagonal matrix. Its norm is therefore ∥g∥ 2 = 1/4 for any system size. With this restriction, we observe in Fig. 9(a) that the bound on the Bures metric in Eq. (D2) is tight as ε S ∝ ∥∆ θ∥ 2 . While all bounds are obeyed, we observe that in particular the bound on ∥ θ∥ 2 is loose, since the bound in Eq. (D10) scales with √ dE max ∝ d 1.5 , but we only observe a d 0.5 scaling. This bound could potentially be further improved by taking into account that the magnitude of the update step is bounded by the change induced of the evolution operator exp(-∆ t H), which is independent of the number of parameters d.

DualQTE

Assume we require K steps to converge. Then the error in the update step δθ is ∥∆(δθ)∥ 2 = ∥∆(δθ (K) )∥ 2 = ∥ δθ (K-1)η ∇L( δθ

(K-1)
)δθ (K-1) + η∇L(δθ (K-1) )∥ 2 ≤ ∥∆(δθ (K-1) )∥ 2 + η∥ ∇L( δθ

(K-1)
) -∇L(δθ (K-1) )∥ 2 ≤ ∥∆(δθ (K-1) )∥ 2 + η∥∆(∇L) max ∥ 2 ≤ ηK∥∆(∇L) max ∥ 2 , loss function is available these criteria become unreliable as the noise in the evaluation might prevent the termination criterion to be fulfilled even though the algorithm converged.

(D18) |0 R Y R Z • R Y R Z |0 R Y R Z • R Y R Z |0 R Y R Z • R Y R Z |0 R Y R Z • R Y R Z . . . . . . . . . . . . . . . |0 R Y R Z R Y R Z
One possible resolution would be to consider a moving average over a past batch of iterations. However, depending on the level of noise, this could require a large batchsize and therefore many iterations until the termination can be checked. Since the dual time evolution only has to compute small corrections, if small timesteps are performed and the optimization are warmstarted, we only expect a few iterations and a moving average is not a resource-efficient solution. Therefore we use a heuristic where the first optimization uses a large number of steps and the subsequent ones perform a fixed number of few iterations.

To demonstrate the effectiveness of warmstarting and to calibrate the number of required steps for noisy evaluations we investigate the dual time evolution in an ideal setting with exact statevector simulations and no finite-sampling statistics. First, we perform the time evolution for a Heisenberg Hamiltonian with periodic boundary conditions with n = 12 sites, J = 1/4, g = -1 and the initial state |+⟩ ⊗n . As circuit model we use the hardware efficient circuit from Fig. 10 with r = 6 repetitions and optimize the update step with a gradient descent routine with a fixed learning rate of η = 0.1. In each timestep we iterate until the change in loss function ∆L is below the threshold of 10 -4 ∆ t = 10 -6 . The results are presented in Fig. 11(a), and we observe that warmstarting drastically reduces the number of required iterations until the convergence criterion is reached.

In a second experiment we analyze how the required number of optimization steps scales with the system size. In Fig. 11(b) we repeat the above experiment for n = 3 to 12 spins and track the number of steps in the first iteration and the mean and standard deviation of the warmstarted iterations. We see that the number of steps scales sublinearly in the number of parameters d and is almost constant for the warmstarted iterations. k ∥ 2 at selected times t, where k indicates the iteration in the optimization within the timestep. Lengths for the fidelity gradients differ as the optimization were performed using a different number of steps, see Table II. In the real-time evolution of the Heisenberg model we use the circuit sketched in Fig. 13, which is the same model used in Ref. [START_REF] Barison | [END_REF]. The dotted box is repeated three times and the rotation layer alternates between Pauli-X rotations (starting from the first layer) and Pauli-Y rotations.

|0

R

X R ZZ R Y |0 R X R ZZ R Y |0 R X R ZZ R Y |0 R X R Y FIG.
13: The circuit model used for the real-time evolution experiments.

FIG. 1 :

 1 FIG. 1: Estimated runtimes of variational imaginary-time evolution (VarQITE) and our proposed dual method (DualQITE) as a function of the number of parameters d of the variational model, for an exemplary Heisenberg model and 200 timesteps. See Appendix A for more details.

FIG. 2 :

 2 FIG. 2: (a) The mean and standard deviation of DualQITE and VarQITE, each averaged over 5 independent experiments for a varying number of shots. (b) The accuracy measured in integrated Bures distance I B (see Eq. (16)) for DualQITE and VarQITE, with mean and standard deviation of 5 experiments. The resources are measured in total number of measurements and are shown for evaluation of gradients (and the QGT, in the case of VarQITE) using parameter-shift rules (dashed) or a LCU approach (solid lines).

FIG. 3 :

 3 FIG. 3: (a) Total number of measurements required to achieve a mean accuracy of I B ≤ 0.1 over an average of 5 experiments. See Table II for the exact algorithm settings. Dotted and dashed lines show fits for the number of measurements. The bumps in the fits are due to the discontinuity of the circuit depth, which depends on ⌈log 2 (n)⌉. VarQITE is not evaluated for n = 14 qubits as it requires too many measurements. (b) Mean accuracy and standard deviation of each point of the top panel. The grey line indicates the infidelity threshold of I B = 0.1.

2 .

 2 Evolve up to imaginary time t = β/2 |ϕ m (β/2)⟩ ∝ e -βH/2 |ϕ m (0)⟩ .

FIG. 4 :

 4 FIG. 4: Energy per site for the Heisenberg model on a 6-spin chain, comparing mean and standard deviation of QMETTS with DualQITE (blue circle and errorbars) with a reference METTS implementation (black line and grey shade).

FIG. 5 :

 5 FIG. 5: Average magnetization in X and Z direction as tracked by different variational algorithms.

FIG. 6 :

 6 FIG. 6: Error of the real-time evolution in Bures distance, plus error bounds obtained with VarQRTE and DualQRTE.

FIG. 7 :

 7 FIG. 7: (a) Values of the loss function L for evaluation with the QGT metric, and with introduced infidelity approximation. (b) Difference of the QGT metric and infidelity as function of the perturbation δθ.

FIG. 8 :

 8 FIG. 8: (a) Loss landscapes and optimal solutions of the original, QGT-based loss function and the dual loss function for different δτ . (b) Error in calculating the parameter derivative θ depending on δτ and the number of measurements N .

FIG. 9 :

 9 FIG. 9: Scaling of different error bound contributions for a product state setting. The labels include the scaling with number of parameters, i.e. ∝ α indicates a scaling with d α , where d is the number of parameters.

FIG. 10 :

 10 FIG. 10: The hardware efficient ansatz for the imaginary-time evolution experiments. In the QMETTS experiments the Pauli-Y rotations is replaced by Pauli-X rotations, if the evolution starts in the Y basis states |±i⟩.

FIG. 11 : 2 FIG. 12 :

 11212 FIG. 11: (a) The number of iterations required per timestep until convergence is reached with different initialization techniques. (b) The number of iterations for different numbers of qubits and the first iteration and warmstarted iterations. The warmstarted points show mean and standard deviation of the number of iteration of all steps after the first.

  Appendix F: Real-time evolution of the Heisenberg model 1. Circuit diagram

We thank Christa Zoufal, Stefano Barison, Almudena Carrera Vazquez, David Sutter, Caroline Tornow, Laurin Fischer and Daniel Egger for insightful conversations on this project.

We acknowledge the use of IBM Quantum services for this work. The views expressed are those of the authors, and do not reflect the official policy or position of IBM or the IBM Quantum team.

IBM, the IBM logo, and ibm.com are trademarks of International Business Machines Corp., registered in many jurisdictions worldwide. Other product and service names might be trade-marks of IBM or other companies. The current list of IBM trademarks is available at https://www.ibm.com/legal/copytrade.

Appendix A: Runtime estimates of variational time evolution

The benchmark in Sec. III C provides a scaling for the total number of measurements N required by VarQITE and DualQITE, which allows a runtime estimation on the algorithms on quantum hardware. In this estimation we neglect the overhead of classical processors and assume a superconducting quantum computer with a basis gate set including √ X, R Z and CX gates, as reported by several IBM Quantum backends, for example. This gate set allows to compile any sequence of single qubit gates into two √ X gates and three virtual R Z gates. For an n-qubit simulation of the Heisenberg model and the considered circuit model (see Fig. 10) with r repetitions, the time for a single measurement can then be approximated as

where t CX is the duration of a CX gate, t √ X the duration of a √ X gate, t meas the time of a measurement and t reset the time to reset the qubits for the next execution. Since the R Z gates are virtual they do not contribute to the runtime. The total runtime is then estimated by N t shot .

Depending on the architecture and the gate decomposition the duration and fidelity of single-and two-qubit operation, as well as measurements, varies on superconducting qubit chips [55]. Here we use gate times of t CX = 451ns, t √ X = 36ns and t meas = 860ns as reported by ibm peekskill (v2.6.5), which is an IBM Quantum Falcon processors [START_REF]IBM Quantum[END_REF]. For shallow circuits in particular, the time to reset qubits for the following execution is a crucial bottleneck. The reset operation can, for example, be implemented by waiting 5-10 T1 times and let the qubits decay to the computational ground-state, but with T1 times of the order of 400µs the reset via relaxation is orders of magnitude slower than the other circuit operations. Active resets instead measure the qubit state and apply an X-operation conditionally if the state |1⟩ is measured. This technique allows to reduce the reset times to typically 50 to 250µs on IBM hardware [START_REF] Tornow | [END_REF], which, however, still dominates the overall runtime for the considered circuits. By using a second excited state it is possible to implement reset schemes with 500ns to 2µs [58,59] and we therefore use t reset = 2µs in our estimation.

Appendix B: Derivation via quantum natural gradient descent

VarQITE is inherently connected to the quantum natural gradient (QNG) algorithm [33]. In fact, this connection is a motivation for the convergence of the QNG as imaginary-time evolution is guaranteed to converge to the ground state, if there is sufficient initial overlap with it.

With a forward Euler integration the VarQITE update rule is

This coincides with the QNG update step for the loss function ℓ(θ) = E(θ)/2 and a learning rate of η = ∆ t ,

The natural gradients step can be expressed in a dual formulation as

with a distance metric d. In this equation we see that the update step is going into the opposite direction of the gradient ∇ℓ, while the magnitude is limited by the distance metric and the learning rate. Standard gradient descent uses the model-agnostic ℓ 2 norm as distance metric. Natural gradients on the other hand limit the update step by the amount of change it induces in the model. To measure the induced change the metric d In this section, we present the derivation on the upper bound of the sample complexity of VarQTE and DualQTE. The target error is measured in integrated Bures distance,

where we used that the error at the initial point is zero, ∥∆(δθ (0) ∥ 2 = 0. The error in the loss function gradient can then be written as

where ∆(∇F (θ, θ + δθ)) = ∇F (θ, θ + δθ) -∇F (θ, θ + δθ) and ∥∆(∇F ) max ∥ 2 is an upper bound on the maximum fidelity gradient error for any parameter.

The error in the gradient of F can be derived via the parameter-shift rule, as

where we used that the variance of the fidelity can be bounded for any state |ψ⟩ as

Hence the total error of the fidelity gradient in ℓ 2 norm is

The bound on ∥∆b max ∥ 2 is already derived in the previous subsection, which gives then a total of

Using the definition θ = δθ/δτ we then obtain

Appendix E: Imaginary-time evolution of the Heisenberg model

Circuit diagram

The circuit used as variational model is schematically presented in Fig. 10. Each Pauli rotation gate has an independent parameter and the dotted box is repeated several times. For r repetitions and n qubits the total number of tunable parameters is thus 2n(r + 1). The CNOT entangling gates are arranged in a pairwise manner to minimize the total depth to 2 per entangling layer.

Termination and warmstarting

Termination criteria for gradient descent algorithms are typically defined as achieving a minimal threshold in the difference in the loss function between update steps or in the gradient norm. However, if only noisy readout of the

Resource requirements for the dual time evolution

This section shows the detailed VarQITE and DualQITE settings for the resource estimations in Sec. III. For VarQITE we only varied the number of shots and in the dual method we additionally allowed to vary the number of iterations in the optimization in each time step. Especially DualQITE has a lot of additional degrees of freedom that could be optimized, such as the kind of optimizer, in addition to settings shared with VarQITE, such as timestep size.

Table I shows the settings for VarQITE and DualQITE for the resource estimation in Fig. 2(b) and Table II the settings for the scaling with system size in Fig. 3.

Gradient benchmark

In this section, we measure how the norm of the loss function gradient, defined as

scales with the number of qubits n in the Heisenberg Hamiltonian. This Hamiltonian is 2-local and since the ansatz depth grows logarithmic with the number of qubits we do not expect exponentially vanishing gradients for the evolution gradient b [42]. Further, as discussed in the main text, the initial point of each timestep δθ is close to 0, which ensures the circuit required to measure to fidelity gradient is close to the identity and we do not expect to encounter a barren plateau [48]. Since both parts of the loss functions are not in a barren plateau setting, we expect the loss function to be measurable efficiently. In Fig. 12, we measure the ℓ 2 norm of both the evolution gradient and fidelity gradients and we find that neither gradient vanishes exponentially. Instead, the evolution gradient increases with system size, which reflects the extensiveness of energy in the Heisenberg model. Since we perform imaginary time evolution the system converges towards the ground state and we expect the energy gradients to vanish, once converged. In this case, we do not infer a barren plateau as the gradient norm does not systematically decrease faster for larger systems. Similarly, the fidelity gradients are expected to decay as the optimal parameter update δθ is found. Note that for system sizes of 12 qubits the exponential decay of the gradients is typically clearly visible [43].