Hamiltonian Assignment for Open Quantum Systems

We investigate the problem of determining the Hamiltonian of a locally interacting open-quantum system. To do so, we construct model estimators based on inverting a set of stationary, or dynamical, Heisenberg-Langevin equations of motion which rely on a polynomial number of measurements and parameters. We validate our Hamiltonian assignment methods by numerically simulating one-dimensional XX-interacting spin chains coupled to thermal reservoirs. We study the effects of systematic noise and discuss connections with alternative Hamiltonian tomography formulations.


I. INTRODUCTION
Fault tolerant quantum computation provides a framework for digitally decomposing unitary operators using a polynomial number of one-and two-qubit operations drawn from a universal gate set [1].For noisy intermediate scale quantum (NISQ) hardware, characterized by fixed gate fidelities and limited coherence times, digitizing a quantum simulation unitary is too costly in terms of the polynomial scaling digital gate complexity and circuit depth.
However, if a programmable quantum device's manybody dynamics are described by an underlying Hamiltonian H, it is prudent to consider digital-analog decompositions [2] leveraging H.It has been proposed that, in such a case, the target unitary can be decomposed as a sequence of native analog unitaries U = exp(−iHt) interleaved with programmable single-qubit operations.For certain applications, such as many body simulation [3], the gate complexity, quantified by the total number of applications of the many body evolution operator U and local rotations, may be significantly smaller than that of the resulting digitized decomposition.
To ensure accuracy the digital-analog quantum simulation's error must be bounded, e.g. in terms of the distance between the target and digital-analog unitaries.In order to upper bound the simulation error, one must first precisely characterize the device Hamiltonian generating the many body-entangler.
While some prior results regarding Hamiltonian estimation exist (e.g.process tomography [4] and Bayesian Hamiltonian learning [5]) of the ones that are efficient (local Hamiltonian tomography [6,7]) the estimation task is complicated by unwanted interactions coupling the system of interest to environmental degrees of freedom.To address this outstanding issue, we study the problem of assigning a Hamiltonian to an open quantum system, provided the principle system and environmental interactions are geometrically local.
This work is organized as follows.We first formulate the task of Hamiltonian inference and summarize previous results.Next, we discuss how the Hamiltonian learning protocols, involving a correlation matrix inversion, may be generalized to the context of open quantum systems.Finally we perform numerical simulations in order to validate the techniques and analyze their behaviour in a noisy setting.Finally, we discus generalizations and future directions for Hamiltonian learning.

II. BACKGROUND
Hamiltonian tomography refers to the task of estimating a Hamiltonian H given access to states whose dynamics is generated by H.While this task is infeasible (exponentially costly) in the most general case, the tomography of local Hamiltonians has recently attracted significant attention [6][7][8][9] due to its scalability.Specifically, we are interested in determining k-local Hamiltonians of the form where each S i is an operator supported on k spatially connected sites.We work in the Pauli basis, such that all k-local operators can be written as where µ j ∈ {I, X, Y, Z}, i.e. as a tensor product over Pauli operators.

A. Stationary Learning
for some the observable O. Inserting Eq. 1, selecting an operator basis {O j }, and measuring the commutators [O j , H] we may write the set of linear equations j,i [O j , S i ] c i = 0.These equations can be concisely expressed in matrix form as where we have introduced the matrix A, with matrix elements A i,j = [O j , S i ] , and the Hamiltonian coefficient Note that in principle A need not be a square matrix as its dimensions are determined by the number of accessible correlation measurements.Since the operators S i are k-local and we have the freedom to choose O j from a local basis, most correlators will vanish due to spatially non-overlapping (O j , S i ) pairs and A will be sparse.
In practice, entries of A arise from noisy measurements which may lead to an erroneous evaluation of eigenvalues.To improve numerical stability of the inversion problem, one could reformulate the problem as a convex optimization problem minimize A c 2 2 which is equivalent to maximizing a Gaussian log likelihood, maximize log e − c T A T A c .The latter formulation is convenient for incorporating Bayesian uncertainty quantification methods e.g., to treat noise in the matrix A. If A is a square matrix than A c = 0 has a unique solution only if A has a non-degenerate zero eigenvalue and associated the eigenvector c is the Hamiltonian estimator.Note that Eq. 3 is actually true for the family of Hamiltonians defined up to a scalar factor, i.e. for a scalar a, |ψ n is also an eigenvector of aH.
In general, we wish to avoid the trivial solution and find a minimizing Hamiltonian with c 2 > 0. The solution to the constrained minimization problem minimize A c 2 2 , subject to c 2 2 = 1 is is the row vector of V T associated with the minimal singular value found by the singular value decomposition (SVD) A = U ΣV T .In Sec.IV A we discuss the the numerical stability of this method.

