Learning the dynamics of open quantum systems from their steady states

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.

Previous works have recovered the dynamics of open quantum systems by tracking their time evolution [41-47, 18, 48-51]. 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 [52,53]: 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 ρ s of  is defined by ( )  r r = =  0 s s . Suppose that we prepare many copies of ρ s and measure expectation values of local observables in the state ρ s . Can  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 | | ñá   i i . In contrast, generic Lindbladians have only a single steady state [54]. 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  r µ , 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 [40], 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.

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 {( )} l l , r s . Expanding the dynamics in this operator basis (see appendix A.1), equation (1) becomes 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 l , r 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 ρ s of . Since ρ s is a steady state, the expectation value ( ) r á ñ = A A Tr s def of any observable must be time-independent where the expectation values are taken with respect to the steady state ρ s . For any operator A, equation (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 equation (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 equation (4). We concatenate the Hamiltonian parameters c j and the dissipative parameters c rs into a single vector  c . In this notation, equation (4) takes the form . Thus,  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 N 1 , we obtain a homogeneous system of linear equations for the coefficients of the true Lindbladian where K is anŃ 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 equation (2). Assuming that we measured K at a steady state of a local Lindbladian, the vector  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  = Kc 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 >  0 4 . For the measured K, the equation  = Kc 0 will likely not have an exact solution. As an approximate solution, we take the normalized coefficient vectorˆ|| ||   = c c c that minimizes ||ˆ|| Kc , i.e. the eigenvector of K K T with smallest eigenvalue. Since the Lindbladian is only recovered up to a multiplicative scalar, we measure the reconstruction error d by the L 2 distance between the normalized coefficient vectors ĉof the recovered Lindbladian and the true Lindbladian, Using perturbation theory, we estimated in [40] the reconstruction error due to independent random noise with standard deviation ò added to each element of K

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 L = 6 spins with random local interactions and dissipation. We consider Lindbladians of the form given in equation (1) with local Hamiltonian terms and on-site jump operators L j given by

We choose open boundary conditions
, 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 a d j, are similarly drawn from a Gaussian distribution, with mean zero and standard deviation a = D 1 2 . We obtain the steady state of each random Lindbladian  by exactly diagonalizing it as a superoperator. We then attempt to recover  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 4 For example, if each observable is measured experimentally using n s copies of ρ s , its expectation value is known up to random noise of order n 1 s . 5 In particular, the reconstruction error is dominated by the gap l 1 of the constraint matrix, since 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 = - 10 4 . We then compute the reconstruction error Δ due to the measurement noise ò.
As soon as the number of constraints approaches the number of unknowns, the reconstruction error Δ drops, and we obtain a good apparoximation of the Lindbladian ( figure 1(a)). The error decreases with the number of constraints, following the estimate of equation (8). We verified numerically that the reconstruction error Δ follows the estimate of equation (8) over several orders of magnitudes of the measurement noise ò, as long as D - 10 2 (figure 1(a) inset).

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 a D of the dissipative terms appearig in equation (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 ( figure 1(b), red). Due to our choice of single-site jump operators, 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 a  ¥; D 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 equation (4) into the non-homogeneous constraint 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). Figure 1(b) 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 ( figure B2).
Next, we study the interplay of different dissipation types. We consider a Lindbladian  which consists of single-site jump operators of two kinds: Figure 1. Reconstruction of Lindbladians from their steady states. We generated steady states of random local Lindbladians on chains of L = 6 spins and measured local observables given by equation (4) for a set of constraint operators { } = A n n N 1 . We then recovered the Lindbladians from these observables by solving equation (6), adding a small random measurement noise of order = - 10 4 to each observable, and computed the error Δ in the recovered Lindbladians (see equation (7)). (a) Reconstruction error (Δ) of random Lindbladians (equations (9), (10)) as a function of the number of constraints (red; shaded area indicates error bars). Recovery succeeded once the number of constraints N approached the number of unknowns M (here M = 117); its accuracy improved as more constraints were added, following the estimate D est (dashed line) of equation (8). Here, the ratio between the magnitudes of the Hamiltonian and dissipation was fixed to a = D 1 2 . Inset: the error-to-noise ratio D  with all 3-local constraints as a function of the measurement noise magnitude ò. The error followed the prediction of equation (8)  Here we used all constraints A n acting on up to 3 consecutive sites. Addition of weak dissipation improved the Lindbladian recovery, which was optimal at a » 0.5 D . A lower reconstruction error was achieved when the Hamiltonian was known (green; D prior ). (c) Dependence of the reconstruction error on the type of dissipation. We used the same ensemble of random Hamiltonians, with dissipation given by equation (12), and α L interpolating between loss and dephasing. When dissipation is almost entirely due to dephasing, a  0 L , the steady state is close to being fully mixed; consequently, recovery improves with increasing loss (increasing α L ). All results were averaged over 300 random Lindbladians, with error bars indicating one standard deviation; means and standard deviations were calculated after taking the log.
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 a   0 1 L to interpolate the relative weights of the loss and dephasing. In addition,  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  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 ( figure 1(c)). This is expected: for  a 1 L , the steady state is close to a fully mixed state, compatible with any Lindbladian with Hermitian jump operators. As the loss intensifies,  || ( )|| a µ  ; L 2 correspondingly, figure 1(c) shows that the reconstruction error decreases as a -L 2 (see also appendix B.4), indicating that the steady state becomes more informative.

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, whose coefficients are drawn from a Gaussian distribution with zero mean and unit variance. Any state ρ 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 so that the dynamics of  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 , 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 equation (13). Due to the small number of unknowns, single-site constraint operators s j y , s j z are sufficient to recover H (s 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 (figure 2).

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 [55][56][57]. In this work, we have used the variational MPO approach of  [56], which iteratively finds the density matrix with the smallest-magnitude eigenvalue of . Using this approach, we obtain steady states of the random Lindbladians considered in equations (9), (10) on chains with L = 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 (» -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 (figure 3, left; see appendix D for analysis). Namely, the error grows as the square root of system size, ( ) L O 1 2 . We find the same square root system-size scaling of the reconstruction error in the second, direct approach ( figure 3, right).
These findings suggest that in order to recover the dynamics of a system of length Λ to a fixed accuracy, each observable should be measured to an accuracy of ( ) L -O

Conclusions
Near-term intermediate-scale quantum devices [1] are invariably subject to noise and coupled to their environments. While tomographic methods can characterize noises acting on a few isolated qubits [58,59], cross-talk between qubits necessitates holistic methods that identify the sources of error in an entire device [60].
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. Figure 3. Reconstruction of Lindbladians on large spin chains: system-size scaling of the reconstruction error. We obtained the steady states of the random Lindbladians described in equations (9), (10) on Λ=100 spins. (left) We recovered the Lindbladians on spatial patches of 6 spins, with overlaps of 2 sites between consecutive patches. We used all constraints supported on up to 3 consecutive sites in the interior of each patch (middle 4 sites for bulk patches). We then stitched consecutive patches to obtain the full Lindbladian on subsystems of increasing length. The reconstruction error increased with system size (red line), following the predicted square-root scaling with the number of patches (dashed line). (right) As a different approach, we built a single large constraint matrix for each subsystem, and obtained the error as a function of subsystem size; this approach yielded a slightly smaller reconstruction error, still scaling as the square root of subsystem size (dashed line).
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 figure 2, addition of controlled terms can allow learning of Hamiltonians consisting of commuting terms, such as those corresponding to topological quantum errorcorrecting codes [61].
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. through the DRINQS program, . Explicitly, the matrix elements K n m , are given by (see also figure A1): ). Let us see how this reflects in equation (4).
If the dissipators L j are real, we can expand them (see equation (A4)) in a basis of Hermitian local operators † = l l r r 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 equation (4) Figure A1. Top: we concatenate the Hamiltonian coefficients  c h and the matrix of dissipative coefficients C d into a long vector of coefficients for the Lindblad evolution. Off-diagonal entries of C d are split into their real and imaginary parts (blue and magenta, correspondingly). Bottom: the constraint matrix is composed of a vertical block corresponding to Hamiltonian terms K h , and a vertical block corresponding to dissipative terms K d . Entries corresponding to Hamiltonian terms are given by their commutators with the constraint operators (red); the formula for the dissipative entries varies between the diagonal entries of the dissipative matrix C d (green), and the real (blue) and imaginary (magenta) entries of C d . A.4. Recovery with prior knowledge If some part of the dynamics is known to high accuracy, equation (4) can be turned into a non-homogenous equation. For instance, if the Hamiltonian is known but the dissipators are not, we obtain equation (11). Using a set of constraint operators { } = A n n N 1 , we obtain the system of equations where K l is the constraint matrix of the dissipative operators alone, and c l are their corresponding coefficients; the vector  b is given by n n Equation (A12) is then solved using least squares.
Appendix B. Error analysis B.1. Recovery of strongly dissipating Lindbladians In figure 1(b), it appears that the recovery error diverges when the relative magnitude of the dissipative terms is large a > 1 D . 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 with all coefficients drawn from a Gaussian distribution with mean zero and standard deviation a ; D for the Hamiltonian terms, we used the same random nearest-neighbor interactions of equation (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 a D (figure B1, blue); thus, the divergence of the reconstruction error is cured when entangling jump operators are added.

B.2. Recovery error: results versus expectation
The recovery error Δ we find in figure 1(a) is slightly higher (by a factor of »1.25) than the estimate of equation (8), derived in [40]. In contrast to our results in this work, the recovery error obtained in [40] 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 [40] we added independent noise to each of the entries of K (even when they contain the same observable). This is because in [40], 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 a = is the reduced density matrix on sites i j , , and the trace distance loc rec is a worst-case measure for the difference between local observables in the two states.
As shown in figure B2, the mean local trace distance peaks at a value below 10 −3 for short times. It then decreases to · » -2 10 4 , which is approximately the measurement accuracy taken for the reconstruction. This is not surprising in retrospect: at long times, r rec is the steady state of the recovered Lindbladian, which was chosen such that the measured local observables would correspond to its steady state. Figure B1. Entangling jump operators facilitate learning of strongly dissipative dynamics. Reconstruction error as a function of dissipation strength a D of Lindbladians with nearest-neighbor Hamiltonian terms (equation (9)); for the dissipation, we took either strictly single-site jump operators (red), as in main text (equation (10)); or both single-site and nearest-neighbor jump operators (blue) (equation (B1)). The divergence of the error at the strong dissipation limit is cured when nearest-neighbor jump operators are added. Here, we added a smaller noise than in the main text ( = - 10 8 rather than 10 −4 ) to probe the behavior at large values of a D . Figure B2. The accuracy of the evolution generated by the reconstructed dynamics as a function of time. We initialize the system in a product state with all spins up. We then measure the deviation between its evolution by the true dynamics ( ) r t and its evolution by the recovered dynamics r rec by the mean local trace distance (red curve, see equation (B3)). For comparison, we also show the mean local trace distance between the true dynamics and the fully mixed state (green curve).

B.4. Scaling of the reconstruction error with relative weight of loss in the dissipation
We argued that figure 1(c) confirms the theoretical expectation that the reconstruction error scale as a -L 2 when the weight of loss relative to dephasing α L (see equation (12)) is small. However, the curve in figure 1(c) did not show a clear power law for small α L , since the reconstruction error approached large values of order 1. We thus repeated these simulations with weaker measurement noise ( = - 10 8 compared to = - 10 4 in the main text), and verified this power law over a wider range of α L (figure B3).