### The VFF algorithm

#### Overview

Given a Hamiltonian *H* on a *d* = 2^{n} dimensional Hilbert space (i.e., on *n* qubits) evolved for a short time Δ*t* with the simulation unitary *e*^{−iHΔt}, the goal is to find an approximation that allows the simulation at later times *T* to be fast-forwarded beyond the coherence time. Figure 2 schematically shows the VFF algorithm, which consists of the following steps:

- 1.
Implement a unitary circuit

*U*(Δ*t*) to approximate*e*^{−iHΔt}, the simulation at a small timestep. - 2.
Compile

*U*(Δ*t*) to a diagonal factorization*V*=*W**D**W*^{†}≈*e*^{−iHΔt}with circuit depth*L*. - 3.
Approximately fast-forward the QS at large time

*T*=*N*Δ*t*using the same circuit of depth*L*:*e*^{−iHT}≈*W**D*^{N}*W*^{†}.

Typically *U*(Δ*t*) will be a single-timestep Trotterized unitary approximating *e*^{−iHΔt}. We variationally search for an approximate diagonalization of *U*(Δ*t*) by compiling it to a unitary with a structure of the form

$$V({boldsymbol{alpha }},Delta t):=W({boldsymbol{theta }})D({boldsymbol{gamma }},Delta t)W{({boldsymbol{theta }})}^{dagger },$$

(1)

with ** α** = (

**,**

*θ***) being a vector of parameters. Here,**

*γ**D*(

**, Δ**

*γ**t*) is a parameterized unitary composed of commuting unitaries that encode the eigenvalues of

*U*(Δ

*t*), while

*W*(

**) is a parameterized unitary matrix consisting of corresponding eigenvectors. In the “Methods” section, we describe layered structures that provide ansätze for the circuits**

*θ**W*(

**) and**

*θ**D*(

**, Δ**

*γ**t*), and we detail our gradient-descent optimization methods for training

**and**

*θ***.**

*γ*To approximately diagonalize *U*(Δ*t*), the parameters ** α** = (

**,**

*θ***) are variationally optimized to minimize a cost function**

*γ**C*

_{LHST}(

*U*(Δ

*t*),

*V*) that can be evaluated using a short-depth quantum circuit called the Local Hilbert–Schmidt Test (LHST)

^{35}shown in Fig. 2c. The compilation procedure we employ to approximate

*U*(Δ

*t*) by

*V*(

**, Δ**

*α**t*) makes use of the quantum-assisted quantum compiling (QAQC) algorithm

^{35}, that was later shown to be robust to quantum hardware noise

^{36}. The “Cost function and cost evaluation” section below elaborates on our cost function.

If we can find such an approximate diagonalization for *U*(Δ*t*) then, for any total simulation time, *T* = *N*Δ*t*, we have:

$${e}^{-iHT}={({e}^{-iHDelta t})}^{N},$$

(2)

$$approx {(U(Delta t))}^{N},$$

(3)

$$approx W({boldsymbol{theta }})D{({boldsymbol{gamma }},Delta t)}^{N}W{({boldsymbol{theta }})}^{dagger },$$

(4)

$$=W({boldsymbol{theta }})D(N{boldsymbol{gamma }},Delta t)W{({boldsymbol{theta }})}^{dagger }.$$

(5)

Hence, a QS for any total time, *T*, may be performed with a fixed-quantum circuit structure as depicted in Fig. 2d. In “Simulation error analysis”, we perform an error analysis to investigate how the approximate equalities in Eqs. (3) and (4) affect the overall simulation error.

#### Cost function and cost evaluation

For the variational compiling step of VFF (shown in Fig. 2c), we employ the cost function *C*_{LHST}(*U*, *V*) introduced in ref. ^{35}. This is defined as

$${C}_{{rm{LHST}}}(U,V)=1-frac{1}{n}mathop{sum }limits_{j = 1}^{n}{F}_{e}^{(j)},$$

(6)

