Quantum process tomography of unitary maps from time-delayed measurements

Quantum process tomography conventionally uses a multitude of initial quantum states and then performs state tomography on the process output. Here we propose and study an alternative approach which requires only a single (or few) known initial states together with time-delayed measurements for reconstructing the unitary map and corresponding Hamiltonian of the time dynamics. The overarching mathematical framework and feasibility guarantee of our method is provided by the Takens embedding theorem. We explain in detail how the reconstruction of a single-qubit Hamiltonian works in this setting and provide numerical methods and experiments for general few-qubit and lattice systems with local interactions. In particular, the method allows to ﬁnd the Hamiltonian of a two qubit system by observing only one of the qubits.


I. INTRODUCTION
System identification refers to the estimation of the dynamics of a system from measurements of its characteristics.Its quantum analogue, quantum process tomography (QPT), is essential for the realization and testing of quantum devices [1][2][3], as a benchmarking tool for quantum algorithms [4,5], and in general for understanding the inner workings of a quantum system [6,7].However, textbook algorithms for QPT [8] might be difficult to realize in laboratory settings due to the required preparation of many initial states and observation of the complete target system.
In this work we focus on the unitary time evolution of a closed quantum system.We propose and investigate an approach based on measurements with different time delays, which should be easily realizable in lab experiments.Our main contribution are algorithms and numerical procedures for identifying the corresponding time evolution operator, and thus indirectly the quantum Hamiltonian.Intriguingly, we can identify the operator for the entire system even if the measurements are restricted to a subsystem (using two known initial quantum states).The Takens embedding theorem, discussed in Sect.III, provides an overarching mathematical foundation and feasibility guarantee for our approach.Concretely, the mathematical framework starts with a manifold M, which we take to be the special unitary group of time evolution operators.Given an initial quantum state, a time step matrix U ∈ M then determines the measurement averages at a sequence of time points, see Fig. 1.The Takens theorem states that the map from M to this measurement vector is actually a smooth embedding (under cer- t/∆t ψ(t)| M |ψ(t) 1: Schematic illustration of the map from a unitary time step operator U to a vector of measurement averages at different time delays.tain conditions), which then allows us to reconstruct U .
Related to the present work, a recent approach for quantum process tomography using tensor networks as a parametrization of the quantum channel was able to achieve high accuracies for systems of up to 10 qubits [9].Furthermore, in [10] the authors derive a lower bound in the number of POVMs required to fully characterize a unitary or near-unitary map.Regarding the related task of quantum state tomography, several methods, some based on machine learning techniques, have been devel-oped [11][12][13][14].However, except for [14], these works do not explicitly take advantage of the time trajectory of the quantum system as envisioned here.Gate Set Tomography [15] does not need pre-calibrated quantum states, but relies on very long time series.
Our work is closely related to earlier research on the identification of quantum Hamiltonians from time series [16,17].The authors obtain all degrees of freedom of the Hamiltonian with known structure by solving a system of equations, mostly in the presence of noise.More recently, neural networks have been used for identification, even from experimental data [18][19][20].The minimal number of observables needed to identify the Hamiltonian is not addressed, and we provide it here with the link to Takens theorem.In our work, the choice of initial state is not as important, as we will demonstrate.
Takens theorem has been fundamental to several system identification methods for general dynamical systems already [21], recently also involving machine learning [22][23][24][25], with a focus on PDE [26,27] or special structure such as (classical) Hamiltonian dynamics [28,29].The identification of a unitary map from measurement data has been discussed by Koopman and von Neumann [30,31], work that has been revived in a data-driven context in the last twenty years and extended to dissipative systems [32][33][34].The connection to quantum systems and their special structure and challenges has not yet been addressed in the work cited above, however.

