Hamiltonian operator approximation for energy measurement and ground state preparation

The Hamiltonian operator plays a central role in quantum theory being a generator of unitary quantum dynamics. Its expectation value describes the energy of a quantum system. Typically being a non-unitary operator, the action of the Hamiltonian is either encoded using complex ancilla-based circuits, or implemented effectively as a sum of Pauli string terms. Here, we show how to approximate the Hamiltonian operator as a sum of propagators using a differential representation. The proposed approach, named Hamiltonian operator approximation (HOA), is designed to benefit analog quantum simulators, where one has direct access to simulation of quantum dynamics, but measuring separate circuits is not possible. We describe how to use this strategy in the hybrid quantum-classical workflow for performing energy measurements. Benchmarking the measurement scheme, we discuss the relevance of the discretization step size, stencil order, number of shots, and noise. We also use HOA to prepare ground states of complex material science models with direct iteration and quantum filter diagonalization, finding the lowest energy for the 12-qubit Hamiltonian of hydrogen chain H$_6$ with $10^{-5}$ Hartree precision using $11$ time-evolved reference states. The approach is compared to the variational quantum eigensolver, proving HOA beneficial for systems at increasing size corresponding to noisy large scale quantum devices. We find that for Heisenberg model with twelve or more spins our approach may outperform variational methods, both in terms of the gate depth and the total number of measurements.

A different strategy to the state preparation relies on applying a non-unitary operation to cool down effectively the higher energy states. This can be done by representing them as a linear combination of unitaries (LCU) [53]. Being a valuable technique used in large-depth protocols for state-of-the-art Hamiltonian dynamics simulation [54] and linear system of equations [55], the value of LCU approaches was also assessed for near-term devices and analog simulators where unitary evolution is an available resource. In this way by measuring wavefunction overlaps one can perform inverse quantum iteration [56], quantum Krylov iteration [57], and quantum filter diagonalization [58] to study low energy properties of correlated materials and molecules. This approach was connected to the quantum version of time grid methods, that were used for the Schrödinger equation in the past [59]. From the variational protocols perspective, the LCU ansatz was also applied to chemistry [60] and linear algebra problems [61], notably showing the ability to avoid vanishing gradient regions.
While overall increasing requirements for the circuit depth and qubit overhead to perform SWAP test [87], the quantum time grid methods provide a way to reach the global energy minimum when starting from a suitable initial state. As these approaches rely on simulation of unitary dynamics for the system, they can be seen as relatives of the quantum phase estimation algorithm (QPEA) [63][64][65][66]. However, quantum time grid methods use the classical post-processing of overlaps and unlike QPEA do not require implementation of the controlled unitary dynamics. This makes them suitable for noisy large scale quantum (NLSQ) devices, where doubling the system size is a less pertinent problem than challenges of the variational search. Time grid methods are particularly favorable for analog quantum simulators [67,68] where one has the access to the simulation of quantum dynamics and overlap measurement [69,70]. Importantly, in this case variational approaches and Hamiltonian averaging are not applicable, and quantum time grid approaches offer the way to study low energy physics.
In this paper we propose a method to approximate a Hamiltonian operator with a sparse sum of unitary propagators. This serves as a building block for many LCU protocols, and is required when measuring the energy of the system in variational approaches, performing direct Hamiltonian iteration [71], implementing Lanczos algorithms, measuring the density of states [72] etc. We use a distinct differential decomposition approach, where Hamiltonian action is simulated with term-by-term quantum evolution, and the resulting energy expectation is read out as a sum of normalized overlaps. This strategy favors analog quantum simulation devices, where energy measurement through string-based Hamiltonian averaging is typically inaccessible. Comparing the proposed approach to the variational quantum eigensolver for the spin-1/2 Heisenberg model at increasing system sizes, we find that Hamiltonian operator approximation combined with quantum filter diagonalization offers better performance for lattices with more than 12 qubits, suggesting that the strategy can be used for ground state search in cases where variational methods suffer from the vanishing gradient problem.

