Hamiltonian tomography for quantum many-body systems with arbitrary couplings

Characterization of qubit couplings in many-body quantum systems is essential for benchmarking quantum computation and simulation. We propose a tomographic measurement scheme to determine all the coupling terms in a general many-body Hamiltonian with arbitrary long-range interactions, provided the energy density of the Hamiltonian remains finite. Different from quantum process tomography, our scheme is fully scalable with the number of qubits as the required rounds of measurements increase only linearly with the number of coupling terms in the Hamiltonian. The scheme makes use of synchronized dynamical decoupling pulses to simplify the many-body dynamics so that the unknown parameters in the Hamiltonian can be retrieved one by one. We simulate the performance of the scheme under the influence of various pulse errors and show that it is robust to typical noise and experimental imperfections.

Characterization of qubit couplings in many-body quantum systems is essential for benchmarking quantum computation and simulation. We propose a tomographic measurement scheme to determine all the coupling terms in a general many-body Hamiltonian with arbitrary long-range interactions, provided the energy density of the Hamiltonian remains finite. Different from quantum process tomography, our scheme is fully scalable with the number of qubits as the required rounds of measurements increase only linearly with the number of coupling terms in the Hamiltonian. The scheme makes use of synchronized dynamical decoupling pulses to simplify the many-body dynamics so that the unknown parameters in the Hamiltonian can be retrieved one by one. We simulate the performance of the scheme under the influence of various pulse errors and show that it is robust to typical noise and experimental imperfections.
Introduction.-Physicists have been striving to understand and harness the power of quantumness since the establishment of the quantum theory. With the flourishing of quantum information science in recent decades [1,2], numerous breakthroughs-both in theory and in experiment-helped to frame a clearer goal: it is the entanglement and the exponentially growing Hilbert space that distinguishes quantum many-body systems from classical systems [3][4][5]. To fully leverage the quantum supremacy, a vital step is to verify and benchmark the quantum device. The standard techniques of quantum state and process tomography [6][7][8][9][10][11], however, are plagued by the same exponential growth of dimensions [12]. A related problem is to directly identify Hamiltonians, the generators of quantum dynamics. They can often be specified by fewer number of parameters that scales polynomially with the system size.
Hamiltonian tomography for generic many-body systems is nevertheless a daunting task. The way to extract information of unknown parameters in a Hamiltonian is by measuring certain features of its generated dynamics. To make this possible, one has to solve the dynamics generated by the Hamiltonian to make a definite connection between its dynamical features and the Hamiltonian parameters. However, for general many-body Hamiltonians, their dynamics are extremely complicated and intractable by numerical simulation as the simulation time increases exponentially with the size of the system. Progress in this direction has mostly be on small systems [13][14][15][16][17] or special many-body systems which are either exactly solvable due to many conserved operators, of limited Hilbert space dimensions amenable to numerical simulation, or short-range interacting systems [18][19][20][21][22][23][24].
In this paper, we propose a scheme to achieve Hamiltonian tomography for general many-body Hamiltonians with arbitrary long-range couplings between the qubits. The key idea is to simplify the dynamics generated by a general many-body Hamiltonian through application of a sequence of dynamical decoupling pulses on individual qubits. Dynamical decoupling (DD) is a powerful technique that uses periodic fast pulses to suppress noise and average out unwanted couplings between the system and the environment [25][26][27][28][29][30][31][32][33][34][35][36][37][38][39][40]. We apply a sequence of synchronized DD pulses on a pair of qubits, which forms a small target system that has coupling with the rest of the qubits in the many-body Hamiltonian, the effective environment. The DD pulses keep the desired couplings within this target system intact while average out its couplings with all the environment qubits. The dynamics under the DD pulses become exactly solvable, from which we can perform a tomographic measurement to determine the coupling parameters within this small target system [13][14][15][16]. We then scan the DD pulses to different pairs of qubits to measure all the other coupling terms in the Hamiltonian. We assume the ability to address individual qubits, which is realistic for many experimental platforms, such as trapped ions [41][42][43], cold atoms [44][45][46], and solid-state qubit systems [47][48][49][50]. Several features make the scheme amenable to experimental implementation. First of all, applying the DD pulse sequence is a standard procedure in many experiments. Post-processing of data is straightforward as it only requires one or two parameter curve fitting. In addition, we demonstrate with explicit numerical simulation that the scheme is robust to various sources of errors in practical implementation, such as the remnant DD coupling error, measurement uncertainties, and different types of pulse errors.
Scheme for Hamiltonian tomography.
-The system we have in mind is the most general Hamiltonian with twobody qubit interactions where J αβ mn characterizes the coupling strength between spins m and n for the α, β components, and b α m represents arXiv:1505.00665v3 [quant-ph] 4 Oct 2015 1. (color online). Schematics for the tomography procedure. (a) To map out the coupling coefficients J αβ ij , a synchronized DD sequence is applied to spins i and j. Both spins will be decoupled from the rest of the system. (b) The XY -8 DD sequence on spins i and j to probe the parameters of the Hamiltonian in Eq. (3). The initial state is for instance prepared to the |00 state for the two spins. (c) To retrieve information about the local fields b α i , XY -8 pulse sequences are applied to the environment spins to decouple spin i from the rest.
the local field on spin m; σ α σ β are the Pauli matrices along the α (β) direction with α, β ∈ (x, y, z). To adopt consistent notations throughout the text, we use m, n to denote a general spin label and i, j to refer to the specific target spins that we are probing with the DD pulses, calling the rest of the spins as environment spins. The terms spin and qubit are used interchangeably. Let the energy unit of the Hamiltonian be J, chosen to be the largest magnitude of all coefficients, so J αβ mn /J and b α m /J are bounded between −1 and 1. In order to map out the coupling coefficient J αβ ij for the target spins, we propose to decouple these two spins from the environment spins by a synchronized DD pulse sequence. A synchronized XY -4 sequence applied to both spins will average out their interactions with other spins while preserving the two-spin coherence (see Fig. 1(a-b) for the schematic and the pulse sequence). Basically, only those interactions that commute with the DD sequence will survive. More rigorously, the evolution operator in one period is where U 0 = e −iHτ , τ is the time interval between two consecutive pulses, and B, the bath, includes all terms of the Hamiltonian that only acts on environment spins. See Supplemental Material for the detailed derivation [51]. To bound the error term to O(J 2 τ 2 ), we assume n J αβ in = O(J), i.e., the interaction strength decays rapidly with spin separation distance so that the energy density of the Hamiltonian is bounded by a constant. This condition is satisfied for any finite systems as in the experiment with arbitrary interactions. In the thermodynamic limit, it is also a reasonable assumption for any physical systems whose energy is extensive. It may also be related to the generalized Lieb-Robinson bound for systems with long-range interactions [52][53][54][55]. The XY -8 pulse sequence, which is the concatenation of XY -4 sequence with its time-reversal, eliminates the error term to the third order O(J 3 τ 3 ). Fig. 1(b) shows the XY -8 DD pulse sequence. Hence, in the Hilbert subspace of the target spins, the effective Hamiltonian is where we use c 1 ≡ J xx ij , c 2 ≡ J yy ij , c 3 ≡ J zz ij to simplify the notation. The effective two-spin unitary evolution after N c cycles of XY -8 sequence is where T = 8N c τ is the total time. From the above expression, one may notice that the Hamiltonian parameters can be retrieved by preparing a particular initial state and measuring its time-evolved output probability in a given basis. In particular, we have where |+ = 1 √ 2 (|0 + |1 ) and |I = 1 √ 2 (|0 + i|1 ) are the rotated basis. The coupling strengths c 1 , c 2 and c 3 can be extracted from the oscillation frequencies of these three sets of measurements at various time points. These particular sets are not the only suite to extract those parameters. They are chosen for the convenience in fitting and in state preparation. Only product states of the two target spins, disentangled from the rest, are required. We also remark that the error incurred is O(N c J 3 τ 3 ), so one needs Jτ 1 for a robust decoupling scheme. In a similar fashion, one can retrieve all other coupling coefficients. Let us denote the synchronized XY -8 DD pulse sequence as X i X j -Y i Y j -8 to show explicitly the particular pulses on specific spins. Replacing the sequence with pulses, we will be able to extract the coefficients J xy ij , J yz ij and J zx ij (J yx ij , J xz ij and J zy ij ), respectively.
By scanning the DD pulses to different target pairs, the above procedure recovers all the coupling coefficients J αβ mn . The retrieval of local field coefficients follows a similar approach. We now need to decouple the particular spin i from the rest without contaminating its own spin term b α i σ α i . Shining a XY -8 DD sequence on spin i removes all information about b α i too. Instead, one could address all the environment spins with XY -8 pulses, and decouple them from spin i (alternative schemes are discussed in the Supplemental Material [51]). This scheme will be very robust to pulse errors, since no laser pulses are directly applied to the target spin [ Fig. 1(c)]. The effective single-spin Hamiltonian is thus i with a unitary evolution U 1-spin = e −iH1-spinT executing a spin rotation on the Bloch sphere. Again, by preparing a particular state and measuring its time evolution, we get P |0 →|0 = 1 + (b z i /b) 2 − 1 sin 2 (bT ) and P |+ →|+ = 1 + is the magnitude of the Bloch vector. These two sets of measurements will determine b x i , b y i and b z i up to a sign. The correct signs from the remaining discrete set can be picked out by measuring P |+ →|0 and P |I →|0 at a single time point [51].
The complete scheme applies to any generic Hamiltonian with interacting qubits. In the most general case, one needs to determine 9N (N − 1)/2 + 3N coefficients. However, in many physical systems, the particular form of the interaction is known and/or the interaction often decays fast enough that one can significantly reduce the number of measurements required. In particular, if J αβ mn can be truncated at some spin separation distance in the case of short-range interactions, the number of measurements will be linear with the system size N . In the following, we numerically simulate the experimental procedure for the most general Hamiltonian, taking into account various sources of errors, including the remnant DD coupling error, measurement uncertainties, and different forms of pulse errors.
Numerical simulation.-We consider the general Hamiltonian given in Eq. (1) with coefficients J αβ ij /J and b α i /J randomly drawn from −1 to 1. In our finite-system simulation, we ignore the decay of J αβ ij with distance, so the system may include unphysically long-range interactions and could simulate Hamiltonians in any dimensions. To retrieve J αα ij , for example, we start with a product state of all spins, and perform time evolution using the entire Hamiltonian from Eq. (1), interspersed with the XY -8 DD pulses on target spins i and j. We would like to emphasize that specific state initialization for the environment spins is not required as long as they are disentangled from the target pair of qubits at the beginning. After N c cycles of the DD sequence, the environment spins are traced out and measurements are made on spins i and j. In the simulation, we do not assume the pure unitary evolution U 2-spin as the remnant coupling to the environment spins may entangle the two spins with the rest. However, any undesired couplings are suppressed to the order of O(J 3 τ 3 ) and we do observe that the two-spin density matrix remains mostly pure (∼ 99.9%) for our chosen parameters.
As the tomography procedure involves measuring the output probability of a certain state, each time point will be measured N m times, which gives an estimate of the probability p m in this state. The measurement uncertainty (standard deviation) will be p m (1 − p m )/N m following the binomial distribution. As discussed above, to map out c 1 , c 2 and c 3 , one needs to measure P |+I →|00 , P |+I →|10 and P |0I →|++ for the target spins at various time points and extract the corresponding oscillation frequencies. Suppose N t different time points are measured for each set. The oscillation frequencies can be found either by Fourier transform or by curve fitting. In general, if data show numerous oscillation periods, Fourier transform will be more robust and reliable [14][15][16]. In our case, however, the long time observations will be undermined by the remnant coupling to the environment spins and possible pulse error accumulation. Simple curving fitting with fewer oscillation periods, therefore, appears to be a better solution. In Fig. 2(a-c), we fit the data with the method of least squares with τ J = 0.01, N m = 100, N t = 50 for spins i = 7 and j = 9 in a N = 12 spin system. The blue solid lines are the best-fit lines, and the red dashed lines are the theoretical lines using the true coupling coefficients. The longest time period requires 800 pulses, which is well within the current experimental technology without significant pulse error accumulation. Table I compares the true values and the estimated ones of J αα 79 . Uncertainties in the estimation stem from the curve fitting due to measurement uncertainties. Corresponding results for b α 6 of spin 6 are shown in Fig. 2(d-e) and Table I. All estimated parameters are accurate within a few percent.
To simulate real experiments, one also needs to include possible pulse errors. One possible source of errors is the finite duration of each control pulse, which limits the minimum cycle time. This is typically not the dominant source of errors and can often be well-controlled [31,35,[56][57][58]. In most experiments, the major cause of errors is the deviation between the control pulses and the ideal X or Y pulses. These can either arise from the amplitude error where the rotation angle differs from the ideal π-pulse or the rotation error where the rotation axis deviates from the x or y axis. In typical experiments, individual pulse errors may be controlled within a percent level. In our simulation, we consider three different forms of pulse errors: Systematic Amplitude pulse Error (SAE), Random Amplitude pulse Error (RAE) and Random Rotation axis Error (RRE). See the caption of Table I for the specific forms of the errors. Moderate systematic errors can be self-compensated by the XY -8 DD sequence. Numerically, we found that 5% of SAE has negligible effect on the parameter estimation. In addition, we also simulated the cases where each pulse experiences a 1% RAE or RRE. Results are summarized in Table I. The average deviation from the true parameters are within 5%. Here, we would like to point out a few features of our scheme that make it inherently robust to errors. First of all, the estimation of the coupling strength J αβ mn only entails frequency estimation, which could endure large deviations of a few measurement points. In addition, the single-parameter curve fitting scheme not only makes the estimation robust but is also more convenient for experiment. Moreover, the retrieval of local fields b α m is remarkably tolerant to pulse errors. Since no pulse is directly applied to the target spin, any pulse errors on the environment spins will only be propagated via the remnant DD coupling error, which is suppressed to the order of O(J 3 τ 3 ). We have numerically tested that a 10% pulse error of any kind would have negligible effects on the estimation of b α m . Alternative schemes to extract the local fields are detailed and discussed in the Supplemental Material [51]. They are less tolerant to pulse errors, but may be easier to implement in some experimental setups.
Discussion and outlook.-We have thus numerically demonstrated that the proposed scheme is robust to various sources of errors present in real experiments. The Random Rotation axis Error; AD: Average Deviation from true values. The last digit in bracket for each number quantifies the estimation error bar due to measurement uncertainties, which is generated by the bootstrapping method. The percentage values in the brackets denote the amount of errors introduced in each pulse. The errors are in the form of: SAE, e i π 2 (1+ )σ ν ; RAE, e i π 2 (1+δ)σ ν ; RRE, e i π 2 (σ ν +ασ x +βσ y +γσ z ) ; where = 5%, δ is randomly chosen from (−1%, 1%), (α, β, γ) is a vector with a random direction but fixed magnitude at 1%, and ν = x, y for the X and Y pulses respectively.

True
Estimated measurement uncertainties can be lowered by increasing N m and the pulse errors can be reduced by limiting the maximum number of pulses needed. The optimal strategy involves a delicate balance between experimental sophistication and error control. For example, by fixing τ J and the total number of measurements for each set, N m × N t , one could devise an optimal estimation procedure. In addition, it is also possible to eliminate the remnant DD coupling error to a higher order with more elaborate pulse sequences such as the concatenated DD sequence [30,31] and reduce pulse errors by designing composite pulses or self-correcting sequences [34,35,59,60]. The scheme can also be extended straightforwardly to qudit systems of higher spins or to bosonic or fermionic systems.
In conclusion, we have proposed a general scheme to achieve full Hamiltonian tomography for generic interacting qubit systems with arbitrary long-range couplings. The required number of measurements scales linearly with the number of terms in the Hamiltonian, and the scheme is robust to typical experimental errors or imperfections.
We would like to thank Z.-X. Gong for discussions. This work is supported by the IARPA MUSIQC program, the ARO and the AFOSR MURI program. In this Supplemental Material, we provide more details on the dynamical decoupling scheme and error estimation, the retrieval of local field parameters and alternative schemes, and include further results taking into account of pulse errors.
where the remnant coupling noise term is In the error term C, the biggest contribution comes from terms like [B, σ α m B α m ]. Our aim is to show that the error does not scale with the system size N , i.e., C = O(J 2 τ 2 ). Let us consider one such term and write it out explicitly (suppressing the α, β summation): Since J αβ mn and J xγ im are rapidly decaying functions of the separation distance, for a fixed site i, m<n J αβ mn J xγ im = O(J 2 ). Note that this differs from the scaling of the Hamiltonian, H ∼ m<n J αβ mn = O(N J). All the other terms in C are either smaller or contribute to the same order as the above term. Therefore, we have C = O(J 2 τ 2 ). To be able to neglect the error terms, one needs to fulfill the condition Jτ 1. The above discussion is pertinent to the XY -4 pulse sequence. We can cancel the second order contribution by using the XY -8 pulse sequence as U 2 = U 1 U R 1 , where U R 1 is just the time-reversed sequence of U 1 . It can be readily seen that the error terms C and C R will cancel each other to the second order O(τ 2 ), since C R contains the same terms as in C only with the role of H 2 and H 3 interchanged. Therefore, the remnant coupling error of the XY -8 pulse sequence is O(J 3 τ 3 ) as discussed in the main text.

LOCAL FIELD RETRIEVAL
In the main text, we proposed a scheme to retrieve the local fields b α i by shining the XY -8 pulse sequences on all the environment spins. Here, we provide more details and outline alternative schemes that may in some experimental setups be easier to implement. By decoupling the environment spins with spin i as illustrated in Fig. 1(c) of the main text, we have the effective single-spin Hamiltonian The time evolution operator is is the magnitude of the Bloch vector. By measuring at various time points, we could determine |b x i |, |b y i |, |b z i |. To pin down the correct signs, one can supplement the above two sets of measurements with another two measurement points: Only one time point is needed to determine the signs. For example, one could take measurements at bT = π/4 and use P |+ →|0 and P |I →|0 to pick out the correct signs. The above procedure requires applying the DD sequences to all spins other than the target spin. In some experimental setting, it may be easier to apply a global DD sequence to all spins and add another individually addressed beam on spin i to cancel the DD sequence on that single spin. See Fig. 3 for illustration. For instance, one could apply synchronized X All Y All -8 global pulses and in addition X i Y i -8 focused pulses on spin i. In this way, spin i effectively 3. Alternative scheme to map out the local fields b α i . A global pulse imposes the XY -8 pulse sequence on all spins and a focused pulse is in addition applied to spin i to cancel the DD sequence on that single spin.
experiences no pulses at all time. The effective Hamiltonian again reduces to the same H 1-spin as above. However, this scheme is not very robust to pulse errors. Any deviation from the ideal pulse will be doubled on spin i and accumulate. The pulse error will affect the single-spin coherence and obscure b α i too. We have tested it numerically that the pulse errors have to be controlled within 0.5% for the scheme to be feasible. So it can be used in some setups where pulse errors are not an issue or the total number of pulses can be reduced. One may also use this scheme and modify the sequence by designing composite pulses or self-correcting sequences to reduce pulse errors.

PULSE ERRORS
In the main text, we discussed different types of pulse errors. In our numerical simulation, we considered Systematic Amplitude pulse Error (SAE), Random Amplitude pulse Error (RAE) and Random Rotation axis Error (RRE). The fitting curves in Fig. 2 of the main text do not take into account of pulse errors. Here, we include the figures (Fig. 4) for the case with a 1% RRE. We can see, for example in Fig. 4(a), that the frequency estimation is still very accurate while some measurement points may have a notable mismatch. We may also notice that the estimation of b α i is exceptionally robust to pulse errors since no pulse is applied to spin i in the scheme. Other pulse errors have similar effects on the estimation of parameters. The RRE is of the form e i π 2 (σ ν +ασ x +βσ y +γσ z ) where (α, β, γ) is a vector with a random direction but fixed magnitude at 1%. (a)-(c) are used to retrieve J xx 79 , J yy 79 and J zz 79 between spins 7 and 9. (d) and (e) are used to extract b x 6 , b y 6 and b z 6 for spin 6. The blue solid lines are the best-fit lines with the simulated experimental data, and red dashed lines are the theoretical ones generated by the true Hamiltonian parameters. Other parameters are the same as in Fig. 2 of the main text.