Daniel Stilck França
email: daniel.stilck_franca@ens-lyon.fr

Liubov A Markovich

V V Dobrovitski

Albert H Werner

Johannes Borregaard

Efficient and robust estimation of many-qubit Hamiltonians

Characterizing the interactions and dynamics of quantum mechanical systems is an essential task in developing quantum technologies. We propose an efficient protocol based on the estimation of the time-derivatives of few qubit observables using polynomial interpolation for characterizing the underlying Hamiltonian dynamics and Markovian noise of a multi-qubit device. For finite range dynamics, our protocol exponentially relaxes the necessary time-resolution of the measurements and quadratically reduces the overall sample complexity compared to previous approaches. Furthermore, we show that our protocol can characterize the dynamics of systems with algebraically decaying interactions. The implementation of the protocol requires only the preparation of product states and single-qubit measurements. Furthermore, we improve a shadow tomography method for quantum channels that is of independent interest and discuss the robustness of the protocol to various errors. This protocol can be used to parallelize the learning of the Hamiltonian, rendering it applicable for the characterization of both current and future quantum devices.

I. INTRODUCTION

Large quantum devices consisting of tens to hundreds of qubits have been realized across various hardware architectures [1][2][3][4] representing a significant step towards the realization of quantum computers and simulators with the potential to solve outstanding problems intractable for classical computers [5,6]. However, continued progress towards this goal requires careful characterization of the underlying Hamiltonians and dissipative dynamics of the hardware to mitigate errors and engineer the desired dynamics. The exponential growth of the dimension of the state space of a quantum device with the number of qubits renders this an outstanding challenge broadly referred to as the Hamiltonian learning problem .

To tackle this challenge, previous approaches make strong assumptions such as the existence of a trusted quantum simulator capable of simulating the unknown Hamiltonian [20,21] or the capability of preparing particular states of the Hamiltonian such as steady states and Gibbs states [23,25,[START_REF] Srinivasan | Sample-efficient learning of interacting quantum systems[END_REF][START_REF] Haah | Optimal learning of quantum Hamiltonians from high-temperature Gibbs states[END_REF][START_REF] Rudinger | Compressed sensing for hamiltonian reconstruction[END_REF][START_REF] Rouze | Learning quantum many-body systems from a few copies[END_REF], which may be difficult for realistic devices subject to various decoherence mechanisms.

Alternatively, several works [START_REF] Da Silva | Practical characterization of quantum devices without tomography[END_REF][START_REF] Bairey | Learning the dynamics of open quantum systems from their steady states[END_REF][START_REF] Zubida | Optimal short-time measurements for hamiltonian learning[END_REF] are built on the observation that a Master equation describes the evolution of any system governed by Markovian dynamics. Through this, one obtains a simple linear relation between time derivatives of expectation values and the parameters of the Hamiltonian, jump operators and decay the time derivatives are estimated using finite difference methods. Obtaining a good precision thus requires high time resolution, which is experimentally challenging given the finite operation time of gates and measurements. To estimate a Lindbladian parameter up to an additive error ϵ, the system has to be probed at times O(ϵ) apart and expectation values of observables have to be estimated up to a precision of O(ϵ 2 ), which translates to an overall O(ϵ -4 ) sample complexity to estimate each parameter.

In this article, we propose a protocol that alleviates these daunting experimental requirements. Our protocol requires only a time resolution of O(polylog(ϵ -1 )) representing an exponential improvement compared to previous protocols and gives an overall sample complexity to recover all parameters of a k-local n qubit Lindbladian up to precision ϵ of O(9 k ϵ -2 polylog(n, ϵ -1 )). We obtain this by estimating time derivatives using multiple temporal sampling points and robust polynomial interpolation [START_REF] Kane | Robust Polynomial Regression up to the Information Theoretic Limit[END_REF]. Furthermore, we show how to use shadow process tomography methods to estimate multiple parameters in parallel. In particular, we improve the results of Refs. [START_REF] Kunjummen | Shadow process tomography of quantum channels[END_REF][START_REF] Levy | Classical Shadows for Quantum Process Tomography on Near-term Quantum Computers[END_REF] in extending the framework of classical shadows to processes and Pauli matrices with an alternative proof, a result that is of independent interest. We also extend our analysis to long-range (algebraically decaying) interactions in the systems, obtaining the first results for such systems to the best of our knowledge. The necessary operations for our protocol are measurements in the Pauli bases on time-evolved product states consisting of Pauli eigenstates. These minimal requirements make our protocol feasible for characterization of both current and future quantum devices.

II. RESULTS

In order to use our protocol for an efficient characterization of a quantum device, two assumptions should be fulfilled:

1. The quantum device implements an (unknown)

Markovian quantum evolution on n qubits described by a time-independent Lindbladian, L.

2. We assume knowledge of the general structure of the interaction graph of the device i.e. which qubits are coupled to each other. Importantly, no assumptions are made regarding the couplings' exact form.

The first assumption ensures that the evolution of a general observable, O(t) is described by the Master equation, i.e. d dt O(t) = L(O(t)). We note that the Lindbladian captures both the Hamiltonian evolution and the dissipative dynamics of the device.

The second assumption bounds the size of the estimation task. It corresponds to making some assumptions about the locality of the generator, as evolutions without some locality assumptions have exponentially many parameters. However, having prior knowledge that, e.g. nearest neighbor couplings dominate in the device, makes the estimation task tractable. For now, we will assume that we know which qubits interact. Later, we will show that a bound on the support of each interaction and some technical assumptions on the evolution suffice to also learn the interaction graph from data.

Using the knowledge of the interaction graph, one can expand the Lindbladian in an operator basis, {P i } constructed from tensor products of single-qubit Pauli matrices and the identity [START_REF] Wolf | Quantum channels & operations: Guided tour[END_REF]:

L(O) = a i [P i , ρ] + D i,j (P † i OP j - 1 2 {P † j P i , O}) (1)
Such an expansion is always possible since this basis amounts to a Hilbert-Schmidt orthogonal set of Hermitian operators spanning the entire vector space. Note that the coefficient matrix D i,j needs to be positive for it to form a valid Lindbladian. Estimating the set of expansion coefficients {a i , D i,j } gives an estimation of L and thus a full characterization of the system. It is well known that the Master equation for the time derivative of the expectation value of a local observable O at time t = 0 for a given initial state ρ of the system gives us a linear equation for the expansion coefficients [START_REF] Da Silva | Practical characterization of quantum devices without tomography[END_REF][START_REF] Bairey | Learning the dynamics of open quantum systems from their steady states[END_REF][START_REF] Zubida | Optimal short-time measurements for hamiltonian learning[END_REF]. We use this to estimate the expansion coefficients going through three stages of classical pre-processing, quantum processing, and classical post-processing (see Fig. I).

Classical pre-processing: After expanding L in an operator basis, the following steps are performed.

1. Find a suitable complete set, {(ρ j , O i )} of multiqubit product states (ρ j ) and observables (O i ) for which the Master equation involves only a few selected expansion parameters of the Lindbladian for each element of the set. The set is complete in the sense that all expansion coefficients can be found by solving the Master equations for all elements in the set. As we show below, such a set can readily be found by considering initial states where only a few qubits are initialized as different eigenstates of the Pauli matrices while the remaining qubits are prepared in the maximally mixed state I/2. Since the initial states and the observables are products, this can be done efficiently.

Calculate the expectation values appearing on the right-hand side of the Master equations

Quantum processing: In order to solve for the expansion coefficients, we also need the values of the timederivatives appearing on the left-hand side of the Master equations, i.e. d dt tr [ρ j O i ]. These can be estimated using the quantum device. The naive approach is the following:

1. The quantum device is prepared in initial state ρ j and evolved for a time t k ∈ {t 0 , t 1 , . . . , t max } after which the observable O i is measured.

2. The above procedure is repeated for each element in the set {(ρ j , O i )} for all evolution times t k to obtain estimates of O i (t k ) j = tr [ρ j (t k )O i ] where ρ j (t k ) is the state of the system having evolved for time t k from the initial state ρ j . We note that the single qubit mixed states can be simulated by sampling eigenstates of the Pauli matrices at random.

For this naive approach, the sample complexity increases linearly with the size of the set {(ρ j , O i )} since the expectation values O i (t k ) j are estimated sequentially. However, we also propose a variation of the classical shadows protocol of Ref. [START_REF] Huang | Predicting many properties of a quantum system from very few measurements[END_REF] for process tomography that can reduce this to a logarithmic scaling. In essence, we can obtain estimates of all elements in the set {(ρ j , O i )} in parallel. This is done by the following steps:

1. Every qubit is prepared in a random Pauli matrix eigenstate and the system is evolved for a time t k after which each qubit is measured in a random single-qubit Pauli basis.

2. The above procedure is repeated O(3 w1+w2 ϵ -2 log(K)) times set by the required precision, ϵ, of the estimates, the size, K, of the set {(ρ j , O i )} and the weights (i.e. maximum number of sites differing from identity) of ρ j (w 1 ) and O i (w 2 ). The whole procedure is repeated for all times t k ∈ {t 0 , t 1 , . . . , t max }.

From the measurement statistics of the above procedure, it is possible to obtain accurate estimates of all O i (t k ) j and thus a parallel estimation is possible. We refer to Sec. V and the Supplementary Note for more details and proof of the method.

Whether to use the sequential approach or the parallel approach depends on the number of qubits and the weight of the states and observables in the set {(ρ j , O i )}. For few qubit processors, the sequential protocol may require fewest samples, however for local Hamiltonians on a lattice, the logarithmic scaling in system size of the parallel method will quickly be advantageous for larger processors. Importantly, both methods only require the preparation of single qubit Pauli states and measurements in single-qubit Pauli bases.

Classical post-processing: The final part of the characterization involves estimating d dt tr [ρ j O i (t)] from the experimentally obtained time trace of O i (t k ) j and solving for the expansion coefficients {a α }. This involves 1. Fit the time trace of O i (t k ) j with a lowdegree polynomial in the time, p i,j (t) and estimate d dt tr [ρ j O i ] as d dt p i,j (t)| t=0 . This is done for each element in the set {(ρ j , O i )}.

2. Solve the set of linear equations from the Master equations d dt tr [ρ j O i ] = tr [ρ j L(O i )] with respect to the expansion coefficients. This is possible since d dt tr [ρ j O i ] has been estimated from the polynomial fits and all expectation values appearing in tr [ρ j L(O i )] have been calculated leaving the expansion coefficients as the only unknown parameters. Following the steps above, a complete characterization of the underlying Hamiltonian and dissipative dynamics of the quantum device as given by the Lindbladian is obtained. The two key steps of the protocol are the choice of the set {(ρ j , O i )} and the polynomial interpolation used to obtain estimates of the time derivatives. Below, we outline the details of both steps and provide rigorous guarantees on the precision of the protocol. Importantly, we show that Lieb-Robinson bounds on the spread of correlations in the system can be used to ensure robust polynomial fitting of the time traces of expectation values allowing for an exponential relaxation of the temporal resolution compared to finite difference methods rendering the protocol feasible for near-term quantum devices.

Choosing the set of initial states and observables:

The first step in the classical pre-processing is to expand L in an operator basis constructed from tensor products of single-qubit Pauli matrices and the identity. The right-hand side (rhs) of the Master equation

d dt tr [ρ j O i (t)] = tr [ρ j L(O i )
] can be expanded as a sum of single Pauli matrices and their products. Our goal is to isolate the unknown expansion coefficients. To this end, we consider an initial state of the form

ρ (i,j) k,l = (I + σ (i) k ) 2 ⊗ (I + σ (j) l ) 2 ⊗ ρ n-2 , ( 2
)
where k, l = x, y, z. Thus, the i'th and j'th qubit are prepared in eigenstates of the Pauli matrices σ k and σ l while the state of the remaining n -2 qubits, ρ n-2 , is assumed to be the maximally mixed state. For a state of the form in Eq. ( 2) the rhs of the Master equation (see above) can be simplified greatly depending on the choice of the observable O. This is due to the properties of the Pauli matrices namely that they have vanishing trace and that

σ k σ l = δ kl I + iε klp σ p , ( 3
)
where δ kl is the Kroenecker delta function and ε klp is the Levi-Civita symbol. As we show explicitly in the Supplementary Note, if a single qubit Pauli observable,

O = σ (i)
l , is chosen, then only the single qubit terms of the rhs of the Master equation involving the i'th qubit will have non-vanishing trace and, using the relation in Eq. ( 3), the different single qubit Pauli expansion coefficients (the coefficients of terms in the expansion that only involves single qubit Pauli matrices) can be isolated.

After isolating the single qubit expansion coefficients, the coefficients related to two-qubit Pauli terms (σ

(i) i ⊗ σ (j) j ) can be isolated by choosing observables of the form O = σ (i) i ⊗ σ (j) j
in a similar manner. This procedure can be iterated to isolate higher order expansion coefficients by considering observables involving more and and more qubits.

In the Supplementary Note, we provide a detailed derivation of how all expansion parameters can be isolated for a general Hamiltonian with terms coupling between two and k qubits and arbitrary single qubit dissipation terms. We note that already for two qubit dissipation terms, deriving linear combinations of initial states and expectation values that allow us to isolate different parameters is quite cumbersome and we do not do this explicitly. However, from a numerical point of view, this is a trivial task. Indeed, as remarked before, each pair of Pauli strings gives us access to a linear equation for the different evolution parameters.

After collecting enough equations to ensure that the linear system is invertible, the precision with which we need to estimate each expectation value to ensure a reliable estimation of the parameters is controlled by the condition number of the matrix describing the system of linear equations. As both estimating the condition number and solving the linear system can be done efficiently, we conclude that estimating dissipative terms acting on a constant number of qubits does not pose a significant challenge from a numerical perspective.

Robust polynomial interpolation: As described above, a key step in our learning algorithm is to obtain information about the time-derivatives of observables at t = 0. For this, we rely on robust polynomial interpolation. Accordingly, based on expectation values O i (t k ) j for a set of times t k we want to extract a polynomial p i,j (t) such that we can estimate d dt tr [ρ j O i ] as d dt p i,j (t)| t=0 . For this approach to work, we have to be able to control the degree of the polynomial p i,j (t) in order to give an upper bound on the number of sampling points t k for which we will have to determine O i (t k ) j experimentally. In the following, we briefly outline how such a guarantee on the degree of p i,j (t) can be obtained and refer to the methods section for a more detailed proof.

Our argument proceeds in two steps: establishing that local expectation values are well approximated by lowdegree polynomials and then showing how to robustly extract the derivative of the polynomial from this information. Before we give an overview of the ideas, we will also need to introduce some notation to deal with Lindbladians acting on different parts of the system. We will denote by Λ the whole system of qubits and for some subset B ⊂ Λ of qubits, we denote by L B , the generator truncated to those qubits.

The first step of our proof is to establish that the expectation value O(t) of a local observable O that evolves under a Lindbladian L B restricted to some sub-region B up to some time t max , can indeed be approximated up to error ε by a degree-d polynomial, where d depends linearly on the size of B, t max and log(ε -1 ).

For the second step of our argument, it remains to show under which circumstances, we can restrict the evolution of the Pauli-strings P m a , that we identified in the previous step, to a local generator. The main insight here is that for finite range (or sufficiently quickly fast decaying) interactions, the dynamics of any local observable O exhibits an effective light cone quantified by a Lieb-Robinson bound (LRB) [START_REF] Lieb | The finite group velocity of quantum spin systems[END_REF][START_REF] Poulin | Lieb-Robinson Bound and Locality for General Markovian Quantum Dynamics[END_REF][START_REF] Kliesch | Lieb-Robinson Bounds and the Simulation of Time-Evolution of Local Observables in Lattice Systems[END_REF][START_REF] Hastings | Locality in Quantum Systems[END_REF][START_REF] Kliesch | Lieb-robinson bounds and the simulation of time-evolution of local observables in lattice systems[END_REF][START_REF] Kuwahara | Strictly Linear Light Cones in Long-Range Interacting Systems of Arbitrary Dimensions[END_REF]. The LR-bound, in turn, allows us to restrict the Lindbladian on the full system to a generator coupling only systems in the vicinity of the support of O, where the size of this shielding region only grows linearly with t max .

Bringing these two arguments together, we can first employ the LR-bound to restrict the dynamics to a sub-region around the support of the Pauli-string, P a m , and then approximate the corresponding evolution on that finite region up to error ε by a polynomial of degree O poly(t max , log(ϵ -1 )) . Now, making use of the techniques from Ref. [START_REF] Kane | Robust Polynomial Regression up to the Information Theoretic Limit[END_REF], we can extract the first derivative of this polynomial from measurements at O polylog(ϵ -1 ) different times t k . Indeed, in that work, the authors show how to perform polynomial interpolation reliably only given approximations of the values and even under the presence of outliers.

Robustness of the protocol to experimental errors:

In the previous sections, we assumed that we can prepare Pauli eigenstates and measure in Pauli eigenbases and did not consider state preparation and measurement (SPAM) errors. Furthermore, we assumed that the dynamics is described exactly be a generator that is local. Such an idealized scenario rarely comes up in practice and it is important to develop protocols that also work when these conditions are only met approximately. Here, we will provide the main arguments showing that our protocol is indeed robust to such imperfections. We refer to Supplementary Note for the more detailed technical derivations and statements.

First, we consider the setting where the SPAM errors are well-characterized and independent of the basis we prepare and measure. In this case, it is possible to adapt the protocol to incorporate this information without changing its performance significantly. For instance, if the SPAM errors are described by local depolarizing noise with depolarizing probability p, we can still recover 2-local Hamiltonians with single qubit noise with a O((1 -p) -4 ) sampling overhead when compared to the noiseless case. In particular, this only depends on the local noise rate. We then show that this is the case more generally: if the SPAM errors are not characterized well or are highly dependent on the particular state we are preparing or the basis we are measuring, then once again the effect of SPAM errors will be independent of the system's size and only depend on the precision we wish to obtain.

Besides SPAM errors, our protocol is also robust to perturbations of the generators, even nonlocal ones.

More precisely, assume that the true generator is of the form L = L ′ + ∆, where L is a Lindbladian that satisfies a LR-bound and ∆ is an arbitrary, potentially global perturbation. A naive bound would imply that for short times, the expectation value of this perturbed evolution could be up to t ∆ ∞→∞ away from the unperturbed value, where ∆ ∞→∞ = sup ∥X∥∞≤1 ∆(X) ∞ . However, we show that only the effect of ∆ on few-qubit observables contributes to the bound. A precise statement is given in Thm. 9.1. of the Supplementary Note 9. As an illustrative example, assume that ∆ is a small all-toall coupling between all pairs of qubits of order τ , e.g.

∆(X) = τ n i,j=1,i̸ =j i[σ (i) x σ (j) x , X]. ( 4
)
Clearly, such an evolution does not satisfy a LR-bound and ∆ ∞→∞ = O(τ n 2 ). However, we show that under such a perturbation, the expectation value is only perturbed by O(tτ n), which corresponds to the local effect of this global perturbation. To show this, we once again resort to LR bounds and the fact that we only need to measure local observables. Taken together, these results show that our protocol is robust to both SPAM noise and deviations from the assumptions we impose on the generators.

Learning the structure of the interactions: So far, we have also assumed that we have knowledge of the interaction graph. However, an astute reader might have remarked that we do not require this explicit knowledge for our protocol to work: indeed, we just use it to decide which parameters to estimate and need to restrict the possible interaction graphs to have an LR bound for the evolution. Thus, a brute-force approach to learning the interaction structure also follows from our results if we wish to estimate all k-body couplings of strength at least η: we just estimate all k body-couplings up to a precision η/4 and discard all of those that we see are smaller than 3 4 η. This can be done with a number of samples that scales like O(9 k η -2 polylog(η -1 ) log(n)).

III. NUMERICAL EXAMPLES

To investigate the performance of our protocol for experimentally relevant parameters, we performed numerical simulations of a multi-qubit superconducting device. We consider a system with tunable couplers similar to the Google Sycamore chip [1]. This design relies on a cancellation of the next-next-nearest coupling between two qubits through the direct coupling with a coupler [START_REF] Yan | Tunable coupling scheme for implementing high-fidelity two-qubit gates[END_REF][START_REF] Sung | Realization of high-fidelity cz and zzfree iswap gates with a tunable coupler[END_REF]. We consider a generic system consisting of a 2D grid of qubits with exchange coupling between nearest neighbors. The dynamics are described through a Lindblad equation with the effective two-qubit Hamiltonian for each neighboring qubit pair (i, j) [START_REF] Yan | Tunable coupling scheme for implementing high-fidelity two-qubit gates[END_REF][START_REF] Sung | Realization of high-fidelity cz and zzfree iswap gates with a tunable coupler[END_REF]

H ij = k=i,j 1 2 ωk σ (k) z + g i g j ∆ ij + g ij (σ (i) + σ (j) -+σ (i) -σ (j) + ) (5)
for i = j = 1, . . . , n and a dissipation term acting on the i'th qubit and having jump operators σ

(i) -, σ (i)
+ (generalised amplitude damping) and σ (k) z (pure dephasing).

Here ωk = ω k + g 2 k ∆ k is the Lamb-shifted qubit frequency, g i is the coupling between the i'th qubit and the coupler, and g ij is the direct two-qubit coupling. We have assumed that ∆ k = ω k -ω c < 0 where ω c (ω k ) is the frequency of the coupler (k'th qubit) and have defined