II. METHODOLOGY
Our goal is to represent a Hamiltonian operatorĤ as a sum of unitary operators, taking as few terms as possible. Typically, this is done by partitioning the Hamiltonian into groups of qubit string operators (mutually commuting inside each group), as described in Hamiltonian averaging technique used to estimate energy of the system [14,73]. However, certain Hamiltonians, for instance in quantum chemistry, have unfavourable scaling for the number of groups (partitions) to be measured, as it grows like ∼ O(N 4 ) (N is a number of qubits) for vanilla Hamiltonian averaging [25]. Instead, we use the differential representation based on the quantum evolution operator (propagator) denoted asÛ(t) = exp(−iĤt), where t is a parameter which defines evolution time,Ĥ is a Hermitian time-independent operator, and we use = 1 throughout the paper. Acting with the time derivative operator on the propagator, we can formally write dÛ/dt = (−iĤ)Û(t). The Hamiltonian operator then readsĤ defined using its derivative at some fixed point of time t 0 , and it is convenient to use t 0 = 0. We approximate the propagator derivative using the finite difference scheme with S stencil points. The accuracy of the approximation scales as O δt S−1 [74], where δt is the distance between neighboring points. The derivative is approximated as a sum of unitary operators where s is a shifting parameter that is arbitrary in general, and we usually choose it to be (S − 1)/2 for symmetry reasons. To find the expansion coefficients q n (s) we decompose our function at stencil points t 0 + nδt (n = −s, −s+1, ..., S −s+1) into Taylor series around t 0 , [75]. Forming equations for each stencil point, we keep only the first S terms in the expansions. Next, we need to compose a linear combination of these equations such that the coefficient before the first derivative is equal to 1, and coefficients before all the other terms on the right are equal to zero. For this we solve the eigenvalue equation M ·q = δ, where the matrix M reads T is a vector of coefficients we want to find, and we set δ = 0, 1, ..., 0 T (as the first derivative term is considered). By inverting the matrix we find the required coefficients q n (s) and by inserting Eq. (2) in Eq. (1) we write the S-stencil differential approximation for the operatorĤ aŝ We refer to this procedure as Hamiltonian operator approximation (HOA).

III. RESULTS
We proceed simulating complex quantum systems and study their low energy properties using the proposed Hamiltonian operator approximation approach. As test systems we choose several paradigmatic material science models. These include Heisenberg model of spin-1/2 particles in the external magnetic field, and stronglycorrelated hydrogen chain as a standard example from quantum chemistry. For the simulation we choose the programming language Julia and use the Yao.jl package [76] as a simulation backend, capable of performing quantum protocols with state-of-the-art efficiency [76][77][78].

