Simplified Quantum Process Tomography

We propose and evaluate experimentally an approach to quantum process tomography that completely removes the scaling problem plaguing the standard approach. The key to this simplification is the incorporation of prior knowledge of the class of physical interactions involved in generating the dynamics, which reduces the problem to one of parameter estimation. This allows part of the problem to be tackled using efficient convex methods, which, when coupled with a constraint on some parameters allows globally optimal estimates for the Kraus operators to be determined from experimental data. Parameterising the maps provides further advantages: it allows the incorporation of mixed states of the environment as well as some initial correlation between the system and environment, both of which are common physical situations following excitation of the system away from thermal equilibrium. Although the approach is not universal, in cases where it is valid it returns a complete set of positive maps for the dynamical evolution of a quantum system at all times.


Introduction
Quantum process tomography (QPT) provides a means to specify the complete map of a set of input states of a quantum system (say the system state at time t = 0) to a set of output states (say the system state at some later time t > 0). It is an essential tool for characterizing the dynamics of quantum systems, and is especially useful for systems undergoing non-unitary evolution [1,2]. Therefore, it is useful not only for understanding the evolution of a quantum system coupled to possibly unknown environments, but also for applications in which quantum systems are manipulated for particular ends. In practical realizations of quantum information processing, for example, QPT is necessary to fully characterize the operation of quantum logic gates, which is, in turn, critical in improving the fidelity of quantum computing devices.
The main challenge in implementing QPT is the large number of required measurements. The positive map that represents the quantum process is specified by N 4 − N 2 parameters for a system Hilbert space of dimension N . The number of experiments and the computational power required to estimate the process therefore scale exponentially with the size of the system specified in qubits. This makes it difficult to realize in systems of even modest dimension.
A common experimental strategy to realize QPT involves the preparation of a complete set of input states that span the Hilbert space of the system. These input states act as probes of the quantum channel described by the positive map. After passing through the channel, the output states are reconstructed by means of quantum state tomography (QST). The channel process is estimated by inversion of the information contained in the difference between the input and the output states. QPT has been implemented in optical systems [3,4], atoms in optical lattices [5], NMR [6,7] and a solid state qubit [8].
A number of approaches exist that may reduce the number of input states that must be prepared. For example, ancilla-assisted process tomography [9,10] and direct characterization of quantum dynamics [11] both use fewer probes than the direct approach, but the probes need to be entangled, which may be either resource-expensive or impossible in some experimental situations. Individual diagonal elements of a process matrix can be estimated efficiently from the average fidelities of appropriately modified channels, while off-diagonal elements require ancillas and controlled quantum gates [12,13].
In this paper, we present an approach to reduce the size of the QPT problem, based on convex optimization. A common assumption in QPT is that the process is completely unknown a priori, it is a 'black box'. The idea introduced in this paper is based on the fact that in most cases in which QPT is wanted some knowledge will be available about what is going on inside the black box and that this prior knowledge could reduce the size of the problem significantly. One of the most elementary non-unitary processes is decoherence, caused by the coupling of a system to its environment. Decoherence manifests itself as phase-or amplitude damping. If it is known that decoherence is present, prior knowledge is available and the problem size can be reduced. The unknowns that need to be estimated are the coupling strength to the bath and the bath distribution function, neither of which can be measured directly. If a large enough and informative set of data is chosen these parameters can be estimated so that not only may the operator mapping the initial state onto a later one be determined, but also the evolution at all times can be estimated. Thus, from two sets of samples of the system state at two different moments in time, the complete evolution can be inferred. We test this approach experimentally by characterizing the rotational dephasing of a vibrational wavepacket in diatomic potassium molecules. In this case, the system size is more than two qubits, so that the feature of problem 3 size reduction is emphasized. This is to our knowledge the first implementation of process tomography in molecules.