II. PHYSICAL MODEL
We assume throughout that the quantum Hamiltonian H is time-independent and will investigate two settings: (i) a few-qubit system and H a general dense Hermitian matrix, and (ii) a (tight binding) Ising-type model on a two-dimensional lattice Λ, as widely studied in condensed matter physics.For concreteness, the physical system for case (ii) consists a local spin degree of freedom at each lattice site.The Hamiltonian is then defined as where J j, ∈ R and h j ∈ R 3 are local parameters defining the interaction strength and the external field, respectively.σ α j for α ∈ {x, y, z} is the α-th Pauli matrix acting on site j ∈ Λ, and σ j = (σ x j , σ y j , σ z j ) the corresponding Pauli vector.The first sum in (1) runs over nearest neighbors on the lattice.We take Λ to contain a finite number n of sites, and assume periodic boundary conditions.The site-dependent parameters allow to simulate disorder.Fig. 2 illustrates the interaction and external field terms of H on a two-dimensional lattice.The precise form of H is not important for the reconstruction as long as the time dynamics it generates can be well approximated by a quantum circuit, for example via Trotterization.We assume that the overall structure of H is known, and our task is to determine the numerical values of the parameters.
FIG. 2: Visualization of the quantum Hamiltonian in Eq. ( 1) on a two-dimensional lattice, with J the interaction strength and h the external field.
For later reference, we state the unitary time evolution operator at time t ∈ R based on the Schrödinger equation (in units of = 1): ( The solution to the Schrödinger equation for an initial wavefunction (statevector) ψ ∈ C N is then ψ(t) = U (t)ψ.

III. TAKENS EMBEDDING FRAMEWORK
The theorems of Takens and Ruelle [35,36], based on the embedding theorems of Whitney [37], are the theoretical foundations of time-delay embedding we use here.
The general idea we employ in this paper is to "embed" the manifold of unitary matrices (describing the dynamics of a quantum system) into the space of measurement trajectories.We first state the mathematical theorem, and then discuss the specialization for quantum time evolution.Let k ≥ d ∈ N, and M ⊂ R k be a d-dimensional, compact, smooth, connected, oriented manifold with Riemannian metric g induced by the embedding in k-dimensional Euclidean space.Note that this setting is sufficient for our presentation, but is more restrictive than allowed by the results cited below.
Together with the results from Packard et al. [38] and Aeyels [39], the definitions and theorems of Takens [36] describe the concept of observability of state spaces of nonlinear dynamical systems.A dynamical system is defined through its state space (here, the manifold M) and a diffeomorphism φ : M → M. Theorem 1. Generic delay embeddings For pairs (φ, y), φ : M → M a smooth diffeomorphism and y : M → R a smooth function, it is a generic property that the map Φ (φ,y) : M → R 2d+1 , defined by is an embedding of M; here, "smooth" means at least C 2 .
Genericity in this context is defined as "an open and dense set of pairs (φ, y)" in the C 2 function space.Open and dense sets can have measure zero, so Sauer et al. [40] later refined this result significantly by introducing the concept of prevalence (a "probability one" analog in infinite dimensional spaces).See [41] for similar results with stochastic systems.
Let N denote the quantum Hilbert space dimension.In our context, M is the special unitary group SU(N ) ⊂ C N ×N when identifying C R 2 .In particular, M (with underlying field R) has dimension d = N 2 −1.In physical terms, the elements of SU(N ) are the unitary time evolution matrices in Eq. ( 2).We assume that the Hamiltonian H is traceless, since adding multiples of the identity to H leads to a global phase factor in the time evolution, which is unobservable in subsequent measurements as envisioned here.We fix a time step ∆t, which we may (without loss of generality) absorb into H, and set U = e −iH in the following.Now define the diffeomorphism φ via a scaling of H by a factor γ > 0, γ = 1, such that Regarding y, fix a randomly chosen initial quantum state ψ ∈ C N and an observable (Hermitian matrix) M .Now let y compute the corresponding expectation value: With these definitions, the output of the map Φ (φ,y) in Eq. (3) becomes physically corresponding to measurements at time points t q = γ q for q = 0, . . ., 2d.
The main point here is the prospect to identify the time step matrix U , and thus indirectly the Hamiltonian, based on a single measurement time trajectory, under the assumption that the initial quantum state is known.As caveat, the relation U = e −iH∆t determines the eigenvalues of H only up to multiples of 2π/∆t.Moreover, the map Φ (φ,y) might not be one-to-one, in the sense that two different unitary matrices give rise to the same measurement trajectory.We discuss such a case in more detail in the following.Such issues can be avoided in practice by additional assumptions on the structure of H.

IV. NUMERICAL METHODS
Throughout this work, we assume that it is feasible to reliably prepare a single (or when indicated several) known initial state(s) ψ ∈ C N .To demonstrate the general applicability of our methods, ψ is chosen at random in the following algorithms and numerical simulations.Specifically, the entries of ψ before normalization are independent and identically distributed (i.i.d.) random numbers sampled from the standard complex normal distribution.
A. Reconstruction algorithm for a single-qubit system We now describe how to reconstruct a single-qubit Hamiltonian from a time series of measurements.The Bloch sphere picture [8] provides a geometric perspective and insight into single-qubit quantum states and operations.The Bloch vector r ∈ R 3 associated with ψ ∈ C 2 is defined via the relation For pure states as considered here, r is a unit vector.We can parametrize any single-qubit Hamiltonian (up to multiples of the identity map, which are irrelevant here) as with h ∈ R 3 .The corresponding time evolution operator is the rotation when representing h = ω v/2 by a unit vector v ∈ R 3 and ω ∈ R. On the Bloch sphere, this operator corresponds to a classical rotation about v, as illustrated in Fig. 3a and described by Rodrigues' rotation formula (θ = ωt): 3) is a rotation matrix (parametrized by θ and v) applied to r.In the above context of Takens embedding with U = U (∆t), we may equivalently work with U = U ω∆t, v and matrix powers thereof.
The time trajectory effected by the rotation applied to r results in a circle embedded within the Bloch sphere.Knowing the circle would allow to determine v, and the dynamics on the circle to determine ω.By construction we also know one point on the circle already, namely r (the initial condition), but we only have access to the expectation value of an observable M for t > 0. In the following, we denote the "measurement direction" by m ∈ R 3 , i.e., M is parametrized as M = m • σ.Algorithm 1 facilitates a recovery of ω and v using the projections of the time trajectory on m.The main idea is to first reconstruct ω based on the time dependence, and then the unit vector v.
As caveat, the solution to this problem is not unique, due to the four possible signs of the coefficients α 2 and α 3 in the algorithm.Fig. 3a visualizes how two rotations can generate the same projection onto the z-axis, with h resulting from (α 2 , α 3 ) → −(α 2 , α 3 ).The two equations for α 2 and α 3 are shown in Fig. 3b.(In the example, α 1 = −0.2051,and hence the radius of the norm constraint circle is close to 1.) In order to decide between these distinct solutions, one could perform one further measurement in a different basis, or may use a-priori knowledge about the Hamiltonian.
Algorithm 1 Reconstruction of a single-qubit Hamiltonian (measurement direction m ∈ R 3 with m ∦ r) Input: Measurement averages yq = m • (U ωtq , v r) with tq = ∆t γ q for q = 0, . . ., 2d Output: Hamiltonian parameters ω and v. 1: Find ω based on the dependency on ωt in Eq. (10): 4: Represent v with respect to the orthonormal basis with to-be determined coefficients α2, α3 ∈ R. Using that u1, u2 and u3 are eigenvectors of Together with the normalization condition Decide between the four possible signs of α2, α3 using apriori information of H or one additional measurement in a different basis.
We remark that the nonlinear optimization in the first step of the algorithm might get trapped in a local minimum, which could be resolved by restarting the optimization.On the other hand, when neglecting uncertainties associated with the measurements, a correct solution is indicated by a zero residual both in steps 1 and 3. We assume that a range of realistic frequencies ω is known beforehand, such that the Nyquist condition (given the non-uniform sampling points t q ) holds [42].

B. Relaxation method
As straightforward approach for determining a unitary time step matrix U which matches the measurement averages, one could start from the natural parametrization U = e −iH (with the time step already absorbed into H), and then optimize the Hermitian matrix H directly.However, this method encounters the difficulty of a complicated optimization landscape with many local minima, in particular when using gradient descent-based approaches.
Here we discuss an alternative numerical approach tailored towards generic (dense) few-qubit Hamiltonians.As discussed in Sect.IV A, for single qubits a transformation ψ = U ψ by a unitary matrix U ∈ SU(2) can equivalently be described by a spatial rotation in three dimensions on the Bloch sphere: r = Ur, with U ∈ SO(3) given by Rodrigues' formula (10).An equivalent mapping from U to U is via for α, β = 1, 2, 3, where we have used the relation (7) together with the orthogonality relation of the Pauli matrices: 1 2 tr[σ α σ β ] = δ αβ .For our purposes, the Bloch picture has the advantage that measurement averages depend linearly on the Bloch vector, and hence on U, e.g., ψ |σ z |ψ = e z • r = e z • (Ur).In practice, we first optimize the entries of U to reproduce the measurement data (under the constraint that U is an orthogonal matrix), and afterwards find U related to U via Eq.(11).
Generalization to a larger number of qubits is feasible via tensor products of Pauli and identity matrices (in other words, Pauli strings).The analogue of (11) for n qubits, with where α and β are now index tuples from the set {0, 1, 2, 3} ⊗n \{(0, . . ., 0)}, using the convention that σ 0 is the 2 × 2 identity matrix.We exclude the tuple (0, . . ., 0) here since it conveys no additional information: the identity matrix is always mapped to itself by unitary conjugation.
Specifying an element of SO(4 n − 1) involves more free parameters than for a unitary matrix from SU(2 n ) in case n ≥ 2, thus the optimization might find an orthogonal U which reproduces the measurement data, but does not originate from a U ∈ SU(2 n ) via (12).We circumvent this representability issue as follows, focusing on the case of two qubits here.Every U ∈ SU(4) can be decomposed as [43,44] with the "entanglement" gate ) and single-qubit unitaries u a , u b , u c , u d ∈ SU (2).The single qubit gates can be handled as before, i.e., represented by orthogonal rotation matrices u a , u b , u c , u d ∈ SO(3).We find the Bloch representation of R via (12) (cf.[45,46]); the parameters θ 1 , θ 2 , θ 3 appear in the matrix entries solely as cos(θ j ) and sin(θ j ) for j = 1, 2, 3.In summary, we express the decomposition (13) in terms of to-be found orthogonal matrices in the Bloch picture.
After switching to the described Bloch representation, we "relax" the condition that an involved (real) matrix U is orthogonal, by admitting any real matrix, but adding the term to the overall cost function, where • F denotes the Frobenius norm.As advantage, we bypass a parametrization of U to enforce strict orthogonality and avoid local minima in the optimization.
In our case, the U q 's are separate matrices, which are set in relation via an additional cost function term For the case of two-qubit Hamiltonians, we additionally substitute two-dimensional vectors (c j , s j ) for (cos(θ j ), sin(θ j )), again to avoid local minima.The condition c 2 j + s 2 j = 1 translates to another penalty term in the overall cost function: C. Partial (subsystem) measurements An intriguing possibility of the time-delay measurements is the identification of the overall Hamiltonian based on measurements restricted to a subsystem.This scenario will actually require the preparation of more than one exactly known initial state.As minimal example, consider a two-qubit system, the choice between two initial states, being able to select a measurement basis via the gate C, and time-delayed measurements on one of the qubits while the other qubit is inaccessible, as illustrated in Fig. 4. For the numerical experiment shown in Fig. 8 below, we use X-, Y -and Z-basis measurements on the top qubit and minimize the mean squared error between the observed measurement averages and the prediction based on a general Ansatz for the (traceless) Hamiltonian, i.e., 15 real parameters.