1/∆ ij = (1/∆ i + 1/∆ j )/2
. By adjusting the frequencies of the coupler and the qubits, the effective qubitqubit interaction can be canceled up to experimental precision. Typical qubit frequencies are around 5 -6 GHz [1]), while ∆ k ∼ -1 GHz, g ij ∼ 10 -20 MHz, and g j ∼ 100 MHz [START_REF] Yan | Tunable coupling scheme for implementing high-fidelity two-qubit gates[END_REF][START_REF] Sung | Realization of high-fidelity cz and zzfree iswap gates with a tunable coupler[END_REF]. In our simulation, we assume that all qubit frequencies and couplings have been characterized up to a precision of 100 kHz using standard characterization techniques [1] and consequently, that all couplers have been tuned off with the same precision i.e. gigj ∆ij + g ij ∼ 100 kHz. Considering a layout of 16 qubits (see below for the interaction graph), we randomly sample all qubit frequencies and qubit-qubit interactions according to Gaussian distributions with zero mean and standard deviation of 100 kHz.

In addition to the Hamiltonian evolution, we also include dissipative dynamics in our numerical simulation. We include quasi-static random frequency shifts of the qubits leading to effective dephasing with a characteristic timescale of T * 2 ∼ 150 µs as well as pure dephasing resulting in a transverse relaxation on a timescale T 2 ∼ 60 µs representing state of the art coherence times [1,[START_REF] Sung | Realization of high-fidelity cz and zzfree iswap gates with a tunable coupler[END_REF]. Finally, we include longitudinal relaxation of the qubits through an amplitude damping channel on the time scale of T 1 ∼ 60 µs. We refer to Supplementary Note 2for a more detailed discussion and Table III for the sampled parameters of our simulation.

In Fig. 2, we plot the average estimation error as a function of the temporal resolution set by the value of the initial time step, t 0 . For this plot, we only included the Hamiltonian evolution in the numerical simulation together with quasi-static random frequency shifts of the qubits. This was to lower the run time of the simulation allowing us to investigate the performance for a broad range of initial times. We assumed the total run time of the experiment was fixed such that t 0 × S is constant, where S is the number of samples. From the figure, we clearly see the improved scaling of our protocol of the estimation error with the time-step size compared to using a finite difference method [START_REF] Da Silva | Practical characterization of quantum devices without tomography[END_REF][START_REF] Bairey | Learning the dynamics of open quantum systems from their steady states[END_REF][START_REF] Zubida | Optimal short-time measurements for hamiltonian learning[END_REF]. Besides already performing better at the time resolution for moderate values of the initial time, we see that after a threshold initial time around 10 -0.7 , the performance is not limited by the initial time, only the shot noise. In contrast, the finite difference method still requires smaller initial times to improve on the error with the same shot noise.

We also investigated the robustness of our method with respect to shot-noise for a fixed time resolution. For these simulations, we again only included the Hamiltonian evolution together with quasi-static random frequency shifts of the qubits to have a practical run time of the simulation. From Fig. 3 we see that for a fixed time resolution of 30 ns our protocol results in an average estimation error that improves linearly with the shot-noise down to an error below 10 -4 . This is in contrast to finite difference methods, where the estimation error plateaus around 10 -3 since it becomes limited by the time resolution. This is a clear effect of the exponential improvement of our protocol w.r.t. the time resolution compared to finite difference methods.

Finally, we performed a numerical simulation that included the pure dephasing and amplitude damping noise as described above and estimated the σ X σ X couplings between the qubits. As shown in Fig. 4, we obtain reliable estimates of all 22 couplings demonstrating how our method allows the estimation of specific terms in the Hamiltonian despite the dynamics being governed by the full (dissipative) Lindbladian. For simplicity, we did not explicitly estimate the single qubit Hamiltonian parameters and the Lindbladian decay rates.

For all estimations above, we fitted to degrees 1 -7 and picked the one with the smallest average error on the sampled points. In Supplementary Note 6 we give explicit rigorous bounds on how to pick the parameters like the degree for a given desired precision, but we believe that a heuristic approach like the one pursued in the numerics performs well in practice: as long as the interpolating polynomial approximates well the observed points and new points we did not fit to, the degree should be adequate. We note that, although the robust interpolation methods of Ref. [START_REF] Kane | Robust Polynomial Regression up to the Information Theoretic Limit[END_REF], in principle, require random times, we performed numerical experiments with deterministic times on systems with 16 qubits.

IV. DISCUSSION

In conclusion, we have proposed a Hamiltonian learning protocol based on robust polynomial interpolation that has rigorous guarantees on the estimation error. Our protocol offers an exponential reduction in the required temporal resolution of the measurements compared to previous methods and a quadratic reduction in the overall sampling complexity for finite-range interactions. Our protocol only requires the preparation of single qubit states and single qubit measurements in the Pauli bases and is robust to various imperfections such as SPAM errors and Hamiltonian perturbations. This makes it suitable for the characterization of both near-term and future quantum devices.

Furthermore, the recovery of multiple parameters can be highly parallelized by resorting to a variation of classical shadows to quantum channels we improve here.

Our method allows for the characterization of a general local Markovian evolution consisting of a unitary Hamiltonian part and a dissipative part. While we have only explicitly considered single-qubit dissipation here, we believe that our protocol is also valid for general multi-qubit dissipation as outlined above but leave the explicit analysis of this to future work. We have also analysed the performance of our protocol for algebraically decaying interactions which we believe to be the first results for Hamiltonian learning of such systems. The convergence of our method can be ensured for interactions decaying faster than the dimension of the system. We note, how- ever, that improved bounds on the locality of such systems might improve this result in the future.

V. METHODS

Here we detail and formalize our results regarding the estimation error guarantees of our protocol. In particular, we detail the use of Lieb-Robinson bounds on the spread of correlations in the system to bound the error. Furthermore, we outline the shadow tomography method for the parallelization of the measurements.

Derivative estimation: Define f (t) = tr e tL (O)ρ and L B to be the Lindbladian truncated to a subregion B of the interaction graph. Our protocol consists of first estimating f (t i ) up to an error O(ϵ) for random times t 1 , . . . , t m . The curve of f (t) is then fitted to a lowdegree polynomial p, and p ′ (0) is taken as an estimate for the derivative f ′ (0) = tr [L(O)ρ]. Below we prove the accuracy and robustness of this method. The first step is Theorem 1, which establishes under what conditions f (t) is indeed well-approximated by a low-degree polynomial.

Theorem 1. Let L be a local Lindbladian on a Ddimensional lattice Λ. Moreover, let t max , ϵ > 0 and O Y be a 2-qubit observable, supported on some region Y ∈ Λ, such that O Y ≤ 1, holds. Then there is a polynomial p of degree d = O poly(t max , log(ϵ -1 )) , (6) such that for all 0 ≤ t ≤ t max :

tr e tLΛ (O Y )ρ -p(t) ≤ ϵ, (7)
and

p ′ (0) = tr [L Λ (O Y )ρ], holds.
The main technical tool required for the proof are Lieb-Robinson bounds (LRB) [START_REF] Lieb | The finite group velocity of quantum spin systems[END_REF][START_REF] Poulin | Lieb-Robinson Bound and Locality for General Markovian Quantum Dynamics[END_REF][START_REF] Kliesch | Lieb-Robinson Bounds and the Simulation of Time-Evolution of Local Observables in Lattice Systems[END_REF][START_REF] Hastings | Locality in Quantum Systems[END_REF][START_REF] Kliesch | Lieb-robinson bounds and the simulation of time-evolution of local observables in lattice systems[END_REF][START_REF] Kuwahara | Strictly Linear Light Cones in Long-Range Interacting Systems of Arbitrary Dimensions[END_REF], which ascertain that the dynamics of local observables under a time evolution with a local Lindbladian have an effective lightcone. More precisely, we need that for regions Y ⊂ B we have

(e tL B -e tLΛ )(O Y ) ≤ c 1 exp(-µ dist(Y, Λ\{B}))(e vt -1), (8)
to hold for constants c 1 , µ and v, where dist() is the distance in the graph.

From the LRB we can show that the dynamics is wellapproximated by a low-degree polynomial. We leave the details of the proof to the Supplementary Note 3 and only discuss the main steps here. The general idea of going from the LRB to the low-degree polynomial is to truncate the Taylor series of the evolution under L B for B large enough and take that as the approximating polynomial. As the derivatives of the evolution under L B only scale with the size of the region B, this allows us to show that the Taylor series converges quickly. To simplify the presentation we did not give explicit numerical constants in the bounds, but in Sec.A of the Supplementary Note 7 we explicitly determine the constants for the polynomial approximation.

Now that we have concluded that the expectation value is well-approximated by a small degree polynomial, we continue to show that we can reliably extract the derivative from approximations of the expectation values for different t. This is formally stated in the following theorem.

Theorem 2. Let L be a Lindbladian on a D-dimensional regular lattice. Suppose we can measure the expectation value of two-body Pauli observables on Pauli eigenstates in the time interval [t 0 , t max ] under L for t 0 as

t -1 0 = O polylog(ϵ -1 ) (9)
and t max = 2 + t 0 . Then, measuring the expectation values at

m = O polylog(ϵ -1 ) (10)
random times up to precision O(ϵ/ polylog(ϵ -1 )), is sufficient to obtain an estimate of the Lindbladian coefficients âi of a i satisfying

|â i -a i | = ϵ. ( 11
)
This yields a total sample complexity of S = O ϵ -2 polylog(ϵ -1 ) .

Importantly, Theorem 2 bypasses both requiring small initial times and O(ϵ -4 ) sample complexities.

To go from Thm. 1 to Thm. 2 we first need to establish that we can robustly infer an approximation of p from finite measurement data subject to shot noise. Subsequently, we need to show that it will also allow us to reliably estimate p ′ (0). Let us start with approximating p.

a. Robust polynomial interpolation: We will resort to the robust polynomial interpolation methods of [START_REF] Kane | Robust Polynomial Regression up to the Information Theoretic Limit[END_REF] to show Thm. 2. We review their methods in more detail in Supplementary Note 6. In our setting, the randomly sampled point x i , i ∈ 1, m correspond to different times ∈ [t 0 , t max ] and the y i ∈ 4 to approximations of the expectation value of the evolution at that time. Furthermore, the y i satisfy the promise that there exists a polynomial p of degree d and some σ > 0 such, that

y i = p(x i ) + w i , |w i | ≤ σ, ( 12
)
hold, for strictly more than half of the y i . The rest might be outliers. In our setting, the magnitude of σ corresponds to the amount of shot noise present in the estimates of the expectation values. The authors of [START_REF] Kane | Robust Polynomial Regression up to the Information Theoretic Limit[END_REF] then show that by sampling m = O(d log(d)) points from the Chebyshev measure on [t 0 , t max ], a combination of ℓ 1 and ℓ ∞ regression allows us to find a polynomial p of degree d that satisfies:

max x∈[t0,tmax] |p(x) -p(x)| = O(σ). ( 13
)
Although the details of the ℓ 1 and ℓ ∞ interpolation are more involved and described in the Supplementary Note 4, a rough simplification of the procedure is the following. First, we find a polynomial p 1 of degree d that minimizes

i |p 1 (x i ) -y i |. After finding p 1 we compute the poly- nomial p ∞ that minimizes max i |p ∞ (x i ) -(y i -p 1 (x 1 ))|.
We then output p = p 1 + p ∞ as our guess polynomial. Note that finding both p 1 and p ∞ can be cast as linear programs and thus can be solved efficiently [START_REF] Vandenberghe | Convex Optimization[END_REF].

By combining this result with Thm. 1, we robustly extract a polynomial that approximates the curve t → tr e tL (O Y )ρ up to O(ϵ) for t ∈ [t 0 , t max ]. Indeed, we only need to estimate the expectation value f (t i ) up to ϵ for enough t i and run the polynomial interpolation.

Note that Eq. ( 13) only allows us to conclude that p -p is small. However, we are ultimately interested in the curve's derivative at t = 0, as the derivative contains information about the parameters of the evolution. For arbitrary smooth functions, two functions being close on an interval does not imply that their derivatives are close as well. Fortunately, for polynomials the picture is simpler. By the Thm. 1 one has to estimate the first derivative of a polynomial at t = 0 but not of the actual function. A classical result from approximation theory, Markov brother's inequality [START_REF] Markoff | über Polynome, die in einem gegebenen Intervalle möglichst wenig von Null ab-weichen[END_REF], allows us to quantify the deviation of the derivatives given a bound on the degree and a bound like Eq. (13). Putting these observations together, we arrive at Thm. 2. The details of the proof are given in the Supplementary Note 4.

Generalizations of Thm. 2:

We also generalize Thm. 2 in two directions. First, we extend the results to interactions acting on k qubits instead of 2. As long as the noise is constrained to act on 1 qubit and k = O(1), this generalization is straightforward. Indeed, we only need to measure an observable that has the same support as the Pauli string and does not commute with it, as it is then always possible to find a product initial state that isolates the parameter. Generalizing to noise acting on more than one qubit makes it more difficult to isolate the parameters of the evolution as described in the main text. In that case, it then becomes necessary to solve a system of linear equations that couples different parameters. Although our method still applies, analysing this scenario would require picking the observables and initial states in a way that the system of equations is well-conditioned and we will not discuss this case in detail here.

Second, another important generalization is to go beyond short-range systems. Although we have only stated our results for short-range systems in Thm. 2, our techniques apply to certain long-range systems. As this generalization is more technical, we leave the details to the Supplementary Note 4 and constrain ourselves to discussing how the statement of Thm. 2 changes for more general interactions.

Only one aspect of the previous discussion changes significantly for long-range interactions: how the r.h.s. of Eq. ( 8) generalizes. More precisely, let us assume that for some injective function h : 4 → 4 with h(x) = o(1), we have

(e tL B -e tLΛ )(O Y ) ≤ h(dist(Y, V \{B}))(e vt -1). (14)
For instance, for short-range or exponentially decaying interactions, h will be an exponentially decaying function. Then we can restate Thm. 2 in terms of h -1 . As we show in Thm. 6.1. of the Supplementary Note 6, for a precision parameter ϵ > 0 and evolution on a Ddimensional lattice, assume that we pick the initial time as

t 0 = O   h -1 ϵ 2(e 2.5v -1) D log(ϵ -1 ) -2   . (15)
Furthermore, assume that we estimate the expectation value of local observables up to precision O(ϵ) at

Õ h -1 ϵ 2(e 2.5v -1) D log(ϵ -1
)) points. Then we can estimate each parameter up to an error of

O   ϵ h -1 ϵ 2(e 2.5v -1) D log(ϵ -1 ) 2   , (16)
through, the same procedure as in the local case. Note that the error in Eq. ( 16) only tends to 0 as ϵ → 0 if

h -1 ϵ 2(e 2.5v -1) D log(ϵ -1 ) = o(ϵ -1
), holds, i.e. the function h must decay fast enough. In the Supplementary Note 5, we discuss examples of systems with algebraically decaying interactions for which this is satisfied. For instance, for potentials that decay like r -α with α > 5D-1

we obtain that h -1 (ϵ) = O(ϵ -1 α-3D
), holds. We summarize the resulting resources in in the Supplementary Table 6.

But the message of bounds like ( 16) is that it is still possible to obtain bounds on the error independent of the system's size beyond short-range systems. However, this comes at the expense of requiring higher precision and sampling from more points.

Another important observation is that the assumption that we know the structure of the interactions exactly is not required. Indeed, our method is robust to Hamiltonian perturbations of the model as long as the resulting evolution still satisfies a LR bound. For instance, suppose that there actually is a non-negligible interaction between qubits i and j that is not accounted by our model. As long as the resulting time evolution still satisfies a LR bound, our results still hold. As the linear equation to isolate any parameter is independent of that parameter, we can still apply our techniques in this setting.

Parallelizing the measurements: To parallelize the measurement procedure and ensure that we can obtain experimental data to estimate all parameters simultaneously, we resort to a classical shadow process tomography method. Although some papers in the literature already discussed classical shadows for process tomography [START_REF] Kunjummen | Shadow process tomography of quantum channels[END_REF][START_REF] Levy | Classical Shadows for Quantum Process Tomography on Near-term Quantum Computers[END_REF], we present a simplified and streamlined proof that also gives an improved sample complexity for the observables relevant to this work in the Supplementary Note.

More precisely, we show that given a quantum channel Φ, Pauli strings P 1 a , . . . , P K1 a that differ from the identity on at most ω a sites and Pauli strings P 1 b , . . . P K2 b that differ from the identity on at most ω b sites, it is possible to obtain estimates êm,l of 2 -n tr P m a Φ(P l b ) satisfying

|2 -n tr P m a Φ(P l b ) -êm,l | ≤ ϵ ( 17
)
for all m, l with probability at least 1 -δ from

O(3 ωa+ω b ϵ -2 log(K 1 K 2 δ -1 )) (18)
samples. More precisely, the protocol of shadow process tomography requires preparing Eq. ( 18) many different random initial product Pauli eigenstates and measuring them in random Pauli bases. This makes it feasible to implement it in the near-term. We discuss it in more detail in the Supplementary Note 7, as this protocol may be of interest beyond the problem at hand. The shadow process tomography protocol is ideally suited for our Hamiltonian learning protocol. Indeed, note that to learn k-body interactions, we only required the preparation of initial states ρ l that differ from the maximally mixed state on k qubits and measure Pauli strings P m supported on at most k qubits. Furthermore, for a system of n qubits in total, there are at most 16 k n k ≤ 16 k n k such states or Pauli strings. We conclude that we can estimate all required expectation values for a given time step using

O(9 k ϵ -2 k log(nδ -1 )) (19)
samples. As our protocol requires estimating expectation values at a total of polylog(ϵ -1 ) time steps, we can gather the data required to recover all the O(n) parameters of the evolution from O(ϵ -2 polylog(n, ϵ -1 )) samples through the shadow process tomography protocol whenever k = O(1).

