Characterization and Verification of Trotterized Digital Quantum Simulation via Hamiltonian and Liouvillian Learning

The goal of digital quantum simulation is to approximate the dynamics of a given target Hamiltonian via a sequence of quantum gates, a procedure known as Trotterization. The quality of this approximation can be controlled by the so called Trotter step, that governs the number of required quantum gates per unit simulation time. The stroboscopic dynamics generated by Trotterization is effectively described by a time-independent Hamiltonian, referred to as the Floquet Hamiltonian. In this work, we propose Floquet Hamiltonian learning to reconstruct the experimentally realized Floquet Hamiltonian order-by-order in the Trotter step. This procedure is efficient, i.e., it requires a number of measurements that scales polynomially in the system size, and can be readily implemented in state-of-the-art experiments. With numerical examples, we propose several applications of our method in the context of verification of quantum devices: from the characterization of the distinct sources of errors in digital quantum simulators to determining the optimal operating regime of the device. We show that our protocol provides the basis for feedback-loop design and calibration of new types of quantum gates. Furthermore it can be extended to the case of non-unitary dynamics and used to learn Floquet Liouvillians, thereby offering a way of characterizing the dissipative processes present in NISQ quantum devices.


I. INTRODUCTION
Quantum simulation aims at 'solving' complex quantum many-body problems efficiently and with controlled errors on quantum devices [1,2]. The development of quantum simulators is presently considered to be one of the most promising near term goals in quantum technologies, where problems of practical relevance from condensed matter physics [3] to high energy physics [4] or quantum chemistry [5,6] can be addressed in regimes potentially beyond the abilities of classical computations. Quantum simulation can be implemented as analog or digital quantum simulators. In analog quantum simulation a many-body Hamiltonian of interest finds a natural implementation on a given quantum device. Examples are ultracold bosonic and fermionic atoms in optical lattices as Hubbard models [7][8][9][10][11], or spin models with trapped ions [12,13], Rydberg tweezer arrays [14][15][16][17][18], and superconducting qubits [19,20]. While analog simulators naturally scale to large particle numbers, there is limited flexibility realizing complex Hamiltonians required, e.g., in high energy physics and quantum chemistry. In contrast, in digital quantum simulation (DQS) the time-evolution is 'programmed' as a sequence of Trotter steps represented by quantum gates [21], experimentally realized in several platforms such as trapped ions [22][23][24][25], Rydberg atoms [26], or superconducting circuits [27][28][29]. While in the future fault tolerant and scalable quantum computers might be available to perform DQS, we are still in an era of NISQ devices, i.e., noisy intermediate scale quantum devices, where the elementary * L.P. and T.O. contributed equally to this work. quantum operations are performed as analog quantum gates without (the possibility of) error correction.
A central challenge faced by analog and digital quantum simulation on NISQ devices is thus the characterization and verification of the quantum device, i.e., a demonstration of its proper functioning as a prerequisite for quantitative predictions (see Refs. [30,31] for recent reviews on quantum verification). This is particularly relevant in regimes where classical simulations become intractable, which induces the need to develop efficient protocols for characterization of quantum simulation experiments that scale to large particle numbers.
Below we address this problem in the context of Trotterized DQS by utilizing recently developed techniques of Hamiltonian learning, which we employ to perform an efficient process tomography of Trotter blocks in DQS.

II. BACKGROUND AND MOTIVATION
In this section we provide background and motivation for our approach. We start with a brief review of Trotterized DQS and its relation to Floquet systems, and elaborate on motivation and techniques for Floquet Hamiltonian learning as efficient process tomography of Trotter blocks.