D. Lattice system with local interactions
Finally, we consider quantum systems defined on a lattice, and assume an Ising-type Hamiltonian as in Eq. ( 1) (instead of a generic Hamiltonian) to avoid an exponential growth of the number of parameters.
For the purpose of reconstructing the Hamiltonian parameters, the Schrödinger time evolution can be well approximated via a second order Strang splitting method [47], i.e., the Time-Evolving Block Decimation (TEBD) algorithm [48].Interpreting the resulting layout of singleand two-site unitary operators as forming a quantum circuit, the reconstruction task amounts to a variational circuit optimization to reproduce the reference measurement averages.Specifically, when considering the interaction and local field parts of the Hamiltonian by themselves, the individual terms in each of the sums pairwise commute.Fig. 5 illustrates the corresponding quantum circuit for a one-dimensional lattice; the constructing works for higher-dimensional lattices as well.19), illustrated for a one-dimensional lattice with n = 5 sites.The single qubit gates on the left and right are rotation gates exp(−i∆t hj • σj/2) for the j-th qubit, and the twoqubit gates are exp(−i∆tJj,j+1σ z j σ z j+1 ).The different shades indicate different parameters at each site.The ordering of the two-qubit interaction gates among each other is not relevant since they pairwise commute.

V. NUMERICAL EXPERIMENTS
Here we present numerical experiments for the methods described in Sect.IV.