VI. DATA AVAILABILITY

The data generated and analyzed during the current study are available from the authors upon request.

VII. CODE AVAILABILITY

The computer codes developed and tested by the authors, and the input files used for producing the presented data, are available upon request. The current version of the codes is not designed for broad distribution and requires substantial hardware resources for reproducing the presented results. Correct installation of the current version of the codes, preparation of correct input files, and correct analysis of the results require substantial expertise, and may need separate instructions in each specific case.

IX. AUTHOR CONTRIBUTION STATEMENTS

DSF, AHW and JB conceived the original project. DSF, AHW, and LM conducted the analytical derivations with input from VVD and JB. The software for the many-qubit simulations was written and tested by VVD. LM and DSF performed the numerical analysis with input from AHW and JB. All authors contributed to the discussion of the results and writing of the manuscript.

X. COMPETING INTERESTS STATEMENT

The authors declare no competing interests.

XI. FIGURE CAPTIONS

Fig. 1: Sketch of the proposed protocol to estimate an unknown Lindbladian, L, of a multi-qubit device. In the first step of classical pre-processing, the interaction graph between qubits is identified from the physical connectivity of the device. Then the unknown Lindbladian is written in a general form using an operator basis of Pauli strings, {P i } and a suitable set of initial states and observables, {(ρ j , O i )} is chosen. In the second step of quantum processing, a time trace (expectation value) of each element of the set is obtained from the preparation and evolution of single qubit Pauli eigenstates on the quantum device followed by measurements in single qubit Pauli basis. In the last step of classical post-processing, each time trace is fitted to a low-degree polynomial to estimate the derivative of the observable. From these, an estimate of the Lindbladian, L est , is obtained from the Master equation.

Fig. 2: The median quality of recovery of of one 2qubit coupling using interpolation methods and those based on numerical derivatives [START_REF] Da Silva | Practical characterization of quantum devices without tomography[END_REF][START_REF] Bairey | Learning the dynamics of open quantum systems from their steady states[END_REF][START_REF] Zubida | Optimal short-time measurements for hamiltonian learning[END_REF] as a function of the initial time. We assumed that the total time of the experiment is fixed. That is, we let the initial time times the total number of samples for each time step to be a constant (10 7 for this plot). We used the expectation value of 40 equally spaced times which had the same distance between each other as the initial time. For each initial time, we simulated 1000 instances of the recovery protocol, always adding shot noise with the same standard deviation to the data. The dots correspond to the median quality of recovery, whereas the lower and upper end correspond to the 25 and 75 percentile. We ran the simulation on a system with 16 qubits.

Fig. 3: Median quality of recovery of one 2-qubit coupling using interpolation methods and those based on numerical derivatives [START_REF] Da Silva | Practical characterization of quantum devices without tomography[END_REF][START_REF] Bairey | Learning the dynamics of open quantum systems from their steady states[END_REF][START_REF] Zubida | Optimal short-time measurements for hamiltonian learning[END_REF] as a function of the standard deviation of the shot noise. The initial time for this estimate is 30 ns and here we also generated 1000 instances of the noise with a given standard deviation.

The plot shows the median quality of the recovery and the 25 and 75 percentiles. We see that the quality of the recovery for the interpolation decays approximately linearly with the shot noise, before plateauing at shot noise -6. For the numerical derivative, we see two regimes: first a linear decay of the error until a shot of noise of order 10 -3 . After that, the error plateaus and does not improve even with smaller shot noise. This is because for numerical derivative methods, at this point the dominant error source comes from the choice of initial time, whereas for interpolation it is at -6. Importantly, we see that interpolation consistently provides estimates that are no worse than the numerical derivatives method.

Fig. 4: Error in the recovery of σ X σ X couplings of a quantum system with a geometry similar to the Sycamore processor using numerical derivatives and interpolation. Note that while we only plot the estimation of the Hamiltonian couplings, the numerical simulation included the full Lindbladian including both dephasing due to quasistatic random frequency shifts of the qubits, pure dephasing and amplitude damping noise. The initial time for each coupling was 0.1 µs in the simulation. Note that interpolation consistently outperforms numerical derivatives, sometimes by several orders of magnitude. We chose the time steps and the number of samples to compare both methods as in Fig. 2.

Supplementary Note 1: Selecting States and Observables to Isolate Parameters

Any Hamiltonian can be written as

H = m1 α1 a (m1) α1 σ (m1) α1 + m1,m2 α1α2 a (m1,m2) α1α2 σ (m1) α1 σ (n) α2 + . . . ( 1
)
≡ m1 α1 a (m1) α1 H (m1) α1 + m1,m2 α1α2 a (m1,m2) α1α2 H (m1,m2) α1α2 + . . . ,
where Roman indices identify the subspace on which the operator acts, and Greek indices identify the Pauli operator, e.g. α = x. No assumption about the dimension or structure of the hermitian Hamiltonian is needed for this expansion to be valid. For a Markovian noise environment, the evolution of a quantum system ρ 0 is described by a Master equation of the form

dρ t dt t=0 = -i[H, ρ 0 ] + n m=1 3 µ,ν=1 L (m) µ,ν (σ (m) µ ρ 0 σ (m) ν † - 1 2 {σ (m) ν † σ (m) µ , ρ 0 }), ( 2
)
where H is the Hamiltonian describing the evolution of the system, L (m) µ,ν are the elements of the Lindblad matrix, expressed in an operator basis consisting of the different combinations of single-qubit Pauli matrices {σ}. Multiplying it from the right hand side on an observable O and taking the trace, we can write

d dt tr [ρ t O] | t=0 = -i m1 α1 a (m1) α1 tr [H (m1) α1 , ρ 0 ]O -i m1m2 α1α2 a (m1m2) α1α2 tr [H (m1m2) α1α2 , ρ 0 ]O -. . . . . . -i m1...m k α1...α k a (m1...m k ) α1...α k tr [H (m1...m k ) α1...α k , ρ 0 ]O . . . (3)
+ m µ,ν L (m) µν tr (σ (m) µ ρ 0 σ (m) ν † - 1 2 {σ (m) ν † σ (m) µ , ρ 0 })O .
Let us introduce the notation

B (m1,...,m k ) α1,...,α k (ρ 0 , O) ≡ -i tr [H (m1...m k ) α1...α k , ρ 0 ]O . ( 4
)
To isolate the coefficients {a

(m1) α1 , a (m1m2) α1α2
, . . . } we observe the n qubit state with the following density matrix

ρ (i,j) τi,τj = ρ (i) τi ⊗ ρ (j) τj ⊗ ρ n-2 , τ = {1, 2, 3} (5)
where the i and j are the Pauli qubits, namely ρ

(i,j) τi,j = (I + σ (i,j)
τi,j )/2, and ρ n-2 is the density matrix of all other qubits, which we set to be maximally mixed. Using σ α σ β = δ αβ I + iε αβγ σ γ , we can find the following relations

B (i) α1 (ρ (i,j) τi,τj , O) = 1 2 ε α1τiγ tr (σ (i) γ ⊗ I)O + tr (σ (i) γ ⊗ σ (j) τj )O , ( 6
)
B (j) α1 (ρ (i,j) τi,τj , O) = 1 2 ε α1τj ζ tr (I ⊗ σ (j) ζ )O + tr (σ (i) τi ⊗ σ (j) ζ )O , ( 7
)
B (ij) α1α2 (ρ (i,j) τi,τj , O) = 1 2 ε α1τiγ δ α2τj tr (σ (i) γ ⊗ I)O + tr (σ (i) γ ⊗ σ (j) α2 )O (8) + ε α2τj η δ α1τi tr (I ⊗ σ (j) η )O + tr (σ (i) α1 ⊗ σ (j) η )O ,
where O is acting on (i, j) qubits. Selecting O = σ

(i) ξi ⊗ σ (j)
ξj , we can rewrite the latter matrix element as follows

B (i) α1 (ρ (i,j) τi,τj , σ (i) ξi ⊗ σ (j) ξj ) = 2ε α1τiγ δ ξiγ δ ξj τj =      2ε α1τiγ δ ητj if ξ i = γ, ξ j = η, 2ε α1τiγ if ξ i = γ, ξ j = η, η = τ j , 0 else , ( 9
)
B (j) α1 (ρ (i,j) τi,τj , σ (i) ξi ⊗ σ (j) ξj ) = 2ε α1τj ζ δ ξj ζ δ ξiτi =      2ε α1τj ζ δ κτi if ξ j = ζ, ξ i = κ, 2ε α1τj ζ if ξ j = ζ, ξ i = κ, κ = τ i , 0 else , B (ij) α1α2 (ρ (i,j) τi,τj , σ (i) ξi ⊗ σ (j) ξj ) = 2 ε α1τiγ δ ξiγ δ ξj α2 + ε α2τj η δ ξiα1 δ ξj η =      2 ε α1τiγ δ ηα2 + ε α2τj η δ γα1 if ξ i = γ, ξ j = η, 2ε α1τiγ δ τj α2 if ξ i = γ, ξ j = η, η = τ j , 0 else , where η, γ ∈ {x, y, z}. Selecting other observable O = σ (i)
ξi , we can rewrite ( 6)-( 8) differently as

B (i) α1 (ρ (i,j) τi,τj , σ (i) ξi ) = ε α1τiγ δ ξiγ = ε α1τiγ if ξ i = γ, 0 else , B (j) α1 (ρ (i,j) τi,τj , σ (i) ξi ) = 0, (10)
B (ij) α1α2 (ρ (i,j) τi,τj , σ (i) ξi ) = 2ε α1τiγ δ α2τj δ ξiγ = 2ε α1τiγ δ α2τj if ξ i = γ, 0 else . Next, selecting O = σ (j)
ξj , we can rewrite ( 6)- (8) as

B (i) α1 (ρ (i,j) τi,τj , O) = 0, B (j) α1 (ρ (i,j) τi,τj , σ (i) ξi ) = ε α1τj ζ δ ξj ζ = ε α1τj ζ if ξ j = ζ, 0 else , ( 11
)
B (ij) α1α2 (ρ (i,j) τi,τj , σ (j) ξj ) = 2ε α2τj η δ α1τi δ ξj η = 2ε α2τj η δ α1τi if ξ j = η, 0 else .
For k ≥ 2 we can write the general matrix element:

B (m1,...,m k ) α1,...,α k (ρ (i,j) τi,τj , O) = 1 2 ε αiτiγ tr (σ (1) α1 ⊗ • • • ⊗ σ (i) γ ⊗ • • • ⊗ σ (j) αj ⊗ • • • ⊗ σ (k) α k )O (12)
+ δ αj τj tr (σ (1)

α1 ⊗ • • • ⊗ σ (i) γ ⊗ • • • ⊗ I (j) ⊗ • • • ⊗ σ (k) α k )O + ε αj τj η tr (σ (1) α1 ⊗ • • • ⊗ σ (i) αi ⊗ • • • ⊗ σ (j) η ⊗ • • • ⊗ σ (k) α k )O + δ αiτi tr (σ (1) α1 ⊗ • • • ⊗ I (i) ⊗ • • • ⊗ σ (j) η ⊗ • • • ⊗ σ (k) α k )O , where O is acting on (1, . . . , k) qubits. Let O = σ (1) ξ1 ⊗ • • • ⊗ σ (k)
ξ k , holds. Then we can rewrite (12) as follows

B (1,...,k) α1,...,α k (ρ (i,j) τi,τj , σ (1) ξ1 ⊗ • • • ⊗ σ (k) ξ k ) ≡ -iT r [σ (1) α1 ⊗ • • • ⊗ σ (k) α k , ρ (i,j) τi,τj ](σ (1) ξ1 ⊗ • • • ⊗ σ (k) ξ k ) = 2 k-1 (ε αiτiγ δ α1ξ1 . . . δ γξi . . . δ αj ξj . . . δ α k ξ k + ε αj τj η δ α1ξ1 . . . δ αiξi . . . δ ηξj . . . δ α k ξ k ) (13) = 2 k-1 (ε αiτiγ δ ηαj + ε αj τj η δ γαi ) if ξ i = γ, ξ j = η, ξ 1,...,k = α 1,...,k , 0 else .
From this result, for k = 3 we get

B (i,j,l) αiαj α l (ρ (i,j) τi,τj , σ (i) ξi ⊗ σ (j) ξj ⊗ σ (l) ξ l ) = 2(ε αiτiγ δ γ,ξi δ αj ξj + ε αj τj η δ αiξi δ ηξj )δ α l ξ l . ( 14
)
To isolate {L (m) µν } we recall the assumption that ρ n-2 is the density matrix of the maximally mixed state. Then the Lindblad part of the equation ( 3) for an observable O = σ

(i) ξi ⊗ σ (j) ξj is L k µν (ρ (i,j) τi,τj , σ (i) ξi ⊗ σ (j) ξj ) ≡ L (k) µν tr σ (k) µ ρ (i,j) τi,τj σ (k) ν - 1 2 {σ (k) ν σ (k) µ , ρ (i,j) τi,τj } (σ (i) ξi ⊗ σ (j) ξj ) (15) =      L (i) µν (2iε µνγ δ γξi + 2δ ντi δ µξi -3 2 δ µν δ ξiτi )δ τj ξj if k = i L (j) µν (2iε µνγ δ γξj + 2δ ντj δ µξj -3 2 δ µν δ ξj τj )δ τiξi if k = j 0 else
.

Let us substitute the conditions ( 9) in (15). We get the following results:

L i µν (ρ (i,j) τi,τj , σ (i) ξi ⊗ σ (j) ξj ) = L (i) µν (2iε µνγ + 2δ ντi δ µγ - 3 2 δ µν δ γτi )δ τj η if ξ i = γ, ξ j = η, ( 16
)
and

L j µν (ρ (i,j) τi,τj , σ (i) ξi ⊗ σ (j) ξj ) = L (j) µν (2iε µνγ δ γη + 2δ ντj δ µη - 3 2 δ µν δ ητj )δ τiγ if ξ i = γ, ξ j = η. ( 17
)
However, for an observable O = σ

(i)
ξi the Lindblad part of the equation ( 3) is the following

L k µν (ρ (i,j) τi,τj , σ (i) ξi ) ≡ L (k) µν tr σ (k) µ ρ (i,j) τi,τj σ (k) ν - 1 2 {σ (k) ν σ (k) µ , ρ (i,j) τi,τj } (σ (i) ξi ) (18)
=      L (i) µν (2iε µνγ δ γξi + δ ντi δ µξi -3 2 δ µν δ ξiτi + δ µξi δ τiν ) if k = i 0 if k = j 0 else .
Substituting the conditions (10) in (18), we get

L i µν (ρ (i,j) τi,τj , σ (i) ξi ) = L (i) µν (2iε µνγ + 2δ ντi δ µγ - 3 2 δ µν δ γτi ), if ξ i = γ. ( 19
)
Next, for an observable O = σ

(i) ξi ⊗ σ (j) ξj ⊗ σ (l)
ξ l we can write

L k µν (ρ (i,j) τi,τj , σ (i) ξi ⊗ σ (j) ξj ⊗ σ (l) ξ l ) ≡ L (k) µν tr σ (k) µ ρ (i,j) τi,τj σ (k) ν - 1 2 {σ (k) ν σ (k) µ , ρ (i,j) τi,τj } (σ (i) ξi ⊗ σ (j) ξj ⊗ σ (l) ξ l ) = 2iL (l) µν ε µνγ δ γξ l δ τiξi δ τj ξj if k = l 0 else. . ( 20
)
Then, according to the results of the previous subsection, we get

L l µν (ρ (i,j) τi,τj , σ (i) ξi ⊗ σ (j) ξj ⊗ σ (l) ξ l ) = 2iL (l) µν ε µνγ δ γα l δ τiγ δ τj η if ξ i = γ, ξ j = η, ξ l = α l . ( 21
)
Finally, for an observable O = σ

(1) ξ1 ⊗ . . . σ (i) ξi • • • ⊗ σ (j) ξj • • • ⊗ σ (k)
ξ k , k > 3, the Lindblad part of the equation ( 3) is the following

L (m) µν tr σ (m) µ ρ (i,j) τi,τj σ (m) ν - 1 2 {σ (m) ν σ (m) µ , ρ (i,j) τi,τj } (σ (1) ξ1 ⊗ • • • ⊗ σ (i) ξi • • • ⊗ σ (j) ξj • • • ⊗ σ (k) ξ k ) = 0. (22)
A. Final Results

After we selected the different observable operators and defined the density matrix ρ 0 = ρ

(i) τi ⊗ρ (j) τj ⊗ I 2n-2 2 2n-2
, where the i and j qubits are in the Pauli states, we are ready to isolate the desired coefficients. For an observable O = σ

(i) γ ⊗ σ (j) η
we can write the equation (3) as follows

d dt tr ρ t (σ (i) γ ⊗ σ (j) η ) | t=0 = -2 α1 a (i) α1 ε α1τiγ δ ητj -2 α1 a (j) α1 ε α1τj η δ γτi (23) -2 α1,α2 a (ij) α1α2 ε α1τiγ δ ηα2 + ε α2τj η δ γα1 + µ,ν L (i) µν (2iε µνγ + 2δ ντi δ µγ - 3 2 δ µν δ γτi )δ τj η + L (j) µν (2iε µνγ δ γη + 2δ ντj δ µη - 3 2 δ µν δ ητj )δ τiγ .
Selecting τ j = η, τ i = γ, we can isolate the coefficients of the type a

(ij)
α1α2 in (23), namely

d dt tr ρ t (σ (i) γ ⊗ σ (j) η ) | t=0 = -2 α1,α2 a (ij) α1α2 ε α1τiγ δ ηα2 + ε α2τj η δ γα1 . ( 24
)
Let us call ρ tτiτj the density matrix evaluated by the Hamiltonian evolution from ρ 0 . From (24) we can find a

(ij)
α1α2 . To this end, we select γ = y, η = y and four pairs τ i , τ j ∈ {z, x; z, z; x, x; x, z} to get the system of equations

d dt tr ρ tzx (σ (i) y ⊗ σ (j) y ) | t=0 = 2(a (ij) xy -a (ij) yz ), d dt tr ρ tzz (σ (i) y ⊗ σ (j) y ) | t=0 = 2(a (ij) xy + a (ij) yx ), (25)
d dt tr ρ txx (σ (i) y ⊗ σ (j) y ) | t=0 = -2(a (ij) yz + a (ij) zy ), d dt tr ρ txz (σ (i) y ⊗ σ (j) y ) | t=0 = -2(a (ij) yx + a (ij) zy ).
Since a

(ij) yz = a (ij)
zy and a

(ij) xy = a (ij)
yx , we can write

a (ij) xy = 1 4 d dt tr ρ tzz (σ (i) y ⊗ σ (j) y ) | t=0 , a (ij) yz = - 1 4 d dt tr ρ txx (σ (i) y ⊗ σ (j) y ) | t=0 . ( 26
)
Selecting γ = x, η = y and four pairs τ i , τ j ∈ {y, z; y, x; z, z; z, x} to get the system of equations

d dt tr ρ tyz (σ (i) x ⊗ σ (j) y ) | t=0 = 2(a (ij) xx + a (ij) zy ), d dt tr ρ tyx (σ (i) x ⊗ σ (j) y ) | t=0 = 2(a (ij) zy -a (ij) xz ), (27)
d dt tr ρ tzz (σ (i) x ⊗ σ (j) y ) | t=0 = 2(a (ij) xx -a (ij) yy ), d dt tr ρ tzx (σ (i) x ⊗ σ (j) y ) | t=0 = -2(a (ij) yy + a (ij) xz ).
Hence