Estimation procedure: incorporating prior knowledge
In order to illustrate the new features of this approach, we begin by developing the general theory of process tomography including an initially mixed quantum state of the environment and incorporating some prior knowledge of the Hamiltonian. We consider the unitary dynamics of the combined state ρ se comprising the system and its environment, where U (t) = e −iH t propagates the state over a time t according to the time-independent Hamiltonian H . We suppose that initially, the system and environment are not correlated, so that ρ se (0) = ρ s (0) ⊗ ρ e (0). Over time, environmental interactions decohere the system, and we examine this decoherence by considering the dynamics of the reduced density matrix ρ s for the system where Tr e indicates a partial trace over the environment. This partial trace is conveniently performed in the basis | j that diagonalizes the environment, taken to have dimension J . In this basis, the initial state of the environment can be written as a thermal ensemble and the system's dynamics reduce to the form where the Kraus operators E jk are given by The problem at hand is to estimate the probability distribution { p j } associated with the environment, as well as the Kraus operators themselves. This expression for the quantum process is known as the operator-sum representation, and is based on the set of Kraus operators E = {E kk |k, k = 1, . . . , J }. The basis can be transformed into a basis for which J n 2 for an n × n dimensional system density matrix. This makes the physical interpretation of the Kraus matrices more difficult, although it simplifies the mathematics. The Kraus operators are not, of course, dependent on the experimental configuration, assuming that the environment itself is not under the control of the experimenter. Nonetheless, knowledge of the Kraus operators implies complete knowledge of the decoherence process.
QPT is a means to estimate the Kraus operators from experimental data. A common procedure is to formulate the estimation algorithm as a convex optimization, that is, a search over all operators satisfying both the experimental data constraints and the mathematical form constraints of the problem in such a way as to guarantee that the solutions are globally optimal. 4 Further, such problems may make use of efficient algorithms to identify such optimal solutions. Unfortunately, the estimation of Kraus operators in the form above is not a convex optimization problem, for two reasons. First, the equality constraint κ k E † jk E jk = I n is not linear in E jk and second the objective function itself is not convex. Fortunately, the problem can be reformulated into a convex form by expanding the Kraus operators in a fixed basis of system space operators. Let {B i } be a set of n 2 operators that span the system Hilbert space. The Kraus operators can be expressed as a convex sum of these operators as Estimating the Kraus operators is then reduced to the estimation of a superoperator X , the representation of which in the fixed operators basis is given in terms of the matrix elements where the superoperator is restricted to be positive and trace-preserving. Estimation of the superoperator from the data set is a convex optimization problem [18]. The problem has now been made convex, but at the expense of an increase of the number of free variables. Since the size of the superoperator is n 2 × n 2 , and the size of the constraints is n 2 , the total number of free variables is n 4 − n 2 . This number needs to be matched by the same number of variables in the measurements to estimate the process. This method is 'expensive', even for experiments in which the collection of large data sets is straightforward.
Adoption of a few physically reasonable assumptions about the system-environment coupling, applied to the form of the Hamiltonian, greatly simplifies this task. In general we have where the first two terms represent the Hamiltonians of the system and environment, and where H se is that describing their interaction. Two particularly simple special cases arise if (i) In the first case, the interaction Hamiltonian H se can be decomposed in the same eigenbasis as the system Hamiltonian H s , so that no transitions between the system's energy eigenstates are induced by the coupling to the environment. The decoherence caused by the environment can then be described as 'pure dephasing', and the Kraus operators can be written in the form where the relation H s |n =hω n |n defines the eigenstates and eigenfrequencies of the system, and where the time-dependent coefficients µ n jk are given by In the second case, the dynamics of the environment are diagonal in the same basis {| j } that diagonalizes the initial state of the environment, meaning that different microstates of the initial thermal ensemble are not coupled. The Kraus operators then take the form where the K j are operators given by K j = j|(H e + H se )| j .
In the present work, we are fortunate that both of the above situations obtain; case (ii) exactly, and case (i) approximately. The resulting expression for the system dynamics can be written in the form where the κ jn are the system eigenvalues of the K j , By inspection, it is clear that (12) describes decoherence, since the off-diagonal elements of ρ s involve summations over exponential phase factors with incommensurate frequencies. In fact, it is straightforward to show that under conditions (i) or (ii) the system evolves into the state Because the final state is a convex sum over system-environment product states, it is not entangled. However, because the final system states are not orthogonal, it does possess nonclassical correlations in the form of nonzero discord. It is the existence of such correlations that causes the decoherence of the system itself.

