{'fraction_non_alphanumeric': 0.03709411414821424, 'fraction_numerical': 0.011760888960986293, 'mean_word_length': 4.390265228796269, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Variational Adiabatic Gauge Transformation on real quantum hardware for effective low-energy Hamiltonians and accurate diagonalization

Effective low-energy theories represent powerful theoretical tools to reduce the complexity in modeling interacting quantum many-particle systems. However, common theoretical methods rely on perturbation theory, which limits their applicability to weak interactions. Here we introduce the Variational Adiabatic Gauge Transformation (VAGT), a non-perturbative hybrid quantum algorithm that can use nowadays quantum computers to learn the variational parameters of the unitary circuit that brings the Hamiltonian to either its block-diagonal or full-diagonal form. If a Hamiltonian can be diagonalized via a shallow quantum circuit, then VAGT can learn the optimal parameters using a polynomial number of runs. The accuracy of VAGT is tested trough numerical simulations, as well as simulations on Rigetti and IonQ quantum computers.


INTRODUCTION
Low-energy approximations permeate many-body physics. Systems as diverse as cold atoms in optical lattices [1], solid state spin systems [2], and even current superconducting quantum computers [3] are accurately described by low energy theories. Powerful theoretical methods have been developed to obtain low-energy Hamiltonians from perturbative expansions, such as the Schrieffer-Wolff transformation [4], or non-perturbative methods [5]. However, due to the exponentially large Hilbert space, numerical calculations rapidly become unfeasible when the dimensionality of the Hilbert space increases, while analytical results are limited to toy models.
Quantum computers and simulators [6,7] are starting to become experimentally available, even in the cloud [8]. It has been shown that quantum computers can accurately approximate the ground state of many-particle systems [9][10][11], estimate molecular energies [12], molecular docking configurations [13], and even some excited states [14]. One of the main challenges in quantum simulation is computing the dynamics of quantum many-particle systems without having to resort to exact diagonalization or conventional perturbation theory. Algorithms based on the Suzuki-Trotter decomposition [7,15] have been adapted to better exploit the capabilities of current noisy hardware [16], while accurate evolutions for longer times can be obtained using variational fast-forwarding [17] or Hamiltonian diagonalization [18,19] methods.
Here we introduce a hybrid variational quantum algorithm for either block-or full-diagonalization of N qubit Hamiltonians, where complex calculations in exponentially large Hilbert spaces are performed on a quantum hardware while, when some assumptions are met, the classical part of the algorithm scales polynomially in N . Our algorithm provides an efficient way of finding a variational circuit that brings the Hamiltonian to the desired diagonal or block-diagonal form, which can be used to extract low-energy interactions or estimate quantum dynamics as in fast forwarding. Our method is based on recent advances [20][21][22][23][24] in the context of adiabatic gauge potentials (AGPs), which are infinitesimal generators of a unitary transformation diagonalizing a given Hamiltonian. The AGP is non-perturbative, it can recover Schrieffer-Wolff transformation in the perturbative limit, and it is tightly connected to the Wegner Hamiltonian flow [25], also called similarity renormalization group [26]. We define the Variational Adiabatic Gauge Transformation (VAGT), as a variational quantum circuit approximation to the unitary generated by the AGP, and show that the parameters of that transformation can be efficiently trained on noisy-intermediate-scale quantum (NISQ) devices [6].
We show that VAGT yields an accurate block diagonalization with few variational parameters (low depth) when the Hamiltonian has some energy separated or symmetry separated blocks, and full-diagonalization when the number of parameters increases. Finally, the feasibility of our method on current noisy quantum hardware is tested with experiments on the Rigetti Aspen-9 and IonQ 11qubit quantum processors. RESULTS We focus on the diagonalization, or blockdiagonalization, of a Hamiltonian H, assuming that there is another Hamiltonian H 0 whose eigenvalues and eigenstates are known, and possibly easy to prepare on a quantum device. We thus split H as where λV = H − H 0 , and λ models the strength of the correction. A good approximation of the ground state of H can be prepared thanks to the adiabatic theorem [27], by starting from the ground state |g 0 of H 0 and then slowly increasing the interaction strength λ. For an evolution time T , the approximate ground state is obtained as |g = T exp −i T 0 H µ(t) dt |g 0 , where µ(t) is a function, typically linear in t, satisfying µ(0) = 0 and µ(T ) = λ. Such adiabatic preparation of the ground A generalization of the adiabatic ground state preparation is given by the adiabatic gauge potential [20,21,28], which defines the infinitesimal generators of a unitary transformation that allows the estimation of more eigenvalues, in some cases even performing full-diagonalization -see also Appendix A for more details. Consider some infinitesimal generators A µ for µ ∈ [0, λ], set the unitary

Energy levels
and the rotated Hamiltoniañ with T ν denoting the ordering with respect to ν. Depending on the problem, we want to reach a diagonal or block-diagonalH µ at the end of the evolution, when µ = λ (see Fig. 1(a)). We assume for the sake of simplicity that H 0 has two blocks which, in the easiest case, correspond to two degenerate eigenvalues h P and h Q of H 0 , with h P h Q , that are possibly split by the term λV , as in Fig. 1(b) -the extension to Hamiltonians with more blocks is straightforward. In such case, we call P and Q, respectively, the projectors on the low-energy and highenergy sectors, but in general P and Q can also be projectors on symmetry sectors of the Hamiltonian H 0 . For block-diagonalization we should impose that all the elements of the off-diagonal blocks are zero, i.e. PH µ Q = 0. Differentiating such equation with respect to µ, we get The latter equation can also be written as [U † µ G µ U µ , H 0 ] = 0 when H 0 = h P P + h Q Q has only two degenerate eigenvalues, as in Fig. 1.
The adiabatic gauge potential is a particular choice of A µ that satisfies the operator equation see Appendix A for more details. Such equation is stronger than Eq. (4) and, when exactly satisfied, the re-sultingH µ has no off-diagonal elements. Approximations of the above exact solution were proposed in [20,21,[29][30][31][32], based on a variational approximation of the A µ , with optimal parameters obtained by variationally minimising, on a classical computer, either [G µ , H µ ] or G µ , where · is the Hilbert-Schmidt norm. With some assumptions, such a variational approximation of the AGP is efficient in suppressing matrix elements between states that belong to different energy sectors or that are wellseparate in the basis defined by the symmetries of H 0 , effectively resulting in a block diagonalization of the Hamiltonian. Therefore, the AGP can also be applied when low-energy and high-energy sectors are a priori unknown. In this work we propose a different variational approach, whose parameters can be optimized in NISQ hardware. Taking inspiration from the success of hybrid variational quantum algorithms [33,34], we consider a variational quantum circuit ansatz for the Adiabatic Gauge Transformation (AGT) U µ defined in Eq. (2), namely we write where α µ are variational parameters, L is the number of layers in the circuit ansatz, and B are local operators. Given the available gates in current quantum hardware we choose B such that e iαB is either a one-or two-qubit gate. If the Hamiltonian H 0 is diagonal in the chosen basis, then U 0 = 1 1, otherwise we assume that U 0 may be efficiently expressed as a known quantum circuit. Each parameter α µ is a continuous function of the running parameter µ ∈ [0, λ]. By dividing such interval in T steps δµ we create a discrete set of T values for µ: As a result we now have a discrete set of LT variational parameters {α t }. Within precision δµ the potential A µ = i(∂ µ U µ )U † µ at step t can be approximated via finite differences as We set α 0 = 0 at step t = 0 for 1 ≤ ≤ L and, starting from this initial configuration, we iteratively impose either Eq. (4) or (5), to get firstly α 1 and then the optimal parameters α t at all steps t. More precisely, as we will clarify in the next sections, setting β t = (α t+1 −α t ) δµ those equations can be written as X t · β t = Y t , for some operators X t and Y t . Callingβ t the solution of such operator equation we get the gradient-like update rule where allβ t are obtained by classical post-processing of quantum measurement results. We notice that, although Eq. (9) resembles a gradient ascent update rule, it was obtained from a completely different route. Parametric quantum circuits like the one in Eq. (6) can give rise to barren plateau in the cost function landscape [35][36][37] when the parameters are randomly initialized or for global cost functions. However, in the VAGT algorithm the cost function is local and all the parameters are initialized to zero and then evolved to the optimal values, a strategy that has been found to address the barren plateau problem [37,38]. An update rule similar to Eq. (9), namely based on the solution of linear system of equations with coefficients estimated via quantum hardware, was discussed in the context of the quantum imaginary time evolution algorithm [33,39], but the resulting circuits are entirely different. In the following sections we will study different applications that can be done efficiently on a quantum hardware.

Variational Quantum Adiabatic Gauge Transformation Algorithm
Since H acts on a Hilbert space whose dimension exponentially increases with the number of qubits, in general the diagonalization or block diagonalization of the Hamiltonian is exponentially hard. Here we show that, provided H can be accurately diagonalized by a shallow circuit, such diagonalizing unitary can be found in polynomial time using a hybrid quantum-classical algorithm that can be run on nowadays NISQ devices. Our algorithm is based on the minimization of the norm G µ that, as we show in Appendix B 1, is equivalent to the solution of the linear system X l, t β t = b l t of L equations, from which we can update the variational parameters following (9). In order to define a quantum circuit to measure the coefficients X l, t and b l t in a quantum computer, we first expand V and H µt in term of Pauli operators where σ j ≡ σ k j 1 ⊗ · · · ⊗ σ k j N are strings of Pauli operators, k j i ∈ {0, x, y, z} and, in principle, the sum index runs up to 4 N where N is the number of qubits. However, in most physical relevant cases the Hamiltonian only contains a limited number of terms, so most coefficients h jt and v j are null. The quantum algorithm we are about to define does not require to run a quantum circuit corresponding to such zero coefficients. In Eq. (10) we also have dropped the dependence on µ t of h j to simplify notation.
We call computational basis the basis in which all the Pauli operators σ z k are diagonal. In terms of such coefficients we find where the sums are restricted to non-null values of h kt and v j . The detailed derivation of the last equations is reported in Appendix B 2, where we also show how such quantities can be estimated on a quantum computer using the circuits given in Fig. 2, where |φ N is the maximally entangled state that can be constructed using O(N ) operations as in Fig. 2(c). Therefore, the number of operations in each circuit is at most O(N + L).
The scaling efficiency of the method depends on the connectivity of the Hamiltonians H 0 and V , namely on the number of terms N V and N H associated to a non-null coefficient in the two quantities in Eq. (10). Suppose that max(N V , N H ) = O(N γ ): for instance, if the Hamiltonians contain just single-qubit terms, then γ = 1; for nearest neighbour interactions γ = 1, too; on the other hand, γ = 2 if H or V contain all possible two qubit interactions. For each step t the number of circuits needed to evaluate all terms (11) and (12) where L is the number of layers in Eq. (6). Therefore, the number of measurements to be performed on the quantum device to calculate all variational parameters through Eq. (9) is while the solution of all linear systems of equations is at most O(L 3 ) for each step t. Therefore, for shallow circuits with L = poly(N ), both the algorithmic and measurement complexities scale polynomially in the number of qubits N and linearly in the number of steps T .

Low-energy approximation
When well-defined energy sectors exist for the problem at hand, as in Fig. 1, the projected block diagonalized Hamiltonian PHP can be interpreted as a low energy effective Hamiltonian, that can be useful in the context of many-body physics, where the original Hamiltonian may be unmanageable for many purposes, such as calculating dynamics. Assume that P is known, and that the lowenergy block can be expanded in the Pauli basis as whereσj are Pauli operators acting on N eff qubits, then the expansion coefficients h ef j can be obtained using a simple quantum circuit. Indeed, using the decomposition (10) we get and such coefficients can be evaluated in-hardware using a simple circuit like the one in Fig. 2. Indeed, suppose that P = 1 1 N eff ⊗ |π N −N eff π| and thatσj nontrivially acts only in the P subspace, then we may write which can be measured using the circuit shown in Fig. 2(d). In the above equations U λ ≡ U λ (α opt ) and α opt are the optimal parameters obtained at the end of the iteration, namely with t = T . Therefore, provided that the number N P of non-null expansion coefficients h ef j in Eq. (14) is suitably small, the effective Hamiltonian can be efficiently obtained. We test our framework using the following three-qubit model Hamiltonian with h = −5. For λ = 0 the Hamiltonian is diagonal in the computational basis with only two degenerate, well separated energy levels. As λ grows, the off-diagonal part of H is designed to completely remove the degeneracy, while keeping the levels in two separate subspaces, as showed in fig 3(b). In the Hamiltonian Eq. 17, qubits 1 and 2 do not interact directly, but can effectively communicate via qubit 3. Qubit 3 is the one that define energy sectors, so we can easily identify the low-energy projector P in the computational basis as where |0 is the eigenstate of σ z 3 with eigenvalue +1. This means that, at the end of the process, one can obtain a effective low energy Hamiltonian that couples qubits 1 and 2, with interactions mediated by qubit 3 without having to take in account its evolution at all. Even if this is just a toy model, it is reminiscent of quantum communication schemes [40][41][42]: if qubit 3 is replaced by a multi-qubit communication channel, forming for example a qubit chain, then this method can be used to find an effective Hamiltonian for the sender and receiver qubits only [43].
In Fig. 3 we present the results obtained with λ = 1 and T = 100. We use a variational ansatz composed by two blocks of layers: the first one is made by three layers of single qubit rotations around the x, y and z axis respectively for each qubit; the second one is made of parametrized two qubit gates, σ x ⊗ σ x , σ y ⊗ σ y and σ z ⊗ σ z , for different pairs of qubits. Since the Hamiltonian Eq. (17) is symmetric with respect to the exchange of qubits 1 and 2, we employed a symmetric ansatz where each operator B in Eq. (6) satisfies [B , S 12 ] = 0, being S 12 the swap operator. Considering such symmetry and alternating and repeating each block three times, we get L = 36 free parameters. Fig. 3(a) shows the absolute value of transformed Hamiltonian components |H jk |, where the block structure due to energy sectors ( Fig. 3(b)) is clearly visible. The resulting effective interaction between qubits 1 and 2 is where only the terms larger than 0.1 have been shown, the full Hamiltonian can be found in Appendix D. Fig. 3(c) shows the state fidelity F 1 between the time evolved state according to the full Hamiltonian (17) and the one obtained from the effective model, together with the fidelity F 2 between the evolved and initial states where ρ 12 (t) = Tr 3 |ψ 123 (t) ψ 123 (t)|, |ψ 123 (t) = e −itH λ |ξ 12 , 0 3 , ψ eff 12 (t) = e −itH eff |ξ 12 , |ξ 12 are randomly generated two qubit state, and t is varied from 1 to 1000. As shown Fig 3(c), F 2 displays a non-trivial behaviour, signaling a non-trivial dynamics. Since F 1 ∼ 1 for t ≤ 1000, such dynamics is accurately reproduced by the effective model for remarkably long times.

Block diagonalization
We now study the performance of the VAGT algorithm with symmetry defined blocks, by focusing on the following spin chain Hamiltonian with open boundary conditions where N is the number of qubits, while h and λ are, respectively, the transverse and longitudinal fields. For λ = 0 the model is exactly solvable, so eigenvalues and eigenvectors can be found in O(poly(N )) time. For this Hamiltonian the blocks can be identified as different magnetization sectors, as H 0 commutes with i σ z i , without necessarily being far away in energy. In Fig. 4 we present the results obtained with N = 4, h = 4.5, λ = 1 and T = 100 by numerical simulation of the quantum circuits. The VAGT in Eq. 6 is composed of two blocks of layers, the first block contains two layers of parametrized single qubit rotation gates around the x and y axes, while the second block contains two layers of parametrized two qubit gates. For the latter we choose only nearest-neighbour σ y ⊗ σ y and σ z ⊗ σ z interactions, and we alternate and repeat both blocks of layers ten times, resulting in L = 140. Even if such ansatz is obviously not universal for a 4-qubit system, the small off-diagonal terms at the end of the optimization, as shown in the left panel of Fig. 4(a), confirm the validity of our algorithm. The solution can also be improved by using deeper circuits and finer slicing, i.e. higher T . In Fig. 4(a), we also show (right panel) the transformed Hamiltonian in the magnetization basis, where the block structure associated with the different magnetization sectors defined by the original symmetry is clearly apparent.

Implementation on NISQ devices
We now discuss the implementation of our algorithm on real quantum hardwares, the Rigetti Aspen-9 quantum processor with 31 qubits, and IonQ quantum processor with 11 qubits, that we access through the cloudbased Amazon Braket service [8]. In order to simplify the experiment, we focus on two-qubit Hamiltonians, as in such case, as shown in Appendix C, we can fully exploit some specific properties to minimize the number of gates employed, and accordingly the simulation cost. With this simplification, valid for N = 2, each circuit requires at most 5 qubits. In order to fully exploit Aspen-9's 31 qubits and reduce cost, we run 4 different experiments in parallel, still guaranteeing the presence of one or two "garbage" qubits between different experiments to reduce possible cross talk. On the IonQ's 11-qubits hardware we run instead 2 experiments in parallel. In numerical simulations the full circuits shown in Fig. 2 is implemented. Using a universal variational ansatz the circuit depth is L = 15, but lower depths are possible by using an ansatz suitably designed for the specific Hamiltonian problem at hand. For this purpose, we choose to test our method on quantum hardware with the highly non-diagonal Hamiltonian defined below: with λ = 1 and where the v k coefficients are randomly chosen between 0 and 1. Due to the lack of well defined sectors, whether they are defined by energy, magnetization or other physical quantities, we do not expect block-diagonalization in any basis, though thanks to the universal variational ansatz, we can expect full diagonalization in the computational basis.
In Fig. 5 we show the results of our numerical simulation and hardware experiments [44]. Fig. 5(a) shows the exact energy levels for different λ and the energy levels obtained via VAGT with T = 10 discretization steps and S = 100 measurement shots. We see that, in spite of the finite discretization steps, finite measurement shots, and imperfect gate implementation, results on the IonQ hardware are very accurate, while simulations on Aspen-9 did not converge. We run different experiments on Aspen-9, always getting similar outcomes, though numerical simulations with Rigetti's decoherence and dephasing times show results comparable with IonQ. We believe that the high errors on Rigetti's hardware are possibly due to the qubit connectivity, that requires extra compilation steps in order to implement the non-local gates required by the VAGT circuits. On the other hand, all qubits in IonQ's hardware are fully connected, so better results are expected. Indeed, we see in Fig. 5 that the accuracy obtained with IonQ hardware is very high.
In order to study the accuracy of the computed diagonal forms, rather than focusing on operator norms or other mathematical distances, we focus on the study of physical quantities, like time-evolved correlation functions. In Fig. 6 we plot the real part of the time-evolved correlation function on the ground state: where α = x, z, |g 0 is the ground state of H 0 , so that U λ |g 0 is our approximation of the ground state of the Hamiltonian;H λ is the Hamiltonian in its diagonal form, reconstructed from either experimental or simulation data. In Fig. 6 we see that all simulations provide an accurate description of the dynamics, with the exception of the simulations on the Rigetti hardware.

DISCUSSION
We have defined the VAGT hybrid quantum algorithm for block-and full-diagonalization of many-body Hamiltonians. It can be used to extract low-energy effective theories in complex many-particle systems or to approximate long-time evolutions, e.g. using fast forwarding.
The VAGT is based on the adiabatic gauge potential (AGP), a non-perturbative method that generalizes the adiabatic theorem to multiple energy levels. The AGP has been successfully used in both analytical calculations with toy models and numerical simulations with classical computers, which are nevertheless limited to few-body operators, because of the exponentially large Hilbert space. The VAGT algorithm on the other hand is specifically made for hybrid quantum-classical simulations, where the complex calculations in exponentially large spaces are efficiently performed by the quantum hardware. It uses a variational quantum circuit approximation of the unitary transformation generated by the AGP, whose optimal parameters are iteratively obtained by merging outcomes from purpose-built quantum measurements with simple classical post-processing routines. When a Hamiltonian can be transformed into a (block)  24), where the AGP is obtained using T steps and L layers, using the same model and parameters of Fig. 5. The Analytic red line is referred to a completely classical simulation of our method, where quantum measurements are replaced with analytically computed expectation values, the Simulation green line is referred to a classical simulation of our algorithm on an ideal quantum processor, where we implemented also a quantum measurement process simulation, using S = 100 shots. Note that an ideal quantum processor is, by definition, unaffected by noise. The Aspen-9 and IonQ purple and brown lines show quantities computed from experimental data, obtained using real quantum hardware. diagonal form using a shallow parametric circuit, then the VAGT algorithm can find the optimal parameters efficiently, using a number of classical and quantum operations that scale polynomially in the number of qubits.
We remind that the algorithm relays on the use of suitable ansatz for the choice of the shallow parametric circuit: How to find such an ansatz, or even understand if it can exist for a given target Hamiltonian, are obviously rather relevant questions, which are out of the scope of the present paper.
To show the performance of the VAGT algorithm, we have considered both random and physically motivated Hamiltonians, where the block structure may come from separated energy bands or may be defined by the symmetries of the problem. We have performed both exact numerical simulations and simulations with realistic error sources (e.g. measurement shots), always obtaining convergence after a few iterations. Moreover, we have also run our algorithm on Rigetti and IonQ quantum computers, finding very accurate results on the latter, possibly thanks to its all-to-all qubit connectivity.
Appendix A: Adiabatic Gauge Potential Consider a system evolving under the Hamiltonian H µ(t) defined in Eq.(1), which is time-dependent through the parameter µ(t) and is written as a matrix in the computational basis where all σ z j are diagonal. We define U µ(t) the instantaneous unitary operator that diagonalizes the Hamiltonian at time t with respect to an eigenbasis of H 0 , and we'll always denote with ∼ the operators in this reference frame. For the moving observer in the instantaneous eigenbasis of H µ(t) the effective Hamiltonian ruling the dynamics is [21]:Ĥ whereÃ µ(t) is the Adiabatic Gauge Potential (AGP) in the moving frame:Ã the AGP in the standard frame is: It is possible to show [21,45] that: where is the adiabatic, or generalized force, operator [21], n µ(t) being the istantaneous n-th eigenstate of H µ(t) .
Let's now suppose we want diagonalize the Hamiltonian H µ(t) at any time t. Instead of calculating directly the unitary operator U µ , we can search for its instantaneous generator, the AGP. We use a variational approach, in the sense that we make an hypothesis of a suitable form of A µ , A µ (α), depending on some variational parameters α. Let's now define: (A7) from equation (A5), if there is a set α * of variational parameters such that G µ = −F ad we also have where D µ is a operator that commute with H µ . In other words, if we find the set of optimal variational parameters α * , we find the AGP apart of diagonal elements in the basis of the instantanous eigenstates of the Hamiltonian. In [45] it is formally demonstrated that searching for the operator A µ that minimize the distance from G µ and F ad , that is searching for the best approximation of the AGP, is equivalent to finding the variational parameters that minimize G µ (α) operator norm: We remind that solving this equation leads to the best approximation of the AGP, except for the diagonal part, that is undetermined by construction, as we explicitly stated by equation (A8).

Appendix B: Our variational ansatz
Taking inspiration from variational hybrid quantumclassical computation [33,34], we propose to use a "quantum circuit"-type ansatz of the operator U µ : where L is the number of layer in the circuit ansatz, the B operators are one-or two-local operators (loosely speaking: e iαB is a one or two-qubit gate); and the arrow over the product sign define the order of the product itself: specifically: Note that in this paper we supposed to know the eigenvalues and eigenstates of the H 0 Hamiltonian defined in the main text in equation (1), meaning we can efficiently construct the quantum circuit realizing the rotation U 0 such that is diagonal, and consequently it is a constant element in the definition of our variational circuit ansatz above. We can compute the generator of U µ , that will be our variational hypotesis for A µ , from its definition in Eq. (A4):  (1) we can assume that µ is a running parameter, µ ∈ 0...λ, and iteratively find the generator A µ for all points µ ∈ 0...λ.
At this level, each parameter α l µ is a continuous function of the running parameter µ. If we now divide the interval 0 . . . λ in T intervals δµ, we create a discrete set of T values for µ: As a result, we now have a discrete set of variational parameters, {α} L,T l,t ; U µt at the step t is expressed as a parametric evolution, depending on variational parameters α t , and the expression of the generator at the step t is now: where we use the finite difference form for the derivative of α wrt µ and we have defined: Since the Hamiltonian H 0 in equation (1) is diagonalized by U 0 , we can take advantage of the fact that U (µ t = 0)| {α}=0 = U 0 : this means that we know that the optimized parameters {α} L,N l,t=0 at the step t = 0 are all zero. Morover, A µt (α t , α t+1 ) is a function of only the subset of {α} L l,t and {α} L l,t+1 at the step t and t + 1 respectively. This two facts leads to an iterative method to solve eq. (A9): • first compute eq. (A9) at the step t = 0 with {α} L l,t=0 = 0, that we know are optimized already ; • this gives a function of only {α} L l,t=1 , the subset of parameters at t = 1, that can be easily optimized (see appendix B 1 for more details); • once the optimal parameters for t = 1 have been obtained, one has to repeat the previous two steps in order to obtain the optimal parameters for t = 2 and so on, until one reaches the last step t = T , that corresponds to µ = λ and the correct generator A(λ) for the unitary operator U (λ) is finally obtained.

Analytic minimization
As we said previously, at the time step t we want optimize the cost: w.r.t. α t+1 , since α t 's are optimized already. Using the form (B7) and the equation (1) for the Hamiltonian, we find: where we use the fact that V is Hermitian and we defined: In order to minimize the cost (B11)we can now calculate the gradient and set it to zero: The equation (B14) is a linear system that we can solve numerically: once the solution for is found, we can write the optimized parameters α t+1 ∀ as

Quantum circuits
For every step t ∈ 1...T the method involves the solution of the linear system (B14) so, at each step, we need to compute: ∀ , l, ∈ 1...L, where we dropped the t step index to simplify the notation.
In order to evaluate them through a quantum computer, let's suppose: where σ k are strings of Pauli operators ∀ k that forms a complete basis of SU (N ) (eventually, some coefficients v j and h j may be zero, depending on the specific model). Our quantities become: where we use also equations (B9) and (B12), again dropping t-labels.
We are now left with the task of evaluating on a quantum computer the following two types of terms: For this aim, we consider the identity below: where A and B are hermitian operators acting on a 2 Ndimensional Hilbert space H, T indicates the transpose and |φ is the maximally entangled state defined as and {|i } is a orthonormal basis for H. Using (B22) we can now express our target terms as: Note that σ j are strings of Pauli operators, so σ T j = (−1) n y j σ j , where n y j is the number of σ y in σ j . Furthermore We can choose a variational ansatz where (B ) T ≡ B ∀l ( in other words, we can use an ansatz in which σ y appear an even number of times in each operator B ), so: where V is our variational ansatz U , defined in equation (B8), taken in the reverse order, i.e.: (of course one can choose an ansatz in which B l too is not symmetric, in this case a factor (−1) n y l , where n y l is the number of σ y in B l , must be inserted in the exponential in the equation above, as well as in equation (B25)).
Our target terms can now be written as: For the first type of terms, we can consider the quantum circuit of Fig. 2(a), where R x is: and the construction of the state |φ from the standard initial state |00...0 is made by applying the Hadamard gate on the first N qubit, followed by N C-NOT gates, controlled by the first N qubit with the second half of the register as a target. In particular, for N = 2 the quantum circuit is the one of Fig. 2(b), that constructs the state |φ and can be used as an input to the circuit above. Once the ancilla qubit is measured in the σ z basis, the probability p A (0) of getting the outcome 0 is linearly related to the target: Similarly, the circuit of Fig. 2(c) can be employed to evaluate the second type of terms, as their value are encoded in the probability of getting the outcome 0 from the first ancilla qubit measurement through the relationship Appendix C: Saving money: an approach for N = 2 Even if the algorithm proposed in the previous section is feasible on NISQ devices within reasonable limits, the real costs of running experiments on real devices can be high if a cost is charged for circuit reconfiguration, that is implicit at any step in our variational approach. Consequently, we develop an alternative method to reduce the number of quantum circuits employed by the algorithm itself: Although this method can lead to lower costs (and we actually used it in our experiments), we remark that it is not efficient from the point of view of scalability,so that it is practically useful for very small values of N only. Dropping again the step index t, let's consider the equation (B11): where we used definitions (B12) and (B15) for Q l and β l . We can expand the operators V and Q l ∀l on the basis of Pauli strings, obtaining expressions with at most 4 N operators (this expansion is exponentially inefficient, but for N = 2 it leads to a 16-terms expansion, that we can afford easily.): where σ i is a Pauli string of two operators and, by definition: After some calculations, recalling that where V is the vector of v i 's, Q is the matrix composed by q kj and β is the vector of β i 's. We recall that the cost function, as it's expressed in equation (C7) is the wellknown least-squares loss function of a multiple linear regression classical problem, that can be solved efficiently. Therefore, our goal is now to estimate efficiently the elements of V and Q defined in eq. (C4) and (C5), respectively.
For what concern y i , they are already known in our setting, since they are coefficients in Pauli decomposition of the operator V , the hard-to-diagonalize part of the Hamiltonian. Let's focus on q kj . Decomposing and taking into account equations (B12) and (B22) the quantity we want to estimate is where n y j is the number of σ y Pauli operators in the string σ j . We already know that we can estimate this quantity via the first type of quantum circuit presented in the section B 2.
Even for N = 2, the method presented above seems to perform worse than the one presented in section B 2. In fact, the number of circuit we have to execute, we have where N γ is the number of terms in the decomposition (C8), L is the length of the variational ansatz, T is the number of steps in the discretization of the parameter µ and the highly inefficient factor 4 N comes from the decomposition of Q k . On the other hand, by this method we only have quantum circuit of the type shown in Fig 2(a), that can be reduced. In fact, consider the figure 7: (a) panel shows the quantum circuit we have to execute in order to calculate q kj , while (b) panel shows a completely equivalent quantum circuit. The probabilities of getting the outcome 0 on the ancilla qubit are respectively: and the equivalence between them follows from identity (B22), and from now on we simply denote both of them with P ljk (0). As we can see also from Eq. (C11) the second circuit gives a probability connected with the commutator where F ljh are the structure coefficients of the algebra. Using the last equation in the (C11) we obtain: From the equation above it is clear that, although in principle one have to run all different N γ · 4 N circuits in fig 7(b) with different Γ lj , since the latter is just a commutator of two sigma strings there are not so many different results for Γ lj , and so there are not so many different quantum circuits one has to really run. Indeed, one can compute the structure coefficients F classically and run only 16 (4 N for N = 2) circuit of the type showed in figure 7(c), where σ h is one of the 16 elements of the SU (4) basis. Keeping trace of the original lj-th term corresponding to a given σ h one can recover the information about the original probability P ljk (0).
Note that, in the quantum circuit in figure 7(c) the rotation of the ancilla qubit right before the measurement process is different that in 7(a) and (b): and the output probability of getting 0 from the ancilla qubit is: and we finally recover the originally searched for probability via Finally, we remark that in quantum circuit in fig 7(c) it is possible to replace the indirect measurement with a direct one: in other words, it is possible to remove the ancilla qubit, replacing the control-B k and control σ h gates with measurements of the expectation value of B k ⊗ σ h on the principal register. This is possible only because both B k and σ h for all k and h's are Pauli strings, so they're observables.