where the ({F}_{e}^{(j)}) are entanglement fidelities and hence satisfy (0le {F}_{e}^{(j)}le 1). Specifically, ({F}_{e}^{(j)}) is the entanglement fidelity for the quantum channel obtained from feeding into the unitary *U**V*^{†}, the maximally mixed state on (overline{j}) and then tracing over (overline{j}) at the output of *U**V*^{†}, where (overline{j}) contains all qubits except for the *j*-qubit. We elaborate on the form of *C*_{LHST}(*U*, *V*) in Supplementary Note 1.

This function has several important properties.

- 1.
It is faithful, vanishing if and only if

*V*=*U*(up to a global phase). - 2.
Nonzero values are operationally meaningful. Namely,

*C*_{LHST}(*U*,*V*) upper bounds the average-case compilation error as follows:$${C}_{{rm{LHST}}}(U,V)ge frac{d+1}{nd}(1-overline{F}(U,V)),$$

(7)

where (overline{F}(U,V)) is the average fidelity of states acted upon by

*V*versus those acted upon by*U*, with the average being over all Haar-measure pure states. - 3.
The cost function appears to be trainable, in the sense that it does not have an obvious barren plateau issue (i.e., exponentially vanishing gradient, see refs.

^{35,37}). - 4.
Estimating the cost function is DQC-1 hard and hence it cannot be efficiently estimated with a classical algorithm

^{35}. - 5.
There exists a short-depth quantum circuit for efficiently estimating the cost and its gradient.

Regarding the last point, each ({F}_{e}^{(j)}) term in Eq. (6) is estimated with a different quantum circuit, and then one classically sums them up to compute *C*_{LHST}(*U*, *V*). An example of such a circuit is depicted in Fig. 2c. It involves 2*n* qubits, with the top (bottom) *n* qubits denoted *A* (*B*). The probability of the 00 measurement outcome on qubits *A*_{j}*B*_{j} in this circuit is precisely the entanglement fidelity ({F}_{e}^{(j)}). Therefore 2*n* single-qubit measurements are required to compute *C*_{LHST}(*U*, *V*), a favorable scaling compared to, for example, the *O*(*n*^{4}) measurements that are naively required for VQE^{15,38,39,40}. We also remark that *C*_{LHST}(*U*, *V*) was recently shown to have noise resilience properties, in that noise acting during the circuit in Fig. 2c tends not to affect the global optimum of this function^{36}.

For simplicity, we will often write our cost function as

$${C}_{{rm{LHST}}}^{{rm{VFF}}}(T):={C}_{{rm{LHST}}}({U}^{frac{T}{Delta t}},{V}^{frac{T}{Delta t}}),$$

(8)

with *U* = *U*(Δ*t*) and *V* = *V*(** α**, Δ

*t*), and note that ({C}_{{rm{LHST}}}^{{rm{VFF}}}(Delta t)) is the quantity that we directly minimize in the optimization loop of VFF.

#### Simulation error analysis

### Linear scaling in *N*

In practice, each of the steps in the VFF algorithm above will generate errors. This includes the algorithmic error from the approximate implementation, *U*(Δ*t*), of the infinitesimal time evolution operator *e*^{−iHΔt} and error from the approximate compilation and diagonalization of *U*(Δ*t*) into *V*(** α**, Δ

*t*). These two error sources bound the overall error via the triangle inequality:

$${epsilon }_{p}^{{rm{FF}}}(Delta t)le {epsilon }_{p}^{{rm{TS}}}(Delta t)+{epsilon }_{p}^{{rm{ML}}}(Delta t).$$

(9)

Here, ({epsilon }_{p}^{{rm{FF}}}(Delta t)) is the overall simulation error for time Δ*t*, ({epsilon }_{p}^{{rm{TS}}}(Delta t)) is the Trotterization error (note that this error may always be reduced using higher-order Trotterizations at the cost of more gates), and ({epsilon }_{p}^{{rm{ML}}}(Delta t)) is the “machine learning” error associated with the variational compilation step. These quantities are defined as

$${epsilon }_{p}^{{rm{FF}}}(Delta t)=| | {e}^{-iHDelta t}-V({boldsymbol{alpha }},Delta t)| {| }_{p},$$

(10)

$${epsilon }_{p}^{{rm{TS}}}(Delta t)=| | {e}^{-iHDelta t}-U(Delta t)| {| }_{p},$$

(11)