A. Exact reconstruction of a single-qubit Hamiltonian
With the algorithm described in Sec.IV A we are able to reconstruct a single-qubit Hamiltonian up to numerical precision.Fig. 6 shows the results obtained for the ω and least squares optimizations.Following the Takens framework, we use 2d + 1 time points, where d = 3 is the number of Hamiltonian parameters to be reconstructed.Specifically, we have used the time points t q = 0.3×(1.3)q for q = 0, . . ., 6 here.With this, we are able to find ω, α 1 and κ up to an error of ∼ 10 −15 .Thus 7 measurements, plus a further measurement in a different basis to distinguish between the possible solutions due to symmetry, are sufficient for the reconstruction.

B. Relaxation method
In this subsection we consider a generic single-or twoqubit Hamiltonian H. Specifically, we draw the coefficients of H with respect to the Pauli basis from the standard normal (Gaussian) distribution.
We first focus on a single-qubit 2 × 2 Hamiltonian.Fig. 7a shows the result of directly fitting the vector h ∈ R 3 (see Eq. ( 8)), with the KL divergence between the predicted and actual measurement probabilities as cost function.This approach is compared to the "relaxation" procedure (described in Sect.IV B) in Fig. 7b, using γ = 2 and ∆t = 0.1.The manifold for the Takens embedding has dimension d = 3 for a single qubit, and hence 2d+1 = 7 time points should be used.However, we face the difficulty that the Hamiltonian is not uniquely specified solely by Pauli-Z measurements according to the discussion in Sect.IV A. For this reason, we include Pauli-X measurement data as well, and reduce the number of time points to 4 (to compensate for this additional source of data).For both versions we use the Adam optimizer [49].The direct fitting method might get trapped in a local minimum and hence not be able to find the ground-truth vector h, which leads to a large variation around the median.This issue is ameliorated by the relaxation method.
In Fig. 7c visualizes the loss function and relative errors when applying the relaxation procedure to reconstruct the unitary time step matrix and corresponding Hamiltonian of a two-qubit system (d = 15 real parameters).Specifically, for parameter optimization we express Eq. ( 13) using the Bloch picture, i.e., in terms of orthogonal rotation matrices for the single-qubit unitaries, and likewise using the Bloch picture analogue of the entanglement gate (14), parametrized by two-dimensional vectors (c j , s j ) for (cos(θ j ), sin(θ j )).The condition c 2 j + s 2 j = 1 translates to another penalty term in the overall cost function, see Eq. ( 18).We set ∆t = 0.05, γ = 2 and use 6 time points for this experiment.The reference measurement data at each time point are the expectation values of the nine observables {I ⊗ σ α , σ α ⊗ I, σ α ⊗ σ α } α∈{x,y,z} .We use 54 measurement data points in total (instead of 2d + 1 = 31) since we found that the additional information improves the reconstruction.As last step (after the main optimization), we find the Hamiltonian entries giving rise to the computed U via the Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm.
According to Fig. 7c, the median relative errors of U and the Hamiltonian are below 10 −3 , but there are still cases where the method cannot find the reference solution at all, even tough the loss value is small.An explanation could be that U and H are not uniquely determined.We leave a clarification of this point for future work.

