Quantum tomography

Using Born's rule, one can derive the simplest form of quantum tomography. Generally, being in a pure state is not known, and a state may be mixed. In this case, many different measurements will have to be performed, many times each. To fully reconstruct the density matrix for a mixed state in a finite-dimensional Hilbert space, the following technique may be used.

Born's rule states ${\displaystyle \mathrm {P} (E_{i}|\rho )=\mathrm {Trace} (E_{i}\rho )}$, where ${\displaystyle E_{i}}$ is a particular measurement outcome projector and ${\displaystyle \rho }$ is the density matrix of the system. Given a histogram of observations for each measurement, one has an approximation ${\displaystyle p_{i}}$ to ${\displaystyle \mathrm {P} (E_{i}|\rho )}$ for each ${\displaystyle E_{i}}$.

Given linear operators ${\displaystyle S}$ and ${\displaystyle T}$, define the inner product

${\displaystyle S\cdot T=\mathrm {Tr} [S^{\dagger }T]={\vec {S}}^{\dagger }{\vec {T}}}$

where ${\displaystyle {\vec {T}}}$ is representation of the ${\displaystyle T}$ operator as a column vector and ${\displaystyle {\vec {S}}^{\dagger }}$ a row vector such that ${\displaystyle {\vec {S}}^{\dagger }{\vec {T}}}$ is the inner product in ${\displaystyle \mathbb {C} ^{d}}$ of the two.

Define the matrix ${\displaystyle A}$ as

${\displaystyle A={\begin{pmatrix}{\vec {E}}_{1}^{\dagger }\\{\vec {E}}_{2}^{\dagger }\\{\vec {E}}_{3}^{\dagger }\\\vdots \end{pmatrix}}}$.

Here E_i is some fixed list of individual measurements (with binary outcomes), and A does all the measurements at once.

Then applying this to ${\displaystyle {\vec {\rho }}}$ yields the probabilities:

${\displaystyle A{\vec {\rho }}={\begin{pmatrix}{\vec {E}}_{1}^{\dagger }{\vec {\rho }}\\{\vec {E}}_{2}^{\dagger }{\vec {\rho }}\\{\vec {E}}_{3}^{\dagger }{\vec {\rho }}\\\vdots \end{pmatrix}}={\begin{pmatrix}E_{1}\cdot \rho \\E_{2}\cdot \rho \\E_{3}\cdot \rho \\\vdots \end{pmatrix}}={\begin{pmatrix}\mathrm {P} (E_{1}|\rho )\\\mathrm {P} (E_{2}|\rho )\\\mathrm {P} (E_{3}|\rho )\\\vdots \end{pmatrix}}\approx {\begin{pmatrix}p_{1}\\p_{2}\\p_{3}\\\vdots \end{pmatrix}}={\vec {p}}}$.

Linear inversion corresponds to inverting this system using the observed relative frequencies ${\displaystyle {\vec {p}}}$ to derive ${\displaystyle {\vec {\rho }}}$ (which is isomorphic to ${\displaystyle \displaystyle \rho }$).

This system is not going to be square in general, as for each measurement being made there will generally be multiple measurement outcome projectors ${\displaystyle E_{i}}$. For example, in a 2-D Hilbert space with 3 measurements ${\displaystyle \sigma _{x},\sigma _{y},\sigma _{z}}$, each measurement has 2 outcomes, each of which has a projector E_i, for 6 projectors, whereas the real dimension of the space of density matrices is (2⋅2²)/2=4, leaving ${\displaystyle A}$ to be 6 x 4. To solve the system, multiply on the left by ${\displaystyle A^{T}}$:

${\displaystyle A^{T}A{\vec {\rho }}=A^{T}{\vec {p}}}$.

Now solving for ${\displaystyle {\vec {\rho }}}$ yields the pseudoinverse:

${\displaystyle {\vec {\rho }}=(A^{T}A)^{-1}A^{T}{\vec {p}}}$.

This works in general only if the measurement list E_i is tomographically complete. Otherwise, the matrix ${\displaystyle A^{T}A}$ will not be invertible.

In infinite dimensional Hilbert spaces, e.g. in measurements of continuous variables such as position, the methodology is somewhat more complex. One notable example is in the tomography of light, known as optical homodyne tomography. Using balanced homodyne measurements, one can derive the Wigner function and a density matrix for the state of the light.