Application: molecular vibrational decoherence by a rotational bath
A prototypical system-environment interaction that serves to test this method of simplified process tomography is that of vibrational-rotational coupling in a diatomic molecule. In this case, a vibrational wavepacket (that is, a coherent superposition of eigenstates of the internuclear potential energy in a particular electronic state) is gradually rendered incoherent by means of dephasing. This arises because of the coupling of the vibrational and rotational degrees of freedom of the molecule, via the moment of inertia. The moment of inertia changes dynamically as the molecule vibrates, leading to a change in the rotational frequency, which, in turn, modifies the vibrational frequency. The net effect on the vibrational wavepacket is to dephase the superposition, leading to a mixed vibrational state. The aim of QPT in this case is to determine by experiment the quantum process that characterizes this dephasing, leading to an estimation of the Kraus operators in the parameterized form described in the previous section. The experimental procedure consists of exciting a vibrational wavepacket in an excited electronic state of a homonuclear dimer and monitoring its evolution by measuring the timeand frequency-resolved fluorescence. This enables QST of the mode as it evolves, and, as we shall show, QST provides sufficient information to implement simplified QPT (SQPT).
The experimental apparatus has been described in detail elsewhere [14,15] and is illustrated in figure 1(a). Briefly, a vibrational wavepacket is prepared in the 1 + u state of K 2 through excitation by an ultrashort pulse, the vibrational period of the wavepacket is around 500 fs. The molecules are kept in a vapour cell at 400 • C. At this temperature around 150 rotational levels are occupied, which creates a particularly detrimental environment for the vibrational mode, and causes rapid dephasing due to the modulation of the centrifugal coupling to the rotational degree of freedom. During the oscillation of the wavepacket in the excited state potential, fluorescence is emitted when the electron makes a spontaneous transition to the ground state. The fluorescence is imaged onto a nonlinear crystal, and, by mixing the fluorescence in the crystal with an ultrashort pulse so that it is gated in time and frequency, a tomographically complete data set for the vibrational mode can be collected, as described in detail in [17].
From the time-frequency fluorescence map, we reconstruct the initial vibrational quantum state and the state at some later time. In order to do this, we assume that there is a separation of timescales, such that data for QST is collected over a short enough period that there is no dephasing. Figure 1(b) illustrates a typical fluorescence quantum beat pattern a particular wavelength, corresponding to the outer turning point of the vibrational wavepacket. The vertical dashed lines indicate the single vibrational periods that form the sampling windows for QST. The initial window occurs immediately following the termination of the exciting pulse, and the second at a sufficiently long time that significant dephasing has taken place.