C. Partial (subsystem) measurements
We study the scenario depicted in Fig. 4, namely performing measurements solely on one out of two qubits.The observables are the three Pauli gates here.For the reconstruction to work, we require that two different, precisely characterized initial states can be prepared, which we choose at random for the numerical experiments.
The ground-truth Hamiltonian is similarly constructed at random, by drawing standard normal-distributed coefficients of Pauli strings, which in sum form the Hamiltonian.We use the time step ∆t = 1  5 and γ = 1.15 here, and have heuristically found 12 time-delayed measure-ments (for each measurement basis) as viable to reliably reconstruct the Hamiltonian.The loss function is the mean squared error between the model prediction for the measurement averages and the ground-truth values.To avoid getting trapped in local minima during the numerical optimization, we use the best out of 10 attempts (in terms of the loss function) with different random starting points.The specific observables are the three Pauli gates applied to the first qubit.Fig. 8a shows exemplary measurement time trajectories, and Fig. 8b the reconstruction error and loss function (median and 25% quantiles based on 20 random realizations of the overall setup).One observes that a reliable and precise reconstruction of the Hamiltonian (even up to numerical rounding errors) is possible in principle, when neglecting inaccuracies of the measurement process.
We remark that we have used only a single initial state and less time points for the "relaxation" method (see Fig. 7c), which can explain the seemingly larger errors there.