a (ij) xx = 1 2 d dt tr ρ tyz (σ (i) x ⊗ σ (j) y ) | t=0 + 1 4 d dt tr ρ txx (σ (i) y ⊗ σ (j) y ) | t=0 , ( 28
)
a (ij) xz = - 1 2 d dt tr ρ tyx (σ (i) x ⊗ σ (j) y ) | t=0 - 1 4
d dt tr ρ txx (σ (i) y ⊗ σ (j) y ) | t=0 , a (ij) yy = - 1 2 d dt tr ρ tzz (σ (i) x ⊗ σ (j) y ) | t=0 + 1 2 d dt tr ρ tyz (σ (i) x ⊗ σ (j) y ) | t=0 + 1 4
d dt tr ρ txx (σ (i) y ⊗ σ (j) y ) | t=0 .
Selecting γ = x, η = z and four pairs τ i , τ j ∈ {y, y; z, x; y, x; z, y}, we get

d dt tr ρ tyy (σ (i) x ⊗ σ (j) z ) | t=0 = -2(a (ij) xx -a (ij) zz ), d dt tr ρ tzx (σ (i) x ⊗ σ (j) z ) | t=0 = 2(a (ij) xy -a (ij) yz ), (29)
d dt tr ρ tyx (σ (i) x ⊗ σ (j) z ) | t=0 = -2(a (ij) zz -a (ij) xy ), d dt tr ρ tzy (σ (i) x ⊗ σ (j) z ) | t=0 = -2(a (ij) xx + a (ij) yz ).
Hence the last coefficient is

a (ij) zz = - 1 2 d dt tr ρ tyx (σ (i) x ⊗ σ (j) z ) | t=0 + 1 4 d dt tr ρ tzz (σ (i) y ⊗ σ (j) y ) | t=0 . ( 30
)
To find the other coefficients we select τ j = η, τ i = γ and rerewrite (23) as

d dt tr ρ t (σ (i) γ ⊗ σ (j) η ) | t=0 = -2 α1 a (i) α1 ε α1τiγ -2 α1,α2 a (ij) α1α2 ε α1τiγ δ ηα2 (31) + 2 µ=x,y,z L (i) µµ δ µτi δ µγ + 2 µ̸ =ν µ,ν=x,y,z L (i) µν (iε µνγ + δ ντi δ µγ ).
Next, for τ j = η, τ i = γ, we can rewrite (23) as

d dt tr ρ t (σ (i) γ ⊗ σ (j) η ) | t=0 = -2 α1 a (j) α1 ε α1τj η -2 α1,α2 a (ij) α1α2 ε α2τj η δ γα1 + 2 µ=x,y,z L (i) µµ δ µτj δ µη + 2 µ̸ =ν µ,ν=x,y,z L (i) µν (iε µνγ δ γη + δ ντj δ µη ). (32)
Selecting an observable O = σ

(i)
γ , we can write

d dt tr ρ t σ (i) γ | t=0 = - α1 a (i) α1 ε α1τiγ -2 α1,α2 a (ij) α1α2 ε α1τiγ δ α2τj ( 33
)
+ 2

µ=x,y,z

L (i) µµ (δ µτi δ µγ - 3 4 δ γτi ) + 2 µ̸ =ν µ,ν=x,y,z L (i) µν (iε µνγ + δ ντi δ µγ ).
For the other observable O = σ

(j)
η , the result is the following

d dt tr ρ t σ (j) η ) | t=0 = - α1 a (j) α1 ε α1τj η -2 α1,α2 a (ij) α1α2 ε α2τj η δ α1τi + 2 µ=x,y,z L (j) µµ (δ µτj δ µη - 3 4 δ ητj ) + 2 µ̸ =ν µ,ν=x,y,z L (j) µν (iε µνγ δ γη + δ ντj δ µη ). (34)
Substituting [START_REF] Gu | Practical black box hamiltonian learning[END_REF] in [START_REF] Bairey | Learning the dynamics of open quantum systems from their steady states[END_REF], we can write

d dt tr ρ t (σ (i) γ ⊗ σ (j) η ) | t=0 - d dt tr ρ t σ (i) γ | t=0 = - α1 a (i) α1 ε α1τiγ (35) + 2 α1,α2 a (ij) α1α2 ε α1τiγ (δ τj α2 -δ ηα2 ) + 3 2 µ=x,y,z L (i) µµ δ γτi .
Since in [START_REF] Bairey | Learning the dynamics of open quantum systems from their steady states[END_REF] the conditions τ j = η, τ i = γ, hold, we can rewrite the latter equation as

d dt tr ρ t (σ (i) γ ⊗ σ (j) η ) | t=0 - d dt tr ρ t σ (i) γ | t=0 = - α1 a (i) α1 ε α1τiγ . ( 36
)
Solving the latter equation, we find a (i) α1 . To this end, we select γ = y, τ i = x, τ j = η ∈ {x, y, z} and get

a (i) x = - d dt tr ρ txη (σ (i) y ⊗ σ (j) η ) | t=0 + d dt tr ρ txη σ (i) y | t=0 . ( 37
)
Next, for γ = x, τ i = z, τ j = η ∈ {x, y, z} we get the solution of (36), namely

a (i) y = - d dt tr ρ tzη (σ (i) x ⊗ σ (j) η ) | t=0 + d dt tr ρ tzη σ (i) x | t=0 . ( 38
)
Finally, for γ = x, τ i = y, τ j = η ∈ {x, y, z} the solution is

a (i) z = d dt tr ρ tyη (σ (i) x ⊗ σ (j) η ) | t=0 - d dt tr ρ tyη σ (i) x | t=0 . ( 39
)
Substituting [START_REF] Yu | Practical and efficient hamiltonian learning[END_REF] in [START_REF] Zubida | Optimal short-time measurements for hamiltonian learning[END_REF], we can write

d dt tr ρ t (σ (i) γ ⊗ σ (j) η ) | t=0 - d dt tr ρ t σ (j) η ) | t=0 = - α1 a (j) α1 ε α1τj η + 2 α1,α2 a (ij) α1α2 ε α2τj η (δ τiα1 -δ γα1 ) + 3 2 µ=x,y,z L (i) µµ δ ητj . ( 40
)
Since τ j = η, τ i = γ, hold, we can rewrite it as follows

d dt tr ρ t (σ (i) γ ⊗ σ (j) η ) | t=0 - d dt tr ρ t σ (j) η ) | t=0 = - α1 a (j) α1 ε α1τj η . ( 41
)
Solving the latter equation, we find a (j) α1 . Selecting η = y, τ j = z and τ i = γ ∈ {x, y, z}, we get

a (j) x = d dt tr ρ tγz (σ (i) γ ⊗ σ (j) y ) | t=0 - d dt tr ρ tγz σ (j) y ) | t=0 . ( 42
)
Selecting η = x, τ j = z and τ i = γ ∈ {x, y, z}, we get

a (j) y = - d dt tr ρ tγz (σ (i) γ ⊗ σ (j) x ) | t=0 + d dt tr ρ tγz σ (j) x ) | t=0 . ( 43
)
Finally, selecting η = y, τ j = x and τ i = γ ∈ {x, y, z}, we get the last coefficient of this type

a (j) z = - d dt tr ρ tγx (σ (i) γ ⊗ σ (j) y ) | t=0 + d dt tr ρ tγx σ (j) y ) | t=0 . ( 44
)
All the coefficients with the corresponding observables and initial states are given in Table I.

a (i) α (O, {τ (i) i , τ (j) j }) Equation a (i) x (σ (i) y ⊗ σ (j) η , {x, η}); (σ (i) y , {x, η}), η ∈ {x, y, z} (37) a
(i) y (σ (i) x ⊗ σ (j) η , {z, η}); (σ (i) x , {z, η}), η ∈ {x, y, z} (38) a
(i) z (σ (i) x ⊗ σ (j) η , {y, η}); (σ (i)
x , {y, η}), η ∈ {x, y, z} (39) a

(j) x (σ (i) γ ⊗ σ (j) y , {γ, z}); (σ (i) y , {γ, z}), γ ∈ {x, y, z} (42) a (j) y (σ (i) γ ⊗ σ (j)
x , {γ, z}); (σ

(i) x , {γ, z}), γ ∈ {x, y, z} (43) a (j) z (σ (i) γ ⊗ σ (j) y , {γ, x}); (σ (i) y , {γ, x}), γ ∈ {x, y, z} (44) a (ij) xx (σ (i) x ⊗ σ (j) y , {y, z}); (σ (i) y ⊗ σ (j) y , {x, x}) (28) a (ij) yy (σ (i) x ⊗ σ (j) y , {z, z}); (σ (i) x ⊗ σ (j) y , {y, z}); (σ (i) y ⊗ σ (j) y , {x, x}) (28) a (ij) zz (σ (i) x ⊗ σ (j) z , {y, x}); (σ (i) y ⊗ σ (j) y , {z, z}) (30) a (ij) xy (σ (i) y ⊗ σ (j) y , {z, z}) (26) a (ij) yz (σ (i) y ⊗ σ (j) y , {x, x}) (26) a (ij) xz (σ (i) x ⊗ σ (j) y , {y, x}); (σ (i) y ⊗ σ (j) y , {x, x}) (28)
Supplementary Table I. The first column represents the type of the estimated Hamiltonian parameters a

(i) α i , a (ij)
α i ,α j , αi, αj ∈ {x, y, z}. In the third column the number of equation for every parameter is provided, depending from the pairs of the observable O and the initial state ρ (i,j)

τ i ,τ j = ρ (i) τ i ⊗ ρ (j)
τ j ⊗ ρn-2, τi, τj = {x, y, z}, given in the second column. For example, to estimate a

(i) x
we use [START_REF] Rouze | Learning quantum many-body systems from a few copies[END_REF] with two pairs of observables and states: {O, ρ (i,j)

τ i ,τ j } = {(σ (i) y ⊗ σ (j) η , ρ (i,j) x,η ); (σ (i) y , ρ (i,j) x,η )}.

For an observable

O = σ (i) γ ⊗ σ (j) η ⊗ σ (l) α l we can write the equation d dt tr ρ t (σ (i) γ ⊗ σ (j) η ⊗ σ (l) α l ) | t=0 = -4 αiαj α l a (ijl) αiαj α l (ε αiτiγ δ ηαj + ε αj τj η δ γαi ) (45)
+ 2i µ,ν L (l) µν ε µνγ δ γα l δ τiγ δ τj η .
Selecting τ j = η, τ i = γ, we can rewrite [START_REF] Kliesch | Lieb-Robinson Bounds and the Simulation of Time-Evolution of Local Observables in Lattice Systems[END_REF] as

d dt tr ρ t (σ (i) γ ⊗ σ (j) η ⊗ σ (l) α l ) | t=0 = -4 αiαj α l a (ijl) αiαj α l ε αiτiγ δ ηαj . ( 46
)
From this equation we can find a (ijl)

αiαj α l . For an observable O = σ (1) α1 ⊗ . . . σ (i) γ • • • ⊗ σ (j) η • • • ⊗ σ (k) α k , k > 3 we can write d dt tr ρ t (σ (1) α1 ⊗ . . . σ (i) γ • • • ⊗ σ (j) η • • • ⊗ σ (k) α k ) | t=0 = -2 k-1 α1,...,α k a (1...k) α1...α k (ε αiτiγ δ ηαj + ε αj τj η δ γαi ). (47)
Selecting τ j = η, τ i = γ, we can rewrite [START_REF] Kliesch | Lieb-robinson bounds and the simulation of time-evolution of local observables in lattice systems[END_REF] as

d dt tr ρ t (σ (1) α1 ⊗ . . . σ (i) γ • • • ⊗ σ (j) η • • • ⊗ σ (k) α k ) | t=0 = -2 k-1 α1,...,α k a (1...k) α1...α k ε αiτiγ δ ηαj . ( 48
)
Solving this equation, we find a (1...k) α1...α k . From (15) we can find three Lindbladian coefficients. Selecting τ i = y, τ j = η, γ = x, we can find

L (i) xy = 1 2 d dt tr ρ tyη (σ (i) x ⊗ σ (j) η ) | t=0 -a (i) z -a (ij) zη . ( 49
)
Next, for τ i = z, τ j = η, γ = x, we deduce

L (i) xz = 1 2 d dt tr ρ tzη (σ (i) x ⊗ σ (j) η ) | t=0 + a (i) y + a (ij) yη . ( 50
)
Finally, for τ i = z, τ j = η, γ = y, the coefficient is

L (i) yz = 1 2 d dt tr ρ tzη (σ (i) y ⊗ σ (j) η ) | t=0 -a (i) x -a (ij) xη . ( 51
)
From (23) we can find three more coefficients. Selecting τ j = y, τ i = z, γ = z, η = x, we get

L (j) xy = 1 2 d dt tr ρ tzy (σ (i) z ⊗ σ (j) x ) | t=0 -a (j) z -a (ij) zz . ( 52
)
For τ j = y, τ i = y, γ = y, η = z, we can deduce

L (j) zy = 1 2 d dt tr ρ tyy (σ (i) y ⊗ σ (j) z ) | t=0 + a (j) x + a (ij) yx . ( 53
)
Selecting τ j = z, τ i = y, γ = y, η = x, we get

L (j) xz = 1 2 d dt tr ρ tyz (σ (i) y ⊗ σ (j) x ) | t=0 + a (j) y + a (ij) yy . ( 54
)
From [START_REF] Gu | Practical black box hamiltonian learning[END_REF] we find the following coefficients:

L (i) xx = 1 10 -3 d dt tr ρ txτj σ (i) x | t=0 -3 d dt tr ρ tyτj σ (i) y | t=0 + 2 d dt tr ρ tzτj σ (i) z | t=0 , ( 55
)
L (i) yy = 1 10 -3 d dt tr ρ txτj σ (i) x | t=0 + 2 d dt tr ρ tyτj σ (i) y | t=0 -3 d dt tr ρ tzτj σ (i) z | t=0 , L (i) zz = 1 10 2 d dt tr ρ txτj σ (i) x | t=0 -3 d dt tr ρ tyτj σ (i) y | t=0 -3 d dt tr ρ tzτj σ (i) z | t=0 .
From [START_REF] Yu | Practical and efficient hamiltonian learning[END_REF] the following coefficients can be found

L (j) xx = 1 10 -3 d dt tr ρ tγ1x σ (i) x | t=0 -3 d dt tr ρ tγ2y σ (i) y | t=0 + 2 d dt tr ρ tγ3z σ (i) z | t=0 , ( 56
)
L (j) yy = 1 10 -3 d dt tr ρ tγ1x σ (i) x | t=0 + 2 d dt tr ρ tγ2y σ (i) y | t=0 -3 d dt tr ρ tγ3z σ (i) z | t=0 , L (j) zz = 1 10 2 d dt tr ρ tγ1x σ (i) x | t=0 -3 d dt tr ρ tγ2y σ (i) y | t=0 -3 d dt tr ρ tγ3z σ (i) z | t=0 .
All the Lindbladian coefficients with the corresponding observables and initial states are given in Table II.

L (i) µν (O, {τ (i) i , τ (j) j }) Equation L (i) xx, L (i) yy , L (i) zz (σ (i)
x , {x, τj1}); (σ

(i) y , {y, τj2}); (σ (i) z , {z, τj3}), τj1,2,3 ∈ {x, y, z} (55) L (j) xx , L (j) yy , L (j) zz (σ (i) x , {γ1, x}); (σ (i) y , {γ2, y}); (σ (i) z , {γ3, z}), γ1 ∈ {y, z}, γ2 ∈ {x, z}, γ3 ∈ {x, y} (56) L (i) xy (σ (i) x ⊗ σ (j) η , {y, η}), η ∈ {x, y, z} (49) L (i) xz (σ (i) x ⊗ σ (j) η , {z, η}), η ∈ {x, y, z} (50) L (i) yz (σ (i) y ⊗ σ (j) η , {z, η}), η ∈ {x, y, z} (51) L (j) xy (σ (i) z ⊗ σ (j) x , {z, y}) (52) L (j) xz (σ (i) y ⊗ σ (j) x , {y, z}) (54) L (j) yz (σ (i) y ⊗ σ (j) z , {y, y}) (53)
Supplementary Table II. The first column represents the type of the estimated Lindbladian parameters

L (i)
µν , µ, ν ∈ {x, y, z}. In the third column the number of equation for every parameter is provided, depending from the pairs of the observable O and the initial state ρ (i,j)

τ i ,τ j = ρ (i) τ i ⊗ ρ (j)
τ j ⊗ ρn-2, τi, τj = {x, y, z}, given in the second column. For example, to estimate L (j) yz we use (53) with one pair of observable and state: {O, ρ (i,j)

τ i ,τ j } = {σ (i) y ⊗ σ (j) z , ρ (i,j) y,y }.

Supplementary Note 2: Numerical simulations

For simulations of the Hamiltonian learning protocol, we employed direct numerical solution of the time-dependent Schrodinger equation for the whole system of n = 16 qubits, whose time-dependent wavefunction |ψ(t) is represented as an array of 2 n complex numbers, normalized to 1. The evolution includes both unitary component, governed by the system's Hamiltonian, and three non-unitary components: the first one, stemming from the quasi-static random frequency shifts, leading to essentially non-Markovian dephasing of the qubits with the characteristic time T * 2 , the second component, described as a set of Lindblad superoperators corresponding to the phase damping channel, leading to Markovian transverse decoherence of the qubits on the timescale T 2 , and the third component, also leading to Markovian evolution of the qubit, and described as a set of Lindblad superoperators corresponding to the amplitude damping channel, which leads to longitudinal relaxation of the qubit on the timescale T 1 . For precise meaning of the terms "Markovian" and "non-Markovian", see explanations below in subsection 2 B.

A. Unitary evolution

The simulation of the unitary (Hamiltonian) evolution was performed using the 2nd order Suzuki-Trotter decomposition of the evolution operator. The total Hamiltonian of the system in question can be written as

H = n m1,m2=1 α1,α2=x,y J m1,m2 σ (m1) α1 σ (m2) α2 + 1 2 n m1=1 Ω m1 σ (m1) z , ( 57
)
For the problem considered in this paper the couplings J m1,m2 are restricted to the nearest neighbor qubits on a 2-D lattice. Note that the actual frequency Ω m1 of the m 1 -th qubit in the Hamiltonian (57) is different from its nominal frequency ωm mentioned in the main text; the reasons for this difference are explained in subsection 2 B below. The Hamiltonian (57) is represented as a sum

H = H X + H Y + H Z (58) H X = n m1,m2=1 J m1,m2 σ (m1) x σ (m2) x , H Y = n m1,m2=1 J m1,m2 σ (m1) y σ (m2) y , H Z = 1 2 n m1=1 Ω m1 σ (m1) z , ( 59
)
and the corresponding Suzuki-Trotter decomposition of the evolution operator U (∆t) for the (small) timestep of duration ∆t has the form

U (∆t) ≡ exp (-iH∆t) ≈ e -iH Z ∆t/2 e -iH Y ∆t/2 e -iH X ∆t e -iH Y ∆t/2 e -iH Z ∆t/2 , ( 60
)
ensuring the overall time discretization error of the order (∆t) 2 . The evolution operator over many time steps is a product of elementary operators U (∆t).

Each term in the sum representing the Hamiltonian H X (and, similarly, H Y and H Z ) commutes with all other terms, therefore

exp (-iH X ∆t) = n m1,m2=1 exp (-iJ m1,m2 ∆t σ (m1) x σ (m2) x ) = n m1,m2=1 cos (J m1,m2 ∆t) -iσ (m1) x σ (m2)
x sin (J m1,m2 ∆t) .

