Single-qubit reaped quantum state tomography.

Quantum state tomography is the experimental procedure of determining an unknown state. It is not only essential for the verification of resources and processors of quantum information but is also important in its own right with regard to the foundation of quantum mechanics. Standard methods have been elusive for large systems because of the enormous number of observables to be measured and the exponential complexity of data post-processing. Here, we propose a new scheme of quantum state tomography that requires the measurement of only three observables (acting jointly on the system and pointer) regardless of the size of the system. The system is coupled to a "pointer" of single qubit, and the wavefunction of the system is "reaped" onto the pointer upon the measurement of the system. Subsequently, standard two-state tomography on the pointer and classical post-processing are used to reconstruct the quantum state of the system. We also developed an efficient and scalable iterative maximum likelihood algorithm to estimate states from statistically incomplete data.

Introduction

To develop reliable quantum resources and devices for quantum information processing, it is crucial to verify their actual performance. This is achieved at various levels, such as quantum process tomography[1] and quantum detector tomography[2,3], and at the most fundamental level is the quantum state tomography, which is the procedure of experimentally determining an unknown quantum state[4]. Quantum state tomography is of great interest in its own right with regards to the foundation of quantum mechanics as well.

In the standard formulation[5], quantum state tomography is accomplished by performing repeated measurements of numerous non-commuting observables on many systems prepared in the same states. As a matter of principle, if the set of non-commuting observables is complete and the measurements are repeated infinitely many times, one can build up a comprehensive description of the quantum state by post-processing the measurement statistics[6,7]. It is recapitulated by the three requirements of the standard quantum state tomography: (i) a complete set of observables to be measured (so-called “quorum”), (ii) accurate measurement statistics, and (iii) efficient post-processing. In practice, the requirement of measuring a complete set of observables causes overwhelming experimental obstacles, which affects the other requirements. Technical reasons and other difficulties may prevent some observables from being measured experimentally. For large systems, the number of required observables is exponentially large and places a serious limit on the number of repetitions of measurements (which is finite anyway in reality). Both issues lead to incomplete measurement statistics and/or limited accuracy of measurement statistics. Furthermore, even if reasonably accurate measurement statistics are attained, the complexity of post-processing itself is exponentially high for large systems. To overcome such difficulties in exact quantum state tomography, various statistical methods have been developed to estimate quantum states, such as the maximum likelihood estimation[8,9] and Bayesian estimation[10-12] methods. Notably, most statistical estimation methods, including the ML and Bayesian approaches, are highly nonlinear procedures and generally suffer from high complexity for large systems.

Here, we propose a new quantum state tomography scheme that requires the measurement of only three observables regardless of the system size[13]. In this scheme, the system is coupled to a “pointer” of a single two-level quantum system (i.e., “qubit”), and the wavefunction of the system is “reaped” onto the pointer upon the measurement of a single observable on the system. The subsequent standard quantum state tomography on the pointer and classical post-processing reconstruct the quantum state of the system, where the classical post-processing requires matrix inversion. We refer to this scheme as single-qubit reaped (or pointer-reaped) quantum state tomography. We have also developed an iterative maximum likelihood (ML) estimation algorithm that is adaptable to the single-qubit reaped scheme. The iterative ML estimation algorithm is demonstrated by numerical simulations with several interesting quantum states, such as the GHZ, W, and Dicke states. Furthermore, by matrix product state (MPS) representations, the iterative ML algorithm is scalable and provides an efficient method to obtain MPS estimates for the mixed states of large systems. The MPS pure state estimate for the mixed state determines the lower bound of the fidelity between the pure and mixed states and can be used to experimentally verify the purity of the laboratory-generated states[14].

Results

Exact tomography

Consider a system of n particles, each of which has dimension d, such that the total dimension of the system is Let be the computational basis of the Hilbert space. Suppose that we have an ensemble of such systems, identically prepared in the unknown state with the “wavefunction” , where is the set of complex numbers. We assume that without a loss of generality (a physical state cannot be a null vector). Our proposed scheme is illustrated in the two equivalent quantum circuits in Fig. 1. We discuss these procedures in the following order:

Figure 1

Two equivalent schematics of the single-qubit reaped quantum state tomography. (a) The system–pointer interaction is described by the p-dependent conditional phase shift on the pointer. (b) The system–pointer interaction is regarded as the pointer-controlled unitary operator on the system. The measurement on the system measures the eigenvalues x of the observable whereas the measurement on the pointer measures the eigenvalues of the Pauli operators , , or .

First, we select a qubit as the “pointer”. The pointer plays a central role in the proposed scheme. Initially, we prepare the pointer in the state , where and are the computational basis states of the pointer such that the initial state of the system plus pointer is given by .