D. Hamiltonian on a lattice with local interactions
We first consider the case of a Hamiltonian (1) with uniform parameters (not depending on the lattice site), i.e., with J ∈ R and h ∈ R 3 .Thus the task consists of reconstructing 4 real numbers.The ground-truth values for the following experiment are J = 1 and h = (0.5, −0.8, 1.1).Λ is a 3 × 4 lattice with periodic boundary conditions, and the initial state is a single wavefunction with complex random entries (independently normally distributed).We compute the time-evolved quantum state ψ(t) via the KrylovKit Julia package [50], and use the squared entries of ψ(t) as reference Born measurement averages.
For the purpose of reconstructing J and h, we approximate a time step via a variational quantum circuit as shown in Fig. 5, where the single-and two-qubit gates share their to-be optimized parameters J and h.For simplicity, we use three uniform time points ∆t, 2∆t, 3∆t with ∆t = 0.2 (instead of time points t q = γ q ∆t), such that the quantum circuit for realizing a single time step can be reused.
We quantify the deviation between the exact Born measurement probabilities, p(t) = |ψ(t)| 2 , and the probabilities p(t) resulting from the trotterized circuit Ansatz, via the Kullback-Leibler divergence: Fig. 9a visualizes the optimization progress, i.e., the relative errors of J and h, and the loss function (21).The darker curves show medians over 100 trials of random initial states, and the lighter shades 25% quantiles.The dashed horizontal line is the loss function evaluated at the ground truth values of J and h; it is non-zero due to the Strang splitting approximation of a time step.Interestingly, the optimization arrives at an even smaller loss value starting around training epoch 70.We interpret this as an artefact of the Strang splitting approximation -note that around this epoch, the relative error of h slightly increases again.We have used the RMSProp optimizer [51] with a learning rate of 0.005 here.In summary, the reachable relative error is around 0.02, and higher accuracy would likely require a smaller time step to reduce the Strang splitting error.
Next, we consider a Hamiltonian (1) with disorder (and external field solely in x-direction): with i.i.d.random coefficients J j, and h j .For the numerical experiment, J j, is uniformly distributed in the interval [0.8, 1.2], and h j ∼ 0.5 N (0, 1) (normal distribution).
The results are shown in Fig. 9b.As before, we evaluate 100 realizations of the experiment, with a set of coefficients and random initial state drawn independently for each trial.We take the maximum over the relative errors of J j, for all j, to arrive at the relative error reported in Fig. 9b, and likewise for h j .The "reference" loss function is the KL divergence evaluated at the ground truth coefficients.It is non-zero due to Strang splitting errors, and fluctuates due to the different random coefficients at each trial.