7
QPT scales poorly with the size of the Hilbert space on which the process acts, and becomes quickly infeasible for larger Hilbert spaces. For the vibrational wavepacket discussed here, the Hilbert space dimension is 6, so that the number of unknowns in the problem is 1260. In order to implement 'blind' QPT, the measurements should consist of 36 orthogonal input states times 35 different settings, which we could do in principle by using differently shaped laser pulses. However, this is technically rather challenging, and becomes more so as we consider exciting more vibrational levels using broader bandwidth pulses. The size of this testbed system provides a useful model by which to test ways in which the size of the problem may be reduced. In the particular case chosen here, a number of simplifications are possible based on knowledge of the Hamiltonian of the system and its environment. In fact, there are more than might normally be afforded. However, this serves to illustrate in a stark manner the extent to which some prior knowledge can render a significant gain in determining the quantum process.
A first simplification is provided by the possible reduction of the number of input states. In standard process tomography experiments a complete set of input states is prepared, subjected to the process, and the output state is reconstructed. In this case the process is dephasing and a single input state is sufficient. Dephasing maps every element of the density matrix onto itself, so a single input state is sufficient as long as all elements of that density matrix are nonzero.
The unitary evolution in equation (1) is determined by the Hamiltonian wherep andq are the internuclear momentum and position operators, µ is the reduced mass, V (q) is the adiabatic vibrational potential of the electronic 1 + u state andĴ is the angular momentum operator of the molecule orthogonal to the internuclear axis. The Hamiltonian can be simplified by approximating the electronic potential as being harmonic and the ro-vibrational centrifugal coupling term can be developed up to first order in the displacement parameter η = q/q for small variations around the equilibrium internuclear separationq [19,20]. Writing the operator describing small displacements around this equilibrium position in terms of canonical creation and annihilation operatorsâ † andâ, [19] the Hamiltonian for a certain rotational level j iŝ H j =hωâ †â +h6Bη 2 j ( j + 1)â †â +h3Bη 2 j ( j + 1)(â †2 +â 2 ), (17) where ω is the vibrational frequency, B the rotational constant. The Heisenberg equation of motion within a rotational subspace can be solved to yield the time-dependent creation operator a † j (t) =â † (0)e iωt+iλ j t with λ j = 6η 2 B j( j + 1). Now the evolution of state | j, m |n within a jsubspace can be written aŝ The Kraus operator arising from coupling to rotational bath mode j is n j, m|U | j, m |n n|.
Since the quantum number m is degenerate, we can sum over m = − j, . . . , j which gives the Kraus operator 8 with N the normalization factor. The state of the system at time t is where p j is the probability of rotational state j being occupied in thermal equilibrium. The estimation of the superoperator is now narrowed down to estimation of the thermal probability distribution and the coupling constants {λ j }. Since at 400 • C around 150 levels are occupied, there are now only 300 free variables, namely the set { p j , λ j }. The number of variables may be reduced further by making use of the explicit form of the interaction Hamiltonian, in which λ j = λj ( j + 1), such that only a single coupling constant λ is required. This is a large reduction of the required number of data required to characterize a completely unknown process, and therefore enables a larger system subspace to be taken into consideration. We solve the estimation of the distribution { p j , λ} as a weighted least-squares problem with uniform weights. That is, a set of measurements is made to estimate the set of probabilities p αγ = Tr(O αγ ρ s (t)), where O αγ is the positive operator-valued measure (POVM) representing the measurement [15]- [17], the experimental configuration of which is labelled γ , and α is the measurement outcome. In our case γ represents the settings that enable detection of fluorescence at a particular wavelengths and time delay after the excitation of the vibrational wavepacket, i.e. λ = (ω fluor , τ ). The measurement outcome α is the strength of the fluorescence signal at that wavelength and delay. The experimentally determined probabilities is the set p emp αγ . Then convex optimization problem becomes subject to p j 0, j p j = 1, ρ s,f satisfies (20).
The procedure is to measure the initial quantum state, and propagate it using the parameterized form of the Kraus operators. The expected probabilities p αγ are subtracted from those empirically determined at the final time window p emp αγ , and the difference minimized by adjusting the parameters of the Kraus operators.
The estimation of the coupling coefficients λ j is not a convex problem, since it is a Hamiltonian parameter [18]. Therefore, we first estimate the optimal distribution of the bath distribution function for a fixed value of λ as a semidefinite program using available interiorpoint method solvers, such as those available via YALMIP [21]. This is repeated for a number of values of the coupling parameter that are between those physically accessible in the experiment. These values are set by the inverse of the minimum temporal resolution of the delay and the range of delays. The minimum of the objective function L over these two optimizations is taken as the global optimal, since the variation over λ is shown to be convex across the feasible range.