B. Dynamical Learning
While the homogeneous operator equations derived from steady states provide a simple inferential formalism, preparing eigen-and thermal-states of an unknown Hamiltonian may be a challenging or time consuming task in and of itself.Earlier work has therefore also explored Hamiltonian estimation in a dynamical context [8,9].In the Heisenberg picture an observable O i evolves as O i (t) = e itH O i (0)e −itH where O i (0) denotes the observable at time t = 0 where it coincides with the Schrodinger picture counterpart.Taking Eq. 2 and integrating over an infinitesimal time interval δt we get, where the trace is taken with respect to an initial state ρ j , which serves as an indexed input degree of freedom.Considering the small time evolution of a set of operators O i , with respect to a set of initial states ρ j , we may construct the a set of heterogeneous Heisenberg equations of motion defined by the measurement settings vector

which can be arranged as,
As before, the assigned Hamiltonian will correspond to the solution vector c optimizing min .

III. OPEN SYSTEMS GENERALIZATION
Unfortunately, the inability to evolve by purely unitary dynamics limits the applicability of the prior methods.In contrast, realistic quantum systems are open and, in the presence of unknown environmental interactions, evolve dissipatively.In order to incorporate environmental couplings in our inferential framework, we consider the generalized master equation dynamics for a density operator ∂ t ρ = L[ρ] generated by the quantum Liouvillian L. As a concrete example, let us consider Lindblad equation given by L Motivated by physical locality, we take the {L n } operators form an orthonormal basis spanning the manifold of J-local superoperators.Furthermore, the coefficient matrix γ is constrained to be positive semi-definite in order to be a valid mapping on the space of positive semi-definite density operators [10].An observable's dynamics will now be given by the Heisenberg-Langevin master equation . Lindbladian fixed points, satisfying L[ρ] = 0, generalize the notion of Hamiltonian steady states.The fixed-point's time-translational invariant dynamics can be used to generalize Eq. 3 into a set of linear equations written as A c + B γ = 0 where B is defined its elements As before, the set of linear equations can be invoked as C x = 0, where C = (A|B) acts on a composite Hamiltonian-Lindbladian model space spanned by the configuration vector x = ( c T , γ T ) T .Without access to Lindbladian fixed points, we again turn our attention to finite-time evolution.Like before, the linearized Heisenberg-Langevin evolution can be constructed at the cost of an approximation error δt 2 .Aside from this approximation the only difference, with respect to the steady state methods, is the appearance of additional input state ρ j degree of freedom, such that the matrix elements are defined as The dynamical equations of motion are simply C x = W where C = (A |B ), and A has been defined in the Hamiltonian context.

IV. NUMERICAL SIMULATIONS AND MEASUREMENT COMPLEXITY
We present numerical simulations in order to quantitatively analyze the efficacy the open-systems learning protocols outlined above.As an illustrative example, we consider 1-dimensional spin chain consisting of N sites whose Hamiltonian is given by the sum H = H 0 + H I .Noisy equilibrium learning for N = 5 using the parameters of Fig. 1 in the presence of N (0, σ) distributed measurement noise.The error in model estimators grows linearly as a function of σ.Note the estimator error is approximately two orders of magnitude larger than the minimal singular value when the gap ∆s is well defined, and the estimate is inaccurate in the regime of a poorly definite gap or degenerate singular spectrum.Results are averaged over 20 disorder configurations.
Here H 0 = i c i • σ i consist of single qubit interactions, with σ i = (σ x i , σ y i , σ z i ), and the nearest neighbor interactions are given by H I = i J i,i+1 σ x i σ x i+1 .We simulate a translationally invariant system with interactions c = (−2.55,0, 0.5), J = 0.25 for all sites.Although such simplifications could dramatically reduce the resource requirements, we work with the full, unconstrained model space which may assign distinct parameters to each region.To mimic the environmental effects, each spin is subjected to thermal noise.We choose this model as its thermal dissipators are not diagonal in the Pauli basis.More precisely, the single qubit thermal excitation and relaxation operators are typically expressed as L + = √ g + σ + and L − = √ g − σ − where g + = gn/2, g − = g(n + 1)/2, here n is a thermal occupation number, and g the coupling to the environmental reservoir.Given an operator O, assumed to be partially supported on a thermally coupled site, O evolves dissipatively as Expanding the ladder operators σ ± = (X ± iY )/2 and re-grouping the terms we see that this expression may be re-written in the Pauli basis as where the coefficients are

