# Efficient and robust estimation of many-qubit Hamiltonians

Daniel Stilck Franca, Liubov A. Markovich, V. V. Dobrovitski, Albert H. Werner, Johannes Borregaard

► To cite this version:

Daniel Stilck Franca, Liubov A. Markovich, V. V. Dobrovitski, Albert H. Werner, Johannes Borregaard. Efficient and robust estimation of many-qubit Hamiltonians. Nature Communications, 2024, 15 (311), 10.1038/s41467-023-44012-5. [https://doi.org/10.1038/s41467-023-44012-5](https://doi.org/10.1038/s41467-023-44012-5)

HAL Id: hal-03675783

https://hal.science/hal-03675783v2

Submitted on 10 Jan 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
Efficient and robust estimation of many-qubit Hamiltonians

Daniel Stilck França,¹,²,* Liubov A. Markovich,³,⁴ V. V.

Dobrovitski,³ Albert H. Werner,¹,⁵ and Johannes Borregaard³,⁶

¹QMATH, Department of Mathematical Sciences, University of Copenhagen,
Universitetsparken 5, 2100 Copenhagen, Denmark

²Univ Lyon, ENS Lyon, UCBL, CNRS, Inria, LIP, F-69342, Lyon Cedex 07, France

³QuTech and Kavli Institute of Nanoscience, Delft University of Technology, 2628 CJ, Delft, The Netherlands

⁴Instituut-Lorentz, Universiteit Leiden, P.O. Box 9506, 2300 RA Leiden, The Netherlands

⁵NBIA, Niels Bohr Institute, University of Copenhagen, Blegdamsvej 17, 2100 Copenhagen, Denmark

⁶Department of Physics, Harvard University, Cambridge, Massachusetts 02138, USA

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–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 [7–35].

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, 26, 29, 36, 37], which may be difficult for realistic devices subject to various decoherence mechanisms.

Alternatively, several works [30–32] 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

<figure><img src="image_2.png" /><figcaption>FIG. 1.</figcaption></figure>

rates (jointly referred to as the parameters of the Lindbladian $\mathcal{L}$) governing the system. Furthermore, for finite range interactions, these approaches can estimate the parameters of the Lindbladian to a given precision from a number of samples that is independent of the system's size [30–32].

A significant drawback of these approaches is that

