Learning the dynamics of open quantum systems from their steady states

TL;DR

This paper demonstrates that the steady states of open quantum systems governed by local Markovian dynamics contain enough information for full reconstruction of their dynamics using only local measurements, even in the presence of dissipation.

Contribution

It introduces a method for reconstructing the Lindbladian of open quantum systems from local steady-state measurements, extending previous work to include dissipative dynamics.

Findings

Reconstruction is feasible for generic local Markovian dynamics.

The method accurately recovers the Lindbladian with fewer measurements for systems with finite-range interactions.

External environment couplings can aid in reconstructing Hamiltonians with commuting terms.

Abstract

Recent works have shown that generic local Hamiltonians can be efficiently inferred from local measurements performed on their eigenstates or thermal states. Realistic quantum systems are often affected by dissipation and decoherence due to coupling to an external environment. This raises the question whether the steady states of such open quantum systems contain sufficient information allowing for full and efficient reconstruction of the system's dynamics. We find that such a reconstruction is possible for generic local Markovian dynamics. We propose a recovery method that uses only local measurements; for systems with finite-range interactions, the method recovers the Lindbladian acting on each spatial domain using only observables within that domain. We numerically study the accuracy of the reconstruction as a function of the number of measurements, type of open-system dynamics and system size. Interestingly, we show that couplings to external environments can in fact facilitate the reconstruction of Hamiltonians composed of commuting terms.

Figures

Equations

$$\dot{\rho }$$

$$=-\text{operatorfont}i\sum _j\left[H_j,\rho \right]+\frac{1}{2}\sum _j\left(\left[L_j\rho ,L_j^{†}\right]+\left[L_j,\rho L_j^{†}\right]\right),$$

$$\dot{\rho }=-\text{operatorfont}i\sum _i{c}_i\left[h_i,\rho \right]+\sum _{r,s}\frac{c_{rs}}{2}\left(\left[l_r\rho ,l_s^{†}\right]+\left[l_r,\rho l_s^{†}\right]\right),$$

$$\text{Tr}\left(A\dot{{\rho }_s}\right)=0.$$

$$-\sum _i{c}_i⟨\text{operatorfont}i\left[A,h_i\right]⟩+\sum _{r,s}\frac{c_{rs}}{2}⟨\left[l_r,A\right]l_s^{†}+l_r\left[A,l_s^{†}\right]⟩=0,$$

$${\vec{k}}^T\vec{c}=0$$

$$K\vec{c}=0,$$

$$\Delta =\text{norm}{{\hat{c}}_{recovered}-{\hat{c}}_{true}}_2.$$

$${\Delta }^{est}=ϵ\sqrt{\sum _{m>0}{\lambda }_m^{-1}},$$

$$H_j=\sum _{\alpha =1}^3{c}_{j,\alpha }{\sigma }_j^{\alpha }+\sum _{\alpha ,\beta =1}^3{c}_{j,\alpha ,\beta }{\sigma }_j^{\alpha }{\sigma }_{j+1}^{\beta },$$

$$L_j=\sum _{\alpha =1}^3{d}_{j,\alpha }{\sigma }_j^{\alpha }.$$

$$\sum _{r,s}\frac{c_{rs}}{2}⟨\left[l_r,A\right]l_s^{†}+l_r\left[A,l_s^{†}\right]⟩=⟨\text{operatorfont}i\left[A,H\right]⟩,$$

$$L_{j,L}={\alpha }_L{\sigma }_j^{-},L_{j,D}=(1-{\alpha }_L){\sigma }_j^z,$$

$$H_{cl}^x=\sum _{j=1}^{\Lambda }{b}_j{\sigma }_j^x+\sum _{j=1}^{\Lambda -1}{J}_j{\sigma }_j^x{\sigma }_{i+1}^x,$$

$$L_j=2{\sigma }_j^{-},$$

$$H_j=\sum _i{c}_i^{(j)}{h}_i,$$

$$\sum _j\left[H_j,\rho \right]=\sum _i{c}_i\left[h_i,\rho \right]$$

$$c_i=\sum _j{c}_i^{(j)}.$$

$$L_j=\sum _i{c}_r^{(j)}{l}_r,$$

$$\frac{1}{2}\sum _j\left(\left[L_j\rho ,L_j^{†}\right]+\left[L_j,\rho L_j^{†}\right]\right)=$$

$$=\sum _{r,s}\frac{c_{rs}}{2}\left(\left[l_r\rho ,l_s^{†}\right]+\left[l_r,\rho l_s^{†}\right]\right),$$

$$c_{rs}=\sum _j{c}_r^{(j)}(c_s^{(j)}{)}^{\ast }$$

$$K_{n,m}=$$