Next, we couple the system and pointer for a certain time, which is assumed to be sufficiently short compared to the typical time scales of the system and pointer. This interaction can be described by a unitary operator of the form[15] where is an observable of the system. For the sake of physical implementation in actual experiments, one can take two different but equivalent views of . One can represent with the phase shift on the pointer conditioned on the system observable . To observe this more explicitly, let be the eigenstate of the observable belonging to the eigenvalue p and rewrite as with the p-dependent phase shift on the pointer. This interpretation is depicted in the quantum circuit representation in Fig. 1a and is analogous to the conventional von Neumann picture of the measurement of the observable . One important difference is that the pointer here is only of two dimensions and is insufficient to directly discriminate the N eigenvalues, p, of . On the other hand, noting that with being the identity operator and , one can regard it as a pointer-controlled unitary operator acting on the system. This picture is illustrated in the quantum circuit in Fig. 1b and is analogous to the quantum phase estimation circuit for a unitary transformation ( in the present case)[16]. Throughout this paper, we will mainly consider the latter interpretation for convenience. After the unitary interaction, the total state becomeswhere V is the matrix representation of in the computational basis,

We then measure the eigenvalues of the observable in the system. When the measurement outcome is x, the (unnormalized) pointer state is reduced to

Equation (3) reveals the key idea of the proposed scheme: the wavefunction appears in the two expansion coefficients and can be determined by the standard quantum state tomography by measuring three independent observables, that is, the Pauli operators , , and in the pointer. One tricky point is that naive two-state tomography does not fix the overall phase, which is necessary to fix the relative phases of for different values of x. We now provide a careful tomographic reconstruction procedure [see Eq. (6)] that is not hindered by this tricky issue.

Physically, the two-step procedure for the measurement of on the system and the subsequent quantum state tomography on the pointer is equivalent to the measurement of the eigenvalues of three observables, , , and . For the purpose of mathematical analysis of measurement outcomes and maximum likelihood estimation process (see below), it is convenient to describe the measurements using the projective POVM elementswhere , , and the index refers to the eigenstates of the Pauli operators , , and , respectively. The joint probabilities determine the ratio between the two coefficients,where Owing to the normalization constraint, the N relations in Eq. (5) are not independent of each other. Instead of directly imposing the normalization constraint, one can just determine the ratio . This casts the relation (5) to the following set of linear equationsfor . Given the experimentally determined measurement statistics , solving the linear equations yields the wavefunction (up to normalization). There are several dangerous cases in which Eq. (6) cannot provide a unique solution. Avoiding or overcoming them is addressed in “Methods”.

One simple example is to select the local basis such that where is the index of p when the eigenvalues are arranged in an ordered sequence. The computational basis and the eigenstates of are related by the quantum Fourier transform[17]. For a system consisting of qubits (), another valuable example is the system operator of the form where denotes the Pauli operator acting on the jth qubit. This leads to a pointer-controlled unitary operator

In this case, and are related to each other via the local Hadamard gates,with

Maximum likelihood estimation algorithm

Above, we have shown that, as a matter of principle, the single-qubit reaped scheme can successfully reconstruct quantum states. It assumes an idealistic situation where the probability distribution corresponding to the POVM elements can be inferred from measurements. It is possible only when the measurements are repeated infinitely many times, apart from other technical imperfections; finite repetitions give rise to statistical errors in the inferred probabilities . Obviously, the statistical errors become more severe as the system size n increases; recall the number of possible measurement outcomes (x, m). A popular method to overcome such an issue is to follow the maximum likelihood (ML) principle and seek the state that is most “likely” given the experimental observations rather than the actual (and impossible-to-infer) wavefunction[6-9,18]. In this section, we develop an iterative ML algorithm that can be combined with the single-qubit reaping scheme discussed above. We note controversies about the physically proper estimation of quantum states from the experimental data[11,18], and it would be valuable to develop other statistical methods, such as Bayesian approaches, that are adaptable to the present tomography scheme.

Consider an ensemble of F systems. Let be the number of experimental observations corresponding to the POVM element , such that The ideal situation corresponds to the limit , where the relative frequency gives the true probability . For finite size (), only estimates approximately. The observation statistics are governed by a multinomial distributionwhereis the probability of obtaining the result (x, m) on the condition that the system plus pointer is prepared in the state

We use the multinomial distribution as the likelihood function. Generally, the likelihood function should depend on the specific measurement apparatus and other experimental conditions. Here, we focus on the generic effects on statistical error, putting aside specific technical issues. The ML approach maximizes(up to irrelevant terms) over all possible states of the system with the normalization constraint. The wavefunction that maximizes the likelihood function satisfies the extremal equation (see “Methods” for details)where the matrix W is defined byand the x-dependent operator on the pointer by