VI. CONCLUSIONS AND OUTLOOK
Our work demonstrates the feasibility of using measurements at different time points to reconstruct the time evolution operator.The approach presented in this work has two main benefits: First, it requires only a single (or few) initial state(s) and reduces the number of different kinds of measurements for the reconstruction of the Hamiltonian.From an experimental point of view, this would be helpful when the experimental setup makes it easier to maintain the time evolution of the system than to implement a variety of measurements.Second, as demonstrated in Sect.V C, the method allows for the reconstruction of a Hamiltonian even when only part of the system can be observed.
A natural extension of our work is QPT in the situation of dissipation and noise processes, i.e., a system governed by a Lindblad equation.This would pose the additional challenge of a limited time window for extracting information and a larger number of free parameters.An interesting alternative approach could be a mapping to a unitary evolution with an unobserved environment [52], which would fit into the present framework.
The relaxation method in Sect.IV B can be regarded as tool to smooth the optimization landscape, but an open question is how to guarantee convergence to the correct solution.This becomes particularly relevant for larger systems and the likewise increasing number of free parameters.
A related task for future work is a sensitivity analysis with respect to inaccuracies and noise in the measurement data.A modified version of Takens theorem still applies in this situation [41], but it is not clear how accurate our algorithm can reconstruct H.A first step could be a gradient calculation of the measurement averages with respect to the Hamiltonian parameters.

conditions on α 2 and α 3 coefficientsFIG. 3 :
FIG. 3: (a)The Bloch sphere picture of the time dynamics effected by a Hamiltonian is a classical rotation.For a single observable, the reconstruction of H is not unique: e.g., the two trajectories for H = h • σ and H = h • σ result in the same Z expectation values.(b) The system of equations to be solved for the reconstruction admits four solutions.

|ψ e − iHt CFIG. 4 :
FIG.4: Minimal example for the time-delayed partial (subsystem) measurement scenario.The goal is to reconstruct the (unknown) Hamiltonian H.We assume that one can prepare two different initial states, and can select a measurement basis via the gate C.

FIG. 5 :
FIG.5:One time step ∆t of the quantum dynamics based on Strang splitting into interaction and local field terms, see Eq.(19), illustrated for a one-dimensional lattice with n = 5 sites.The single qubit gates on the left and right are rotation gates exp(−i∆t hj • σj/2) for the j-th qubit, and the twoqubit gates are exp(−i∆tJj,j+1σ z j σ z j+1 ).The different shades indicate different parameters at each site.The ordering of the two-qubit interaction gates among each other is not relevant since they pairwise commute.
FIG. 6: (a) step 1 and (b) step 3 of Algorithm 1 using Z-basis measurements, for the time evolution shown in Fig. 3a.
FIG.7:(a) and (b): Single-qubit Hamiltonian reconstruction based on four Z-and four X-basis measurements, (a) by direct fitting of the Hamiltonian parameters to the measurement data, and (b) using the "relaxation" procedure (Sect.IV B), with the actual Hamiltonian parameters determined from U at the end (blue line at the rightmost section of the plot).Solid lines show the median, and lighter shades are 25% quantiles.(c) "Relaxation procedure" for two-qubit Hamiltonian reconstruction, working with the Bloch picture of the representation(13) to find the time step matrix U, and then determining a two-qubit Hamiltonian giving rise to this U.

FIG. 8 :
FIG. 8: Numerical reconstruction of a generic two-qubit Hamiltonian based on time-delayed measurements of only one of the qubits, cf.Fig. 4. (a) Exemplary measurement time trajectory of the operators σ x ⊗ I, σ y ⊗ I and σ z ⊗ I. Two random initial states and their respective trajectories are used as input to the reconstruction algorithm.(b) Reconstruction error and loss function based on 20 random realizations.

FIG. 9 :
FIG.9: Relative reconstruction error of the Hamiltonian parameters for the 3 × 4 lattice model in (a) Eq. (20) and (b) Eq. (22), and loss function during training, when fitting a quantum circuit representation of a time step (Fig.5) to the exact Born measurement averages.