A. Energy estimation
As one of useful applications forĤ dynamical approximation we consider the measurement of the expected energy. Usually this is done through the procedure of Hamiltonian averaging, as commonly used in VQE [14]. We propose to use a hybrid quantumclassical approach, where energy is measured as a combination of wavefunction overlaps. Using HOA [Eq. (4)] we can write the Hamiltonian expectation as Ĥ ≈ (i/δt) S−s−1 n=−s q n (s) ψ|ψ(nδt) , where energy estimation requires calculating the overlap between |ψ and the evolved state |ψ(nδt) = e −iĤ(nδt) |ψ . Importantly, since we use the approximated derivative of the evolution operator, HOA approach favors short time evolution (see discussion below).
First, we consider the Heisenberg Hamiltonian of an N -qubit chain with open boundaries,Ĥ = −J whereX j ,Ŷ j ,Ẑ j are Pauli operators acting at site j. The results of simulation are shown in Fig. 1(a,c) for N = 8 and h/J = 0.1. We consider a uniform quantum state created by the string of Hadamard operators as |ψ = (|0 + |1 ) ⊗N /2 N/2 , and measure energy expectation by evaluating terms using noiseless statevector simulation and analog unitary evolution. In Fig. 1(a) we plot the difference between true expected energy and HOA expectation ∆E as a function of δt, plotted for increasing number of stencil points. We observe that the difference decreases as O(δt S−1 ), until it reaches numerical precision-impacted region (see the full discussion on the error scaling in Appendix A). In Fig. 1(c) we plot the energy difference for varying S and observe exponential improvement in energy difference, which also requires smaller δt (and correspondingly smaller circuit depth). We note that this monotonous dependence holds for the infinite number of measurement shots, and the finite shots case is discussed later.
The number of required overlaps to estimate remains small even for increasing S, and the complexity of HOA is defined by the maximal propagation phase T = Sδt. While an analog simulation of dynamics highly benefits the proposed approach by construction, we can also use digital simulation of dynamics. Asymptotically optimal simulation of quantum dynamics can be achieved with qubitization [79] or LCU-based approaches [54], at the expense of increased system size due to ancillary register. A resource-frugal alternative for many near-and midterm applications can be the Trotterization approach, which was used previously for quantum time grid methods [56,58]. Recent advances in Trotterized simulation of quantum dynamics suggest that scaling can be improved to O(N T /r + N T 3 /r 2 ) for r steps [80], and for local Hamiltonians there is a strong evidence that circuits can be further optimized [81,82]. For the Heisenberg model Trotterization can be readily performed by alternating couplings on even/odd sublattices [83,84], and implemented in various hardware platforms.
To understand the scaling of HOA with Trotterized quantum evolution we study the contribution of errors coming from the product formula truncation and differentiation (see Appendix B for the full discussion). As HOA favors small propagation time to reduce differentiation error we observe that the number of required Trotter steps remains small for increasing system size. For N = 13 spin-1/2 Heisenberg model, δtJ = 0.1, and S = 25, this allows for negligible Trotterization error already at r = 5 steps. The gate count scales as 23N r + 3N and gives 1534 gates. For a hundred-qubit system with similar settings we get ∼ 10 4 gates, which may become possible for improved NLSQ devices. Next, we proceed studying quantum chemistry problems, and consider molecular hydrogen (open) chains H n with homogeneous bond distance d between n atoms. Specifically, we use examples from QunaSys competition obtained using the STO-3G minimal basis set in the spinful case [85]. Qubit-encoded Hamiltonians are then obtained using Jordan-Wigner transformation from OpenFermion [86]. The chain of four hydrogen atoms H 4 at d = 1.0Å bond distance is encoded using 8qubit Hamiltonian. We evaluate the energy expectation of the Hartree-Fock state, providing benchmarking in Fig. 1(b,d). We find that the overall results are similar to the Heisenberg case, suggesting the internal structure of the Hamiltonian does not have a major impact on HOA accuracy.