Experimental results
The quantum state of the system ρ s is reconstructed by means of QST [17] based on an optimal experiment design protocol, using 250 homogeneously distributed experimental settings γ . The initial (t = 0) and final (t of the Wigner function, whereas the darker, turquoise fringes are negative regions, indicative of nonclassical character. The initially angularly localized Wigner distribution diffuses around the classical trajectory as a result of rotational dephasing. This corresponds to a decay of the off-diagonal elements in the density matrix. The purity of the initial state is 0.4 and that of the final state is 0.2, which is still big enough to observe quantum beats in the signal, since the corresponding Wigner distribution is not yet spread out over the entire trajectory. The inner product between the estimated initial density matrix ρ est,i and the estimated final density matrix ρ est,f is Tr(ρ est,i ρ est,f ) = 0.15, which provides a simple quantification of the dephasing process.
Since the state tomography problem scales as N 2 , for this case the number of unknown variables in the state tomography problem is 36. The outcome of the estimation is shown in figure 3. The bath distribution function has been estimated for several different values of the coupling parameter. The values span the physically feasible range, from 6η 2 B = 2.1 × 10 7 s −1 to 2.9 × 10 7 s −1 , including the expected value, based on spectroscopic data, of 2.73 × 10 7 s −1 . Over this range, the objective L( p j , D) is convex, with a minimum value of 0.0021, as shown in the inset. At the minimum of this curve, the bath distribution function is very similar in shape to the rotational thermal distribution. The thermal bath distribution function is shown as the dashed line in figure 3. The optimal distribution is plotted in green, and the distributions for the extremal values of the coupling parameter as red and cyan. The peak of the optimal distribution is close to that of the rotational distribution at the experimental temperature of 400 • C.
On the one hand, the outcome is not surprising, because the protocol returns the thermal rotational distribution expected for a system in equilibrium. On the other hand, the outcome is nontrivial, because there is no direct access to the bath and the bath distribution function is estimated by measurements on the system alone.
The values of the objective L, as shown in the inset of figure 3, quantify the errors in the reconstructed measurement statistics. The differences between the empirical and estimated probabilities are so small that it is difficult to distinguish them by eye in a plot of the two sets of probabilities. We used a large data set for the estimation, and therefore the empirical   probabilities are a good approximation to the true probabilities. These facts together suggest that our reconstructed statistics are a faithful rendering of the true probabilities. Now it is possible to combine the estimated bath distribution { p j } and our prior knowledge-encapsulated by the form of the map in equation (20)-to construct the superoperator X . The real and imaginary parts of the 36 × 36 matrix are shown in figures 4(b) and (c). The superoperator is positive semi-definite and trace preserving. Although the visual information from the figures is limited, one can recognize two general characteristics of the dephasing process: the process leaves the diagonal density matrix elements unchanged, whereas the dephasing is larger for off-diagonal density matrix elements further away from the diagonal.
To evaluate the accuracy of the current approach, we should ideally compare this estimated process with the true process. Of course, we do not have access to the true process, but we can surmise that SQPT is reasonably accurate by comparing the final density matrix ρ s,f -as generated by propagating the initial reconstructed state ρ est,i according to our reconstructed process-with the final state ρ est,f as reconstructed directly via QST. We find the value of Tr{ρ est,f ρ s,f } to be 0.19, which is close to the purity Tr{ρ 2 est,f } = 0.2.

Conclusions
We have formulated and demonstrated a method to reconstruct completely the decoherence process that is caused by the coupling of a vibrational system to its rotational environment in a diatomic molecule. This method avoids the main problem of process tomography, which is the scaling of the problem. The way in which our approach avoids this problem is to apply some prior knowledge of the quantum process. This information consists of an assumption about the form of the coupling of the system to its environment, and provides a general means for simplification. The particular conditions under which it can be applied are broadly applicable-the linearity of the system-environment coupling, and the commutativity of the bath or system operators with the coupling Hamiltonian-though they are not universal. The problem is thereby reduced to one of parameter estimation, in which the number of parameters is vastly smaller than the number of elements of the process operators. Further, it provides a route for including realistic properties of the environment, such as a nonzero temperature, as well as constructing the quantum channel for any particular time evolution. The problem of estimating the environment distribution function at finite temperature has the form of a convex optimization problem, which makes the problem easy to solve numerically, and guarantees that the solution is the global optimum. Estimating the coupling parameters, which are part of the Hamiltonian, is not strictly convex, but there are not enough parameters to make searching the space straightforward and not computationally onerous. We anticipate that the formalism can, for example, be generalized for the case of a quantum logic gate in the presence of dephasing and to other systems in which nonunitary dynamics plays a role.