(61) Each term in this direct product acts on the wavefunction |ψ(t) in a straightforward manner: the entries of the array that represents the wavefunction turn into linear combinations of themselves. Similar direct-product representation holds for H Y and H Z as well, such that the action of the total evolution operator U (∆t) is easy to compute, without the need to calculate or store 2 N × 2 N matrices.

In order to represent the situation where the non-initialized part of the system is in the completely mixed state, but avoid using the density matrix explicitly (which would imply dealing with 2 n × 2 n matrix instead of the single array of the size 2 n ), we represent the completely mixed state as a wavefunction with random entries [1,2]. Specifically, we sampled the real and the imaginary parts of each entry of the corresponding wavefunction independently from Gaussian distribution with zero mean and unit variance, and then normalized the resulting wavefunction to one. In this way, for instance, the situation where the first and the second qubit are both initialized in the state |0 , while the rest of the system is in completely mixed state, i.e. when the system's density matrix is

ρ = |0 0| ⊗ |0 0| ⊗ I n-2 2 n-2 , ( 62
)
where I n-2 is an identity matrix of the size 2 n-2 × 2 n-2 , is represented using the total wavefunction in the form

|ψ = |0 ⊗ |0 ⊗ |ψ (r) n-2 , ( 63
)
where the random state |ψ (r)

n-2 of the remaining n -2 qubits is generated as described above. Such an approximation provides high accuracy, of the order of exp (-n/2), due to the measure concentration phenomenon [3].

Further improvement in accuracy was achieved by averaging the values of the relevant observables over M = 189 independent realizations of the random wavefunction (as well as other random quantities, see below), which reduced the error by an additional factor of the order ∼ 1/ √ M ≈ 0.07. The accuracy was also independently controlled by estimating the variance in the calculated values of the observables, and ensuring that this variance remains much smaller than the statistical error caused by the shot noise produced by sampling the relevant observables for each qubit.

B. Non-unitary evolution

The first non-unitary component of the system's evolution, dephasing of the m 1 -th qubit on the timescale T * 2,m1 , caused by its random static frequency shift, is modeled by directly reproducing the underlying physical picture. Namely, we assumed that the actual frequency Ω m1 of the m 1 -th qubit, see Eq. 57, is a sum of two contributions: the nominal value ωm1 , and a random shift β m1 that remains constant during the system's evolution. The values of β m1 were independently sampled from Gaussian distributions with zero mean and variance b 2 m1 , which can be different for different qubits. The parameter b m1 determines the dephasing time T * 2,m1 of the m 1 -th qubit: if this qubit were uncoupled from the rest of the system, then, after averaging over β m1 , its transverse (x-and y-) components would undergo Gaussian decay with time dependence exp (-b 2 m1 t 2 /2), i.e. b m1 = √ 2/T * 2,m1 . As mentioned above in subsection 2 A, the evolution of the system was repeated M = 189 times; each time we used different realizations of the set of the random frequency shifts β m1 (as well as other random quantities, such as e.g. different realizations of the random wavefunction, also see below). Within this approach, for each particular realization of the parameters β j , the evolution of the system is unitary, and can be simulated using the system's wavefunction as described in subsection 2 A above, while the non-unitary decay occurs due to averaging of the relevant observables over different realizations of the random parameters β m1 (along with other random quantities).

Note that the dephasing caused by averaging over the static random frequency shifts, with its characteristic Gaussian-like decay of the density matrix elements, cannot be described via Lindblad operators. It is an example of non-Markovian evolution, in the sense that it cannot be described by a set of first-order differential equations (with respect to time t), which would include only current values of the (averaged) elements of the system's density matrix ρ(t); in other words, the future values of the (averaged) density matrix elements, at times t + s (s > 0), are not completely determined by their current values at the moment of time t. At the same time, the static noise processes β m1 (t) representing the random frequency shifts are, of course, Markovian random processes, sastisfying the Chapman-Kolmogorov equation.

The two other components of the non-unitary evolution, addressed below, are Markovian, and can be described using the Lindblad operators. However, in order to avoid dealing with the density matrix, these components were also modeled by employing the random processes and calculating the averages of the relevant observables.

The second non-unitary component of the evolution corresponds to the Markovian dephasing, and can be described via the set of Lindblad superoperators corresponding to the phase damping channel. For an isolated qubit, this would lead to exponential decay of the transverse components of the m 1 -th qubit, having the form exp (-t/T 2,j ). This kind of dephasing, being Markovian, can be described by a set of first-order differential equations, generalizing the well-known Bloch-Redfield equations [4], which include only the current values of the elements of the system's density matrix ρ(t), such that the future values of the density matrix elements, at times t + s (s > 0), are completely determined by their current values at the moment of time t.

This decay was modeled by taking the z-rotation of the m 1 -th qubit produced by Ω m1 (i.e., produced by the action of the operator exp (-iH Z ∆t) in Eq. 60), and adding to it another time-dependent rotation around the z-axis by the angle Λ m1 (t). For each time step of duration ∆t, the values Λ m1 (t) were sampled randomly, indepedently of each other and of their previous values, from Gaussian distribution with zero mean and variance 4g 2 m1 τ m1 ∆t. This choice for the quantity Λ(t) can be visualized as a rotation induced by a time-dependent frequency shift δω m1 (t), which is represented by an Ornstein-Uhlenbeck noise process with the correlation function δω m1 (t) δω m1 (t + s) = g 2 m1 exp (-|s|/τ m1 ), in the limit where the correlation time τ m1 is much smaller than ∆t, while the magnitude g m1 is large (formally, τ m1 → 0 and g m1 → ∞), but the combination g 2 m1 τ m1 = 1/T 2,m1 remains finite. For an isolated qubit, the average evolution under the influence of such Ornstein-Uhlenbeck noise δω m1 (t) is known [5] to produce exponential decay of the qubit's transverse (x-and y-) components with the decay time T 2 . Again, for each particular realization of the time-dependent random process Λ m1 (t), the evolution of the system is unitary, and can be simulated using the system's wavefunction as described in subsection 2 A (provided, of course, that ∆t T 2,m1 , to ensure accuracy of the Suzuki-Trotter decomposition), while the non-unitary decay occurs due to averaging over different realizations of the noise.

The third non-unitary component, describing exponential relaxation of the m 1 -th qubit towards the state |0 on a timescale T 1,m1 , was simulated in a similar manner, by representing the non-unitary evolution via averaging over many realizations of a random unitary evolution, employing the approach described in Ref. [6], with some modifications improving the accuracy. Namely, at each time step, we calculated the probability p m1 for the m 1 -th qubit to make a transition ("quantum jump") from the state |1 to the state |0 ; the corresponding value is p m1 = w m1,1 µ 2 m1 , where µ m1 = 1 -exp (-∆t/T 1,m1 ), and w m1,1 is the total probability of the system to be in the subspace corresponding to the j-th qubit in the state |1 . This transition was implemented with the probability p m1 at each time step: the part of the system's wavefunction corresponding to the m 1 -th qubit in the state |0 was replaced by its complement, i.e. by the part corresponding to the m 1 -th qubit in the state |1 , multiplied by the factor µ m1 , and the part of the wavefunction corresponding to the j-th qubit in the state |1 was set to zero. Alternatively, with the probability 1-p m1 at each time step, the part of the wavefunction corresponding to the m 1 -th qubit in the state |1 was multiplied by the factor exp [-(1/2) ∆t/T 1,m1 ], while its complement was left unchanged. These transformations were applied to the wavefunction in succession, for all qubits (for all m 1 = 1, . . . n), and the resulting modified wavefunction was normalized back to 1. Since all these transformations commute with the action of the operator exp (-iH Z ∆t) in Eq. 60, they were applied at the end of each unitary time-step evolution, after application of the operator U (∆t) given by Eq. 60, in parallel with the action of the operators exp (-iH Z ∆t) or exp (-iH Z ∆t/2).

Note that this implementation corresponds to the application to the wavefunction of the Krauss operators E 0 or E 1 (see Ref. 7), describing the amplitude damping quantum channel, with the corresponding probabilities, where E 1 corresponds to the event of the "quantum jump", and E 0 corresponds to the absence of it.

Supplementary Note 3: Numerical simulation of superconducting qubit platform

From the discussion in the main text we are simulating a 2D grid of qubits that interact only with the nearest neighbours. The coupling between two neighbouring qubits through a coupler can be described by a Hamiltonian. In our notations it can be rewritten as

H = k=i,j a (k) z H (k) z + a (i,j) xx H (i,j) xx + a (i,j) yy H (i,j) yy , ( 64
)
where we introduce the notations

a (i) z = 1 2 ω(i) , a (j) z = 1 2 ω(j) , a (i,j) xx = a (i,j) yy = 1 2
g i g j ∆ + g ij , ( 65
)
H (i) z = σ (i) z , H (j) z = σ (j) z , H (i,j) xx = σ (i) x σ (j)
x , H (i,j) yy = σ (i) y σ (j) y .

Thus, we define the 16 qubits 2D grid, where we generate ω 1 , ω 2 , ω c , g 1 , g 2 , g 12 from Gaussian distribution with mean and variance N (0, 1). The parameters a

(i) z , a (j) z , a (ij) xx , a (ij)
yy we estimate in our simulation and the rates of decay of T 1 , T 2 and T ⋆ 2 are given in the Table III. The observables and initial states, isolating the desired coefficients a 16. The rates of decay of T1, T2 and T ⋆ 2 , corresponding to i's qubit are given in the fifth and six's columns, respectively.

(i) z , a (j) z , a (i,j) xx (ij) a (ij) xx , a (ij) yy [kHz] (i) a (i) z [kHz] T1 [µs] T2 [µs] T ⋆ 2 [µs] 1, 2 1.
α i , αi ∈ {x, y, z}, i = 1, . . . ,
and a (i,j) yy , are given in Table IV. One can see, that we need three starting states, namely ρ 01 = ρ

(i) x ⊗ ρ (j) x ⊗ I 2n-2 2 2n-2 , ρ 02 = ρ (i) y ⊗ ρ (j) z ⊗ I 2n-2 2 2n-2 and ρ 03 = ρ (i) z ⊗ ρ (j)
z ⊗ I 2n-2 2 2n-2 to isolate all four unknown coefficients. We measure the expectation values of observables in different times. Next, the time traces of these expectation values are fitted, using the polynomial interpolation method, and the derivatives estimation is preceded. Finally, using (39), ( 44), ( 28) and ( 28), the estimates of the coefficients a ,j) yy for the pair of (i, j) are obtained. We repeat this process for all pairs of interacting qubits to obtain all coefficients of the Hamiltonian of the 2D grid.

(i) z , a (j) z , a (i,j) xx , a (i
In the presence of the noise, the observables and initial states required to isolate the Lindbladian coefficients are given in Table V. One can see, that we need three extra starting states in the presence of the Lindbladian noise, namely ρ 04 = ρ

(i) z ⊗ ρ (j) y ⊗ I 2n-2 2 2n-2 , ρ 05 = ρ (i) y ⊗ ρ (j) y ⊗ I 2n-2 2 2n-2 and ρ 06 = ρ (i) z ⊗ ρ (j) x ⊗ I 2n-2 2 2n-2 to find L (i) µν , µ, ν ∈ {x, y, z}. a (i) α i {O, ρ (i,j) τ i ,τ j } Equation a (i) z {σ (i) x ⊗ σ (j) z , ρ (i,j) y,z }; {σ (i)
x , ρ

(i,j) y,z } (39) a (j) z {σ (i) x ⊗ σ (j) y , ρ (i,j) x,x }; {σ (i) y , ρ (i,j) x,x } (44) a (ij) xx {σ (i) x ⊗ σ (j) y , ρ (i,j) y,z }; {σ (i) y ⊗ σ (j) y , ρ (i,j) x,x } (28) a (ij) yy {σ (i) x ⊗ σ (j) y , ρ (i,j) z,z }; {σ (i) x ⊗ σ (j) y , ρ (i,j) y,z }; {σ (i) y ⊗ σ (j) y , ρ (i,j) x,x } (28)

Supplementary Table IV. In this table the minimal selection of the pairs {O, ρ

(i,j) τ i ,τ j } is presented for our specific example. The first column represents the type of the estimated Hamiltonian (64) parameters a (i)

α i , a (ij)
α i ,α j , αi, αj ∈ {x, y, z}. In the third column the number of equation for every parameter is provided, depending from the pairs of the observable O and the initial state ρ (i,j)

τ i ,τ j = ρ (i) τ i ⊗ ρ (j)
τ j ⊗ ρN-2, τi, τj = {x, y, z}, given in the second column. To estimate all four parameters we need only three initial states: ρ (i,j)

x,x , ρ (i,j) z,z and ρ

(i,j) y,z . L (i) µν {O, ρ (i,j) τ i ,τ j } Equation L (i) xx, L (i) yy , L (i) zz {σ (i)
x , ρ (i,j)

x,x }; {σ

(i) y , ρ (i,j) y,z }; {σ (i) z , ρ (i,j) z,z } (55) L (j) xx , L (j) yy , L (j) zz {σ (i) x , ρ (i,j) z,x }; {σ (i) y , ρ (i,j) z,y }; {σ (i) z , ρ (i,j) y,z } (56) L (i) xy {σ (i) x ⊗ σ (j) z , ρ (i,j) y,z } (49) L (i) xz {σ (i) x ⊗ σ (j) η , ρ (i,j) z,z } (50) L (i) yz {σ (i) y ⊗ σ (j) η , ρ (i,j) z,z } (51) L (j) xy {σ (i) z ⊗ σ (j) x , ρ (i,j) z,y } (52) L (j) xz {σ (i) y ⊗ σ (j) x , ρ (i,j) y,z } (54) L (j) yz {σ (i) y ⊗ σ (j) z , ρ (i,j) y,y } (53)
Supplementary Table V. In this table the minimal selection of the pairs {O, ρ (i,j) τ i ,τ j } is presented for our specific example. The first column represents the type of the estimated Lindbladian parameters L (i) µν , µ, ν ∈ {x, y, z}. In the third column the number of equation for every parameter is provided, depending from the pairs of the observable O and the initial state ρ (i,j)

τ i ,τ j = ρ (i) τ i ⊗ ρ (j)
τ j ⊗ ρN-2, τi, τj = {x, y, z}, given in the second column. To estimate all parameters we need only three extra states to the ones given by the previous table, namely: ρ

(i,j) z,x , ρ (i,j)
z,y and ρ (i,j) y,y .

Supplementary Note 4: Approximating local time evolutions by polynomials

One of the main points behind our method is the fact that the time evolution of local observables at constant times is well-approximated by polynomials. The purpose of this section is to make this assertion precise.

Before we do that, let us set some notation. Given a system of n qubits on a D dimensional lattice Γ, we let L Λ : M 2 n → M 2 n be a Lindbladian which models the time evolution of the system in the Heisenberg picture. Note that in the supplementary information we consider a slightly more general class of evolutions than in the main text. There, we restricted to evolutions whose Hamiltonians were short range with two-body interactions and the noise acted on at most one qubit at a time. Here, in contrast, we will also consider k-local evolutions with long range.

We will assume that this Lindbladian can be written as:

L Λ = A⊂Λ L A , ( 66
)
where L A is a Lindbladian only acting on the qubits in A. Given some graph G = (V, E) on n vertices, we will say that L is k-local if L A = 0 only if A is a subset of vertices of G containing at most k vertices. Furthermore, we will say that L Λ is locally bounded if there is a constant g > 0 such that for all B ⊂ Λ we have that:

A⊂Λ:A∩B̸ =∅ L A ≤ g|B|. ( 67
)
This condition is satisfied if e.g. L is a local Lindbladian on a D-dimensional lattice. In that case, we have g = O(D). However, this condition is also fulfilled for generators with algebraically decaying tails that decay at least like x - (D+δ) for some δ > 0. Moreover, for ease of notation we will let for a region B ⊂ Λ

L B = A⊂B L A ( 68
)
be the generator restricted or truncated to a subregion B. Furthermore, given the lattice Λ, some region X ⊂ V and r > 0, we will denote by Λ r (X) the set of vertices that are a distance at most r from X:

Λ r (X) = {v ∈ V : ∃x ∈ X s.t. dist(x, v) ≤ r}. ( 69
)
Here the distance is according to the lattice. We will also require some norms for superoperators. Given a superoperator Φ :

M 2 n → M 2 n we define for p, q ≥ 1 Φ p→q = sup X∈M 2 n Φ(X) q X p , ( 70
)
where • p corresponds to the Schatten p-norm. Also note that p = ∞ corresponds to the operator norm, for which we will often drop the +∞ and only write as • . We will also consider completely bounded versions of these norms, which are given by:

Φ p→q,cb = sup d≥1 sup X∈M 2 n ⊗M d Φ ⊗ id d (X) q X p . ( 71
)
It is then simple to see that derivatives of locally bounded evolutions can only increases with the size of the region they are defined on: Then for all t ≥ 0:

f (k) B (t) ≤ (g|B|) k . ( 72
)
In particular, for any 0 < t < t max we have that:

f B (t) - d k=0 f (k) B (0) k! t k ≤ (t max g|B|) d+1 (d + 1)! . ( 73
)
Proof. The proof is elementary. Note that:

f (k) B (t) = tr e L B ((tL B ) k (O))ρ . ( 74
)
Now, by Hölder's inequality we have that:

f (k) B (t) ≤ e tL B ((tL B ) k (O)) ∞ ρ 1 (1) ≤ t k e tL B ∞→∞ L B k ∞→∞ (2) ≤ g k |B| k ,
where in (1) we used the submultiplicativity of the operator norm, i.e. Φ 1 Φ 2 ∞→∞ ≤ Φ 1 ∞→∞ Φ 2 ∞→∞ for all linear maps Φ 1 , Φ 2 . In (2) we used the fact that for any quantum channel e tL ∞→∞ = 1 and the fact that the Lindbladian is locally bounded with constant g. The estimate in Eq. ( 73) then immediately follows from Taylor's remainder theorem.

We then immediately have: Corollary 4.1. In the same setting as Lemma 4.1 it holds that for any given ϵ > 0 and t max > 0 s.t.

log log(ϵ -1 )2et max g|B| ≥ 1 (75)
there is a polynomial p of degree

d = 2et max g|B| log(ϵ -1 ) -1 (76)
such that for all 0 ≤ t ≤ t max we have that

|f B (t) -p(t)| ≤ ϵ. ( 77
)
Proof. It follows from Sitrling's approximation that the error in Eq. ( 73) is bounded by

f (t) - d k=0 f (k) (0) k! t k ≤ 1 √ 2πd et max g|B| d + 1 d+1 . ( 78
)
It is then easy to see that picking d = 2et max g|B| log(ϵ -1 ) -1 is sufficient to ensure that the error in (73) is at most ϵ. Indeed, plugging in the value of d into Eq. (78), we get:

f (t) - d k=0 f (k) (0) k! t k ≤ 1 √ 2πd 1 log(ϵ -1 ) 2etmaxg|B| log(ϵ -1 ) = (79) 1 √ 2πd exp -log(log(ϵ -1 )) log(ϵ -1 )2et max g|B| = 1 √ 2πd ϵ log log(ϵ -1 )2etmaxg|B| ≤ ϵ. ( 80
)
Thus, the truncated Taylor expansion yields the desired polynomial.

We conclude from Lemma 4.1 and Cor. 4.1 that local time evolutions are well-approximated by polynomials whose degree grows like the size of the region times the maximal time of evolution.

Also note that the estimate in Eq. ( 79) is quite loose and shows that for d as in Eq. ( 76) the error decays like a polynomial of high-degree in ϵ. But that rough approximation will be sufficient for our purposes.

Cor. 4.1 is an important step to prove our Thm. 4.1, but still does not correspond to the exact statement we wish to prove. This is because Cor. 4.1 is a statement about the local evolution, whereas Thm. 4.1 is a statement about the global evolution being well-approximated by a polynomial of small degree. The strategy to go from the local to the global evolution, is to show that for the (local) observables required for our protocl, the local evolution approximates the global one well.

Our main tool to show this approximatability of expectation values are Lieb-Robinson bounds [8][9][10][11], which exactly give conditions under which the local time evolution and the global one are close for small enough times and local observables. In order to provide a self-contained presentation, we include a brief introduction to Lieb-Robinson bounds in Sec. 7 of this appendix.

In fact, there are various ways of quantifying that this idea of local approximability and, thus, LR-bounds come in various forms. The version that we are going to work with here and which is discussed in detail in Sec. 7, considers an observable O Y initially supported on in a region Y . For a region B ⊃ Y we set R = dist(Λ\{B}, Y )/k, where k is the locality of the generator, then it is shown in Lemma 7.1 that indeed

(e tLΛ -e tL Λr (Y ) )(O Y ) ∞ ≤ (e vt -1)h(R), (81)
where h : R + → R + is a monotonically increasing function such that lim R→+∞ h(R) = 0 and v is some constant that depends on the generator, usually called the LR-velocity. On the other hand, the decay of the function h typically depends on how fast the interactions in the system decays spatially (i.e. if it is strictly local, exponentially decaying in the distance or even algebraically decaying) and the geometry of the underlying lattice. However, the important point for our purposes is that it does not depend on the system's size. For the specific case of short-range Hamiltonians discussed in the main text, we have that h(R) = e -µR for some constant µ > 0. For algebraically decaying evolutions we usually have h(x) = x k for some k ∈ R + . We refer again to Sec. 7 for a discussion of various LR-bounds available in the literature. But from Eq. ( 81) we immediately conclude that the values of the expectation values of global and local evolutions are well-approximated by each other. More precisely: Proposition 4.1. Let O Y be an observable supported on some region Y , ϵ > 0 and t max be given. Assume Eq. (81) holds for the time evolution L Λ and a function h. Let l > 0 be given by

l = h -1 ϵ e vt max -1 . ( 82
)
Then we have for Λ l (Y ) and all 0 < t < t max and any initial state ρ that

| tr e tLΛ -e tL Λ l (Y ) (O Y )ρ | ≤ ϵ. ( 83
)
Proof. The claim follows directly from Eq. (81) or Lemma 7.1 and a Hölder inequality. Indeed, for the value of l in Eq. ( 83), we obtain from Eq. (81) after some simplification that

(e tLΛ -e tL Λr (Y ) )(O Y ) ∞ ≤ ϵ. ( 84
)
Note that in the case of short-range systems we have that h -1 (x) = µ -1 log(x). From now on we will suppress the terms of order log(log(ϵ -1 )) or higher from the equations and denote bounds where we do this with Õ. Thus, combining 4.1 with Prop. 4.1 we conclude that: Theorem 4.1. Let L Λ be a locally bounded Lindbladian on a D-dimensional regular lattice with constant g. Moreover, let t max , ϵ > 0 be given and O Y is an observable such that O Y ≤ 1 and O Y is supported on a constant number of qubits. Assume that L Λ satisfies Eq. (81). Then there is a polynomial p of degree

d = Õ h -1 ϵ e vt max -1 D t max log(ϵ -1 ) (85)
such that for all 0 ≤ t ≤ t max :

tr e tLΛ (O Y )ρ -p(t) ≤ ϵ, ( 86
)
and

p ′ (0) = tr [L Λ (O Y )ρ] . ( 87
)
Proof. It follows from Prop. 4.1 that a region of radius

l = Õ h -1 ϵ 2(e t max -1) (88)
is enough to approximate the time evolution of e tLΛ (O Y ) up to ϵ/2. If the original region Y has a constant number of qubits, then for a D-dimensional lattice we have |Λ l (Y )| = O(l D ). It then follows from Cor. 4.1 that for regions of this size, it is sufficient to pick a degree that is

O h -1 ϵ 2(e t max -1) D t max log(ϵ -1 )) (89)
to approximate the expectation value of the local evolution up to an error ϵ/2. This concludes the proof by a triangle inequality. Eq. ( 87) is clear from properties of the truncated Taylor series.