A. Equilibrium learning
In order to the simulate steady state learning protocol we numerically solve for the fixed point density operator, i.e. the solution to L[ρ] = 0 where the Lindblad operator L has been defined above [11].Next, we select a basis for the Hamiltonian and Lindbladian model spaces, as well as for the set of input operators {O i }.Our choice of candidate learning models corresponds to K = 2, J = 1local Hamiltonians/Linbladians.Applying this approach in general, one should begin with a small model space and increase it until adequate convergence.In order to fully explore the maximal quality of learning model, and to deal with finite size scaling corrections (which are especially important for small system sizes) we take {O i } to be the overcomplete union of all 1, 2, 3, 4-local operators.
As illustrated by the dashed lines in Fig. 1 we sweep through the {O i } operators at each level of locality, where they are enumerated beginning with those having support on the left end of the chain until the end of the chain.The top panel in Fig. 1 Our normalized estimator is constructed by performing a singular value decomposition on the correlation matrix C = U ΣV T and selecting as a solution the row vector of V T associated with the minimal singular value of Σ = diag(s m−1 , • • • , s 0 ).Fig. 1 shows how both estimator errors rapidly converges to zero as the dimensionality of the input {O i }space surpasses the model dimensionality, i.e. when C becomes left invertible.At this point the minimal singular value s 0 simultaneously approaches zero, such that Eq. 3 is well approximated by the associated right singular vector.Our solution is unique as indicated by the large singular gap ∆ s = s 1 − s 0 , in the sense that ∆ s /s 0 >> 1.
In a realistic setting stochastic, measurement, and modeling errors will deteriorate the estimator quality.In order to quantify these effects, we add normally distributed random noise with mean zero and variance σ 2 to the elements of the correlation matrix.Fig. 2 illustrates the noise's impact on the correlation matrix properties and estimator errors.Overall, the estimator error grows linearly in σ, up to a critical point ∆ s ≈ s 0 where the solution is degenerate.Before this point, ∆ s remains roughly constant while s 0 rises linearly, and that the estimator infidelity rises approximately proportionately to the characteristic ratio ∆ s /s 0 .

B. Dynamical learning
Lastly, we simulate the dynamical learning scenario whose solution satisfies min x ||C x − W | 2 2 .Refs. 8, 9 previously suggest a set of input product states from which the state may evolve and its dynamics computed.Since product states may not actually be available as a re-  source, especially for near term hardware, we consider a slightly modified scenario.That is, let us take the fixed-point density operator as a resource and conjugate it with respect to a set of unitary operators {U j }.In this way, we generate the set of approximate product states {ρ j } = {U j ρ ss U † j }.Choosing {U j } to consist of all 1 − local Pauli operators, and beginning a time step δt/c z ∼ 1e − 3, we find that the finite time protocol behaves overall similarly to that of the steady states.That is, the estimator error dramatically vanishes when the set of equations is complete with respect to the model space.This is followed by a slower, less dramatic convergence as more information is collected.The normalized minimum estimator error for both the Hamiltonian and the total Hamiltonian-Lindbladian estimates is plotted in Fig. 3 as a function of system size.Here we illustrate the errors for both components of the model system using i) a naive least squares fit and ii) imposing positivity constraints on the Lindbladian process.Interestingly Fig. 3 shows that the Hamiltonian estimator error is insensitive to this distinction whereas the environmental error is greatly reduced with the physical constraints.
One potential drawback of the current approach is that the small time resolution needed for the linear finite time approximation may be difficult to achieve in practice.In order to overcome this obstacle, we consider a slight modification and replace the first order finite difference approximation with a higher order approximation ) such that the error in the equations of motion is O(δt 3 ).The bottom panel of Fig. 3 the Hamiltonian and composite estimator errors as a function of δt for both the linear and quadratic finite difference approximations.Note one major drawback of this modification lies in the increased variance for the derivative estimator.Denoting V 1(2) as the variance in evaluating the first (second) order time derivative of an operator O, and assuming independence and that Var[O(mδt)] is constant for all m, we have V 2 /V 1 = 13/4 which corresponds to a ∼ 10X increase in the number of samples to reduce the second order estimator variance to that of the first.It is also of interest to note that in the large time regime the least squares residual error for the Hamiltonian estimator is smaller than that of the convex solution.One possible explanation for this phenomena is that imposing physicality in the presence of a large approximation error biases the estimator away from the true model while the least squares minimization is afforded more freedom in obtaining a Hamiltonian estimator.