B. Measurement and noise
To assess the full performance of the Hamiltonian operator approximation, we account for other imprecision sources coming from the measurement (finite number of shots) and effects of the environment (qubit decay).
The overlap measurement can be performed using several methods. The most popular choices correspond to SWAP and Hadamard tests [87,88]. The SWAP test for N -qubit system requires additional N + 1 qubits, where one qubit is used for the read-out [60]. A destructive version of the SWAP test can be performed using 2N qubits [69,88]. The Hadamard test uses a single ancilla added to the register, but doubles the circuit depth [58,88]. We also note that in certain cases direct ancilla-free measurement is possible, for instance when circuits possess a specific structure [89,90] or one of the eigenstates is known [56]. Further advances in the area include overlap estimation through randomized measurements [91], atom interferometry [70], or recent projected kernel techniques [92] and shadow tomography [93]. The latter was proven to be optimal for predicting fidelities and entanglement entropy. We consider the improvement of overlap measurement as an important direction for advancing Hamiltonian approximation techniques.
We simulate the full measurement scheme using the SWAP test to estimate the overlap between the initial state |ψ and the propagated state. Following the steps in [60] we prepare a superposition 1 √ 2 |ψ |0 + iÛ |ψ |1 for the system register and an ancilla. We apply Hadamard gate to the ancilla (last qubit), and measure the expectation forẐ N +1 , getting Ẑ N +1 = −Im ψ|Û |ψ . The results are shown in Fig. 2 for the case of the Heisenberg model with h/J = 0.1, S = 5, and δtJ = 0.5. We plot the difference between true energy of the uniform state in the system with HamiltonianĤ and the energy of the uniform state calculated with HOA, ∆E, shown as a function of the number of measurement shots N meas for which Ẑ N +1 is averaged. We observe that the upper bound of error decreases approximately as ∼ N −1/2 meas (Fig. 2a). Fixing N meas = 10 9 we calculate the dependence of the error on the time step δt. In the absence of noise the results are presented by blue dots in Fig. 2b), in line with the variance scaling ∝ c/N meas (c being a system-dependent constant) found in [60]. While for relatively large δt the error behaves similarly to ideal statevector simulator, the sampling noise increases the error at δtJ < 0.3, meaning that optimal δt depends on N meas . This suggests that resolving a small difference between two states (initial and propagated) becomes difficult at very small δt.
Next, we simulate the effect of noise using the wave function Monte-Carlo approach [94]. For this we run N iter trajectories of unitary evolution, subject to probabilistic action of jump operators √ γ(X j + iŶ j )/2 that can de-excite each qubit with decay rate γ. We simulate the Heisenberg model example and plot the energy deviation as a function of δt (Fig. 2b, magenta bars). The effect of decay is pronounced at large times δtJ > 1, where a significant number of jumps leads to erroneous overlap estimates at impacted trajectories (we consider large decay rate γ/J = 0.1, and simulate 10 9 trajectories). As δt and maximal propagation time T decrease, so does the jump probability, leading to improved energy precision ∆E. However, for small δtJ < 0.1 we find that even a small number of jumps increases the estimate imprecision, given that each collapsed trajectory leading to corrupted overlap estimate enters with δt −1 weight. We observe that optimal δt ≈ 0.2J −1 lies close to N measdefined case, depending on decay rate γ, and we note that a small time step region is the overall preferred choice for HOA.

C. Direct Hamiltonian iteration
We generalize our consideration to approximate higher powers of the Hamiltonian operatorĤ k . This can be defined through the k-th order derivative of the propagator, d kÛ /dt k = (−iĤ) kÛ (t), and following the numerical differentiation at S stencil points we express it aŝ The first nontrivial and useful example corresponds to the approximation ofĤ 2 operator, allowing for the straightforward measurement of the variance. This is often required for variational schemes, as together with the energy it may help to detect the convergence of the algorithm [95]. Hamiltonian averaging through Pauli string measurement has daunting scaling in this case. We note that to get the same order of accuracy as for the expansion ofĤ, we just need to increase the number of stencil points by one. This holds when going beyond k = 2, where each power increment requires at least one additional stencil point, and we generally prefer using an odd number of stencil points for the interval to be symmetric. However, the scaling becomes unfavorable for small δt, with error growing as O(δt S−1 ), thus shifting the optimal δt region.
We proceed using HOA as a part of the direct iteration algorithm to prepare the dominant eigenstate of the Hamiltonian matrix. We start from an initial state |Ψ 0 having a non-vanishing overlap with the ground state.
We simulate the repeated application ofĤ such that at iteration step k we get |Ψ k =Ĥ |Ψ k−1 / Ĥ |Ψ k−1 , where the denominator accounts for the state normalization. At sufficiently large k → K we get |Ψ K = H K |Ψ 0 / Ĥ K |Ψ 0 ≈ |ground state , and corresponding expected energy approaches the ground state energy. We apply the described procedure to prepare a low energy state for the Heisenberg model at the critical point h/J = 1 (Fig. 3).
Starting from the product state corresponding to a meanfield solution, we lower the energy by one order of magnitude using just four iterations. We also observe that when h/J > 1 the convergence further improves. At the same time, we note that direct iteration typically is an unstable procedure with critical dependence on the condition number, and works best for diagonally dominant matrices.

