Indirect Quantum Tomography of Quadratic Hamiltonians

A number of many-body problems can be formulated using Hamiltonians that are quadratic in the creation and annihilation operators. Here, we show how such quadratic Hamiltonians can be efficiently estimated indirectly, employing very few resources. We find that almost all properties of the Hamiltonian are determined by its surface, and that these properties can be measured even if the system can only be initialised to a mixed state. Therefore our method can be applied to various physical models, with important examples including coupled nano-mechanical oscillators, hopping fermions in optical lattices, and transverse Ising chains.


I. INTRODUCTION
There has been considerable interest in the problem of Hamiltonian identification through indirect probing , thereby developing various quantum mechanical versions of classical system tomography or classical 'inverse scattering' problems [1].For certain types of interactions, it was found [2][3][4][5][6][7] that only few resources are required to obtain an accurate model of the system.Indirect Hamiltonian estimation is therefore an interesting problem for both pragmatic purposes and fundamental insights.We are interested in the following questions: How can we obtain precise information about a Hamiltonian under restricted access?What can we learn about the 'inside' of a large system by only looking at a subsystem of it?Under which conditions is such indirect probing possible?
Recent studies have focused on this problem for cases of chains and networks of spin-1/2 particles.The common question addressed can be formulated as follows: Can we estimate all parameters, such as coupling strengths and local fields, by accessing only one or a few spins?It should be emphasised that even direct Hamiltonian estimation, or more generally, process tomography, is hard, because the required number of measurements and the complexity of the post-processing both scale exponentially with the system size.However, in realistic situations, we usually have a priori knowledge based on the underlying physics.It has been shown that such knowledge can be used to develop compressed sensing protocols [8,9], which greatly reduce the complexity of process tomography.Various works on indirect Hamiltonian estimation have relied on similar assumptions; namely, that the dynamics is restricted to a subspace of polynomial dimension [2][3][4].In [2], the efficiency of the estimation in terms of the required time and the number of measurements is discussed.An interesting example that does not rely on a subspace was analysed by Di Franco et al. [5].We will see here that this is a special case of the generic estimation of quadratic Hamiltonians, which can be estimated efficiently due to a simple description of their dynamics in the Heisenberg picture.Di Franco also found that the estimation is quite robust against noise [5].In [3] the 1D methods were generalised to arbitrary graphs, and the possible elimination of degeneracies was discussed.Also, Wiesniak and Markiewicz [4] went beyond the simplest subspace in order to study quasi-1D systems.Table I summarises the results obtained so far in terms of the settings and assumptions considered.However, the analysis of physically important cases, such as the transverse Ising model and the XY model with a magnetic field, have remained open.Solutions to both cases will be presented in this work.
Our main goal in this paper is to develop a method to perform indirect quantum tomography for many-body systems of identical particles.Even though the method is analogous to the spin case, the Hamiltonians considered here have a higher number of parameters, and it is surprising that they can still be estimated in a similar fashion.The class of Hamiltonians we study here is those of quadratic form in bosonic or fermionic operators.There has been tremendous Figure 1: Indirect classical system tomography of a quadratic Hamiltonian.In this example, the spring constants ki and masses mi of a chain of coupled harmonic oscillators can be determined by monitoring the dynamics of a single particle at the end (in red).See [1] for details.

Interaction type
Needs preparation Geometry Obtain Reference progress in experiments of quantum random walks [10], optical lattices [11], coupled cavities [12], nano-mechanical oscillators [13], etc., which can be modelled by such quadratic Hamiltonians.Thus, the indirect estimation scheme we present here will be of use in reducing the necessary resources for modelling such systems.For the case of bosons, it is the most direct translation of the work by Gladwell [1] to the quantum case.Gladwell studied how the spring constants and masses of coupled classical harmonic oscillator chains can be estimated by looking at the movement of only one particle (see also Fig. 1).Furthermore, our estimation protocol gives a natural generalisation of the problems on spin chains, since quadratic Hamiltonians of fermions also describe a certain class of spin systems, such as the transverse Ising model.
The main method of indirect estimation is summarised as follows.First, the system is initialised to an arbitrary but fixed state.This can be, for example, even a thermal state, and can occur on slow time-scales via relaxation.Then, some simple single-particle properties are initialised and measured at a later time.Finally, the accumulated data is Fourier transformed, and the parameters are extracted through a set of linear equations.This simple method is outlined in Fig. 2 for the 1D case.The procedure is analogous to an 'inverse scattering' problem because the perturbation introduced in one edge of the sample (e.g., rotation of the first qubit) propagates through the sample, 'scattering' with the inner structure of the Hamiltonian, and then this information encodes the structure of the system.While it is obvious that this procedure provides some information on the system, the surprising result here is that all information can be uniquely identified.Since we obtain information on anisotropies as well, not only the coupling strengths but also the type of interaction can be determined by our method.
The paper is structured as follows.First, we introduce the necessary notation for quadratic Hamiltonians and some techniques for their diagonalisations.Since these are well established methods, we will only discuss them as a 'recipe' for the estimation procedure.We then discuss the simplest case of estimation, namely when the system is a chain of hopping particles, and later generalise it to arbitrary graphs.Finally we discuss how the results apply to 1D chains of spins and conclude.