$${epsilon }_{p}^{{rm{MS}}}(Delta t)=| | U(Delta t)-V({boldsymbol{alpha }},Delta t)| {| }_{p},$$

(12)

where (| | M| {| }_{p}={({sum }_{j}{m}_{j}^{p})}^{1/p}) is the Schatten *p*-norm, with {*m*_{j}} the singular values of *M*.

Ultimately, we are interested in fast-forwarding and hence we want to bound ({epsilon }_{p}^{{rm{FF}}}(T)) with *T* = *N*Δ*t*. For this purpose, we prove a lemma in Supplementary Note 2 stating that

$$| | {U}_{1}^{N}-{U}_{2}^{N}| {| }_{p}le N| | {U}_{1}-{U}_{2}| {| }_{p},$$

(13)

for any two unitaries *U*_{1} and *U*_{2}. Combining this lemma with the triangle inequality in Eq. (9) gives

$${epsilon }_{p}^{{rm{FF}}}(T)le N({epsilon }_{p}^{{rm{TS}}}(Delta t)+{epsilon }_{p}^{{rm{ML}}}(Delta t)).$$

(14)

Equation (14) implies that the overall simulation error scales at worst linearly with the number of timesteps, *N*.

We remark that, for the special case of *p* = 2, Eq. (13) can be reformulated in terms of our cost function as:

$${C}_{{rm{LHST}}}^{{rm{VFF}}}(T)mathop{<}limits_{approx} n {N}^{2} {C}_{{rm{LHST}}}^{{rm{VFF}}}(Delta t),$$

(15)

with ({C}_{{rm{LHST}}}^{{rm{VFF}}}) given by Eq. (8). The approximation in Eq. (15) holds when the cost function ({C}_{{rm{LHST}}}^{{rm{VFF}}}(Delta t)) is small, which is the case after a successful optimization procedure. See Supplementary Note 2 for the non-approximate version of Eq. (15). Thus, we find that the VFF cost function scales at worst quadratically in *N* under fast-forwarding.

### Certifiable error and a termination condition

Equation (14) holds for all Schatten norms, but of particular interest for our purposes is the Hilbert–Schmidt norm, *p* = 2, from which we can derive certifiable error bounds on the average-case error. In addition, the operator norm, *p* = *∞*, quantifies the worst-case error and is often used in the QS literature^{41,42}. For our numerical implementations (“Implementations” section), we will consider both worst-case and average-case error. On the other hand, for our analytical results presented here, we will focus on average-case error since it is naturally suited to providing a termination condition for the optimization in VFF.

As an operationally meaningful measure of average-case error, we consider the average gate fidelity between the target unitary *e*^{−iHT} and the approximate simulation *V*(** α**,

*T*), arising from the VFF algorithm:

$$overline{F}(T)={int}_{psi }| leftlangle psi right|V{({boldsymbol{alpha }},T)}^{dagger }{e}^{-iHT}left|psi rightrangle {| }^{2}dpsi ,$$

(16)

where the integral is over all states (left|psi rightrangle) chosen according to the Haar measure.

In Supplementary Note 2, we show that one can lower bound (overline{F}(T)) based on the value of the VFF cost function,

$$overline{F}(T)mathop{>}limits_{approx} 1-frac{d}{d+1}{N}^{2}{left({epsilon }_{infty }^{{rm{TS}}}(Delta t)+sqrt{n{C}_{{rm{LHST}}}^{{rm{VFF}}}(Delta t)}right)}^{2}.$$

(17)

This inequality holds to a good approximation in the limit that ({C}_{{rm{LHST}}}^{{rm{VFF}}}(Delta t)) is small, as is the case after a successful optimization procedure. See Supplementary Note 2 for the exact lower bound on (overline{F}(T)), from which Eq. (17) is derived.

In addition, Eq. (17) provides a termination condition for the variational portion of VFF. If one has a desired threshold for (overline{F}(T)), then this threshold can be guaranteed provided that ({C}_{{rm{LHST}}}^{{rm{VFF}}}(Delta t)) is below a certain value. Once ({C}_{{rm{LHST}}}^{{rm{VFF}}}(Delta t)) dips below this value, then the variational portion of VFF can be terminated. Specifically, the termination condition is ({C}_{{rm{LHST}}}^{{rm{VFF}}}(Delta t)le {C}_{text{Threshold}}), where

