Process Tomography for Systems in a Thermal State

We propose a new method for implementing process tomography that is based on the information extracted from temporal correlations between observables, rather than on state preparation and state tomography. As such, the approach is applicable to systems that are in a mixed state, and in particular to thermal states. We illustrate the method for an arbitrary evolution described by Kraus operators, as well as for simpler cases such as a general Gaussian channels, and qubit dynamics.


Introduction
Quantum process tomography deals with estimating the dynamics of unknown systems. It has long been studied due to its fundamental importance in the fields of quantum communication and quantum computation.
The common approach for implementing quantum process tomography is based on applying the dynamics on each element of a complete set of input states, and then performing tomographic measurement of the output states. This procedure allows one to completely reconstruct the superoperator representing the dynamics [1]. In recent years several improvements have been proposed for this method [2,3,4], which reduce its complexity in some cases.
Another approach is based on applying the dynamics on random states and comparing the results with the theoretical output state which would have been received had the transformation been purely unitary [5]. The major advantage of this method is that it is efficient (in the sense that it scales polynomially with the number of particles in the system). However, it does not reconstruct the superoperator representing the channel, but rather estimates the strength of the noise.
In the present work we propose a different method for performing process tomography of an unknown evolution, which is based on the information extracted from temporal correlations between observables, and as such does not require state preparation. Therefore, it would be useful in situations were one is restricted to a small set of physical input states and so standard preparation-evolution-measurement scheme cannot be applied. In fact, it can be used with unknown mixed states, and in particular with thermal states. Therefore it is potentially useful in "hot" systems whose state cannot be controlled with current technology, or cooled to the ground state.
The proposed method is based on two key points: (a) we utilize a proper set of temporal correlations between observables at two consecutive instances of time, say t 2 and t 1 , t 2 > t 1 , that encodes the dynamical evolution of the system, and (b) we employ weak measurements, rather than ordinary disturbing measurements, in order to measure such correlations. We cannot use ordinary measurements to observe the above temporal correlations because the interaction with the system at the earlier time will disrupt the correlations. Nevertheless, as we show, by using weak measurements which barely affect the system, such correlations can be measured.
Our method is applicable for general discrete systems whose evolution is described by a set of Kraus operators, as well as continuous systems whose evolution is linear with respect to a given set of operators. The rest of the paper is organized as follows: In the next section we define the temporal covariance matrix which will be sufficient for determining the above type of general evolution. We then show how the temporal covariance can be measured using weak measurements. In section 3 we show how to find the Kraus operators for general discrete systems, in section 4 we show how to find the dynamics for Gaussian channels, in section 5 we demonstrate the proposed method on qubit dynamics and in section 6 we present a numerical study.

Measurement of a temporal covariance matrix
Let us begin by defining the two-point temporal covariance matrix: where {B i } are Hermitian operators. This unusual definition of the covariance matrix, where the measurements are taken at different times, proves to be essential for our method.
In order to measure the covariance matrix one has to measure the observables B i (t) and the correlations {B i (t 1 ) , B j (t 2 )}.
While the measurement of the observables is straightforward, the measurement of the correlations is not trivial and the rest of this section is devoted to it.
The general scheme for the measurement of a single correlation {B i (t 1 ) , B j (t 2 )} involves two measurement devices: one measures B i at time t 1 , and the other measures B j at time t 2 . In order to measure the correlation one has to measure a joint operator of both measurement devices. Our method can be implemented in numerous ways. For simplicity, we shall demonstrate it with spin pointers and assume unitary time evolution. In this case, the Hamiltonian of the coupling between the system and the n'th pointer measuring Assuming both pointers initially point at |↑ x ≡ 1 √ 2 (|↓ z + |↑ z ), we denote the initial state as ρ in ≡ a p a |ψ a,sys , ↑ x , ↑ x ψ a,sys , ↑ x , ↑ x |, where the label 'sys' stands for 'system'. Setting ρ in = ρ (t 1 ), the final state, ρ fi , is given by where U (t 2 ) is the time propagator. The pointers' expectation value and variance are therefore where |f (B i (t 1 ) , B j (t 2 ))| ≤ 1 3 B 4 (the notion · stands for operator norm [6]), for more details see appendix III. This procedure can be easily generalized to non-unitary time evolutions and pointers of general dimension. Similar calculations for continuous pointers are found in [7]. Note that while the error described in Eq. (5) is random, the error described in Eq. (4) is systematic. This means that in principle, for a given system, one can calculate f (B i (t 1 ) , B j (t 2 )) (either exactly or perturbatively), fix the systematic error and thus dramatically improve the efficiency of this process. This direction will be discussed further in section 6.
It can be shown that following a single correlation measurement procedure (for short times), i.e., the state is hardly influenced by the measurement, hence we refer to this measurement as "weak" [8,9] (Note that we do not use post selection). This is the key feature that enables this measurement method to work. As shown in equations (4) and (5) the "weakness" of the measurement comes with a price, every single weak measurement is highly inaccurate. This can be compensated by a large number of measurements. The error after N such measurements of a single correlation is δ = max )| 2 , hence in order to reach this error the optimal measurement strength is = δ |f (B i (t 1 ) , B j (t 2 ))| −1 and the number of In appendix II we present an alternative method for measuring the correlations. This alternative method involves only single pointer measurements, and so we believe this new method may be easier to implement.