<footnote>* daniel.stilck_franca@ens-lyon.fr</footnote>
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 $\epsilon$, the system has to be probed at times $O(\epsilon)$ apart and expectation values of observables have to be estimated up to a precision of $O(\epsilon^2)$, which translates to an overall $O(\epsilon^{-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(\text{polylog}(\epsilon^{-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 $\epsilon$ of $O(9^k \epsilon^{-2} \text{polylog}(n, \epsilon^{-1}))$. We obtain this by estimating time derivatives using multiple temporal sampling points and robust polynomial interpolation [38]. Furthermore, we show how to use shadow process tomography methods to estimate multiple parameters in parallel. In particular, we improve the results of Refs. [39, 40] 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, $\mathcal{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. $\frac{d}{dt}O(t) = \mathcal{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 [41]:

$$ \mathcal{L}(O) = \sum a_i [P_i, \rho] + \sum D_{i,j} (P_i^\dagger O P_j - \frac{1}{2} \{P_j^\dagger P_i, O\}) \quad (1) $$

Such an expansion is always possible since this basis amounts to a Hilbert-Schmidt orthogonal set of Hermitean 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 $\mathcal{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 $\rho$ of the system gives us a linear equation for the expansion coefficients [30–32]. 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 $\mathcal{L}$ in an operator basis, the following steps are performed.

1. Find a suitable complete set, $\{(\rho_j, O_i)\}$ of multi-qubit product states $(\rho_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$.

2. Calculate the expectation values appearing on the right-hand side of the Master equations $\frac{d}{dt}\text{tr}[\rho_j O_i(t)] = \text{tr}[\rho_j \mathcal{L}(O_i)]$ for all elements in the set $\{(\rho_j, O_i)\}$ in terms of the expansion coefficients. Since the initial states and the observables are products, this can be done efficiently.

**Quantum processing:** In order to solve for the expansion coefficients, we also need the values of the time-derivatives appearing on the left-hand side of the Master equations, i.e. $\frac{d}{dt}\text{tr}[\rho_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 $\rho_j$ and evolved for a time $t_k \in \{t_0, t_1, \dots, t_{\max}\}$ after which the observable $O_i$ is measured.

2. The above procedure is repeated for each element in the set $\{(\rho_j, O_i)\}$ for all evolution times $t_k$ to obtain estimates of $\langle O_i(t_k) \rangle_j = \text{tr}[\rho_j(t_k)O_i]$ where $\rho_j(t_k)$ is the state of the system having evolved for time $t_k$ from the initial state $\rho_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 $\{(\rho_j, O_i)\}$ since the expectation values $\langle O_i(t_k) \rangle_j$ are estimated sequentially. However, we also propose a variation of the classical shadows protocol of Ref. [42] for process tomography that can reduce this to a logarithmic scaling. In essence, we can obtain estimates of all elements in the set $\{(\rho_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^{w_1+w_2}\epsilon^{-2}\log(K))$ times set by the required precision, $\epsilon$, of the estimates, the size, $K$, of the set $\{(\rho_j, O_i)\}$ and the weights (i.e. maximum number of sites differing from identity) of $\rho_j$ ($w_1$) and $O_i$ ($w_2$). The whole procedure is repeated for all times $t_k \in \{t_0, t_1, \dots, t_{\max}\}$.

From the measurement statistics of the above procedure, it is possible to obtain accurate estimates of all $\langle O_i(t_k) \rangle_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 $\{(\rho_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 $\frac{d}{dt}\text{tr}[\rho_j O_i(t)]$ from the experimentally obtained time trace of $\langle O_i(t_k) \rangle_j$ and solving for the expansion coefficients $\{a_\alpha\}$. This involves

1. Fit the time trace of $\langle O_i(t_k) \rangle_j$ with a low-degree polynomial in the time, $p_{i,j}(t)$ and estimate $\frac{d}{dt}\text{tr}[\rho_j O_i]$ as $\frac{d}{dt}p_{i,j}(t)|_{t=0}$. This is done for each element in the set $\{(\rho_j, O_i)\}$.

2. Solve the set of linear equations from the Master equations $\frac{d}{dt}\text{tr}[\rho_j O_i] = \text{tr}[\rho_j \mathcal{L}(O_i)]$ with respect to the expansion coefficients. This is possible since $\frac{d}{dt}\text{tr}[\rho_j O_i]$ has been estimated from the polynomial fits and all expectation values appearing in $\text{tr}[\rho_j \mathcal{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 $\{(\rho_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 $\mathcal{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 $\frac{d}{dt}\text{tr}[\rho_j O_i(t)] = \text{tr}[\rho_j \mathcal{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

$$ \rho_{k,l}^{(i,j)} = \frac{(I + \sigma_k^{(i)})}{2} \otimes \frac{(I + \sigma_l^{(j)})}{2} \otimes \rho_{n-2}, \quad (2) $$

where $k, l = x, y, z$. Thus, the i'th and j'th qubit are prepared in eigenstates of the Pauli matrices $\sigma_k$ and $\sigma_l$ while the state of the remaining $n-2$ qubits, $\rho_{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

$$ \sigma_k \sigma_l = \delta_{kl} I + i \epsilon_{klp} \sigma_p, \quad (3) $$

where $\delta_{kl}$ is the Kroenecker delta function and $\epsilon_{klp}$ is the Levi-Civita symbol. As we show explicitly in the Supplementary Note, if a single qubit Pauli observable, $O = \sigma_l^{(i)}$, 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 $(\sigma_i^{(i)} \otimes \sigma_j^{(j)})$
$\sigma_j^{(j)}$) can be isolated by choosing observables of the form
$O = \sigma_i^{(i)} \otimes \sigma_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 iso-
lated for a general Hamiltonian with terms coupling be-
tween two and *k* qubits and arbitrary single qubit dissi-
pation terms. We note that already for two qubit dissipa-
tion 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 $\langle O_i(t_k) \rangle_j$ for a set of times $t_k$ we want to extract a polynomial $p_{i,j}(t)$ such that we can estimate $\frac{d}{dt}\text{tr}[\rho_j O_i]$ as $\frac{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 $\langle O_i(t_k) \rangle_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 low-
degree polynomials and then showing how to robustly
extract the derivative of the polynomial from this infor-
mation. Before we give an overview of the ideas, we will
also need to introduce some notation to deal with Lind-
bladians acting on different parts of the system. We will
denote by $\Lambda$ the whole system of qubits and for some
subset $B \subset \Lambda$ of qubits, we denote by $\mathcal{L}_B$, the generator
truncated to those qubits.

The first step of our proof is to establish that the expec-
tation value $\langle O(t) \rangle$ of a local observable $O$ that evolves
under a Lindbladian $\mathcal{L}_B$ restricted to some sub-region $B$
up to some time $t_{max}$, can indeed be approximated up
to error $\varepsilon$ by a degree-$d$ polynomial, where $d$ depends

linearly on the size of $B$, $t_{max}$ and $\log(\varepsilon^{-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_a^m$, that we identified in the previ-
ous step, to a local generator. The main insight here is
that for finite range (or sufficiently quickly fast decay-
ing) interactions, the dynamics of any local observable
$O$ exhibits an effective light cone quantified by a Lieb-
Robinson bound (LRB) [43–48]. 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 $\varepsilon$ by a polynomial of
degree $\mathcal{O}[\text{poly}(t_{max}, \log(\varepsilon^{-1}))]$. Now, making use of
the techniques from Ref. [38], we can extract the first
derivative of this polynomial from measurements at
$\mathcal{O}[\text{polylog}(\varepsilon^{-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 pre-
pare Pauli eigenstates and measure in Pauli eigenbases
and did not consider state preparation and measurement
(SPAM) errors. Furthermore, we assumed that the dy-
namics 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 er-
rors 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 with-
out changing its performance significantly. For instance,
if the SPAM errors are described by local depolarizing
noise with depolarizing probability $p$, we can still re-
cover 2-local Hamiltonians with single qubit noise with
a $\mathcal{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 sys-
tem'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 $\mathcal{L} = \mathcal{L}' + \Delta$, where $\mathcal{L}$ is a Lindbladian that satisfies a LR-bound and $\Delta$ 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\|\Delta\|_{\infty\to\infty}$ away from the unperturbed value, where $\|\Delta\|_{\infty\to\infty} = \sup_{\|X\|_{\infty}\le1} \|\Delta(X)\|_{\infty}$. However, we show that only the effect of $\Delta$ 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 $\Delta$ is a small all-to-all coupling between all pairs of qubits of order $\tau$, e.g.

$$ \Delta(X) = \tau \sum_{i,j=1, i \neq j}^{n} i[\sigma_x^{(i)} \sigma_x^{(j)}, X]. \quad (4) $$

Clearly, such an evolution does not satisfy a LR-bound and $\|\Delta\|_{\infty\to\infty} = \mathcal{O}(\tau n^2)$. However, we show that under such a perturbation, the expectation value is only perturbed by $\mathcal{O}(t\tau 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 $\eta$: we just estimate all k body-couplings up to a precision $\eta/4$ and discard all of those that we see are smaller than $\frac{3}{4}\eta$. This can be done with a number of samples that scales like $\mathcal{O}(9^k \eta^{-2} \text{polylog}(\eta^{-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 [49, 50]. 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)$ [49, 50]

$$ H_{ij} = \sum_{k=i,j} \frac{1}{2} \tilde{\omega}_k \sigma_z^{(k)} + \left[ \frac{g_i g_j}{\Delta_{ij}} + g_{ij} \right] (\sigma_+^{(i)} \sigma_-^{(j)} + \sigma_-^{(i)} \sigma_+^{(j)}) \quad (5) $$

for $i \neq j = 1, \dots, n$ and a dissipation term acting on the $i$'th qubit and having jump operators $\sigma_-^{(i)}, \sigma_+^{(i)}$ (generalised amplitude damping) and $\sigma_z^{(k)}$ (pure dephasing). Here $\tilde{\omega}_k = \omega_k + \frac{g_k^2}{\Delta_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 $\Delta_k = \omega_k - \omega_c < 0$ where $\omega_c$ ($\omega_k$) is the frequency of the coupler (k'th qubit) and have defined $1/\Delta_{ij} = (1/\Delta_i + 1/\Delta_j)/2$. By adjusting the frequencies of the coupler and the qubits, the effective qubit-qubit interaction can be canceled up to experimental precision. Typical qubit frequencies are around 5–6 GHz [1]), while $\Delta_k \sim -1$ GHz, $g_{ij} \sim 10-20$ MHz, and $g_j \sim 100$ MHz [49, 50]. 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. $\frac{g_i g_j}{\Delta_{ij}} + g_{ij} \sim 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^* \sim 150 \ \mu\text{s}$ as well as pure dephasing resulting in a transverse relaxation on a timescale $T_2 \sim 60 \ \mu\text{s}$ representing state of the art coherence times [1, 50]. Finally, we include longitudinal relaxation of the qubits through an amplitude damping channel on the time scale of $T_1 \sim 60 \ \mu\text{s}$. We refer to Supplementary Note 2 for 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 \times 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 [30–32]. Besides already
<figure><img src="image_3.png" /><figcaption>FIG. 2.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>FIG. 3.</figcaption></figure>

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 $\sigma_X\sigma_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. [38], 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-
<figure><img src="image_5.png" /><figcaption>Recovery error using interpolation (log$_{10}$)</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>FIG. 4.</figcaption></figure>

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) = \text{tr}[e^{t\mathcal{L}(O)}\rho]$ and $\mathcal{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 $\mathcal{O}(\epsilon)$ for random times $t_1, \dots, t_m$. The curve of $f(t)$ is then fitted to a low-degree polynomial $p$, and $p'(0)$ is taken as an estimate for the derivative $f'(0) = \text{tr}[\mathcal{L}(O)\rho]$. 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 $\mathcal{L}$ be a local Lindbladian on a $D$-dimensional lattice $\Lambda$. Moreover, let $t_{\text{max}}$, $\epsilon > 0$ and $O_Y$ be a 2-qubit observable, supported on some region $Y \in \Lambda$, such that $\|O_Y\| \le 1$, holds. Then there is a polynomial $p$ of degree

$$d = \mathcal{O}[\text{poly}(t_{\text{max}}, \log(\epsilon^{-1}))], \quad (6)$$

such that for all $0 \le t \le t_{\text{max}}$:

$$|\text{tr}[e^{t\mathcal{L}_\Lambda}(O_Y)\rho] - p(t)| \le \epsilon, \quad (7)$$

and $p'(0) = \text{tr}[\mathcal{L}_\Lambda(O_Y)\rho]$, holds.

The main technical tool required for the proof are Lieb-Robinson bounds (LRB) [43–48], 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 \subset B$ we have

$$\| (e^{t\mathcal{L}_B} - e^{t\mathcal{L}_\Lambda})(O_Y) \| \le c_1 \exp(-\mu \text{dist}(Y, \Lambda \setminus \{B\})) (e^{vt} - 1), \quad (8)$$

to hold for constants $c_1$, $\mu$ and $v$, where $\text{dist}()$ is the distance in the graph.

From the LRB we can show that the dynamics is well-approximated 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 $\mathcal{L}_B$ for $B$ large enough and take that as the approximating polynomial. As the derivatives of the evolution under $\mathcal{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 $\mathcal{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_{\text{max}}]$ under $\mathcal{L}$ for $t_0$ as

$$t_0^{-1} = \mathcal{O}[\text{polylog}(\epsilon^{-1})] \quad (9)$$

and $t_{\text{max}} = 2 + t_0$. Then, measuring the expectation values at

$$m = \mathcal{O}[\text{polylog}(\epsilon^{-1})] \quad (10)$$

random times up to precision $\mathcal{O}(\epsilon/\text{polylog}(\epsilon^{-1}))$, is sufficient to obtain an estimate of the Lindbladian coefficients $\hat{a}_i$ of $a_i$ satisfying

$$|\hat{a}_i - a_i| = \epsilon. \quad (11)$$
This yields a total sample complexity of $S = O(\epsilon^{-2} \text{polylog}(\epsilon^{-1}))$.

Importantly, Theorem 2 bypasses both requiring small initial times and $O(\epsilon^{-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 [38] 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 \in \overline{1, m}$ correspond to different times $\in [t_0, t_{\text{max}}]$ and the $y_i \in \mathbb{R}$ 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 $\sigma > 0$ such that

$$y_i = p(x_i) + w_i, \quad |w_i| \le \sigma, \qquad (12)$$

hold, for strictly more than half of the $y_i$. The rest might be outliers. In our setting, the magnitude of $\sigma$ corresponds to the amount of shot noise present in the estimates of the expectation values.

The authors of [38] then show that by sampling $m = O(d \log(d))$ points from the Chebyshev measure on $[t_0, t_{\text{max}}]$, a combination of $\ell_1$ and $\ell_\infty$ regression allows us to find a polynomial $\hat{p}$ of degree $d$ that satisfies:

$$\max_{x \in [t_0, t_{\text{max}}]} |p(x) - \hat{p}(x)| = O(\sigma). \qquad (13)$$

Although the details of the $\ell_1$ and $\ell_\infty$ 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 $\sum_i |p_1(x_i) - y_i|$. After finding $p_1$ we compute the polynomial $p_\infty$ that minimizes $\max_i |p_\infty(x_i) - (y_i - p_1(x_1))|$. We then output $\hat{p} = p_1 + p_\infty$ as our guess polynomial. Note that finding both $p_1$ and $p_\infty$ can be cast as linear programs and thus can be solved efficiently [51].

By combining this result with Thm. 1, we robustly extract a polynomial that approximates the curve $t \mapsto \text{tr}[e^{t\mathcal{L}}(O_Y)\rho]$ up to $O(\epsilon)$ for $t \in [t_0, t_{\text{max}}]$. Indeed, we only need to estimate the expectation value $f(t_i)$ up to $\epsilon$ for enough $t_i$ and run the polynomial interpolation.

Note that Eq. (13) only allows us to conclude that $p - \hat{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 [52], 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: \mathbb{R} \to \mathbb{R}$ with $h(x) = o(1)$, we have

$$\| (e^{t\mathcal{L}_B} - e^{t\mathcal{L}_\Lambda})(O_Y) \| \le h(\text{dist}(Y, V\setminus\{B\}))(e^{vt} - 1). \qquad (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 $\epsilon > 0$ and evolution on a $D$-dimensional lattice, assume that we pick the initial time as

$$t_0 = O \left[ \left( h^{-1} \left( \frac{\epsilon}{2(e^{2.5\epsilon} - 1)} \right)^D \log(\epsilon^{-1}) \right)^{-2} \right]. \qquad (15)$$

Furthermore, assume that we estimate the expectation value of local observables up to precision $O(\epsilon)$ at $\tilde{O}\left[ \left( h^{-1} \left( \frac{\epsilon}{2(e^{2.5\epsilon} - 1)} \right)^D \log(\epsilon^{-1}) \right)^2 \right]$ points. Then we can
estimate each parameter up to an error of

$$ \mathcal{O}\left[\epsilon\left(h^{-1}\left(\frac{\epsilon}{2(e^{2.5\nu}-1)}\right)^D \log(\epsilon^{-1})\right)^2\right], \quad (16) $$

through, the same procedure as in the local case. Note that the error in Eq. (16) only tends to 0 as $\epsilon \to 0$ if $h^{-1}\left(\frac{\epsilon}{2(e^{2.5\nu}-1)}\right)^D \log(\epsilon^{-1}) = o(\epsilon^{-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^{-\alpha}$ with $\alpha > 5D-1$ we obtain that $h^{-1}(\epsilon) = \mathcal{O}(\epsilon^{-\frac{1}{\alpha-3D}})$, holds. We summarize the resulting resources 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 [39, 40], 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 $\Phi$, Pauli strings $P_a^1, \dots, P_a^{K_1}$ that differ from the identity on at most $\omega_a$ sites and Pauli strings $P_b^1, \dots, P_b^{K_2}$ that differ from the identity on at most $\omega_b$ sites, it is possible to obtain estimates $\hat{e}_{m,l}$ of $2^{-n} \text{tr}[P_a^m \Phi(P_b^l)]$ satisfying

$$ |2^{-n} \text{tr}[P_a^m \Phi(P_b^l)] - \hat{e}_{m,l}| \le \epsilon \quad (17) $$

for all $m, l$ with probability at least $1 - \delta$ from

$$ \mathcal{O}(3^{\omega_a + \omega_b} \epsilon^{-2} \log(K_1 K_2 \delta^{-1})) \quad (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 $\rho_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 \binom{n}{k} \le 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

$$ \mathcal{O}(9^k \epsilon^{-2} k \log(n\delta^{-1})) \quad (19) $$

samples. As our protocol requires estimating expectation values at a total of polylog($\epsilon^{-1}$) time steps, we can gather the data required to recover all the $\mathcal{O}(n)$ parameters of the evolution from $\mathcal{O}(\epsilon^{-2}\text{polylog}(n, \epsilon^{-1}))$ samples through the shadow process tomography protocol whenever $k = \mathcal{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.

[1] F. Arute, K. Arya, R. Babbush, D. Bacon, J. C. Bardin, R. Barends, R. Biswas, S. Boixo, F. G. Brando, D. A. Buell, B. Burkett, Y. Chen, Z. Chen, B. Chiaro, R. Collins, W. Courtney, A. Dunsworth,

E. Farhi, B. Foxen, A. Fowler, C. Gidney, M. Giustina, R. Graff, K. Guerin, S. Habegger, M. P. Harrigan, M. J. Hartmann, A. Ho, M. Hoffmann, T. Huang, Travis S. Humble, S. V. Isakov, E. Jeffrey, Z. Jiang,
D. Kafri, K. Kechedzhi, J. Kelly, P. V. Klimov, S. Knysch, A. Korotkov, F. Kostritsa, D. Landhuis, M. Lindmark, E. Lucero, D. Lyakh, S. Mandrà, J. R. McClean, M. McEwen, A. Megrant, X. Mi, K. Michielsen, M. Mohseni, J. Mutus, O. Naaman, M. Neeley, C. Neill, M. Y. Niu, E. Ostby, A. Petukhov, J. C. Platt, C. Quintana, E. G. Rieffel, P. Roushan, N. C. Rubin, D. Sank, K. J. Satzinger, V. Smelyanskiy, K. J. Sung, M. D. Trevithick, A. Vainsencher, B. Villalonga, T. White, Z. J. Yao, P. Yeh, A. Zalcman, H. Neven, and J. M. Martinis. Quantum supremacy using a programmable superconducting processor. *Nature*, 574(7779):505–510, 2019.

[2] H. S. Zhong, H. Wang, Y. H. Deng, M. C. Chen, L. C. Peng, Y. H. Luo, J. Qin, D. Wu, X. Ding, Y. Hu, P. Hu, X. Y. Yang, W. J. Zhang, H. Li, Y. Li, X. Jiang, L. Gan, G. Yang, L. You, Z. Wang, L. Li, N. L. Liu, C. Y. Lu, and J. W. Pan. Quantum computational advantage using photons. *Science*, 370(6523):1460–1463, dec 2020.

[3] P. Scholl, M. Schuler, H. J. Williams, A. A. Eberharter, D. Barredo, K. N. Schymik, V. Lienhard, L. P. Henry, T. C. Lang, T. Lahaye, A. M. Lauchli, and A. Browaeys. Quantum simulation of 2d antiferromagnets with hundreds of rydberg atoms. *Nature*, 595(7866):233–238, jul 2021.

[4] S. Ebadi, T. T. Wang, H. Levine, A. Keesling, G. Semeghini, A. Omran, D. Bluvstein, R. Samajdar, H. Pichler, W. W. Ho, S. Choi, S. Sachdev, M. Greiner, V. Vuletić, and M. D. Lukin. Quantum phases of matter on a 256-atom programmable quantum simulator. *Nature*, 595(7866):227–232, jul 2021.

[5] J. I. Cirac and P. Zoller. Goals and opportunities in quantum simulation. *Nat. Phys.*, 8(4):264–266, apr 2012.

[6] I. L. Chuang M. A. Nielsen. *Quantum Computation and Quantum Information*. Cambridge University Pr., December 2010.

[7] J. Eisert, D. Hangleiter, N. Walk, I. Roth, D. Markham, R. Parekh, U. Chabaud, and E. Kashefi. Quantum certification and benchmarking. *Nat. Rev. Phys.*, 2(7):382–390, jun 2020.

[8] M. D. Shulman, S. P. Harvey, J. M. Nichol, S. D. Bartlett, A. C. Doherty, V. Umansky, and A. Yacoby. Suppressing qubit dephasing using real-time hamiltonian estimation. *Nat. Commun.*, 5(1), oct 2014.

[9] J. Zhang and M. Sarovar. Quantum hamiltonian identification from measurement time traces. *Phys. Rev. Lett.*, 113(8):080401, aug 2014.

[10] J. Zhang and M. Sarovar. Identification of open quantum systems from observable time traces. *Phys. Rev. A*, 91(5):052121, may 2015.

[11] C. Di Franco, M. Paternostro, and M. S. Kim. Hamiltonian tomography in an access-limited setting without state initialization. *Phys. Rev. Lett.*, 102:187203, May 2009.

[12] S. Sheldon, E. Magesan, J. M. Chow, and J. M. Gambetta. Procedure for systematically tuning up cross-talk in the cross-resonance gate. *Phys. Rev. A*, 93(6):060302, jun 2016.

[13] A. Sone and P. Cappellaro. Hamiltonian identifiability assisted by a single-probe measurement. *Phys. Rev. A*, 95(2):022335, feb 2017.

[14] G. O. Samach, A. Greene, J. Borregaard, M. Christandl, D. K. Kim, C. M. McNally, A. Melville, B. M. Niedzielski, Y. Sung, D. Rosenberg, M. E. Schwartz, J. L. Yoder, T. P. Orlando, J. I-J. Wang, S. Gustavsson, M. Kjaergaard, and W. D. Oliver. Lindblad tomography of a superconducting quantum processor. May 2021.

[15] S. T. Wang, D. L. Deng, and L. M. Duan. Hamiltonian tomography for quantum many-body systems with arbitrary couplings. *New J. Phys.*, 17(9):093017, sep 2015.

[16] J. Wang, S. Paesani, R. Santagati, S. Knauer, A. A. Gentile, N. Wiebe, M. Petruzzella, J. L. O'Brien, J. G. Rarity, A. Laing, and M. G. Thompson. Experimental quantum hamiltonian learning. *Nat. Phys.*, 13(6):551–555, mar 2017.

[17] A. Valenti, E. van Nieuwenburg, S. Huber, and E. Grøplova. Hamiltonian learning for quantum error correction. *Phys. Rev. Research*, 1(3):033092, nov 2019.

[18] A. Valenti, G. Jin, J. Léonard, S. D. Huber, and E. Grøplova. Scalable hamiltonian learning for large-scale out-of-equilibrium quantum dynamics. *Physical Review A*, 105(2):023302, 2022.

[19] C. E. Granade, C. Ferrie, N. Wiebe, and D. G. Cory. Robust online hamiltonian learning. *New J. Phys*, 14(10):103013, oct 2012.

[20] N. Wiebe, C. Granade, C. Ferrie, and D. G. Cory. Hamiltonian learning and certification using quantum resources. *Phys. Rev. Lett.*, 112:190501, May 2014.

[21] N. Wiebe, C. Granade, C. Ferrie, and D. Cory. Quantum hamiltonian learning using imperfect quantum resources. *Phys. Rev. A*, 89:042314, Apr 2014.

[22] A. Shabani, M. Mohseni, S. Lloyd, R. L. Kosut, and H. Rabitz. Estimation of many-body quantum hamiltonians via compressive sensing. *Phys. Rev. A*, 84:012107, Jul 2011.

[23] X. L. Qi and D. Ranard. Determining a local hamiltonian from a single eigenstate. *Quantum*, 3:159, jul 2019.

[24] E. Chertkov and B. K. Clark. Computational inverse method for constructing spaces of quantum models from wave functions. *Phys. Rev. X*, 8(3):031029, jul 2018.

[25] E. Bairey, I. Arad, and N. H. Lindner. Learning a local hamiltonian from local measurements. *Phys. Rev. Lett.*, 122:020504, Jan 2019.

[26] Anurag A., Srinivasan A., Tomotaka K., and Mehdi S. Sample-efficient learning of interacting quantum systems. *Nat. Phys.*, 17(8):931–935, may 2021.

[27] T.J. Evans, R. Harper, and S.T. Flammia. Scalable bayesian hamiltonian learning. December 2019.

[28] L. Zhi, Z. Liujun, and H.H. Timothy. Hamiltonian tomography via quantum quench. *Phys. Rev. Lett.*, 124(16):160502, apr 2020.

[29] J. Haah, R. Kothari, and E. Tang. Optimal learning of quantum Hamiltonians from high-temperature Gibbs states, 2021. arXiv:2108. 04842v1.

[30] M. P. da Silva, O. Landon-Cardinal, and D. Poulin. Practical characterization of quantum devices without tomography. *Phys. Rev. Lett.*, 107:210404, Nov 2011.

[31] E. Bairey, C. Guo, D. Poletti, N. H. Lindner, and I. Arad. Learning the dynamics of open quantum systems from their steady states. *N. J. Phys.*, 22(3):032001, mar 2020.

[32] A. Zubida, E. Yitzhaki, N. H. Lindner, and E. Bairey. Optimal short-time measurements for hamiltonian learning. August 2021. arXiv:2108. 08824.

[33] A. Gu, L. Cincio, and P.J. Coles. Practical black box hamiltonian learning. June 2022.

[34] W. Yu, J. Sun, Z. Han, and X. Yuan. Practical and efficient hamiltonian learning. January 2022.

[35] D. Rattacaso, G. Passarelli, and P. Lucignano. High-accuracy hamiltonian learning via delocalized quantum
state evolutions. *Quantum*, 7:905, 2023.

[36] K. Rudinger and R. Joynt. Compressed sensing for hamiltonian reconstruction. *Phys. Rev. A*, 92:052322, Nov 2015.

[37] C. Rouze and D. Stilck França. Learning quantum many-body systems from a few copies, 2021. arXiv:2107.03333v2.

[38] D. Kane, S. Karmalkar, and E. Price. Robust Polynomial Regression up to the Information Theoretic Limit. In *2017 IEEE 58th Annual Symposium on Foundations of Computer Science (FOCS)*, pages 391–402, Berkeley, CA, October 2017. IEEE.

[39] J. Kunjummen, M.C. Tran, D. Carney, and J.M. Taylor. Shadow process tomography of quantum channels, 2021. arXiv:2110.03629v2.

[40] R. Levy, D. Luo, and B.K. Clark. Classical Shadows for Quantum Process Tomography on Near-term Quantum Computers, 2021. arXiv:2110.02965v1.

[41] M. M. Wolf. Quantum channels & operations: Guided tour. *Lecture notes available at http://www-m5.ma.tum.de/foswiki/pub M*, 5, 2012.

[42] H.-Y. Huang, R. Kueng, and J. Preskill. Predicting many properties of a quantum system from very few measurements. *Nature Physics*, 16(10):1050–1057, jun 2020.

[43] E. H. Lieb and D. W. Robinson. The finite group velocity of quantum spin systems. *Commun. Math. Phys*, 28(3):251–257, sep 1972.

[44] D. Poulin. Lieb-Robinson Bound and Locality for General Markovian Quantum Dynamics. *Phys. Rev. Let.*, 104(19):190401, May 2010.

[45] M. Kliesch, C. Gogolin, and J. Eisert. Lieb-Robinson Bounds and the Simulation of Time-Evolution of Local Observables in Lattice Systems. In Volker Bach and Luigi Delle Site, editors, *Many-Electron Approaches in Physics, Chemistry and Mathematics*, pages 301–318. Springer International Publishing, Cham, 2014. Series Title: Mathematical Physics Studies.

[46] M. B. Hastings. Locality in Quantum Systems. arXiv:1008. 5137 [math-ph, physics:quant-ph], August 2010. arXiv: 1008. 5137.

[47] M. Kliesch, C. Gogolin, and J. Eisert. Lieb-robinson bounds and the simulation of time-evolution of local observables in lattice systems. In *Many-Electron Approaches in Physics, Chemistry and Mathematics*, pages 301–318. Springer International Publishing, 2014.

[48] T. Kuwahara and K. Saito. Strictly Linear Light Cones in Long-Range Interacting Systems of Arbitrary Dimensions. *Phys. Rev. X*, 10(3):031010, July 2020.

[49] F. Yan, P. Krantz, Y. Sung, M. Kjaergaard, D. L. Campbell, T. P. Orlando, S. Gustavsson, and W. D. Oliver. Tunable coupling scheme for implementing high-fidelity two-qubit gates. *Phys. Rev. Applied*, 10:054062, Nov 2018.

[50] Y. Sung, L. Ding, J. Braumüller, A. Vepsäläinen, B. Kanan, M. Kjaergaard, A. Greene, G. O. Samach, C. McNally, D. Kim, A. Melville, B. M. Niedzielski, M. E. Schwartz, J. L. Yoder, T. P. Orlando, S. Gustavsson, and W. D. Oliver. Realization of high-fidelity cz and zz-free iswap gates with a tunable coupler. *Phys. Rev. X*, 11:021058, Jun 2021.

[51] L. Vandenberghe S. Boyd. *Convex Optimization*. Cambridge University Press, March 2004.

[52] W. Markoff and J. Grossmann. über Polynome, die in einem gegebenen Intervalle möglichst wenig von Null ab-

weichen. *Mathematische Annalen*, 77(2):213–258, June 1916.

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).

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*ᵢ} and a suitable set of initial states and observables, {{(*ρ*ᵢ, *O*ᵢ)} 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, $\mathcal{L}_{est}$, is obtained from the Master equation.

**Fig. 2:** The median quality of recovery of one 2-qubit coupling using interpolation methods and those based on numerical derivatives [30–32] 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 [30–32] 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 $\sigma_X\sigma_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 quasi-static 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.
# Supplemental Material to Efficient and robust estimation of many-qubit Hamiltonians

Daniel Stilck França,¹,²,* Liubov A. Markovich,³,⁴ V. V.

Dobrovitski,³ Albert H. Werner,¹,⁵ and Johannes Borregaard³

¹QMATH, Department of Mathematical Sciences, University of Copenhagen,
Universitetsparken 5, 2100 Copenhagen, Denmark

²Univ Lyon, ENS Lyon, UCBL, CNRS, Inria, LIP, F-69342, Lyon Cedex 07, France

³QuTech and Kavli Institute of Nanoscience, Delft University of Technology, 2628 CJ, Delft, The Netherlands

⁴Instituut-Lorentz, Universiteit Leiden, P.O. Box 9506, 2300 RA Leiden, The Netherlands

⁵NBIA, Niels Bohr Institute, University of Copenhagen, Blegdamsvej 17, 2100 Copenhagen, Denmark

## Supplementary Note 1: Selecting States and Observables to Isolate Parameters

Any Hamiltonian can be written as

$$ 
\begin{aligned}
H &= \sum_{m_1} \sum_{\alpha_1} a_{\alpha_1}^{(m_1)} \sigma_{\alpha_1}^{(m_1)} + \sum_{m_1, m_2} \sum_{\alpha_1 \alpha_2} a_{\alpha_1 \alpha_2}^{(m_1, m_2)} \sigma_{\alpha_1}^{(m_1)} \sigma_{\alpha_2}^{(m_2)} + \dots \\
&\equiv \sum_{m_1} \sum_{\alpha_1} a_{\alpha_1}^{(m_1)} \mathcal{H}_{\alpha_1}^{(m_1)} + \sum_{m_1, m_2} \sum_{\alpha_1 \alpha_2} a_{\alpha_1 \alpha_2}^{(m_1, m_2)} \mathcal{H}_{\alpha_1 \alpha_2}^{(m_1, m_2)} + \dots,
\end{aligned} 
\quad (1) $$

where Roman indices identify the subspace on which the operator acts, and Greek indices identify the Pauli operator, e.g. $\alpha = 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 $\rho_0$ is described by a Master equation of the form

$$ \left. \frac{d\rho_t}{dt} \right|_{t=0} = -i[H, \rho_0] + \sum_{m=1}^{n} \sum_{\mu, \nu=1}^{3} L_{\mu, \nu}^{(m)} (\sigma_{\mu}^{(m)} \rho_0 \sigma_{\nu}^{(m)\dagger} - \frac{1}{2} \{\sigma_{\nu}^{(m)\dagger} \sigma_{\mu}^{(m)}, \rho_0\}), \quad (2) $$

where $H$ is the Hamiltonian describing the evolution of the system, $L_{\mu, \nu}^{(m)}$ are the elements of the Lindblad matrix, expressed in an operator basis consisting of the different combinations of single-qubit Pauli matrices $\{\sigma\}$. Multiplying it from the right hand side on an observable $O$ and taking the trace, we can write

$$ 
\begin{aligned}
\left. \frac{d}{dt} \operatorname{tr}[\rho_t O] \right|_{t=0} &= -i \sum_{m_1} \sum_{\alpha_1} a_{\alpha_1}^{(m_1)} \operatorname{tr}[[H_{\alpha_1}^{(m_1)}, \rho_0]O] - i \sum_{m_1 m_2} \sum_{\alpha_1 \alpha_2} a_{\alpha_1 \alpha_2}^{(m_1 m_2)} \operatorname{tr}[[H_{\alpha_1 \alpha_2}^{(m_1 m_2)}, \rho_0]O] - \dots \\
&\quad -i \sum_{m_1 \dots m_k} \sum_{\alpha_1 \dots \alpha_k} a_{\alpha_1 \dots \alpha_k}^{(m_1 \dots m_k)} \operatorname{tr}[[H_{\alpha_1 \dots \alpha_k}^{(m_1 \dots m_k)}, \rho_0]O] \dots \\
&\quad + \sum_m \sum_{\mu, \nu} L_{\mu \nu}^{(m)} \operatorname{tr}[(\sigma_{\mu}^{(m)} \rho_0 \sigma_{\nu}^{(m)\dagger} - \frac{1}{2} \{\sigma_{\nu}^{(m)\dagger} \sigma_{\mu}^{(m)}, \rho_0\})O].
\end{aligned} 
\quad (3) $$

Let us introduce the notation

$$ B_{\alpha_1, \dots, \alpha_k}^{(m_1, \dots, m_k)}(\rho_0, O) \equiv -i \operatorname{tr}[[H_{\alpha_1, \dots, \alpha_k}^{(m_1, \dots, m_k)}, \rho_0]O]. \quad (4) $$

To isolate the coefficients $\{a_{\alpha_1}^{(m_1)}, a_{\alpha_1 \alpha_2}^{(m_1 m_2)}, \dots\}$ we observe the $n$ qubit state with the following density matrix

$$ \rho_{\tau_i, \tau_j}^{(i,j)} = \rho_{\tau_i}^{(i)} \otimes \rho_{\tau_j}^{(j)} \otimes \rho_{n-2}, \quad \tau = \{1, 2, 3\} \quad (5) $$

where the $i$ and $j$ are the Pauli qubits, namely $\rho_{\tau_i, \tau_j}^{(i,j)} = (I + \sigma_{\tau_i, \tau_j}^{(i,j)})/2$, and $\rho_{n-2}$ is the density matrix of all other qubits, which we set to be maximally mixed. Using $\sigma_\alpha \sigma_\beta = \delta_{\alpha\beta}I + i\epsilon_{\alpha\beta\gamma}\sigma_\gamma$, we can find the following relations

$$ B_{\alpha_1}^{(i)}(\rho_{\tau_i, \tau_j}^{(i,j)}, O) = \frac{1}{2} \left( \operatorname{tr}[(\sigma_\gamma^{(i)} \otimes I)O] + \operatorname{tr}[(\sigma_\gamma^{(i)} \otimes \sigma_{\tau_j}^{(j)})O] \right), \quad (6) $$

<footnote>* daniel.stilck_franca@ens-lyon.fr</footnote>
$$B_{\alpha_1}^{(j)}(\rho_{\tau_i, \tau_j}^{(i,j)}, O) = \frac{1}{2} \left( \varepsilon_{\alpha_1 \tau_j \zeta} \left( \mathrm{tr} \left[ (I \otimes \sigma_{\zeta}^{(j)}) O \right] + \mathrm{tr} \left[ (\sigma_{\tau_i}^{(i)} \otimes \sigma_{\zeta}^{(j)}) O \right] \right) \right), \quad (7)$$

$$B_{\alpha_1 \alpha_2}^{(ij)}(\rho_{\tau_i, \tau_j}^{(i,j)}, O) = \frac{1}{2} \left( \varepsilon_{\alpha_1 \tau_i \gamma} \left( \delta_{\alpha_2 \tau_j} \mathrm{tr} \left[ (\sigma_{\gamma}^{(i)} \otimes I) O \right] + \mathrm{tr} \left[ (\sigma_{\gamma}^{(i)} \otimes \sigma_{\alpha_2}^{(j)}) O \right] \right) \right. \\ \left. + \varepsilon_{\alpha_2 \tau_j \eta} \left( \delta_{\alpha_1 \tau_i} \mathrm{tr} \left[ (I \otimes \sigma_{\eta}^{(j)}) O \right] + \mathrm{tr} \left[ (\sigma_{\alpha_1}^{(i)} \otimes \sigma_{\eta}^{(j)}) O \right] \right) \right), \quad (8)$$

where O is acting on (i, j) qubits. Selecting $O = \sigma_{\xi_i}^{(i)} \otimes \sigma_{\xi_j}^{(j)}$, we can rewrite the latter matrix element as follows

$$B_{\alpha_1}^{(i)}(\rho_{\tau_i, \tau_j}^{(i,j)}, \sigma_{\xi_i}^{(i)} \otimes \sigma_{\xi_j}^{(j)}) = 2\varepsilon_{\alpha_1 \tau_i \gamma} \delta_{\xi_i \gamma} \delta_{\xi_j \tau_j} = \begin{cases} 2\varepsilon_{\alpha_1 \tau_i \gamma} \delta_{\eta \tau_j} & \text{if } \xi_i = \gamma, \xi_j = \eta, \\ 2\varepsilon_{\alpha_1 \tau_i \gamma} & \text{if } \xi_i = \gamma, \xi_j = \eta, \eta = \tau_j, \\ 0 & \text{else} \end{cases}, \quad (9)$$

$$B_{\alpha_1}^{(j)}(\rho_{\tau_i, \tau_j}^{(i,j)}, \sigma_{\xi_i}^{(i)} \otimes \sigma_{\xi_j}^{(j)}) = 2\varepsilon_{\alpha_1 \tau_j \zeta} \delta_{\xi_j \zeta} \delta_{\xi_i \tau_i} = \begin{cases} 2\varepsilon_{\alpha_1 \tau_j \zeta} \delta_{\kappa \tau_i} & \text{if } \xi_j = \zeta, \xi_i = \kappa, \\ 2\varepsilon_{\alpha_1 \tau_j \zeta} & \text{if } \xi_j = \zeta, \xi_i = \kappa, \kappa = \tau_i, \\ 0 & \text{else} \end{cases}$$

$$B_{\alpha_1 \alpha_2}^{(ij)}(\rho_{\tau_i, \tau_j}^{(i,j)}, \sigma_{\xi_i}^{(i)} \otimes \sigma_{\xi_j}^{(j)}) = 2(\varepsilon_{\alpha_1 \tau_i \gamma} \delta_{\xi_i \gamma} \delta_{\xi_j \alpha_2} + \varepsilon_{\alpha_2 \tau_j \eta} \delta_{\xi_i \alpha_1} \delta_{\xi_j \eta}) \\ = \begin{cases} 2(\varepsilon_{\alpha_1 \tau_i \gamma} \delta_{\eta \alpha_2} + \varepsilon_{\alpha_2 \tau_j \eta} \delta_{\gamma \alpha_1}) & \text{if } \xi_i = \gamma, \xi_j = \eta, \\ 2\varepsilon_{\alpha_1 \tau_i \gamma} \delta_{\tau_j \alpha_2} & \text{if } \xi_i = \gamma, \xi_j = \eta, \eta = \tau_j, \\ 0 & \text{else} \end{cases}$$

where $\eta, \gamma \in \{x, y, z\}$. Selecting other observable $O = \sigma_{\xi_i}^{(i)}$, we can rewrite (6)-(8) differently as

$$B_{\alpha_1}^{(i)}(\rho_{\tau_i, \tau_j}^{(i,j)}, \sigma_{\xi_i}^{(i)}) = \varepsilon_{\alpha_1 \tau_i \gamma} \delta_{\xi_i \gamma} = \begin{cases} \varepsilon_{\alpha_1 \tau_i \gamma} & \text{if } \xi_i = \gamma, \\ 0 & \text{else} \end{cases}, \quad B_{\alpha_1}^{(j)}(\rho_{\tau_i, \tau_j}^{(i,j)}, \sigma_{\xi_i}^{(i)}) = 0, \quad (10)$$

$$B_{\alpha_1 \alpha_2}^{(ij)}(\rho_{\tau_i, \tau_j}^{(i,j)}, \sigma_{\xi_i}^{(i)}) = 2\varepsilon_{\alpha_1 \tau_i \gamma} \delta_{\alpha_2 \tau_j} \delta_{\xi_i \gamma} = \begin{cases} 2\varepsilon_{\alpha_1 \tau_i \gamma} \delta_{\alpha_2 \tau_j} & \text{if } \xi_i = \gamma, \\ 0 & \text{else} \end{cases}.$$

Next, selecting $O = \sigma_{\xi_j}^{(j)}$, we can rewrite (6)-(8) as

$$B_{\alpha_1}^{(i)}(\rho_{\tau_i, \tau_j}^{(i,j)}, O) = 0, \quad B_{\alpha_1}^{(j)}(\rho_{\tau_i, \tau_j}^{(i,j)}, \sigma_{\xi_i}^{(i)}) = \varepsilon_{\alpha_1 \tau_j \zeta} \delta_{\xi_j \zeta} = \begin{cases} \varepsilon_{\alpha_1 \tau_j \zeta} & \text{if } \xi_j = \zeta, \\ 0 & \text{else} \end{cases}, \quad (11)$$

$$B_{\alpha_1 \alpha_2}^{(ij)}(\rho_{\tau_i, \tau_j}^{(i,j)}, \sigma_{\xi_j}^{(j)}) = 2\varepsilon_{\alpha_2 \tau_j \eta} \delta_{\alpha_1 \tau_i} \delta_{\xi_j \eta} = \begin{cases} 2\varepsilon_{\alpha_2 \tau_j \eta} \delta_{\alpha_1 \tau_i} & \text{if } \xi_j = \eta, \\ 0 & \text{else} \end{cases}.$$

For $k \ge 2$ we can write the general matrix element:

$$B_{\alpha_1, \dots, \alpha_k}^{(m_1, \dots, m_k)}(\rho_{\tau_i, \tau_j}^{(i,j)}, O) = \frac{1}{2} \left( \varepsilon_{\alpha_i \tau_i \gamma} \left( \mathrm{tr} \left[ (\sigma_{\alpha_1}^{(1)} \otimes \dots \otimes \sigma_{\gamma}^{(i)} \otimes \dots \otimes \sigma_{\alpha_j}^{(j)} \otimes \dots \otimes \sigma_{\alpha_k}^{(k)} ) O \right] \right. \right. \\ \left. + \delta_{\alpha_j \tau_j} \mathrm{tr} \left[ (\sigma_{\alpha_1}^{(1)} \otimes \dots \otimes \sigma_{\gamma}^{(i)} \otimes \dots \otimes I^{(j)} \otimes \dots \otimes \sigma_{\alpha_k}^{(k)} ) O \right] \right) \\ + \varepsilon_{\alpha_j \tau_j \eta} \left( \mathrm{tr} \left[ (\sigma_{\alpha_1}^{(1)} \otimes \dots \otimes \sigma_{\alpha_i}^{(i)} \otimes \dots \otimes \sigma_{\eta}^{(j)} \otimes \dots \otimes \sigma_{\alpha_k}^{(k)} ) O \right] \right. \\ \left. + \delta_{\alpha_i \tau_i} \mathrm{tr} \left[ (\sigma_{\alpha_1}^{(1)} \otimes \dots \otimes I^{(i)} \otimes \dots \otimes \sigma_{\eta}^{(j)} \otimes \dots \otimes \sigma_{\alpha_k}^{(k)} ) O \right] \right) \right), \quad (12)$$

where O is acting on $(1, \dots, k)$ qubits. Let $O = \sigma_{\xi_1}^{(1)} \otimes \dots \otimes \sigma_{\xi_k}^{(k)}$ holds. Then we can rewrite (12) as follows

$$\begin{align} B_{\alpha_1, \dots, \alpha_k}^{(1, \dots, k)}(\rho_{\tau_i, \tau_j}^{(i,j)}, \sigma_{\xi_1}^{(1)} \otimes \dots \otimes \sigma_{\xi_k}^{(k)}) &\equiv -i \mathrm{Tr} \left( [\sigma_{\alpha_1}^{(1)} \otimes \dots \otimes \sigma_{\alpha_k}^{(k)}, \rho_{\tau_i, \tau_j}^{(i,j)}] (\sigma_{\xi_1}^{(1)} \otimes \dots \otimes \sigma_{\xi_k}^{(k)}) \right) \nonumber \\ &= 2^{k-1} (\varepsilon_{\alpha_i \tau_i \gamma} \delta_{\alpha_1 \xi_1} \dots \delta_{\gamma \xi_i} \dots \delta_{\alpha_j \xi_j} \dots \delta_{\alpha_k \xi_k} + \varepsilon_{\alpha_j \tau_j \eta} \delta_{\alpha_1 \xi_1} \dots \delta_{\alpha_i \xi_i} \dots \delta_{\eta \xi_j} \dots \delta_{\alpha_k \xi_k}) \tag{13} \\ &= \begin{cases} 2^{k-1} (\varepsilon_{\alpha_i \tau_i \gamma} \delta_{\eta \alpha_j} + \varepsilon_{\alpha_j \tau_j \eta} \delta_{\gamma \alpha_i}) & \text{if } \xi_i = \gamma, \xi_j = \eta, \xi_1, \dots, k = \alpha_1, \dots, k, \\ 0 & \text{else} \end{cases}. \nonumber \end{align}$$
From this result, for $k=3$ we get

$$B_{\alpha_i \alpha_j \alpha_l}^{(i,j,l)} (\rho_{\tau_i, \tau_j}^{(i,j)}, \sigma_{\xi_i}^{(i)} \otimes \sigma_{\xi_j}^{(j)} \otimes \sigma_{\xi_l}^{(l)}) = 2(\varepsilon_{\alpha_i \tau_i \gamma} \delta_{\gamma, \xi_i} \delta_{\alpha_j \xi_j} + \varepsilon_{\alpha_j \tau_j \eta} \delta_{\alpha_i \xi_i} \delta_{\eta \xi_j}) \delta_{\alpha_l \xi_l} \quad (14)$$

To isolate $\{L_{\mu\nu}^{(m)}\}$ we recall the assumption that $\rho_{n-2}$ is the density matrix of the maximally mixed state. Then the Lindblad part of the equation (3) for an observable $O = \sigma_{\xi_i}^{(i)} \otimes \sigma_{\xi_j}^{(j)}$ is

$$ \begin{aligned} \mathcal{L}_{\mu\nu}^{k}(\rho_{\tau_i, \tau_j}^{(i,j)}, \sigma_{\xi_i}^{(i)} \otimes \sigma_{\xi_j}^{(j)}) &\equiv L_{\mu\nu}^{(k)} \operatorname{tr} \left[ \left( \left( \sigma_{\mu}^{(k)} \rho_{\tau_i, \tau_j}^{(i,j)} \sigma_{\nu}^{(k)} - \frac{1}{2} \{ \sigma_{\nu}^{(k)} \sigma_{\mu}^{(k)}, \rho_{\tau_i, \tau_j}^{(i,j)} \} \right) (\sigma_{\xi_i}^{(i)} \otimes \sigma_{\xi_j}^{(j)}) \right] \right. \\ &= \left. \begin{cases} L_{\mu\nu}^{(i)}(2i\varepsilon_{\mu\nu\gamma}\delta_{\gamma\xi_i} + 2\delta_{\nu\tau_i}\delta_{\mu\xi_i} - \frac{3}{2}\delta_{\mu\nu}\delta_{\xi_i\tau_i})\delta_{\tau_j\xi_j} & \text{if } k=i \\ L_{\mu\nu}^{(j)}(2i\varepsilon_{\mu\nu\gamma}\delta_{\gamma\xi_j} + 2\delta_{\nu\tau_j}\delta_{\mu\xi_j} - \frac{3}{2}\delta_{\mu\nu}\delta_{\xi_j\tau_j})\delta_{\tau_i\xi_i} & \text{if } k=j \\ 0 & \text{else} \end{cases} \right] \end{aligned} \quad (15) $$

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

$$ \mathcal{L}_{\mu\nu}^{i}(\rho_{\tau_i, \tau_j}^{(i,j)}, \sigma_{\xi_i}^{(i)} \otimes \sigma_{\xi_j}^{(j)}) = L_{\mu\nu}^{(i)}(2i\varepsilon_{\mu\nu\gamma} + 2\delta_{\nu\tau_i}\delta_{\mu\gamma} - \frac{3}{2}\delta_{\mu\nu}\delta_{\gamma\tau_i})\delta_{\tau_j\eta} \quad \text{if } \xi_i = \gamma, \xi_j = \eta, \quad (16) $$

and

$$ \mathcal{L}_{\mu\nu}^{j}(\rho_{\tau_i, \tau_j}^{(i,j)}, \sigma_{\xi_i}^{(i)} \otimes \sigma_{\xi_j}^{(j)}) = L_{\mu\nu}^{(j)}(2i\varepsilon_{\mu\nu\gamma}\delta_{\gamma\eta} + 2\delta_{\nu\tau_j}\delta_{\mu\eta} - \frac{3}{2}\delta_{\mu\nu}\delta_{\eta\tau_j})\delta_{\tau_i\gamma} \quad \text{if } \xi_i = \gamma, \xi_j = \eta. \quad (17) $$

However, for an observable $O = \sigma_{\xi_i}^{(i)}$ the Lindblad part of the equation (3) is the following

$$ \begin{aligned} \mathcal{L}_{\mu\nu}^{k}(\rho_{\tau_i, \tau_j}^{(i,j)}, \sigma_{\xi_i}^{(i)} \otimes \sigma_{\xi_j}^{(j)}) &\equiv L_{\mu\nu}^{(k)} \operatorname{tr} \left[ \left( \left( \sigma_{\mu}^{(k)} \rho_{\tau_i, \tau_j}^{(i,j)} \sigma_{\nu}^{(k)} - \frac{1}{2} \{ \sigma_{\nu}^{(k)} \sigma_{\mu}^{(k)}, \rho_{\tau_i, \tau_j}^{(i,j)} \} \right) (\sigma_{\xi_i}^{(i)}) \right] \right. \\ &= \left. \begin{cases} L_{\mu\nu}^{(i)}(2i\varepsilon_{\mu\nu\gamma}\delta_{\gamma\xi_i} + \delta_{\nu\tau_i}\delta_{\mu\xi_i} - \frac{3}{2}\delta_{\mu\nu}\delta_{\xi_i\tau_i} + \delta_{\mu\xi_i}\delta_{\tau_i\nu}) & \text{if } k=i \\ 0 & \text{if } k=j \\ 0 & \text{else} \end{cases} \right] \end{aligned} \quad (18) $$

Substituting the conditions (10) in (18), we get

$$ \mathcal{L}_{\mu\nu}^{i}(\rho_{\tau_i, \tau_j}^{(i,j)}, \sigma_{\xi_i}^{(i)}) = L_{\mu\nu}^{(i)}(2i\varepsilon_{\mu\nu\gamma} + 2\delta_{\nu\tau_i}\delta_{\mu\gamma} - \frac{3}{2}\delta_{\mu\nu}\delta_{\gamma\tau_i}), \quad \text{if } \xi_i = \gamma. \quad (19) $$

Next, for an observable $O = \sigma_{\xi_i}^{(i)} \otimes \sigma_{\xi_j}^{(j)} \otimes \sigma_{\xi_l}^{(l)}$ we can write

$$ \begin{aligned} \mathcal{L}_{\mu\nu}^{k}(\rho_{\tau_i, \tau_j}^{(i,j)}, \sigma_{\xi_i}^{(i)} \otimes \sigma_{\xi_j}^{(j)} \otimes \sigma_{\xi_l}^{(l)}) &\equiv L_{\mu\nu}^{(k)} \operatorname{tr} \left[ \left( \left( \sigma_{\mu}^{(k)} \rho_{\tau_i, \tau_j}^{(i,j)} \sigma_{\nu}^{(k)} - \frac{1}{2} \{ \sigma_{\nu}^{(k)} \sigma_{\mu}^{(k)}, \rho_{\tau_i, \tau_j}^{(i,j)} \} \right) (\sigma_{\xi_i}^{(i)} \otimes \sigma_{\xi_j}^{(j)} \otimes \sigma_{\xi_l}^{(l)}) \right] \\ &= \left. \begin{cases} 2iL_{\mu\nu}^{(l)}\varepsilon_{\mu\nu\gamma}\delta_{\gamma\xi_i}\delta_{\tau_i\xi_i}\delta_{\tau_j\xi_j} & \text{if } k=l \\ 0 & \text{else.} \end{cases} \right] \end{aligned} \quad (20) $$

Then, according to the results of the previous subsection, we get

$$ \mathcal{L}_{\mu\nu}^{l}(\rho_{\tau_i, \tau_j}^{(i,j)}, \sigma_{\xi_i}^{(i)} \otimes \sigma_{\xi_j}^{(j)} \otimes \sigma_{\xi_l}^{(l)}) = 2iL_{\mu\nu}^{(l)}\varepsilon_{\mu\nu\gamma}\delta_{\gamma\alpha_l}\delta_{\tau_i\gamma}\delta_{\tau_j\eta} \quad \text{if } \xi_i = \gamma, \xi_j = \eta, \xi_l = \alpha_l. \quad (21) $$

Finally, for an observable $O = \sigma_{\xi_1}^{(1)} \otimes \dots \sigma_{\xi_i}^{(i)} \otimes \dots \otimes \sigma_{\xi_j}^{(j)} \otimes \dots \otimes \sigma_{\xi_k}^{(k)}$, $k>3$, the Lindblad part of the equation (3) is the following

$$ L_{\mu\nu}^{(m)} \operatorname{tr} \left[ \left( \left( \sigma_{\mu}^{(m)} \rho_{\tau_i, \tau_j}^{(i,j)} \sigma_{\nu}^{(m)} - \frac{1}{2} \{ \sigma_{\nu}^{(m)} \sigma_{\mu}^{(m)}, \rho_{\tau_i, \tau_j}^{(i,j)} \} \right) (\sigma_{\xi_1}^{(1)} \otimes \dots \otimes \sigma_{\xi_i}^{(i)} \otimes \dots \otimes \sigma_{\xi_j}^{(j)} \otimes \dots \otimes \sigma_{\xi_k}^{(k)}) \right] = 0. \quad (22) $$
A. Final Results

After we selected the different observable operators and defined the density matrix $\rho_0 = \rho_{\tau_i}^{(i)} \otimes \rho_{\tau_j}^{(j)} \otimes \frac{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 = \sigma_{\gamma}^{(i)} \otimes \sigma_{\eta}^{(j)}$ we can write the equation (3) as follows

$$ \frac{d}{dt} \mathrm{tr} \left[ \rho_t (\sigma_{\gamma}^{(i)} \otimes \sigma_{\eta}^{(j)}) \right] \bigg|_{t=0} = -2 \sum_{\alpha_1} a_{\alpha_1}^{(i)} \varepsilon_{\alpha_1 \tau_i \gamma} \delta_{\eta \tau_j} - 2 \sum_{\alpha_1} a_{\alpha_1}^{(j)} \varepsilon_{\alpha_1 \tau_j \eta} \delta_{\gamma \tau_i} \tag{23} \\ -2 \sum_{\alpha_1, \alpha_2} a_{\alpha_1 \alpha_2}^{(ij)} (\varepsilon_{\alpha_1 \tau_i \gamma} \delta_{\eta \alpha_2} + \varepsilon_{\alpha_2 \tau_j \eta} \delta_{\gamma \alpha_1}) + \sum_{\mu, \nu} \left( L_{\mu\nu}^{(i)} (2i\varepsilon_{\mu\nu\gamma} + 2\delta_{\nu\tau_i}\delta_{\mu\gamma} - \frac{3}{2}\delta_{\mu\nu}\delta_{\gamma\tau_i}) \delta_{\tau_j\eta} \right. \\ \left. + L_{\mu\nu}^{(j)} (2i\varepsilon_{\mu\nu\gamma} \delta_{\gamma\eta} + 2\delta_{\nu\tau_j}\delta_{\mu\eta} - \frac{3}{2}\delta_{\mu\nu}\delta_{\eta\tau_j}) \delta_{\tau_i\gamma} \right). $$

Selecting $\tau_j \neq \eta$, $\tau_i \neq \gamma$, we can isolate the coefficients of the type $a_{\alpha_1 \alpha_2}^{(ij)}$ in (23), namely

$$ \frac{d}{dt} \mathrm{tr} \left[ \rho_t (\sigma_{\gamma}^{(i)} \otimes \sigma_{\eta}^{(j)}) \right] \bigg|_{t=0} = -2 \sum_{\alpha_1, \alpha_2} a_{\alpha_1 \alpha_2}^{(ij)} (\varepsilon_{\alpha_1 \tau_i \gamma} \delta_{\eta \alpha_2} + \varepsilon_{\alpha_2 \tau_j \eta} \delta_{\gamma \alpha_1}). \tag{24} $$

Let us call $\rho_{t\tau_i\tau_j}$ the density matrix evaluated by the Hamiltonian evolution from $\rho_0$. From (24) we can find $a_{\alpha_1\alpha_2}^{(ij)}$. To this end, we select $\gamma = y$, $\eta = y$ and four pairs $\tau_i, \tau_j \in \{z, x; z, z; x, x; x, z\}$ to get the system of equations

$$ \frac{d}{dt} \mathrm{tr} \left[ \rho_{txz} (\sigma_y^{(i)} \otimes \sigma_y^{(j)}) \right] \bigg|_{t=0} = 2(a_{xy}^{(ij)} - a_{yz}^{(ij)}), \quad \frac{d}{dt} \mathrm{tr} \left[ \rho_{tzz} (\sigma_y^{(i)} \otimes \sigma_y^{(j)}) \right] \bigg|_{t=0} = 2(a_{xy}^{(ij)} + a_{yz}^{(ij)}), \tag{25} $$

$$ \frac{d}{dt} \mathrm{tr} \left[ \rho_{txx} (\sigma_y^{(i)} \otimes \sigma_y^{(j)}) \right] \bigg|_{t=0} = -2(a_{yz}^{(ij)} + a_{zy}^{(ij)}), \quad \frac{d}{dt} \mathrm{tr} \left[ \rho_{txz} (\sigma_y^{(i)} \otimes \sigma_y^{(j)}) \right] \bigg|_{t=0} = -2(a_{yx}^{(ij)} + a_{zy}^{(ij)}). $$

Since $a_{yz}^{(ij)} = a_{zy}^{(ij)}$ and $a_{xy}^{(ij)} = a_{yx}^{(ij)}$, we can write

$$ a_{xy}^{(ij)} = \frac{1}{4} \frac{d}{dt} \mathrm{tr} \left[ \rho_{tzz} (\sigma_y^{(i)} \otimes \sigma_y^{(j)}) \right] \bigg|_{t=0}, \quad a_{yz}^{(ij)} = -\frac{1}{4} \frac{d}{dt} \mathrm{tr} \left[ \rho_{txx} (\sigma_y^{(i)} \otimes \sigma_y^{(j)}) \right] \bigg|_{t=0}. \tag{26} $$

Selecting $\gamma = x$, $\eta = y$ and four pairs $\tau_i, \tau_j \in \{y, z; y, x; z, z; z, x\}$ to get the system of equations

$$ \frac{d}{dt} \mathrm{tr} \left[ \rho_{tyz} (\sigma_x^{(i)} \otimes \sigma_y^{(j)}) \right] \bigg|_{t=0} = 2(a_{xx}^{(ij)} + a_{zy}^{(ij)}), \quad \frac{d}{dt} \mathrm{tr} \left[ \rho_{tyx} (\sigma_x^{(i)} \otimes \sigma_y^{(j)}) \right] \bigg|_{t=0} = 2(a_{zy}^{(ij)} - a_{xz}^{(ij)}), \tag{27} $$

$$ \frac{d}{dt} \mathrm{tr} \left[ \rho_{tzz} (\sigma_x^{(i)} \otimes \sigma_y^{(j)}) \right] \bigg|_{t=0} = 2(a_{xx}^{(ij)} - a_{yy}^{(ij)}), \quad \frac{d}{dt} \mathrm{tr} \left[ \rho_{txz} (\sigma_x^{(i)} \otimes \sigma_y^{(j)}) \right] \bigg|_{t=0} = -2(a_{yy}^{(ij)} + a_{xz}^{(ij)}). $$

Hence

$$ a_{xx}^{(ij)} = \frac{1}{2} \frac{d}{dt} \mathrm{tr} \left[ \rho_{tyz} (\sigma_x^{(i)} \otimes \sigma_y^{(j)}) \right] \bigg|_{t=0} + \frac{1}{4} \frac{d}{dt} \mathrm{tr} \left[ \rho_{txx} (\sigma_y^{(i)} \otimes \sigma_y^{(j)}) \right] \bigg|_{t=0}, \tag{28} $$

$$ a_{xz}^{(ij)} = -\frac{1}{2} \frac{d}{dt} \mathrm{tr} \left[ \rho_{tyx} (\sigma_x^{(i)} \otimes \sigma_y^{(j)}) \right] \bigg|_{t=0} - \frac{1}{4} \frac{d}{dt} \mathrm{tr} \left[ \rho_{txx} (\sigma_y^{(i)} \otimes \sigma_y^{(j)}) \right] \bigg|_{t=0}, $$

$$ a_{yy}^{(ij)} = -\frac{1}{2} \frac{d}{dt} \mathrm{tr} \left[ \rho_{tzz} (\sigma_x^{(i)} \otimes \sigma_y^{(j)}) \right] \bigg|_{t=0} + \frac{1}{2} \frac{d}{dt} \mathrm{tr} \left[ \rho_{tyz} (\sigma_x^{(i)} \otimes \sigma_y^{(j)}) \right] \bigg|_{t=0} \\ + \frac{1}{4} \frac{d}{dt} \mathrm{tr} \left[ \rho_{txx} (\sigma_y^{(i)} \otimes \sigma_y^{(j)}) \right] \bigg|_{t=0}. $$

Selecting $\gamma = x$, $\eta = z$ and four pairs $\tau_i, \tau_j \in \{y, y; z, x; y, x; z, y\}$, we get

$$ \frac{d}{dt} \mathrm{tr} \left[ \rho_{tyy} (\sigma_x^{(i)} \otimes \sigma_z^{(j)}) \right] \bigg|_{t=0} = -2(a_{xx}^{(ij)} - a_{zz}^{(ij)}), \quad \frac{d}{dt} \mathrm{tr} \left[ \rho_{txz} (\sigma_x^{(i)} \otimes \sigma_z^{(j)}) \right] \bigg|_{t=0} = 2(a_{xy}^{(ij)} - a_{yz}^{(ij)}), \tag{29} $$

$$ \frac{d}{dt} \mathrm{tr} \left[ \rho_{tyx} (\sigma_x^{(i)} \otimes \sigma_z^{(j)}) \right] \bigg|_{t=0} = -2(a_{zz}^{(ij)} - a_{xy}^{(ij)}), \quad \frac{d}{dt} \mathrm{tr} \left[ \rho_{tzy} (\sigma_x^{(i)} \otimes \sigma_z^{(j)}) \right] \bigg|_{t=0} = -2(a_{xx}^{(ij)} + a_{yz}^{(ij)}). $$

Hence the last coefficient is

$$ a_{zz}^{(ij)} = -\frac{1}{2} \frac{d}{dt} \mathrm{tr} \left[ \rho_{tyx} (\sigma_x^{(i)} \otimes \sigma_z^{(j)}) \right] \bigg|_{t=0} + \frac{1}{4} \frac{d}{dt} \mathrm{tr} \left[ \rho_{tzz} (\sigma_y^{(i)} \otimes \sigma_y^{(j)}) \right] \bigg|_{t=0}. \tag{30} $$
To find the other coefficients we select $\tau_j = \eta$, $\tau_i \neq \gamma$ and rereswrite (23) as

$$ \frac{d}{dt} \mathrm{tr} [\rho_t(\sigma_{\gamma}^{(i)} \otimes \sigma_{\eta}^{(j)})] |_{t=0} = -2 \sum_{\alpha_1} a_{\alpha_1}^{(i)} \varepsilon_{\alpha_1 \tau_i \gamma} - 2 \sum_{\alpha_1, \alpha_2} a_{\alpha_1 \alpha_2}^{(ij)} \varepsilon_{\alpha_1 \tau_i \gamma} \delta_{\eta \alpha_2} \\ + 2 \sum_{\mu=x,y,z} L_{\mu\mu}^{(i)} \delta_{\mu\tau_i} \delta_{\mu\gamma} + 2 \sum_{\mu,\nu=x,y,z} L_{\mu\nu}^{(i)} (i\varepsilon_{\mu\nu\gamma} + \delta_{\nu\tau_i} \delta_{\mu\gamma}). \quad (31) $$

Next, for $\tau_j \neq \eta$, $\tau_i = \gamma$, we can rewrite (23) as

$$ \frac{d}{dt} \mathrm{tr} [\rho_t(\sigma_{\gamma}^{(i)} \otimes \sigma_{\eta}^{(j)})] |_{t=0} = -2 \sum_{\alpha_1} a_{\alpha_1}^{(j)} \varepsilon_{\alpha_1 \tau_j \eta} - 2 \sum_{\alpha_1, \alpha_2} a_{\alpha_1 \alpha_2}^{(ij)} \varepsilon_{\alpha_2 \tau_j \eta} \delta_{\gamma \alpha_1} \\ + 2 \sum_{\mu=x,y,z} L_{\mu\mu}^{(i)} \delta_{\mu\tau_j} \delta_{\mu\eta} + 2 \sum_{\mu,\nu=x,y,z} L_{\mu\nu}^{(i)} (i\varepsilon_{\mu\nu\gamma} \delta_{\gamma\eta} + \delta_{\nu\tau_j} \delta_{\mu\eta}). \quad (32) $$

Selecting an observable $O = \sigma_{\gamma}^{(i)}$, we can write

$$ \frac{d}{dt} \mathrm{tr} [\rho_t \sigma_{\gamma}^{(i)}] |_{t=0} = - \sum_{\alpha_1} a_{\alpha_1}^{(i)} \varepsilon_{\alpha_1 \tau_i \gamma} - 2 \sum_{\alpha_1, \alpha_2} a_{\alpha_1 \alpha_2}^{(ij)} \varepsilon_{\alpha_1 \tau_i \gamma} \delta_{\alpha_2 \tau_j} \\ + 2 \sum_{\mu=x,y,z} L_{\mu\mu}^{(i)} (\delta_{\mu\tau_i} \delta_{\mu\gamma} - \frac{3}{4} \delta_{\gamma\tau_i}) + 2 \sum_{\mu,\nu=x,y,z} L_{\mu\nu}^{(i)} (i\varepsilon_{\mu\nu\gamma} + \delta_{\nu\tau_i} \delta_{\mu\gamma}). \quad (33) $$

For the other observable $O = \sigma_{\eta}^{(j)}$, the result is the following

$$ \frac{d}{dt} \mathrm{tr} [\rho_t \sigma_{\eta}^{(j)}] |_{t=0} = - \sum_{\alpha_1} a_{\alpha_1}^{(j)} \varepsilon_{\alpha_1 \tau_j \eta} - 2 \sum_{\alpha_1, \alpha_2} a_{\alpha_1 \alpha_2}^{(ij)} \varepsilon_{\alpha_2 \tau_j \eta} \delta_{\alpha_1 \tau_i} \\ + 2 \sum_{\mu=x,y,z} L_{\mu\mu}^{(j)} (\delta_{\mu\tau_j} \delta_{\mu\eta} - \frac{3}{4} \delta_{\eta\tau_j}) + 2 \sum_{\mu,\nu=x,y,z} L_{\mu\nu}^{(j)} (i\varepsilon_{\mu\nu\gamma} \delta_{\gamma\eta} + \delta_{\nu\tau_j} \delta_{\mu\eta}). \quad (34) $$

Substituting (33) in (31), we can write

$$ \frac{d}{dt} \mathrm{tr} [\rho_t(\sigma_{\gamma}^{(i)} \otimes \sigma_{\eta}^{(j)})] |_{t=0} - \frac{d}{dt} \mathrm{tr} [\rho_t \sigma_{\gamma}^{(i)}] |_{t=0} = - \sum_{\alpha_1} a_{\alpha_1}^{(i)} \varepsilon_{\alpha_1 \tau_i \gamma} \\ + 2 \sum_{\alpha_1, \alpha_2} a_{\alpha_1 \alpha_2}^{(ij)} \varepsilon_{\alpha_1 \tau_i \gamma} (\delta_{\tau_j \alpha_2} - \delta_{\eta \alpha_2}) + \frac{3}{2} \sum_{\mu=x,y,z} L_{\mu\mu}^{(i)} \delta_{\gamma \tau_i}. \quad (35) $$

Since in (31) the conditions $\tau_j = \eta$, $\tau_i \neq \gamma$, hold, we can rewrite the latter equation as

$$ \frac{d}{dt} \mathrm{tr} [\rho_t(\sigma_{\gamma}^{(i)} \otimes \sigma_{\eta}^{(j)})] |_{t=0} - \frac{d}{dt} \mathrm{tr} [\rho_t \sigma_{\gamma}^{(i)}] |_{t=0} = - \sum_{\alpha_1} a_{\alpha_1}^{(i)} \varepsilon_{\alpha_1 \tau_i \gamma}. \quad (36) $$

Solving the latter equation, we find $a_{\alpha_1}^{(i)}$. To this end, we select $\gamma = y$, $\tau_i = x$, $\tau_j = \eta \in \{x, y, z\}$ and get

$$ a_x^{(i)} = -\frac{d}{dt} \mathrm{tr} [\rho_{tx\eta}(\sigma_y^{(i)} \otimes \sigma_{\eta}^{(j)})] |_{t=0} + \frac{d}{dt} \mathrm{tr} [\rho_{tx\eta} \sigma_y^{(i)}] |_{t=0}. \quad (37) $$

Next, for $\gamma = x$, $\tau_i = z$, $\tau_j = \eta \in \{x, y, z\}$ we get the solution of (36), namely

$$ a_y^{(i)} = -\frac{d}{dt} \mathrm{tr} [\rho_{tz\eta}(\sigma_x^{(i)} \otimes \sigma_{\eta}^{(j)})] |_{t=0} + \frac{d}{dt} \mathrm{tr} [\rho_{tz\eta} \sigma_x^{(i)}] |_{t=0}. \quad (38) $$

Finally, for $\gamma = x$, $\tau_i = y$, $\tau_j = \eta \in \{x, y, z\}$ the solution is

$$ a_z^{(i)} = \frac{d}{dt} \mathrm{tr} [\rho_{ty\eta}(\sigma_x^{(i)} \otimes \sigma_{\eta}^{(j)})] |_{t=0} - \frac{d}{dt} \mathrm{tr} [\rho_{ty\eta} \sigma_x^{(i)}] |_{t=0}. \quad (39) $$
Substituting (34) in (32), we can write

$$ \frac{d}{dt} \mathrm{tr} \left[ \rho_t (\sigma_\gamma^{(i)} \otimes \sigma_\eta^{(j)}) \right] \Big|_{t=0} - \frac{d}{dt} \mathrm{tr} \left[ \rho_t \sigma_\eta^{(j)} \right] \Big|_{t=0} = - \sum_{\alpha_1} a_{\alpha_1}^{(j)} \varepsilon_{\alpha_1 \tau_j \eta} \\ + 2 \sum_{\alpha_1, \alpha_2} a_{\alpha_1 \alpha_2}^{(ij)} \varepsilon_{\alpha_2 \tau_j \eta} (\delta_{\tau_i \alpha_1} - \delta_{\gamma \alpha_1}) + \frac{3}{2} \sum_{\mu=x,y,z} L_{\mu\mu}^{(i)} \delta_{\eta \tau_j}. \quad (40) $$

Since $\tau_j \neq \eta$, $\tau_i = \gamma$, hold, we can rewrite it as follows

$$ \frac{d}{dt} \mathrm{tr} \left[ \rho_t (\sigma_\gamma^{(i)} \otimes \sigma_\eta^{(j)}) \right] \Big|_{t=0} - \frac{d}{dt} \mathrm{tr} \left[ \rho_t \sigma_\eta^{(j)} \right] \Big|_{t=0} = - \sum_{\alpha_1} a_{\alpha_1}^{(j)} \varepsilon_{\alpha_1 \tau_j \eta}. \quad (41) $$

Solving the latter equation, we find $a_{\alpha_1}^{(j)}$. Selecting $\eta = y$, $\tau_j = z$ and $\tau_i = \gamma \in \{x, y, z\}$, we get

$$ a_x^{(j)} = \frac{d}{dt} \mathrm{tr} \left[ \rho_t \gamma z (\sigma_\gamma^{(i)} \otimes \sigma_y^{(j)}) \right] \Big|_{t=0} - \frac{d}{dt} \mathrm{tr} \left[ \rho_t \gamma z \sigma_y^{(j)} \right] \Big|_{t=0}. \quad (42) $$

Selecting $\eta = x$, $\tau_j = z$ and $\tau_i = \gamma \in \{x, y, z\}$, we get

$$ a_y^{(j)} = -\frac{d}{dt} \mathrm{tr} \left[ \rho_t \gamma z (\sigma_\gamma^{(i)} \otimes \sigma_x^{(j)}) \right] \Big|_{t=0} + \frac{d}{dt} \mathrm{tr} \left[ \rho_t \gamma z \sigma_x^{(j)} \right] \Big|_{t=0}. \quad (43) $$

Finally, selecting $\eta = y$, $\tau_j = x$ and $\tau_i = \gamma \in \{x, y, z\}$, we get the last coefficient of this type

$$ a_z^{(j)} = -\frac{d}{dt} \mathrm{tr} \left[ \rho_t \gamma x (\sigma_\gamma^{(i)} \otimes \sigma_y^{(j)}) \right] \Big|_{t=0} + \frac{d}{dt} \mathrm{tr} \left[ \rho_t \gamma x \sigma_y^{(j)} \right] \Big|_{t=0}. \quad (44) $$

All the coefficients with the corresponding observables and initial states are given in Table I.

<table><thead><tr><th>a<sub>&alpha;</sub><sup>(i)</sup></th><th>(O, {&tau;<sub>i</sub><sup>(i)</sup>, &tau;<sub>j</sub><sup>(j)</sup>}</th><th>Equation</th></tr></thead><tbody><tr><td>a<sub>x</sub><sup>(i)</sup></td><td>(&sigma;<sub>y</sub><sup>(i)</sup> &otimes; &sigma;<sub>&eta;</sub><sup>(j)</sup>, {x, &eta;}); (&sigma;<sub>y</sub><sup>(i)</sup>, {x, &eta;}), &eta; &isin; {x, y, z}</td><td>(37)</td></tr><tr><td>a<sub>y</sub><sup>(i)</sup></td><td>(&sigma;<sub>x</sub><sup>(i)</sup> &otimes; &sigma;<sub>&eta;</sub><sup>(j)</sup>, {z, &eta;}); (&sigma;<sub>x</sub><sup>(i)</sup>, {z, &eta;}), &eta; &isin; {x, y, z}</td><td>(38)</td></tr><tr><td>a<sub>z</sub><sup>(i)</sup></td><td>(&sigma;<sub>x</sub><sup>(i)</sup> &otimes; &sigma;<sub>&eta;</sub><sup>(j)</sup>, {y, &eta;}); (&sigma;<sub>x</sub><sup>(i)</sup>, {y, &eta;}), &eta; &isin; {x, y, z}</td><td>(39)</td></tr><tr><td>a<sub>x</sub><sup>(j)</sup></td><td>(&sigma;<sub>&gamma;</sub><sup>(i)</sup> &otimes; &sigma;<sub>y</sub><sup>(j)</sup>, {&gamma;, z}); (&sigma;<sub>y</sub><sup>(i)</sup>, {&gamma;, z}), &gamma; &isin; {x, y, z}</td><td>(42)</td></tr><tr><td>a<sub>y</sub><sup>(j)</sup></td><td>(&sigma;<sub>&gamma;</sub><sup>(i)</sup> &otimes; &sigma;<sub>x</sub><sup>(j)</sup>, {&gamma;, z}); (&sigma;<sub>x</sub><sup>(i)</sup>, {&gamma;, z}), &gamma; &isin; {x, y, z}</td><td>(43)</td></tr><tr><td>a<sub>z</sub><sup>(j)</sup></td><td>(&sigma;<sub>&gamma;</sub><sup>(i)</sup> &otimes; &sigma;<sub>y</sub><sup>(j)</sup>, {&gamma;, x}); (&sigma;<sub>y</sub><sup>(i)</sup>, {&gamma;, x}), &gamma; &isin; {x, y, z}</td><td>(44)</td></tr><tr><td>a<sub>xx</sub><sup>(ij)</sup></td><td>(&sigma;<sub>x</sub><sup>(i)</sup> &otimes; &sigma;<sub>y</sub><sup>(j)</sup>, {y, z}); (&sigma;<sub>y</sub><sup>(i)</sup> &otimes; &sigma;<sub>y</sub><sup>(j)</sup>, {x, x})</td><td>(28)</td></tr><tr><td>a<sub>yy</sub><sup>(ij)</sup></td><td>(&sigma;<sub>x</sub><sup>(i)</sup> &otimes; &sigma;<sub>y</sub><sup>(j)</sup>, {z, z}); (&sigma;<sub>x</sub><sup>(i)</sup> &otimes; &sigma;<sub>y</sub><sup>(j)</sup>, {y, z}); (&sigma;<sub>y</sub><sup>(i)</sup> &otimes; &sigma;<sub>y</sub><sup>(j)</sup>, {x, x})</td><td>(28)</td></tr><tr><td>a<sub>zz</sub><sup>(ij)</sup></td><td>(&sigma;<sub>x</sub><sup>(i)</sup> &otimes; &sigma;<sub>z</sub><sup>(j)</sup>, {y, x}); (&sigma;<sub>y</sub><sup>(i)</sup> &otimes; &sigma;<sub>y</sub><sup>(j)</sup>, {z, z})</td><td>(30)</td></tr><tr><td>a<sub>xy</sub><sup>(ij)</sup></td><td>(&sigma;<sub>y</sub><sup>(i)</sup> &otimes; &sigma;<sub>y</sub><sup>(j)</sup>, {z, z})</td><td>(26)</td></tr><tr><td>a<sub>yz</sub><sup>(ij)</sup></td><td>(&sigma;<sub>y</sub><sup>(i)</sup> &otimes; &sigma;<sub>y</sub><sup>(j)</sup>, {x, x})</td><td>(26)</td></tr><tr><td>a<sub>zx</sub><sup>(ij)</sup></td><td>(&sigma;<sub>x</sub><sup>(i)</sup> &otimes; &sigma;<sub>y</sub><sup>(j)</sup>, {y, x}); (&sigma;<sub>y</sub><sup>(i)</sup> &otimes; &sigma;<sub>y</sub><sup>(j)</sup>, {x, x})</td><td>(28)</td></tr></tbody></table>

**Supplementary Table I.** The first column represents the type of the estimated Hamiltonian parameters $a_{\alpha_i}^{(i)}$, $a_{\alpha_i, \alpha_j}^{(ij)}$, $\alpha_i$, $\alpha_j \in \{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 $\rho_{\tau_i, \tau_j}^{(i,j)} = \rho_{\tau_i}^{(i)} \otimes \rho_{\tau_j}^{(j)} \otimes \rho_{n-2}$, $\tau_i, \tau_j = \{x, y, z\}$, given in the second column. For example, to estimate $a_x^{(i)}$ we use (37) with two pairs of observables and states: $\{O, \rho_{\tau_i, \tau_j}^{(i,j)}\} = \{(\sigma_y^{(i)} \otimes \sigma_\eta^{(j)}, \rho_{x,\eta}^{(i,j)}); (\sigma_y^{(i)}, \rho_{x,\eta}^{(i,j)})\}$.

For an observable $O = \sigma_\gamma^{(i)} \otimes \sigma_\eta^{(j)} \otimes \sigma_{\alpha_l}^{(l)}$ we can write the equation

$$ \frac{d}{dt} \mathrm{tr} \left[ \rho_t (\sigma_\gamma^{(i)} \otimes \sigma_\eta^{(j)} \otimes \sigma_{\alpha_l}^{(l)}) \right] \Big|_{t=0} = -4 \sum_{\alpha_i \alpha_j \alpha_l} a_{\alpha_i \alpha_j \alpha_l}^{(ijl)} (\varepsilon_{\alpha_i \tau_i \gamma} \delta_{\eta \alpha_j} + \varepsilon_{\alpha_j \tau_j \eta} \delta_{\gamma \alpha_i}) \\ + 2i \sum_{\mu, \nu} L_{\mu\nu}^{(l)} \varepsilon_{\mu\nu\gamma} \delta_{\eta\alpha_l} \delta_{\tau_i\gamma} \delta_{\tau_j\eta}. \quad (45) $$

Selecting $\tau_j = \eta$, $\tau_i \neq \gamma$, we can rewrite (45) as

$$ \frac{d}{dt} \mathrm{tr} \left[ \rho_t (\sigma_\gamma^{(i)} \otimes \sigma_\eta^{(j)} \otimes \sigma_{\alpha_l}^{(l)}) \right] \Big|_{t=0} = -4 \sum_{\alpha_i \alpha_j \alpha_l} a_{\alpha_i \alpha_j \alpha_l}^{(ijl)} \varepsilon_{\alpha_i \tau_i \gamma} \delta_{\eta \alpha_j}. \quad (46) $$
From this equation we can find $a_{\alpha_i \alpha_j \alpha_l}^{(ijl)}$. For an observable $O = \sigma_{\alpha_1}^{(1)} \otimes \dots \otimes \sigma_{\gamma}^{(i)} \dots \otimes \sigma_{\eta}^{(j)} \dots \otimes \sigma_{\alpha_k}^{(k)}$, $k > 3$ we can write

$$ \frac{d}{dt} \mathrm{tr} \left[ \rho_t (\sigma_{\alpha_1}^{(1)} \otimes \dots \otimes \sigma_{\gamma}^{(i)} \dots \otimes \sigma_{\eta}^{(j)} \dots \otimes \sigma_{\alpha_k}^{(k)}) \right] \Big|_{t=0} = -2^{k-1} \sum_{\alpha_1, \dots, \alpha_k} a_{\alpha_1 \dots \alpha_k}^{(1 \dots k)} (\varepsilon_{\alpha_i \tau_i \gamma} \delta_{\eta \alpha_j} + \varepsilon_{\alpha_j \tau_j \eta} \delta_{\gamma \alpha_i}). \quad (47) $$

Selecting $\tau_j = \eta$, $\tau_i \neq \gamma$, we can rewrite (47) as

$$ \frac{d}{dt} \mathrm{tr} \left[ \rho_t (\sigma_{\alpha_1}^{(1)} \otimes \dots \otimes \sigma_{\gamma}^{(i)} \dots \otimes \sigma_{\eta}^{(j)} \dots \otimes \sigma_{\alpha_k}^{(k)}) \right] \Big|_{t=0} = -2^{k-1} \sum_{\alpha_1, \dots, \alpha_k} a_{\alpha_1 \dots \alpha_k}^{(1 \dots k)} \varepsilon_{\alpha_i \tau_i \gamma} \delta_{\eta \alpha_j}. \quad (48) $$

Solving this equation, we find $a_{\alpha_1 \dots \alpha_k}^{(1 \dots k)}$.

From (15) we can find three Lindbladian coefficients. Selecting $\tau_i = y$, $\tau_j = \eta$, $\gamma = x$, we can find

$$ L_{xy}^{(i)} = \frac{1}{2} \frac{d}{dt} \mathrm{tr} \left[ \rho_{ty\eta} (\sigma_x^{(i)} \otimes \sigma_\eta^{(j)}) \right] \Big|_{t=0} - a_z^{(i)} - a_{z\eta}^{(ij)}. \quad (49) $$

Next, for $\tau_i = z$, $\tau_j = \eta$, $\gamma = x$, we deduce

$$ L_{xz}^{(i)} = \frac{1}{2} \frac{d}{dt} \mathrm{tr} \left[ \rho_{tz\eta} (\sigma_x^{(i)} \otimes \sigma_\eta^{(j)}) \right] \Big|_{t=0} + a_y^{(i)} + a_{y\eta}^{(ij)}. \quad (50) $$

Finally, for $\tau_i = z$, $\tau_j = \eta$, $\gamma = y$, the coefficient is

$$ L_{yz}^{(i)} = \frac{1}{2} \frac{d}{dt} \mathrm{tr} \left[ \rho_{tz\eta} (\sigma_y^{(i)} \otimes \sigma_\eta^{(j)}) \right] \Big|_{t=0} - a_x^{(i)} - a_{x\eta}^{(ij)}. \quad (51) $$

From (23) we can find three more coefficients. Selecting $\tau_j = y$, $\tau_i = z$, $\gamma = z$, $\eta = x$, we get

$$ L_{xy}^{(j)} = \frac{1}{2} \frac{d}{dt} \mathrm{tr} \left[ \rho_{tzy} (\sigma_z^{(i)} \otimes \sigma_x^{(j)}) \right] \Big|_{t=0} - a_z^{(j)} - a_{zz}^{(ij)}. \quad (52) $$

For $\tau_j = y$, $\tau_i = y$, $\gamma = y$, $\eta = z$, we can deduce

$$ L_{zy}^{(j)} = \frac{1}{2} \frac{d}{dt} \mathrm{tr} \left[ \rho_{tyy} (\sigma_y^{(i)} \otimes \sigma_z^{(j)}) \right] \Big|_{t=0} + a_x^{(j)} + a_{yx}^{(ij)}. \quad (53) $$

Selecting $\tau_j = z$, $\tau_i = y$, $\gamma = y$, $\eta = x$, we get

$$ L_{xz}^{(j)} = \frac{1}{2} \frac{d}{dt} \mathrm{tr} \left[ \rho_{tyz} (\sigma_y^{(i)} \otimes \sigma_x^{(j)}) \right] \Big|_{t=0} + a_y^{(j)} + a_{yy}^{(ij)}. \quad (54) $$

From (33) we find the following coefficients:

$$ L_{xx}^{(i)} = \frac{1}{10} \left( -3 \frac{d}{dt} \mathrm{tr} \left[ \rho_{tx\tau_j} \sigma_x^{(i)} \right] \Big|_{t=0} - 3 \frac{d}{dt} \mathrm{tr} \left[ \rho_{ty\tau_j} \sigma_y^{(i)} \right] \Big|_{t=0} + 2 \frac{d}{dt} \mathrm{tr} \left[ \rho_{tz\tau_j} \sigma_z^{(i)} \right] \Big|_{t=0} \right), \quad (55) $$

$$ L_{yy}^{(i)} = \frac{1}{10} \left( -3 \frac{d}{dt} \mathrm{tr} \left[ \rho_{tx\tau_j} \sigma_x^{(i)} \right] \Big|_{t=0} + 2 \frac{d}{dt} \mathrm{tr} \left[ \rho_{ty\tau_j} \sigma_y^{(i)} \right] \Big|_{t=0} - 3 \frac{d}{dt} \mathrm{tr} \left[ \rho_{tz\tau_j} \sigma_z^{(i)} \right] \Big|_{t=0} \right), $$

$$ L_{zz}^{(i)} = \frac{1}{10} \left( 2 \frac{d}{dt} \mathrm{tr} \left[ \rho_{tx\tau_j} \sigma_x^{(i)} \right] \Big|_{t=0} - 3 \frac{d}{dt} \mathrm{tr} \left[ \rho_{ty\tau_j} \sigma_y^{(i)} \right] \Big|_{t=0} - 3 \frac{d}{dt} \mathrm{tr} \left[ \rho_{tz\tau_j} \sigma_z^{(i)} \right] \Big|_{t=0} \right). $$

From (34) the following coefficients can be found

$$ L_{xx}^{(j)} = \frac{1}{10} \left( -3 \frac{d}{dt} \mathrm{tr} \left[ \rho_{t\gamma_1 x} \sigma_x^{(i)} \right] \Big|_{t=0} - 3 \frac{d}{dt} \mathrm{tr} \left[ \rho_{t\gamma_2 y} \sigma_y^{(i)} \right] \Big|_{t=0} + 2 \frac{d}{dt} \mathrm{tr} \left[ \rho_{t\gamma_3 z} \sigma_z^{(i)} \right] \Big|_{t=0} \right), \quad (56) $$

$$ L_{yy}^{(j)} = \frac{1}{10} \left( -3 \frac{d}{dt} \mathrm{tr} \left[ \rho_{t\gamma_1 x} \sigma_x^{(i)} \right] \Big|_{t=0} + 2 \frac{d}{dt} \mathrm{tr} \left[ \rho_{t\gamma_2 y} \sigma_y^{(i)} \right] \Big|_{t=0} - 3 \frac{d}{dt} \mathrm{tr} \left[ \rho_{t\gamma_3 z} \sigma_z^{(i)} \right] \Big|_{t=0} \right), $$

$$ L_{zz}^{(j)} = \frac{1}{10} \left( 2 \frac{d}{dt} \mathrm{tr} \left[ \rho_{t\gamma_1 x} \sigma_x^{(i)} \right] \Big|_{t=0} - 3 \frac{d}{dt} \mathrm{tr} \left[ \rho_{t\gamma_2 y} \sigma_y^{(i)} \right] \Big|_{t=0} - 3 \frac{d}{dt} \mathrm{tr} \left[ \rho_{t\gamma_3 z} \sigma_z^{(i)} \right] \Big|_{t=0} \right). $$

All the Lindbladian coefficients with the corresponding observables and initial states are given in Table II.
<table><thead><tr><th>L<sub>μν</sub><sup>(i)</sup></th><th>(O, {τ<sub>i</sub><sup>(i)</sup>, τ<sub>j</sub><sup>(j)</sup>})</th><th>Equation</th></tr></thead><tbody><tr><td>L<sub>xx</sub><sup>(i)</sup>, L<sub>yy</sub><sup>(i)</sup>, L<sub>zz</sub><sup>(i)</sup></td><td>(σ<sub>x</sub><sup>(i)</sup>, {x, τ<sub>j1</sub>}); (σ<sub>y</sub><sup>(i)</sup>, {y, τ<sub>j2</sub>}); (σ<sub>z</sub><sup>(i)</sup>, {z, τ<sub>j3</sub>}), τ<sub>j1,2,3</sub> ∈ {x, y, z}</td><td>(55)</td></tr><tr><td>L<sub>xx</sub><sup>(j)</sup>, L<sub>yy</sub><sup>(j)</sup>, L<sub>zz</sub><sup>(j)</sup></td><td>(σ<sub>x</sub><sup>(i)</sup>, {γ<sub>1</sub>, x}); (σ<sub>y</sub><sup>(i)</sup>, {γ<sub>2</sub>, y}); (σ<sub>z</sub><sup>(i)</sup>, {γ<sub>3</sub>, z}), γ<sub>1</sub> ∈ {y, z}, γ<sub>2</sub> ∈ {x, z}, γ<sub>3</sub> ∈ {x, y}</td><td>(56)</td></tr><tr><td>L<sub>xy</sub><sup>(i)</sup></td><td>(σ<sub>x</sub><sup>(i)</sup> ⊗ σ<sub>η</sub><sup>(j)</sup>, {y, η}), η ∈ {x, y, z}</td><td>(49)</td></tr><tr><td>L<sub>xz</sub><sup>(i)</sup></td><td>(σ<sub>x</sub><sup>(i)</sup> ⊗ σ<sub>η</sub><sup>(j)</sup>, {z, η}), η ∈ {x, y, z}</td><td>(50)</td></tr><tr><td>L<sub>yz</sub><sup>(i)</sup></td><td>(σ<sub>y</sub><sup>(i)</sup> ⊗ σ<sub>η</sub><sup>(j)</sup>, {z, η}), η ∈ {x, y, z}</td><td>(51)</td></tr><tr><td>L<sub>xy</sub><sup>(j)</sup></td><td>(σ<sub>z</sub><sup>(i)</sup> ⊗ σ<sub>x</sub><sup>(j)</sup>, {z, y})</td><td>(52)</td></tr><tr><td>L<sub>xz</sub><sup>(j)</sup></td><td>(σ<sub>y</sub><sup>(i)</sup> ⊗ σ<sub>x</sub><sup>(j)</sup>, {y, z})</td><td>(54)</td></tr><tr><td>L<sub>yz</sub><sup>(j)</sup></td><td>(σ<sub>y</sub><sup>(i)</sup> ⊗ σ<sub>z</sub><sup>(j)</sup>, {y, y})</td><td>(53)</td></tr></tbody></table>

**Supplementary Table II.** The first column represents the type of the estimated Lindbladian parameters $L_{\mu\nu}^{(i)}$, $\mu, \nu \in \{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 $\rho_{\tau_i, \tau_j}^{(i,j)} = \rho_{\tau_i}^{(i)} \otimes \rho_{\tau_j}^{(j)} \otimes \rho_{n-2}$, $\tau_i, \tau_j = \{x, y, z\}$, given in the second column. For example, to estimate $L_{yz}^{(j)}$ we use (53) with one pair of observable and state: $\{O, \rho_{\tau_i, \tau_j}^{(i,j)}\} = \{\sigma_y^{(i)} \otimes \sigma_z^{(j)}, \rho_{y,y}^{(i,j)}\}$.

## 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 $|\psi(t)\rangle$ 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 = \sum_{m_1,m_2=1}^{n} \sum_{\alpha_1, \alpha_2=x,y} J_{m_1,m_2} \sigma_{\alpha_1}^{(m_1)} \sigma_{\alpha_2}^{(m_2)} + \frac{1}{2} \sum_{m_1=1}^{n} \Omega_{m_1} \sigma_z^{(m_1)}, \quad (57)$$

For the problem considered in this paper the couplings $J_{m_1,m_2}$ are restricted to the nearest neighbor qubits on a 2-D lattice. Note that the actual frequency $\Omega_{m_1}$ of the $m_1$-th qubit in the Hamiltonian (57) is different from its nominal frequency $\tilde{\omega}_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 \quad (58)$$

$$H_X = \sum_{m_1, m_2=1}^{n} J_{m_1, m_2} \sigma_x^{(m_1)} \sigma_x^{(m_2)}, \quad H_Y = \sum_{m_1, m_2=1}^{n} J_{m_1, m_2} \sigma_y^{(m_1)} \sigma_y^{(m_2)}, \quad H_Z = \frac{1}{2} \sum_{m_1=1}^{n} \Omega_{m_1} \sigma_z^{(m_1)}, \quad (59)$$

and the corresponding Suzuki-Trotter decomposition of the evolution operator $U(\Delta t)$ for the (small) timestep of duration $\Delta t$ has the form

$$U(\Delta t) \equiv \exp(-iH\Delta t) \approx e^{-iH_Z\Delta t/2} e^{-iH_Y\Delta t/2} e^{-iH_X\Delta t} e^{-iH_Y\Delta t/2} e^{-iH_Z\Delta t/2}, \quad (60)$$

ensuring the overall time discretization error of the order $(\Delta t)^2$. The evolution operator over many time steps is a product of elementary operators $U(\Delta 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\Delta t) = \prod_{m_1,m_2=1}^{n} \exp(-iJ_{m_1,m_2}\Delta t \sigma_x^{(m_1)}\sigma_x^{(m_2)}) = \prod_{m_1,m_2=1}^{n} \left[ \cos(J_{m_1,m_2}\Delta t) - i\sigma_x^{(m_1)}\sigma_x^{(m_2)} \sin(J_{m_1,m_2}\Delta t) \right]. \quad (61) $$

Each term in this direct product acts on the wavefunction $|\psi(t)\rangle$ 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(\Delta t)$ is easy to compute, without the need to calculate or store $2^N \times 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 \times 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\rangle$, while the rest of the system is in completely mixed state, i.e. when the system's density matrix is

$$ \rho = |0\rangle\langle 0| \otimes |0\rangle\langle 0| \otimes \frac{I_{n-2}}{2^{n-2}}, \quad (62) $$

where $I_{n-2}$ is an identity matrix of the size $2^{n-2} \times 2^{n-2}$, is represented using the total wavefunction in the form

$$ |\psi\rangle = |0\rangle \otimes |0\rangle \otimes |\psi_{n-2}^{(r)}\rangle, \quad (63) $$

where the random state $|\psi_{n-2}^{(r)}\rangle$ 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 $\sim 1/\sqrt{M} \approx 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,m_1}^*$, caused by its random static frequency shift, is modeled by directly reproducing the underlying physical picture. Namely, we assumed that the actual frequency $\Omega_{m_1}$ of the $m_1$-th qubit, see Eq. 57, is a sum of two contributions: the nominal value $\tilde{\omega}_{m_1}$, and a random shift $\beta_{m_1}$ that remains constant during the system's evolution. The values of $\beta_{m_1}$ were independently sampled from Gaussian distributions with zero mean and variance $b_{m_1}^2$, which can be different for different qubits. The parameter $b_{m_1}$ determines the dephasing time $T_{2,m_1}^*$ of the $m_1$-th qubit: if this qubit were uncoupled from the rest of the system, then, after averaging over $\beta_{m_1}$, its transverse (x- and y-) components would undergo Gaussian decay with time dependence $\exp(-b_{m_1}^2 t^2/2)$, i.e. $b_{m_1} = \sqrt{2}/T_{2,m_1}^*$.

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 $\beta_{m_1}$ (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 $\beta_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 $\beta_{m_1}$ (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 $\rho(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 $\beta_{m_1}(t)$ representing the random frequency shifts are, of course, Markovian random processes, satisfying 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 $\rho(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 $\Omega_{m_1}$ (i.e., produced by the action of the operator $\exp(-iH_Z\Delta t)$ in Eq. 60), and adding to it another time-dependent rotation around the z-axis by the angle $\Lambda_{m_1}(t)$. For each time step of duration $\Delta t$, the values $\Lambda_{m_1}(t)$ were sampled randomly, indepedently of each other and of their previous values, from Gaussian distribution with zero mean and variance $4g_{m_1}^2\tau_{m_1}\Delta t$. This choice for the quantity $\Lambda(t)$ can be visualized as a rotation induced by a time-dependent frequency shift $\delta\omega_{m_1}(t)$, which is represented by an Ornstein-Uhlenbeck noise process with the correlation function $\langle\delta\omega_{m_1}(t)\delta\omega_{m_1}(t+s)\rangle = g_{m_1}^2 \exp(-|s|/\tau_{m_1})$, in the limit where the correlation time $\tau_{m_1}$ is much smaller than $\Delta t$, while the magnitude $g_{m_1}$ is large (formally, $\tau_{m_1} \to 0$ and $g_{m_1} \to \infty$), but the combination $g_{m_1}^2\tau_{m_1} = 1/T_{2,m_1}$ remains finite. For an isolated qubit, the average evolution under the influence of such Ornstein-Uhlenbeck noise $\delta\omega_{m_1}(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 $\Lambda_{m_1}(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 $\Delta t \ll T_{2,m_1}$, 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\rangle$ on a timescale $T_{1,m_1}$, 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_{m_1}$ for the $m_1$-th qubit to make a transition (“quantum jump”) from the state $|1\rangle$ to the state $|0\rangle$; the corresponding value is $p_{m_1} = w_{m_1,1}\mu_{m_1}^2$, where $\mu_{m_1} = \sqrt{1 - \exp(-\Delta t/T_{1,m_1})}$, and $w_{m_1,1}$ is the total probability of the system to be in the subspace corresponding to the j-th qubit in the state $|1\rangle$. This transition was implemented with the probability $p_{m_1}$ at each time step: the part of the system's wavefunction corresponding to the $m_1$-th qubit in the state $|0\rangle$ was replaced by its complement, i.e. by the part corresponding to the $m_1$-th qubit in the state $|1\rangle$, multiplied by the factor $\mu_{m_1}$, and the part of the wavefunction corresponding to the j-th qubit in the state $|1\rangle$ was set to zero. Alternatively, with the probability $1-p_{m_1}$ at each time step, the part of the wavefunction corresponding to the $m_1$-th qubit in the state $|1\rangle$ was multiplied by the factor $\exp[-(1/2)\Delta t/T_{1,m_1}]$, while its complement was left unchanged. These transformations were applied to the wavefunction in succession, for all qubits (for all $m_1 = 1, \dots, 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\Delta t)$ in Eq. 60, they were applied at the end of each unitary time-step evolution, after application of the operator $U(\Delta t)$ given by Eq. 60, in parallel with the action of the operators $\exp(-iH_Z\Delta t)$ or $\exp(-iH_Z\Delta 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 = \sum_{k=i,j} a_z^{(k)} H_z^{(k)} + a_{xx}^{(i,j)} H_{xx}^{(i,j)} + a_{yy}^{(i,j)} H_{yy}^{(i,j)}, \qquad (64)$$

where we introduce the notations

$$a_z^{(i)} = \frac{1}{2} \tilde{\omega}^{(i)}, \quad a_z^{(j)} = \frac{1}{2} \tilde{\omega}^{(j)}, \quad a_{xx}^{(i,j)} = a_{yy}^{(i,j)} = \frac{1}{2} \left[ \frac{g_i g_j}{\Delta} + g_{ij} \right], \qquad (65)$$

$$H_z^{(i)} = \sigma_z^{(i)}, \quad H_z^{(j)} = \sigma_z^{(j)}, \quad H_{xx}^{(i,j)} = \sigma_x^{(i)} \sigma_x^{(j)}, \quad H_{yy}^{(i,j)} = \sigma_y^{(i)} \sigma_y^{(j)}.$$
Thus, we define the 16 qubits 2D grid, where we generate $\omega_1, \omega_2, \omega_c, g_1, g_2, g_{12}$ from Gaussian distribution with mean and variance $N(0, 1)$. The parameters $a_z^{(i)}, a_z^{(j)}, a_{xx}^{(ij)}, a_{yy}^{(ij)}$ 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_z^{(i)}, a_z^{(j)}, a_{xx}^{(i,j)}$

<table><thead><tr><th>(ij)</th><th>a<sub>xx</sub><sup>(ij)</sup>, a<sub>yy</sub><sup>(ij)</sup> [kHz]</th><th>(i)</th><th>a<sub>z</sub><sup>(i)</sup> [kHz]</th><th>T<sub>1</sub> [µs]</th><th>T<sub>2</sub> [µs]</th><th>T<sub>2</sub><sup>*</sup> [µs]</th></tr></thead><tbody><tr><td>1, 2</td><td>1.28112</td><td>1</td><td>1.73807</td><td>58.5227</td><td>65.9752</td><td>151.515</td></tr><tr><td>2, 3</td><td>-0.716875</td><td>2</td><td>-0.816877</td><td>60.0269</td><td>65.1704</td><td>166.667</td></tr><tr><td>3, 4</td><td>-0.956949</td><td>3</td><td>-1.0602</td><td>59.2424</td><td>64.6375</td><td>163.934</td></tr><tr><td>4, 5</td><td>-0.819328</td><td>4</td><td>-0.913223</td><td>61.0255</td><td>65.7397</td><td>149.254</td></tr><tr><td>1, 6</td><td>-1.1682</td><td>5</td><td>-1.23118</td><td>59.0545</td><td>66.0886</td><td>147.059</td></tr><tr><td>2, 7</td><td>-0.213057</td><td>6</td><td>-0.654699</td><td>60.0915</td><td>66.1118</td><td>151.515</td></tr><tr><td>3, 8</td><td>-0.563789</td><td>7</td><td>-0.514756</td><td>59.8856</td><td>65.1432</td><td>153.846</td></tr><tr><td>4, 9</td><td>1.74022</td><td>8</td><td>2.0817</td><td>61.0389</td><td>64.8252</td><td>158.73</td></tr><tr><td>5, 10</td><td>1.68348</td><td>9</td><td>-0.568581</td><td>60.5375</td><td>66.2155</td><td>158.73</td></tr><tr><td>6, 7</td><td>-1.51535</td><td>10</td><td>-0.710498</td><td>61.5036</td><td>65.389</td><td>149.254</td></tr><tr><td>7, 8</td><td>-0.729672</td><td>11</td><td>1.86153</td><td>59.8949</td><td>65.825</td><td>147.059</td></tr><tr><td>8, 9</td><td>1.6622</td><td>12</td><td>-2.03725</td><td>60.3777</td><td>65.1203</td><td>153.846</td></tr><tr><td>9, 10</td><td>-0.314438</td><td>13</td><td>-1.31695</td><td>57.5781</td><td>65.6052</td><td>156.25</td></tr><tr><td>11, 12</td><td>-0.475787</td><td>14</td><td>-0.902159</td><td>59.1881</td><td>65.8892</td><td>158.73</td></tr><tr><td>6, 11</td><td>1.3663</td><td>15</td><td>-0.202118</td><td>60.0283</td><td>66.2967</td><td>144.928</td></tr><tr><td>7, 12</td><td>-2.03531</td><td>16</td><td>0.136975</td><td>58.9397</td><td>66.0541</td><td>149.254</td></tr><tr><td>12, 13</td><td>-1.22632</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>13, 8</td><td>-0.717182</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>13, 14</td><td>-0.546421</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>14, 9</td><td>1.90836</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>14, 15</td><td>-0.781306</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>11, 16</td><td>-0.358714</td><td></td><td></td><td></td><td></td><td></td></tr></tbody></table>

**Supplementary Table III.** The first column represents the numbers of qubits whose coefficients $a_{\alpha_i, \alpha_j}^{(ij)}$, $\alpha_i, \alpha_j \in \{x, y, z\}$ are not zero, given in the second column. The fourth column contains $a_{\alpha_i}^{(i)}, \alpha_i \in \{x, y, z\}, i = 1, \dots, 16$. The rates of decay of $T_1, T_2$ and $T_2^*$, corresponding to *i*'s qubit are given in the fifth and six's columns, respectively.

and $a_{yy}^{(i,j)}$, are given in Table IV. One can see, that we need three starting states, namely $\rho_{01} = \rho_x^{(i)} \otimes \rho_x^{(j)} \otimes \frac{I^{2n-2}}{2^{2n-2}}$, $\rho_{02} = \rho_y^{(i)} \otimes \rho_z^{(j)} \otimes \frac{I^{2n-2}}{2^{2n-2}}$ and $\rho_{03} = \rho_z^{(i)} \otimes \rho_z^{(j)} \otimes \frac{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_z^{(i)}, a_z^{(j)}, a_{xx}^{(i,j)}, a_{yy}^{(i,j)}$ 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.

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 $\rho_{04} = \rho_z^{(i)} \otimes \rho_y^{(j)} \otimes \frac{I^{2n-2}}{2^{2n-2}}$, $\rho_{05} = \rho_y^{(i)} \otimes \rho_y^{(j)} \otimes \frac{I^{2n-2}}{2^{2n-2}}$ and $\rho_{06} = \rho_z^{(i)} \otimes \rho_x^{(j)} \otimes \frac{I^{2n-2}}{2^{2n-2}}$ to find $L_{\mu\nu}^{(i)}, \mu, \nu \in \{x, y, z\}$.

<table><thead><tr><th>a<sub>α<sub>i</sub></sub><sup>(i)</sup></th><th>{O, ρ<sub>τ<sub>i</sub>,τ<sub>j</sub></sub><sup>(i,j)</sup>}</th><th>Equation</th></tr></thead><tbody><tr><td>a<sub>z</sub><sup>(i)</sup></td><td>{σ<sub>x</sub><sup>(i)</sup> ⊗ σ<sub>z</sub><sup>(j)</sup>, ρ<sub>y,z</sub><sup>(i,j)</sup>}; {σ<sub>x</sub><sup>(i)</sup>, ρ<sub>y,z</sub><sup>(i,j)</sup>}</td><td>(39)</td></tr><tr><td>a<sub>z</sub><sup>(j)</sup></td><td>{σ<sub>x</sub><sup>(i)</sup> ⊗ σ<sub>y</sub><sup>(j)</sup>, ρ<sub>x,x</sub><sup>(i,j)</sup>}; {σ<sub>y</sub><sup>(i)</sup>, ρ<sub>x,x</sub><sup>(i,j)</sup>}</td><td>(44)</td></tr><tr><td>a<sub>xx</sub><sup>(ij)</sup></td><td>{σ<sub>x</sub><sup>(i)</sup> ⊗ σ<sub>y</sub><sup>(j)</sup>, ρ<sub>y,z</sub><sup>(i,j)</sup>}; {σ<sub>y</sub><sup>(i)</sup> ⊗ σ<sub>y</sub><sup>(j)</sup>, ρ<sub>x,x</sub><sup>(i,j)</sup>}</td><td>(28)</td></tr><tr><td>a<sub>yy</sub><sup>(ij)</sup></td><td>{σ<sub>x</sub><sup>(i)</sup> ⊗ σ<sub>y</sub><sup>(j)</sup>, ρ<sub>z,z</sub><sup>(i,j)</sup>}; {σ<sub>x</sub><sup>(i)</sup> ⊗ σ<sub>y</sub><sup>(j)</sup>, ρ<sub>y,z</sub><sup>(i,j)</sup>}; {σ<sub>y</sub><sup>(i)</sup> ⊗ σ<sub>y</sub><sup>(j)</sup>, ρ<sub>x,x</sub><sup>(i,j)</sup>}</td><td>(28)</td></tr></tbody></table>

**Supplementary Table IV.** In this table the minimal selection of the pairs $\{O, \rho_{\tau_i, \tau_j}^{(i,j)}\}$ is presented for our specific example. The first column represents the type of the estimated Hamiltonian (64) parameters $a_{\alpha_i}^{(i)}, a_{\alpha_i, \alpha_j}^{(ij)}, \alpha_i, \alpha_j \in \{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 $\rho_{\tau_i, \tau_j}^{(i,j)} = \rho_{\tau_i}^{(i)} \otimes \rho_{\tau_j}^{(j)} \otimes \rho_{N-2, N-2}$, $\tau_i, \tau_j = \{x, y, z\}$, given in the second column. To estimate all four parameters we need only three initial states: $\rho_{x,x}^{(i,j)}, \rho_{z,z}^{(i,j)}$ and $\rho_{y,z}^{(i,j)}$.
<table><thead><tr><th>L<sub>μν</sub><sup>(i)</sup></th><th>{O, ρ<sub>τ<sub>i</sub>,τ<sub>j</sub></sub><sup>(i,j)</sup>}</th><th>Equation</th></tr></thead><tbody><tr><td>L<sub>xx</sub><sup>(i)</sup>, L<sub>yy</sub><sup>(i)</sup>, L<sub>zz</sub><sup>(i)</sup></td><td>{σ<sub>x</sub><sup>(i)</sup>, ρ<sub>x,x</sub><sup>(i,j)</sup>}; {σ<sub>y</sub><sup>(i)</sup>, ρ<sub>y,z</sub><sup>(i,j)</sup>}; {σ<sub>z</sub><sup>(i)</sup>, ρ<sub>z,z</sub><sup>(i,j)</sup>}</td><td>(55)</td></tr><tr><td>L<sub>xx</sub><sup>(j)</sup>, L<sub>yy</sub><sup>(j)</sup>, L<sub>zz</sub><sup>(j)</sup></td><td>{σ<sub>x</sub><sup>(i)</sup>, ρ<sub>z,x</sub><sup>(i,j)</sup>}; {σ<sub>y</sub><sup>(i)</sup>, ρ<sub>z,y</sub><sup>(i,j)</sup>}; {σ<sub>z</sub><sup>(i)</sup>, ρ<sub>y,z</sub><sup>(i,j)</sup>}</td><td>(56)</td></tr><tr><td>L<sub>xy</sub><sup>(i)</sup></td><td>{σ<sub>x</sub><sup>(i)</sup> ⊗ σ<sub>z</sub><sup>(j)</sup>, ρ<sub>y,z</sub><sup>(i,j)</sup>}</td><td>(49)</td></tr><tr><td>L<sub>xz</sub><sup>(i)</sup></td><td>{σ<sub>x</sub><sup>(i)</sup> ⊗ σ<sub>η</sub><sup>(j)</sup>, ρ<sub>z,z</sub><sup>(i,j)</sup>}</td><td>(50)</td></tr><tr><td>L<sub>yz</sub><sup>(i)</sup></td><td>{σ<sub>y</sub><sup>(i)</sup> ⊗ σ<sub>η</sub><sup>(j)</sup>, ρ<sub>z,z</sub><sup>(i,j)</sup>}</td><td>(51)</td></tr><tr><td>L<sub>xy</sub><sup>(j)</sup></td><td>{σ<sub>z</sub><sup>(i)</sup> ⊗ σ<sub>x</sub><sup>(j)</sup>, ρ<sub>z,y</sub><sup>(i,j)</sup>}</td><td>(52)</td></tr><tr><td>L<sub>xz</sub><sup>(j)</sup></td><td>{σ<sub>y</sub><sup>(i)</sup> ⊗ σ<sub>x</sub><sup>(j)</sup>, ρ<sub>y,z</sub><sup>(i,j)</sup>}</td><td>(54)</td></tr><tr><td>L<sub>yz</sub><sup>(j)</sup></td><td>{σ<sub>y</sub><sup>(i)</sup> ⊗ σ<sub>z</sub><sup>(j)</sup>, ρ<sub>y,y</sub><sup>(i,j)</sup>}</td><td>(53)</td></tr></tbody></table>

**Supplementary Table V.** In this table the minimal selection of the pairs {$O, \rho_{\tau_i, \tau_j}^{(i,j)}$} is presented for our specific example. The first column represents the type of the estimated Lindbladian parameters $L_{\mu\nu}^{(i)}$, $\mu, \nu \in \{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 $\rho_{\tau_i, \tau_j}^{(i,j)} = \rho_{\tau_i}^{(i)} \otimes \rho_{\tau_j}^{(j)} \otimes \rho_{N-2}$, $\tau_i, \tau_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: $\rho_{z,x}^{(i,j)}, \rho_{z,y}^{(i,j)}$ and $\rho_{y,y}^{(i,j)}$.

**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 $\Gamma$, we let $\mathcal{L}_{\Lambda}: \mathcal{M}_{2^n} \to \mathcal{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:

$$ \mathcal{L}_{\Lambda} = \sum_{A \subset \Lambda} \mathcal{L}_{A}, \qquad (66) $$

where $\mathcal{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 $\mathcal{L}$ is $k$-local if $\mathcal{L}_A \neq 0$ only if $A$ is a subset of vertices of $G$ containing at most $k$ vertices. Furthermore, we will say that $\mathcal{L}_{\Lambda}$ is locally bounded if there is a constant $g > 0$ such that for all $B \subset \Lambda$ we have that:

$$ \| \sum_{A \subset \Lambda: A \cap B \neq \emptyset} \mathcal{L}_A \| \leq g|B|. \qquad (67) $$

This condition is satisfied if e.g. $\mathcal{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+\delta)}$ for some $\delta > 0$. Moreover, for ease of notation we will let for a region $B \subset \Lambda$

$$ \mathcal{L}_B = \sum_{A \subset B} \mathcal{L}_A \qquad (68) $$

be the generator restricted or truncated to a subregion $B$.

Furthermore, given the lattice $\Lambda$, some region $X \subset V$ and $r > 0$, we will denote by $\Lambda_r(X)$ the set of vertices that are a distance at most $r$ from $X$:

$$ \Lambda_r(X) = \{v \in V : \exists x \in X \text{ s.t. } \mathrm{dist}(x, v) \leq r\}. \qquad (69) $$

Here the distance is according to the lattice.

We will also require some norms for superoperators. Given a superoperator $\Phi: \mathcal{M}_{2^n} \to \mathcal{M}_{2^n}$ we define for $p, q \geq 1$

$$ \|\Phi\|_{p \to q} = \sup_{X \in \mathcal{M}_{2^n}} \frac{\|\Phi(X)\|_q}{\|X\|_p}, \qquad (70) $$

where $\|\cdot\|_p$ corresponds to the Schatten $p$-norm. Also note that $p = \infty$ corresponds to the operator norm, for which we will often drop the $+\infty$ and only write as $\|\cdot\|$. We will also consider completely bounded versions of these norms,
which are given by:

$$ \|\Phi\|_{p \to q, cb} = \sup_{d \ge 1} \sup_{X \in \mathcal{M}_{2^n} \otimes \mathcal{M}_d} \frac{\|\Phi \otimes \mathrm{id}_d(X)\|_q}{\|X\|_p}. \quad (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:

**Lemma 4.1** (Derivatives of local evolutions). Let $\mathcal{L}_\Lambda$ be a locally bounded Lindbladian with constant $g$. For an observable $O$ such that $\|O\| \le 1$, an initial state $\rho$ and a region $B \subset \Lambda$ defines $f_B : \mathbb{R}^+ \to \mathbb{R}$ as $f_B(t) = \mathrm{tr}[e^{t\mathcal{L}_B}(O)\rho]$. Then for all $t \ge 0$:

$$ |f_B^{(k)}(t)| \le (g|B|)^k. \quad (72) $$

In particular, for any $0 < t < t_{\max}$ we have that:

$$ \left| f_B(t) - \sum_{k=0}^{d} \frac{f_B^{(k)}(0)}{k!} t^k \right| \le \frac{(t_{\max}g|B|)^{d+1}}{(d+1)!}. \quad (73) $$

*Proof.* The proof is elementary. Note that:

$$ f_B^{(k)}(t) = \mathrm{tr}[e^{\mathcal{L}_B}((t\mathcal{L}_B)^k(O))\rho]. \quad (74) $$

Now, by Hölder's inequality we have that:

$$ |f_B^{(k)}(t)| \le \|e^{t\mathcal{L}_B}((t\mathcal{L}_B)^k(O))\|_{\infty} \|\rho\|_1 \stackrel{(1)}{\le} t^k \|e^{t\mathcal{L}_B}\|_{\infty \to \infty} \|\mathcal{L}_B\|_{\infty \to \infty} \stackrel{(2)}{\le} g^k |B|^k, $$

where in (1) we used the submultiplicativity of the operator norm, i.e. $\|\Phi_1\Phi_2\|_{\infty \to \infty} \le \|\Phi_1\|_{\infty \to \infty} \|\Phi_2\|_{\infty \to \infty}$ for all linear maps $\Phi_1, \Phi_2$. In (2) we used the fact that for any quantum channel $\|e^{t\mathcal{L}}\|_{\infty \to \infty} = 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 $\epsilon > 0$ and $t_{\max} > 0$ s.t.

$$ \log \log(\epsilon^{-1}) 2et_{\max}g|B| \ge 1 \quad (75) $$

there is a polynomial $p$ of degree

$$ d = 2et_{\max}g|B|\log(\epsilon^{-1}) - 1 \quad (76) $$

such that for all $0 \le t \le t_{\max}$ we have that

$$ |f_B(t) - p(t)| \le \epsilon. \quad (77) $$

*Proof.* It follows from Sitrling's approximation that the error in Eq. (73) is bounded by

$$ \left| f(t) - \sum_{k=0}^{d} \frac{f^{(k)}(0)}{k!} t^k \right| \le \frac{1}{\sqrt{2\pi d}} \left( \frac{et_{\max}g|B|}{d+1} \right)^{d+1}. \quad (78) $$

It is then easy to see that picking $d = 2et_{\max}g|B|\log(\epsilon^{-1}) - 1$ is sufficient to ensure that the error in (73) is at most $\epsilon$. Indeed, plugging in the value of $d$ into Eq. (78), we get:

$$ \left| f(t) - \sum_{k=0}^{d} \frac{f^{(k)}(0)}{k!} t^k \right| \le \frac{1}{\sqrt{2\pi d}} \left( \frac{1}{\log(\epsilon^{-1})} \right)^{2et_{\max}g|B|\log(\epsilon^{-1})} = \quad (79) $$

$$ \frac{1}{\sqrt{2\pi d}} \exp[-\log(\log(\epsilon^{-1})) \log(\epsilon^{-1}) 2et_{\max}g|B|] = \frac{1}{\sqrt{2\pi d}} \epsilon^{\log(\log(\epsilon^{-1})) 2et_{\max}g|B|} \le \epsilon. \quad (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 $\epsilon$. 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 protocol, the local evolution approximates the global one well.

Our main tool to show this approximability of expectation values are Lieb-Robinson bounds [8–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 \supset Y$ we set $R = \text{dist}(\Lambda\setminus\{B\}, Y)/k$, where $k$ is the locality of the generator, then it is shown in Lemma 7.1 that indeed

$$ \|(e^{t\mathcal{L}_\Lambda} - e^{t\mathcal{L}_{\Lambda_r(Y)}})(O_Y)\|_{\infty} \le (e^{vt} - 1)h(R), \quad (81) $$

where $h: \mathbb{R}^+ \to \mathbb{R}^+$ is a monotonically increasing function such that $\lim_{R \to +\infty} 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^{-\mu R}$ for some constant $\mu > 0$. For algebraically decaying evolutions we usually have $h(x) = x^k$ for some $k \in \mathbb{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$, $\epsilon > 0$ and $t_{\max}$ be given. Assume Eq. (81) holds for the time evolution $\mathcal{L}_\Lambda$ and a function $h$. Let $l > 0$ be given by

$$ l = h^{-1} \left( \frac{\epsilon}{e^{vt_{\max}} - 1} \right). \quad (82) $$

Then we have for $\Lambda_l(Y)$ and all $0 < t < t_{\max}$ and any initial state $\rho$ that

$$ | \text{tr} \left[ (e^{t\mathcal{L}_\Lambda} - e^{t\mathcal{L}_{\Lambda_l(Y)}})(O_Y)\rho \right] | \le \epsilon. \quad (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^{t\mathcal{L}_\Lambda} - e^{t\mathcal{L}_{\Lambda_r(Y)}})(O_Y)\|_{\infty} \le \epsilon. \quad (84) $$

□

Note that in the case of short-range systems we have that $h^{-1}(x) = \mu^{-1}\log(x)$. From now on we will suppress the terms of order $\log(\log(\epsilon^{-1}))$ or higher from the equations and denote bounds where we do this with $\tilde{O}$. Thus, combining 4.1 with Prop. 4.1 we conclude that:

**Theorem 4.1.** Let $\mathcal{L}_\Lambda$ be a locally bounded Lindbladian on a $D$-dimensional regular lattice with constant $g$. Moreover, let $t_{\max}, \epsilon > 0$ be given and $O_Y$ is an observable such that $\|O_Y\| \le 1$ and $O_Y$ is supported on a constant number of qubits. Assume that $\mathcal{L}_\Lambda$ satisfies Eq. (81). Then there is a polynomial $p$ of degree

$$ d = \tilde{O} \left[ \left( h^{-1} \left( \frac{\epsilon}{e^{vt_{\max}} - 1} \right) \right)^D t_{\max} \log(\epsilon^{-1}) \right] \quad (85) $$

such that for all $0 \le t \le t_{\max}$:

$$ |\text{tr}[e^{t\mathcal{L}_\Lambda}(O_Y)\rho] - p(t)| \le \epsilon, \quad (86) $$
and

$$p'(0) = \mathrm{tr}[\mathcal{L}_{\Lambda}(O_Y)\rho]. \qquad (87)$$

*Proof.* It follows from Prop. 4.1 that a region of radius

$$l = \tilde{O}\left[h^{-1}\left(\frac{\epsilon}{2(e^{t\max} - 1)}\right)\right] \qquad (88)$$

is enough to approximate the time evolution of $e^{t\mathcal{L}_\Lambda(O_Y)}$ up to $\epsilon/2$. If the original region $Y$ has a constant number of qubits, then for a $D$-dimensional lattice we have $|\Lambda_l(Y)| = \mathcal{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

$$\mathcal{O}\left[h^{-1}\left(\frac{\epsilon}{2(e^{t\max} - 1)}\right)^D t_{\max} \log(\epsilon^{-1})\right] \qquad (89)$$

to approximate the expectation value of the local evolution up to an error $\epsilon/2$. This concludes the proof by a triangle inequality. Eq. (87) is clear from properties of the truncated Taylor series. $\square$

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 $\epsilon^{-1} = \mathcal{O}(1)$ and time $t_{\max} = \mathcal{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 \mapsto \mathrm{tr}[e^{t\mathcal{L}_\Lambda(O_Y)\rho}]$ 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) = \mathrm{tr}[\rho e^{t\mathcal{L}_\Lambda(O)}]$ for some suitably-picked initial state $\rho$ 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 \in [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), \dots, f(t_{\tilde{d}})$ points for $\tilde{d} \sim 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 $\epsilon_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 $\epsilon_s$ by sampling from the output of the device at time $t \mathcal{O}(\epsilon_s^{-2})$ 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, \hat{p}$ that are $\epsilon$ close in some interval $[a, b]$, how well can we infer the derivative of $p$ at 0 from that of $\hat{p}$? First, note that if the interval $[a, b]$ is very small, then two functions can differ by $\epsilon$ and their derivatives can still differ by $\epsilon/(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] \to \mathbb{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 \ge 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 $\alpha > 1/2$ of them satisfies for some $\sigma > 0$ that:

$$p(x_i) = y_i + w_i, \quad |w_i| \le \sigma. \qquad (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

$$\frac{1}{\pi\sqrt{1-x^2}} \qquad (91)$$

on the interval $[-1, 1]$. We then have:

**Theorem 5.1 (Robust polynomial interpolation).** Let $p: [-1, 1] \to \mathbb{R}$ be a polynomial of degree $d$ and assume we are given $m$ samples $(x_i, y_i)$ such that a fraction $\alpha > 1/2$ of them satisfies Eq. (90) for some $\sigma > 0$. Moreover, suppose that the $x_i$ were sampled independently and at random from the Chebyshev measure. Then for any $\delta > 0$

$$m = \mathcal{O}\left(d \log\left(\frac{d}{\delta}\right)\right) \qquad (92)$$

samples suffice to with probability of success at least $1 - \delta$ recover a polynomial $\hat{p}$ that satisfies:

$$\max_{x \in [-1, 1]} |p(x) - \hat{p}(x)| \le 3\sigma. \qquad (93)$$

Moreover, $\hat{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 $\epsilon < 1/4$ in their statement. $\square$

We note that the same result holds for random points picked from the uniform measure with $m = \mathcal{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 $\sigma$, then we also recover the whole polynomial up to some error proportional to $\sigma$. 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 \in \mathbb{N}$ define the constant $C_M(d, k)$ to be given by

$$C_M(d, k) = \frac{d^2 (d^2 - 1^2) (d^2 - 2^2) \cdots (d^2 - (k-1)^2)}{1 \cdot 3 \cdot 5 \cdots (2k-1)} \quad (94)$$

Then for any polynomial $p$ of degree $d$ we have that:

$$\max_{x \in [-1,1]} |p^{(k)}(x)| \le C_M(d, k) \max_{x \in [-1,1]} |p(x)|. \quad (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 $\tilde{p}(x) = p(\frac{b-a}{2}x + \frac{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 $\tilde{p}$, it follows from a straightforward application of the chain rule that

$$\max_{x \in [a,b]} |p^{(k)}(x)| \le \left| \frac{2}{(b-a)} \right|^k C_M(d,k) \max_{x \in [a,b]} |p(x)|. \quad (96)$$

From this we conclude that:

**Lemma 5.2** (Extrapolating the derivative at 0). Let $p: [0, b] \to \mathbb{R}$ be a polynomial of degree $d$ such that for some $\epsilon > 0$ and $0 < a < b$:

$$\max_{x \in [a,b]} |p(x)| \le \epsilon. \quad (97)$$

Then

$$\max_{x \in [0,a]} |p(x)| \le \epsilon \left( \sum_{k=0}^{d} \left| \frac{2}{(b-a)} \right|^k a^k \frac{C_M(d,k)}{k!} \right) \quad (98)$$

and

$$|p'(0)| \le \epsilon \left( \sum_{k=1}^{d} \left| \frac{2}{b-a} \right|^k a^{k-1} \frac{C_M(d,k)}{(k-1)!} \right). \quad (99)$$

*Proof.* It follows from Eq. (97) and Markov brothers' inequality that

$$|p^{(k)}(a)| \le \left| \frac{2}{b-a} \right|^k C_M(d, k) \epsilon. \quad (100)$$

By a Taylor expansion we know that for $x \in [0, a]$:

$$p(x) = \sum_{k=0}^{d} p^{(k)}(a) \frac{(x-a)^k}{k!}. \quad (101)$$

The claim in Eq. (98) then follows by combining this expansion with Eq. (100) and a triangle inequality. Similarly we have

$$p'(0) = \sum_{k=1}^{d} p^{(k)}(a) (-1)^k \frac{a^{k-1}}{(k-1)!}, \quad (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 - \hat{p}$, as in Eq. (93) to control the error we make by estimating the derivative at 0 by evaluating $\hat{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) = \left( \sum_{k=1}^{d} \left| \frac{2}{b-a} \right|^k a^{k-1} \frac{C_M(d, k)}{(k-1)!} \right). \qquad (103)$$

Then for $\delta > 0$ sampling

$$m = \mathcal{O}(d \log(d\delta^{-1})) \qquad (104)$$

i.i.d. points $(x_i, y_i)$ from the Chebyshev measure on $[a,b]$ satisfying

$$p(x_i) = y_i + w_i, \quad |w_i| \le \sigma$$

for at least a fraction $\alpha > \frac{1}{2}$ of the points is sufficient to obtain a polynomial $\hat{p}$ satisfying

$$|p'(0) - (\hat{p})'(0)| \le 3\sigma E(a, b, d). \qquad (105)$$

with probability of success at least $1 - \delta$.

*Proof.* From Thm. 5.1 we know that this number of samples suffices to obtain the polynomial $\hat{p}$ satisfying

$$\max_{x \in [a,b]} |p(x) - \hat{p}(x)| \le 3\sigma. \qquad (106)$$

Applying Lemma 5.2 to $p - \hat{p}$ yields the claim.

<figure><img src="image_7.png" /></figure>

This then yields a simple condition on how small $\sigma$ has to be in the regime of interest to us:

**Corollary 5.1.** In the same setting as in Prop. 5.1 for some $\epsilon > 0$ let $a \le d^{-2}$, $b = 2 + a$ and

$$\sigma = \epsilon d^{-2}. \qquad (107)$$

Then

$$|p'(0) - \hat{p}'(0)| = \mathcal{O}(\epsilon). \qquad (108)$$

*Proof.* It is easy to see that we have:

$$C_M(d, k) \le \frac{d^{2k}}{k!!}, \qquad (109)$$

where $k!! = 1 \times 3 \times 5 \times \dots \times (2k - 1)$ is the double factorial.

Thus, we see from this and Eq. (99) that by our choice

$$a = \frac{1}{d^2}, \quad b = 2 + a \qquad (110)$$

we have from Eq. (105) that the estimated polynomial $\hat{p}(0)$ satisfies

$$|\hat{p}'(0) - p'(0)| \le 3\sigma \left( \sum_{k=1}^{d} \left| \frac{2}{(a-b)} \right|^k a^{-1} \frac{(ad^2)^k}{k!!(k-1)!} \right) \le 3\sigma \left( \sum_{k=1}^{d} a^{-1} \frac{1}{k!!(k-1)!} \right), \qquad (111)$$

where we used the fact that $a \le d^{-2}$.

Thus, as

$$\left( \sum_{k=1}^{d} \frac{1}{k!!(k-1)!} \right) \le e, \qquad (112)$$

we conclude that with this choice of parameters we have

$$|\hat{p}'(0) - p'(0)| \le 3e\sigma a^{-1}, \qquad (113)$$

which gives the claim by our choice of $\sigma$.

<figure><img src="image_8.png" /></figure>

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_\alpha$ of $\mathcal{L}_\Lambda$ we will be interested in estimating the sample complexity of obtaining an estimate $\hat{a}_\alpha$ satisfying

$$ |a_{\alpha} - \hat{a}_{\alpha}| \le \epsilon \qquad (114) $$

with high probability for some given error $\epsilon$. 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 $\mathcal{L}_\Lambda$ be a locally bounded Lindbladian on a $D$-dimensional regular lattice with growth constant $g = \mathcal{O}(1)$. Let $O$ be an observable of constant support and $\rho$ and arbitrary quantum state. Let $\epsilon > 0$ be given. Assume that $\mathcal{L}_\Lambda$ satisfies Eq. (81) for some function $h$. Then picking $t_0$ as

$$ t_0 = \mathcal{O} \left[ \left( h^{-1} \left( \frac{\epsilon}{2(e^{2.5v} - 1)} \right)^D \log(\epsilon^{-1}) \right)^{-2} \right] \qquad (115) $$

and $t_{max} = 2 + t_0$ and measuring the expectation value $f(t) = \text{tr}[e^{t\mathcal{L}_\Lambda}(O)\rho]$ for

$$ m = \tilde{O} \left[ \left( h^{-1} \left( \frac{\epsilon}{2(e^{2.5v} - 1)} \right)^D \log(\epsilon^{-1}) \right) \right] \qquad (116) $$

random times $t_i \in [t_0, t_{max}]$ up to precision $\mathcal{O}(\epsilon)$ is sufficient to obtain an estimate $(\hat{f})'(0)$ satisfying

$$ \left| (\hat{f})'(0) - \text{tr}[\mathcal{L}_\Lambda(O)\rho] \right| \le 3e\epsilon t_0^{-1} = \mathcal{O} \left[ \epsilon \left( h^{-1} \left( \frac{\epsilon}{2(e^{2.5v} - 1)} \right)^D \log(\epsilon^{-1}) \right)^2 \right]. \qquad (117) $$

In particular, this estimate can be obtained from

$$ \tilde{O}(\epsilon^{-2} \log(\delta^{-1})) \qquad (118) $$

samples from the time evolved state $\rho$ with probability of success at least $1 - \delta$.

*Proof.* First, note that by Lemma. 4.1, if we pick $t_{max}$ as described above, then a polynomial $p$ of degree

$$ d = \mathcal{O} \left[ \left( h^{-1} \left( \frac{\epsilon}{2(e^{2.5v} - 1)} \right)^D \log(\epsilon^{-1}) \right) \right] \qquad (119) $$

is sufficient to approximate the expectation value in the interval $[0, t_{max}]$ up to an error $\epsilon/2$, assuming that $\epsilon$ is small enough to ensure that $t_{max} \le 2.5$. We will estimate the value of the polynomial at each point up to an error $\sigma > 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 $\sigma = \mathcal{O}(\epsilon)$ to obtain an overall error of

$$ \mathcal{O} \left[ \epsilon \left( h^{-1} \left( \frac{\epsilon}{2(e^{2.5v} - 1)} \right)^D \log(\epsilon^{-1}) \right)^2 \right] \qquad (120) $$

As we imposed that the precision with which the polynomial approximates the expectation values is $\epsilon/2$, we can estimate the value of the polynomial for a given time up to an error $\mathcal{O}(\epsilon)$ from $\mathcal{O}(\epsilon^{-2})$ samples.

As we have to sample

$$ \mathcal{O}(d \log(d)) = \tilde{O} \left[ \left( h^{-1} \left( \frac{\epsilon}{2(e^{2.5v} - 1)} \right)^D \log(\epsilon^{-1}) \right) \right] \qquad (121) $$

points to perform the stable interpolation, we obtain the advertised sample complexity. □
<table><thead><tr><th></th><th>Sample Complexity</th><th>Number of points</th><th>Initial time</th></tr></thead><tbody><tr><td>Finite range</td><td>ε<sup>-2</sup></td><td>polylog(ε<sup>-1</sup>)</td><td>polylog(ε)</td></tr><tr><td>Exponentially</td><td>ε<sup>-2</sup></td><td>polylog(ε<sup>-1</sup>)</td><td>polylog(ε)</td></tr><tr><td>Algebraically (α ≥ 5D - 1)</td><td>ε<sup>-2</sup> α<sup>-3D</sup> α<sup>-5D</sup></td><td>ε<sup>-D</sup> α<sup>-5D</sup></td><td>ε<sup>-2D</sup> α<sup>-5D</sup></td></tr></tbody></table>

**Supplementary Table VI.** Scaling of different resources required to obtain a recovery up to additive error $\epsilon$ of a parameters of the evolution. We have only included the leading order term and $\alpha$ denotes the decay of the potential in space, whereas $D$ the dimension of the lattice.

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

$$h^{-1} \left( \frac{\epsilon}{2(e^{2.5\nu} - 1)} \right) = \text{poly log}(\epsilon^{-1}). \qquad (122)$$

In that case the sample complexity is of order $\tilde{O}(\epsilon^{-2})$. Thus, in this case we see that the inverse initial time $t_0^{-1}$ and the number of points we need to sample from is polylogarithmic in $\epsilon$. Furthermore, the sample complexity to obtain an error $\epsilon$ is also $\tilde{O}(\epsilon^{-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

$$\left( h^{-1} \left( \frac{\epsilon}{2(e^{2.5\nu} - 1)} \right)^D \log(\epsilon^{-1}) \right)^{-2} = o(\epsilon^{-1}). \qquad (123)$$

Indeed, in this case we have the property that it is possible to suitably re-scale the error $\epsilon$ to ensure that the total precision is at some desired precision $\tilde{\epsilon}$. For instance, let us assume that

$$h^{-1} \left( \frac{\epsilon}{2(e^{2.5\nu} - 1)} \right) = \mathcal{O}(\epsilon^{-r}), \qquad (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

$$\mathcal{O}(\epsilon^{1-2Dr} \log(\epsilon^{-1})^2). \qquad (125)$$

Ignoring the $\log(\epsilon^{-1})$ term, we see that by picking $\epsilon = \tilde{\epsilon}^{1-1/2Dr}$ we can ensure an error of order $\tilde{\epsilon}$ for the estimate. Thus, the growth of $h^{-1}$ has to be at most $r \le 1/2D$ and the sample complexity would also grow like $\tilde{O}(\epsilon^{-2-2Dr-r})$, as we would need to sample $\tilde{O}(d) = \tilde{O}(\epsilon^{-r})$ points up to precision $\tilde{\epsilon}^{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 $\epsilon^{-2}$ scaling for short range 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 $l_1$ regression and an iteration of $l_\infty$ 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. $l_1$ regression: before we define the $l_1$ regression, we need to define the Chebyshev partitions:

**Definition 6.1** (Chebyshev partitions). Let $m \in \mathbb{N}$ be given. The size $m$ Chebyshev partitions of $[-1, 1]$ is the set of intervals $I_j = [\cos \frac{\pi j}{m}, \cos \frac{\pi (j-1)}{m}]$ for $1 \le j \le m$.
We also define $\mathcal{P}_d$ to be the space of polynomials of degree at most $d$.

With these definitions at hand, we define the $\ell_1$ regression solution as follows:

**Definition 6.2.** Given a set of $n$ points $(x_i, y_i)$ and $m \in \mathbb{N}$, we define the result of the degree $d$ $\ell_1$ regression with $m$ Chebyshev partitions $\hat{p}$ to be the polynomial

$$ \operatorname{argmin}_{\hat{p} \in \mathcal{P}_d} \sum_{i=1}^{n} |I_j| \operatorname{mean}_{x_i \in I_j} |y_i - \hat{p}(x_i)|, \qquad (126) $$

where $\mathcal{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 $\ell_1$ regression problem with $n = \mathcal{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 - \hat{p}\|_{\ell_1} = \mathcal{O}(\sigma), \qquad (127) $$

where as usual $\sigma$ is the error in each estimate $y_i$ and

$$ \|p - \hat{p}\|_{\ell_1} = \int_{-1}^{1} |p(x) - \hat{p}(x)| dx. \qquad (128) $$

However, the results of the previous sections required us to obtain a good solution in the $\| \cdot \|_{\ell_\infty}$ distance, and in general

$$ \|p - \hat{p}\|_{\ell_1} = \mathcal{O}(d^2 \|p - \hat{p}\|_{\infty}). \qquad (129) $$

Although, as commented in the last section, we are interested in the regime of polynomial of relatively small degree, by adding a $\ell_\infty$ regression iteration on top of the $\ell_1$ regression, it is possible to get rid of this $d^2$ prefactor.

b. $\ell_\infty$ regression: besides getting rid of the unwanted $d^2$ factor on the promise for the error of the $\ell_1$ regression, adding a $\ell_\infty$ regression step also has the favourable feature of making the whole procedure more robust to outliers in the data.

**Definition 6.3** ($\ell_\infty$ regression). Given a set of $n$ points $(x_i, y_i)$ and $m \in \mathbb{N}$ given. For the $m$ Chebyshev partitions $I_j$, choose $\tilde{x}_j \in I_j$ arbitrarily and let

$$ \tilde{y}_j = \operatorname{median}_{x_i \in I_j} y_i. \qquad (130) $$

We define the result of the degree $d$ $\ell_\infty$ regression with $m$ Chebyshev partitions $\hat{p}$ to be the polynomial $\hat{p} \in \mathcal{P}_d$

$$ \operatorname{argmin}_{\hat{p} \in \mathcal{P}_d} \max_{j \in [m]} |\hat{p}(\tilde{x}_j) - \tilde{y}_j|. \qquad (131) $$

Note that the problem in Eq. (131) also corresponds to a linear program and, thus, can be solved efficiently. The output of the $\ell_\infty$ regression algorithm is guaranteed to satisfy

$$ \|\hat{p} - p\|_{\infty} \le 2.5\sigma + \frac{1}{2}\|p\|_{\infty} \qquad (132) $$

as long as $m = \mathcal{O}(d)$ and we pick $n = \mathcal{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 $\infty$-norm up to the unwanted $\|p\|_{\infty}$ term. This can be solved by iterating the $\ell_\infty$ regression step.

c. Iterating the $\ell_\infty$ step: the last step to obtain the desired robust polynomial interpolation is to iteratively apply the $\ell_\infty$ iteration step to the residual. More precisely, we first perform the $\ell_1$-regression on our data, obtaining a polynomial $\hat{p}_0$. We can then define the new data points

$$ (x_i, \tilde{y}_i^0 = y_i - \hat{p}_0(x_i)) \qquad (133) $$

and run the $\ell_\infty$ interpolation on this residual error, obtaining a polynomial $\hat{p}_1$. From Eq. (129) and our promise on the output of the $\ell_\infty$ interpolation, we know that the result of the interpolation will satisfy

$$ \|p - \hat{p}_1\|_{\infty} \le 2.5\sigma + \frac{1}{2}\mathcal{O}(d^2\sigma). \qquad (134) $$
But then we can iterate this procedure by just running the $\ell_\infty$ regression on

$$ (x_i, \tilde{y}_i^1 = y_i - \hat{p}_1(x_i)). \qquad (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(\epsilon)$, 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.

## 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–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 $\text{dist}(X, Y)$ between the regions $X, Y$, i.e. a bound of the form

$$ \| [A_X, O_Y(t)] \| \leq Ch(\text{dist}(X, Y))(e^{vt} - 1), \qquad (136) $$

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, \cdot]$ for an arbitrary bounded super-operator $K_X: M_{2^n} \to M_{2^n}$ supported on $X$ leading to a Lieb-Robinson-bound of the form

$$ \| K_X(O_Y(t)) \| \leq Ch(\text{dist}(X, Y))(e^{vt} - 1), \qquad (137) $$

with $C$ depending on $\|K_X\|_{\infty \to \infty,cb}$. However, if $K_X$ is of the form $K_X = [A_X, \cdot]$, we have $\|K_X\|_{\infty \to \infty,cb} \leq 2 \|A_X\|_{\infty}$, which allows us to recover the commutator [10, 19]. In the following, we consider a regular lattice $\Lambda$ and assume that the dynamics is generated by a Lindbladian that decomposes according to

$$ \mathcal{L} = \sum_{X \subset \Lambda} \mathcal{L}_X. \qquad (138) $$

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

$$ J = \sup_{X \subset \Lambda} \| \mathcal{L}_X \|_{1 \to 1,cb} \qquad (139) $$

as well as the decay behaviour of the interactions

$$ \mu(r) = \sup_{X \subset \Lambda: \text{diam}(X)=r} \frac{\|\mathcal{L}_X\|_{1 \to 1,cb}}{J} \qquad (140) $$
in terms of the stabilized 1-to-1-norm $||T||_{1→1,cb} = \sup_n ||T ⊗ \text{id}_n||_{1→1}$. We can then characterize $\mathcal{L}$ as finite range if $\mu(r) = 0$ for $r ≥ R > 0$, exponentially decaying if $\mu(r) ≤ e^{-\mu r}$ and algebraically decaying if $\mu(r) ≤ (1 + r)^{-\alpha}$ for $\alpha > 0$ and state the following Lieb-Robinson-bound for Lindbladians:

**Theorem 7.1** (dissipative LR-bound [19]). Let $\mathcal{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(\text{dist}(X, Y))(e^{vt} - 1), \quad (141) $$

with $h(r) = e^{-\nu r}$ for $\mathcal{L}$ exponentially decaying or finite range and $h(r) = (1 + r)^\nu$ if $\mathcal{L}$ is algebraically decaying with $\alpha > 2D + 1$ with $\nu < \alpha - (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^{t\mathcal{L}_\Lambda} - e^{t\mathcal{L}_{\Lambda_r(Y)}})(O_Y)|| ≤ c_1h(\text{diam}(\Lambda_r(Y)))(e^{vt} - 1), \quad (142) $$

which reflects directly that the dynamics of the system can already be described by a generator $\mathcal{L}_{\Lambda_r(Y)}$ restricted to a region $\Lambda_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 $\mathcal{L}_\Lambda$ as compared to a restriction $\mathcal{L}_{\Lambda_r}$ to the subset $\Lambda_r ⊂ Λ$ according to

$$ (e^{t\mathcal{L}_\Lambda} - e^{t\mathcal{L}_{\Lambda_r}})O_Y = - \int_0^t ds \partial_s (e^{s\mathcal{L}_{\Lambda_r}} e^{(t-s)\mathcal{L}_\Lambda}) O_Y = \int_0^t ds e^{s\mathcal{L}_{\Lambda_r}} (\mathcal{L}_\Lambda - \mathcal{L}_{\Lambda_r}) e^{(t-s)\mathcal{L}_\Lambda} (O_Y) \quad (143) $$

Taking norms on both sides, we therefore obtain an upper bound of the form

$$ ||(e^{t\mathcal{L}_\Lambda} - e^{t\mathcal{L}_{\Lambda_r}})O_Y|| ≤ \sum_{X ⊂ Λ_r} \int_0^t ds ||\mathcal{L}_X e^{(t-s)\mathcal{L}_\Lambda}(O_Y)||. \quad (144) $$

We notice, that the term inside the integral is exactly of the form of the left-hand side of (137) with $K_X = \mathcal{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 $\mathcal{L}$ be a Lindbladian of the form (138), $O_Y$ an observable supported on $Y ⊂ Λ$ and $r > 0$ hold. Then for $\mathcal{L}_{\Lambda_r(y)} = \sum_{X ⊂ Λ_r(y)} Λ_X$ we have

$$ ||(e^{t\mathcal{L}} - e^{t\mathcal{L}_{\Lambda_r(y)}}) O_Y|| ≤ ||O_Y|| |Y| J \frac{e^{vt} - 1 - vt}{v} h(r) \quad (145) $$

with $h(r)$ exponentially decaying in $r$ for $\mathcal{L}$ finite range or exponentially decaying and $h(r)$ decaying as $(1 + r)^{-\beta}$ if $\mathcal{L}$ is algebraically decaying with $\alpha > 2D + 1$ and $\beta = \alpha - 3D$ for $\alpha ≥ 5D - 1$ and $\beta = \frac{1}{2}(\alpha - D - 1)$ if $\alpha ≤ 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. Explicit 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\Lambda} \sum_{n \ge 0} \frac{|b_x(n)|}{|b_x(n-1)|} \sum_{r \ge n} \mu(r) |b_y(r)| h(r) \le v, \quad (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 = \sup_{X \in E} \|\mathcal{L}_X\| \le \infty$

2. finite neighbourhood: $Z = \max_{X \in E} |\{Y \in E : X \cap Y \ne \emptyset\}| \le \infty$

3. bounded growth of neighbours: $\max_{X \in E} |\{Y : \text{dist}(X, Y) = n\}| \le Mn^\eta$

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 by

$$\frac{2M}{Z} \|O_Y\| r^{\eta} e^{(egZ) \cdot t-r} = \frac{2M}{Z} e^{v \cdot t} h(r) \quad (147)$$

with $h(r) = r^{\eta}e^{-r}$ and the LR-velocity $v = egZ$. Hence, following the reasoning in section 4, we choose $r$ as $h^{-1}(\frac{Z\varepsilon}{4M\|O_Y\|\exp(v \cdot t_{max})})$, which for $\|O_Y\| \le 1$ and for $r$ large enough such that $(1 - \eta \frac{\log r}{r}) > \frac{2}{3}$ means that $r$ should be chosen as $\frac{3}{2} \log(\frac{4Me^{(egZ) \cdot t_{max}}}{Z\varepsilon})$. This in turn results for the truncated Taylor-expansion for the corresponding generator $\mathcal{L}_{\Lambda_y(r)}$ restricted to the region $\Lambda_y(r)$ in a polynomial degree of $d = 2et_{max}g|\Lambda_y(r)| - 1$, which in the case of a $D$-dimensional lattice results in

$$d_{exp} = 2et_{max}g \left( \frac{3}{2} \log \left( \frac{4Me^{(egZ) \cdot t_{max}}}{Z\varepsilon} \right) + \text{diam}(\text{supp}(O_Y)) \right)^D \log(\varepsilon^{-1}) - 1. \quad (148)$$

Following the reasoning in section 6, this implies a required minimal starting time $t_0 = d_{exp}^2$ and maximal evolution time $t_{max} = 2 + t_0$ as well as sampling $m = O(d \log(d\delta^{-1})$ evolution times within this interval according to the Chebychev distribution in order to succeed with the polynomial interpolation.

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 \sim \sqrt{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.

## 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 $\Phi$ and two collections of Pauli strings $P_a^1, \dots, P_a^{K_1}$ and $P_b^1, \dots, P_b^{K_2}$ that have combined weight (i.e. number of Paulis that are nonidentity) $\omega_a + \omega_b$ at most $\omega$ it allows us to estimate all the overlaps of the form

$$2^{-n} \text{tr} [P_a^j \Phi(P_b^k)] \quad (149)$$

up to an error $\epsilon$ with probability at least $1 - \delta$ from a number of samples that grows like $O(3^w \log(K_1 K_2 \delta^{-1}) \epsilon^{-2})$. 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:
1. Draw two Pauli strings $B_i, S_i \in \{\sigma_x, \sigma_y, \sigma_z\}^n$ uniformly at random and also a sign vector $E_i = \{-1, +1\}^n$ uniformly at random.

2. Prepare the quantum state $\rho_i = \otimes_{j=1}^n \phi_j(S_i, E_i)$, where $\phi_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 $\rho_i$ by $\Phi$.

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)$.

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 = \sigma_x \otimes \sigma_y \otimes I$ and $B_i = \{\sigma_x, \sigma_y, \sigma_z\}$ overlap. However, $P_a = \sigma_x \otimes \sigma_y \otimes I$ and $B_i = \{\sigma_x, \sigma_z, \sigma_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 $\rho_i$ is a positive eigenstate of $P_b$. We say we overlap negatively if it is a negative eigenstate. We will also define $\omega(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) = \begin{cases} 0 & \text{if } B_i \text{ does not overlap with } P_a \text{ or } S_i \text{ does not overlap with } P_b, \\ \frac{3^{\omega(P_a)+\omega(P_b)}}{2} & \text{if } (B_i, S_i) \text{ overlap with } (P_a, P_b), \text{ and both do so positively,} \\ \frac{3^{\omega(P_a)+\omega(P_b)}}{2} & \text{if } (B_i, S_i) \text{ overlap with } (P_a, P_b), \text{ and both do so negatively,} \\ \frac{-3^{\omega(P_a)+\omega(P_b)}}{2} & \text{if } (B_i, S_i) \text{ overlap with } (P_a, P_b), \text{ one positively, the other negatively.} \end{cases} \quad (150)$$

We will now show that:

$$\mathbb{E}(X_{a,b}) = 2^{-n} \text{tr} (P_a \Phi(P_b)). \quad (151)$$

and

$$\mathbb{E}(X_{a,b}^2) \le 3^{w(P_a)+w(P_b)} \quad (152)$$

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 $\hat{\mu}_i$, with $1 \le i \le K$. We then set our estimator of the mean to be:

$$\hat{\mu}_{\text{MoM}} = \text{median}(\hat{\mu}_1, \dots, \hat{\mu}_K). \quad (153)$$

The main property of this estimator is that we have that if the variance of $X_{a,b}$ is $\sigma^2$, then:

$$\mathbb{P}(|\hat{\mu}_{\text{MoM}} - \mathbb{E}(X_{a,b})| \ge \epsilon) \le e^{-2K\left(\frac{1}{2} - \frac{\sigma^2}{B\epsilon^2}\right)}. \quad (154)$$

In particular, if we pick $B = 4\sigma^2\epsilon^{-2}$ and $K = 2\log(\delta^{-1})$, then:

$$\mathbb{P}(|\hat{\mu}_{\text{MoM}} - \mathbb{E}(X_{a,b})| \ge \epsilon) \le \delta. \quad (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):
**Corollary 8.1.** Let {$P_a^l$}$_{l=1}^{K_1}$ and {$P_b^j$}$_{j=1}^{K_2}$ be two collections of Pauli matrices on n qubits such that for any l, j the following condition

$$ \omega(P_a^l) + \omega(P_b^j) \le \omega, \qquad (156) $$

holds. Then

$$ O(3^{\omega \log (K_1 K_2 \delta^{-1}) \epsilon^{-2}}) \qquad (157) $$

runs of the protocol above suffice to obtain an estimate $e_{a,b}^{l,j}$ satisfying

$$ \left| |2^{-n} \operatorname{tr} [P_b^l \Phi(P_a^j)] - e_{a,b}^{l,j}| \right| \le \epsilon \qquad (158) $$

for all pairs l, j with probability at least $1 - \delta$.

*Proof.* Let $X_{a,b}^{l,j}$ be the random variable we have defined above for a pair of Paulis $P_a^l$ and $P_b^j$. By the bound in (152) we have that $\mathbb{E}((X_{a,b}^{l,j})^2) \le 3^\omega$, holds. Thus, for the median of means estimator with $B = 3^\omega \epsilon^{-2}$ samples per group and $K = 2 \log(\delta^{-1} K_1 K_2)$ we have an estimate satisfying (158) with probability of failure at most $\delta/K_1 K_2$.

By a union bound, the median of means estimator of all $K_1 K_2$ combinations of $P_a^l$ and $P_b^j$ satisfy (158) with probability of failure at most $\delta$. The total number of samples required for this is

$$ S = KB = O(3^{\omega \log (K_1 K_2 \delta^{-1}) \epsilon^{-2}}), \qquad (159) $$

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

$$ \mathbb{E}(X_{a,b}) = \mathbb{E}(X_{a,b}|B \text{ and } S \text{ overlap})\mathbb{P}(B \text{ and } S \text{ overlap}), \qquad (160) $$

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

$$ \mathbb{P}(B \text{ and } S \text{ overlap}) = 3^{-\omega(P_a)} 3^{-\omega(P_b)}. \qquad (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

$$ \mathbb{E}(X_{a,b}|B \text{ and } S \text{ overlap}). \qquad (162) $$

Here we distinguish four cases:

1. $S$ overlaps positively with $P_b$ and $B$ overlaps positively with $P_a$.

2. $S$ overlaps negatively with $P_b$ and $B$ overlaps negatively with $P_a$,

3. $S$ overlaps positively with $P_b$ and $B$ overlaps negatively with $P_a$,

4. $S$ overlaps negatively with $P_b$ and $B$ overlaps positively with $P_a$.

We can then break down the expectation of (160) down into the cases:

$$ \mathbb{E}(X_{a,b}|B \text{ and } S \text{ overlap})\mathbb{P}(B \text{ and } S \text{ overlap}) = \sum_{i=1}^{4} \mathbb{E}(X_{a,b}|B \text{ and } S \text{ overlap}, c=i)\mathbb{P}(c=i), \qquad (163) $$

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

$$ \mathbb{E}(X_{a,b}|B \text{ and } S \text{ overlap}, c = 1, 2) = 3^{\omega(P_a) + \omega(P_b)/2}. \qquad (164) $$

If $c = 3, 4$:

$$ \mathbb{E}(X_{a,b}|B \text{ and } S \text{ overlap}, c = 3, 4) = -3^{\omega(P_a) + \omega(P_b)/2}. \qquad (165) $$

Note that the $3^{\omega(P_a)+\omega(P_b)}$ counteracts the $\mathbb{P}(B \text{ and } S \text{ overlap})$ term, up to an additional 1/2 factor.
Thus, all that is left is to estimate

$$ \mathbb{P}(c = i), \quad i = 1, 2, 3, 4. \qquad (166) $$

Let us estimate $\mathbb{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 $Q_b^+$ and $Q_b^-$. Clearly, $P_a = (Q_a^+ - Q_a^-) \otimes I^{\otimes n-\omega(P_a)}$, holds. 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

$$ \sigma = \frac{Q_b^+}{2\omega(P_b)-1} \otimes \left(\frac{\mathbb{I}}{2}\right)^{\otimes n-\omega(P_b)}. \qquad (167) $$

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 $\text{tr}(Q_a^+ \otimes \mathbb{I}^{n-\omega(P_a)} \Phi(\sigma))$.

We conclude that

$$ \mathbb{P}(c = 1) = \frac{1}{2} \text{tr}(Q_a^+ \otimes \mathbb{I}^{n-\omega(P_a)} \Phi(\sigma)) = \frac{1}{2^{n-1}} \text{tr}(Q_a^+ \otimes \mathbb{I}^{n-\omega(P_a)} \Phi(Q_b^+ \otimes \mathbb{I}^{n-\omega(P_b)})). \qquad (168) $$

Similarly

$$ \mathbb{P}(c = 2) = \frac{1}{2^{n-1}} \text{tr}(Q_a^- \otimes 1^{n-\omega(P_a)} \Phi(Q_b^- \otimes 1^{n-\omega(P_b)})), \qquad (169) $$

$$ \mathbb{P}(c = 3) = \frac{1}{2^{n-1}} \text{tr}(Q_a^+ \otimes 1^{n-\omega(P_a)} \Phi(Q_b^- \otimes 1^{n-\omega(P_b)})), \qquad (170) $$

$$ \mathbb{P}(c = 4) = \frac{1}{2^{n-1}} \text{tr}(Q_a^- \otimes 1^{n-\omega(P_a)} \Phi(Q_b^+ \otimes 1^{n-\omega(P_b)})), \qquad (171) $$

As

$$ 2^{-n} \text{tr}(P_a \Phi(P_b)) = 2^{-n}(\text{tr}(Q_a^+ \Phi(Q_b^+)) + \text{tr}(Q_a^- \Phi(Q_b^-)) - \text{tr}(Q_a^+ \Phi(Q_b^-)) - \text{tr}(Q_a^- \Phi(Q_b^+))), \qquad (172) $$

we conclude that

$$ \mathbb{E}(X_{a,b}) = 2^{-n} \text{tr}(P_a \Phi(P_b)). \qquad (173) $$

Computing the second moment of $X_{a,b}$ turns out to be quite simple. Note that the only nonzero value the random variable $X_{a,b}^2$ takes is $3^{2\omega(P_a)+2\omega(P_b)}/4$ with probability $3^{-\omega(P_a)-\omega(P_b)}$. Thus, we clearly have that:

$$ \mathbb{E}(X_{a,b}^2) = 3^{\omega(P_a)+\omega(P_b)}/4, \qquad (174) $$

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 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:
1. perturbations of the locality of the interaction.

2. SPAM errors.

3. 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 $\mathcal{L} = \mathcal{L}' + \Delta$, where $\mathcal{L}$ is a Lindbladian that satisfies our locality assumptions and $\Delta$ 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:

**Lemma 9.1.** Let $\mathcal{L}$ be a Lindbladian on a $D$-dimensional regular lattice and $\epsilon > 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 $\{\Phi_t\}_{t_0 \le t \le t_{\max}}$ s.t. for all $t \in [t_0, t_{\max}]$ we have for 2-body Paulis P:

$$ \| (\Phi_t^* - e^{t\mathcal{L}})(P) \| \le \mathcal{O}(\epsilon/\log(\epsilon^{-1})). \qquad (175) $$

for

$$ t_0^{-1} = \mathcal{O}[\text{polylog}(\epsilon^{-1})] \qquad (176) $$

and $t_{\max} = 2 + t_0$. Assume further that we can prepare Pauli eigenstates and perform Pauli measurements up to an error $\mathcal{O}(\epsilon/\log(\epsilon^{-1}))$. Then, measuring such Pauli measurements on Pauli eigenstates under the evolution of $\Phi_t^*$ for

$$ m = \mathcal{O}[\text{polylog}(\epsilon^{-1})] \qquad (177) $$

random times up to precision $\mathcal{O}(\epsilon/\text{polylog}(\epsilon^{-1}))$, is sufficient to obtain with probability of success at least 2/3 an estimate $\hat{a}_{\alpha}$ of $a_{\alpha}$ satisfying

$$ |\hat{a}_{\alpha} - a_{\alpha}| = \epsilon. \qquad (178) $$

*Proof.* By Eq. (175), measuring the expectation value of $P$ under the time evolution will yield values that are within $\mathcal{O}(\epsilon/\log(\epsilon^{-1}))$ of the expectation value for the evolution under $\mathcal{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. $\square$

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

$$ \mathcal{L} = \mathcal{L}' + \Delta, \qquad (179) $$

where $\mathcal{L}'$ is a time evolution that satisfies a LR bound and $\Delta$ is an arbitrary superoperator that is supposed to model the non-local parts of the evolution. Furthermore, we will define the following expansion quantity $\kappa(\Delta, k)$ that measures the impact of the non-local part on observables supported on at most $k$ sites:

$$ \kappa(\Delta, k) = \sup_{O:\text{supp}(O) \le k, \|O\| \le 1} \|\Delta(O)\|. \qquad (180) $$
To see why this quantity measures the impact of global perturbations on local observables, let us consider the case where we have

$$ \Delta(X) = \tau \sum_{i,j=1, i \neq j}^{n} i[H_{i,j}, X], \qquad (181) $$

for some arbitrary Hamiltonian evolution terms $||H_{i,j}|| \le 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 $\sim n^2$ in general. In contrast, $\kappa(\Delta, k) \sim \tau 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 $\sim nk$ terms acting on that site.

We then have:

**Theorem 9.1.** Let $\mathcal{L}$ be a Lindbladian and suppose it is of the form

$$ \mathcal{L} = \mathcal{L}' + \Delta, \qquad (182) $$

where $\mathcal{L}'$ satisfies a LR bound as in Eq. (142). Then we have for an observables $O$ s.t. $||O|| \le 1$ and

$$ ||e^{t\mathcal{L}} - e^{t\mathcal{L}'}(O)|| \le \inf_{\delta>0} t\delta ||\Delta||_{+\infty \to +\infty} + \int_0^t ds\kappa(\Delta, |B(\text{supp}(O), \log(\delta^{-1}) + vs)|). \qquad (183) $$

*Proof.* We have that:

$$ (e^{t\mathcal{L}} - e^{t\mathcal{L}'})(O) = \int_0^t dse^{(t-s)\mathcal{L}}\Delta e^{s\mathcal{L}'}(O). \qquad (184) $$

Now for any $\delta > 0$ we can approximate $e^{s\mathcal{L}'}(O)$ up to $\delta$ in operator norm by an observable supported on a ball of radius $r_{s,\delta} = vs + \log(\delta^{-1})$ around the original support of $O$, i.e.:

$$ e^{s\mathcal{L}'}(O) = O_{r_s} + \delta\tilde{O}_s \qquad (185) $$

where $||O_{r_s}||, ||\tilde{O}_s|| \le 1$ and $\tilde{O}_s$ does not necessarily have any further locality. Now:

$$ ||(e^{t\mathcal{L}} - e^{t\mathcal{L}'})(O)|| \le \int_0^t ds ||e^{(t-s)\mathcal{L}}\Delta e^{s\mathcal{L}'}(O)|| \le \int_0^t ds (||e^{(t-s)\mathcal{L}}\Delta O_{r_{s,\delta}}|| + ||e^{(t-s)\mathcal{L}}\Delta(\delta\tilde{O}_s)||). \qquad (186) $$

We can estimate the first summand of the integral by

$$ \int_0^t ds ||e^{(t-s)\mathcal{L}}\Delta(O_{r_{s,\delta}})|| \le \int_0^t ds\kappa(\Delta, |B(\text{supp}(O), \log(\delta^{-1}) + vs)|) \qquad (187) $$

as, by definition, $||\Delta(O_{r_{s,\delta}})|| \le \kappa(\Delta, r(s, \delta))$ and $e^{(t-s)\mathcal{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:

$$ \int_0^t ds ||e^{(t-s)\mathcal{L}}\Delta(\delta\tilde{O}_s)|| \le \delta t ||\Delta||_{+\infty \to +\infty}. \qquad (188) $$

As $\delta > 0$ here was arbitrary, we can further optimize the bound over all $\delta > 0$, which gives the claim. $\square$

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, $||\Delta|| \le \tau \frac{n(n-1)}{2}$. Thus, for an observable initially supported on a constant number of sites, we have that:

$$ \kappa(\Delta, r(s, \delta)) = \tau n (vs + C \log(\delta^{-1}))^D, \qquad (189) $$