D. Quantum filter diagonalization with Hamiltonian operator approximation
Another application where our algorithm provides significant improvement is the quantum filter diagonalization (QFD) proposed by Parrish and McMahon [58]. The goal of the QFD procedure is to find the low energy spectrum of the system. Importantly, the procedure allows estimating both the ground state energy and excited state energies using the same resources (overlap measurements), thus being an efficient alternative to constrained variational methods [46,88,96]. In QFD one uses the trial wavefunction as a sum of 2k max +1 propagated states where {|ψ j } is a set of initial states and κ is a spectral width parameter, which generally shall be greater than the difference of maximal and minimal eigenvalues |E max − E min |. Next, the variational Rayleigh-Ritz approach is applied to find coefficients c j,k such that the energy functional E = Ψ|Ĥ |Ψ / Ψ|Ψ is minimized. This corresponds to solving the generalized eigenvalue problem where overlaps are calculated using a quantum circuit. SubstitutingĤ with the differential approximation (4) we avoid the Hamiltonian averaging procedure performed through Pauli string measurements, and the only difference for overlaps in LHS and RHS is the shift in evolution time by (k −k )/κ. As setting the time step for QFD requires the knowledge of the spectral width κ, in the following we use the Gershgorin circle theorem to provide its estimate in a scalable way [58]. The results of QFD using Hamiltonian operator approximation are shown in Fig. 4. First, we consider H 4 hydrogen chain (N = 8 qubits) for different bond distances d = 0.5Å, 1.0Å, 2.0Å and benchmark the difference between the true ground state and the variational state as a function of k max (Fig. 4a). We observe fast convergence to the true ground state, where starting from the Hartree-Fock initial state we can reach high precision (< 10 −5 Ha) with only few components in the ansatz. We find that QFD results with the ideal Hamiltonian (dashed curves) and HOA (solid curves) only deviate at larger k max , and strong correlations at larger bond distance lead to a slower convergence rate. As k max increases the expressibility of the ansatz improves, both due to larger number of states and longer propagation time. However, at large k max the solution of generalized eigenvalue problem may be challenging due to instabilities, and care must be taken to choose it in the optimal way. In Fig. 4b) we confirm that for fixed QFD parameters the quality of HOA improves with growing S and decreasing δt. Going to the larger scale example of N = 12-qubit H 6 hydrogen chain (d = 1.0Å), we again see exponential convergence to the groundstate energy with k max increase, reaching chemical precision ∆E = 10 −3 Ha already with five propagated states (Fig. 4c).
Finally, we consider the challenging example of a spinful fermion lattice described by the Fermi-Hubbard model. This is described by the Hamiltonian that includes on-site Coulomb repulsion U between opposite spins, nearest-neighbour hopping J, chemical potential µ (see the definitions and description in OpenFermion package [86]). We consider a minimal two-dimensional lattice with four sites, and use the Jordan-Wigner transformation to write the N = 8 Hamiltonian, where we additionally break the symmetry between up and down spin components with a weak effective magnetic field h. Choosing our initial state to be the uniform state, and staying at the half-filling, we observe that QFD+HOA approach can gradually bring the system towards the low energy state (Fig. 4d). Importantly, HOA performs so well that no significant difference between approximation and the ideal Hamiltonian can be spotted, as the two curves overlay. Given the huge progress in analog quantum simulation of Fermi-Hubbard lattices with cold atoms [68] and the possibility to perform interferometric measurements [69], this suggests a promising route towards studying exotic fermion phases.
We note that the problem of preparing exact ground state of the many-body Hamiltonian is QMAcomplete [97]. This is defined by the overlap between initial state and the ground state, that may be exponentially small. However, the same concern holds for quantum phase estimation algorithms. Choosing physically motivated initial states that allow for efficient ground state preparation is an important goal for future studies in the field.