Constructing the dynamics for discrete systems
In this section we show how the dynamics can be estimated using the two-point temporal covariance matrix. Given a D level system, we wish to define a complete basis of operators {B a } ≡ {B 0 , B i } (i ∈ 1..D 2 − 1 ), such that B 0 = 1 D×D and {B i } are chosen to be Hermitian matrices. In the Heisenberg picture, since the set of operators is complete, and superoperators representing quantum channels are linear, the following equation holds Substituting this equation in Eq. (1) we obtain This equation defines the time evolution of the covariance matrix. M could be retrieved from the covariance matrices simply by In general, σ (t 0 , t 0 ) might be singular. This happens if and only if the initial density matrix is singular (see appendix I for a proof). In this case the density operator doesn't sample all the states in the Hilbert space, and so it is impossible to gain a complete knowledge of the evolution. However, it is worth mentioning that when picking a random initial state over a continuous distribution, the chance of producing a singular density matrix is zero.
From Eq. (7) we get Note that the estimation of χ (t, t 0 ) does not require measurements in addition to the ones needed in order to determine σ (t, t 0 ). Since M and χ entirely encode the dynamics, one can fully estimate the dynamics by measuring the two-point temporal covariance matrix.
Since this method relies on the measurement of σ ij (t 0 , t 0 ) and σ ij (t, t 0 ) for every i, j ∈ 1..D 2 − 1 its complexity grows with D as O D 4 . For n qudits (d level systems) D = d n and thus the complexity grows with the number of qudits as O d 4n , the same as in [1].
Note that the described method works for general non-singular states, and in particular maximally mixed states and thermal states. An alternative procedure, based on covariance matrices defined at one time only, would not work with general states.
While M and χ entirely encode the dynamics, it would be convenient to describe the channel using the Kraus representation. The Kraus representation of a superoperator is where the following constraint holds By transferring the time dependence from the Schrodinger picture to the Heisenberg picture, demanding that the expectation value of a general operatorÔ remain unaffected by the picture, one can deduce the equivalent equation for operators: Since the set of operators {B a } ≡ {B 0 , B i } (B a ≡ B a (t 0 )) is complete, the Kraus operators can be represented as (summation on double indices is assumed) where u µ are complex vectors of coefficients and the following relations hold: where f abc and g abc determine the structure of the basis. Substituting Eq. (14)- (16) in Eq. (13) we obtain: In order to bring this equation to the form of Eq. (7) we set M ij (t, t 0 ) ≡ Next, we use additional constraints on u aµ 's which are derived by substituting Eq. (14) in Eq. (12) which means The symmetrical part of M , anti-symmetrical part of M and the trace of M , together with χ and the two constraints of Eq. (19) and (20) form 6 equations for 6 unknown objects which are u * 0µ u 0µ , Re u * iµ u 0µ , Im u * iµ u 0µ , Re u * iµ u jµ , Im u * iµ u jµ and u * iµ u iµ . The solution of this set of equations fully determines the matrix u µ ⊗ u † µ . Since this matrix is Hermitian it can be diagonalized. Using the resulting sets of D orthonormal eigenvectors {v µ } and corresponding eigenvalues {λ µ }, we obtain u µ = λ µ v µ (no summation on µ), and thus a possible set of D 2 Kraus operators that govern the system's dynamics is given by Given the ability to measure M and χ at small time intervals it is possible to estimate the time derivative of the Kraus operators. For Markovian systems this allows one to estimate the Lindblad equation as described in [10].