V. CONCLUSION
In this work we have studied the task of assigning a local Hamiltonian to open-quantum systems.By restrict-ing ourselves to a Lindbladian formulation, with a polynomial number of model parameters describing the evolution, we are able to generalize and implement previous local model estimation techniques in both a steady state and dynamical setting.
To validate our constructions, we have performed numerical simulations of open-systems Hamiltonian assignment in the context of an XX-interacting spin chain subject to thermal relaxation.Our results verify how in the clean limit and, contingent upon an appropriately chosen model space, an approximate Hamiltonian may inferred both in-and out-of-equilibrium.Furthermore, we have considered the effects of noise and simple modifications which may increase estimator accuracy in the dynamical context.
Our work paves the way for the determination of manybody Hamiltonians in open quantum systems.The extension of our results to a greater diversity of openquantum systems remains an open research direction.For example, it will be of great interest to understand how our current work can be enhanced by incorporating existing techniques such as Bayesian learning [5], and quantum process identification [12].
The authors acknowledge DOE ASCR funding under the Quantum Computing Application Teams program, FWP number ERKJ347.This research used quantum computing system resources supported by the U.S. Department of Energy, Office of Science, Office of Advanced Scientific Computing Research program office.The authors thank A. Seif, P. T. Bhattacharjee, and R. S. Bennink for many helpful conversations.
Suppose we have access to either i) eigenstates |ψ n of H, with H|ψ n = E n |ψ n or ii) thermal states ρ = arXiv:1911.11092v1[quant-ph] 25 Nov 2019 exp(−βH)/Z where Z = Tr[exp(−βH)] is the partition function.Aside from the trivial phase factors acquired, both classes of states are stationary under the parent Hamiltonian's dynamics.Transforming to the Heisenberg picture, where observable quantities are likewise stationary, we may write

FIG. 1 .
FIG.1.Top panel: Steady state estimator error and singular value properties as a function of the cardinality of the set |{Oi}|, i.e. the number of rows in Eq. 3).The normalized Hamiltonian(-Lindbladian) estimator error ∆H (∆H−L) is given by the blue x's (orange triangles) and the minimal singular value and gap to the next largest singular value are denoted by the green solid and red dashed lines respectively.Note the abrupt decrease in estimator accuracy and singular spectrum upon the correlation matrix becoming left invertible (our 5-site model contains 141 free parameters).Bottom panel: Measurement complexity |{Mi}|, in terms of the cardinality of the set of measurement operators which generated as the union of unique operators contained in the matrix elements of C. The the local, XX, and thermal interactions given parameters c = (−2.55,0, 0.5), J = 0.25, g = 0.05 respectively and the gray dashed lines indicate the positions at which 1,2,3, and 4-local input operators are exhausted.
FIG. 2.Noisy equilibrium learning for N = 5 using the parameters of Fig.1in the presence of N (0, σ) distributed measurement noise.The error in model estimators grows linearly as a function of σ.Note the estimator error is approximately two orders of magnitude larger than the minimal singular value when the gap ∆s is well defined, and the estimate is inaccurate in the regime of a poorly definite gap or degenerate singular spectrum.Results are averaged over 20 disorder configurations.
illustrates an evaluation of the fixed-point model estimation protocol for a N = 5 site chain with open boundary conditions.The key quantities of interest are the error in the Hamiltonian estimator ∆H = || c T − c E || 2 /|| c T || 2 , where T, E refer to the true and estimator models, and, recalling that x = ( c T , γ T ) T , the error in the total estimator ∆

1 FIG. 3 .
FIG. 3. Top panel: Normalized Hamiltonian and Hamiltonian-Lindbladian reconstruction errors as function of chain length using least squares minimization (x's and open circles) and minimization by convex optimization with the Lindbladian model subject to positive semi-definite constraints (open triangles and squares).We re-use the earlier model parameters use evolve for a time δt/cz ∼ 1e−3.Bottom panel: Normalized reconstruction errors for first and second order finite difference derivative approximations.Here the evolution time is varied and the Hamiltonian (open symbols) and total (lines) estimator errors are given.The least squares (lst) and convex minimixation (cvx) Hamiltonians agree at small times as illustrated in the log-linear inset.