$$\{\begin{array}{ll}-⟨i\left[A_n,h_j\right]⟩& c_j\\ \frac{1}{2}⟨\left[l_r,A_n\right]l_s^{†}+l_r\left[A_n,l_s^{†}\right]⟩& c_{rr}\\ \frac{1}{2}\left(⟨\left[l_r,A_n\right]l_s^{†}+l_r\left[A_n,l_s^{†}\right]⟩+\left\{r\leftrightarrow s\right\}\right)& \Re c_{rs};r>s\\ \frac{i}{2}\left(⟨\left[l_r,A_n\right]l_s^{†}+l_r\left[A_n,l_s^{†}\right]⟩-\left\{r\leftrightarrow s\right\}\right)& \Im c_{r,s};r>s.\end{array}$$

$$-\sum _i{c}_i\text{Tr}\left(\text{operatorfont}i\left[A,h_i\right]\right)+\sum _{r,s}\frac{c_{rs}}{2}\text{Tr}\left(\left[l_r,A\right]l_s+l_r\left[A,l_s\right]\right)=0.$$

$$\text{Tr}\left(\left[l_r,A\right]l_s\right)+\text{Tr}\left(l_r\left[A,l_s\right]\right)$$

$$=\text{Tr}\left(\left[l_s,l_r\right]A\right)+\text{Tr}\left(A\left[l_s,l_r\right]\right)$$

$$=2\text{Tr}\left(\left[l_s,l_r\right]A\right),$$

$$\sum _{r,s}{c}_{rs}\text{Tr}\left(\left[l_s,l_r\right]A\right)=0.$$

$$K_l\overrightarrow{c_l}=\vec{b},$$

$$b_n=⟨\text{operatorfont}i\left[A_n,H\right]⟩.$$

$$L_j=\sum _{\alpha =1}^3{d}_{j,\alpha }{\sigma }_j^{\alpha }+d_{j,x,x}{\sigma }_j^x{\sigma }_{j+1}^x+d_{j,y,y}{\sigma }_j^y{\sigma }_{j+1}^y,$$

$$\rho (0)=\text{outerproduct}\uparrow \uparrow \cdots \uparrow \uparrow \uparrow \cdots \uparrow ,$$

Full text

I Introduction

The development of quantum simulators and computation devices has rapidly progressed over the last few years Preskill (2018). These developments span a multitude of physical platforms, including ultracold atoms Monroe (2002); Bloch et al. (2008); Gross and Bloch (2017); Bernien et al. (2017), trapped ions Cirac and Zoller (1995); Blatt and Roos (2012); Monroe and Kim (2013), photonic circuits Kok et al. (2007); Aspuru-Guzik and Walther (2012); Flamini et al. (2019); Takeda and Furusawa (2019), Josephson junction arrays You and Nori (2005); Houck et al. (2012); Devoret and Schoelkopf (2013); Wendin (2017); Neill et al. (2018) and more, reaching ever larger complexity. The growth in the complexity of these systems calls for efficient methods to characterize and verify their dynamics. The resources required by these methods, whether classical computations or quantum measurements, should scale polynomially with the number of degrees of freedom in the system.

An isolated quantum system can be characterized by learning its underlying Hamiltonian. This can be achieved by monitoring the dynamics that the Hamiltonian generates Da Silva et al. (2011); Burgarth et al. (2009); Di Franco et al. (2009); Shabani et al. (2011); Zhang and Sarovar (2014); Wang et al. (2015); De Clercq et al. (2016); Sone and Cappellaro (2017); Wang et al. (2018); Granade et al. (2012); Wiebe et al. (2014a, b, 2015); Wang et al. (2017), or by measuring local observables in one of its eigenstates or thermal states Rudinger and Joynt (2015); Huber and Gühne (2016); Kieferová and Wiebe (2017); Qi and Ranard (2019); Chertkov and Clark (2018); Greiter et al. (2018); Kappen (2018); Turkeshi et al. (2019); Bairey et al. (2019). However, realistic quantum systems are never fully isolated. This raises the need for methods to characterize the dynamics of open quantum systems which are coupled to external environments.

Previous works have recovered the dynamics of open quantum systems by tracking their time evolution Chuang and Nielsen (1997); Childs et al. (2001); Boulant et al. (2003); Mohseni and Lidar (2006); Howard et al. (2006); Bellomo et al. (2009, 2010); Da Silva et al. (2011); Akerman et al. (2012); Glickman et al. (2013); Zhang and Sarovar (2015); Ficheux et al. (2018). However, the possibility of recovering open system dynamics from their steady states has not been addressed. We focus on open quantum systems evolving under Markovian and local dynamics, for which the evolution can be described by the Lindblad master equation formalism Gorini (1976); Lindblad (1976):

[TABLE]

where each $H_{j}$, $L_{j}$ is a local operator. Throughout this paper, a local operator will be defined as acting on at most $k$ spatially contiguous degrees of freedom (e.g. spins). While the Hamiltonian terms $H_{j}$ are Hermitian, the $L_{j}$ operators, known as the ‘jump operators’, are generically not. A steady state $\rho_{s}$ of $\mathcal{L}$ is defined by $\dot{\rho}_{s}=\mathcal{L}(\rho_{s})=0$. Suppose that we prepare many copies of $\rho_{s}$ and measure expectation values of local observables in the state $\rho_{s}$. Can $\mathcal{L}$ be recovered using the data obtained from these measurements?

Parameter counting suggests this should be possible. The number of parameters describing a local Lindbladian scales polynomially with the system size, similarly to a Hamiltonian. On the other hand, a quantum state is described by exponentially many parameters. Thus, the steady state of a local Lindbladian may potentially contain sufficient information for inferring the dynamics that generated it.

However, steady states of Lindbladians differ from eigenstates and thermal states of local Hamiltonians. Every Hamiltonian commutes with the density matrix corresponding to each of its eigenstates $\outerproduct{\epsilon_{i}}{\epsilon_{i}}$. In contrast, generic Lindbladians have only a single steady state Evans (1977). Dissipation can cause this unique steady state to be highly mixed, possibly reducing its information content. As an extreme example, the steady state of any Lindbladian whose jump operators $L_{j}$ are all Hermitian is the fully mixed state $\rho\propto\mathbbm{1}$, from which there is no hope to recover the Lindbladian. Does this impose a fundamental difficulty to Lindbladian reconstruction? Or do the steady states of local many-body dissipative dynamics generically bear clear signatures of the preceding dynamics? Can these dynamics be extracted efficiently and accurately?

In this work, we study this question by providing an efficient algorithm for learning the dynamics of local Lindbladians from their steady states. Extending the methods of Ref. Bairey et al., 2019, our algorithm exploits strong constraints that locality imprints on the steady states of generic local Lindbladians. Using this algorithm, we (i) explore which types of Lindbladians can be accurately reconstructed from their steady states, (ii) study numerically and analytically the system-size scaling of the reconstruction accuracy, and (iii) show that coupling to a bath can in fact facilitate the reconstruction of certain classes of Hamiltonians, which pose a challenge for methods based on their eigenstates or Gibbs states.

II Algorithm

We begin by choosing a basis of local Hermitian operators for the unitary dynamics $\{{h_{i}\}}$, and a basis of local operator pairs for the dissipative dynamics $\{{\left(l_{r},l_{s}\right)\}}$. Expanding the dynamics in this operator basis (see Appendix A.1), Eq. (1) becomes

[TABLE]

with real coefficients $c_{j}$, and $c_{rs}$ forming a complex-valued positive semidefinite matrix. The locality of the Lindbladian restricts the pairs of non-zero elements of $c_{rs}$; for instance, if the jump operators $L_{j}$ are on-site, $c_{rs}$ vanishes whenever $l_{r},l_{s}$ act on different sites. Our goal is to infer the values of the non-zero coefficients $c_{j}$, $c_{rs}$.

To this end, we identify a set of local constraints that apply to any steady state $\rho_{s}$ of $\mathcal{L}$. Since $\rho_{s}$ is a steady state, the expectation value $\left\langle A\right\rangle\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\Tr\left({A\rho_{s}}\right)$ of any observable must be time-independent,

[TABLE]

Plugging in Eq. (2) and using the cyclic properties of the trace, $\Tr\left({ABC}\right)=\Tr\left({CAB}\right)$ and $\Tr\left({A\left[B,C\right]}\right)=\Tr\left({C\left[A,B\right]}\right)$, we obtain the linear constraint

[TABLE]

where the expectation values are taken with respect to the steady state $\rho_{s}$. For any operator $A$, Eq. (4) yields a linear equation for the parameters $c_{j}$ and $c_{rs}$. We will use a set of constraint operators $\{{A_{n}\}}$ to obtain a system of linear equations for the Lindbladian coefficients.

Importantly, assuming that local $A_{n}$ operators are chosen, the constraints derived from Eq. (4) are local in two ways. First, these constraints involve only local observables, which are easier to measure in most experimental settings. Second, if the $A_{n}$ operators act only within a given region, they commute with all the Lindblad terms that are supported outside that region. This allows to recover the Lindbladian of a region from measurements of that region alone.

We now introduce a convenient notation for representing the constraints derived from Eq. (4). We concatenate the Hamiltonian parameters $c_{j}$ and the dissipative parameters $c_{rs}$ into a single vector $\vec{c}$. In this notation, Eq. (4) takes the form

[TABLE]

for a corresponding vector of expectation values $\vec{k}$. Since $c_{rs}$ is Hermitian, its upper and lower parts are redundant; each pair of off-diagonal elements contributes only a single pair of real parameters, $\Re{c_{rs}}=\frac{1}{2}\left(c_{rs}+c_{sr}\right)$ and $\Im{c_{rs}}=\frac{1}{2{\mkern 1.5mu\operatorfont i\mkern 1.5mu}}\left(c_{rs}-c_{sr}\right)$. Thus, $\vec{c}$ is a real vector with four types of elements: Hamiltonian coefficients $c_{j}$, diagonal dissipative coefficients $c_{r,r}$, and the real and imaginary parts of the off-diagonal dissipative coefficients $c_{rs}$ for $r>s$.

Repeating this procedure for a set of constraints $\{{A_{n}\}}_{n=1}^{N}$, we obtain a homogeneous system of linear equations for the coefficients of the true Lindbladian,

[TABLE]

where $K$ is an ${N\times M}$ matrix of expectation values (see Appendix A.2), with $N$ the number of constraints and $M$ the number of unknown parameters. Each of its rows corresponds to a constraint operator $A_{n}$, and each column to a different Hamiltonian term or jump operator appearing in Eq. (2).

Assuming that we measured $K$ at a steady state of a local Lindbladian, the vector $\vec{c}$ corresponding to that Lindbladian must lie in the kernel of $K$. If the steady state is shared by a family of Lindbladians, the kernel will be spanned by the whole family (see Appendix A.3 for the example of the fully mixed state). If the steady state corresponds to a unique local Lindbladian, the kernel of $K$ will become one-dimensional once sufficiently many constraints are used. We expect this to occur when the number of equations reaches the number of unknowns, revealing the true Lindbladian parameters up to an overall multiplicative constant. When a Lindbladian has multiple steady states, any of them may be used for the reconstruction; however, the reconstruction quality may depend on the steady state used.

Thus, if the elements of $K$ are known exactly, our method recovers a unique Lindbladian whenever the equation $K\vec{c}=0$ has a unique solution. Put differently, the spectrum of singular values of $K$ must contain a single zero. In practice, the elements of $K$ are only known to a finite precision due to measurement noise. The spectrum of $K$ determines the difficulty, or noise sensitivity, of the Lindbladian reconstruction.

Suppose that each observable is only measured to an additive error $\epsilon>0$ 111For example, if each observable is measured experimentally using $n_{s}$ copies of $\rho_{s}$, its expectation value is known up to random noise of order $\epsilon\sim 1/\sqrt{n_{s}}$. For the measured $K$, the equation $K\vec{c}=0$ will likely not have an exact solution. As an approximate solution, we take the normalized coefficient vector $\hat{c}=\vec{c}/\norm{\vec{c}}$ that minimizes $\norm{K\hat{c}}$, i.e. the eigenvector of $K^{T}K$ with smallest eigenvalue. Since the Lindbladian is only recovered up to a multiplicative scalar, we measure the reconstruction error $\delta$ by the $L_{2}$ distance between the normalized coefficient vectors $\hat{c}$ of the recovered Lindbladian and the true Lindbladian,

[TABLE]

Using perturbation theory, we estimated in Ref. Bairey et al., 2019 the reconstruction error due to independent random noise with standard deviation $\epsilon$ added to each element of $K$,

[TABLE]

where $\lambda_{m}$ are the eigenvalues of $K^{T}K$ 222In particular, the reconstruction error is dominated by the gap $\lambda_{1}$ of the constraint matrix, since $\Delta_{est}\leq\epsilon\sqrt{M\lambda_{1}^{-1}}$ (i.e., the squared singular values of $K$).

III Results

III.1 Recovery of random local Lindbladians

We apply our method for the reconstruction of random local Lindbladians from their respective steady states. We start by focusing on chains of $\Lambda=6$ spins with random local interactions and dissipation. We consider Lindbladians of the form given in Eq.(1) with local Hamiltonian terms

[TABLE]

and on-site jump operators $L_{j}$ given by

[TABLE]

We choose open boundary conditions $c_{L,\alpha,\beta}=0$, and draw the remaining Hamiltonian coefficients from a Gaussian distribution with zero mean and unit variance, setting the energy scale for what follows. The real and imaginary parts of the dissipative coefficients $d_{j,\alpha}$ are similarly drawn from a Gaussian distribution, with mean zero and standard deviation $\alpha_{D}=\frac{1}{\sqrt{2}}$.

We obtain the steady state of each random Lindbladian $\mathcal{L}$ by exactly diagonalizing it as a superoperator. We then attempt to recover $\mathcal{L}$ using an increasing number $N$ of constraints $A_{n}$. We start with all the constraints $A_{n}$ acting on single sites and nearest neighbors, and add constraints supported on three consecutive sites in random order. To assess the reconstruction difficulty in practical settings, we add to each measured observable a small, independent, Gaussian noise with mean zero and standard deviation $\epsilon=10^{-4}$. We then compute the reconstruction error $\Delta$ due to the measurement noise $\epsilon$.

As soon as the number of constraints approaches the number of unknowns, the reconstruction error $\Delta$ drops, and we obtain a good approximation of the Lindbladian (Fig. 1a). The error decreases with the number of constraints, following the estimate of Eq. (8). We verified numerically that the reconstruction error $\Delta$ follows the estimate of Eq. (8) over several orders of magnitudes of the measurement noise $\epsilon$, as long as $\Delta\lesssim 10^{-2}$ (Fig. 1a inset).

III.2 Effect of dissipation type and strength

Next, we study how the accuracy of the method depends on the type and strength of the dissipative terms appearing in the Lindbladian. First, we vary the magnitude $\alpha_{D}$ of the dissipative terms appearig in Eq. (10) relative to the Hamiltonian terms. We repeat the recovery experiment on the steady states of these different dynamics, using all 3-local constraints $A_{n}$. We find that the accuracy of the method improves upon adding weak dissipation to a Hamiltonian; the recovery is optimal when the dissipative terms are comparable in magnitude to the Hamiltonian terms (Fig. 1b, red). Due to our choice of single-site jump operators in Eq. (10), steady states at the strong dissipation limit approach product states. Since any product state is a steady state of many different Lindbladians, the reconstruction error diverges for $\alpha_{D}\rightarrow\infty$; this divergence of the error is cured when two-site nearest-neighbor jump operators are added (see Appendix B.1).

In practical situations, the jump operators $L_{j}$ may be unknown even if the Hamiltonian is well-characterized. We can incorporate prior knowledge about the Hamiltonian by turning Eq. (4) into the non-homogeneous constraint

[TABLE]

where the RHS is directly obtained by measurements. The dissipative coefficients $c_{rs}$ are then obtained by solving a system of non-homogeneous linear equations (see Appendix A.4). Fig. 1b shows that recovery with such prior knowledge of the Hamiltonian achieves a lower reconstruction error of the Lindbladian (green curve). Since the recovery with prior knowledge leaves no ambiguity in the magnitude of the Lindbladian, we can also compare the dynamics generated by the true and recovered Lindbladians starting from a fixed initial state; indeed, we find an excellent agreement (Fig. S3).

Next, we study the interplay of different dissipation types. We consider a Lindbladian $\mathcal{L}$ which consists of single-site jump operators of two kinds:

[TABLE]

where $\sigma^{-}\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\frac{1}{2}\left(\sigma^{x}-{\mkern 1.5mu\operatorfont i\mkern 1.5mu}\sigma^{y}\right)$. The “loss” $L_{j,L}$ relaxes the system towards a pure steady state, e.g. due to loss of particles; the “dephasing” $L_{j,D}$ scrambles relative phases between pure states in a specific basis. We tune the parameter $0\leq\alpha_{L}\leq 1$ to interpolate the relative weights of the loss and dephasing. In addition, $\mathcal{L}$ contains Hamiltonian terms of the form (9), with coefficients drawn from a Gaussian distribution with zero mean and unit variance. We then attempt to recover both the Hamiltonian and the jump operators from the steady state of $\mathcal{L}$ using all 3-local constraints $A_{n}$, without assuming that the form of the on-site jump operators is known.

We find that reconstruction of strongly dephasing Lindbladians is hard (Fig. 1c). This is expected: for $\alpha_{L}\ll 1$, the steady state is close to a fully mixed state, compatible with any Lindbladian with Hermitian jump operators. As the loss intensifies, $\norm{\mathcal{L}(\mathbbm{1})}\propto\alpha_{L}^{2}$; correspondingly, Fig. 1c shows that the reconstruction error decreases as $\alpha_{L}^{-2}$ (see also Appendix B.4), indicating that the steady state becomes more informative.

III.3 Loss facilitates learning of commuting Hamiltonians

Motivated by the insight that loss can lead to non-trivial steady states, we investigate whether dissipation can aid in learning Hamiltonians that could not be recovered from their own steady states. In particular, we consider classical Hamiltonians with random nearest-neighbor interactions in the X-basis alone,

[TABLE]

whose coefficients are drawn from a Gaussian distribution with zero mean and unit variance. Any state $\rho$ diagonal in the X-basis is a steady state of $H_{cl}^{x}$, revealing no information about its coefficients. We therefore add on-site jump operators

[TABLE]

so that the dynamics of $\mathcal{L}$ are comprised of Hamiltonian dynamics in the $X$ basis and loss in the $Z$ basis. We then attempt to recover $H$ from the steady state of $\mathcal{L}$, assuming that the jump operators $L_{j}$ are known.

We find that the addition of controlled loss facilitates efficient learning of the classical Hamiltonians of Eq. (13). Due to the small number of unknowns, single-site constraint operators $\sigma_{j}^{y}$, $\sigma_{j}^{z}$ are sufficient to recover $H$ ($\sigma_{j}^{x}$ are not required as they commute with $H$). Moreover, the reconstruction is very robust: when nearest-neighbor constraints are added, the accuracy of the recovered Hamiltonian approaches the measurement accuracy (Fig. 2).

III.4 System-size scaling

Finally, we demonstrate that our method can recover Lindbladians on long spin chains. Various approaches have been proposed for computing steady states of large-scale open quantum systems using matrix product operators  Cui et al. (2015); Mascarenhas et al. (2015); Werner et al. (2016). In this work, we have used the variational MPO approach of  Ref. Mascarenhas et al., 2015, which iteratively finds the density matrix with the smallest-magnitude eigenvalue of $\mathcal{L}$. Using this approach, we obtain steady states of the random Lindbladians considered in Eqs. (9), (10) on chains with $\Lambda=100$ spins (see Appendix C for details).

To study the system-size scaling of our method, we focus on subsystems of increasing sizes. We begin with the $6$ leftmost spins and add $4$ spins in each step, eventually covering the whole chain. We then attempt to recover the Lindbladian of each of these subsystems from observables within that subsystem only, using all 3-local constraints.

We employ two different approaches for recovering the full Lindbladians of these increasingly large subsystems. In the first approach, we partition the subsystem to overlapping patches of $6$ spins, and recover the Lindbladian on each patch independently. The recovery does not determine the overall scale factor of the Lindbladian on the patch; we therefore re-scale the coefficients of neighboring patches according to the coefficients of their shared terms (see Appendix D). In the second approach, we apply our method directly on the whole subsystem, forming a large constraint matrix $K$ which grows with the subsystem size.

Both approaches successfully recover the full-system Lindbladian using the same set of measurements. Here we do not add measurement noise; the error in a single patch ($\approx 10^{-6}$) is controlled by the numerical precision of the MPO steady state. Due to the uncertainty in the coefficients shared between each pair of patches, the norm of the recovered Lindbladian performs a random walk, leading to a total error growing as the square root of the number of patches (Fig. 3, left; see Appendix D for analysis). Namely, the error grows as the square root of system size, $O(\Lambda^{\frac{1}{2}})$. We find the same square root system-size scaling of the reconstruction error in the second, direct approach (Fig. 3, right).

These findings suggest that in order to recover the dynamics of a system of length $\Lambda$ to a fixed accuracy, each observable should be measured to an accuracy of $O(\Lambda^{-\frac{1}{2}})$. In other words, each observable should be measured $n_{s}=O(\Lambda)$ times. The number of observables required scales also as $O(\Lambda)$; however, since they are all local, each copy of the steady state $\rho_{s}$ can be used to measure $O(\Lambda)$ observables. Thus, we expect that $n_{s}=O(\Lambda)$ copies of $\rho_{s}$ overall suffice.

IV Conclusions

Near-term intermediate-scale quantum devices Preskill (2018) are invariably subject to noise and coupled to their environments. While tomographic methods can characterize noises acting on a few isolated qubits Merkel et al. (2013); Blume-Kohout et al. (2017), cross-talk between qubits necessitates holistic methods that identify the sources of error in an entire device Proctor et al. (2019).

Our results suggest that the noises acting on quantum devices may be efficiently characterized from measurements of their steady states. Left to themselves, quantum devices naturally reach their steady states at times longer than their typical relaxation and decoherence timescales. If in addition to single-qubit dissipation, the qubits are also coupled by a Hamiltonian or affected by correlated dissipation, we find that their steady state would be informative enough to recover both the Hamiltonian and the dissipative processes.

In addition to scalability, our approach to characterizing dynamics through their steady states offers a few advantages. It does not require precise control of either state initialization or measurement time. It is independent on the dimensionality of the local Hilbert space, and is effective also for bosonic systems with an infinite-dimensional local Hilbert space. As shown in Fig. 2, addition of controlled terms can allow learning of Hamiltonians consisting of commuting terms, such as those corresponding to topological quantum error-correcting codes Valenti et al. (2019).

Having demonstrated that open quantum system dynamics can generically be learned from their steady states, it is important to obtain rigorous bounds on the number of measurements required for the learning process. Such bounds could be obtained by identifying conditions under which our constraint matrix is guaranteed to be gapped. It could also be interesting to study our method as a means to certify quantum states prepared as the steady states of given quantum dynamics. Finally, adapting our method to the setting of quantum circuits may yield means to certify, characterize and benchmark quantum devices.

Acknowledgements.

We thank Yotam Shapira for useful discussions. E. B. and N. L. acknowledge financial support from the European Research Council (ERC) under the European Union Horizon 2020 Research and Innovation Programme (Grant Agreement No. 639172). D.P. acknowledges support from the Singapore Ministry of Education, Singapore Academic Research Fund Tier-II (project MOE2018-T2-2-142). N. L. acknowledges support from the People Programme (Marie Curie Actions) of the European Union’s Seventh Framework Programme (No. FP7/2007–2013) under REA Grant Agreement No. 631696 and the Defense Advanced Research Projects Agency through the DRINQS program, grant No. D18AC00025. The content of the information presented here does not necessarily reflect the position or the policy of the U.S. government, and no official endorsement should be inferred. I.A. acknowledges the support of the Israel Science Foundation (ISF) under the Individual Research Grant No. 1778/17.

Appendix A Details of the recovery algorithm

A.1 Expanding the Lindblad dynamics in a fixed set of operators: derivation of Eq. (2)

Formally, to derive Eq. (2) from Eq. (1), we first expand each local Hamiltonian term in a fixed basis of local operators

[TABLE]

so that the unitary evolution term becomes

[TABLE]

with

[TABLE]

Similarly, we expand each jump operator in a fixed basis of local operators

[TABLE]

so that the dissipative dynamics may be rewritten as

[TABLE]

where

[TABLE]

forms a positive semi-definite matrix by definition.

A.2 Exact form of the constraint matrix

As derived in Eqs (4-6), the elements of the constraint matrix $K$ are expectation values of different observables. The explicit form of the element $K_{n,m}$ varies, depending on the term in the expansion of the Lindbladian in Eq. (2) which corresponds to the index $m$: (i) coefficients $c_{j}$ of Hamiltonian terms; (ii) diagonal entries of the matrix of dissipative coefficients $c_{rr}$; (iii) the real part of the off-diagonal dissipative coefficients $\frac{1}{2}\Re{c_{rs}+c_{sr}}$; (iiii) the imaginary part of the off-diagonal dissipative coefficients $\frac{1}{2{\mkern 1.5mu\operatorfont i\mkern 1.5mu}}\Re{c_{rs}-c_{sr}}$. Explicitly, the matrix elements $K_{n,m}$ are given by (see also Fig. S1):

[TABLE]

A.3 Eq. (4) for a fully mixed state

From Eq. (1) it is clear that the fully mixed steady state $\rho\propto\mathbbm{1}$ is a steady state of any Lindbladian with Hermitian jump operators $L_{j}=L_{j}^{\dagger}$ (in fact, it is sufficient that the jump operators are normal, $\left[L_{j},L_{j}^{\dagger}\right]=0$). Let us see how this reflects in Eq. (4).

If the dissipators $L_{j}$ are real, we can expand them (see Eq. (S4)) in a basis of Hermitian local operators $l_{r}=l_{r}^{\dagger}$ using real coefficients; subsequently, the coefficient matrix $c_{rs}$ will be real and symmetric. At a fully mixed state, the expectation value of any operator is proportional to its trace, and Eq. (4) becomes

[TABLE]

Since commutators are traceless $\Tr\left[A,B\right]=\Tr AB-\Tr BA=0$, the first part vanishes; in other words, the fully mixed state is a steady state of any Hamiltonian. We now note that the second term is antisymmetric in $r\leftrightarrow s$: using the cyclic properties of the trace $\Tr AB=\Tr BA$ and $\Tr\left(A\left[B,C\right]\right)=\Tr\left(C\left[A,B\right]\right)$,

[TABLE]

which is antisymmmetric to $r\leftrightarrow s$ due to the commutator. On the other hand, $c_{rs}$ is symmetric, so the sum over $r,s$ vanishes:

[TABLE]

Thus, the fully mixed state obeys Eq. (4) for any constraint operator $A$ if the jump operators $L_{j}$ are Hermitian.

A.4 Recovery with prior knowledge

If some part of the dynamics is known to high accuracy, Eq. (4) can be turned into a non-homogenous equation. For instance, if the Hamiltonian is known but the dissipators are not, we obtain Eq. (11). Using a set of constraint operators $\{{A_{n}\}}_{n=1}^{N}$, we obtain the system of equations

[TABLE]

where $K_{l}$ is the constraint matrix of the dissipative operators alone, and $c_{l}$ are their corresponding coefficients; the vector $\vec{b}$ is given by

[TABLE]

Eq. (S12) is then solved using least squares.

Appendix B Error analysis

B.1 Recovery of strongly dissipating Lindbladians

In Fig. 1b, it appears that the recovery error diverges when the relative magnitude of the dissipative terms is large $\alpha_{D}>1$. We conjectured that this divergence does not indicate that recovery is generically impossible in the limit of strong dissipation; rather, it is an artifact of the choice of strictly single-site dissipation we simulated.

To verify this conjecture, we added nearest-neighbor jump operators to our random Lindbladians

[TABLE]

with all coefficients drawn from a Gaussian distribution with mean zero and standard deviation $\alpha_{D}$; for the Hamiltonian terms, we used the same random nearest-neighbor interactions of Eq. (9). We then recovered these Lindbladians from their steady states, assuming that the form of the jump operators is known but their coefficients are not. We found that the reconstruction error of these Lindbladians saturates at large $\alpha_{D}$ (Fig. S2, blue); thus, the divergence of the reconstruction error is cured when entangling jump operators are added.

B.2 Recovery error: results vs. expectation

The recovery error $\Delta$ we find in Fig. 1a is slightly higher (by a factor of $\approx 1.25$) than the estimate of Eq. (8), derived in Ref. Bairey et al., 2019. In contrast to our results in this work, the recovery error obtained in Ref. Bairey et al., 2019 was lower than the prediction of the same estimate, which is indeed expected to be pessimistic due to the use of Jensen’s inequality.

We believe the difference is due to the different noise model used in both papers: here we add noise to each measured observable, while in Ref. Bairey et al., 2019 we added independent noise to each of the entries of $K$ (even when they contain the same observable). This is because in Ref. Bairey et al., 2019, we wished to test the theoretical validity of the error estimate. The estimate assumes that the noise in each entry of the constraint matrix $K$ is independent, and we thus added an independent random noise to each of its entries. Realistically though, noise is incurred in each measured observable. Since many different entries of $K$ feature the same observable, this introduces correlations between the noise in different entries.

B.3 Accuracy of the reconstructed dynamics

To assess how well the recovered dynamics approximate the true dynamics, we compared the time evolution generated by the recovered and true Lindbladians starting from a fixed initial state. We focused on random Lindbladians with a relative dissipation magnitude $\alpha_{D}=\frac{1}{\sqrt{2}}$ and a known Hamiltonian, exactly as in Fig. 1b (green curve). The knowledge of the Hamiltonian allows to recover the Lindbladian exactly (including its overall magnitude), allowing a meaningful comparison of time dynamics.

We initialized the system in a product state with all spins up,

[TABLE]

and computed its evolution under the true Lindbladian $\rho(t)$ and under the recovered Lindbladian $\rho_{rec}(t)$. At each point in time, we compared these two states by the average trace distance between their reduced density matrices on pairs of consecutive sites,

[TABLE]

where $\rho^{(i,j)}=\Tr_{\Lambda\setminus\{{i,j\}}}\left(\rho\right)$ is the reduced density matrix on sites $i,j$, and the trace distance

[TABLE]

bounds the difference in the expectation value of any POVM element. Thus, $D_{loc}(\rho,\rho_{rec})$ is a worst-case measure for the difference between local observables in the two states.

As shown in Fig. S3, the mean local trace distance peaks at a value below $10^{-3}$ for short times. It then decreases to $\approx 2\cdot 10^{-4}$, which is approximately the measurement accuracy taken for the reconstruction. This is not surprising in retrospect: at long times, $\rho_{rec}$ is the steady state of the recovered Lindbladian, which was chosen such that the measured local observables would correspond to its steady state.

B.4 Scaling of the reconstruction error with the relative weight of loss in the dissipation

We argued that Fig. 1c confirms the theoretical expectation that the reconstruction error scale as $\alpha_{L}^{-2}$ when the weight of loss relative to dephasing $\alpha_{L}$ [see Eq. (12)] is small. However, the curve in Fig. 1c did not show a clear power law for small $\alpha_{L}$, since the reconstruction error approached large values of order $1$. We thus repeated these simulations with weaker measurement noise ($\epsilon=10^{-8}$ compared to $\epsilon=10^{-4}$ in the main text), and verified this power law over a wider range of $\alpha_{L}$ (Fig. S4).

Appendix C Computing the steady state using variational matrix product operators

The steady state of the Lindbladian can be obtained by computing the eigenstate of the Lindblad operator $\mathcal{L}$ corresponding to eigenvalue [math]Cui et al. (2015); Mascarenhas et al. (2015) (the system we studied has no degeneracy). Internally the density operator $\rho$ is reorganized into a long vector and treated similarly to the state vector of a unitary system. We use the variational matrix product operator algorithm proposed in Ref. Mascarenhas et al., 2015, where the iterative procedure to search for the steady state is done in the same way as the unitary case, except that one keeps the eigenstate corresponding to the eigenvalue with smallest magnitude instead of smallest algebraic value. For a system with $100$ spins, we have used $D=100$ number of states (bond dimension) for our simulations. We obtained a steady state with eigenvalue of the order $10^{-8}$, and residual $\norm{\mathcal{L}\rho}\approx 10^{-5}$. To check the convergence against different $D$s, we have done another simulation with $D=150$, and compared the local observables $\langle\sigma^{z}_{j}\rangle$, obtaining a mean error $\sum_{j=1}^{\Lambda}\absolutevalue{\langle\sigma^{z}_{j}\rangle_{D=100}-\langle\sigma^{z}_{j}\rangle_{D=150}}/\Lambda\approx 10^{-9}$. We also compared the distances between the reduced density matrices with a patching size $6$, and a patching spacing $4$, and obtained a mean error of the order of $10^{-8}$.

Appendix D Stitching up recovered patches

Recall that the Lindbladian on each patch is only recovered up to a multiplicative scalar. Suppose we recover the Lindbladian of two overlapping patches and wish to “stitch” them together into one Linbladian acting on the joint patch. In the absence of noise, the recovered Lindbladians of the first two patches would be given by

[TABLE]

where $\mathcal{L}_{m}$ is the vector of terms [$h_{j}$ and pairs $(l_{r},l_{s})$] acting on the overlapping region of the two patches; for the analysis below, we assume that each individual recovered Lindbladian is normalized: $\norm{\vec{c}_{l}}^{2}+\norm{\vec{c}_{m}}^{2}=\norm{\vec{c}\mkern 2.0mu\vphantom{c}^{\prime}_{m}}^{2}+\norm{\vec{c}_{r}}^{2}=1$. The coefficients $c_{m},\vec{c}\mkern 2.0mu\vphantom{c}^{\prime}_{m}$ of the overlapping region will generically differ since the Lindbladian on each patch is only recovered up to a multiplicative scalar. We therefore use these overlapping coefficients to determine the relative scale of the two patches, by multiplying the Lindbladian of the second patch by a factor of $\frac{\norm{\vec{c}_{m}}}{\norm{\vec{c}\mkern 2.0mu\vphantom{c}^{\prime}_{m}}}$:

[TABLE]

In fact, we also need to fix the relative signs of the two patches using a similar factor of $\frac{\text{sign}({\vec{c}_{m}})}{\text{sign}({\vec{c}\mkern 2.0mu\vphantom{c}^{\prime}_{m}})}$, where the sign can be determined e.g. according to the coefficient of a fixed shared term. While this last detail is crucial for the stitching process, it does not contribute to the recovery error due to noise, as long as the error in each patch is small relative to its size, so that no coefficient flips its sign.

To recover the Lindbladian of a sequence of patches $1,\dots,n$, we repeat this procedure iteratively and obtain

[TABLE]

where

[TABLE]

with $\mathcal{L}_{1,2}$ denoting the terms acting on the overlapping region of the first two patches. For any $j>1$,

[TABLE]

If each individual patch is recovered perfectly up to a corresponding multiplicative scalar, this procedure yields the full system Lindbladian up to a single overall multiplicative scalar. However, noise introduces error in the recovered Lindbladian of each individual patch: $\vec{c}_{j}\mapsto\vec{c}_{j}+\vec{\delta}_{j}$.

Error in each individual patch affects the overall stitched Lindbladian in two ways. One effect is a rotation of each $\mathcal{L}_{patch}$ with respect to its true value, the $\mathcal{L}_{j}$ component pointing to $\vec{c}_{j}+\vec{\delta}_{j}$ rather than $\vec{c}_{j}$. Since this error is additive, it is absorbed in the normalization of $\mathcal{L}_{stitched}$; assuming that the error is approximately uniform across patches, $\norm{\vec{\delta}_{j}}\approx\delta$, it leads to an overall error of order $\delta$ in the total $\mathcal{L}_{stitched}$, which is independent of the number of patches.

A second effect caused by the errors in the recovery of individual patches is a stretch of each $\mathcal{L}_{patch}$. This effect is induced through the errors’ effect on the relative scale factor $\prod_{i=1}^{j-1}\frac{\norm{\vec{c}_{i,i+1}}}{\norm{\vec{c}\mkern 2.0mu\vphantom{c}^{\prime}_{i,i+1}}}$. Assuming that the errors of the different patches $\vec{\delta}_{j}$ are independent, this scale factor performs a multiplicative random walk, fluctuating from its true value by a deviation of order $\sqrt{j}\delta$. This is most easily seen by taking a log:

[TABLE]

To first order in $\delta$, each of these is an independent random variable with zero mean and standard deviation of order $\delta$:

[TABLE]

where $\hat{c}_{i,i+1}=\vec{c}_{i,i+1}/\norm{\vec{c}_{i,i+1}}$. Therefore, the ratio between the true scale factor and its noisy version is given by $e^{\tilde{\delta}}$, where $\tilde{\delta}$ is the random variable given by Eq. (S24). Its standard deviation scales as $\sqrt{j}\delta\leq\sqrt{n}\delta$, where $n$ is the total number of patches. While the order $\delta^{2}$ correction is always positive, resulting in a drift, it sums up across the patches to $O(n\delta^{2})$, and is therefore higher order in $\sqrt{n}\delta$. Thus, as long as $\sqrt{n}\delta\ll 1$, the Lindbladian on each patch is stretched by a factor of at most $\approx 1\pm\sqrt{n}\delta$, leading to a total recovery error of order $\sqrt{n}\delta$. This explains the square root scaling of the error with system size seen in Fig. 3.