A. Trotterized Digital Quantum Simulation
DQS approximates the time-evolution operator e −iĤt generated by a time-independent many-body Hamilto-nianĤ with quantum gate operations available on a NISQ device. This can be achieved via Trotterization according to a Suzuki-Trotter decomposition [32,33], as Trotter blocks implemented in experiments are affected by control errors: they are treated as 'black boxes' to be characterized by our protocol. The protocol aims at learning the FHĤ exp F (τ ) from measurements performed on a set of initial states, after repeated application of the Trotter blocks. The measurements are used to fit an ansatz for a truncation ofĤ exp F (τ ) to a given order L in τ . The procedure is repeated for several values of τ . (c) The corresponding learning error is measured as a function of τ : if it scales as τ L+1 (for τ < τ * ), then the reconstruction was successful. The data presented in (c) are obtained from the protocol applied to a DQS of the Heisenberg model Eqs. (14)- (15), for learning up to order L = 1 of the FH. (d) The procedure can be generalized to learning Floquet Liouvillians L exp F (τ ) in presence of dissipation (see Sec. V). explained below. We are interested in DQS of generic many-body Hamiltonians with a local structure [34], i.e., we writeĤ with the Hamiltonian componentsĤ j being few-body operators. This structure implies thatĤ depends on a number of parameters that scales only polynomially in the system size. We emphasize that this notion of local-ity includes also Hamiltonians with long-range, few-body interactions. The time-evolution for a sufficiently small time step τ is carried out via Trotter blocksÛ τ . In their simplest form, these Trotter blocks are assembled from the Hamiltonian componentsĤ j [35] according tô as illustrated in Fig. 1(a). The overall coherent timeevolution e −iĤt at times t = nτ (n = 1, 2, . . .) is thus approximated by a sequence of Trotter blocksÛ n τ . This becomes exact when τ → 0, i.e., e −iĤt = lim n→∞Û n t/n . While the Trotter block in Eq. (2) approximatesĤ up to deviations of order τ , we refer to [36,37] for decompositions with higher order Trotter errors. The following discussion, and the methods developed in the present work, will apply to any Trotterization.
Trotterized DQS has the structure of a Floquet problem [38][39][40][41], i.e., a periodically driven many-body problem [42] with stroboscopic time evolution generated bŷ U τ . Writing the Trotter block aŝ we can identifyĤ F (τ ) with the Floquet Hamiltonian (FH) of the driven system. The stroboscopic time-evolution generated by DQS is thusÛ n τ = e −iĤF(τ )nτ , to be compared with the targeted time-evolution e −iĤnτ . Therefore we can interpretĤ F (τ ) as the physical or 'actual' Hamiltonian implemented for a given Trotter step τ on the quantum device.
The FH plays a key role in our experimental protocol for efficient process tomography of a Trotter block below, and thus we will briefly elaborate on its properties. For sufficiently small τ we can write the FH as a Magnus expansion [42]Ĥ and for reference we write out the first few terms explicitly for the Trotter decomposition Eq. (2) in Appendix A. HereΩ 0 =Ĥ is the target Hamiltonian of the DQS, while higher order termsΩ l (l > 0) represent Trotter errors, written as perturbation to the original Hamiltonian H.
In Eq. (4) we distinguish the regime where the Magnus expansion converges, i.e., the regime where DQS is a controlled approximation to the target Hamiltonian. On the other hand, divergence of the series leads to proliferation of Trotter errors at a critical step size τ * [43]. The existence of this Trotter threshold [39][40][41] is intrinsically connected to the implemented dynamics turning quantum chaotic for τ > τ * , where the spectrum ofÛ τ exhibits universal properties described by random matrix theory. This is a phenomenon extensively studied with simple model systems as the kicked top [40,44,45], and in the context of periodically driven many-body systems [44,[46][47][48].

B. Process Tomography of Trotter Blocks via Floquet Hamiltonian and Liouvillian Learning
Below we will be interested in efficient process tomography of Trotter blocks in an experimental setting, where the Trotter block is considered as a 'black box' to be characterized in an experiment. To this end we will apply techniques of Hamiltonian learning to infer the FH parametrizing the Trotter block.
Hamiltonian learning (HL) techniques were originally developed as protocols to infer local Hamiltonians of closed systems from steady states [49][50][51][52][53][54][55][56], as well as methods for reconstructing Hamiltonians governing the time-evolution in quench dynamics [34,[57][58][59][60][61][62] from measurements performed on the quantum device. The central assumption underlying HL is that typically implemented many-body Hamiltonians have an operator structure consisting of a small number of terms that scales polynomial in the system size. The prototypical HL protocol starts with an ansatz for the operator content of the Hamiltonian to be learned, and infers a best fit from comparison with experimental data. These protocols also have a built-in toolbox for error assessment in the learned Hamiltonian, allowing for a systematic exploration by enlarging the ansatz to include missing terms.
In the present context, we note that the target Hamiltonian Eq. (1) was assumed to have a local structure, and this property is inherited by any truncation of the FH (4) for τ < τ * (see Appendix A). This allows HL developed for quench dynamics to be applied to learn the FH implemented by a DQS.
In contrast to the ideal case of Eq. (4), in implementing Trotter blocks in an experiment there will be errors, e.g., control errors so that the experimentally realized Trotter blockÛ exp τ will differ from the idealÛ τ . More generally, quantum gates will suffer from decoherence, i.e., the experimental Trotter block must be described by open system dynamics as a Kraus map, instead of a unitary operator. This is schematically shown in Fig. 1(b), where the Trotter blocks are treated as 'black boxes' to be characterized.
Starting with the unitary case, our goal is thus to learn the experimental FHĤ exp F (τ ) defined viâ This is achieved by inferring the operator content of order-by-order in τ (for τ < τ * ). Specifically, we formulate an ansatz for a truncation of Eq. (6) to a chosen order τ L , and fit it to measurement data obtained from the output states |ψ(nτ ) = (Û exp τ ) n |ψ 0 for chosen initial states |ψ 0 , where experiments are performed for various n = 1, 2, . . . and τ [63]. As discussed in detail in Sec. III, the behavior of the learning error as a function of τ reveals whether a L-order truncation ofĤ exp F (τ ) has been learned: a learning error scaling with τ L+1 confirms a successful reconstruction, while the scaling with a lower power suggests a systematic extension of the ansatz by additional terms. This is illustrated in Fig. 1(c), where we show numerical data obtained from the application of our protocol to a DQS of the Heisenberg model (see Sec. IV A): an incomplete ansatz leads to a plateau of the learning error (pink curve) while scaling (red curve) is observed if and only if all the terms up to order τ L are included. The number of samples required for such a reconstruction scales only polynomially with the system size, as long as τ < τ * . Furthermore, monitoring the learning error as a function of τ provides a mechanism to detect the Trotter threshold τ * , as the point at which learningĤ exp F (τ ) fails: for τ > τ * , an ansatz consisting of a polynomial number of terms generally cannot approximate the FH, which is reflected in a sharp increase of the learning error (marked by the grey vertical line in Fig. 1(c)).
Note that by learning the operator content ofĤ exp F (τ ), and comparing it with the target HamiltonianĤ, both control and Trotter errors are expressed as perturbations to the original (target) HamiltonianĤ, and are thus directly amenable to physical interpretation; i.e., by learning the operator structure of these perturbations we expect direct insight regarding their physical effects in simulating quantum phases and quench dynamics [64]. We contrast this to measuring Trotter block or quantum gate fidelities by comparing the experimentalÛ exp τ with the targeted e −iĤτ , which provides only a global assessment of errors [65]. Finally, in an experiment the FH learning protocol can be run in a quantum feedback loop to finetune parameters as a way of engineering desired Trotter blocks, which will even work in a regime where the quantum device cannot be accurately calibrated [66].
While the above discussion has focused on coherent evolution, these ideas are readily extended to include dissipative processes (as schematically shown in Fig. 1(d)). Assuming that the dissipative errors can be modelled as a quantum Markov process, i.e., described by a master equation, we can learn in a similar way the operator structure of the corresponding Floquet Liouvillian L exp F (τ ), which also admits a Magnus expansion in powers of τ in analogy to Eq. (6) (see Eq. (24) in Sec. V).
The remainder of this paper is organized as follows. In Sec. III we introduce the technical details of our procedure for process tomography of Trotter blocks including order-by-order HL of Floquet Hamiltonians. Sec. IV is dedicated to several possible applications of this scheme: We illustrate the characterization of Trotter errors in Sec. IV A, the detection of imperfect implementation of Trotter blocks in Sec IV B, a verification of a DQS of the lattice Schwinger model in Sec. IV C as well as the design of many-body gates in Sec. IV D. Finally in Sec. V we extend our order-by-order learning scheme to the characterization of dissipative processes in the context of Floquet Liouvillian learning.

III. ORDER-BY-ORDER LEARNING OF THE
FLOQUET HAMILTONIANĤ exp F (τ ) In this section we describe the protocol for process tomography of Trotter blocks in DQS, based on quantifying the individual terms in the experimental Floquet Hamiltonian Eq. (6) order-by-order in τ . We refer to this as Floquet Hamiltonian learning (FHL). Before going into the details, we provide a brief overview of Hamiltonian tomography in quench dynamics, which serves as the basis for our protocol. While in the following we focus on Hamiltonian learning, i.e., assuming unitary dynamics, we will include dissipative processes in the context of Floquet Liouvillian learning in later sections.

A. Hamiltonian learning in quench dynamics
Hamiltonian learning techniques, as developed in the context of non-equilibrium many-body dynamics, form an important sub-routine of our protocol. In this paper we exploit different methods for learning local Hamiltonians from projective measurements as introduced in Ref. [34] and Ref. [52]. We start by discussing the scheme of Ref. [34], which outlines a method to recover the structure of a time-independent HamiltonianĤ that governs the unitary time-evolution in a set of quantum quenches |ψ i (t) = e −iĤt |ψ i (0) , by measuring local observables on the output states |ψ i (t) . Here, {|ψ i (0) } Ncon i=1 denote a set of (in general) arbitrarily chosen initial states (see below). The basic idea is to reconstruct the HamiltonianĤ implemented on the quantum device from a parametrized ansatz HamiltonianĤ which is composed of N A few-body operatorsĥ j chosen from a set A whose number scales polynomially with the system size. Within this ansatz class, the goal is now to find optimal parameters c rec which yield the best possible approximation to the physically implemented HamiltonianĤ from the condition of energy conservation. In particular, for a time-independent Hamiltonian we have the conditions Ĥ i,t = Ĥ i,0 for all i, in which Ĥ i,t = ψ i (t)|Ĥ|ψ i (t) . As proposed in Ref. [34], such an optimal approximation for the ansatz Hamiltonian can be determined by minimizing the energy differences between initial and final states over the parameters c, i.e.
This is a quadratically constrained problem where we optimize under the condition that c T c = 1. Note that the parameters are learned only up to an overall scale, which can be determined separately via time-resolved measurement of a single observable (see App. B). Every quench experiment starting from an initial state |ψ i (0) provides an independent constraint on the coefficients c, implying that in general N con > N A constraints are required for a unique reconstruction. The cost function (8) has the advantage that it can be readily expressed in matrix form using the expansion Eq. (7). In particular we can rewrite Eq. (8) as where the matrix elements M i,j = ĥ j i,t − ĥ j i,0 correspond to local observables to be measured on the output states |ψ i (t) . Minimizing (8) thus amounts to identifying c rec with the right singular vector corresponding to the smallest singular value λ 1 of M [34] Here a non-zero value of λ 1 is caused by the ansatz set A being insufficient to account for the operator content of the implementedĤ. Thus, we refer to λ 1 as the learning error, as its value becomes an indicator for the quality of the reconstructed c rec . The learning error will be the central quantity for order-by-order learning ofĤ exp F (τ ). exp l τ l , that is the generator of an experimental Trotter blockÛ exp τ . We will start from an ansatz for a truncation ofĤ exp F (τ ) in the form of Eq. (7), and attempt the reconstruction of the parameters c using the method above.
The central idea of our approach is to repeat the reconstruction for several values of τ with the chosen ansatz components {ĥ j }, and to measure the learning error λ 1 as a function of the Trotter step τ . The behavior of λ 1 (τ ) with τ is used to rigorously assess whether the given ansatz captures the operator content of allΩ exp l up to the chosen truncation order L. Specifically, the scaling λ 1 (τ ) ∼ τ L+1 ensures that the chosen ansatz captures all components ofĤ exp F (τ ) up to order τ L , thus verifying the reconstruction. The key element enabling this rigorous error assessment is the fact that a scaling λ 1 (τ ) ∼ τ L+1 implies an equivalent scaling for the distance between the implementedĤ exp F (τ ) and the reconstructed Hamiltonian H rec (τ ) ≡Ĥ(c rec (τ )), i.e., which we prove in Appendix C. Importantly, this result enables us to distinguish between Trotter errors and control errors that will appear in an experimental setting. For example, if an ansatz A encompassing the target HamiltonianĤ of the DQS leads to a constant value of λ 1 (τ ), we can conclude that additional terms, corresponding to experimental imperfections, need to be added to A. Once these terms have been added and λ 1 (τ ) ∼ τ , one can further improve the ansatz by adding (nested) commutators to learn higher order termsΩ exp l , thereby achieving an order-by-order learning ofĤ exp F (τ ). Our protocol thus comprises the following steps: (1) Expand the ansatz HamiltonianĤ(c) in a local operator basis according to Eq. (7). The operator basis A = A 0 ∪ A 1 ∪ · · · ∪ A L is first of all motivated from the structure of the terms appearing in the (theoretical) Magnus expansion, i.e., A l contains the l-nested commutators appearing inΩ l (see App. A). In an experimental setting, A should then include additional terms whose presence we want to test.
The Magnus expansion forĤ exp F (τ ) is guaranteed to converge only up to a critical value of τ [43], that is connected to the breakdown of Trotterization at the Trotter threshold at τ * [39][40][41]. In general, for τ > τ * the Floquet HamiltonianĤ exp F (τ ) does not exhibit a local structure, meaning that process tomography based on a local ansatz will fail which is reflected in a sharp increase of λ 1 . Therefore our protocol provides an experimentally feasible way of detecting the Trotter threshold in DQS, and thus probing the regime of validity of the Trotter approximation (see Sec. IV A).
In practice, the matrix elements in M (τ ) will be noisy due to the fact that they are estimated from a finite number of projective measurements. As we will show, the effect of a finite number of samples on λ 1 can be analytically estimated, and does not spoil the validity of the protocol. We will discuss this in more detail in the examples presented in the next section.
Among the several types of imperfections, static deviations from the desired resources/target Hamiltonians can be detected with our protocol by including the corresponding terms in the ansatz set A. Moreover, Eq. (8) is agnostic to state preparation errors as long as sufficiently many independent constraints are provided.
Note that while the presentation in Sec. III puts emphasis on learning the Hamiltonian in the context of coherent dynamics, the entire procedure can be extended to dissipative dynamics which we present in Sec. V. In such a more general Liouvillian learning procedure the Hamiltonian as well as the dissipative terms are parametrized by an ansatz in Lindblad form and the corresponding parameters are reconstructed from the minimization of an experimentally measurable cost function in analogy to Sec. III. Similar to Eq. (11) the scaling of the minimum of the cost function can be used to verify the quality of the reconstruction.

IV. FLOQUET HAMILTONIAN LEARNING IN DQS: NUMERICAL EXAMPLES
In this section, we present numerical results on some proposed applications of Floquet Hamiltonian learning. We do so in the context of a DQS of several experimentally relevant systems, i.e., the XXZ/Heisenberg model, the lattice Schwinger model and the 1D cluster model. All implemented Trotter gate sequences are built from the following set of resources for any pair of qubits i, j and with tunable parameters J, θ and rotation axis n, withR Here and in the following,σ µ j (µ = x, y, z) denote Pauli operators acting on the j-th qubit.

A. Disordered XXZ-Heisenberg Model
As a first example we consider process tomography of Trotter blocks for a spin-1/2 XXZ-model. In particular we study the scaling relation for λ 1 in Eq. (11) and show how one can detect the breakdown of Trotterization from the behavior of λ 1 as a function of τ . In a second step, we discuss the effects that a finite measurement budget has on λ 1 . Finally, we show that our protocol can be employed to learn disorder patterns in quantum manybody systems, which is important in the context of, e.g., many-body localization [67].
The target Hamiltonian considered here iŝ A possible Trotterization for this model readŝ where the individual terms can be generated from the resources (13) aŝ µ (θ) = j e −iθσ µ j being rotations on all qubits (µ = x, y, z). In this section, for FHL we choose a zeroth order ansatz and a first order ansatz comprising all the terms inΩ 0 andΩ 1 in the implemented Floquet Hamiltonian. The case of an 'incomplete' ansatz set will be treated in Sec. IV B. In Fig. 2(a) we show the behavior of λ 1 as a function of τ for A 0 (blue data points) and A 0 ∪ A 1 (orange data points), plotted for different measurement budgets, represented here by the number n s of measurements used for the estimation of each of the expectation values that appear in the constraint matix M . Let us first focus on the data for n s → ∞. There, for τ smaller than a critical value τ * marking the Trotter threshold, we can observe scaling of λ 1 with τ . A zeroth order ansatz leads to a scaling λ 1 ∼ τ , confirming the successful reconstruction of all zeroth order terms inΩ 0 , and the presence of (Trotter) errors proportional to τ . Similarly, a first order ansatz leads to a scaling λ 1 ∼ τ 2 indicating the presence of Trotter errors proportional to τ 2 . Analogous scaling can be observed in the parameter distance c ex (τ ) − c rec (τ ) between the vectors of reconstructed (c rec (τ )) and exact (c ex (τ )) Hamiltonian coefficients, shown in Fig. 2 this is a consequence of the scaling relation (11).
Sharp threshold behavior of λ 1 at τ ≈ τ * signals the Trotter threshold, marked by a vertical grey line in Figs. 2(a)-(b). For τ > τ * , the FH cannot be approximated by a local ansatz, which is signaled by the absence of scaling of λ 1 (τ ) and c ex (τ ) − c rec (τ ) . In this sense the behavior of λ 1 (τ ) can be used as an experimentally measurable signature for the proliferation of Trotter errors. For a deeper discussion of FHL in the context of the Trotter threshold and the transition to quantum chaos, we refer the reader to Appendix D and to our recent work on collective spin systems [45].
In contrast to the limit n s → ∞, a finite measurement budget leads to a plateau of λ 1 at small values of τ , as shown in Fig. 2(a). This is due to the presence of (simulated) measurement noise on the elements of the constraint matrix M defined in Eq. (12), that limits the resolution of λ 1 . The value of this plateau for a given measurement budget can be estimated as as shown in Appendix E. The quality of the reconstruction in presence of measurement noise, as measured by λ 1 or the parameter distance, is proportional to n −1/2 s , as shown in the inset of Fig. 2(b). The estimate in Eq. (16) is valid when errors due to an incomplete ansatz are smaller than the noise introduced by a finite measurement budget. Thus, comparing λ 1 (τ ) to the above estimate allows one to decide whether the ansatz needs to be complemented with additional terms, as we discuss in Sec. IV B below.
In Fig. 2(c)-(d) we show the reconstructed Hamiltonian parameters to zeroth (Ω exp 0 ) and first order (Ω exp 1 ), respectively, for a given value of the Trotter step τ . There, the spatial dependence of the coefficients is accurately reproduced for all operators appearing in the first orders of the FH (the dashed lines show the exact values). This demonstrates the ability of FHL in resolving the spatial structure of the Hamiltonian parameters, thereby learning the disorder pattern characterizing a given manybody system. The overall scale of the coefficients can be efficiently extracted using the method described in App. B.

B. Detection of Static Control Errors
In this section, we show how our protocol can be used to detect, and learn, possible errors affecting the resource gates in a DQS. The errors considered here are of two types: (i) deviations of the resource Hamiltonians from the expected ones, and (ii) static errors in the singlequbit rotations. In the latter case, our protocol also allows for estimating the optimal value of the Trotter step for which the DQS most faithfully reproduces the target Hamiltonian.
To understand the origin of these two types of errors in a specific example we consider the Trotterization for the XXZ/Heisenberg model used in the previous section, which we rewrite here for a translation invariant Hamiltonian asÛ The entangling gates e −iτ J µĤ XX are implemented as the time-evolution with the analog resource HamiltonianĤ XX over 'times' τ J µ : in realistic implementations, the operator content of H XX may slightly deviate from the expected one (i.e., jσ x jσ x j+1 ) by additional unknown terms. We refer to these terms as deviations from the expected resource Hamiltonians. Additionally, a DQS can be affected by errors in the qubit rotationsR µ (θ), resulting in generally small and site-dependent over-or under-rotations. We now address these two cases separately, showing how our protocol can be used to detect and learn the corresponding imperfections.

Static deviations from expected resource Hamiltonians
We consider an imperfect implementation of the XX resource Hamiltonian that appears in Eq. (17) of the form withV j,j+1 being an unknown two-body operator acting on qubits j and j + 1. To learn the termsV j,j+1 , we employ FHL with the aim of constructing an ansatz set A for which scaling λ 1 ∼ τ is observed down to the noise threshold of Eq. (16). We do so by repeatedly adding to A new operators to be measured, until this condition is met. Errors or unknown terms generated by physical processes typically correspond to few-body operators: this restricts the search for new operators to a space whose dimension scales only polynomially with the system size.
We illustrate this adaptive protocol in Fig. 3(a) and further apply it to a DQS of the translation invariant XXZ/Heisenberg Hamiltonian Eq. (14) in Fig. 3(b). The results are obtained for a system of N = 10 spins, with uniform couplings J x = J y = 1, J z = 0.75, B x = 0.5 in Eq. (14). We use 55 initial product states in random bases, each evolved to 6 different final times up to a total simulation time of 16. The error bars stem from measurement noise and from different initial states realizations.
While an ansatz based on the target Hamiltonian leads to a plateau of λ 1 at a value which is significantly higher than the noise threshold (the red dashed line in the figure), the addition of ansatz terms encompassing the unknown operators lowers this value. When the ansatz is such that the plateau value of λ 1 falls below the noise threshold, we can conclude that if additional unknown terms are present, they are small and cannot be resolved with the current measurement budget.

Rotation errors, and determining the optimal Trotter step
We consider now the case in which the rotations in Eq. (17) are affected by qubit-dependent over-or under-  (14) (yellow data), and a full zeroth order ansatz set A0 (blue data, which leads to scaling down to the noise threshold. The results are obtained for a system of N = 10 spins, with target couplings J x = J y = 1, J z = 0.7, B x = 0.75 in Eq. (14). Here ns = 5000. We use 55 initial product states in random bases, each evolved to 6 different final times up to a total simulation time of 16. The error bars stem from measurement noise and from different initial states realizations.
rotations with strength δθ, aŝ with δθ j uniformly sampled in [−δθ, δθ]. Such rotation errors become more important for smaller Trotter steps: the optimal Trotter step at which the DQS most accurately simulates the target Hamiltonian is determined by a tradeoff between rotation errors and Trotter errors, dominant for large τ . Here, we show that such optimal Trotter step can be determined using FHL.
In Fig. 4 we show an example of our scheme applied in presence of rotation errors. In Fig. 4(a) we show the parameter distance c targ −c rec between the reconstructed coefficients and the target ones c targ (the ones corresponding to the target Hamiltonian), as a measure of the accuracy of the DQS. Different colors correspond to different strengths δθ of the rotation errors. As can be seen in Fig. 4(a), for each δθ, there is a clearly visible minimum separating the regimes of dominating Trotter errors (large τ ) and dominating rotation errors (small τ ): this minimum can be used as an estimate for the optimal Trotter step for the DQS. As shown in Fig. 4(b), the same behavior is observed in the learning error λ 1 (τ ) for an ansatz encompassing the target Hamiltonian (14), i.e., A = {X, XX, Y Y, ZZ} (yellow data). An ansatz capturing all operators generated by the imperfect rotations leads to scaling down to the noise threshold, as expected.

C. Verifying a DQS of the Lattice Schwinger Model
Simulating the time-evolution of lattice gauge theories on digital quantum devices has recently attracted considerable experimental [24] and theoretical [68][69][70] interest. Here we numerically illustrate that FHL can be applied to verify the DQS of the lattice Schwinger model as implemented in [24]. In particular we demonstrate reconstruction ofĤ F (τ ) up to first order in τ in presence of the long-range interactions, that emerge from the elimination of the gauge degrees of freedom (see [24]).
The Hamiltonian of the lattice Schwinger model, acting as the target Hamiltonian of the DQS readŝ wherê In Fig. 5 we show the reconstruction error c ex − c rec and the corresponding reconstructed J zz rec (i, j), with J = w = m = 1 setting the units of energy. As shown in Fig. 5(a), for small τ ≈ 0.1 a zeroth order ansatz yields a low reconstruction error, while for larger τ Trotter errors become more dominant and require a first order ansatz to achieve a faithful reconstruction. When using a first order ansatz, the structure of the reconstructed J zz (i, j) resembles that of the target Schwinger model up to τ ≈ 0.6 as can be seen in Fig. 5(c). Moreover, the reconstructed parameters of the first order terms characterize deviations from the target Hamiltonian in the form of Trotter errors. At τ ≈ 0.8 Trotterization breaks down which is reflected in the learned coefficients in Fig. 5(d).

D. Design and Verification of Many-Body Gates:
Example of the 1D Cluster Hamiltonian In this section, we employ our order-by-order learning protocol as a tool to design multi-body interaction gates, i.e., gates coupling more than two qubits, from the resources in (13). Such interactions can be generated from higher order terms in the Magnus expansion, and their strength can be 'enhanced' by choosing a Trotterization canceling the lower orders. The cancellation of lower order terms generally requires the ability of reversing the sign of the time-evolution implemented by the resource gates. We consider here the situation in which this inversion can be achieved only approximately, leading to 'spurious' low-order terms which we want to correct for. We show that FHL can be integrated in a feedback loop optimization scheme for tuning the experimental control parameters to achieve this correction, thereby enhancing the higher order multi-body target terms.
We numerically exemplify such a protocol in a DQS targeting the one-dimensional cluster Hamiltonian, a prototypical example featuring three-body interactions [71]. Specifically, our target Hamiltonian readŝ A Trotter sequence implementing this model with the resources given in (13) readŝ wherê  Ncon for different values of the difference ∆J between forward and backward time-evolution, for an ansatz A1 containing only the first order terms. Here, ns = 2000. As long as ∆J > 0, λ1/ √ Ncon cannot show scaling for arbitrarily small values of τ (as shown by its increase for smaller τ ). Thus, λ1 can be used as a cost function to be minimized by tuning the experimental parameters controlling ∆J. As shown in (c), minimizing λ1 corresponds to minimizing the norm of the zeroth order terms c rec,0 ord. , which can be reconstructed by adding A0 to the ansatz. Once λ1/ √ Ncon falls below the noise threshold in Eq. (16), we conclude that within our resolution set by ns, ∆J ≈ 0. The results are obtained for a system of N = 10 spins, using 55 initial product states in random bases, each evolved to 6 different final times up to a total simulation time of 10/τ . The error bars stem from measurement noise and from different initial states realizations.
where the zeroth order terms in τ vanish. We point out that the addition of theσ z fields is a consequence of the Trotterization and is not essential: it can be easily canceled by addingσ z fields with opposite signs.
We consider the situation in which the coefficients J of the forward (Û XX τ andÛ Y X τ ) and backward (Û XX −τ and U Y X −τ ) evolutions slightly differ from each others, leading to a imperfect cancellation of the zeroth order terms XX and Y X. We denote with J b the coefficients of the backward evolution operators, and define ∆J = |J − J b |. While here this difference is the same for each gate, the following results also apply to the more general case of non-uniform deviations. As before, the ansatz set to zeroth order is denoted by A 0 = {XX, Y X}, and the ansatz set for the first order terms by A 1 = {XZX, Z}, i.e., containing the target terms. In order to assess whether the zeroth order terms in A 0 are canceled and A 1 captures the leading terms inĤ F (τ ), we measure the learning error λ 1 (τ ) and we show that its value can be used as a cost function to be minimized by optimizing the experimental parameters controlling ∆J. The protocol is summarized in Fig. 6(a), and reads as follows: (1) Start from an initial guess for the backward evolution parameters (resulting in a generally non-zero ∆J).
(2) Measure λ 1 as a function of τ for a first order ansatz A 1 .
(3) If the scaling is not satisfied down to the noise threshold (16), add A 0 to the ansatz and measure the corresponding zeroth order coefficients c rec 0 ord. . (4) Use the result to optimize the backward evolution parameters in order to minimize λ 1 , and thereby the norm c rec 0 ord. . (5) Repeat from 2. until scaling down to the noise threshold is observed.
The results of this procedure for the DQS of the onedimensional cluster model are shown in Fig. 6(b)-(c) for a system of 10 qubits, with J = 1 setting the units of energy. In Fig. 6(b) we show the behavior of λ 1 with τ for different values of ∆J. As expected, for a finite ∆J scaling cannot be observed down to arbitrarily small values of τ . Instead, as the Trotter step decreases, one observes (i) an initial decrease of λ 1 for larger values of τ , due to the first order terms in the Floquet Hamiltonian being larger than the (small) spurious zeroth order ones, and (ii) an increase of λ 1 for smaller values of τ due to the zeroth order terms becoming dominant. When ∆J becomes sufficiently small, the effect of the non-canceled zeroth order terms becomes too small to be resolved in the chosen range of Trotter steps. Additional details on this example, in particular the linear scaling of λ 1 (τ ), can be found in Appendix F. In Fig. 6(c) we show the behavior of λ 1 as a function of ∆J for a fixed value of τ , where it is shown to reflect the behavior of the norm of the zeroth order coefficients c rec 0 ord. . This confirms that the learning error λ 1 can be used as an experimentally measurable cost function for verifying and optimizing the implementation of target Hamiltonians.

V. FLOQUET LIOUVILLIAN LEARNING IN TROTTERIZED DQS
In this section, we extend our order-by-order FHL to the verification of DQS subject to dissipative processes, by characterizing the generator of Trotter blocks in terms of a Floquet-Liouvillian. In order to describe this approach, we follow the logic of section III by first introducing Liouvillian learning as the underlying routine of Floquet Liouvillian learning (FLL), where coherent as well as incoherent processes are parametrized by an ansatz in Lindblad form. Second, we demonstrate how the quality of the reconstruction can be verified by scaling relations analogous to Eq. (11) which we illustrate via numerical examples.

A. Order-by-order FLL
The situation we have in mind is a DQS that implements a Trotter sequence in analogy to Eq. (2), where the generator of each gate, instead of being a Hamiltonian H i , becomes a time-independent Lindblad superoperator L i whose action on a stateρ is described by with theL k being the jump operators describing the dissipative channels, and γ k the associated rates. TheĤ i andL k are generally few-body operators, meaning that physically implemented Lindbladians L i can be described by a polynomial number of parameters. In presence of dissipation, the role of the evolution op-eratorÛ τ is taken by a superoperator which generates stroboscopic time-evolution of any initial density matrixρ(0) according toρ(nτ ) = P n τρ (0). In the following, we describe a method for achieving the order-by-order reconstruction of the generator L F (τ ) of the stroboscopic dynamics.
For sufficiently small τ , L F (τ ) can be expressed as a Magnus expansion in powers of τ with Q 0 = i L i and Q l (l > 0) given by l-nested commutators of the L i in complete analogy with the unitary case. As a consequence, a finite truncation of the above expansion for L F (τ ) is still a local operator, and can be learned efficiently as we explain below. For more details on the Floquet Liouvillian, including an explicit formula for the commutator of any two Lindblad superoperators, see Ref. [72]. As in the FHL case, the starting point of the FLL protocol is an ansatz for the experimentally realized L exp F (τ ), which can be chosen of the form [72] specified by a set of operators A = A H ∪ A D = {ĥ j } j ∪ {ˆ k } k and depending on parameters c H and c D to be determined. To generate constraints for reconstructing the coefficients c H and c D we impose the validity of the Ehrenfest theorem, in the spirit of Refs. [52,73]. Specifically, for any observableÂ, initial stateρ(0) and final time nτ , it holds that with · t = tr(·ρ(t)), where the above integral is necessarily discretized since the integrand is measured only at stroboscopic times. Thus, the Liouvillian parameters can be determined by choosing N con > N A constraint operatorsÂ k and by minimizing the sum of the residuals This is a linear optimization problem which can be recast in matrix form as where the matrices G H and G D and the vector b are determined from measured expectation values as The role of λ 1 as learning error is taken by the minimum of the cost function in Eq. (28), which we denote by  As for FHL in Sec. III, we study ∆ as a function of the Trotter step for a fixed ansatz. We obtain a scaling relation analogous to Eq. (11) for the reconstructed Liouvil- which we prove in App. G. In analogy to order-by-order FHL in Sec. III, Eq. (33) can be used to assess whether the chosen ansatz encompasses the experimental Floquet Liouvillian.

B. Numerical Illustration of FLL
In this section we present a numerical illustration of FLL applied to the DQS of the Heisenberg model as given in Eq. (15), in presence of dissipation. Specifically, we replace the generators of each individual gateÛ i τ , with i = {X, XX, Y Y, ZZ} in Eq. (15), by a Lindblad superoperator L i , defined as in Eq. (22) with jump operatorŝ and corresponding dissipation rates γ for each spin j, with Eq. (23) being the associated evolution superoperator over one Trotter cycle. The dissipative processes in this example are chosen in order to illustrate the capability of FLL to learn multi-qubit jump operators. The characterization of such dissipative processes might give insight into the underlying physical processes and help improve the quantum device. (See for example Ref. [74] on how cosmic radiation affects the coherence time of superconducting qubits.) Our numerical results are presented in Fig. 7. The blue and red data in Fig. 7(a) show the behavior of ∆(τ ) for an ansatz missing and including the dissipative terms, respectively. In complete analogy to Sec. IV B, the missing ansatz terms result in a plateau of ∆(τ ) for small τ . Only when the dissipative terms are included in the ansatz the scaling is observed down to small values of τ , hence verifying the complete reconstruction to zeroth order. The red dashed line in Fig. 7(a) corresponds to the noise threshold, providing an upper bound on ∆(τ ) in presence of measurement noise (see Appendix E for a derivation). The blue data being larger than this bound implies the presence of terms not accounted for in the corresponding ansatz. In Fig. 7(b) we plot the reconstructed parameters for both Ansätze (blue and red dots), together with the exact parameters (purple crosses). Note that adding dissipative terms to the ansatz not only allows for reconstructing the corresponding coefficients, but also improves the accuracy of the learned Hamiltonian parameters.

VI. SUMMARY AND OUTLOOK
In this work we proposed an efficient method for process tomography of Trotter blocks in DQS as an experimental protocol for the characterization and verification of NISQ devices. The protocol is based on 'learning' either the Floquet HamiltonianĤ exp F (τ ) as generator of unitary evolution for a Trotter step τ , or the Floquet Liouvillian L exp F (τ ) in the case of dissipative dynamics. The Floquet Hamiltonian, and similarly the Liouvillian, admits an expansion in powers of the Trotter step aŝ exp l τ l . Owing to this structure, our learning protocol starts from an ansatz for the operator content of a truncation ofĤ exp F (τ ) or L exp F (τ ) to a chosen order τ L . The ansatz contains theoretically expected operators (e.g., the target Hamiltonian of the DQS) together with terms associated with implementation errors or dissipative channels that one wants to test for. The corresponding coefficients are reconstructed by fitting the ansatz to measurement data from a set of time-evolved states, for several values of τ . This procedure comes with an associated, measurable learning error λ 1 (τ ), whose behavior as a function of τ (i) reveals up to which order in τ the chosen ansatz reconstructs the implemented FloquetĤ exp F (τ ) or L exp F (τ ), (ii) provides means of systematically improving the ansatz if needed, and (iii) enables the detection of the Trotter threshold τ * , i.e., the point at which the Trotter approximation breaks down andĤ exp F (τ ) cannot be described by a power-expansion in τ . This protocol is scalable, with a number of required measurement runs scaling polynomially with system size (see also Appendix H).
With our procedure, Trotter errors, control errors as well as dissipative processes affecting an experimental DQS can be learned and distinguished from one another. These direct insights into the microscopic processes underlying the implemented dynamics offer a potential route for compensation of errors via fine tuning, e.g., via a feedback loop.
Furthermore, the present techniques provide ingredients for the design of complex quantum gates or manybody interactions (see Sec. IV D) and subsequent verification of the elementary Trotter blocks. A natural extension of our approach for tailoring N -body interactions in DQS is the integration of the protocol into a variational feedback-loop scheme for quantum gate design. In contrast to traditional approaches to quantum compiling [70,[75][76][77][78], which optimize the fidelity with respect to some target unitary, here one exploits Hamiltonian learning to directly minimize the Hamiltonian distance to design a desired target Hamiltonian.
We note that the sample complexity of our protocol, i.e., the number of required experimental runs, can be reduced significantly with efficient parametrization schemes of Hamiltonians or Liouvillians. In particular, one can devise adaptive schemes similar to the one of Sec. IV B, where the ansatz parametrization is improved on the fly as new measurements are collected. Other promising strategies for reducing the sample complexity involve the preparation of (optimized) entangled states as inputs to the Hamiltonian or Liouvillian learning routine (see for instance Ref. [79]).
Finally, we emphasize that while the examples considered in the present work have focused on DQS of quantum many-body models motivated by condensed matter physics, the technique finds immediate application to the simulation of quantum chemistry problems [80], or the real time dynamics of lattice gauge theories [24,68,81,82]. In this appendix we explain how to determine the overall scale of a Hamiltonian that has been reconstructed using FHL as described in Sec. III. As discussed in the main text, the protocol of Ref. [34] can determine the Hamiltonian parameters only up to an undetermined multiplicative factor α(τ ), since c rec (τ ) and α(τ ) c rec (τ ) are equivalent solutions of the minimization problem (8). The overall-scale reconstruction derives from the HL protocol proposed in Ref. [52], which determines the Hamiltonian from measurements of time-resolved local observables. Likewise, here we determine α(τ ) based on the fact that, for any observableÂ it holds that for any time nτ , where · nτ = ψ(nτ )| · |ψ(nτ ) with |ψ(nτ ) = (Û exp τ ) n |ψ(0) for any chosen |ψ(0) . In the above equation, there is an inherent discretization error associated with the integral, which comes from the fact that expectation values are measured at integer multiples of τ . These errors can be suppressed by the choice of the integration routine, and become negligible compared to measurement noise and Trotter errors for small τ . For a reconstruction up to order L the Floquet Hamiltonian with c rec j (τ ) being the parameters learned from the protocol in Sec. III, and α(τ ) the prefactor we want to reconstruct. Using Eq. (B1) we can obtain α as follows We illustrate an application of this method in Fig. 8, where we plot the reconstructed α(τ ) as a function of τ for the first order learning shown in Fig. 2. In this example we used a fourth order integration routine. As a result, the discretization errors in evaluating the integral in Eq. (B3) are negligible for all values of τ up to approximately the Trotterization threshold. Thus, with this method we can accurately reconstruct α for each τ < τ * (the shaded vertical region in the figure). We expect this behavior to be rather generic: denoting with J the dominant energy scale inĤ exp F (τ ), in general we would have τ * ≤ J −1 . Since J −1 is the fastest timescale at which the expectation values −i[Â,ĥ j ] t change, we expect the τ discretization for τ < τ * ≤ J −1 , together with a high-order integration routine, to be sufficient to suppress integration errors. In this appendix we prove the scaling relation Eq. (11), that forms the basis of our order-by-order learning procedure, in the context of unitary dynamics. Its extension to dissipative DQS is proven in App. G.
Consider a DQS that implements a Trotter blockÛ exp τ . Based on an ansatz A and a set of initial states {|ψ m } m one obtains a constraint matrix M A and reconstructŝ H rec (τ ). In the following, we prove the following scaling relation where λ 1 is the lowest singular value of M A . Here O(τ k ) = {f (τ ) = ∞ l=k a l τ l |a l ∈ C}. Note that for any functions f (τ ), g(τ ) it holds that At the end of this section we comment on how Eq. (C1) relates to Eq. (11) in the main text. During the proof we will drop the superscript 'exp' inÛ exp τ to improve readability of the equations.
To prove Eq. (C1) we first note that because λ 1 = M A c . Moreover we find that where |ψ m the initial state that defines the constraint. For simplicity we set n = 1 in the following. We proceed by showing that for any operatorX one has which, together with Eq. (C2) and (C3), will prove the scaling relation (C1) for the commutator. We do this in two steps.
Step 1: Given an operator X we show that . The unitaryÛ τ is analytic in τ for all τ and hence each element can be written (Û τ ) i,j = l u i,j,l τ l for u i,j,l ∈ C. BecauseÛ τ acts linearly it follows that Step 2: Given an operator X we show that Direction "⇒" follows because taking the expectation value is a linear operation that cannot introduce lower orders in τ .
To prove direction "⇐" we splitX =X 0 +X τ witĥ X 0 = k−1 j=0x j τ j a polynomial in τ andX τ ∈ O(τ k ). From the assumed scaling ∀m : ψ m |X|ψ m ∈ O(τ k ) it follows ∀m : ψ m |X 0 |ψ m ∈ O(τ k ). Moreover, because taking the expectation value is a linear operation, it cannot increase the scaling beyond the maximal power inX 0 , which is τ k−1 and hence ∀m : ψ m |X 0 |ψ m = 0. Then ifX 0 = 0 we can choose an eigenstateX 0 |u = u|u with u = 0. This state fulfills u|X 0 |u = u = 0 and choosing sufficient constraints such that ∃m : |u = |ψm proves thatX 0 = 0 via contradiction. This means that X =X τ ∈ O(τ k ). Note that the proof assumes that |u is not an eigenstate ofÛ τ , since such a state would result in a row of zeros in the constraint matrix, and hence does not constitute a valid constraint. We note here that the choice of optimal constraints is a difficult task in general. However, in all the many-body examples presented in the main text, product states, obviously different from eigenstates ofÛ τ , were always a good choice as constraints.
In the main text we refer to the scaling relation Eq. (C1) in the context of the operator distanceĤ rec (τ )− H exp F (τ ) between the reconstructed HamiltonianĤ rec (τ ) and the experimental Floquet HamiltonianĤ exp F (τ ) (see Eq. (11)). This simplification applies in the case, where the ansatz does not contain any conserved quantity of any truncation of the Magnus expansion ofĤ exp F (τ ) (except the truncation ofĤ exp F (τ ) itself). (Note that if these conserved quantities are known, they can be 'projected out' from the ansatz.) In that case, using [Ĥ rec (τ ), Direction "⇐" in Eq. (C5) follows because the commutation withÛ exp τ , being a linear operation, cannot intro- In this appendix, we provide further analysis of the Trotterization threshold which can be detected via FHL. Additionally, we show how λ 1 [Eq. (10)] can be used as an indicator of the onset of quantum chaos in Trotterized dynamics [39][40][41]45]. For large values of τ (in units of the Hamiltonian timescales) no approximately local effective Hamiltonian can be reconstructed, which is signaled by λ 1 saturating to an approximately constant value in τ after the scaling regime (see, e.g., Fig. 2(a)). This corresponds to the transition to quantum chaos, where the properties of the operatorÛ τ can be approximately described by random matrix theory (RMT) [44,[46][47][48]. This notion of chaos can be probed through λ 1 , extracted from FHL, thanks to the following observation. When the dynamics are governed by unitaries from a given random matrix ensemble, and when the chosen initial states for the HL protocol consist of product states, we can analytically estimate the form of the squared constraint matrix Q, defined as and thus the value of λ 1 . Given these estimates, denoted with Q RMT and λ RMT , respectively, we can compute the differences between them and the measured Q and λ 1 : a decrease of these differences after the threshold τ * constitutes an indicator of the onset of chaos. Here, · for a matrix denotes its Frobenius norm. For the details on the analytical derivations of Q RMT and λ RMT we refer the reader to our recent work [45] in the context of collective spin systems. The basic idea is to calculate the elements of Q RMT as the expectation value of the elements of Q over the relevant matrix ensemble, which is the CUE in the case of the XXZ/Heisenberg DQS in the main text, i.e., where M i,j = tr(ρ iĥj ) − tr(ρ iÛ †ĥ jÛ ) withρ i denoting the chosen initial states for the protocol. This can be done using the techniques developed in [84]. Then, λ RMT can be estimated as the square root of the smallest eigenvalue of Q RMT .
Using these results, we compute the relative deviations (D2) and (D3) for the XXZ/Heisenberg DQS described by the Trotter cycle in Eq. (15). The results are shown in Fig. 9, for a zeroth order ansatz A 0 . Since in general one needs to perform a number n of Trotter cycles which is large enough in order forÛ n τ to resemble a random unitary [40], we first study how ∆ RMT λ and ∆ RMT Q change as a function of the number of cycles n -or the total simulation time-for fixed τ , by repeating the learning procedure described in Sec. III with increasing value of final time. This is shown in Fig. 9(a)-(b). For large values of τ (beyond the Trotterization threshold) these deviations decay to a small steady value, which is consistently smaller than the respective value for small τ . Then, we compute the steady value of ∆ RMT λ as a function of τ , expecting it to be large for small τ and close to zero in the chaotic regime. This is confirmed by the results shown in Fig. 9(c), where we see a sharp decrease of ∆ RMT λ in correspondence of the Trotterization threshold. This analysis further corroborates the intuitive connection between the Trotterization threshold in DQS and the onset of quantum chaos in the properties ofÛ τ , identifying in λ 1 a measurable indicator of this transition. In this appendix, we derive the estimates for the cost functions λ 1 (for Hamiltonian learning) and ∆ (for Liouvillian learning) in presence of measurement noise. We start with the derivation of the noise threshold in Eq (16) for Hamiltonian learning. We denote the measured constraint matrix with M = M + E, with M being the exact constraint matrix (in absence of measurement noise), E and additive error matrix and ∼ n −1/2 s . In the limit of large number of shots n s , we can approximate the error matrix elements E i,j as being independent and identically distributed normal random variables, i.e., E i,j ∼ N (0, 1). In the following, we calculate the expectation value of the smallest singular value of M , that is λ 1 ( M ), over the noise realizations of E. To this end, we consider the N A × N A squared constraint matrix where Q = M † M and From this definition, we have that where ε 1 ( Q) is the smallest eigenvalue of Q. We calculate ε 1 ( Q) using standard perturbation theory, taking V ( ) as a perturbation on the exact squared constraint matrix Q. We work under the assumption that λ 1 (M ) 1, and hence ε 1 (Q) 1, which corresponds to the case of a complete ansatz, or an ansatz missing only small terms, e.g., the higher order terms inĤ exp F (τ ) for small τ . Then, denoting with v n the eigenvectors of Q with eigenvalues ε n (Q), using standard perturbation theory up to second order in one has where (v|w) denotes the scalar product between two vectors. Using Eq. (E2), the identities and keeping terms up to order 2 , we obtain Using Jensen's inequality, we finally have which is the equation in the main text. We now show how this perturbative calculation can be extended to the case of Liouvillian learning explained in Sec. V. As before, we denote with G = G + G E the measured constraint matrix, with G = (G H G D ), and with b = b + e the measured constraint vector, with e an error vector with elements being i.i.d. normal random variables. Here ≈ n −1/2 s while G can be estimated by bounding the variance of the integrals (29)-(30) defining the matrix elements of G as with n being the number of Trotter steps used to evaluate the integral. The factor 4 comes from the integrands in Eq. (30) being commutators of Pauli operators. Assuming G † G being non-singular, we can write the solution of the noisy linear system Gc = b as where c ex = (G † G) −1 G † b is the exact coefficient vector.
Assuming that a perturbative expansion in the small parameter is valid for δc, i.e., δc = ∞ m=1 m δ m , we can insert it in the r.h.s. of Eq. (E11) obtaining The noisy residual vector ∆ = Gc rec − b becomes thus where R = I − G (G † G) −1 G † = α|λα(G)=0 |u α )(u α |, with |u α ) being the right singular vectors of G. Now we can calculate the expectation values over the random variables e and E obtaining E e,E∼N (0,1) ∆ ≤ (N con − N A )(1 + 4nτ 2 c ex 2 ) . (E15) Note that this estimate depends on c ex , which is in principle unknown. Denoting with J the typical energy scales in the Liouvillian, c ex 2 ∼ J 2 N , hence the unknown term in the above equation is suppressed for τ J −1 nN . When this is not the case, the above estimate is useful only if one has access to an estimate for c ex . main text. Here, the role of λ 1 (τ ) as 'learning error' is taken by ∆(τ ) = G(τ ) c rec (τ ) − b(τ ) , with G(τ ) = (G H (τ ) G D (τ )) and c rec (τ ) = (c H rec (τ ) c D rec (τ )) defined in Sec. V. In order to show the relation ∆(τ ) ∈ O(τ K ) ⇔ L exp F (τ ) − L rec (τ ) ∈ O(τ K ) , (G1) we rewrite the components of the residuals vector ∆(τ ) = G(τ ) c rec (τ ) − b(τ ) using the superoperator formalism. Each constraint operatorÂ k or density matrixρ(t) is represented by a vector |A k or |ρ(t) in the Liouville space H ⊗ H, and the Liouvillian L exp F (τ ) becomes a linear (super)operator acting on them. Denoting with A k |ρ(t) = tr Â † kρ (t) the inner product in the Liouville space, we can write which holds for all constraints A k |, initial states |ρ(0) and times nτ . We can therefore write which directly relates the components of the residuals vector with the (super)operator difference between the exact and the reconstructed Liouvillian. To show that (G1) holds, it is sufficient to prove that for anyÂ k and ρ(t) being analytic functions of τ one has Direction "⇒" is straightforward. Given that V(τ ) ∈ O(τ K ) element-wise, and for |A k and |ρ(t) analytic functions of τ , taking the expectation value A k |V(τ )|ρ(t) is a linear operation that cannot lower the power K. Let us now focus on direction "⇐". We start by noticing that the constraint operators |A k and the density operators |ρ(t) can be chosen as elements of an (hermitian) basis of the Liouville space. An example of such basis of density operators is given by the tensor product of bases of single-qubit density matrices, for instance asρ l ∈  H). Hence, since A k |V(τ )|ρ(t) ∈ O(τ K ) for any |A k and |ρ(t) , we have that ∀ k, l , V(τ ) k,l = ρ k |V(τ )|ρ l ∈ O(τ K ) , meaning that V(τ ) ∈ O(τ K ) element-wise.
The observation that V(τ ) ∈ O(τ K ) ⇔ L exp F (τ ) − L rec (τ ) ∈ O(τ K ), with · denoting the Frobenius matrix norm, concludes the proof of (G1). As a final remark, we note that the integrand A k |L exp F (τ )|ρ(t) in Eq. (G2) must be measured at integer multiples of τ , which becomes the minimal timestep that can be used for evaluating the integral. Hence, the evaluation of the integral in Eq. (G2) comes with an inherent integration error due to the Trotterization, which scales as O(τ m ) with an exponent m controlled by the integration routine. Hence, in order to observe scaling ∆(τ ) ∈ O(τ K ) due to Trotter errors, one has to choose an integration routine with m ≥ K.
tion of the zeroth order terms in the Floquet Hamiltonian of the DQS of a disordered Heisenberg/XXZ model, for several values of N and a fixed value of τ (below the Trotter threshold). To access larger system sizes, we used matrix product states (MPS) techniques implemented via the ITensor [85] library. The Trotterized time-evolution is very naturally implemented on MPS using the timeevolving block decimation (TEBD) algorithm. In Fig. 10 we show our results, where for each value of N we added simulated measurements noise from a finite measurement budget N meas. chosen to achieve a fixed parameter distance ĉ rec −ĉ ex ≈ 0.1, with the hat denoting a normalized coefficient vector.
As it is visible from Fig. 10(a), the total number N meas. of measurement shots required to reach a fixed accuracy in the reconstructed Hamiltonian scales polynomially with the system size N (note the log-log scale). The exponent of this power-law scaling is governed by the locality of the underlying Hamiltonian ansatz, and by the scaling of the gap of the constraint matrix with N . More specifically, in the present case we observe a scaling N meas. ∼ N 2ξ+1 , where the '+1' contribution comes from the fact that the number of ansatz terms to be measured scales linearly with N , while the 2ξ contribution comes from the spectral gap of M , i.e., ∆λ = λ 2 − λ 1 , as shown in the blue data of Fig. 10(b). There, we observe that the gap decreases as ∆λ ∼ N −ξ for the present example. To understand how the gap influences the number of measurements, we note that a perturbative result for the parameter distance ĉ rec −ĉ ex in presence of measurement noise parametrized by ≈ 1/ √ n s can be easily derived as ĉ rec −ĉ ex ≈ ∆λ [34,52], which implies that the number of shots n s to achieve a fixed parameter distance must scale as ∆λ 2 . Hence the observed scaling N meas. ∼ N 2ξ+1 . Finally, we note that the polynomial decay of ∆λ with the system size is a consequence of our choice of the model, which is here a disordered XXZ chain exhibiting MBL. For systems obeying the eigenstate thermalization hypothesis (ETH), the existence of a gap in the thermodynamic limit was proven in Ref. [34]. However, as already mentioned, the polynomial decay of ∆λ only worsen the scaling of the measurement budget for HL, without invalidating its scalability.