II. NOTATION AND DIAGONALISATION
The most general quadratic Hamiltonian of N indistinguishable particles is written as where the a † n and a m are creation and annihilation operators and A, B are matrices describing the parameters we would like to estimate.For H to be Hermitian we must have A = A † and B T = −ǫB, where we introduced the parameter ǫ = 1 for fermions and ǫ = −1 for bosons.We will mostly follow the notation of [14], though we shall write all vectors in Dirac notation.At first, we put all Hamiltonian parameters into the Hermitian 2N × 2N matrix and introduce the column vector operator , so that Eq. ( 1) can be expressed up to a constant as Throughout this paper, we make two technical assumptions: first, all the phases of M ij are assumed to be known.
Although some phases might be easy to determine, this requires complicated studies of gauge invariance that do not seem to be worthwhile, as in many practical cases all elements are real and positive.Second, for the bosonic case we assume that the matrix M is positive definite.Again, in principle this can be generalised, but this way we avoid technical difficulties of symplectic transformations [14].
As in [3] the efficiency of our method depends on how many entries of M are a priori known to be zero; that is, on knowledge of the coupling graph.If such knowledge is not available, we have to perform measurements on all but one of the qubits.If the graph is known to be highly sparse (for instance a chain) we only need to access a single qubit.But before going into the details of Hamiltonian identification, let us briefly review the diagonalisation of the Hamiltonian Eq. ( 1) and thus the dynamics, introducing some notations.For more detailed descriptions on the diagonalisation procedure, see, e.g., [14].
For quadratic Hamiltonians of the form in Eq. ( 1), there exist quasi-particle creation and annihilation operators b † k and b ℓ , with which the Hamiltonian can be represented by the simple form of non-interacting modes, For this reason, quadratic Hamiltonians are also referred to as 'quasi-free' interactions.We need to know the transformation T that maps the operators a and a † for particles to b and b † for quasi-particles, i.e., β = T α, where β is defined by .
In order to ensure the canonical commutation relations for the operators b k and b † k , T must satisfy T −1 = ηT † η, where The Hamiltonian is now written as It can then be shown that ηM is diagonalised by T as and the inverse of T is given by For bosons, the matrix ηM is not Hermitian and the distinction between right and left eigenvectors is necessary.This gives rise to a few further peculiarities, such as the modified normalisation and completeness relationship (see below).Fortunately, in this work we are solving an inverse problem, and do not have to discuss how to find these vectors and how numerically stable the corresponding algorithms are.It is worth pointing out that the |E k are not representing physical states but just introduced here as a part of solving the Heisenberg equation of motion for the creation and annihilation operators.The completeness relationship is given by and the vectors |E k are chosen to fulfil the normalisation relationship For convenience, let us also introduce vectors |n as the canonical basis vectors: Due to the structure of the matrix M, the upper and lower eigenvectors of ηM are related as where ⊕ is the addition modulo 2N.The dynamics of the original operators α can be found from β n (t) = e −iEmt β n (0) as where we have introduced a sign function s(m, k) through

A. Experimental requirements
Let us first consider a 1D chain of interacting particles, meaning that A and B are tridiagonal.Also assume that we can initialise the chain in a fixed state ρ 0 .This state could, for instance, be a thermal state, but we do not require to have the exact form of ρ 0 : we just have to be able to repeatedly initialise the chain to the same state ρ 0 .As before [2,5], we perform initialisations followed by measurements at the first site.The quantity we need to measure at the first site is a 1 ± a † 1 for different times up to N 2 [2].In order to eradicate the dependence on the initial state ρ 0 , we measure two sequences of a 1 (t) after preparing the first site to give two different initial values, i.e., a 1 (0) = c 1 and a 1 (0) = c 2 .Using a n = a † n , and thus a † 1 (0) = c * i , and subtracting the measurement results, we obtain a quantity that only depends on ∆c ≡ c 1 − c 2 .It is given through Eq. ( 6) by This initialisation can be performed by a von Neumann measurement, or, as long as the reduced density matrix at site one is not maximally mixed, by applying different single qubit rotations (in some experiments von Neumann measurements are hard).As we see in Eq. ( 7), the dependence on the initial state is completely removed.This is thanks to the absence of the interactions between particles: they (almost) do not see each other, so the information on the 'injected' particle can be extracted by subtracting the influence from others.
For the spin chain case the eigenfrequencies are non-degenerate and T −1 1k = 1|E k = 0 (∀k) [2,3].For the present case of quadratic Hamiltonians we were unable to prove this, but could confirm it numerically.Hence, a Fourier ..... ..... analysis provides us with the frequencies E k and the amplitudes ∆c . Summing these amplitudes gives the value of ∆c : where we used the completeness relationship Eq. ( 4).Equation ( 7) still contains mixtures of the coefficients We can separate them by measuring another pair of initialisations c ′ i : as long as ∆c ′ = r∆c (r ∈ R), we can solve the linear equation for |T −1 k1 |.Without loss of generality, we choose by arranging the global phase of each eigenstate |E k .In conclusion, a few random rotations or initialisations of the first qubit, followed by measurements, provide us with E k and 1|E k .
Let us now describe how to obtain the parameters of M from this observed data.We have to distinguish between the generic case where the off-diagonal couplings A n,n+1 and B n,n+1 are distinct, and the special case where they are equal.