$${C}_{text{Threshold}}approx frac{1}{n}{left(frac{1}{N}sqrt{frac{d+1}{d}(1-overline{F}(T))}-{epsilon }_{infty }^{{rm{TS}}}(Delta t)right)}^{2},$$

(18)

with the approximation holding when ({C}_{{rm{LHST}}}^{{rm{VFF}}}(Delta t)) is small. Again, for the exact expression for *C*_{Threshold}, see Supplementary Note 2.

### Implementations

Here, we present results simulated classically and on quantum hardware. We refer the reader to the “Methods” section for details about our ansatz and optimization methods.

For our results below, it will be convenient, for a given simulation error tolerance *δ*, to define fast-forwarding as the case whenever ({R}_{delta }^{{rm{FF}}}>1), where

$${R}_{delta }^{{rm{FF}}}={T}_{delta }^{{rm{FF}}}/{T}_{delta }^{{rm{Trot}}}$$

(19)

is the ratio of the simulation time achievable with VFF, ({T}_{delta }^{{rm{FF}}}), to the simulation time achievable with standard Trotterization, ({T}_{delta }^{{rm{Trot}}}). We note that ({T}_{delta }^{{rm{Trot}}}) is a good empirical measure of the coherence time, since it can account for both decoherence and gate infidelity, and hence the condition ({R}_{delta }^{{rm{FF}}}>1) intuitively captures the idea of simulation beyond the coherence time.

#### Comparing VFF to Trotterization and compiled Trotterizations

Here, we compare the performance of VFF to that of two other simulation strategies. One of these strategies is the standard Trotterization approach depicted in Fig. 1a. Another strategy is to first optimally compile the Trotterization step to a short-depth gate sequence, and use this optimal circuit (in place of the Trotterization step) for the approach in Fig. 1a. We refer to this as the QAQC strategy, since one can use the QAQC algorithm^{35} to obtain the optimal compilation of the Trotterization step. With the QAQC strategy, when finding the optimal short-depth compilation, we make no assumptions about the structure of the compiled circuit, which is in contrast with Trotterization, where the structure of the circuit is dictated by interactions in the Hamiltonian.

For concreteness, we consider the task of simulating the XY model, defined by the Hamiltonian

$${H}_{{rm{XY}}}=-mathop{sum }limits_{i}{X}^{i}{X}^{i+1}+{Y}^{i}{Y}^{i+1},$$

(20)

where *X*^{i} and *Y*^{i} are Pauli operators on qubit *i*. We consider a five-qubit system with open boundary conditions. From the analytical diagonalization of the XY model^{43}, it follows that the ansatz for the diagonal matrix *D* can be truncated at the first nontrivial term, as described in the “Methods” section (see Eq. (26)).

Figure 3 summarizes our results. Figure 3a shows how the ({C}_{{rm{LHST}}}^{{rm{VFF}}}) cost function is iteratively minimized during the optimization procedure. We selected four approximate diagonalizations corresponding to different cost values denoted by VFF_{n}, *n* = 1, …, 4, depicted by colored circles in Fig. 3a. Note that the colors match those used in other panels. Figure 3b compares these diagonalizations with different QS methods: Trotterization (dashed line), QAQC-compiled (dashed-dotted line), and VFF at different stages of optimization. We compare simulated time evolution governed by the XY model and observe how the fidelity decays with evolution time. The fidelity is given by ({rm{Tr}}(left|psi (T)rightrangle leftlangle psi (T)right|rho (T))), where (left|psi (T)rightrangle) is the exact evolved state and *ρ*(*T*) is the state obtained with a noisy simulator. The initial state (left|psi (0)rightrangle) was chosen randomly such that it has nonvanishing overlap with every eigenstate of the Hamiltonian. The circuits were simulated using a noisy trapped-ion simulator with error model from ref. ^{44}. We took error rates, as specified in their Fig. 14 in ref. ^{44} and reduced them by a factor of five for clearer demonstration of VFF’s capabilities.