Constructing the dynamics for Gaussian channels
In general, our proposed method requires the use of a complete set of operators. Nevertheless, in special cases it can be applied using a limited set of observables. One important example for this statement is the class of systems described by quadratic Hamiltonians H = λ ij η i η j + α i η i (η ≡ (x,p) wherex andp are arrays of coordinate and conjugate momenta respectively). The evolution of a subset of the system is then described as a Gaussian channel [11,12] that dictates for the conjugate coordinates an evolution of the type This equation is identical to Eq. (7). Therefore the process of estimating the matrix M is the same as described in the previous section. For Gaussian states M is sufficient in order to calculate the density matrix evolution in time. While in the case of discrete systems the density matrix must be non singular in order for σ (t 0 , t 0 ) to be non singular, in the case of the Gaussian channels, the method would work for every state, even for states that do not sample the whole basis of states. The proof goes as follows: according to a variant of the uncertainty principle [13,14] for every covariance matrix of canonical operators there exists a symplectic matrix S and a diagonal matrix W where W ii ≥ 1 2 which satisfies W = Sσ (t, t) S T . Since S is symplectic it is invertible and so σ (t, t) = S −1 W S T −1 . In this form its clear that σ (t, t) −1 = S T W −1 S always exists, and in particular for t = t 0 . The meaning of this feature is that in the case of Gaussian channels, our method is truly state independent.
Since this method relies on the measurement of σ ij (t 0 , t 0 ) and σ ij (t, t 0 ) for every i, j ∈ {1..2n} (where n is the number of particles) its complexity grows with the number of particles as O n 2 .
It is interesting to remark that in this Gaussian case, there is an analogy between the proposed method and Gaussian state tomography. In the "spatial problem" it is well known that the spatial correlations encoded in the covariance matrix fully determine a Gaussian state [15]. In our case we see that the Gaussian dynamics is encoded in the temporal correlations which form a temporal correlation matrix.
We believe Gaussian channel could be among the most practical channels to estimate using our method. A simple realization of a single particle Gaussian channel would be an ion inside a trap. The spatial degrees of freedom will be regarded as the system and the spin will be used as the measurement device. Before the beginning of the estimation process the system can conveniently be prepared in a thermal state, while the measurement device has to be in a pure state. The interaction between the system and the measurement device could be implemented using lasers. Following each measurement the spin would have to be brought back to the pure state, however no initialization is required for the system. The only time consuming process in this procedure is the initialization of the spin.

Example: Constructing the Kraus operators for qubit dynamics
For qubit dynamics we choose the natural operator basis B 0 = 1 The basis structure is therefore and the rest of the coefficients vanish. Using the explicit form of f abc and g abc we calculate Recalling that in this case D = 2 we obtain From the constraints of Eq. (19) and (20) we obtain two additional equations: The solution of equations (26)-(31) is Re u * iµ u jµ = Equations (32)-(36) construct the matrix u µ ⊗ u † µ which is used to calculate a possible set of Kraus operators as explained above. This result can be easily generalized to n interacting qubits.

Numerical study
We have simulated the channel estimation process for several channels and initial states. In Figs. (1-3) we show simulations of a phase-damping channel estimation process. Fig. (1) presents a single channel estimation process, Fig. (2) presents histograms of estimated values of M 11 , M 12 and M 33 obtained by 1, 000 different simulations and Fig. (3) presents ∆M ≡ M estimated − M theoretical as a function of the number of measurements, N , for numerous measurement couplings.
As shown in Section 2, for the method presented here, the number of measurements required to reach an error of δ for each element of M scales as O δ −4 . For the standard method [1], on the other hand, the number of measurements scales as O δ −2 . In the particular example, analysed numerically here and presented in Fig. (1), we find that the parameters of the phase-damping channel can be estimated with an error of 0.1, using 2, 500 measurements per correlation, and 22, 500 measurements in total. This should then be compared with the standard estimation method, that by the scaling argument above is expected to be more efficient. For the same channel and up to the same error, we find that in total only 1, 000 measurements are required.

Discussion
In the present work we have proposed a method for implementing process tomography which is based on the information encoded in temporal correlations. We have shown that such correlations, embodied in the temporal correlations matrix, provide a general state-independent method to reconstruct the dynamical evolution in terms of Kraus operators. The complexity of our method grows exponentially with the number of particles in the discrete case (as in the usual approach), and quadratically with the number of particles in the Gaussian channel case.
Compared with the standard method, our proposed approach is clearly less efficient. Nevertheless, since it does not require state preparation, that is essential for ordinary process tomography, it is applicable in a wider class of cases. This includes systems in some general mixed states and in particular in thermal states. In this respect, we hope that the present method could be beneficial in various experimental systems, such as NEMS [16,17] and mesoscopic systems.
In addition, for systems where the time scales are very short, and the coupling constant between the system and the measuring device is small, it would be difficult to implement strong measurements. Since our method utilizes weak measurement, we believe our method would be suitable for such systems as well. Weak measurements have been recently realized in various systems [18,19].
To further improve the efficiency of our approach, it would be necessary to develop methods to reduce the error discussed in section 2. A significant part of this error is systematic, and may be further reduced by using either analytical or perturbative methods. For example, when dealing with qubit dynamics and the evolution is known to be unitary, one can calculate f (B i (t 1 ) , B j (t 2 )) and the higher order corrections in terms of the measured correlation and thus eliminate the systematic error. For general, non-unitary, dynamics we believe that similar results can achieved perturbatively. Finally, it would be interesting to investigate whether the present method can be combined with existing efficient methods for selective estimation of dynamics [3,4].