One approach involves measurements along different rotated directions in phase space. For each direction ${\displaystyle \theta }$, one can find a probability distribution ${\displaystyle w(q,\theta )}$ for the probability density of measurements in the ${\displaystyle \theta }$ direction of phase space yielding the value ${\displaystyle q}$. Using an inverse Radon transformation (the filtered back projection) on ${\displaystyle w(q,\theta )}$ leads to the Wigner function, ${\displaystyle \mathrm {W} (x,p)}$,[5] which can be converted by an inverse Fourier transform into the density matrix for the state in any basis.[4] A similar technique is often used in medical tomography.

Field amplitudes or quadratures with high efficiencies can be measured with photodetectors together with temporal mode selectivity. Balanced homodyne tomography is a reliable technique of reconstructing quantum states in the optical domain. This technique combines the advantages of the high efficiencies of photodiodes in measuring the intensity or photon number of light, together with measuring the quantum features of light by a clever set-up called the homodyne tomography detector. This is explained by the following example. A laser is directed onto a 50-50% beamsplitter, splitting the laserbeam into two beams. One is used as local oscillator (LO) and the other is used to generate photons with a particular quantum state. The generation of quantum states can be realized, e.g. by directing the laser beam through a frequency doubling crystal [6] and then onto a parametric down-conversion crystal. This crystal generates two photons in a certain quantum state. One of the photons is used as a trigger signal used to trigger (start) the readout event of the homodyne tomography detector. The other photon is directed into the homodyne tomography detector, in order to reconstruct its quantum state. Since the trigger and signal photons are entangled (this is explained by the Spontaneous parametric down-conversion article), it is important to realize, that the optical mode of the signal state is created nonlocal only when the trigger photon impinges the photodetector (of the trigger event readout module) and is actually measured. More simply said, it is only when the trigger photon is measured, that the signal photon can be measured by the homodyne detector.

Now consider the homodyne tomography detector as depicted in figure 4 (figure missing). The signal photon (this is the quantum state we want to reconstruct) interferes with the local oscillator, when they are directed onto a 50-50% beamsplitter. Since the two beams originate from the same so called master laser, they have the same fixed phase relation. The local oscillator must be intense, compared to the signal so it provides a precise phase reference. The local oscillator is so intense, that we can treat it classically (a = α) and neglect the quantum fluctuations. The signal field is spatially and temporally controlled by the local oscillator, which has a controlled shape. Where the local oscillator is zero, the signal is rejected. Therefore, we have temporal-spatial mode selectivity of the signal. The beamsplitter redirects the two beams to two photodetectors. The photodetectors generate an electric current proportional to the photon number. The two detector currents are subtracted and the resulting current is proportional to the electric field operator in the signal mode, depended on relative optical phase of signal and local oscillator.

Since the electric field amplitude of the local oscillator is much higher than that of the signal the intensity or fluctuations in the signal field can be seen. The homodyne tomography system functions as an amplifier. The system can be seen as an interferometer with such a high intensity reference beam (the local oscillator) that unbalancing the interference by a single photon in the signal is measurable. This amplification is well above the photodetectors noise floor.

The measurement is reproduced a large number of times. Then the phase difference between the signal and local oscillator is changed in order to ‘scan’ a different angle in the phase space. This can be seen from figure 4. The measurement is repeated again a large number of times and a marginal distribution is retrieved from the current difference. The marginal distribution can be transformed into the density matrix and/or the Wigner function. Since the density matrix and the Wigner function give information about the quantum state of the photon, we have reconstructed the quantum state of the photon.

The advantage of this method is that this arrangement is insensitive to fluctuations in the frequency of the laser.

The quantum computations for retrieving the quadrature component from the current difference are performed as follows.

The photon number operator for the beams striking the photodetectors after the beamsplitter is given by:

${\displaystyle {\hat {n}}_{i}={\hat {a}}_{i}^{\dagger }{\hat {a}}_{i}}$,

where i is 1 and 2, for respectively beam one and two. The mode operators of the field emerging the beamsplitters are given by:

${\displaystyle {\hat {a}}_{1}=2^{-1/2}({\hat {a}}-\alpha _{LO})}$

${\displaystyle {\hat {a}}_{2}=2^{-1/2}({\hat {a}}+\alpha _{LO})}$

The ${\displaystyle {\hat {a}}}$ denotes the annihilation operator of the signal and alpha the complex amplitude of the local oscillator. The number of photon difference is eventually proportional to the quadrature and given by:

${\displaystyle {\hat {n}}_{21}={\hat {n}}_{2}-{\hat {n}}_{1}=\alpha _{LO}^{*}{\hat {a}}+\alpha _{LO}{\hat {a}}^{\dagger }}$,

Rewriting this with the relation:

${\displaystyle {\hat {q}}=2^{-1/2}({\hat {a}}^{\dagger }+{\hat {a}})}$

Results in the following relation:

${\displaystyle {\hat {n}}_{21}=2^{1/2}|{\hat {\alpha }}_{LO}|{\hat {q}}_{\theta }}$,

where we see clear relation between the photon number difference and the quadrature component ${\displaystyle {\hat {q}}_{\theta }}$. By keeping track of the sum current, one can recover information about the local oscillator's intensity, since this is usually an unknown quantity, but an important quantity for calculating the quadrature component ${\displaystyle {\hat {q}}_{\theta }}$.

One of the primary problems with using linear inversion to solve for the density matrix is that in general the computed solution will not be a valid density matrix. For example, it could give negative probabilities or probabilities greater than 1 to certain measurement outcomes. This is particularly an issue when fewer measurements are made.

Another issue is that in infinite dimensional Hilbert spaces, an infinite number of measurement outcomes would be required. Making assumptions about the structure and using a finite measurement basis leads to artifacts in the phase space density.[4]

Maximum likelihood estimation (also known as MLE or MaxLik) is a popular technique for dealing with the problems of linear inversion. By restricting the domain of density matrices to the proper space, and searching for the density matrix which maximizes the likelihood of giving the experimental results, it guarantees the state to be theoretically valid while giving a close fit to the data. The likelihood of a state is the probability that would be assigned to the observed results had the system been in that state.

Suppose the measurements ${\displaystyle \{|y_{j}\rangle \langle y_{j}|\}}$ have been observed with frequencies ${\displaystyle f_{j}}$. Then the likelihood associated with a state ${\displaystyle {\hat {\rho }}}$ is

${\displaystyle L({\hat {\rho }})=\prod _{j}\langle y_{j}|{\hat {\rho }}|y_{j}\rangle ^{f_{j}}}$

where ${\displaystyle \langle y_{j}|{\hat {\rho }}|y_{j}\rangle }$ is the probability of outcome ${\displaystyle y_{j}}$ for the state ${\displaystyle {\hat {\rho }}}$.

Finding the maximum of this function is non-trivial and generally involves iterative methods.[7][8] The methods are an active topic of research.

Maximum likelihood estimation suffers from some less obvious problems than linear inversion. One problem is that it makes predictions about probabilities that cannot be justified by the data. This is seen most easily by looking at the problem of zero eigenvalues. The computed solution using MLE often contains eigenvalues which are 0, i.e. it is rank deficient. In these cases, the solution then lies on the boundary of the n-dimensional Bloch sphere. This can be seen as related to linear inversion giving states which lie outside the valid space (the Bloch sphere). MLE in these cases picks a nearby point that is valid, and the nearest points are generally on the boundary.[3]

This is not physically a problem, the real state might have zero eigenvalues. However, since no value may be less than 0, an estimate of an eigenvalue being 0 implies that the estimator is certain the value is 0, otherwise they would have estimated some ${\displaystyle \epsilon }$ greater than 0 with a small degree of uncertainty as the best estimate. This is where the problem arises, in that it is not logical to conclude with absolute certainty after a finite number of measurements that any eigenvalue (that is, the probability of a particular outcome) is 0. For example, if a coin is flipped 5 times and each time heads was observed, it does not mean there is 0 probability of getting tails, despite that being the most likely description of the coin.[3]

Bayesian mean estimation (BME) is a relatively new approach which addresses the problems of maximum likelihood estimation. It focuses on finding optimal solutions which are also honest in that they include error bars in the estimate. The general idea is to start with a likelihood function and a function describing the experimenter's prior knowledge (which might be a constant function), then integrate over all density matrices using the product of the likelihood function and prior knowledge function as a weight.

Given a reasonable prior knowledge function, BME will yield a state strictly within the n-dimensional bloch sphere. In the case of a coin flipped N times to get N heads described above, with a constant prior knowledge function, BME would assign ${\displaystyle \scriptstyle {\frac {1}{N+2}}}$ as the probability for tails.[3]

BME provides a high degree of accuracy in that it minimizes the operational divergences of the estimate from the actual state.[3]