Formally, is reminiscent of a similar operator (denoted by ) that appears in the iterative maximization algorithm adapted to the standard tomography scheme[9]. In our case, acts on the pointer and not on the system itself. In an ideal experiment where , the true wavefunction indeed gives the extremum solution, , as . In a realistic experiment with a finite-size ensemble (), in general , but is simply the wavefunction most likely for the given measurements data.

It should be noted that the operator depends functionally on the state through the probability , and hence the extremum Eq. (13) is nonlinear. Solving such a nonlinear equation is unviable, particularly for large systems (involving a large number of variables ). Instead, we have developed an iterative algorithm[9,18-20]. First, we need to choose an initial trial wavefunction. From the pointer state in Eq. (3) upon the measurement readout x, it follows that the probability is directly proportional to . This implies that is a reasonable choice. At each iterative step k, the wavefunction is updated using the mappingwhere the iteration generator is constructed from the matrix W in Eq. (14). Interestingly, the iteration procedure can be represented by the quantum circuit shown in Fig. 2, which illustrates the crucial role of the pointer from another perspective. The quantum circuit itself is not advantageous when one evaluates the iterations directly. However, as we will observe later, it clearly reveals the simple mathematical structure of the iteration generator , which permits the scalability of the iterative algorithm.

Figure 2

Quantum circuit interpretation of the maximum likelihood iteration. The solid dot indicates the “controlled”- (or ) acting only when the pointer is in the state whereas the open circle indicates “conditional”- on the pointer conditioned on the state of the system. Despite the quantum circuit interpretation, the iteration procedure is not linear as the operator depends on the trial state .

The convergence of no iterative ML algorithm has been analytically proven[18]. However, in standard ML approaches[21,22], numerical tests have demonstrated convergence for physically interesting states, and a diluted iterative algorithm is available when the convergence is critical[18]. Here, we demonstrate the algorithm numerically using several examples for a system of six qubits ( and ). The first example is the symmetric Dicke state where refers to the summation over all permutations of the qubits. We simulated the measurements for an ensemble of 24,000 systems () all prepared in the same state . The resulting relative frequencies, , of the measurement readouts (x, m) are shown in Fig. 3a. We then obtained the ML estimate for the measurement data () through 500 iterations in accordance with (16). As shown in Fig. 3 (b, blue curve), the infidelity between the states from consecutive iterations was already less than after 150 iterations. The fidelity, with the true wavefunction is larger then 0.997.

Figure 3

(a) Relative frequencies of the measurement readouts (x, m) from the simulation with an ensemble of 24,000 systems in the symmetric Dicke state with six qubits (, , ). (b) Convergence behaviors of the iterative maximization procedure for different system states (the W state, Dicke state, GHZ state, and the ground state of the transverse field Ising model in the ordered phase) exhibited by the fidelity between the states from consecutive iterations. For all the four cases, the fidelities between the resulting states and the true wavefunctions, respectively, are better than 0.99.

We performed similar simulations and made the ML estimates for the simulation results for the W-state the GHZ state and the ground state of the transverse-field Ising model in the ordered phase. Figure 3b corroborates the excellent convergence for all those cases. The fidelities between the ML estimates and the respective true wavefunctions were also as good as 0.99 or larger.

Scalability and mixed states

Each ML iteration in Eq. (16) involves the multiplication of exponentially large matrices and vectors, and the computational cost of many iterations for the desired accuracy may still be high for large systems. This can be overcome by means of matrix product state (MPS) and matrix product operator (MPO) representations (see “Methods”). We first examine the quantum circuit shown in Fig. 2 more closely to better understand the MPO structure of the iteration generator, . Let be the extended operator acting on the system and pointer, which results in when averaging over the pointer with the state . consists of the controlled-unitary operator and the conditional-unitary operator

The former is an MPO with a bond dimension of 2 when the coupling observable (and hence ) is local [Eq. (7) is an example]. The latter is also an MPO with a finite bond dimension provided that the input state is an MPS with a finite bond dimension because an MPS only has finite correlations[23,24]; see “Methods”. Therefore, , the product of three MPOs, should be an MPO with a finite bond dimension, and so is as it corresponds to a partial trace of an MPO. Currently, the operation of an MPO on an MPS can be efficiently evaluated[23,24]. In summary, if the laboratory states are MPS, the iteration generator is represented by an MPO, and the ML iterations in Eq. (16) can be updated efficiently. Recently, a formally similar iterative algorithm (from a different tomography scheme) powered by MPO and MPS representations has been demonstrated in detail[22].

Because only a polynomial number of parameters is required for the MPS representations, they span only a small portion of the entire Hilbert space. However, it is well known that many states relevant to quantum information processing, condensed matter physics, and other areas of physics exist in the MPS form. The ground states of the strongly correlated many-body Hamiltonians as well as the cluster states are notable examples.