IV. DISCUSSION
In this section we discuss the prospects of using Hamiltonian operator approximation with near-and mid-term devices, and compare them to state-of-the-art techniques in the field.
VQE comparison. First, we compare the QFD+HEA approach to the variational quantum eigensolver, being the standard NISQ tool. For this task we choose the problem of the ground state energy search for the Heisenberg ring at increasing system size. Specifically, we con-  sider an additional magnetic field, such that the system is at the critical point (h/J = 1) with highly entangled eigenstates. For the VQE we use the state-of-the-art optimization tool-set relying on analytic gradients (automatic differentiation) [98,99] and adaptive gradient descent [100]. It is known that automatic differentiation helps to mitigate the noise as compared to numerical differentiation techniques for updating the variational parameters [37]. Specifically we used the Adam optimizer [100] widely used in machine learning, and enables the efficient search for optimal variational angles.
We also note that unlike chemistry problems that have a possibility to use the unitary coupled cluster ansatz [15], material science problems typically rely on the hardware efficient strategies, as employed for Heisenberg model in [13,101]. This corresponds to the ansatz choice as layers of generic rotations (concatenated parametrized rotationsR Z -R X -R Z ), followed by the network of CNOTs. Another possibility is to employ alternating operator ansatz, and specifically the Hamiltonian variational ansatz (HVA) with individual tunable blocks based on separate Hamiltonian terms [102]. This can be seen as a physically motivated ansatz for spin systems.
The results for the statevector VQE simulation are shown in Fig. 5(a,b). We consider Heisenberg rings with the number of spins ranging from 6 to 14, and use the hardware efficient ansatz [13,101]. We optimize variational angles for maximum 3000 iterations, setting Adam's learning rate to the highest-performing value of 0.001. We observe that for smaller system sizes (N = 6, 8) VQE can prepare the state being close to the true ground state, also staying within relatively shallow depth [ Fig. 5(a)], and correspondingly small budget [ Fig. 5(b)]. As N grows the convergence to ground state requires ever increasing depth, and the corresponding gate count for reaching an approximate ground state at N = 14 surpasses 6000 gates (depth of ∼ 100). While we do not observe yet the zero gradients regime, the ground state preparation requires a largely overparametrized circuit. Previously VQE with the depth of 73 VQE was simulated for the same Heisenberg model at criticality for N = 7 qubits [101]. We perform additional tests with HVA (see details in Appendix C), and observe similar convergence.
Now, let us check the cost of the quantum time grid approach, corresponding to HOA+VQE, for the same system. We use the Trotterization approach detailed in Appendix B. We find that for small systems the gate count from Trotterization dominates the budget of QFD+HOA, making it more costly than VQE for N smaller than 10 qubits. However, at N = 12 it converges to the ground state with just 9 propagated states and below 2000 gates. For N = 14 the convergence holds, needing 2250 gates for the evolution. The weak scaling suggests that our approach has a larger constant as a starting budget, but outperforms VQE as the system grows to relevant sizes. We also note that our approach does not rely on optimization, and naturally avoids barren plateaus problem expected as the system scales to ∼ 20 qubit size [21]. These factors make the proposed approach suitable for NLSQ devices, while the point where it becomes beneficial depends on the task and hardware properties.
Comparing the measurement cost we see that for increased system sizes VQE needs over 3000 iterations for convergence. Taking N = 14 as an example, for each iteration the gradient measurement requires the number of independent runs being equal to six times the number of variational parameters (∼ 6 × 2000), times the number of measurement shots (needs to be > 10 5 for Adam, as shown in Ref. [37]). This leads to the measurement budget of 3 × 10 12 shots. Our approach requires in general N overlaps = (4k max + 1)S overlaps, being the number of time points of the QFD method multiplied by the number of points in the HOA method. This number decreases further if we adjust the time step of the HOA method to be equal to the time step in the QFD method. In such case many overlaps become degenerate, and the procedure only needs the calculation of (4k max + 1) + 4 unique overlaps. For N = 14 we calculate 21 unique overlaps, meaning that for the same measurement budget 10 11 shots can be used to reach extra-high precision on each run. Conducting additional studies to compare the Hamiltonian averaging with HOA, we note similar ∝ 1/N meas drop in the variance for observables and overlaps, with the latter requiring a constant overhead. The point at which the total measurement budget for HOA becomes favorable then depends on the number of variational parameters and number of VQE iterations, ultimately defined by the efficiency of the variational workflow.
QPEA comparison. Next, we compare the proposed approach to quantum phase estimation following [16]. The goal of QPEA is to determine eigenvalues of a unitary operatorÛ, when acting on the corresponding eigenstate. While the full QPEA circuit requires an extended ancillary register and quantum Fourier transform [65], the closest QPEA variation to quantum time grid methods is represented by Kitaev's phase estimation [63,66]. This requires only a single ancillary qubit and uses the adaptive protocol. However, to estimate the eigenvalue it uses the controlled unitary operation for the evolution operator, unlike separate overlap measurements. While in general the implementation of controlled unitaries is complex [106], for the simulation of dynamics we can rewrite it as a simulation of a modified Hamiltonian with the reduced locality. Taking Heisenberg model as an example, this requires simulation of an effective three-body terms involving the ancilla coupled to other spins. Using the Trotterized evolution, we can decompose each Trotter step into single qubit rotations, Hadamards, phase gates, and CNOT operations. For QPEA this translates to controlled two-qubit gates and Toffoli gates, adding a significant (∼ 10) constant overhead, which increases further if restricted connectivity is considered. Also, QPEA is not applicable to analog quantum simulators. The relation between approaches thus depends on the platform and available quantum resources.
Scaling. Finally, we ask the question: can Hamiltonian operator approximation become a viable strategy for the task of energy estimation in near-and mid-term future where larger systems are available but remain noisy? The issue of performing efficient energy measurement has recently gained attention [25,38,[41][42][43][44], and the advances are nicely summarized in Ref. [25], Table I. In particular, considering vanilla Hamiltonian averaging with commuting Pauli heuristics one gets a simple measurement circuit (constant depth of Pauli rotations), but pays for it with O(N 4 ) scaling for the number of partitions [14] (we call it depth-frugal methods). On the other side, the methods based on the basis rotation grouping have much better scaling for the number of partitions being O(N ), while requiring measurement circuits with the gatecount of N 2 /4 [42] (we say they correspond to partition-frugal methods). Another related partition-frugal technique corresponds to the unitary partitioning approach proposed in [39]. Hamiltonian operator approximation thus takes the ultimate position in partition-frugal methods, where the number of independent terms (overlaps) to measure is minimized, at the expense of increased depth if a digital Hamiltonian simulation is used. The practical utility of HOA then depends on the maximal propagation time T and the number of stencil points S being used. Numerically we find that this number does not depend on the system size or has weak, at most log(N ), dependence. The depth of the corresponding measurement circuit increases with T translating to gatecount O(N 2 T ) considering Trotterization [82]. However, the propagation time T = Sδt in HOA favors small values, and in many cases few Trotter steps suffice to reduce errors coming from the product formula (Appendix B). Searching for examples where Hamiltonian operator approximation offers the advantage is an important task for the future research.

V. CONCLUSIONS
We proposed the Hamiltonian operator approximation technique that allows representing a HamiltonianĤ as a linear combination of unitary operators. Using numerical differentiation rules we rewriteĤ as a sparse sum of quantum propagators, and benchmark it for relevant problems of energy estimation and ground state preparation. We found that the expected energy of the quantum system can be measured with high precision once we have access to the simulation of its dynamics, as for instance available in analog and digital quantum simulators. This holds in the presence of imperfections, including shot noise and decoherence. We also showed how the Hamiltonian operator approximation incorporates naturally in quantum Krylov-type approaches, and prepared the ground state of 12-qubit H 6 hydrogen chain using the quantum filter diagonalization. Comparing the proposed approach to variational quantum eigensolver and Hamiltonian averaging we see it can become beneficial both in terms of gatecount and total number of shots for the increasing system size. However, the cross-over point depends on the task and quantum platform.
Note added.-During completion of this work, we became aware of the independent work [107] that has been carried out in parallel.

Appendix A: Approximation errors
Hamiltonian operator approximation relies on differentiation of the evolution operator, and introduces errors that depend on the finite differencing procedure. Our goal is to bound an approximation error as Ĥ − H diff (S, δt) < , where is a pre-defined constant. This shall be achieved for a minimal product of the step size and the total number of points in the differentiation grid, Sδt. Additionally, we note that small S expansion is generally favoured due to its simplicity, and larger δt helps at the stage where physical implementation errors are introduced. The approximation error is a function of S, δt, and the Hamiltonian structure. Using the expansion procedure discussed before, for the infinite-precision arithmetic case the truncation error reads where C S is a coefficient depending on the expansion. This suggests that the approximation improves as S increases, and small δt is beneficial. Similarly, we can write the bound for an expectation value of the Hamiltonian operator, and introduce a required expected energy precision as ψ|Ĥ |ψ − ψ|Ĥ diff (S, δt) |ψ < energy . We further note that the scaling for HOA is in fact more involved as we approach the limit δt → 0 [74] and finite precision arithmetic is used. In this case the round-off error becomes important, and difference estimation is limited by minimal tolerance * energy / Ĥ , typically of the order of 10 −16 . For instance, in the case of quadratic approximation of derivative the total scaling reads as Ĥ −Ĥ diff (1, δt) < ε appr + ε num , where ε appr = 1 6 (δt) 2 max t d S (exp[−itĤ])/dt S is the approximation error as in Eq. (A1), and ε num = ( * /δt) max t exp[−itĤ] comes from the finite precision. We find that at small δt the finite-precision rounding error ε num starts to dominate, and there exists δt * such that the sum of the two contributions is minimized. However, we note that this limit is physically infeasible when the full measurement procedure is considered (see the discussion in the main text).

Appendix B: Trotterization errors
We aim to understand the scaling of Hamiltonian operator approximation for systems with Trotterized simulation of dynamics at increasing system size. For this we compare errors coming from the differential representation and Trotterization.
Specifically, we choose to simulate quantum dynamics digitally using the second-order Trotterization approach [80], and taking the Heisenberg Hamiltonian as an example. The Hamiltonian evolution for time τ can be  where equality is valid for Jτ 1 and r 1. For the chains with the periodic boundary we set indices as N + 1 = 1.
In Fig. 6 we show how the error caused by Trotterization (Fig. 6a) compares to the overall error of the HOA (Fig. 6b). We see that as the propagation times for HOA itself are relatively small it is only a small number of Trotter steps needed for the error caused by Trotterization to become much less than the HOA error. For 5and 8-qubit Heisenberg rings the error caused by Trotterization is small already at r = 1. For a larger system of 13 qubits with the same HOA parameters it is enough to have 5 Trotter steps for the Trotter error to become negligible. In the simulation for Fig. 6 we use values δtJ = 0.1 and S = 25, but note that typically these can be smaller, leading to even more pronounced decrease of the Trotterization error.
When considering the energy measurement as a part of the ground state preparation process we see that the maximal evolution time increases. In QFD+HOA this is posed by larger k max as the system grows. In this case we expect increasing importance of the Trotterization error, that translates to larger gatecounts.

Appendix C: VQE with Hamiltonian variational ansatz
In the main text for the VQE comparison we use the hardware-efficient ansatz based on layers of arbitrary single-qubit rotations and CNOTs. Additionally, we check the performance of the Hamiltonian variational ansatz for improving VQE convergence. The HVA circuit includes layers of Hamiltonian terms corresponding toX jXj+1 ,Ŷ jŶj+1 ,Ẑ jẐj+1 , andẐ j evolution. Addi-tionalR X andR Z rotations are included to improve expressivity. We used stochastic gradient descent in the Adam form with the learning rate of 0.001 and 3000 iterations to achieve convergence. The results are shown and discussed in Fig. 7. We observe the initial decrease of variational state energy as increasing depth (d ∼ N ), but note that circuits are not expressive enough to prepare the ground state with high fidelity. After the depth of twenty layers VQE reaches the barren plateaus, and halts the efficient search. With further increase of depth the overparametrization of ansatz improves the search, reaching high quality solution at large depth. We note that small gradients in this regime require increased number of shots to navigate the derivative-based optimization.