The number of measurements needed for a full quantum state tomography for a multi-particle system scales exponentially with the number of particles, which makes such a procedure impossible even for modest system sizes. Hence, several methods have been developed to realize quantum tomography with fewer measurements.

The concept of matrix completion and compressed sensing have been applied to reconstruct density matrices from an incomplete set of measurements (that is, a set of measurements which is not a quorum). In general, this is impossible, but under assumptions (for example, if the density matrix is a pure state, or a combination of just a few pure states) then the density matrix has fewer degrees of freedom, and it may be possible to reconstruct the state from the incomplete measurements.[9]

Permutationally Invariant Quantum Tomography[10] is a procedure that has been developed mostly for states that are close to being permutationally symmetric, which is typical in nowadays experiments. For two-state particles, the number of measurements needed scales only quadratically with the number of particles. [11] Besides the modest measurement effort, the processing of the measured data can also be done efficiently: It is possible to carry out the fitting of a physical density matrix on the measured data even for large systems. [12] Permutationally Invariant Quantum Tomography has been combined with compressed sensing in a six-qubit photonic experiment.[13]

A quantum process, also known as a quantum dynamical map, ${\displaystyle {\mathcal {E}}(\rho )}$, can be described by a completely positive map

${\displaystyle {\mathcal {E}}(\rho )=\sum _{i}A_{i}\rho A_{i}^{\dagger }}$,

where ${\displaystyle \rho \in {\mathcal {B(H)}}}$, the bounded operators on Hilbert space; with operation elements ${\displaystyle \displaystyle A_{i}}$ satisfying ${\displaystyle \textstyle \sum _{i}A_{i}^{\dagger }A_{i}\leq I}$ so that ${\displaystyle \mathrm {Tr} [{\mathcal {E}}(\rho )]\leq 1}$.

Let ${\displaystyle \displaystyle \{E_{i}\}}$ be an orthogonal basis for ${\displaystyle {\mathcal {B(H)}}}$. Write the ${\displaystyle \displaystyle A_{i}}$ operators in this basis

${\displaystyle \displaystyle A_{i}=\sum _{m}a_{im}E_{m}}$.

This leads to

${\displaystyle {\mathcal {E}}(\rho )=\sum _{m,n}\chi _{mn}E_{m}\rho E_{n}^{\dagger }}$,

where ${\displaystyle \chi _{mn}=\sum _{i}a_{mi}a_{ni}^{*}}$.

The goal is then to solve for ${\displaystyle \displaystyle \chi }$, which is a positive superoperator and completely characterizes ${\displaystyle {\mathcal {E}}}$ with respect to the ${\displaystyle \displaystyle \{E_{i}\}}$ basis.[16][17]

SQPT approaches this using ${\displaystyle d^{2}}$ linearly independent inputs ${\displaystyle \rho _{j}}$, where ${\displaystyle d}$ is the dimension of the Hilbert space ${\displaystyle {\mathcal {H}}}$. For each of these input states ${\displaystyle \rho _{j}}$, sending it through the process gives an output state ${\displaystyle {\mathcal {E}}(\rho )}$ which can be written as a linear combination of the ${\displaystyle \rho _{k}}$, i.e. ${\displaystyle \textstyle {\mathcal {E}}(\rho _{j})=\sum _{k}c_{jk}\rho _{k}}$. By sending each ${\displaystyle \rho _{j}}$ through many times, quantum state tomography can be used to determine the coefficients ${\displaystyle c_{jk}}$ experimentally.

Write

${\displaystyle E_{m}\rho _{j}E_{n}^{\dagger }=\sum _{k}B_{m,n,j,k}\rho _{k}}$,

where ${\displaystyle B}$ is a matrix of coefficients. Then

${\displaystyle \sum _{k}c_{jk}\rho _{k}={\mathcal {E}}(\rho _{j})=\sum _{m,n}\chi _{m,n}E_{m}\rho _{j}E_{n}^{\dagger }=\sum _{m,n}\sum _{k}\chi _{m,n}B_{m,n,j,k}\rho _{k}}$.

Since ${\displaystyle \rho _{k}}$ form a linearly independent basis,

${\displaystyle \displaystyle c_{jk}=\sum _{m,n}\chi _{m,n}B_{m,n,j,k}}$.

Inverting ${\displaystyle B}$ gives ${\displaystyle \chi }$:

${\displaystyle \chi _{m,n}=\sum _{j,k}B_{m,n,j,k}^{-1}c_{jk}}$.