Results shown in Fig. 3b indicate that QAQC performed better than Trotterization, which is expected as QAQC optimizes in a circuit space larger than that given by Trotterization. Both approaches give better results than VFF_{1} (red curve). This confirms the intuition that at early stages of the VFF optimization (large values of ({C}_{{rm{LHST}}}^{{rm{VFF}}})), the error of the diagonalization is too big to outperform other methods. As the cost function is decreased, the length of time one can simulate using VFF increases. Indeed, as one can see from the green, blue, and purple curves (which are associated with cost values ⪅10^{−2}), VFF dramatically outperforms both Trotterization and QAQC.

One can see another important feature in Fig. 3b. For short simulations, Trotterization and QAQC are always more accurate than VFF, no matter how accurate the diagonalization is. That is because for small *T*, there are just a few timesteps taken by Trotterization and QAQC implementations, and the resulting circuits are shorter than VFF circuits implementing *W*, *D*, and *W*^{†}. This disadvantage rapidly diminishes since VFF circuits do not grow with *T* and the only error that impacts the fidelity comes from imperfect diagonalization. On the other hand, Trotter and QAQC circuits grow linearly with *T* and as a result, fidelity is dominated by noise (and not imperfections in the decomposition).

Figure 3c shows how the fast-forwarding factor ({R}_{delta }^{{rm{FF}}}) and the error in the eigenvalue approximation (pertinent if VFF is used for eigenvalue estimation as discussed in the “Estimating energy eigenvalue“ section) depend on the cost function ({C}_{{rm{LHST}}}^{{rm{VFF}}}). The data suggest power-law dependence in both cases. Bringing the cost function down to 10^{−3} allows us to reduce the eigenvalue error <0.1 and reach a fast-forwarding factor of ~30. Note that for VFF to be more efficient than Trotterization (({R}_{delta }^{{rm{FF}}},>,1)), one has to lower the cost function below ~0.04. For this case, *δ* is defined as (1-{rm{Tr}}(left|psi (T)rightrangle leftlangle psi (T)right|rho (T))) and we considered *δ* = 0.2. Figure 3d presents a more detailed analysis of the eigenvalue error, showing how the error of individual eigenvalues is reduced as the cost function is minimized.

#### Using VFF to fast-forward models across a range of parameters

Here, we show how to use VFF to efficiently find diagonalizations for new models that are nearby in parameter space, from previously diagonalized models.

### Hubbard model

We applied VFF to Trotterized QS unitaries, (U(Delta t)approx {e}^{-i{H}_{{rm{Hub}}}Delta t}), of the Fermi–Hubbard model

$${H}_{{rm{Hub}}}=-tau mathop{sum }limits_{i,j,sigma }left({c}_{i,sigma }^{dagger }{c}_{j,sigma }+{c}_{j,sigma }^{dagger }{c}_{i,sigma }right)+umathop{sum }limits_{i = 1}^{N}{n}_{i,uparrow }{n}_{i,downarrow }.$$

(21)

Here, ({c}_{i,sigma }^{dagger }) and *c*_{i,σ} are electron creation and annihilation operators (resp.) for spin *σ* ∈ {*↓*, *↑*} at site *i* and ({n}_{i,sigma }={c}_{i,sigma }^{dagger }{c}_{i,sigma }) is the electron number operator. The parameters *τ* and *u* are the hopping strength and on-site interaction (resp.). We studied a two-site, two-qubit Fermi–Hubbard model^{45}, which, after translation via the Jordan-Wigner transform, takes the form

$${H}_{{rm{Hub}},2}=-tau (Xotimes I+Iotimes X)+uZotimes Z.$$

(22)

We took *τ* = 1 for our initial diagonalization, then perturbatively increased *u* from 0 to 0.1 in increments of 0.01. For *U*(Δ*t*), we used a first-order Trotterization of (exp (-i{H}_{{rm{Hub}},2}Delta t)). We set a threshold for optimization of 10^{−6}. We used a three-layer ansatz for *W* and a two-layer ansatz for *D*, which we describe in the “Ansatz” section.