Thus, we see that as long as the time evolution of the system satisfies a Lieb-Robinson bound, we can approximate the expectation value of a local observable as a function of time by a polynomial whose degree is dictated by how fast the function h decays and the maximal time of the evolution. In particular, for ϵ -1 = O(1) and time t max = O(1), we conclude that the degree of the polynomial is independent of the system's size.

In Thm. 4.1 we established that we can approximate the function f : t → tr e tLΛ (O Y )ρ well by a polynomial for constant times. However, to estimate the parameters of the Hamiltonian we are ultimately interested in the derivative of f at time 0. We will later show in Sec. 5 that for the special case of polynomials of bounded degree, a good recovery of the polynomial also implies a good recovery of the derivative. In particular, as for the polynomial in Lemma 4.1 we have that p ′ (0) = f ′ (0), it is sufficient to argue that any two polynomials that approximates the curve f (t) up to sufficiently large precision in a sufficiently large number of points must have close derivatives at 0 as well. This will be the subject of the next section and proved in Prop. 5.1.

Supplementary Note 5: (Robust) polynomial interpolation and derivative estimation

In this section of the appendix, we are going to review a result in the literature [12] that shows how to perform polynomial interpolation in a robust way even in the presence of outliers. Furthermore, we will show that good polynomial interpolation also implies a good approximation of derivatives of the polynomial, which is our end goal. We will use and review the results and algorithms of [12] for the robust polynomial interpolation and resort to Markov brothers' inequality [13] for estimating the error on the derivatives.

Let us start by briefly recalling the technical problems we wish to overcome. We assume we are able to approximate the expectation value of f (t) = tr ρe tL (O) for some suitably-picked initial state ρ and time-evolved observable O. As argued in Sec. 4, the function f is well-approximated by a low degree polynomial p whenever the time evolution is generated by a local Lindbladian. Moreover, as shown in Sec. 1, by suitably choosing the observable and initial state, we can easily read off the value of the coupling of the Hamiltonian from the value of f ′ (0). Thus, our goal is to find a polynomial p that approximates f from values of f (t i ) for some t i ∈ [a, b] for then use p to infer f ′ (0).

It is well-known that if p is a polynomial of degree d, then it is uniquely determined by its values at d + 1 points. Thus, one could naively expect that having access to f (t 1 ), . . . , f (t d) points for d ∼ d times is sufficient to reconstruct f . However, the present situation exhibits three challenges that need to be overcome to ensure that we can reliably apply polynomial interpolation methods and recover p from points f (t i ):

1. we can only estimate f (t), and not p(t). And the value of f only approximates that of p up to some error ϵ a , as discussed in Thm. 4.1.

2. we do not have access to the value of f (t) directly, but can only approximate it to a precision ϵ s by sampling from the output of the device at time t O(ϵ -2 s ) times. 3. we are interested in the value of p ′ (0) and not in the polynomial p itself. Thus, we need to ensure a small error in estimating the derivative.

To deal with the first two problems the polynomial interpolation technique we use has to be robust to the noise stemming from both the approximation error from the polynomial approximation and the statistical noise. To deal with the third issue, we will show that we need to pick the final and initial time of the interpolation in a judicious manner.

To obtain some intuition about how to pick the times, let us consider the case of estimating the derivative of a quadratic polynomial at 0. I.e., if we have two linear functions p, p that are ϵ close in some interval [a, b], how well can we infer the derivative of p at 0 from that of p? First, note that if the interval [a, b] is very small, then two functions can differ by ϵ and their derivatives can still differ by ϵ/(b -a) even for linear functions. This indicates we probably do not want to pick b -a too small. However, as we know that increasing b also implies that, in our setting, we need to increase the degree of the polynomial for the fitting, which in turns increases the number of points we need to estimate, this hints at the fact that picking b -a of constant order will be optimal.

On the other hand, it is also clear that the closer a is to 0, the more information about the value of p ′ (0) we can infer from the interpolation. Thus, this discussion suggests that picking a as close to 0 as possible and b of constant order should give the best results. We will prove this intuition later in this section, but first will discuss robust interpolation.

Directly interpolating through the noisy data can be an unstable procedure if we do not pick the interpolating points wisely and perform a suitable regression. Recent results have shown how to perform polynomial interpolation in an essentially optimal fashion in a robust way even with a fraction of the points being outliers [12]. Let us now review the results of [12], but note that we expect that less sophisticated interpolation techniques, i.e. based on least squares regressions, should have a similar performance in practice. The reason we use the interpolation described in [12] is that it comes with very strong theoretical guarantees under minal assumptions.

We will now assume we wish to estimate a polynomial p : [-1, 1] → R of degree d, as this corresponds to the setting of [12]. Note that for Hamiltonian learning, we will be interested in the case where the domain is of the form [a, b] for a, b ≥ 0. However, we can simply shift and rescale the domain to [-1, 1]. When we summarize our results later, we will dicuss the effect of this rescaling explicitly. We will assume we are given access to m random samples (x i , y i ) of points such that a fraction of at least α > 1/2 of them satisfies for some σ > 0 that:

p(x i ) = y i + w i , |w i | ≤ σ. ( 90
)
There are results available for various different ways of sampling the points x i . However, the best available sample complexity is given by sampling from the Chebyshev measure, which has density

1 π √ 1 -x 2 (91)
on the interval [-1, 1]. We then have:

Theorem 5.1 (Robust polynomial interpolation). Let p : [-1, 1] → R be a polynomial of degree d and assume we are given m samples (x i , y i ) such that a fraction α > 1/2 of them satisfies Eq. (90) for some σ > 0. Moreover, suppose that the x i were sampled independently and at random from the Chebyshev measure. Then for any δ > 0

m = O d log d δ ( 92
)
samples suffice to with probability of success at least 1 -δ recover a polynomial p that satisfies:

max x∈[-1,1] |p(x) -p(x)| ≤ 3σ. ( 93
)
Moreover, p can be computed in time polynomial in the number of samples m.

Proof. We refer to [12,Corollary 1.5] for a proof and note that we obtain the statement by setting the parameter ϵ < 1/4 in their statement.

We note that the same result holds for random points picked from the uniform measure with m = O(d 2 ).

The result above solves our problem of robust polynomial interpolation outlined in points 1 and 2. It shows that it if we can ensure that we can approximate sufficiently many points of the polynomial up to some σ, then we also recover the whole polynomial up to some error proportional to σ. Moreover, the number of samples required only has a logarithmic overhead in d when compared with the case where we know the points exactly. As we will see later, for our puproses it will be important to choose d to be small. Thus, in a nutshell, we see that Thm. 5.1 ensures that we can reliably and robustly perform polynomial interpolation by only a small overhead when compared to when we know the points exactly.

We will later describe in more detail the algorithm given in [12] whose output satisfies the promises of Thm. 5.1. However, before that we will show how the condition in Eq. ( 93) ensures that we can also recover the derivative of the polynomial as long as the degree d is small.

To do that, we will resort to Markov brothers' inequality, which we restate now for completeness.

Lemma 5.1 (Markov brothers' inequality). For d, k ∈ N define the constant C M (d, k) to be given by

C M (d, k) = d 2 d 2 -1 2 d 2 -2 2 • • • d 2 -(k -1) 2 1 • 3 • 5 • • • (2k -1) . ( 94
)
Then for any polynomial p of degree d we have that:

max x∈[-1,1] p (k) (x) ≤ C M (d, k) max x∈[-1,1] |p(x)| . ( 95
)
Proof. We refer to [14, Theorem 1.2] for a proof and discussion of this result.

Note that the value of CM (d, k) increases exponentially with d for k constant. We remark that having further promises on the structure of the polynomial, such as the location of its zeros, can greatly improve this estimate. We once again refer to [14,Chapter 1] for a discussion on this. It would be interesting to see if recent results on the analyticity of the partition function [15] could be used in our context to also improve the error estimate. However, this general bound will suffice for our purposes.

Let us now discuss the rescaling. It is easy to see that for polynomials defined on some interval [a, b], the polynomial p(x) = p( b-a 2 x+ a+b 2 ) is defined on [-1, 1] and we can use this simple transformation to obtain a variation of Eq. ( 95) for polynomials defined on general intervals. Indeed, applying Eq. ( 95) to p, it follows from a straightforward application of the chain rule that

max x∈[a,b] p (k) (x) ≤ 2 (b -a) k C M (d, k) max x∈[a,b] |p(x)| . ( 96
)
From this we conclude that: and

|p ′ (0)| ≤ ϵ d k=1 2 b -a k a k-1 C M (d, k) (k -1)! . ( 99
)
Proof. It follows from Eq. ( 97) and Markov brothers' inequality that

p (k) (a) ≤ 2 b -a k C M (d, k)ϵ. ( 100
)
By a Taylor expansion we know that for x ∈ [0, a]:

p(x) = d k=0 p (k) (a) (x -a) k k! . ( 101
)
The claim in Eq. (98) then follows by combining this expansion with Eq. (100) and a triangle inequality. Similarly we have

p ′ (0) = d k=1 p (k) (a)(-1) k a k-1 (k -1)! , ( 102
)
for which a similar argument yields Eq. (99).

The proposition above essentially allows us to control to what extent the derivative of a polynomial at 0 can deviate from 0 given that the polynomial is small on another interval [a, b]. We can then apply it to the polynomial p -p, as in Eq. ( 93) to control the error we make by estimating the derivative at 0 by evaluating p(0).

By combining the arguments above we conclude that:

Proposition 5.1 (Precision and number of samples for robust interpolation). Let p be a polynomial of degree d. For some 0 < a < b define

E(a, b, d) = d k=1 2 b -a k a k-1 C M (d, k) (k -1)! . ( 103
)
Then for δ > 0 sampling

m = O(d log(dδ -1 )) (104) i.i.d. points (x i , y i ) from the Chebyshev measure on [a.b] satisfying p(x i ) = y i + w i , |w i | ≤ σ
for at least a fraction α > 1 2 of the points is sufficient to obtain a polynomial p satisfying

p ′ (0) -(p) ′ (0) ≤ 3σE(a, b, d). ( 105
)
with probability of success at least 1 -δ.

Proof. From Thm. 5.1 we know that this number of samples suffices to obtain the polynomial p satisfying Proof. It is easy to see that we have:

C M (d, k) ≤ d 2k k!! , ( 109
)
where

k!! = 1 × 3 × 5 × • • • × (2k -1)
is the double factorial. Thus, we see from this and Eq. (99) that by our choice

a = 1 d 2 , b = 2 + a (110)
we have from Eq. ( 105) that the estimated polynomial p(0) satisfies

|p ′ (0) -p ′ (0)| ≤ 3σ d k=1 2 (a -b) k a -1 (ad 2 ) k k!!(k -1)! ≤ 3σ d k=1 a -1 1 k!!(k -1)! , ( 111
)
where we used the fact that a ≤ d -2 . Thus, as

d k=1 1 k!!(k -1)! ≤ e, ( 112
)
we conclude that with this choice of parameters we have

|p ′ (0) -p ′ (0)| ≤ 3eσa -1 , ( 113
)
which gives the claim by our choice of σ.

We will discuss in Sec. 6 how to specialize the discussion and results above to the scenario of Hamiltonian learning.

Supplementary Note 6: Choice of parameters and performance guarantee of the protocol

Let us now combine the results from Sections 4 and 5 to see how to pick the various parameters of the algorithm to ensure a good recovery of the couplings of the Hamiltonian.

More precisely, given a coupling parameter a α of L Λ we will be interested in estimating the sample complexity of obtaining an estimate âα satisfying

|a α -âα | ≤ ϵ (114)
with high probability for some given error ϵ. As extensively discussed by now, we can easily reduce estimating the couplings to estimating derivatives of time evolutions of local observables.

As expected, we will see that the main parameters we need to control are the maximal observation time t max and the initial time of measurement t 0 . This is showcased in the following Theorem: Theorem 6.1 (Choice of final and initial time). Let L Λ be a locally bounded Lindbladian on a D-dimensional regular lattice with growth constant g = O (1). Let O be an observable of constant support and ρ and arbitrary quantum state. Let ϵ > 0 be given. Assume that L Λ satisfies Eq. (81) for some function h. Then picking t 0 as

t 0 = O   h -1 ϵ 2(e 2.5v -1) D log(ϵ -1 ) -2   ( 115
)
and t max = 2 + t 0 and measuring the expectation value f (t) = tr e tLΛ (O)ρ for

m = Õ h -1 ϵ 2(e 2.5v -1) D log(ϵ -1 )) ( 116
)
random times t i ∈ [t 0 , t max ] up to precision O(ϵ) is sufficient to obtain an estimate ( f ) ′ (0) satisfying ( f ) ′ (0) -tr [L Λ (O)ρ] ≤ 3eϵt -1 0 = O   ϵ h -1 ϵ 2(e 2.5v -1) D log(ϵ -1 ) 2   . ( 117
)
In particular, this estimate can be obtained from

Õ(ϵ -2 log(δ -1 )) (118)
samples from the time evolved state ρ with probability of success at least 1 -δ.

Proof. First, note that by Lemma. 4.1, if we pick t max as described above, then a polynomial p of degree

d = O h -1 ϵ 2(e 2.5v -1) D log(ϵ -1 )) (119)
is sufficient to approximate the expectation value in the interval [0, t max ] up to an error ϵ/2, assuming that ϵ is small enough to ensure that t max ≤ 2.5. We will estimate the value of the polynomial at each point up to an error σ > 0, which is to be determined later. Thus, by inserting the bound on the degree d in Eq. ( 107), we need to estimate each value of the polynomial up to a precision σ = O(ϵ) to obtain an overall error of

O   ϵ h -1 ϵ 2(e 2.5v -1) D log(ϵ -1 ) 2   (120)
As we imposed that the precision with which the polynomial approximates the expectation values is ϵ/2, we can estimate the value of the polynomial for a given time up to an error O(ϵ) from O(ϵ -2 ) samples.

As we have to sample

O(d log(d)) = Õ h -1 ϵ 2(e 2.5v -1) D log(ϵ -1 ) (121)
points to perform the stable interpolation, we obtain the advertised sample complexity. VI. Scaling of different resources required to obtain a recovery up to additive error ϵ of a parameters of the evolution. We have only included the leading order term and α denotes the decay of the potential in space, whereas D the dimension of the lattice.

Algebraically (α ≥ 5D -1) ϵ -2 α-3D α-5D ϵ -D α-5D ϵ 2D α-5D

Supplementary Table

For the case of strictly local or exponentially decaying interactions we have that

h -1 ϵ 2(e 2.5v -1)
= poly log(ϵ -1 ).

In that case the sample complexity is of order Õ(ϵ -2 ). Thus, in this case we see that the inverse initial time t -1 0 and the number of points we need to sample from is polylogarithmic in ϵ. Furthermore, the sample complexity to obtain an error ϵ is also Õ(ϵ -2 ) up to polylogarithmic corrections.

For the sake of completeness, let us now discuss the conditions under which our protocol works beyond the setting of exponentially decaying or short-range interactions. From Eq. (117) the condition for our procedure to work becomes transparent: we need that

h -1 ϵ 2(e 2.5v -1) D log(ϵ -1 ) -2 = o(ϵ -1 ). (123)
Indeed, in this case we have the property that it is possible to suitably re-scale the error ϵ to ensure that the total precision is at some desired precision ε. For instance, let us assume that