B. Different off-diagonal couplings
As we have seen above, what we diagonalised is the 2N × 2N matrix ηM , so it is helpful to regard M as a representation of a graph consisting of 2N nodes (see Fig. 3).Its off-diagonal entries correspond to the coupling strengths between nodes, whereas the diagonal elements represent the intensity of the 'field' at each node.We can then start with a recursive algorithm similar to [1][2][3] by applying M to the localised states at sites 1 and N + 1; Using these equations for 1|ηM |E k and N + 1|ηM |E k we arrive at Note that B 11 = 0 for fermions.Among these parameters, we can immediately obtain A 11 and B 11 from the quantities that are estimated by measuring a 1 (t) as Therefore we can collect all the known terms of the above equations on the left-hand-side as Taking linear combinations, we obtain It may appear as if the right-hand-side contains too many unknowns to solve these equations.However, similar to the original work by Gladwell [1] and the spin chain case, [2], we can use the normalisation of the eigenstates.In this case, it is given by Eq. ( 5).Summing up the mod squares of the above equations, the dependencies on 2|E k and N + 2|E k vanishes and we can infer the absolute value of each coefficient, i.e., .By induction, we then see that all matrix elements of A and B can be obtained, as desired.It is worth pointing out that the scheme works even though the graph in Fig. 3 is not infecting [3].
We will now look at the cases with equal off-diagonal couplings in more detail, because such physical systems are often encountered, e.g., transverse Ising for fermions, coupled harmonic oscillators for bosons.

C. Equal off-diagonal couplings
When |A n+1,n | = |B n+1,n |, the above method fails.This is the case for interacting harmonic oscillators without the rotating wave approximation [15] and for quantum Ising models, and therefore of interest in a number of practical situations.We focus on the case where A and B are real, e.g., A n,n+1 = B n,n+1 = −ǫB n+1,n .The diagonal elements A nn and B nn are always different if there is a transverse field (fermions) or if the masses are finite (bosons), so it is reasonable to assume A nn = B nn (the Ising model without transverse field cannot be estimated using our method because excitations do not propagate).
It is convenient to introduce A simple calculation then shows that where we set A 01 = A N,(N +1) = 0.In some sense, this is similar to a 1D chain case.As the elements 1 ± |E k are already known, we learn A 11 and B 11 from which can also be written in terms of the known variables by inserting the completeness relation Eq. (4).Using and finally We found a simple and efficient method to identify the Hamiltonian of a system of coupled bosons or fermions.While the methods are completely analogous to the spin case [2,3], it is surprising that the higher number of parameters in the Hamiltonian that arises from the non-conservation of excitations can still be estimated using the same resources.Similarly to [5], we can deal with very weak system initialisation, such as thermal states.Therefore our methods can drastically reduce the required resources for system identification.Since we allow for site-dependent anisotropies, not only the coupling strengths but also the type of interaction is determined along the way.From the theory side, our work once more confirmed a type of 'area law' for estimation: looking only at the surface of short-range interacting systems can completely determine their Hamiltonian.It would be intriguing to see if this has direct connections with area laws of entanglement [19].While the efficiency of our method relied on the quadratic form of the Hamiltonian, we conjecture that even for models with true interaction terms, e.g., quartic terms in the Hamiltonian, all system parameters remain discoverable on the surface in principle.

Figure 2 :
Figure 2: Schematic overview of our estimation scheme for a chain of length N. The main requirements are (i) relatively fast single-qubit rotations and measurements on the first site, and (ii) a decoherence time of at least N 2 .
where E = diag {E 1 , ..., E N }, to have the desired form of Eq. (3).Note that the energy eigenvalues appear in pairs of positive E k and negative E k+N ≡ −E k values (k = 1, ..., N ).The matrix T consists of the right eigenvectors |E k of ηM as

Figure 3 :
Figure 3: The matrix M as a coupling graph on the nodes |n .

Table I :
Overview of indirect Hamiltonian estimation schemes.The interaction types represent the Pauli matrices involved in the spin-spin coupling, e.g., XX + Y Y + ∆ZZ + Z stands for a Hamiltonian of the form n,m Anm (XX + Y Y + ∆ZZ) n,m + n BnZn.
|B12| |A 21 | 2 − |B 21 | 2 .If |A 21 | = |B 21 |, we obtain |B 21 /A 21 | ≡ g by dividing the above two equations and then through |A 21 |(1 − g 2 ), both |A 21 | and |B 21 | are obtained.Because the phases of A and B are known, then we learn A 21 and B 21 .We can then express a similar set of equations for the next site 2|ηM |E k and N + 2|ηM |E k