In representative results shown in Fig. 4a, b, we see that, after an initial optimization taking a number of iteration steps, VFF reached the optimization threshold. Then, as we perturbed away from *u* = 0, VFF rapidly found new parameters that diagonalized (exp (-i{H}_{{rm{Hub}},2}Delta t)) to below the cost threshold. For all approximate diagonalizations, for an error tolerance of *δ* = 10^{−2}, the simulation error remained below this tolerance for *T* = 30Δ*t*. The diagonalization used nine single-qubit gates and seven two-qubit gates. The Trotterization used two single-qubit gates and one two-qubit gate. Thus, the fast-forwarded simulations used nine single-qubit layers and seven two-qubit gates, but the equivalent Trotterized simulations used 60 single-qubit gates and 30 two-qubit gates. Thus, VFF gave significant depth compression versus the Trotterized simulations, particularly with respect to entangling gates.

### Heisenberg model

Next, we applied VFF to the Heisenberg model,

$${H}_{{rm{Heis}}}=mathop{sum }limits_{i}{J}_{z}{Z}^{i}{Z}^{i+1}+{J}_{x}{X}^{i}{X}^{i+1}+{J}_{y}{Y}^{i}{Y}^{i+1}+h{Z}^{i},$$

(23)

where *X*^{j}, *Y*^{j}, and *Z*^{j} are Pauli spin matrices acting on qubit *j*, and *h*, *J*_{x}, *J*_{y}, and *J*_{z} are parameters.

Here, we took *h* = 1.0 and investigated the model acting on three qubits (whose Hamiltonian we denote *H*_{Heis,3}). We used a first-order Trotterization of (exp (-i{H}_{{rm{Heis,3}}}Delta t)). We set an optimization threshold of 10^{−6} and used a ten-layer ansatz for *W* and a two-layer ansatz for *D*. From *J*_{z} = 1.0 (a noninteracting Hamiltonian), we increased *J*_{z} to 5.0 in increments of 1.0. For these parameter values, *H*_{Heis} is an antiferromagnetic classical Ising model.

Next, we kept *h* = 1.0 and *J*_{z} = 5.0 fixed, and increased *J*_{x} = *J*_{y} from 0.0 to 8.0 in increments of 2.0. When *J*_{x} = *J*_{y}, these are often called XXZ Heisenberg models.

Finally, we kept *h* = 1.0, *J*_{z} = 5.0, *J*_{x} = 8.0 and varied *J*_{y} from 0.0 to 10.0 in increments of 1.0 (XYZ Heisenberg models).

As may be seen in the representative results plotted in Fig. 4c, d, VFF rapidly found new diagonalizations (WD{W}^{dagger }approx exp (-i{H}_{{rm{Heis,3}}}Delta t)) for all models considered. We performed additional searches for diagonalizations of ferromagnetic models (*J*_{z}, *J*_{x}, and *J*_{y} < 0) with similar results. For all approximate diagonalizations, the simulation error remained below an error tolerance of *δ* = 10^{−2}, up to *T* ≈ 100Δ*t*. For this simulation time, each diagonalization used 40 two-qubit gates and 71 single-qubit gates (111 total), whereas each Trotterization used 1200 two-qubit gates and 2500 single-qubit gates (3700 total).

#### VFF implemented on quantum hardware

We implemented VFF on 1 + 1 qubits (i.e., diagonalizing a random single-qubit unitary) on the Rigetti Aspen-4 QC (Fig. 5). Here, we considered the first-order Trotterization of the Hamiltonian *H* = *α*_{x}*σ*_{x} + *α*_{y}*σ*_{y} + *α*_{z}*σ*_{z}, where *α* was a randomly chosen unit vector, at the time Δ*t* = 0.5. We used *W* = *R*_{z}(*θ*_{z})*R*_{x}(*θ*_{x}) and *D* = *R*_{z}(*γ*_{z}). The VFF cost function, as evaluated on the QC with *n*_{samp} = 10^{4}, was optimized to ({C}_{{rm{LHST}}}^{{rm{VFF}}}(Delta t)approx 1{0}^{-1}).

With this system, we investigated how well VFF performed by classically computing the true, noiseless, cost for the parameters found on the Rigetti QC. This true cost converged to two orders of magnitude below the QC-evaluated cost, demonstrating significant robustness of VFF to the noise on the Rigetti QC.