h -1 ϵ 2(e 2.5v -1) = O(ϵ -r ), (124)
for some r > 0. As we discuss later, this is typically the case for algebraically decaying interactions. For such a LR-bound, we see that the resulting error in Eq. ( 117) is

O(ϵ 1-2Dr log(ϵ -1 ) 2 ). (125)
Ignoring the log(ϵ -1 ) term, we see that by picking ϵ = ε 1 1-2Dr we can ensure an error of order ε for the estimate. Thus, the growth of h -1 has to be at most r ≤ 1 2D and the sample complexity would also grow like Õ(ϵ -2-2Dr-r ), as we would need to sample Õ(d) = Õ(ϵ -r ) points up to precision ε 1 1-2Dr . Thus, we see that our protocol has a sample complexity that is independent of the system size to estimate one parameter and the expected ϵ -2 scaling for short range evolutions evolutions, up to log factors. For algebraically decaying interactions, however, the sample complexity has a worse scaling that depends on the exact decay of the potential, but still independent of system size.

We summarize the sample complexities, smallest initial time t 0 and number of different times steps we need for various different potentials in Table VI.

A. Algorithm for robust polynomial interpolation

Now that we have established that the results of [12] indeed allow us to estimate the derivative at 0, let us now describe their polynomial interpolation algorithm in more detail for completeness. The algorithm consists of two parts, one ℓ 1 regression and an iteration of ℓ ∞ regressions. Following [12], we will only consider the case in which we interpolate over [-1, 1]. But it is straightforward to also interpolate over other intervals by a suitable affine transformation of the domain, as discussed before.

a. ℓ 1 regression: before we define the ℓ 1 regression, we need to define the Chebyshev partitions: We also define P d to be the space of polynomials of degree at most d.

With these definitions at hand, we define the ℓ 1 regression solution as follows: Definition 6.2. Given a set of n points (x i , y i ) and m ∈ N, we define the result of the degree d ℓ 1 regression with m Chebyshev partitions p to be the polynomial

argmin p∈P d n i=1 |I j | mean xi∈Ij |y i -p(x i )| , ( 126
)
where P d is the set of polynomials of degree d.

Note that the optimization problem above is a linear program and, thus, can be solved efficiently. Solving the ℓ 1 regression problem with n = O(d log(d)) samples from the Chebyshev measure is guaranteed to give us a good solution on average. More precisely, as shown in [12, Lemma 1.2], the solution is guaranteed to satisfy

p -p ℓ1 = O(σ), ( 127
)
where as usual σ is the error in each estimate y i and

p -p ℓ1 = 1 -1 |p(x) -p(x)| dx. ( 128
)
However, the results of the previous sections required us to obtain a good solution in the • ℓ∞ distance, and in general

p -p ℓ1 = O(d 2 p -p ∞ ). (129)
Although, as commented in the last section, we are interested in the regime of polynomial of relatively small degree, by adding a ℓ ∞ regression iteration on top of the ℓ 1 regression, it is possible to get rid of this d 2 prefactor. b. ℓ ∞ regression: besides getting rid of the unwanted d 2 factor on the promise for the error of the ℓ 1 regression, adding a ℓ ∞ regression step also has the favourable feature of making the whole procedure more robust to outliers in the data. (130)

We define the result of the degree d ℓ ∞ regression with m Chebyshev partitions p to be the polynomial p ∈ P d

argmin p∈P d max j∈[m] | p ( x j ) -y j | . ( 131
)
Note that the problem in Eq. ( 131) also corresponds to a linear program and, thus, can be solved efficiently. The output of the ℓ ∞ regression algorithm is guaranteed to satisfy c. Iterating the ℓ ∞ step: the last step to obtain the desired robust polynomial interpolation is to iteratively apply the ℓ ∞ iteration step to the residual. More precisely, we first perform the ℓ 1 -regression on our data, obtaining a polynomial p0 . We can then define the new data points

p -p ∞ ≤ 2.5σ + 1 2 p ∞ ( 132
(x i , ỹ0 i = y i -p0 (x i )) (133)
and run the ℓ ∞ interpolation on this residual error, obtaining a polynomial p1 . From Eq. (129) and our promise on the output of the ℓ ∞ interpolation, we know that the result of the interpolation will satisfy

p -p1 ∞ ≤ 2.5σ + 1 2 O(d 2 σ). ( 134
)
But then we can iterate this procedure by just running the ℓ ∞ regression on

(x i , ỹ1 i = y i -p1 (x i )). (135)
Each time we run the interpolation on the residual, we exponentially reduce the error. By repeating the procedure O(log 2 (d)) times, we then arrive at a polynomial satisfying the promises of Thm. 5.1. Note, however, that the procedure used in Sec. III of the main text to demonstrate the viability of our method differs slightly from the ones discussed here. The main difference is that we used equally spaced time steps that were not random. However, in spite of this difference, we still obtained high quality solutions.

B. Certifying the degree of the polynomial in the fitting

Our approach requires us to fit the expectation values to a low-degree polynomial. However, it is natural to ask how to pick the degree of the polynomial we fit the data to and to certify that the degree was chosen correctly. A heuristic approach to this problem worked remarkably well in our numerics: we increased the degree d until the maximal difference between the interpolating polynomial we obtained and the data we collected would fall below O(ϵ), the target accuracy.

In Sec. 7 A we give a version of our polynomial approximation results with explicit constants. The exact required degree is given in Eq. ( 148). However, our numerics suggests that a polynomial of degree orders of magnitude smaller than given by the bound suffices to obtain a good approximation and offers greater numerical stability. We leave to future work to obtain more refined estimated on the degree required. in terms of the stabilized 1-to-1-norm T 1→1,cb = sup n T ⊗ id n 1→1 . We can then characterize L as finite range if µ(r) = 0 for r ≥ R > 0, exponentially decaying if µ(r) ≤ e -µr and algebraically decaying if µ(r) ≤ (1 + r) -α for α > 0 and state the following Lieb-Robinson-bound for Lindbladians: Theorem 7.1 (dissipative LR-bound [19]). Let L be a Lindbladian of the form (138), O Y an observable supported on Y ⊂ Λ and K

X : M 2 n → M 2 n with K X (id X ) = 0. Then K X (O Y (t)) ≤ K X ∞→∞,cb O Y min(|X|, |Y |)h(dist(X, Y ))(e vt -1), (141)
with h(r) = e -νr for L exponentially decaying or finite range and h(r) = (1 + r) ν if L is algebraically decaying with α > 2D + 1 with ν < α -(2D + 1), where D is the dimension of the lattice.

As stated above, in this work, we require a slightly different formulation of the LR-bound as given in Formula (7) of the main text, namely

(e tLΛ -e tL Λr (Y ) )(O Y ) ≤ c 1 h(diam(Λ r(Y ) ))(e vt -1), (142)
which reflects directly that the dynamics of the system can already be described by a generator L Λ r(Y ) restricted to a region Λ r(Y ) of diameter r around the initial support Y of the observable O Y . To convert a bound of the form (137), we follow the reasoning given in [10,19]. We can express the difference of the dynamics generated by the full Lindblad generator L Λ as compared to a restriction L Λr to the subset Λ r ⊂ Λ according to

(e tLΛ -e tLΛ r )O Y = - t 0 ds ∂ s e sLΛ r e (t-s)LΛ O Y = t 0 ds e sLΛ r (L Λ -L Λr ) e (t-s)LΛ (O Y ) (143)
Taking norms on both sides, we therefore obtain an upper bound of the form

(e tLΛ -e tLΛ r )O Y ≤ X̸ ⊂Λr t 0 ds L X e (t-s)LΛ (O Y ) . ( 144
)
We notice, that the term inside the integral is exactly of the form of the left-hand side of (137) with K X = L X . Hence, we can insert the standard LR-bound for dissipative dynamics from (137) here and are left with a combinatorial problem to obtain a bound on the approximation. This can be done explicitly for several standard interaction decays, such as finite range, exponentially decaying or algebraically decaying interactions [10,19]. In particular, based on the Lieb-Robinson bound in Thm. 7.1, we obtain Lemma 7.1 ([19]). Let L be a Lindbladian of the form (138), O Y an observable supported on Y ⊂ Λ and r > 0 hold. Then for L Λr(y) = X⊂Λr(y) Λ X we have

e tL -e tL Λr (y) O Y ≤ O Y |Y |J e vt -1 -vt v h(r) (145)
with h(r) exponentially decaying in r for L finite range or exponentially decaying and h(r) decaying as

(1 + r) -β if L is algebraically decaying with α > 2D + 1 and β = α -3D for α ≥ 5D -1 and β = 1 2 (α -D -1) if α ≤ 5D -1
, where again D is the dimension of the lattice.

We remark that all these bound give us the required independence of the right-hand side from the overall system size. We expect that with the help of recent more stringent estimates on Hamiltonians with algebraic decay, it will most likely be possible to extend and strengthen these bounds for other algebraic decays.

A. Explict Lieb-Robinson velocities

As discussed in section 10, making a priori estimates about the sampling complexity beyond the principle scaling, requires explicit knowledge about the constants in the LR-bounds we are using, in order to obtain a bound on the degree of the involved polynomials. This includes in particular the Lieb-Robinson velocity that governs the decay behaviour of the equations in time. Hence, within this section, we will discuss some models were this dependence can be made more explicit.

First of all, we note that for a given anticipated interaction graph and assuming a bound on the interaction strength and decay behaviour as defined in (139) and (140), according to [19] the Lieb-Robinson velocity can always be chosen according to

2J sup

x,yΛ n≥0

|b x (n)| |b x (n -1)| r≥n µ(r)|b y (r)|h(r) ≤ v, ( 146
)
where b x (n) denotes the ball of size n around lattice side x and the function h(r) is specified in Thm. 7.1 depending on the decay of the interactions. This will in principle allow us to obtain an explicit bound for any anticipated model our protocol can be applied to. However, more explicit bounds can be derived in the case of stronger assumptions on the locality of the interactions. Indeed, let us consider the case of finite-range interactions specified on an interaction graph with (hyper)edge set E that satisfies the following requirements [11] 1. finite interaction strength: g Following the reasoning in section 6, this implies a required minimal starting time t 0 = d 2 exp and maximal evolution time t max = 2 + t 0 as well as sampling m = O d log(dδ -1 evolution times within this interval according to the Chebychev distribution in order to succeed with the polynomial interpolation.

= sup X∈E L X ≤ ∞ 2. finite neighbourhood: Z = max X∈E |{Y ∈ E : X ∩ Y = ∅} ≤ ∞
Finally, let us mention that by assuming even more structure on the interactions and using the method developed in [22] based on the so-called commutativity graph, even stronger bounds on the LR-velocity can be obtained for certain models. Indeed, for example for the translation invariant transverse field Ising model in two-dimensions, [22] shows that v ∼ √ g instead of the linear dependence discussed before and we would expect that this scaling will also be an upper bound in the case of non-translation invariance.

1. Draw two Pauli strings B i , S i ∈ {σ x , σ y , σ z } n uniformly at random and also a sign vector E i = {-1, +1} n uniformly at random.

2. Prepare the quantum state ρ i = ⊗ n j=1 ϕ j (S i , E i ), where ϕ j (S i , E i ) is an eigenstate of the j-th Pauli on the string S i corresponding to the eigenvalue in the j-th entry of E i .

3. Evolve ρ i by Φ. Let us now introduce some notation to explain how to postprocess the samples obtained from the protocol above. Let P a and P b be the given Pauli operators on n qubits. We say that a basis B i overlaps with P a if all qubits on which P a acts non trivially are also measured in the same basis. For example, P a = σ x ⊗ σ y ⊗ I and B i = {σ x , σ y , σ z } overlap. However, P a = σ x ⊗ σ y ⊗ I and B i = {σ x , σ z , σ z } do not overlap. Moreover, if the basis and P a overlap, we say that the measurement outcome M i overlaps positively if we measure a positive eigenstate of P a . Otherwise, we say it is negative.

We will also define a similar notion of the state overlap given the Pauli P b , the basis S i and sign E i . We say that P b and (S i , E i ) overlap positively if S i coincides with P b on all qubits it acts non trivially and the state ρ i is a positive eigenstate of P b . We say we overlap negatively if it is a negative eigenstate. We will also define ω(P a ) to be the number of qubits on which P a acts non trivially.

We can now finally introduce a random variable given the Paulis P a , P b and the data from the experiment B i , S i , E i and M i . Let us define a function in terms of the outcomes and inputs of one round of the protocol:

X a,b (B i , S i , E i , M i ) =          0 if B i does not
Before we prove that, let us discuss how these estimates on moments on X a,b suffice to obtain the claimed sample complexity. As in other works on classical shadows, it will be crucial to use the method of median of means estimator [25] to estimate the expectation value of X a,b . The method of medians of means works as follows. We take a sample of size S and divide it into K subsets of size B, i.e. S = KB. We then compute the empirical mean on each of the K subsamples. Denote them by μi , with 1 ≤ i ≤ K. We then set our estimator of the mean to be: μMoM = median(μ 1 , . . . , μK ).

(

) 153
The main property of this estimator is that we have that if the variance of X a,b is σ 2 , then:

P(|μ MoM -E(X a,b )| ≥ ϵ) ≤ e -2K
(

1 2 -σ 2 Bϵ 2 ) . ( 154
)
In particular, if we pick B = 4σ 2 ϵ -2 and K = 2 log(δ -1 ), then:

P(|μ MoM -E(X a,b )| ≥ ϵ) ≤ δ. ( 155
)
We see that the median of means method allows to have a logarithmic scaling of the error probability from a bound on the variance. On the other hand, just using the empirical mean directly combined with an estimate on the variance gives a polynomial dependence only.

Armed with these facts about the median of means estimator, the following result immediately follows from (151) and (152):

1. perturbations of the locality of the interaction.

2. SPAM errors.

unknown structure of interactions.

By unknown structure of interactions we mean the scenario where we do not know exactly which qubits are actually coupled to each other, i.e. the precise form of the generator is not given, but we assume such a local structure exists. And by approximately local generators we mean the case of Lindbladians that are of the form L = L ′ + ∆, where L is a Lindbladian that satisfies our locality assumptions and ∆ is an arbitrary perturbation. Note that although we could extend all of our results to dynamics that have k-body interactions, we will here restrict to at most 2-body interactions for simplicity.

But before discussing more sophisticated bounds, let us make the general observation that the polynomial interpolation method we are using is robust to noise in the data. This is both in the sense of there being outliers in the data and the estimates we obtain for the points being close to that of a polynomial.

Thus, we have the following general robustness statement: Proof. By Eq. (175), measuring the expectation value of P under the time evolution will yield values that are within O(ϵ/ log(ϵ -1 )) of the expectation value for the evolution under L. As the polynomial interpolation method we use is robust to noise on the values of the polynomial and we only need the expectation value of local Paulis for our recovery, the statement follows.

Thus, we see that, in general, as long as the perturbation to the idealized, local Lindbladian is of the same order as the estimate's precision required at each step for local Paulis, we can still recover the parameters of the idealized Lindbladian.

A. Robustness to perturbations

Our results strongly rely on the fact that local quantum dynamics satisfy LR bounds. Thus, it is natural to what extend our protocol is robust against global perturbations, i.e. what happens if we have small couplings between distant parts of the system. Let us now show that our protocol is robust against such perturbations.

To model such perturbations, we will assume that the true generator is of the form

L = L ′ + ∆, (179)
where L ′ is a time evolution that satisfies a LR bound and ∆ is an arbitrary superoperator that is supposed to model the non-local parts of the evolution. Furthermore, we will define the following expansion quantity κ(∆, k) that measures the impact of the non-local part on observables supported on at most k sites:

κ(∆, k) = sup O:supp(O)≤k,∥O∥≤1 ∆(O) . ( 180
)
To see why this quantity measures the impact of global perturbations on local observables, let us consider the case where we have

∆(X) = τ n i,j=1,i̸ =j i[H i,j , X], (181)
for some arbitrary Hamiltonian evolution terms H i,j ≤ 1 supported on qubits i, j. That is, this additional term models the situation where we have a small interactions between all pairs of qubits. Note that the norm of this perturbation scales like ∼ n 2 in general. In contrast, κ(∆, k) ∼ τ kn. This because for an observable supported on k sites, only the terms that also interact with that site will act nontrivially on it. So we have a total of ∼ nk terms acting on that site. We then have:

Theorem 9.1. Let L be a Lindbladian and suppose it is of the form

L = L ′ + ∆, ( 182
)
where L ′ satisfies a LR bound as in Eq. (142). ds( e (t-s)L ∆O r s,δ + e (t-s)L ∆(δ Õs ) ).

We can estimate the first summand of the integral by ) and e (t-s)L is a quantum channel in the Heisenberg picture and cannot increase the operator norm. We can estimate the second summand in a similar way and get: As δ > 0 here was arbitrary, we can further optimize the bound over all δ > 0, which gives the claim.

In a nutshell, the bound above says that our protocol is only sensitive to the extent at which global couplings affect local observables, which makes it significantly more robust than one would expect naively. Let us illustrate the scaling of the bound above with the example above of all-to-all couplings in Eq. (181). By a triangle inequality, ∆ ≤ τ n(n-1)

2

. Thus, for an observable initially supported on a constant number of sites, we have that: κ(∆, r(s, δ)) = τ n(vs + C log(δ -1 )) D , (189)

d

  dt tr [ρ j O i (t)] = tr [ρ j L(O i )] for all elements in the set {(ρ j , O i )} in terms of the expansion coefficients.

  FIG. 2.

  FIG. 3.

FIG. 4 .

 4 FIG. 4.

Lemma 4 . 1 (

 41 Derivatives of local evolutions). Let L Λ be a locally bounded Lindbladian with constant g. For an observable O such that O ≤ 1, an initial state ρ and a region B ⊂ Λ defines f B : R + → R as f B (t) = tr e tL B (O)ρ .

Lemma 5 . 2 (

 52 Extrapolating the derivative at 0). Let p : [0, b] → R be a polynomial of degree d such that for some ϵ > 0 and 0 < a < b: max x∈[a,b] |p(x)| ≤ ϵ.

Corollary 5 . 1 .

 51 max x∈[a,b] |p(x) -p(x)| ≤ 3σ. (106) Applying Lemma 5.2 to p -p yields the claim. This then yields a simple condition on how small σ has to be in the regime of interest to us: In the same setting as in Prop. 5.1 for some ϵ > 0 let a ≤ d -2 , b = 2 + a and σ = ϵd -2 . (107) Then |p ′ (0) -p′ (0)| = O(ϵ). (108)

Definition 6 . 1 (

 61 Chebyshev partitions). Let m ∈ N be given. The size m Chebyshev partitions of [-1, 1] is the set of intervals I j = cos πj m , cos π(j-1) m for 1 ≤ j ≤ m.

Definition 6 . 3

 63 (ℓ ∞ regression). Given a set of n points (x i , y i ) and m ∈ N given. For the m Chebyshev partitions I j , choose xj ∈ I j arbitrarily and let y j = median xi∈Ij y i .

  ) as long as m = O(d) and the we pick n = O(d log(d)) samples from the Chebyshev measure, as shown in [12, Lemma 1.3]. Thus, the procedure gives us a promise of recovery in ∞-norm up to the unwanted p ∞ term. This can be solved by iterating the ℓ ∞ regression step.