Moreover, as was pointed out recently[14], the tomographic estimation of MPS pure states is valuable even when the system is in a mixed state. That is, it allows us to determine a lower bound on the fidelity between the pure state estimate and mixed states compatible with the experimental observations, thereby certifying the purity of the laboratory state via experiments. A scalable ML method has been proposed to directly reconstruct mixed states via local measurements[21,22], assuming that the states are close to a MPS. For their method, however, experimenters are required to measure many non-commuting observables whereas our scheme requires the measurement of only three observables , and regardless of the system size[13].

Discussion

A seemingly similar idea to couple the system with an ancillary system and measure only one observable (over the entire system plus ancilla) has been previously proposed[25]; this is the so-called ancilla-assisted quantum state tomography and has been demonstrated in recent experiments[26,27]. However, their scheme required the ancilla to be as large as or even larger than the system (one obvious advantage is that it can directly estimate the density matrix of the system). Moreover, no ML algorithm has been developed for their scheme.

The convergence of the ML iterations varies for different states. For example, it is noted in Fig. 3b that the convergence of the ML iterations is slower for the GHZ state (approximately 500 iterations are required for similar accuracy) than for other states. Recalling the massive and long distance entanglement in the GHZ state, this fact raises an interesting question about the relation between the convergence behavior of our ML iterations and the properties (such as multi-partite entanglement) of the state. We leave the relation as an inspiring open question for future works.

Methods

State-reconstruction equation

Here, we derive the state-reconstruction Eq. (5). We begin with the (unnormalized) pointer state in Eq. (1)where we have defined and for notational simplicity. We want to express the ratio in terms of the joint probabilities . The joint probabilities satisfy the following relationship:

Using the last two relations, one can obtain

This implies that the relative phase between and , which is the essential part for quantum coherence effects, can be extracted by combining the join probabilities on the left-hand side. More explicitly, we express it asand observe thatwhich is identical to Eq. (5). The physical implication of the above relation is that the probabilities and in the computational basis of the pointer give the relative magnitudes of and , whereas the probabilities and give the relative phases between them.

Dangerous cases

There are three dangerous cases where the wavefunction extraction scheme in Eq. (13) may not give a unique solution:

In the first case, is compatible with the computational basis, (). Then, are essentially eigenstates of , and the pointer state upon the measurement of becomes Because is an overall factor, it cannot be extracted.

In the second case, the unitary is block diagonal (possibly after simultaneous permutations of rows and columns) in a given basis. Suppose that with and operating on orthogonal subspaces and , respectively, of . Accordingly, any state is decomposed into . Upon the measurement of , the pointer is cast to for (). Therefore, in this case, one can assess by applying the wavefunction extraction scheme (6) for each sector . However, it is impossible to extract the phase relations between different sectors.

The third case is a special case where happens to be an eigenstate of (i.e., ) belonging to a degenerate eigenvalue p. Suppose that the pointer is in the state after the measurement of on the system. The two-state tomography can successfully extract the relative phase factor , and hence p. If p is non-degenerate, the eigenvalue itself uniquely identifies as its eigenstate. However, it is impossible if p is degenerate. Fortunately, this special case can be discerned experimentally because is independent of x, and for all x.

The first two cases can be avoided simply by properly choosing either the coupling operator or the computational basis .

Iterative ML algorithm

Here, we detail the maximization of the likelihood function over the entire Hilbert space. Because of the normalization constraint, it is more convenient to maximizewhere is the Lagrange multiplier. Suppose that the system was initially in a definite state and went through the unitary interaction with the pointer. Let be the pointer state upon the measurement outcome x on the system. Explicitly, it can be expressed as

The pointer state resulting from the general initial state of the system is related to by .

In terms of , the joint probability can be expressed as

For later use, it should be noted that its derivative with respect to has the form

Then, the extremal equation for the maximization problem (23) is given by

We define an x-dependent operator on the pointer by

Then, the extremal equation (27) is

Putting (24) into the above equation, we obtainwhich is identical to the matrix equation (13).

Matrix product states and operators

Consider a system of n particles, each of which has Hilbert space dimension d. We denote the computational basis state for as , where are the base d digits in x, .

An open boundary matrix product state (MPS)[23,24] is represented bywhere are the complex matrices, depending on the local state , and . Similarly, an open boundary matrix product operator (MPO) takes the formwhere are the basis operators of the Hilbert space of all linear operators acting on particle j; and are complex matrices ().

One can observe that the conditional operator is an MPO with a finite bond dimension provided that the state is an MPS with a finite bond dimension. Because an MPS has finite correlations, the probabilities are factorized as they are statistically independent of the uncorrelated parts[23,24]; we recall the base-d digits representation of x. This is also the case for the experimental observed frequencies . Therefore, the conditional operator is an MPO with a finite bond dimension.