We next simulated single-qubit evolution on the QC by (1) iterating the original Trotterization, *U*(Δ*t*)^{N}, and (2) using the VFF diagonalization (5). We then used process tomography to compare the resultant noisy process resulting from the Trotterization, and the process resulting from VFF to the exact process *U*(Δ*t*)^{N} calculated classically.

In this single-qubit case, the Trotterized simulation unitary could have been compiled to a circuit with many fewer gates; however, this would not be true for higher-dimensional unitaries and for this reason we did not compile the iterated gate sequence here.

In Fig. 5b, we show that VFF performed much better than the iterated Trotterization, giving a high-fidelity simulation. In these results, the entanglement fidelity between the process implemented using VFF and the exact process remained high until at least *N*_{VFF} = 150 and never reached a value <0.7. On the other hand, the fidelity of the iterated Trotterization approach was already 0.586 by *N* = 25.

These results demonstrate that VFF on current QCs can allow for simulations beyond the coherence time. For example, taking an entanglement infidelity of 0.3 as our error tolerance *δ*, it follows from the table in Fig. 5b that we obtained a fast-forwarding beyond the coherence time of at least ({R}_{delta }^{{rm{FF}}}=6).

#### Estimating energy eigenvalues

We primarily foresee VFF being used to study the long-time evolution of the observables of a system. But one may also use VFF to reduce the gate complexity of eigenvalue estimation algorithms, such as quantum phase estimation^{46} or time series analyses^{47,48}. Such algorithms require simulating a Hamiltonian up to time (T={mathcal{O}}left(frac{1}{sigma }right)) to obtain eigenvalue estimates of accuracy *σ*. Due to the constant depth circuits produced by VFF, we can therefore reduce the number of gates required for these algorithms by a factor of ({mathcal{O}}left(frac{1}{sigma }right)), increasing the viability of eigenvalue estimation on noisy QCs.

The eigenvalues of the Hamiltonian simulated via VFF are not directly accessible from the diagonal unitary *D*, since they are encoded in a set of Pauli operators. However, these can be extracted using the time series analysis in ref. ^{48}. This method does not require large ancillary systems nor large numbers of controlled-unitary operations, and thus is a promising avenue for eigenvalue estimation in the NISQ era.

To demonstrate the practical utility of VFF for eigenvalue estimation, we numerically computed the spectrum of a two-site Hubbard model in Fig. 6a–c and in Fig. 6d, we show eigenenergy estimates for a five-qubit XY model. Specifically, we used a one-clean-qubit (DQC-1) quantum circuit to discretely sample the function

$$g(t)={rm{Tr}}left(left|psi rightrangle leftlangle psi right| {e}^{-iDt}right)=mathop{sum }limits_{j}{e}^{-i{lambda }_{j}t},$$

(24)

where (left|psi rightrangle :=frac{1}{sqrt{2}}{(left|0rightrangle +left|1rightrangle)}^{otimes n}), and then used classical time series analysis to estimate the eigenvalues. This is achieved by computing each spectral estimate *S*(*λ*) with respect to a discrete number of values for time variable, ({t}_{j}={0,ldots ,{t}_{max }}) in increments of 0.2Δ*t*. A higher number of discrete points results in a better resolution of *S*(*λ*). The signal processing uncertainty principle constrains the spectral widths (variance in the estimate of *λ*) to obey ({sigma }_{{lambda }_{j}}{t}_{max }ge c), where *c* is a constant of order 1. In Fig. 6a–c, we show three examples with successively better optimization and, hence, longer integration times, ({t}_{max }). We plot eigenenergy estimates of diagonalized two-site Hubbard models with parameters *τ* = 1 and *u* ∈ {0.0, 0.2, …, 1.0} ranging from weakly to strongly coupled models.

The time series analysis extracts an estimate for the spectrum of energies corresponding to *V*, the approximate unitary given by VFF of the target unitary *U*, up to a global phase. The Hoffman–Wielandt theorem^{49} gives a bound on the total variation distance between spectra of *U* and *V*, in terms of the two-norm that in turn is directly related to the VFF cost function. For the concrete bounds we refer to Supplementary Note 3. This illustrates that the estimated spectral differences resulting from classical time series analysis give better approximations to the energy differences of the target Hamiltonian, when decreasing our VFF cost function. In Fig. 3c, d, we provide additional numerical analysis supporting this feature.