3 .

 3 bounded growth of neighbours: max X∈E |{Y : dist(X, Y ) = n}| ≤ M n η For such a systems, [11, Thm.2] derives a Lieb-Robinson bound according to Lemma 7.1 with the right-hand side of (145) given by2M Z O Y r η e (egZ)•t-r = 2M Z e v•t h(r) (147)with h(r) = r η e -r and the LR-velocity v = egZ. Hence, following the reasoning in section 4, we choose r ash -1 ( Zε 4M ∥O Y ∥ exp(v•tmax) ), which for O Y ≤ 1 and for r large enough such that (1 -η log r r ) > 2 3 means that r should be chosen as 3 2 log 4M e (egZ)•tmax Zε . This in turn results for the truncated Taylor-expansion for the corresponding generator L Λy(r) restricted to the region Λ y (r) in a polynomial degree of d = 2et max g|Λ y (r)| -1, which in the case of a D-dimensional lattice results in d exp = 2et max g 3 2 log 4M e (egZ)•tmax Zε + diam(supp(O Y )) D log(ε -1 ) -1 . (148)

4 .

 4 Measure in the Pauli basis defined by B i . Denote the measurement outcome by M i . 5. Output (B i , S i , E i , M i ).

  overlap with P a or S i does not overlap with P b ,3 ω(Pa)+ω(P b ) /2if (B i , S i ) overlap with (P a , P b ), and both do so positively,3 ω(Pa)+ω(P b ) /2 if (B i , S i ) overlapwith (P a , P b ), and both do so negatively, -3 ω(Pa)+ω(P b ) /2 if (B i , S i ) overlap with (P a , P b ), one positively, the other negatively. . (150) We will now show that: E(X a,b ) = 2 -n tr (P a Φ(P b )). (151) and E(X 2 a,b ) ≤ 3 w(Pa)+w(P b )

  t-s)L ∆(O r s,δ ) ≤ t 0 dsκ(∆, |B(supp(O), log(δ -1 ) + vs)|) (187)as, by definition, ∆(O r s,δ ) ≤ κ(∆, r(s, δ)

  t-s)L ∆(δ Õs ) ≤ δt ∆ +∞→+∞ . (188)

  VIII. ACKNOWLEDGEMENTS D.S.F. was supported by VILLUM FONDEN via the QMATH Centre of Excellence under Grant No. 10059 and the European Research Council (Grant agreement No. 818761). This work is part of HQI initiative (www.hqi.fr) and DSF was supported by France 2030 under the French National Research Agency award number "ANR-22-PNCQ-0002". A.H.W. thanks the VILLUM FONDEN for its support with a Villum Young Investigator Grant (Grant No. 25452). J.B. acknowledges support from The AWS Quantum Discovery Fund at the Harvard Quantum Initiative and funding from the NWO Gravitation Program Quantum Software Consortium. V.V.D. work is a part of the research program NWO QuTech Physics Funding (QTECH, programme 172) with project number 16QTECH02, which is (partly) financed by the Dutch Research Council (NWO); the work was partially supported by the Kavli Institute of Nanoscience Delft. L.A.M. was supported by the Netherlands Organisation for Scientific Research (NWO/OCW), as part of the Quantum Software Consortium program (project number 024.003.037 / 3368).

  Lemma 9.1. Let L be a Lindbladian on a D-dimensional regular lattice and ϵ > 0 be given. Suppose we can measure the expectation value of two-body Pauli observables in the time interval [t 0 , t max ] under a family of quantum channels {Φ t } t0≤t≤tmax s.t. for all t ∈ [t 0 , t max ] we have for 2-body Paulis P :

	(Φ *	(175)
	for
	t -1 0 = O polylog(ϵ -1 )	(176)

t -e tL )(P ) ≤ O(ϵ/ log(ϵ -1 )). and t max = 2 + t 0 . Assume further that we can prepare Pauli eigenstates and perform Pauli measurements up to an error O(ϵ/ log(ϵ -1 )). Then, measuring such Pauli measurements on Pauli eigenstates under the evolution of Φ * t for m = O polylog(ϵ -1 ) (177) random times up to precision O(ϵ/ polylog(ϵ -1 )), is sufficient to obtain with probability of success at least 2/3 an estimate âα of a α satisfying |â α -a α | = ϵ. (178)

  Then we have for an observables O s.t. O ≤ 1 and e tL -e tL ′ (O) ≤ inf

δ>0 tδ ∆ +∞→+∞ + t 0 dsκ(∆, |B(supp(O), log(δ -1 ) + vs))|). (

183

)

Proof. We have that:

(e tL -e tL ′ )(O) = t 0 dse (t-s)L ∆e sL ′ (O).

(184)

Now for any δ > 0 we can approximate e sL ′ (O) up to δ in operator norm by an observable supported on a ball of radius r s,δ = vs + log(δ -1 ) around the original support of O, i.e.:

e sL ′ (O) = O rs + δ Õs

(185)

where O rs , Õs ≤ 1 and Õs does not necessarily have any further locality. Now:

(e tL -e tL ′ )(O) ≤ t 0 ds e (t-s)L ∆e sL ′ (O) ≤ t 0

Supplemental Material to Efficient and robust estimation of many-qubit Hamiltonians

Daniel Stilck França, 1, 2, * Liubov A. Markovich, 3, 4 V. V. Dobrovitski, 3 Albert H. Werner, 1, 5 and Johannes Borregaard 3

Supplementary Note 7: Lieb-Robinson bounds

This section gives a brief overview of Lieb-Robinson bounds. In particular, we give more explicit formulas for the functions h (see Formula (13) in the main text) in terms of the decay of the interactions and the dimension of the lattice. Lieb-Robinson bounds are by now a standard tool in quantum many-body systems and quantum information theory and we refer to [8,9,11,[16][17][18][19][20] for a more general overview over the mathematical background and some latest bounds for algebraically decaying interactions.

At the heart of any Lieb-Robinson bound is the intuitive idea that if interactions in a system happen locally this should imply a bound on how fast information can be transmitted. The usual way to codify this property for Hamiltonian systems in the Heisenberg picture is to give a bound on the operator norm of the commutator between a time-evolved observable O Y (t) initially located in region Y and a second observable A X located in a region X in the distance dist(X, Y ) between the regions X, Y , i.e. a bound of the form

where C will typically depend on the size of the regions |X| and |Y | as well as on the operator norms of A and O(t). However, in the context of Markovian dynamics and master equations, the bound is usually generalized by substituting the super-operator [A X , •] for an arbitrary bounded super-operator K X : M 2 n → M 2 n supported on X leading to a Lieb-Robinson-bound of the form

with C depending on K X ∞→∞,cb .However, if K X is of the form K X = [A X , •], we have K X ∞→∞,cb ≤ 2 A X ∞ , which allows us to recover the commutator [10,19]. In the following, we consider a regular lattice Λ and assume that the dynamics is generated by a Lindbladian that decomposes according to

Following [9,19,21], we define the maximal interaction strength

as well as the decay behaviour of the interactions

Supplementary Note 8: Parallelizing the Measurements: shadow process tomography

In this section we introduce a method to parallelize the estimation of the Pauli overlaps required for our protocol. For a quantum channel on n qubits Φ and two collections of Pauli strings P 1 a , . . . , P K1 a and P 1 b , . . . , P K2 b that have combined weight (i.e. number of Paulis that are nonidentity) ω a +ω b at most ω it allows us to estimate all the overlaps of the form

up to an error ϵ with probability at least 1 -δ from a number of samples that grows like O(3

Moreover, it only requires us to prepare simple Pauli eigenstates and measure in Pauli eigenbases. Recently two works [23,24] have considered how to generalize the shadows protocol to the setting of process tomography. Here we give an alternative analysis that has a better sample complexity than the previous works when specialized to Pauli strings. Thus, we believe that this section may be of independent interest. The protocol we will introduce now allows us to estimate all the data required for our Hamiltonian learning protocol in parallel. Every round i of the protocol is performed as follows:

Corollary 8.1. Let {P l a } K1 l=1 and {P j b } K2 j=1 be two collections of Pauli matrices on n qubits such that for any l, j the following condition

holds. Then

runs of the protocol above suffice to obtain an estimate e l,j a,b satisfying

for all pairs l, j with probability at least 1 -δ.

Proof. Let X l,j a,b be the random variable we have defined above for a pair of Paulis P l a and P j b . By the bound in (152) we have that E((X l,j a,b ) 2 ) ≤ 3 ω , holds. Thus, for the median of means estimator with B = 3 ω ϵ -2 samples per group and K = 2 log(δ -1 K 1 K 2 ) we have an estimate satisfying (158) with probability of failure at most δ/K 1 K 2 .

By a union bound, the median of means estimator of all K 1 K 2 combinations of P l a and P j b satisfy (158) with probability of failure at most δ. The total number of samples required for this is

which yields the claim.

To conclude we only need to show that (151) and (152) hold. Let us start with the expectation value. First, note that

holds. Here we say that B and S overlap if B overlaps with P a and S overlaps with P b . The latter formula holds, because if we do not have overlap then X a,b = 0 by the case 1 of the definition (150). Now observe that P(B and S overlap) = 3 -ω(Pa) 3 -ω(P b ) .

(161)

This holds because we have a 1/3 chance of "hitting the right Pauli" at each point of the support of either P a or P b and they are all independent. Thus, all we need is to determine

Here we distinguish four cases: We can then break down the expectation of (160) down into the cases:

where c is a random variable which keeps track of which case we have. By definition, if c = 1, 2, then

Note that the 3 ω(Pa)+ω(P b ) counteracts the P(B and S overlap) term, up to an additional 1/2 factor.

Thus, all that is left is to estimate

Let us estimate P(c = 1), the other cases will be analogous. To this end, we will introduce Q + a and Q - a , which are the projectors onto the positive and negative eigenvalues, respectively, of the support of P a . We will use analogous notation for

Here and in what follows we will always assume for simplicity that the Pauli strings are supported on the first qubits.

Let us estimate the expected initial state for the case c = 1. Conditioned on being a state with positive overlap with P b , by construction we know that the state is uniformly distributed on the positive eigenspace of P b . This corresponds to the state

Given that this is the initial state and that we are measuring in the eigenbasis of P a (recall that we have overlap), the probability of measuring a positive outcome is tr (Q + a ⊗ I n-ω(Pa) Φ(σ)). We conclude that

Similarly

As

we conclude that

Computing the second moment of X a,b turns out to be quite simple. Note that the only nonzero value the random variable X 2 (a,b) takes is 3 2ω(Pa)+2ω(P b ) /4 with probability 3 -ω(Pa)-ω(P b ) . Thus, we clearly have that:

which yields (152). This concludes all the computations required for Cor. (8.1). Note that the protocol only requires the preparation of Pauli states and Pauli measurements. Thus, it should be feasible to implement it on near term devices. Additionally, the required postprocessing is efficient, as evaluating the value of X a,b can be efficiently done given a sample.

Furthermore, note that if we further have the information that we do not to wish to recover certain bases (i.e. we do not wish to recover Pauli strings with Y terms), it is possible to adapt the protocol and do not prepare initial states or measure in that basis. This will reduce the sample complexity accordingly.

Supplementary Note 9: Robustness of the protocol to various errors

It is natural to inquire to what extent our protocol is robust against state preparation and measurement errors (SPAM) and to what extent it still outputs reliable estimates whenever the underlying assumptions are only valid approximately. The goal of this subsection is to address all of these concerns and show that our protocol is highly robust.

In the following subsections we will discuss the following forms of robustness:

From this we conclude that for such all-to-all perturbations, we have that

where for the last bound we picked δ = n -1 . Thus, we see that as long as ntτ ϵ, our protocol still produces reliable results even under such global perturbations. Also note that our more refined perturbation bound using the LR-bound is quadratically better than the naive bound n 2 tτ ϵ we directly obtain from observing that ∆ +∞→+∞ ≤ n 2 .

B. Robustness to characterized SPAM errors

Let us consider how to adapt our protocol to the case where we have SPAM errors. We will model such errors by assuming we know quantum channels T P and T M that model the noise when preparing the state and the measurement and we assume that these channels are independent of the Pauli basis we are measuring or preparing.

In that case, the process shadow tomography would be estimating the value of

If the observable T * M (P a ) is sufficiently local, i.e. it does not increase the support of P a signigicantly, then our previous results hold and we know that the function f : t → 2 -n tr (T * M (P a ))e tL (T P (P b )) will be well-approximated by a polynomial of low-degree. Thus, we can see that we can evaluate the derivative of the function at 0 exactly as before. However, now will get access to the equation:

If the channels T M , T P are sufficiently simple (i.e. approximately product), then it is possible to expand the expectation value in Eq. ( 192) efficiently and obtain a linear equation for the parameters of L. Note, however, that these SPAM errors might decrease the condition number of the system of equations. Thus, it would be necessary to see how the condition number changes with the noise and then adjust the precision with which we estimate the derivative accordingly.

Let us exemplify this in a simple form of SPAM: local depolarizing noise. Let us assume that both T M and T P consist of D ⊗n p , where for some parameter p D p (ρ) = pρ + (1 -p) I 2 . In this case, we see that T P (P b ) = p w(P b ) P b , T * M (P a ) = p w(Pa) P a . Thus, in this case we see that the expectation values are rescaled by a factor p w(Pa)+w(P b ) and we can run exactly the same protocol, but now with a sampling overhead of p 2k to learn k-local Hamiltonians. In conclusion, we see that for noise models like local depolarizing (or in fact any Pauli noise), it is easy to recover our original results with a mild sampling overhead by adjusting the linear systems.

Developing an analytical understanding of how the condition number changes for other important noise models like amplitude damping is beyond the scope of this work. But we once again emphasize that this is straightforward and efficient to do numerically.

C. Robustness to unknown SPAM errors

In the last subsection we discussed how to deal with SPAM errors if they are well-characterized and independent of the measurement and state being prepared. Now we are going to discuss the case in which they are unknown and can potentially depend on the measurement and state we wish to prepare. We will show that the effect of SPAM errors will be local, in the sense that as long as the local SPAM rate is of order at most ϵ, then the protocol still yields reliable results. For simplicity, we will consider SPAM errors that are modelled by single-qubit product quantum channels. However, our proof works analogously for channels with an approximately finite light-cone. That is, they map local observables to observables that are approximately local. More formally, in the SPAM-free setting, our protocol entails measuring states of the form:

where U i , V i are single-qubit Pauli rotations. We then estimate local observables on such states. We will now model that instead we prepare the state

where T i,Ui and R i,Vi are single-qubit quantum channels that model the SPAM errors and can potentially depend on the gates we are implementing. We will now see that only local SPAM errors affect the performance of our protocol. But before that, let us discuss the sensitivity of the shadows protocol of Sec. 8 to noise. Suppose we wish to use the shadow protocol to estimate a pair of Paulis P a , P b of respective weight w a , w b , i.e. we wish to estimate

Then the protocol discussed in Sec. 8 consists of measuring an observable of O a,b,{Vi},{Ui} that depends on the state preparation and measurement unitaries, is supported only on the support of P b and such that O a,b,{Vi},{Ui} = 3 wa+w b /4 over an ensemble of random states. Denoting by B the support of P b , we conclude that if for all possible U i , V i we have

then by an application of Hölder's inequality, we have that the expectation value of the shadows protocol will deviate by at most τ 3 wa+w b /4 from the ideal value. Although this exponential scaling of the sensitivy in the locality is undesirable, if the SPAM errors are very strong, we can also adopt the less efficient strategy of preparing each state for the isolation separately. As for that case we only need to measure a local observable of operator norm 1, we see that we avoid the exponential factor and the overall error in the estimation for the observable directly follows from Eq. ( 196) and corresponds to τ .

Let us now show that τ does not scale with the system's size and only depends on the local noise rate: Proposition 9.1. In the notations above, assume that for all i, U i , V i we have that R i,Vi -id ⋄ , T i,Vi -id ⋄ ≤ τ . Furthermore, assume that L satisfies a LR-bound like that of Eq. (81). Then for t, w a , w b = O(1) we have that:

In particular, under such a SPAM error we have that the process shadow protocol for Paulis P a , P b will have an expectation value that is O(3 wa+w b poly(w a , w b )τ polylog(τ -1 )) close to the SPAM-free setting.

Proof. the statement is equivalent to having for all

First, note that we have that

where |ϵ| ≤ w b τ and O ′ B is again an observable of operator norm at most 1 on B. Here we have used that all the unitaries and channels acting outside of B will have a trivial effect. For those acting on B we use our bound on the diamond distance to the identity channel. Now, by our LR bound we have that for times t = O(1), we can write e tL * (O ′ B ) = O t,B(t,r) + X r,t , where O t,B(t,r) is an observable supported on a ball of size tv + r around B and X r,t ≤ O(e-Cr). In particular, picking r = O(log(τ -1 )) and as t = O(1) we get that e tL * (O ′ B ) is localized on a number of sites that is at most w b (vt + C log(τ -1 )) D up to an error ϵ, where D is the dimension of the lattice.

Thus:

where ϵ 2 ≤ τ by our choice of r. We can then follow the same procedure as before for the T i,Vi to conclude that only those channels on the support of O ′ t,B(t,r) contribute to the noise. Thus, we get an additional error of order O(τ polylog(τ -1 )). This concludes the proof.

Recall that for the protocols discussed here we have w(P a ), w(P b ) = O (1). Thus, we see that our protocol is robust to local, gate-dependent SPAM noise as long as the local noise rates are smaller than the target accuracy and the SPAM errors do not need to decrease with the system's size.

Supplementary Note 10: Learning the structure of the interactions

In the main text we assumed that we know which qubits interact with which qubits, i.e. the interaction graph is known, and we only wish to recover the parameters of this interaction. However, a close inspection of the proof shows that we do not really require exact knowledge of the structure; all we need is the promise that the evolution satisfies a LR bound. Indeed, the only structure required for our proofs to go through is that the evolution of local observables can be well-approximated by that of truncated Lindbladian, and this remains true if the structure of the interactions is unknown, but the true evolution still satisfies a LR bound.

Thus, in this section we will discuss how to apply our protocol in the setting in which we expect the evolution to satisfy a LR bound, but we do not necessarily know all qubits that are coupled to each other. Discussing this setting in full detail and generality goes beyond the scope of this paper, so we will restrict to the case where there is no Lindblad term and we restrict to interaction acting on only two qubits.

Furthermore, we will make the assumption that we have a lower-bound on the strength of all the interactions that are non-zero. This is necessary to obtain bounds on the performance of the protocol, as if we don't have such a lower-bound it will in principle be necessary to measure expectation values with a very high precision to see their effect.

We will thus assume that our generator is of the form

where the (unknown) interaction graph G = (V, E) is such that the evolution induced by L satisfies a LR bound as in Eq. ( 81), Proof. The protocol is as follows: we first perform the process shadow tomography protocol for m random times t i and measure the output O(η -2 log(n)/ polylog(η -1 )) times for each t i . Then, by a union bound, we will obtain the expectation value of all the O(n 4 ) pairs input and output 2-body Pauli strings up to precision η/ polylog(η -1 ) each with probability at least 2/3. Conditioning on the event that all of the estimates were correct, we can use our previous results to estimate the interaction strength for all 2-body Pauli terms up to precision η/4. We then set an edge between all qubits such that they have a 2-body interaction at least 3 4 η in absolute value. By our assumption that all the interactions are at least η, we will then identify all interactions correctly. Thus, we see that our protocol can also correctly identify the interactions as long as some LR-bound is given. To the best of our knowledge, this is the first protocol that does not necessarily require knowledge of the graph, just some assumptions on the locality. Furthermore, note that the sample complexity of this version of the protocol requires approximately the same number of samples as to estimate the values themselves. However, note that learning the graph is computationally more expensive. This is because we need to run O(n 2 ) polynomial interpolations to estimate all couplings, whereas when we are given the graph, the number of interpolation steps required is O(n).