format: text

text:
# Quantum Tomography

### Acknowledgments

The writing of the present Review has been co-sponsored by: the Italian Ministero dell’Istruzione, dell’Universita’ e della Ricerca (MIUR) under the Cofinanziamento 2002 Entanglement assisted high precision measurements, the Istituto Nazionale di Fisica della Materia under the project PRA-2002-CLON, and by the European Community programs ATESIT (Contract No. IST-2000-29681) and EQUIP (Contract No. IST-1999-11053). G. M. D. acknowledges partial support by the Department of Defense Multidisciplinary University Research Initiative (MURI) program administered by the Army Research Office under Grant No. DAAD19-00-1-0177. M. G. A. P. is research fellow at Collegio Alessandro Volta.

## Chapter 1 Introduction

The state of a physical system is the mathematical description that provides a complete information on the system. Its knowledge is equivalent to know the result of any possible measurement on the system. In Classical Mechanics it is always possible, at least in principle, to devise a procedure made of multiple measurements which fully recovers the state of the system. In Quantum Mechanics, on the contrary, this is not possible, due to the fundamental limitations related to the Heisenberg uncertainty principle  and the no-cloning theorem . In fact, on one hand one cannot perform an arbitrary sequence of measurements on a single system without inducing on it a back-action of some sort. On the other hand, the no-cloning theorem forbids to create a perfect copy of the system without already knowing its state in advance. Thus, there is no way out, not even in principle, to infer the quantum state of a single system without having some prior knowledge on it .

It is possible to estimate the unknown quantum state of a system when many identical copies are available in the same state, so that a different measurement can be performed on each copy. A procedure of such kind is called quantum tomography. The problem of finding a procedure to determine the state of a system from multiple copies was first addressed in 1957 by Fano , who called quorum a set of observables sufficient for a complete determination of the density matrix. However, since for a particle it is difficult to devise concretely measurable observables other than position, momentum and energy, the fundamental problem of measuring the quantum state has remained at the level of mere speculation up to almost ten years ago, when the issue finally entered the realm of experimental physics with the pioneering experiments by Raymer’s group  in the domain of quantum optics. In quantum optics, in fact, using a balanced homodyne detector one has the unique opportunity of measuring all possible linear combinations of position and momentum of a harmonic oscillator, which here represents a single mode of the electromagnetic field.

The first technique to reconstruct the density matrix from homodyne measurements — so called homodyne tomography — originated from the observation by Vogel and Risken  that the collection of probability distributions achieved by homodyne detection is just the Radon transform of the Wigner function $`W`$. Therefore, as in classical imaging, by Radon transform inversion one can obtain $`W`$, and then from $`W`$ the matrix elements of the density operator. This first method, however, was affected by uncontrollable approximations, since arbitrary smoothing parameters are needed for the inverse Radon transform. In Ref.  the first exact technique was given for measuring experimentally the matrix elements of the density operator in the photon-number representation, by simply averaging functions of homodyne data. After that, the method was further simplified , and the feasibility for non-unit quantum efficiency of detectors—above some bounds—was established.

The exact homodyne method has been implemented experimentally to measure the photon statistics of a semiconductor laser , and the density matrix of a squeezed vacuum . The success of optical homodyne tomography has then stimulated the development of state-reconstruction procedures for atomic beams , the experimental determination of the vibrational state of a molecule , of an ensemble of helium atoms , and of a single ion in a Paul trap .

Through quantum tomography the state is perfectly recovered in the limit of infinite number of measurements, while in the practical finite-measurements case, one can always estimate the statistical error that affects the reconstruction. For infinite dimensions the propagation of statistical errors of the density matrix elements make them useless for estimating the ensemble average of unbounded operators, and a method for estimating the ensemble average of arbitrary observable of the field without using the density matrix elements has been derived . Further insight on the general method of state reconstruction has lead to generalize homodyne tomography to any number of modes , and then to extend the tomographic method from the harmonic oscillator to an arbitrary quantum system using group theory . A general data analysis method has been designed in order to unbias the estimation procedure from any known instrumental noise . Moreover, algorithms have been engineered to improve the statistical errors on a given sample of experimental data—the so-called adaptive tomography —and then max-likelihood strategies  have been used that improved dramatically statistical errors, however, at the expense of some bias in the infinite dimensional case, and of exponential complexity versus $`N`$ for the joint tomography of $`N`$ quantum systems. The latest technical developments  derive the general tomographic method from spanning sets of operators, the previous group theoretical approaches  being just a particular case of this general method, where the group representation is just a device to find suitable operator “orthogonality” and “completeness” relations in the linear algebra of operators. Finally, very recently, a method for tomographic estimation of the unknown quantum operation of a quantum device has been derived , which uses a single fixed input entangled state, which plays the role of all possible input states in quantum parallel on the tested device, making finally the method a true “quantum radiography” of the functioning of a device.

In this Review we will give a self-contained and unified derivation of the methods of quantum tomography, with examples of applications to different kinds of quantum systems, and with particular focus on quantum optics, where also some results from experiments are reexamined. The Review is organized as follows.

In Chapter 2 we introduce the generalized Wigner functions  and we provide the basic elements of detection theory in quantum optics, by giving the description of photodetection, homodyne detection, and heterodyne detection. As we will see, heterodyne detection also provides a method for estimating the ensemble average of polynomials in the field operators, however, it is unsuitable for the density matrix elements in the photon-number representation. The effect of non unit quantum efficiency is taken into account for all such detection schemes.

In Chapter 3 we give a brief history of quantum tomography, starting with the first proposal of Vogel and Risken  as the extension to the domain of quantum optics of the conventional tomographic imaging. As already mentioned, this method indirectly recovers the state of the system through the reconstruction of the Wigner function, and is affected by uncontrollable bias. The exact homodyne tomography method of Ref.  (successively simplified in Ref. ) is here presented on the basis of the general tomographic method of spanning sets of operators of Ref. . As another application of the general method, the tomography of spin systems  is provided from the group theoretical method of Refs. . In this chapter we include also further developments to improve the method, such as the deconvolution techniques of  to correct the effects of experimental noise by data processing, and the adaptive tomography  to reduce the statistical fluctuations of tomographic estimators.

Chapter 4 is devoted to the evaluation from Ref.  of the expectation value of arbitrary operators of a single-mode radiation field via homodyne tomography. Here we also report from Ref.  the estimation of the added noise with respect to the perfect measurement of field observables, for some relevant observables, along with a comparison with the noise that would have been obtained using heterodyne detection.

The generalization of Ref.  of homodyne tomography to many modes of radiation is reviewed in Chapter 5, where it is shown how tomography of a multimode field can be performed by using only a single local oscillator with a tunable field mode. Some results of Monte Carlo simulations from Ref.  are also shown for the state that describes light from parametric downconversion.

Chapter 6 reviews some applications of quantum homodyne tomography to perform fundamental test of quantum mechanics. The first is the proposal of Ref.  to measure the nonclassicality of radiation field. The second is the scheme of Ref.  to test the state reduction rule using light from parametric downconversion. Finally, we review some experimental results about tomography of coherent signals with applications to the estimation of losses introduced by simple optical components .

Chapter 7 reviews the tomographic method of Ref.  to reconstruct the quantum operation of a device, such as an amplifier or a measuring device, using a single fixed input entangled state, which plays the role of all possible input states in a quantum parallel fashion.

Chapter 8 is devoted to the reconstruction technique of Ref.  based on the maximum likelihood principle. As mentioned, for infinite dimensions this method is necessarily biased, however, it is more suited to the estimation of a finite number of parameters, as proposed in Ref. , or to the state determination in the presence of very low number of experimental data . Unfortunately, the algorithm of this method has exponential complexity versus the number of quantum systems for a joint tomography of many systems.

Finally, in Chapter 9 we briefly review Ref. , showing how quantum tomography could be profitably used as a tool for reconstruction and compression in classical imaging.

## Chapter 2 Wigner functions and elements of detection theory

In this chapter we review some simple formulas from Ref.  that connect the generalized Wigner functions for $`s`$-ordering with the density matrix, and vice-versa. These formulas prove very useful for quantum mechanical applications as, for example, for connecting master equations with Fokker-Planck equations, or for evaluating the quantum state from Monte Carlo simulations of Fokker-Planck equations, and finally for studying positivity of the generalized Wigner functions in the complex plane. Moreover, as we will show in Chapter 3, the first proposal of quantum state reconstruction  used the Wigner function as an intermediate step.

In the second part of the chapter, we evaluate the probability distribution of the photocurrent of photodetectors, balanced homodyne detectors, and heterodyne detectors. We show that under suitable limits the respective photocurrents provide the measurement of the photon number distribution, of the quadrature, and of the complex amplitude of a single mode of the electromagnetic field. When the effect of non-unit quantum efficiency is taken into account an additional noise affects the measurement, giving a Bernoulli convolution for photo-detection, and a Gaussian convolution for homodyne and heterodyne detection. Extensive use of the results in this chapter will be made in the next chapters devoted to quantum homodyne tomography.

### 2.1 Wigner functions

Since Wigner’s pioneering work , generalized phase-space techniques have proved very useful in various branches of physics . As a method to express the density operator in terms of c-number functions, the Wigner functions often lead to considerable simplification of the quantum equations of motion, as for example, for transforming master equations in operator form into more amenable Fokker-Planck differential equations (see, for example, Ref. ). Using the Wigner function one can express quantum-mechanical expectation values in form of averages over the complex plane (the classical phase-space), the Wigner function playing the role of a c-number quasi-probability distribution, which generally can also have negative values. More precisely, the original Wigner function allows to easily evaluate expectations of symmetrically ordered products of the field operators, corresponding to the Weyl’s quantization procedure . However, with a slight change of the original definition, one defines generalized $`s`$-ordered Wigner function $`W_s(\alpha ,\alpha ^{})`$, as follows

$`W_s(\alpha ,\alpha ^{})={\displaystyle _{}}{\displaystyle \frac{d^2\lambda }{\pi ^2}}e^{\alpha \lambda ^{}\alpha ^{}\lambda +\frac{s}{2}|\lambda |^2}\text{Tr}[D(\lambda )\rho ],`$ (2.1)
where $`\alpha ^{}`$ denotes the complex conjugate of $`\alpha `$, the integral is performed on the complex plane with measure $`d^2\lambda =d\text{Re}\lambda d\text{Im}\lambda `$, $`\rho `$ represents the density operator, and

$`D(\alpha )\mathrm{exp}(\alpha a^{}\alpha ^{}a)`$ (2.2)
denotes the displacement operator, where $`a`$ and $`a^{}`$ ($`[a,a^{}]=1`$) are the annihilation and creation operators of the field mode of interest. The Wigner functions in Eq. (2.1) allow to evaluate $`s`$-ordered expectation values of the field operators through the following relation

$`\text{Tr}[\text{:}(a^{})^na^m\text{:}_s\rho ]={\displaystyle _{}}d^2\alpha W_s(\alpha ,\alpha ^{})\alpha ^n\alpha ^m.`$ (2.3)
The particular cases $`s=1,0,1`$ correspond to anti-normal, symmetrical, and normal ordering, respectively. In these cases the generalized Wigner function $`W_s(\alpha ,\alpha ^{})`$ are usually denoted by the following symbols and names

$`\begin{array}{cc}\frac{1}{\pi }Q(\alpha ,\alpha ^{})\hfill & \text{ for }s=1\text{ “}Q\text{ function”}\hfill \\ W(\alpha ,\alpha ^{})\hfill & \text{ for }s=0\text{ (usual Wigner function)}\hfill \\ P(\alpha ,\alpha ^{})\hfill & \text{ for }s=1\text{ “}P\text{ function”}\hfill \end{array}`$ (2.7)
For the normal ($`s=1`$) and anti-normal ($`s=1`$) orderings, the following simple relations with the density matrix are well known

$`Q(\alpha ,\alpha ^{})\alpha |\rho |\alpha ,`$ (2.8)
$`\rho ={\displaystyle _{}}d^2\alpha P(\alpha ,\alpha ^{})|\alpha \alpha |,`$ (2.9)
where $`|\alpha `$ denotes the customary coherent state $`|\alpha =D(\alpha )|0`$, $`|0`$ being the vacuum state of the field. Among the three particular representations (2.7), the $`Q`$ function is positively definite and infinitely differentiable (it actually represents the probability distribution for ideal joint measurements of position and momentum of the harmonic oscillator: see Sec. 2.4). On the other hand, the $`P`$ function is known to be possibly highly singular, and the only pure states for which it is positive are the coherent states . Finally, the usual Wigner function has the remarkable property of providing the probability distribution of the quadratures of the field in form of marginal distribution, namely

$`{\displaystyle _{\mathrm{}}^{\mathrm{}}}d\text{Im}\alpha W(\alpha e^{i\phi },\alpha ^{}e^{i\phi })={}_{\phi }{}^{}\text{Re}\alpha |\rho |\text{Re}\alpha _{\phi }^{},`$ (2.10)
where $`|x_\phi `$ denotes the (unnormalizable) eigenstate of the field quadrature

$`X_\phi ={\displaystyle \frac{a^{}e^{i\phi }+ae^{i\phi }}{2}}`$ (2.11)
with real eigenvalue $`x`$. Notice that any couple of quadratures $`X_\phi `$, $`X_{\phi +\pi /2}`$ is canonically conjugate, namely $`[X_\phi ,X_{\phi +\pi /2}]=i/2`$, and it is equivalent to position and momentum of a harmonic oscillator. Usually, negative values of the Wigner function are viewed as signature of a non-classical state, the most eloquent example being the Schrödinger-cat state , whose Wigner function is characterized by rapid oscillations around the origin of the complex plane. From Eq. (2.1) one can notice that all $`s`$-ordered Wigner functions are related to each other through Gaussian convolution

$`W_s(\alpha ,\alpha ^{})`$ $`=`$ $`{\displaystyle _{}}d^2\beta W_s^{}(\beta ,\beta ^{}){\displaystyle \frac{2}{\pi (s^{}s)}}\mathrm{exp}\left({\displaystyle \frac{2}{s^{}s}}|\alpha \beta |^2\right)`$ (2.12)
$`=`$ $`\mathrm{exp}\left({\displaystyle \frac{s^{}s}{2}}{\displaystyle \frac{^2}{\alpha \alpha ^{}}}\right)W_s^{}(\alpha ,\alpha ^{}),(s^{}>s).`$ (2.13)
Equation (2.12) shows the positivity of the generalized Wigner function for $`s<1`$, as a consequence of the positivity of the $`Q`$ function. From a qualitative point of view, the maximum value of $`s`$ keeping the generalized Wigner functions as positive can be considered as an indication of the classical nature of the physical state .

An equivalent expression for $`W_s(\alpha ,\alpha ^{})`$ can be derived as follows . Eq. (2.1) can be rewritten as

$`W_s(\alpha ,\alpha ^{})=\text{Tr}[\rho D(\alpha )\widehat{W}_sD^{}(\alpha )],`$ (2.14)
where

$`\widehat{W}_s={\displaystyle _{}}{\displaystyle \frac{d^2\lambda }{\pi ^2}}e^{\frac{s}{2}|\lambda |^2}D(\lambda ).`$ (2.15)
Through the customary Baker-Campbell-Hausdorff (BCH) formula

$`\mathrm{exp}A\mathrm{exp}B=\mathrm{exp}\left(A+B+{\displaystyle \frac{1}{2}}[A,B]\right),`$ (2.16)
which holds when $`[A,[A,B]]=[B,[A,B]]=0`$, one writes the displacement in normal order, and integrating on $`\mathrm{arg}(\lambda )`$ and $`|\lambda |`$ one obtains

$`\widehat{W}_s={\displaystyle \frac{2}{\pi (1s)}}{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{1}{n!}}\left({\displaystyle \frac{2}{s1}}\right)^na^na^n={\displaystyle \frac{2}{\pi (1s)}}\left({\displaystyle \frac{s+1}{s1}}\right)^{a^{}a},`$ (2.17)
where we used the normal-ordered forms

$`\text{:}(a^{}a)^n\text{:}=(a^{})^na^n=a^{}a(a^{}a1)\mathrm{}(a^{}an+1),`$ (2.18)
and the identity

$`\text{:}e^{xa^{}a}\text{:}={\displaystyle \underset{l=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{(x)^l}{l!}}(a^{})^la^l=(1x)^{a^{}a}.`$ (2.19)
The density matrix can be recovered from the generalized Wigner functions using the following expression

$`\rho ={\displaystyle \frac{2}{1+s}}{\displaystyle _{}}d^2\alpha W_s(\alpha ,\alpha ^{})e^{\frac{2}{1+s}|\alpha |^2}e^{\frac{2\alpha }{1+s}a^{}}\left({\displaystyle \frac{s1}{s+1}}\right)^{a^{}a}e^{\frac{2\alpha ^{}}{1+s}a}.`$ (2.20)
For the proof of Eq. (2.20) the reader is referred to Ref. . In particular, for $`s=0`$ one has the inverse of the Glauber formula

$`\rho =2{\displaystyle _{}}d^2\alpha W(\alpha ,\alpha ^{})D(2\alpha )()^{a^{}a},`$ (2.21)
whereas for $`s=1`$ one recovers Eq. (2.9) that defines the $`P`$ function.

### 2.2 Photodetection

Light is revealed by exploiting its interaction with atoms/molecules or electrons in a solid, and, essentially, each photon ionizes a single atom or promotes an electron to a conduction band, and the resulting charge is then amplified to produce a measurable pulse. In practice, however, available photodetectors are not ideally counting all photons, and their performances is limited by a non-unit quantum efficiency $`\zeta `$, namely only a fraction $`\zeta `$ of the incoming photons lead to an electric signal, and ultimately to a count: some photons are either reflected from the surface of the detector, or are absorbed without being transformed into electric pulses.

Let us consider a light beam entering a photodetector of quantum efficiency $`\zeta `$, i.e. a detector that transforms just a fraction $`\zeta `$ of the incoming light pulse into electric signal. If the detector is small with respect to the coherence length of radiation and its window is open for a time interval $`T`$, then the Poissonian process of counting gives a probability $`p(m;T)`$ of revealing $`m`$ photons that writes

$`p(m;T)=\text{Tr}\left[\rho \text{:}{\displaystyle \frac{[\zeta I(T)T]^m}{m!}}\mathrm{exp}[\zeta I(T)T]\text{:}\right],`$ (2.22)
where $`\rho `$ is the quantum state of light, : : denotes the normal ordering of field operators, and $`I(T)`$ is the beam intensity

$`I(T)={\displaystyle \frac{2ϵ_0c}{T}}{\displaystyle _0^T}\text{E}^{()}(\text{r},t)\text{E}^{(+)}(\text{r},t)𝑑t,`$ (2.23)
given in terms of the positive (negative) frequency part of the electric field operator $`\text{E}^{(+)}(\text{r},t)`$ ($`\text{E}^{()}(\text{r},t)`$). The quantity $`p(t)=\zeta \text{Tr}\left[\rho I(T)\right]`$ equals the probability of a single count during the time interval $`(t,t+dt)`$. Let us now focus our attention to the case of the radiation field excited in a stationary state of a single mode at frequency $`\omega `$. Eq. (2.22) can be rewritten as

$`p_\eta (m)=\text{Tr}\left[\rho \text{:}{\displaystyle \frac{(\eta a^{}a)^m}{m!}}\mathrm{exp}(\eta a^{}a)\text{:}\right],`$ (2.24)
where the parameter $`\eta =\zeta c\mathrm{}\omega /V`$ denotes the overall quantum efficiency of the photodetector. Using Eqs. (2.18) and (2.19) one obtains

$`p_\eta (m)={\displaystyle \underset{n=m}{\overset{\mathrm{}}{}}}\rho _{nn}\left(\begin{array}{c}n\\ m\end{array}\right)\eta ^m(1\eta )^{nm},`$ (2.27)
where

$`\rho _{nn}n|\rho |n=p_{\eta =1}(n).`$ (2.28)
Hence, for unit quantum efficiency a photodetector measures the photon number distribution of the state, whereas for non unit quantum efficiency the output distribution of counts is given by a Bernoulli convolution of the ideal distribution.

The effects of non unit quantum efficiency on the statistics of a photodetector, i.e. Eq. (2.27) for the output distribution, can be also described by means of a simple model in which the realistic photodetector is replaced with an ideal photodetector preceded by a beam splitter of transmissivity $`\tau \eta `$. The reflected mode is absorbed, whereas the transmitted mode is photo-detected with unit quantum efficiency. In order to obtain the probability of measuring $`m`$ clicks, notice that, apart from trivial phase changes, a beam splitter of transmissivity $`\tau `$ affects the unitary transformation of fields

$`\left(\begin{array}{c}c\\ d\end{array}\right)U_\tau ^{}\left(\begin{array}{c}a\\ b\end{array}\right)U_\tau =\left(\begin{array}{cc}\sqrt{\tau }& \sqrt{1\tau }\\ \sqrt{1\tau }& \sqrt{\tau }\end{array}\right)\left(\begin{array}{c}a\\ b\end{array}\right),`$ (2.37)
where all field modes are considered at the same frequency. Hence, the output mode $`c`$ hitting the detector is given by the linear combination

$`c=\sqrt{\tau }a\sqrt{1\tau }b,`$ (2.38)
and the probability of counts reads

$`p_\tau (m)`$ $`=`$ $`\text{Tr}\left[U_\tau \rho |00|U_\tau ^{}|mm|1\right]`$ (2.41)
$`=`$ $`{\displaystyle \underset{n=m}{\overset{\mathrm{}}{}}}\rho _{nn}\left(\begin{array}{c}n\\ m\end{array}\right)(1\tau )^{nm}\tau ^m.`$
Eq. (2.41) reproduces the probability distribution of Eq. (2.27) with $`\tau =\eta `$. We conclude that a photo-detector of quantum efficiency $`\eta `$ is equivalent to a perfect photo-detector preceded by a beam splitter of transmissivity $`\eta `$ which accounts for the overall losses of the detection process.

### 2.3 Balanced homodyne detection

The balanced homodyne detector provides the measurement of the quadrature of the field $`X_\phi `$ in Eq. (2.11). It was proposed by Yuen and Chan , and subsequently demonstrated by Abbas, Chan and Yee .

The scheme of a balanced homodyne detector is depicted in Fig. 2.1. The signal mode $`a`$ interferes with a strong laser beam mode $`b`$ in a balanced 50/50 beam splitter. The mode $`b`$ is the so-called the local oscillator (LO) mode of the detector. It operates at the same frequency of $`a`$, and is excited by the laser in a strong coherent state $`|z`$. Since in all experiments that use homodyne detectors the signal and the LO beams are generated by a common source, we assume that they have a fixed phase relation. In this case the LO phase provides a reference for the quadrature measurement, namely we identify the phase of the LO with the phase difference between the two modes. As we will see, by tuning $`\phi =\mathrm{arg}z`$ we can measure the quadrature $`X_\phi `$ at different phases.

After the beam splitter the two modes are detected by two identical photodetectors (usually linear avalanche photodiodes), and finally the difference of photocurrents at zero frequency is electronically processed and rescaled by $`2|z|`$. According to Eqs. (2.37), the modes at the output of the $`50/50`$ beam splitter ($`\tau =1/2`$) write

$`c={\displaystyle \frac{ab}{\sqrt{2}}},d={\displaystyle \frac{a+b}{\sqrt{2}}},`$ (2.42)
hence the difference of photocurrents is given by the following operator

$`I={\displaystyle \frac{d^{}dc^{}c}{2|z|}}={\displaystyle \frac{a^{}b+b^{}a}{2|z|}}.`$ (2.43)
Let us now proceed to evaluate the probability distribution of the output photocurrent $`I`$ for a generic state $`\rho `$ of the signal mode $`a`$. In the following treatment we will follow Refs. .

Let us consider the moments generating function of the photocurrent $`I`$

$`\chi (\lambda )=\text{Tr}\left[\rho |zz|e^{i\lambda I}\right],`$ (2.44)
which provides the probability distribution of $`I`$ as the Fourier transform

$`P(I)={\displaystyle _{\mathrm{}}^+\mathrm{}}{\displaystyle \frac{d\lambda }{2\pi }}e^{i\lambda I}\chi (\lambda ).`$ (2.45)
Using the BCH formula  for the $`SU(2)`$ group, namely

$`\mathrm{exp}\left(\xi ab^{}\xi ^{}a^{}b\right)=e^{\zeta b^{}a}\left(1+|\zeta |^2\right)^{\frac{1}{2}\left(b^{}ba^{}a\right)}e^{\zeta ^{}a^{}b},\zeta ={\displaystyle \frac{\xi }{|\xi |}}\mathrm{tan}|\xi |,`$ (2.46)
one can write the exponential in Eq. (2.44) in normal-ordered form with respect to mode $`b`$ as follows

$`\chi (\lambda )=e^{i\mathrm{tan}\left(\frac{\lambda }{2|z|}\right)b^{}a}\left[\mathrm{cos}\left({\displaystyle \frac{\lambda }{2|z|}}\right)\right]^{a^{}ab^{}b}e^{i\mathrm{tan}\left(\frac{\lambda }{2|z|}\right)a^{}b}_{ab}.`$ (2.47)
Since mode $`b`$ is in a coherent state $`|z`$ the partial trace over $`b`$ can be evaluated as follows

$`\chi (\lambda )`$ $`=`$ $`e^{i\mathrm{tan}\left(\frac{\lambda }{2|z|}\right)z^{}a}\left[\mathrm{cos}\left({\displaystyle \frac{\lambda }{2|z|}}\right)\right]^{a^{}a}e^{i\mathrm{tan}\left(\frac{\lambda }{2|z|}\right)za^{}}_a`$ (2.48)
$`\times `$ $`z\left|\left[\mathrm{cos}\left({\displaystyle \frac{\lambda }{2|z|}}\right)\right]^{b^{}b}\right|z.`$
Using now Eq. (2.16), one can rewrite Eq. (2.48) in normal order with respect to $`a`$, namely

$`\chi (\lambda )=e^{iz\mathrm{sin}\left(\frac{\lambda }{2|z|}\right)a^{}}\mathrm{exp}\left[2\mathrm{sin}^2\left({\displaystyle \frac{\lambda }{4|z|}}\right)(a^{}a+|z|^2)\right]e^{iz^{}\mathrm{sin}\left(\frac{\lambda }{2|z|}\right)a}_a,`$ (2.49)
In the strong-LO limit $`z\mathrm{}`$, only the lowest order terms in $`\lambda /|z|`$ are retained, $`a^{}a`$ is neglected with respect to $`|z|^2`$, and Eq. (2.49) simplifies as follows

$`\underset{z\mathrm{}}{lim}\chi (\lambda )=e^{i\frac{\lambda }{2}e^{i\phi }a^{}}\mathrm{exp}\left[{\displaystyle \frac{\lambda ^2}{8}}\right]e^{i\frac{\lambda }{2}e^{i\phi }a}_a=\mathrm{exp}[i\lambda X_\phi ]_a,`$ (2.50)
where $`\phi =\text{arg}z`$. The generating function in Eq. (2.50) is then equivalent to the POVM

$`\mathrm{\Pi }(x)={\displaystyle _{\mathrm{}}^+\mathrm{}}{\displaystyle \frac{d\lambda }{2\pi }}\mathrm{exp}[i\lambda (X_\phi x)]=\delta (X_\phi x)|x_\phi _\phi x|,`$ (2.51)
namely the projector on the eigenstate of the quadrature $`X_\phi `$ with eigenvalue $`x`$. In conclusion, the balanced homodyne detector achieves the ideal measurement of the quadrature $`X_\phi `$ in the strong LO limit. In this limit, the probability distribution of the output photocurrent $`I`$ approaches exactly the probability distribution $`p(x,\phi )={}_{\phi }{}^{}x|\rho |x_{\phi }^{}`$ of the quadrature $`X_\phi `$, and this for any state $`\rho `$ of the signal mode $`a`$.

It is easy to take into account non-unit quantum efficiency at detectors. According to Eq. (2.38) one has the replacements

$`c`$ $``$ $`\sqrt{\eta }c\sqrt{1\eta }u,u,v\text{ vacuum modes}`$ (2.52)
$`d`$ $``$ $`\sqrt{\eta }d\sqrt{1\eta }v,`$ (2.53)
and now the output current is rescaled by $`2|z|\eta `$, namely

$`I_\eta {\displaystyle \frac{1}{2|z|}}\{[a+\sqrt{{\displaystyle \frac{1\eta }{2\eta }}}(u+v)]b^{}+h.c\},`$ (2.54)
where only terms containing the strong LO mode $`b`$ are retained. The POVM is then obtained by replacing

$`X_\phi X_\phi +\sqrt{{\displaystyle \frac{1\eta }{2\eta }}}(u_\phi +v_\phi )`$ (2.55)
in Eq. (2.51), with $`w_\phi =(w^{}e^{i\phi }+we^{i\phi })/2`$, $`w=u,v`$, and tracing the vacuum modes $`u`$ and $`v`$. One then obtains

$`\mathrm{\Pi }_\eta (x)`$ $`=`$ $`{\displaystyle _{\mathrm{}}^+\mathrm{}}{\displaystyle \frac{d\lambda }{2\pi }}e^{i\lambda (X_\phi x)}|0|e^{i\lambda \sqrt{\frac{1\eta }{2\eta }}u_\phi }|0|^2={\displaystyle _{\mathrm{}}^+\mathrm{}}{\displaystyle \frac{d\lambda }{2\pi }}e^{i\lambda (X_\phi x)}e^{\lambda ^2\frac{1\eta }{8\eta }}`$ (2.56)
$`=`$ $`{\displaystyle \frac{1}{\sqrt{2\pi \mathrm{\Delta }_\eta ^2}}}\mathrm{exp}\left[{\displaystyle \frac{(xX_\phi )^2}{2\mathrm{\Delta }_\eta ^2}}\right]`$
$`=`$ $`{\displaystyle \frac{1}{\sqrt{2\pi \mathrm{\Delta }_\eta ^2}}}{\displaystyle _{\mathrm{}}^+\mathrm{}}dx^{}e^{\frac{1}{2\mathrm{\Delta }_\eta ^2}(xx^{})^2}|x^{}_\phi _\phi x^{}|,`$
where

$`\mathrm{\Delta }_\eta ^2={\displaystyle \frac{1\eta }{4\eta }}.`$ (2.57)
Thus the POVM, and in turn the probability distribution of the output photocurrent, are just the Gaussian convolution of the ideal ones with rms $`\mathrm{\Delta }_\eta =\sqrt{(1\eta )/(4\eta )}`$.

### 2.4 Heterodyne detection

Heterodyne detection allows to perform the joint measurement of two conjugated quadratures of the field . The scheme of the heterodyne detector is depicted in Fig. 2.2.

A strong local oscillator at frequency $`\omega `$ in a coherent state $`|\alpha `$ hits a beam splitter with transmissivity $`\tau 1`$, and with the coherent amplitude $`\alpha `$ such that $`\gamma |\alpha |\sqrt{\tau (1\tau )}`$ is kept constant. If the output photocurrent is sampled at the intermediate frequency $`\omega _{IF}`$, just the field modes $`a`$ and $`b`$ at frequency $`\omega \pm \omega _{IF}`$ are selected by the detector. Modes $`a`$ and $`b`$ are usually referred to as signal band and image band modes, respectively. In the strong LO limit, upon tracing the LO mode, the output photocurrent $`I(\omega _{IF})`$ rescaled by $`\gamma `$ is equivalent to the complex operator

$`Z={\displaystyle \frac{I(\omega _{IF})}{\gamma }}=ab^{},`$ (2.58)
where the arbitrary phases of modes have been suitably chosen. The heterodyne photocurrent $`Z`$ is a normal operator, equivalent to a couple of commuting selfadjoint operators

$`Z=\text{Re}Z+i\text{Im}Z,[Z,Z^{}]=[\text{Re}Z,\text{Im}Z]=0.`$ (2.59)
The POVM of the detector is then given by the orthogonal eigenvectors of $`Z`$. It is here convenient to introduce the notation of Ref.  for vectors in the tensor product of Hilbert spaces $``$

$`|A={\displaystyle \underset{nm}{}}A_{nm}|n|m(AI)|I(IA^\tau )|I,`$ (2.60)
where $`A^\tau `$ denotes the transposed operator with respect to some pre-chosen orthonormal basis. Eq. (2.60) exploits the isomorphism between the Hilbert space of the Hilbert-Schmidt operators $`A,B\mathrm{𝖧𝖲}()`$ with scalar product $`A,B=\text{Tr}[A^{}B]`$, and the Hilbert space of bipartite vectors $`|A,|B`$, where one has $`A|BA,B`$.

Using the above notation it is easy to write the eigenvectors of $`Z`$ with eigenvalue $`z`$ as $`\frac{1}{\sqrt{\pi }}|D(z)`$. In fact one has

$`Z|D(z)`$ $`=`$ $`(ab^{})(D_a(z)I_b)|I=(D_a(z)I_b)(ab^{}+z){\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}|n|n`$ (2.61)
$`=`$ $`z(D_a(z)I_b)|I=z|D(z).`$
The orthogonality of such eigenvectors can be verified through the relation

$`D(z)|D(z^{})=\text{Tr}[D^{}(z)D(z^{})]=\pi \delta ^{(2)}(zz^{}),`$ (2.62)
where $`\delta ^{(2)}(\alpha )`$ denotes the Dirac delta function over the complex plane

$`\delta ^{(2)}(\alpha )={\displaystyle _{}}{\displaystyle \frac{d^2\gamma }{\pi ^2}}\mathrm{exp}(\gamma \alpha ^{}\gamma ^{}\alpha ).`$ (2.63)
In conventional heterodyne detection the image band mode is in the vacuum state, and one is just interested in measuring the field mode $`a`$. In this case we can evaluate the POVM upon tracing on mode $`b`$. One has

$`\mathrm{\Pi }(z,z^{})`$ $`=`$ $`{\displaystyle \frac{1}{\pi }}\text{Tr}_b[|D(z)D(z)|I_a|00|]`$ (2.64)
$`=`$ $`{\displaystyle \frac{1}{\pi }}D(z)|00|D^{}(z)={\displaystyle \frac{1}{\pi }}|zz|,`$
namely one obtain the projectors on coherent states. The coherent-state POVM provides the optimal joint measurement of conjugated quadratures of the field . In fact, heterodyne detection allows to measure the Q-function in Eq. (2.7). According to Eq. (2.3) then it provides the expectation value of anti-normal ordered field operator. For a state $`\rho `$ the expectation value of any quadrature $`X_\phi `$ is obtained as

$`X_\phi =\text{Tr}[\rho X_\phi ]={\displaystyle _{}}{\displaystyle \frac{d^2\alpha }{\pi }}\text{Re}(\alpha e^{i\phi })Q(\alpha ,\alpha ^{}).`$ (2.65)
The price to pay for jointly measuring noncommuting observables is an additional noise. The rms fluctuation is evaluated as follows

$`{\displaystyle _{}}{\displaystyle \frac{d^2\alpha }{\pi }}[\text{Re}(\alpha e^{i\phi })]^2Q(\alpha ,\alpha ^{})X_\phi ^2=\mathrm{\Delta }X_\phi ^2+{\displaystyle \frac{1}{4}},`$ (2.66)
where $`\mathrm{\Delta }X_\phi ^2`$ is the intrinsic noise, and the additional term is usually referred to as “the additional 3dB noise due to the joint measure” .

The effect of non-unit quantum efficiency can be taken into account in analogous way as in Sec. 2.3 for homodyne detection. The heterodyne photocurrent is rescaled by an additional factor $`\eta ^{1/2}`$, and vacuum modes $`u`$ and $`v`$ are introduced, thus giving

$`Z_\eta =ab^{}+\sqrt{{\displaystyle \frac{1\eta }{\eta }}}u\sqrt{{\displaystyle \frac{1\eta }{\eta }}}v^{}.`$ (2.67)
Upon tracing over modes $`u`$ and $`v`$, one obtain the POVM

$`\mathrm{\Pi }_\eta (z,z^{})`$ $`=`$ $`{\displaystyle _{}}{\displaystyle \frac{d^2\gamma }{\pi ^2}}{}_{u}{}^{}0|{}_{v}{}^{}0|e^{\gamma (Z_\eta ^{}z^{})\gamma ^{}(Z_\eta z)}|0_{u}^{}|0_{v}^{}`$
$`=`$ $`{\displaystyle _{}}{\displaystyle \frac{d^2\gamma }{\pi ^2}}e^{\gamma (Z^{}z^{})\gamma ^{}(Zz)}e^{\frac{1\eta }{\eta }|\gamma |^2}`$
$`=`$ $`{\displaystyle \frac{\eta }{\pi (1\eta )}}e^{\frac{\eta }{1\eta }|Zz|^2}={\displaystyle _{}}{\displaystyle \frac{d^2z^{}}{\pi \mathrm{\Delta }_\eta ^2}}e^{\frac{|z^{}z|^2}{\mathrm{\Delta }_\eta ^2}}|D(z^{})D(z^{})|.`$
The probability distribution is then a Gaussian convolution on the complex plane of the ideal probability with rms $`\mathrm{\Delta }_\eta ^2=(1\eta )/\eta `$.

Analogously, the coherent-state POVM for conventional heterodyne detection with vacuum image band mode is replaced with

$`\mathrm{\Pi }_\eta (z,z^{})={\displaystyle _{}}{\displaystyle \frac{d^2z^{}}{\pi \mathrm{\Delta }_\eta ^2}}e^{\frac{|z^{}z|^2}{\mathrm{\Delta }_\eta ^2}}|z^{}z^{}|.`$ (2.69)
From Eqs. (2.12) we can equivalently say that the heterodyne detection probability density is given by the generalized Wigner function $`W_s(\alpha ,\alpha ^{})`$, with $`s=1\frac{2}{\eta }`$. Notice that for $`\eta <1`$, the average of functions $`\alpha ^n\alpha ^m`$ is related to the expectation value of a different ordering of field operators. However, one has the relevant identity

$`\text{:}(a^{})^na^m\text{:}_s={\displaystyle \underset{k=0}{\overset{(n,m)}{}}}k!\left({\displaystyle \genfrac{}{}{0pt}{}{n}{k}}\right)\left({\displaystyle \genfrac{}{}{0pt}{}{m}{k}}\right)\left({\displaystyle \frac{ts}{2}}\right)^k\text{:}(a^{})^{nk}a^{mk}\text{:}_t,`$ (2.70)
where $`(n,m)=\mathrm{min}(n,m)`$, and then

$`{\displaystyle _{}}d^2\alpha W_{1\frac{2}{\eta }}(\alpha ,\alpha ^{})\alpha ^m\alpha ^n`$
$`={\displaystyle \underset{k=0}{\overset{(n,m)}{}}}k!\left({\displaystyle \genfrac{}{}{0pt}{}{n}{k}}\right)\left({\displaystyle \genfrac{}{}{0pt}{}{m}{k}}\right)\left({\displaystyle \frac{1\eta }{\eta }}\right)^ka^{mk}(a^{})^{nk}.`$ (2.71)
Notice that the measure of the Q-function (or any smoothed version for $`\eta <1`$) does not allow to recover the expectation value of any operator through an average over heterodyne outcomes. In fact, one needs the admissibility of anti-normal ordered expansion  and the convergence of the integral in Eq. (2.71). In particular, the matrix elements of the density operator cannot be recovered. For some operators in which heterodyne measurement is allowed, a comparison with quantum homodyne tomography will be given in Sec. 4.3.

Finally, it is worth mentioning that all results of this section are valid also for an image-band mode with the same frequency of the signal. In this case a measurement scheme based on multiport homodyne detection should be used .

## Chapter 3 General tomographic method

In this chapter we review the general tomographic method of spanning sets of operators of Ref. , and re-derive in this general framework the exact homodyne tomography method of Ref. . In the first section, we first give a brief history of quantum tomography, starting with the original proposal of Vogel and Risken , that extended the conventional tomographic imaging to the domain of quantum optics. Here we will briefly sketch the conventional imaging tomographic method, and show the analogy with the method of Ref. . The latter achieves the quantum state via the Wigner function, which in turn is obtained by inverse Radon transform of the homodyne probability distributions for varying phase with respect to the LO. As already mentioned, the Radon transform inversion is affected by uncontrollable bias: such limitations and the intrinsic unreliability of this method are thoroughly explained in the same section.

As opposite to the Radon transform method, the first exact method of Ref.  (and successively refined in Ref. ) allows the reconstruction of the density matrix $`\rho `$, bypassing the step of the Wigner function, and achieving the matrix elements of $`\rho `$—or the expectation of any arbitrary operator—by just averaging the pertaining estimators (also called Kernel functions or pattern functions), evaluated on the experimental homodyne data. This method will be re-derived in Subsec. 3.3.3, as a special case of the general tomographic method of Ref.  here reviewed in Sect. 3.3, where we introduce the concept of “quorum”, which is the complete set of observables whose measurement provides the expectation value of any desired operator. Here we also show how some “orthogonality” and “completeness” relations in the linear algebra of operators are sufficient to individuate a quorum. As another application of the general method, in Subsect. 3.3.5 the tomography of spin systems  is reviewed, which was originally derived from the group theoretical methods of Refs. . Another application is the quantum tomography of a free particle state, given in Subsect. 3.3.6.

In this same chapter, in Sec. 3.4 we include some further developments to improve the tomographic method, such as the deconvolution techniques of Ref.  to correct the imperfections of detectors and experimental apparatus with a suitable data processing, and the adaptive tomography of Ref.  to reduce the statistical fluctuations of tomographic estimators, by adapting the averaged estimators to the given sample of experimental data.

The other relevant topics of homodyning observables, multimode tomography, and tomography of quantum operations will be given a separate treatment in the following chapters of the Review.

### 3.1 Brief historical excursus

The problem of quantum state determination through repeated measurements on identically prepared systems was originally stated by Fano in 1957 , who first recognized the need of measuring more that two noncommuting observables to achieve such purpose. However, it was only with the proposal by Vogel and Risken  that quantum tomography was born. The first experiments, which already showed reconstructions of coherent and squeezed states were performed by Michael Raymer’s and his group at the University of Oregon . The main idea at the basis of the first proposal is that it is possible to extend to the quantum domain the algorithms that are conventionally used in medical imaging to recover two dimensional (mass) distributions from unidimensional projections in different directions. This first tomographic method, however, was unreliable for the reconstruction of an unknown quantum state, since arbitrary smoothing parameters were needed in the Radon-transform based imaging procedure. The first exact unbiased tomographic method was proposed in Ref. , and successively simplified in Ref. . Since then, the new exact method has been practically implemented in many experiments, such as the measurement of the photon statistics of a semiconductor laser , and the reconstruction of the density matrix of a squeezed vacuum . The success of optical homodyne tomography has then stimulated the development of state-reconstruction procedures in other quantum harmonic oscillator systems, such as for atomic beams , and the vibrational state of a molecule , of an ensemble of helium atoms , and of a single ion in a Paul trap .

After the original exact method, quantum tomography has been generalized to the estimation of arbitrary observable of the field , to any number of modes , and, finally, to arbitrary quantum systems via group theory , with further improvements such as noise deconvolution , adaptive tomographic methods , and the use of max-likelihood strategies , which has made possible to reduce dramatically the number of experimental data, up to a factor $`10^3÷10^5`$, with negligible bias for most practical cases of interest. Finally, more recently, a method for tomographic estimation of the unknown quantum operation of a quantum device has been proposed , where a fixed input entangled state plays the role of all input states in a sort of quantum parallel fashion. Moreover, as another manifestation of such a quantum parallelism, one can also estimate the ensemble average of all operators by measuring only one fixed ”universal” observable on an extended Hilbert space in a sort of quantum hologram . This latest development is based on the general tomographic method of Ref. , where the tomographic reconstruction is based on the existence of spanning sets of operators, of which the irreducible unitary group representations of the group-methods of Refs.  are just a special case.

### 3.2 Conventional tomographic imaging

In conventional medical tomography, one collects data in the form of marginal distributions of the mass function $`m(x,y)`$. In the complex plane the marginal $`r(x,\phi )`$ is a projection of the complex function $`m(x,y)`$ on the direction indicated by the angle $`\phi [0,\pi ]`$, namely

$`r(x,\phi )={\displaystyle _{\mathrm{}}^+\mathrm{}}{\displaystyle \frac{dy}{\pi }}m((x+iy)e^{i\phi },(xiy)e^{i\phi }).`$ (3.1)
The collection of marginals for different $`\phi `$ is called “Radon transform”. The tomography process essentially consists in the inversion of the Radon transform (3.1), in order to recover the mass function $`m(x,y)`$ from the marginals $`r(x,\phi )`$.

Here we derive inversion of Eq. (3.1). Consider the identity

$`m(\alpha ,\alpha ^{})={\displaystyle _{}}d^2\beta \delta ^{(2)}(\alpha \beta )m(\beta ,\beta ^{}),`$ (3.2)
where $`\delta ^{(2)}(\alpha )`$ denotes the Dirac delta function of Eq. (2.63), and $`m(\alpha ,\alpha ^{})m(x,y)`$, with $`\alpha =x+iy`$ and $`\alpha ^{}=xiy`$. It is convenient to rewrite Eq. (2.63) as follows

$`\delta ^{(2)}(\alpha )={\displaystyle _0^+\mathrm{}}{\displaystyle \frac{dk}{4}}k{\displaystyle _0^{2\pi }}{\displaystyle \frac{d\phi }{\pi ^2}}e^{ik\alpha _\phi }={\displaystyle _{\mathrm{}}^+\mathrm{}}{\displaystyle \frac{dk}{4}}|k|{\displaystyle _0^\pi }{\displaystyle \frac{d\phi }{\pi ^2}}e^{ik\alpha _\phi },`$ (3.3)
with $`\alpha _\phi \text{Re}(\alpha e^{i\phi })=\alpha _{\phi +\pi }`$. Then, from Eqs. (3.2) and (3.3) the inverse Radon transform is obtained as follows

$`m(x,y)={\displaystyle _0^\pi }{\displaystyle \frac{d\phi }{\pi }}{\displaystyle _{\mathrm{}}^+\mathrm{}}𝑑x^{}r(x^{},\phi ){\displaystyle _{\mathrm{}}^+\mathrm{}}{\displaystyle \frac{dk}{4}}|k|e^{ik(x^{}\alpha _\phi )}.`$ (3.4)
Eq. (3.4) is conventionally written as

$`m(x,y)={\displaystyle _0^\pi }{\displaystyle \frac{d\phi }{\pi }}{\displaystyle _{\mathrm{}}^+\mathrm{}}𝑑x^{}r(x^{},\phi )K(x^{}\alpha _\phi ),`$ (3.5)
where $`K(x)`$ is given by

$`K(x){\displaystyle _{\mathrm{}}^+\mathrm{}}{\displaystyle \frac{dk}{4}}|k|e^{ikx}={\displaystyle \frac{1}{2}}\text{Re}{\displaystyle _0^+\mathrm{}}𝑑kke^{ikx}={\displaystyle \frac{1}{2}}𝒫{\displaystyle \frac{1}{x^2}},`$ (3.6)
with $`𝒫`$ denoting the Cauchy principal value. Integrating by parts Eq. (3.5) one obtains the tomographic formula that is usually found in medical imaging, i.e.

$`m(x,y)={\displaystyle \frac{1}{2\pi }}{\displaystyle _0^\pi }𝑑\phi 𝒫{\displaystyle _{\mathrm{}}^+\mathrm{}}𝑑x^{}{\displaystyle \frac{1}{x^{}\alpha _\phi }}{\displaystyle \frac{}{x^{}}}r(x^{},\phi ),`$ (3.7)
which allows the reconstruction of the mass distribution $`m(x,y)`$ from its projections along different directions $`r(x,\phi )`$.

#### 3.2.1 Extension to the quantum domain

In the “quantum imaging” process the goal is to reconstruct a quantum state in the form of its Wigner function starting from its marginal probability distributions. As shown in Sec. 2.1, the Wigner function is a real normalized function that is in one-to-one correspondence with the state density operator $`\rho `$. As noticed in Eq. (2.10), the probability distributions of the quadrature operators $`X_\phi =(a^{}e^{i\phi }+ae^{i\phi })/2`$ are the marginal probabilities of the Wigner function for the state $`\rho `$. Thus, by applying the same procedure outlined in the previous subsection, Vogel and Risken  proposed a method to recover the Wigner function via an inverse Radon transform from the quadrature probability distributions $`p(x,\phi )`$, namely

$`W(x,y)={\displaystyle _0^\pi }{\displaystyle \frac{d\phi }{\pi }}{\displaystyle _{\mathrm{}}^+\mathrm{}}𝑑x^{}p(x^{},\phi ){\displaystyle _{\mathrm{}}^+\mathrm{}}{\displaystyle \frac{dk}{4}}|k|e^{ik(x^{}x\mathrm{cos}\phi y\mathrm{sin}\phi )}.`$ (3.8)
\[Surprisingly, in the original paper  the connection to the tomographic imaging method was never mentioned\]. As shown in Sec. 2.3 the experimental measurement of the quadratures of the field is obtained using the homodyne detector. The method proposed by Vogel and Risken, namely the inversion of the Radon transform, was the one which has been used in the first experiments .

This first method is, however, not reliable for the reconstruction of an unknown quantum state, due to the intrinsic unavoidable systematic error related to the fact that the integral on $`k`$ in Eq. (3.8) is unbounded. In fact, in order to evaluate the inverse Radon transform, one would need the analytical form of the marginal distribution of the quadrature $`p(x,\phi )`$, which, in turn, can only be obtained by collecting the experimental data into histograms, and thence “spline-ing” them. This, of course, is not an unbiased procedure since the degree of spline-ing, the width and the number of the histogram bins, and finally the number of different phases used to collect the experimental data sample introduce systematic errors if they are not set above some minimal values, which actually depend on the unknown quantum state that one wants to reconstruct. Typically, an over-spline-ing will wash-out the quantum features of the state, whereas, vice-versa, an under-spline-ing will create negative photon probabilities in the reconstruction (see Ref.  for details).

A new exact method was then proposed in Ref. , alternative to the Radon-transform technique. This approach, referred to as quantum homodyne tomography, allows to recover the quantum state of the field $`\rho `$—along with any ensemble average of arbitrary operators—by directly averaging functions of the homodyne data, abolishing the intermediate step of the Wigner function, which is the source of all systematic errors. Only statistical errors are present, and they can be reduced arbitrarily by collecting more experimental data. This exact method will be re-derived from the general tomographic theory in Sec. 3.3.3.

### 3.3 General method of quantum tomography

In this section the general method of quantum tomography is explained in detail. First, we give the basics of Monte Carlo integral theory which are needed to implement the tomographic algorithms in actual experiments and in numerical simulations. Then, we derive the formulas on which all schemes of state reconstruction are based.

#### 3.3.1 Basic statistics

The aim of quantum tomography is to estimate, for arbitrary quantum system, the mean value $`O`$ of a system operator $`O`$ using only the results of the measurements on a set of observables $`\{Q_\lambda ,\lambda \mathrm{\Lambda }\}`$, called “quorum”. The procedure by which this can be obtained needs the estimator or “Kernel function” $`[O](x,\lambda )`$ which is a function of the eigenvalues $`x`$ of the quorum operators. Integrating the estimator with the probability $`p(x,\lambda )`$ of having outcome $`x`$ when measuring $`Q_\lambda `$, the mean value of $`O`$ is obtained as follows

$`O={\displaystyle _\mathrm{\Lambda }}𝑑\lambda {\displaystyle 𝑑\mu _\lambda (x)p(x,\lambda )[O](x,\lambda )},`$ (3.9)
where the first integral is performed on the values of $`\lambda `$ that designate all quorum observables, and the second on the eigenvalues of the quorum observable $`Q_\lambda `$ determined by the $`\lambda `$ variable of the outer integral. For discrete set $`\mathrm{\Lambda }`$ and/or discrete spectrum of the quorum, both integrals in (3.9) can suitably replaced by sums.

The algorithm to estimate $`O`$ with Eq. (3.9) is the following. One chooses a quorum operator $`Q_\lambda `$ by drawing $`\lambda `$ with uniform probability in $`\mathrm{\Lambda }`$ and performs a measurement, obtaining the result $`x_i`$. By repeating the procedure $`N`$ times, one collects the set of experimental data $`\{(x_i,\lambda _i)`$, with $`i=1,\mathrm{},N\}`$, where $`\lambda _i`$ identifies the quorum observable used for the $`i`$th measurement, and $`x_i`$ its result. From the same set of data the mean value of any operator $`O`$ can be obtained. In fact, one evaluates the estimator of $`O`$ and the quorum $`Q_\lambda `$, and then samples the double integral of (3.9) using the limit

$`O=\underset{N\mathrm{}}{lim}{\displaystyle \frac{1}{N}}{\displaystyle \underset{i=1}{\overset{N}{}}}[O](x_i,\lambda _i).`$ (3.10)
Of course the finite sum

$`F_N={\displaystyle \frac{1}{N}}{\displaystyle \underset{i=1}{\overset{N}{}}}[O](x_i,\lambda _i).`$ (3.11)
gives an approximation of $`O`$. To estimate the error in the approximation one applies the central limit theorem that we recall here.

Central limit theorem. Consider $`N`$ statistically uncorrelated random variables $`\{z_i,`$ $`i=1,\mathrm{},N\}`$, with mean values $`\mu (z_i)`$, variances $`\sigma ^2(z_i)`$ and bounded third order moments. If the variances $`\sigma ^2(z_i)`$ are all of the same order then the statistical variable “average” $`y`$ defined as

$$y_N=\frac{1}{N}\underset{i=1}{\overset{N}{}}z_i$$
(3.12)
has mean and variance

$$\mu (y_N)=\frac{1}{N}\underset{i=1}{\overset{N}{}}\mu (z_i),\sigma ^2(y_N)=\frac{1}{N^2}\underset{i=1}{\overset{N}{}}\sigma ^2(z_i).$$
(3.13)
The distribution of $`y_N`$ approaches asymptotically a Gaussian for $`N\mathrm{}`$. In practical cases, the distribution of $`y`$ can be considered Gaussian already for $`N`$ as low as $`N10`$.

For our needs the hypotheses are met if the estimator $`[O](x_i,\lambda _i)`$ in Eq. (3.11) has limited moments up to the third order, since, even though $`x_i`$ have different probability densities depending on $`\lambda _i`$, nevertheless, since $`\lambda _i`$ is also random all $`z_i`$ here given by

$$z_i=[O](x_i,\lambda _i)$$
(3.14)
have common mean

$`\mu (z_i)=O`$ (3.15)
and variance

$`\sigma ^2(z_i)={\displaystyle _\mathrm{\Lambda }}𝑑\lambda {\displaystyle 𝑑\mu _\lambda (x)p(x,\lambda )^2[O](x,\lambda )}O^2.`$ (3.16)
Using the central limit theorem, we can conclude that the experimental average $`yF_N`$ in Eq. (3.11) is a statistical variable distributed as a Gaussian with mean value $`\mu (y_N)\mu (z_i)`$ and variance $`\sigma ^2(y_N)\frac{1}{N}\sigma ^2(z_i)`$. Then the tomographic estimation converges with statistical error that decreases as $`N^{1/2}`$. A statistically precise estimate of the confidence interval is given by

$$ϵ_N=\sqrt{\frac{_{i=1}^N[z_iy_N]^2}{N(N1)}},$$
(3.17)
with $`z_i`$ given by Eq. (3.14) and $`y_N`$ by Eq. (3.12). In order to test that the confidence intervals are estimated correctly, one can check that the $`F_N`$ distribution is actually Gaussian. This can be done by comparing the histogram of the block data with a Gaussian, or by using the $`\chi ^2`$ test.

#### 3.3.2 Characterization of the quorum

Different estimation techniques have been proposed tailored to different quantum systems, such as the radiation field , trapped ions and molecular vibrational states , spin systems , etc. All the known quantum estimation techniques can be embodied in the following approach.

The tomographic reconstruction of an operator $`O`$ is possible when there exists a resolution of the form

$`O={\displaystyle _\mathrm{\Lambda }}𝑑\lambda \text{Tr}\left[OB^{}(\lambda )\right]C(\lambda ),`$ (3.18)
where $`\lambda `$ is a (possibly multidimensional) parameter living on a (discrete or continuous) manifold $`\mathrm{\Lambda }`$. The only hypothesis in (3.18) is the existence of the trace. If, for example, $`O`$ is a trace–class operator, then we do not need to require $`B(\lambda )`$ to be of Hilbert-Schmidt class, since it is sufficient to require $`B(\lambda )`$ bounded. The operators $`C(\lambda )`$ are functions of the quorum of observables measured for the reconstruction, whereas the operators $`B(\lambda )`$ form the dual basis of the set $`C(\lambda )`$. The term

$`[O](\lambda )=\text{Tr}\left[OB^{}(\lambda )\right]C(\lambda )`$ (3.19)
represents the quantum estimator for the operator $`O`$. The expectation value of $`O`$ is given by the ensemble average

$`O\text{Tr}\left[O\rho \right]={\displaystyle _\mathrm{\Lambda }}𝑑\lambda \text{Tr}\left[OB^{}(\lambda )\right]\text{Tr}\left[C(\lambda )\rho \right]{\displaystyle _\mathrm{\Lambda }}𝑑\lambda [O](\lambda ),`$ (3.20)
where $`\rho `$ is the density matrix of the quantum system under investigation. Notice that the quantity $`\text{Tr}\left[C(\lambda )\rho \right]`$ depends only on the quantum state, and it is related to the probability distribution of the measurement outcomes, whereas the term $`\text{Tr}\left[OB^{}(\lambda )\right]`$ depends only on the quantity to be measured. In particular, the tomography of the quantum state of a system corresponds to writing Eq. (3.18) for the operators $`O=|kn|`$, $`\{|n\}`$ being a given Hilbert space basis. For a given system, the existence of a set of operators $`C(\lambda )`$, together with its dual basis $`B(\lambda )`$ allows universal quantum estimation, i. e. the reconstruction of any operator.

We now give two characterizations of the sets $`B(\lambda )`$ and $`C(\lambda )`$ that are necessary and sufficient conditions for writing Eq. (3.18).

Condition 1: bi-orthogonality

Let us consider a complete orthonormal basis of vectors $`|n(n=0,1,\mathrm{})`$. Formula (3.18) is equivalent to the bi-orthogonality condition

$`{\displaystyle _\mathrm{\Lambda }}𝑑\lambda q|B^{}(\lambda )|pm|C(\lambda )|l=\delta _{mp}\delta _{lq},`$ (3.21)
where $`\delta _{ij}`$ is the Kronecker delta. Eq. (3.21) can be straightforwardly generalized to a continuous basis.

Condition 2: completeness

If the set of operators $`C(\lambda )`$ is complete, namely if any operator can be written as a linear combination of the $`C(\lambda )`$ as

$`O={\displaystyle _\mathrm{\Lambda }}𝑑\lambda a(\lambda )C(\lambda ),`$ (3.22)
then Eq. (3.18) is also equivalent to the trace condition

$`\text{Tr}\left[B^{}(\lambda )C(\mu )\right]=\delta (\lambda ,\mu ),`$ (3.23)
where $`\delta (\lambda ,\mu )`$ is a reproducing kernel for the set $`B(\lambda )`$, namely it is a function or a tempered distribution which satisfies

$`{\displaystyle _\mathrm{\Lambda }}𝑑\lambda B(\lambda )\delta (\lambda ,\mu )=B(\mu ).`$ (3.24)
An analogous identity holds for the set of $`C(\lambda )`$

$`{\displaystyle _\mathrm{\Lambda }}𝑑\lambda C(\lambda )\delta (\lambda ,\mu )=C(\mu ).`$ (3.25)
The proofs are straightforward. The completeness condition on the operators $`C(\lambda )`$ is essential for the equivalence of (3.18) and (3.23). A simple counterexample is provided by the set of projectors $`P(\lambda )=|\lambda \lambda |`$ over the eigenstates of a self-adjoint operator $`L`$. In fact, Eq. (3.23) is satisfied by $`C(\lambda )=B(\lambda )P(\lambda )`$. However, since they do not form a complete set in the sense of Eq. (3.22), it is not possible to express a generic operator in the form $`X=_\mathrm{\Lambda }𝑑\lambda \lambda |O|\lambda |\lambda \lambda |`$.

If either the set $`B(\lambda )`$ or the set $`C(\lambda )`$ satisfy the additional trace condition

$`\text{Tr}\left[B^{}(\mu )B(\lambda )\right]=\delta (\lambda ,\mu ),`$ (3.26)
$`\text{Tr}\left[C^{}(\mu )C(\lambda )\right]=\delta (\lambda ,\mu ),`$ (3.27)
then we have $`C(\lambda )=B(\lambda )`$ (notice that neither $`B(\lambda )`$ nor $`C(\lambda )`$ need to be unitary). In this case, Eq. (3.18) can be rewritten as

$`O={\displaystyle _\mathrm{\Lambda }}𝑑\lambda \text{Tr}\left[OC^{}(\lambda )\right]C(\lambda ).`$ (3.28)
A certain number of observables $`Q_\lambda `$ constitute a quorum when there are functions $`f_\lambda (Q_\lambda )=C(\lambda )`$ such that $`C(\lambda )`$ form an irreducible set. The quantum estimator for $`O`$ in Eq. (3.19) then writes as a function of the quorum operators

$`[O](\lambda )_\lambda [O](Q_\lambda ).`$ (3.29)
Notice that if a set of observables $`Q_\lambda `$ constitutes a quorum, than the set of projectors $`|q_\lambda _\lambda q|`$ over their eigenvectors provides a quorum too, with the measure $`d\lambda `$ in Eq. (3.18) including the measure $`d\mu _\lambda (q)`$. Notice also that, even once the quorum has been fixed, the unbiased estimator for an operator $`O`$ will not in general be unique, since there can exist functions $`𝒩(Q_\lambda )`$ that satisfies

$`{\displaystyle _\mathrm{\Lambda }}𝑑\lambda 𝒩(Q_\lambda )=0,`$ (3.30)
and that will be called ‘null estimators’. Two unbiased estimators that differ by a null estimator yield the same results when estimating the operator mean value. We will see in Sec. 3.4.2 how the null estimators can be used to reduce the statistical noise.

In terms of the quorum observables $`Q_\lambda `$ Eq. (3.20) rewrites

$`O`$ $`=`$ $`{\displaystyle _\mathrm{\Lambda }}𝑑\lambda \text{Tr}\left[OB^{}(\lambda )\right]\text{Tr}\left[f_\lambda (Q_\lambda )\rho \right]`$ (3.31)
$`=`$ $`{\displaystyle _\mathrm{\Lambda }}𝑑\lambda {\displaystyle 𝑑\mu _\lambda (q)p(q,\lambda )\text{Tr}\left[OB^{}(\lambda )\right]f_\lambda (q)},`$
where $`p(q,\lambda )={}_{\lambda }{}^{}q|\rho |q_{\lambda }^{}`$ is the probability density of getting the outcome $`q`$ from the measurement of $`Q_\lambda `$ on the state $`\rho `$. Eq. (3.31) is equivalent to the expression (3.9), with estimator

$`[O](q,\lambda )=\text{Tr}\left[OB^{}(\lambda )\right]f_\lambda (q).`$ (3.32)
Of course it is of interest to connect a quorum of observables to a resolution of the form (3.18), since only in this case there can be a feasible reconstruction scheme. If a resolution formula is written in terms of a set of selfadjoint operators, the set itself constitutes the desired quorum. However, in general a quorum of observables is functionally connected to the corresponding resolution formula. If the operators $`C(\lambda )`$ are unitary, then they can always be taken as the exponential map of a set of selfadjoint operators, which then are identified with our quorum $`Q_\lambda `$. The quantity $`\text{Tr}\left[C(\lambda )\rho \right]`$ is thus connected with the moment generating function of the set $`Q_\lambda `$, and hence to the probability density $`p(q,\lambda )`$ of the measurement outcomes, which play the role of the Radon transform in the quantum tomography of the harmonic oscillator. In general, the operators $`C(\lambda )`$ can be any function (neither self-adjoint nor unitary) of observables and, even more generally, they may be connected to POVMs rather than observables.

The dual set $`B(\lambda )`$ can be obtained from the set $`C(\lambda )`$ by solving Eq. (3.23). For finite quorums, this resorts to a matrix inversion. An alternative procedure uses the Gram-Schmidt orthogonalization procedure . No such a general procedure exists for a continuous spanning set. Many cases, however, satisfy conditions (3.26) and (3.27), and thus we can write $`B(\lambda )=C(\lambda )^{}`$.

#### 3.3.3 Quantum estimation for harmonic-oscillator systems

The harmonic oscillator models several systems of interest in quantum mechanics, as the vibrational states of molecules, the motion of an ion in a Paul trap, and a single mode radiation field. Different proposals have been suggested in order to reconstruct the quantum state of a harmonic system, which all fit the framework of the previous section, which is also useful for devising novel estimation techniques. Here, the basic resolution formula involves the set of displacement operators $`D(\alpha )=\mathrm{exp}(\alpha a^{}\alpha ^{}a)`$, which can be viewed as exponentials of the field-quadrature operators $`X_\phi =(a^{}e^{i\phi }+ae^{i\phi })/2`$. We have shown in Sec. 2.3 that for a single-mode radiation field $`X_\phi `$ is measured through homodyne detection. For the vibrational tomography of a molecule or a trapped ion $`X_\phi `$ corresponds to a time-evolved position or momentum. The set of displacement operators satisfies Eqs. (3.23) and (3.27), since

$`\text{Tr}[D(\alpha )D^{}(\beta )]=\pi \delta ^{(2)}(\alpha \beta ),`$ (3.33)
whereas Eq. (3.28) reduces to the Glauber formula

$`O={\displaystyle _{}}{\displaystyle \frac{d^2\alpha }{\pi }}\text{Tr}\left[OD^{}(\alpha )\right]D(\alpha ).`$ (3.34)
Changing to polar variables $`\alpha =(i/2)ke^{i\phi }`$, Eq. (3.34) becomes

$$O=_0^\pi \frac{d\phi }{\pi }_{\mathrm{}}^+\mathrm{}\frac{dk|k|}{4}\text{Tr}[Oe^{ikX_\phi }]e^{ikX_\phi },$$
(3.35)
which shows explicitly the dependence on the quorum $`X_\phi `$. Taking the ensemble average of both members and evaluating the trace over the set of eigenvectors of $`X_\phi `$, one obtains

$`O={\displaystyle _0^\pi }{\displaystyle \frac{d\phi }{\pi }}{\displaystyle _{\mathrm{}}^+\mathrm{}}𝑑xp(x,\phi )[O](x,\phi ),`$ (3.36)
where $`p(x;\phi )={}_{\phi }{}^{}x|\rho |x_{\phi }^{}`$ is the probability distribution of quadratures outcome. The estimator of the operator ensemble average $`O`$ is given by

$`[O](x,\phi )=\text{Tr}[OK(X_\phi x)],`$ (3.37)
where $`K(x)`$ is the same as in Eq. (3.6).

Eq. (3.36) is the basis of quantum homodyne tomography. Notice that even though $`K(x)`$ is unbounded, however, the matrix element $`\psi |K(X_\phi x)|\varphi `$ can be bounded, whence it can be used to sample the matrix element $`\psi |\rho |\varphi `$ of the state $`\rho `$, which, according to Sec. 3.3.1, is directly obtained by averaging the estimator (3.37) over homodyne experimental values. In fact, for bounded $`\psi |K(X_\phi x)|\varphi `$, the central limit theorem guarantees that

$`\psi |\rho |\varphi `$ $`=`$ $`{\displaystyle _0^\pi }{\displaystyle \frac{d\phi }{\pi }}{\displaystyle _{\mathrm{}}^+\mathrm{}}𝑑xp(x,\phi )\psi |K(X_\phi x)|\varphi `$ (3.38)
$`=`$ $`\underset{N\mathrm{}}{lim}{\displaystyle \frac{1}{N}}{\displaystyle \underset{n=0}{\overset{N}{}}}\psi |K(x_{\phi _n}x_n))|\varphi ,`$ (3.39)
where $`x_n`$ is the homodyne outcome measured at phase $`\phi _n`$ and distributed with probability $`p(x,\phi )`$. Systematic errors are eliminated by choosing randomly each phase $`\phi _n`$ at which homodyne measurement is performed. As shown in Sec. 3.3.1, for finite number of measurements $`N`$, the estimate (3.39) of the integral in Eq. (3.38) is Gaussian distributed around the true value $`\psi |\rho |\varphi `$, with statistical error decreasing as $`N^{1/2}`$. Notice that the measurability of the density operator matrix element depends only on the boundedness of the matrix element of the estimator, and that no adjustable parameters are needed in the procedure, which thus is unbiased.

The general procedure for noise deconvolution is presented in Sec. 3.4.1. However, we give here the main result for the density matrix reconstruction. As shown in Sec. 2.3, the effect of the efficiency in homodyne detectors is a Gaussian convolution of the ideal probability $`p(x,\phi )`$, as

$`p_\eta (x,\phi )=\sqrt{{\displaystyle \frac{2\eta }{\pi (1\eta )}}}{\displaystyle _{\mathrm{}}^+\mathrm{}}𝑑x^{}e^{\frac{2\eta }{1\eta }(xx^{})^2}p(x^{},\phi ).`$ (3.40)
The tomographic reconstruction procedure still holds upon replacing $`p(x,\phi )`$ with $`p_\eta (x,\phi )`$, so that

$`\rho ={\displaystyle _0^\pi }{\displaystyle \frac{d\phi }{\pi }}{\displaystyle _{\mathrm{}}^+\mathrm{}}𝑑xp_\eta (x,\phi )K_\eta (X_\phi x),`$ (3.41)
where now the estimator is

$`K_\eta (x)={\displaystyle \frac{1}{2}}\text{Re}{\displaystyle _0^+\mathrm{}}k𝑑ke^{\frac{1\eta }{8\eta }k^2+ikx}.`$ (3.42)
In fact, by taking the Fourier transform of both members of Eq. (3.40), one can easily check that

$`\rho `$ $`=`$ $`{\displaystyle _0^\pi }{\displaystyle \frac{d\phi }{\pi }}{\displaystyle _{\mathrm{}}^+\mathrm{}}𝑑xp_\eta (x,\phi )K_\eta (X_\phi x)`$ (3.43)
$`=`$ $`{\displaystyle _0^\pi }{\displaystyle \frac{d\phi }{\pi }}{\displaystyle _{\mathrm{}}^+\mathrm{}}𝑑xp(x,\phi )K(X_\phi x).`$
Notice that the anti-Gaussian in Eq. (3.42) causes a much slower convergence of the Monte Carlo integral (3.41): the statistical fluctuation will increase exponentially for decreasing detector efficiency $`\eta `$. In order to achieve good reconstructions with non-ideal detectors, then one has to collect a larger number of data.

It is clear from Eq. (3.39) that the measurability of the density matrix depends on the chosen representation and on the quantum efficiency of the detectors. For example, for the reconstruction of the density matrix in the Fock basis the estimators are given by

$`_\eta [|nn+d|](x,\phi )={\displaystyle _{\mathrm{}}^+\mathrm{}}{\displaystyle \frac{dk|k|}{4}}e^{\frac{1\eta }{8\eta }k^2ikx}n+d|e^{ikX_\phi }|n`$ (3.44)
$`=e^{id(\phi +\frac{\pi }{2})}\sqrt{{\displaystyle \frac{n!}{(n+d)!}}}{\displaystyle _{\mathrm{}}^+\mathrm{}}𝑑k|k|e^{\frac{12\eta }{2\eta }k^2i2kx}k^dL_n^d(k^2),`$
where $`L_n^d(x)`$ denotes the generalized Laguerre polynomials. Notice that the estimator is bounded only for $`\eta >1/2`$, and below the method would give unbounded statistical errors. However, this bound is well below the values that are reasonably achieved in the lab, where actual homodyne detectors have efficiency ranging between $`70\%`$ and $`90\%`$ . Moreover, a more efficient algorithm is available, that uses the factorization formulas that hold for $`\eta =1`$

$`[|nn+d|](x,\phi )=e^{id\phi }[4xu_n(x)v_{n+d}(x)`$ (3.45)
$`2\sqrt{n+1}u_{n+1}(x)v_{n+d}(x)2\sqrt{n+d+1}u_n(x)v_{n+d+1}(x)],`$
where $`u_j(x)`$ and $`v_j(x)`$ are the normalizable and unnormalizable eigenfunctions of the harmonic oscillator with eigenvalue $`j`$, respectively. The noise from quantum efficiency can be unbiased via the inversion of the Bernoulli convolution, which holds for $`\eta >1/2`$ .

The use of Eq. (3.36) to estimate arbitrary operators through homodyne tomography will be the subject of the following Chapter. Notice that Eq. (3.34) cannot be used for unbounded operators, however the estimators also for some unbounded operators will be derived in Sec. 4.1.

#### 3.3.4 Some generalizations

Using condition (3.23) one can see that the Glauber formula can be generalized to

$`O={\displaystyle _{}}{\displaystyle \frac{d^2\alpha }{\pi }}\text{Tr}\left[OF_1D(\alpha )F_2\right]F_2^1D^{}(\alpha )F_1^1,`$ (3.46)
where $`F_1`$ and $`F_2`$ are two generic invertible operators. By choosing $`F_1^{}=F_2=S(\zeta )`$, where $`S(\zeta )`$ is the squeezing operator

$`S(\zeta )=\mathrm{exp}\left[{\displaystyle \frac{1}{2}}\left(\zeta ^2a^2\zeta _{}^{}{}_{}{}^{2}a^2\right)\right],\zeta ,`$ (3.47)
we obtain the tomographic resolution

$`O={\displaystyle _0^\pi }{\displaystyle \frac{d\phi }{\pi }}{\displaystyle _{\mathrm{}}^+\mathrm{}}𝑑xp_\zeta (x,\phi )\text{Tr}\left[OK(X_{\phi \zeta }x)\right],`$ (3.48)
in terms of the probability distribution of the generalized squeezed quadrature operators

$`X_{\phi \zeta }=S^{}(\zeta )X_\phi S(\zeta )={\displaystyle \frac{1}{2}}\left[(\mu e^{i\phi }+\nu e^{i\phi })a^{}+(\mu e^{i\phi }+\nu ^{}e^{i\phi })a\right],`$ (3.49)
with $`\mu =\mathrm{cosh}|\zeta |`$ and $`\nu =\mathrm{sinh}|\zeta |\mathrm{exp}[2i\mathrm{arg}(\zeta )]`$. Such an estimation technique has been investigated in detail in Ref. .

A different estimation technique can be obtained by choosing in Eq. (3.46) $`F_1=I`$, the identity operator, and $`F_2=()^{a^{}a}`$, the parity operator. In this case one gets

$`O={\displaystyle _{}}{\displaystyle \frac{d^2\alpha }{\pi }}\text{Tr}\left[OD^{}(\alpha )()^{a^{}a}\right]()^{a^{}a}D(\alpha ).`$ (3.50)
Changing variable to $`\alpha =2\beta `$ and using the relation

$`()^{a^{}a}D(2\beta )=D^{}(\beta )()^{a^{}a}D(\beta )`$ (3.51)
it follows

$`O={\displaystyle _{}}{\displaystyle \frac{d^2\beta }{\pi }}\text{Tr}\left[O4D^{}(\beta )()^{a^{}a}D(\beta )\right]\text{Tr}\left[D(\beta )\rho D^{}(\beta )()^{a^{}a}\right].`$ (3.52)
Hence, it is possible to estimate $`O`$ by repeated measurement of the parity operator on displaced versions of the state under investigation. An approximated implementation of this technique for a single mode radiation field has been suggested in Refs.  through the measurement of the photon number probability on states displaced by a beam splitter. A similar scheme has been used for the experimental determination of the motional quantum state of a trapped atom . In comparison with the approximated methods, Eq. (3.52) allows to directly obtain the estimator $`[O](\alpha )`$ for any operator $`O`$ for which the trace exists. For instance, the reconstruction of the density matrix in the Fock representation is obtained by averaging the estimator

$`[|nn+d||](\alpha )`$ $`=`$ $`4n+d|D^{}(\alpha )()^{a^{}a}D(\alpha )|n`$
$`=`$ $`4()^{n+d}\sqrt{{\displaystyle \frac{n!}{(n+d)!}}}(2\alpha )^de^{2|\alpha |^2}L_n^d(4|\alpha |^2),`$
without the need of artificial cut-off in the Fock space .

#### 3.3.5 Quantum estimation for spin systems

The spin tomographic methods of Refs.  allow the reconstruction of the quantum state of a spin system. These methods utilize measurements of the spin in different directions, i.e. the quorum is the set of operators of the form $`\stackrel{}{S}\stackrel{}{n}`$, where $`\stackrel{}{S}`$ is the spin operator and $`\stackrel{}{n}(\mathrm{cos}\phi \mathrm{sin}\vartheta ,\mathrm{sin}\phi \mathrm{sin}\vartheta ,\mathrm{cos}\vartheta )`$ is a varying unit vector. Different quorums can be used, that exploit different sets of directions.

The easiest choice for the set of directions $`\stackrel{}{n}`$ is to consider all possible directions. The procedure to derive the tomographic formulas for this quorum is analogous to the one employed in Sec. 3.3.3 for homodyne tomography. The reconstruction formula for spin tomography for the estimation of an arbitrary operator $`O`$ writes

$`O={\displaystyle \underset{m=s}{\overset{s}{}}}{\displaystyle _\mathrm{\Omega }}{\displaystyle \frac{d\stackrel{}{n}}{4\pi }}p(m,\stackrel{}{n})[O](m,\stackrel{}{n}),`$ (3.54)
where $`p(m,\stackrel{}{n})`$ is the probability of obtaining the eigenvalue $`m`$ when measuring the spin along direction $`\stackrel{}{n}`$, $`[O](m,\stackrel{}{n})`$ is the tomographic estimator for the operator $`O`$, and $`\mathrm{\Omega }`$ is the unit sphere. In this case the operators $`C(\lambda )`$ of Eq. (3.18) are given by the set of projectors over the eigenstates $`|m,\stackrel{}{n}`$ of the operators $`\stackrel{}{S}\stackrel{}{n}`$. Notice that this is a complete set of operators on the system Hilbert space $``$. In order to find the dual basis $`B`$, one must consider the unitary operators obtained by exponentiating the quorum, i.e. $`D(\psi ,\stackrel{}{n})=\mathrm{exp}(i\psi \stackrel{}{S}\stackrel{}{n})`$, which satisfy the bi-orthogonality condition (3.21). In fact, $`D(\psi ,\stackrel{}{n})`$ constitutes a unitary irreducible representation of the group $`𝔾=SU(2)`$, and the bi-orthogonality condition is just the orthogonality relations between the matrix elements of the group representation , i.e.

$`{\displaystyle _𝔾}𝑑gD_{jr}(g)D_{tk}^{}(g)={\displaystyle \frac{V}{d}}\delta _{jk}\delta _{tr},`$ (3.55)
where $`D`$ is a unitary irreducible representation of dimension $`d`$, $`dg`$ is the group Haar invariant measure, and $`V=_𝔾𝑑g`$. For $`𝔾=SU(2)`$, with the $`2s+1`$-dimensional unitary irreducible representation $`D(\psi ,\stackrel{}{n})`$ ($`\stackrel{}{n}S^2`$ unit vector on the sphere, and $`\psi [0,4\pi ]`$ the rotation angle around $`\stackrel{}{n}`$) the Haar’s invariant measure is $`\mathrm{sin}^2\frac{\psi }{2}\mathrm{sin}\vartheta d\vartheta d\phi d\psi `$, and $`\frac{V}{d}=\frac{8\pi ^2}{2s+1}`$. We need however to integrate only for $`\psi [0,2\pi ]`$ (the change of sign for $`2\pi `$ rotation is irrelevant), whence the bi-orthogonality condition writes

$`{\displaystyle \frac{2s+1}{4\pi ^2}}{\displaystyle _\mathrm{\Omega }}𝑑\stackrel{}{n}{\displaystyle _0^{2\pi }}𝑑\psi \mathrm{sin}^2{\displaystyle \frac{\psi }{2}}j|e^{i\psi \stackrel{}{n}\stackrel{}{S}}|rt|e^{i\psi \stackrel{}{n}\stackrel{}{S}}|k=\delta _{jk}\delta _{tr},`$ (3.56)
and hence the spin tomography identity is given by

$`O={\displaystyle \frac{2s+1}{4\pi ^2}}{\displaystyle _\mathrm{\Omega }}𝑑\stackrel{}{n}{\displaystyle _0^{2\pi }}𝑑\psi \mathrm{sin}^2{\displaystyle \frac{\psi }{2}}\text{Tr}\left[OD^{}(\psi ,\stackrel{}{n})\right]D(\psi ,\stackrel{}{n}).`$ (3.57)
Notice the analogy between Eq. (3.57) and Glauber’s formula (3.34). In fact, both homodyne and spin tomography can be derived using the method of Group Tomography , and the underlying groups are the Weyl-Heisenberg group and the $`SU(2)`$ group, respectively. Formula (3.54) is obtained from Eq. (3.57) through the expectation value calculated on the eigenstates of $`\stackrel{}{S}\stackrel{}{n}`$. Thus, the explicit form of the tomographic estimator is obtained as

$`[O](m,\stackrel{}{n})={\displaystyle \frac{2s+1}{\pi }}{\displaystyle _0^{2\pi }}𝑑\psi \mathrm{sin}^2{\displaystyle \frac{\psi }{2}}\text{Tr}\left[Oe^{i\psi \stackrel{}{S}\stackrel{}{n}}\right]e^{i\psi m}.`$ (3.58)
As already noticed, there are other possible quorums for spin tomography. For example, for spin $`s=\frac{1}{2}`$ systems, a self-dual basis for the operator space is given by the identity and the Pauli matrices. In fact, from the properties $`\text{Tr}[\sigma _\alpha ]=0`$ and $`\sigma _\alpha \sigma _\beta =i_\gamma ϵ_{\alpha \beta \gamma }\sigma _\alpha \sigma _\beta `$ ($`\alpha ,\beta ,\gamma =x,y,z`$), both the bi-orthogonality relation (3.21) and the trace condition (3.23) follow. In this case the reconstruction formula writes

$`O={\displaystyle \frac{1}{2}}\text{Tr}\left[O\right]+{\displaystyle \frac{1}{2}}{\displaystyle \underset{\alpha =x,y,z}{}}{\displaystyle \underset{m=\pm \frac{1}{2}}{}}mp(m,\stackrel{}{n}_\alpha )\text{Tr}\left[O\sigma _\alpha \right].`$ (3.59)
In the case of generic $`s`$ spin system, Weigert has also shown  that by choosing $`(2s+1)^2`$ arbitrary directions for $`\stackrel{}{n}`$, it is possible to obtain (in almost all cases) a quorum of projectors $`|s,\stackrel{}{n}_js,\stackrel{}{n}_j|`$ ($`j=1,\mathrm{},(2s+1)^2`$), where $`|s,\stackrel{}{n}_j`$ is the eigenstate pertaining to the maximum eigenvalue $`s`$ of $`\stackrel{}{S}\stackrel{}{n}_j`$.

#### 3.3.6 Quantum estimation for a free particle

The state of a moving packet can be inferred from position measurement at different times . Assuming a particle with unit mass and using normalized unit $`\mathrm{}/2=1`$, the free Hamiltonian is given by the square of momentum operator $`H_F=p^2`$. In terms of the eigenvectors $`|x`$ of the position operator and of the selfadjoint operator

$`R(x,\tau )=e^{ip^2\tau }|xx|e^{ip^2\tau },`$ (3.60)
the probability density of the position of the free particle at time $`\tau `$ is obtained as $`p(x,\tau )=\text{Tr}[\rho R(x,\tau )]`$. The operators $`R(x,\tau )`$ provide a self-dual basis, and an arbitrary particle state can be written as

$`\rho ={\displaystyle _{}}{\displaystyle _{}}𝑑x𝑑\tau p(x,\tau )R(x,\tau ).`$ (3.61)
### 3.4 Noise deconvolution and adaptive tomography

In this section we will analyze: 1) the noise deconvolution scheme of Refs. , that allows to eliminate the experimental noise that arises from imperfect detection and lossy devices; 2) the adaptive tomography technique of Ref.  that allows to tune the unbiased tomographic estimators to a specific sample of experimental data, in order to reduce the statistical noise.

#### 3.4.1 Noise deconvolution

In short, it is possible to eliminate detection noise when it is possible to invert the noise map. A noise process is described by a trace preserving completely positive map $`\mathrm{\Gamma }`$. The noise can be deconvolved at the data analysis if

* the inverse of $`\mathrm{\Gamma }`$ exists, namely $`\mathrm{\Gamma }^1:()()`$, with $`\mathrm{\Gamma }^1\left[\mathrm{\Gamma }\left[O\right]\right]=O`$, for $`O()`$.
* the estimator $`_\lambda [O](Q_\lambda )`$ is in the domain of $`\mathrm{\Gamma }^1`$
* the map $`\mathrm{\Gamma }^1\left[_\lambda [O](Q_\lambda )\right]`$ is a function of $`Q_\lambda `$.

If the above conditions are met, we can recover the “ideal” expectation value $`O`$ that we would get without noise. This is achieved by replacing $`_\lambda [O](Q_\lambda )`$ with $`\mathrm{\Gamma }^1[_\lambda [O](Q_\lambda )]`$, and evaluating the ensemble average with the state $`\mathrm{\Gamma }^\tau (\rho )`$, namely the state affected by the noise ($`\mathrm{\Gamma }^\tau `$ represents the dual map, that provides the evolution in the Schroedinger picture). Hence, one has

$`O={\displaystyle _\mathrm{\Lambda }}𝑑\lambda \text{Tr}[\mathrm{\Gamma }^1[_\lambda [O](Q_\lambda )]\mathrm{\Gamma }^\tau (\rho )]{\displaystyle _\mathrm{\Lambda }}𝑑\lambda \mathrm{\Gamma }^1[_\lambda [O](Q_\lambda )]_\mathrm{\Gamma }.`$ (3.62)
Consider for example the noise arising from non-unity quantum efficiency $`\eta `$ of homodyne detectors. Recall that the ideal probability density is replaced by a Gaussian convolution with rms $`\mathrm{\Delta }_\eta ^2=(1\eta )/(4\eta )`$. Then, the map $`\mathrm{\Gamma }_\eta `$ acts on the quorum as follows

$`\mathrm{\Gamma }_\eta [e^{ikX_\phi }]`$ $`=`$ $`{\displaystyle _{\mathrm{}}^+\mathrm{}}𝑑xe^{ikx}\mathrm{\Gamma }_\eta [|xx|]`$
$`=`$ $`{\displaystyle _{\mathrm{}}^+\mathrm{}}𝑑x{\displaystyle _{\mathrm{}}^+\mathrm{}}𝑑x^{}e^{ikx}e^{\frac{(xx^{})^2}{2\mathrm{\Delta }^2}}[|x^{}x^{}|]=e^{\frac{1}{2}\mathrm{\Delta }^2k^2}e^{ikX_\phi }.`$
Of course one has

$`\mathrm{\Gamma }_\eta ^1[e^{ikX_\phi }]=e^{\frac{1}{2}\mathrm{\Delta }^2k^2}e^{ikX_\phi }.`$ (3.64)
In terms of the Fourier transform of the estimator

$`\stackrel{~}{}[O](y,\phi )={\displaystyle _{\mathrm{}}^+\mathrm{}}{\displaystyle \frac{dx}{2\pi }}e^{ixy}[O](x,\phi ),`$ (3.65)
one has

$`\stackrel{~}{}_\eta [O](y,\phi )=e^{\frac{1}{2}\mathrm{\Delta }^2y^2}\stackrel{~}{}[O](y,\phi ).`$ (3.66)
We applied the above result in Sec. 3.3.3, where the effect of non-unity quantum efficiency for reconstructing the density matrix elements was discussed. The use of the estimator in Eq. (3.42) and the origin of the bound $`\eta >1/2`$ is now more clear.

Another simple example of noise deconvolution is given here for a spin $`1/2`$ system. Consider the map that describes the “depolarizing channel”

$`\mathrm{\Gamma }_p[O]=(1p)O+{\displaystyle \frac{p}{2}}\text{Tr}[O]I,0p1.`$ (3.67)
This map can be inverted for $`p1`$ as follows

$`\mathrm{\Gamma }_p^1[O]={\displaystyle \frac{1}{1p}}\left(O{\displaystyle \frac{p}{2}}\text{Tr}[O]I\right).`$ (3.68)
Then Eq. (3.59) can be replaced with

$`O={\displaystyle \frac{1}{2}}\text{Tr}\left[O\right]+{\displaystyle \frac{1}{2(1p)}}{\displaystyle \underset{m=\pm \frac{1}{2}}{}}{\displaystyle \underset{\alpha =x,y,z}{}}mp_p(m,\stackrel{}{n}_\alpha )\text{Tr}\left[O\sigma _\alpha \right],`$ (3.69)
where now $`p_p(m,\stackrel{}{n}_\alpha )`$ represents the probability of outcome $`m`$ when measuring $`\sigma _\alpha `$ on the noisy state $`\mathrm{\Gamma }_p[\rho ]`$.

#### 3.4.2 Adaptive tomography

The idea of adaptive tomography is that the tomographic null estimators of Eq. (3.30) can be used to reduce statistical errors. In fact, the addition of a null estimator in the ideal case of infinite statistics does not change the average since its mean value is zero, but can change the variance. Thus, one can look for a procedure to reduce the variance by adding suitable null functions. Consider the class of equivalent estimators for $`O`$

$`_\lambda ^{}[O](Q_\lambda )=_\lambda [O](Q_\lambda )+{\displaystyle \underset{i=1}{\overset{M}{}}}\nu _i𝒩_i(Q_\lambda ).`$ (3.70)
Each estimator in the class $`^{}`$ is identified by the coefficient vector $`\stackrel{}{\nu }`$. The variance of the tomographic averages can be evaluated as

$`\overline{\mathrm{\Delta }^2^{}[O]}=\overline{\mathrm{\Delta }^2[O]}+2{\displaystyle \underset{i=1}{\overset{M}{}}}\nu _i\overline{𝒩_i[O]}+{\displaystyle \underset{i,j=1}{\overset{M}{}}}\nu _i\nu _j\overline{𝒩_i𝒩_j},`$ (3.71)
where $`\overline{F}_\mathrm{\Lambda }𝑑\lambda F(Q_\lambda )`$, and

$`\overline{\mathrm{\Delta }^2[O]}=\overline{^2[O]}\overline{[O]}^2.`$ (3.72)
Minimizing $`\overline{\mathrm{\Delta }^2^{}[O]}`$ with respect to the coefficients $`\nu _i`$, one obtains the equation

$`{\displaystyle \underset{j=1}{\overset{M}{}}}\nu _j\overline{𝒩_i𝒩_j}=\overline{[O]𝒩_i},`$ (3.73)
which can be solved starting from the estimated mean values, with the vector $`\stackrel{}{\nu }`$ as unknown. Notice that the obtained vector $`\stackrel{}{\nu }`$ will depend on the experimental data, and has to be calculated with the above procedure for any new set of data.

In this way we obtain an adaptive tomographic algorithm, which consists in the following steps:

* Find the null estimators $`𝒩_i(Q_\lambda )(i=1,\mathrm{},M)`$ for the quorum which is being used in the experiment.
* Execute the experiment and collect the input data.
* Calculate, using the obtained data, the mean values $`\overline{𝒩_i𝒩_j}`$ and $`\overline{[O]𝒩_i}`$, and solve the linear system (3.73), to obtain $`\stackrel{}{\nu }`$.
* Use the vector $`\stackrel{}{\nu }`$ obtained in the previous step to build the ‘optimized estimator’ $`^{}[O](Q_\lambda )=[O](Q_\lambda )+_i\nu _i𝒩_i(Q_\lambda )`$. Using the data collected in the first step, the mean value $`O`$ is now evaluated as

$`O={\displaystyle _\mathrm{\Lambda }}𝑑\lambda _{}^{}{}_{\lambda }{}^{}[O](Q_\lambda ),`$ (3.74)
where the optimized estimator has been used.
* For each new set of data the whole procedure must be repeated, as $`\stackrel{}{\nu }`$ is dependent on the data.

Notice that also the experimental mean values are slightly modified in the adaptive tomographic process, since null estimators do not change mean values only in the limiting case of infinite statistics. Examples of simulations of the adaptive technique that efficiently reduce statistical noise of homodyne tomographic reconstructions can be found in Ref. . In homodyne tomography null estimators are obtained as linear combinations of the following functions

$`𝒩_{k,n}(X_\phi )=X_\phi ^ke^{\pm i(k+2+2n)\phi },k,n0.`$ (3.75)
One can easily check that such functions have zero average over $`\phi `$, independently on $`\rho `$. Hence, for every operator $`O`$ one actually has an equivalence class of infinitely many unbiased estimators, which differ by a linear combination of functions $`𝒩_{k,n}(X_\phi )`$. It is then possible to minimize the rms error in the equivalence class by the least-squares method, obtaining in this way an optimal estimator that is adapted to the particular set of experimental data.

## Chapter 4 Universal homodyning

As shown in Ref. , homodyne tomography can be used as a kind of universal detector, for measuring generic field operators, however at expense of some additional noise. In this chapter, the general class of field operators that can be measured in this way is reviewed, which includes also operators that are inaccessible to heterodyne detection. In Ref.  the most relevant observables were analyzed—such as the intensity, the real, the complex field, and the phase—showing how their tomographic measurements are affected by noise that is always larger than the intrinsic noise of the direct detection of the considered observables. On the other hand, by comparing the noise from homodyne tomography with that from heterodyning (for those operators that can be measured in both ways), in Ref.  it was shown that for some operators homodyning is better than heterodyning when the mean photon number is sufficiently small, i.e. in the quantum regime, and in this chapter such comparison will be also reviewed.

### 4.1 Homodyning observables

Homodyne tomography provides the maximum achievable information on the quantum state of a single-mode radiation field through the use of the estimators in Sec. 3.3.3. In principle, the knowledge of the density matrix should allow one to calculate the expectation value also for unbounded operators. However, this is generally true only when one has an analytic knowledge of the density matrix, but it is not true when the matrix has been obtained experimentally. In fact, the Hilbert space is actually infinite dimensional, whereas experimentally one can achieve only a finite matrix, each element being affected by an experimental error. Notice that, even though the method allows one to extract any matrix element in the Hilbert space from the same bunch of experimental data, it is the way in which errors converge in the Hilbert space that determines the actual possibility of estimating the trace $`O=\text{Tr}[O\rho ]`$ for an arbitrary operator $`O`$. This issue has been debated in the set of papers of Ref. . Consider for example the number representation, and suppose that we want to estimate the average photon number $`a^{}a`$. In Ref.  it has been shown that for non-unit quantum efficiency the statistical error for the diagonal matrix element $`n|\rho |n`$ diverges faster than exponentially versus $`n`$, whereas for $`\eta =1`$ the error saturates for large $`n`$ to the universal value $`\epsilon _n=\sqrt{2/N}`$ that depends only on the number $`N`$ of experimental data, but is independent on both $`n`$ and on the quantum state. Even for the unrealistic case $`\eta =1`$, one can see immediately that the estimated expectation value $`a^{}a=_{n=0}^{H1}n\rho _{nn}`$ based on the measured matrix elements $`\rho _{nn}`$, will exhibit an unbounded error versus the truncated-space dimension $`H`$, because the non-vanishing error of $`\rho _{nn}`$ versus $`n`$ multiplies the increasing eigenvalue $`n`$.

Here, we report the estimators valid for any operator that admits a normal ordered expansion, giving the general class of operators that can be measured in this way, also as a function of the quantum efficiency $`\eta `$. Hence, from the same tomographic experiment, one can obtain not only the density matrix, but also the expectation value of various field operators, also unbounded, and including some operators that are inaccessible to heterodyne detection. However, the price to pay for such detection flexibility is that all measured quantities will be affected by noise. If one compares this noise with that from heterodyning (for those operators that can be measured in both ways), it turns out that for some operators homodyning is anyway less noisy than heterodyning, at least for small mean photon numbers. The procedure for estimating the expectation $`O`$ will be referred to as homodyning the observable $`O`$.

By homodyning the observable $`O`$ we mean averaging an appropriate estimator $`[O](x,\phi )`$, independent on the state $`\rho `$, over the experimental homodyne data, achieving in this way the expectation value $`O`$ for every state $`\rho `$, as in Eq. (3.36). For unbounded operators one can obtain the explicit form of the estimator $`[O](x,\phi )`$ in a different way. Starting from the identity involving trilinear products of Hermite polynomials

$`{\displaystyle _{\mathrm{}}^+\mathrm{}}𝑑xe^{x^2}H_k(x)H_m(x)H_n(x)={\displaystyle \frac{2^{\frac{m+n+k}{2}}\pi ^{\frac{1}{2}}k!m!n!}{(sk)!(sm)!(sn)!}},`$ (4.1)
for $`k+m+n=2s`$ even, Richter proved the following nontrivial formula for the expectation value of the normally ordered field operators

$`a^{}{}_{}{}^{n}a_{}^{m}={\displaystyle _0^\pi }{\displaystyle \frac{d\phi }{\pi }}{\displaystyle _{\mathrm{}}^+\mathrm{}}𝑑xp(x,\phi )e^{i(mn)\phi }{\displaystyle \frac{H_{n+m}(\sqrt{2}x)}{\sqrt{2^{n+m}}\left(\genfrac{}{}{0pt}{}{n+m}{n}\right)}},`$ (4.2)
which corresponds to the estimator

$`[a^{}{}_{}{}^{n}a_{}^{m}](x,\phi )=e^{i(mn)\phi }{\displaystyle \frac{H_{n+m}(\sqrt{2}x)}{\sqrt{2^{n+m}}\left(\genfrac{}{}{0pt}{}{n+m}{n}\right)}}.`$ (4.3)
This result can be easily extended to the case of non-unit quantum efficiency $`\eta <1`$. Using Eq. (3.66) one obtains

$`_\eta [a^{}{}_{}{}^{n}a_{}^{m}](x,\phi )=e^{i(mn)\phi }{\displaystyle \frac{H_{n+m}(\sqrt{2\eta }x)}{\sqrt{(2\eta )^{n+m}}\left(\genfrac{}{}{0pt}{}{n+m}{n}\right)}}.`$ (4.4)
From Eq. (4.4) by linearity one can obtain the estimator $`_\eta [f](x,\phi )`$ for any operator function $`f`$ that has normal ordered expansion

$`ff(a,a^{})={\displaystyle \underset{nm=0}{\overset{\mathrm{}}{}}}f_{nm}^{(N)}a^{}{}_{}{}^{n}a_{}^{m}.`$ (4.5)
From Eq. (4.4) one obtains

$`_\eta [f](x,\phi )`$ $`=`$ $`{\displaystyle \underset{s=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{H_s(\sqrt{2\eta }x)}{s!(2\eta )^{s/2}}}{\displaystyle \underset{nm=0}{\overset{\mathrm{}}{}}}f_{nm}^{(N)}e^{i(mn)\phi }n!m!\delta _{n+m,s}`$ (4.6)
$`=`$ $`{\displaystyle \underset{s=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{H_s(\sqrt{2\eta }x)i^s}{s!(2\eta )^{s/2}}}{\displaystyle \frac{d^s}{dv^s}}|_{v=0}[f](v,\phi ),`$
where

$`[f](v,\phi )={\displaystyle \underset{nm=0}{\overset{\mathrm{}}{}}}f_{nm}^{(N)}\left({\displaystyle \genfrac{}{}{0pt}{}{n+m}{m}}\right)^1(iv)^{n+m}e^{i(mn)\phi }.`$ (4.7)
Continuing from Eq. (4.6) one has

$`_\eta [f](x,\phi )=\mathrm{exp}\left({\displaystyle \frac{1}{2\eta }}{\displaystyle \frac{d^2}{dv^2}}+{\displaystyle \frac{2ix}{\sqrt{\eta }}}{\displaystyle \frac{d}{dv}}\right)|_{v=0}[f](v,\phi ),`$ (4.8)
and finally

$`_\eta [f](x,\phi )={\displaystyle _{\mathrm{}}^+\mathrm{}}{\displaystyle \frac{dw}{\sqrt{2\pi \eta ^1}}}e^{\frac{\eta }{2}w^2}[f](w+2ix/\sqrt{\eta },\phi ).`$ (4.9)
Hence one concludes that the operator $`f`$ can be measured by homodyne tomography if the function $`[f](v,\phi )`$ in Eq. (4.7) grows slower than $`\mathrm{exp}(\eta v^2/2)`$ for $`v\mathrm{}`$, and the integral in Eq. (4.9) grows at most exponentially for $`x\mathrm{}`$ (assuming $`p(x,\phi )`$ goes to zero faster than exponentially at $`x\mathrm{}`$).

The robustness to additive phase-insensitive noise of this method of homodyning observables has also been analyzed in Ref. , where it was shown that just half photon of thermal noise would spoil completely the measurement of the density matrix elements in the Fock representation.

In Table 4.1 we report the estimator $`_\eta [O](x,\phi )`$ for some operators $`O`$. The operator $`\widehat{W}_s`$ gives the generalized Wigner function $`W_s(\alpha ,\alpha ^{})`$ for ordering parameter $`s`$ through the relation in Eq. (2.14). From the expression of $`_\eta [\widehat{W}_s](x,\phi )`$ it follows that by homodyning with quantum efficiency $`\eta `$ one can measure the generalized Wigner function only for $`s<1\eta ^1`$: in particular the usual Wigner function for $`s=0`$ cannot be measured for any quantum efficiency.

### 4.2 Noise in tomographic measurements

In this section we will review the analysis of Ref. , where the tomographic measurement of following four relevant field quantities have been studied: the field intensity, the real field or quadrature, the complex field, and the phase. For all these quantities the conditions given after Eq. (4.9) are fulfilled.

The tomographic measurement of the observable $`O`$ is provided in terms of the average $`\overline{w_\eta }`$ of the estimator $`w_\eta _\eta [O](x,\phi )`$ over the homodyne data. The precision of the measurement is given by the confidence interval $`\sqrt{\overline{\mathrm{\Delta }w_\eta ^2}}`$. When $`w_\eta `$ is a real quantity, one has

$`\overline{\mathrm{\Delta }w_\eta ^2}=\overline{w_\eta ^2}\overline{w_\eta }^2,`$ (4.10)
where

$`\overline{w_\eta ^2}\overline{_\eta ^2[O](x,\phi )}={\displaystyle _0^\pi }{\displaystyle \frac{d\phi }{\pi }}{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑xp_\eta (x,\phi )_\eta ^2[O](x,\phi ).`$ (4.11)
When $`w_\eta `$ is complex, one has to consider the eigenvalues of the covariance matrix, namely

$`\overline{\mathrm{\Delta }w_\eta ^2}={\displaystyle \frac{1}{2}}\left[\overline{|w|_\eta ^2}|\overline{w_\eta }|^2\pm |\overline{w_\eta ^2}\overline{w_\eta }^2|\right].`$ (4.12)
When the observable $`O`$ can also be directly measured by a specific setup we can compare the tomographic precision $`\overline{\mathrm{\Delta }w^2}`$ with $`\mathrm{\Delta }O^2=O^2O^2`$. Notice that, when we deal with $`\eta <1`$ the noise $`\mathrm{\Delta }O_\eta ^2`$ is larger that the quantum fluctuations due to smearing effect of non-unit quantum efficiency. As we will see, the tomographic measurement is always more noisy than the corresponding direct measurement for any observable at any quantum efficiency $`\eta `$. This is not surprising, in view of the larger amount of information retrieved in the tomographic measurement as compared to the direct measurement of a single quantity.

According to Eq. (4.11), the evaluation of the added noise requires the average of the squared estimator. For the estimators in Eq. (4.4) it is very useful the following identity for the Hermite polynomials

$`H_n^2(x)=2^nn!^2{\displaystyle \underset{k=0}{\overset{n}{}}}{\displaystyle \frac{H_{2k}(x)}{k!^2\mathrm{\hspace{0.33em}2}^k(nk)!}},`$ (4.13)
that allows to write

$`_\eta ^2[a^na^m](x,\phi )=`$
$`e^{2i\phi (mn)}{\displaystyle \frac{n!^2m!^2}{\eta ^{m+n}}}{\displaystyle \underset{k=0}{\overset{m+n}{}}}{\displaystyle \frac{(2k)!\eta ^k}{k!^4(n+mk)!}}_\eta [a^ka^k](x,\phi ),`$ (4.14)
namely the squared estimator $`_\eta ^2[a^na^m](x,\phi )`$ can be written just in terms of ”diagonal” estimators $`_\eta [a^ka^k](x,\phi )`$.

#### 4.2.1 Field-Intensity

Photodetection is the direct measurement of the field-intensity. For non-unit quantum efficiency $`\eta `$, the probability of detecting $`m`$ photons is given by the Bernoulli convolution in Eq. (2.27). Let us consider the rescaled photocurrent

$`I_\eta ={\displaystyle \frac{1}{\eta }}a^{}a,`$ (4.15)
which traces the photon number, namely

$`I_\eta ={\displaystyle \frac{1}{\eta }}{\displaystyle \underset{m=0}{\overset{\mathrm{}}{}}}mp_\eta (m)=a^{}a\overline{n}.`$ (4.16)
The variance of $`I_\eta `$ is given by

$`\mathrm{\Delta }I_\eta ^2={\displaystyle \frac{1}{\eta ^2}}{\displaystyle \underset{m=0}{\overset{\mathrm{}}{}}}m^2p(m)\overline{n}^2=\mathrm{\Delta }n^2+\overline{n}\left({\displaystyle \frac{1}{\eta }}1\right),`$ (4.17)
where $`\mathrm{\Delta }n^2`$ denotes the intrinsic photon number variance, and $`\overline{n}(\eta ^11)`$ represents the noise introduced by inefficient detection. The tomographic estimator that traces the photon number is given by the phase-independent function $`w_\eta 2x^2(2\eta )^1`$. Using Eq. (4.14) we can evaluate its variance as follows

$`\overline{\mathrm{\Delta }w_\eta ^2}=\mathrm{\Delta }n^2+{\displaystyle \frac{1}{2}}n^2+\overline{n}\left({\displaystyle \frac{2}{\eta }}{\displaystyle \frac{3}{2}}\right)+{\displaystyle \frac{1}{2\eta ^2}}.`$ (4.18)
The noise $`N[n]`$ added by tomography in the measurement of the field intensity $`n`$ is then given by

$`N[n]=\overline{\mathrm{\Delta }w_\eta ^2}\mathrm{\Delta }I^2_\eta ={\displaystyle \frac{1}{2}}\left[n^2+\overline{n}\left({\displaystyle \frac{2}{\eta }}1\right)+{\displaystyle \frac{1}{\eta ^2}}\right].`$ (4.19)
Notice that $`N[n]`$ is always positive, and largely depends on state under examination. For coherent states we have the noise-ratio

$`\delta n_\eta =\sqrt{{\displaystyle \frac{\overline{\mathrm{\Delta }w_\eta ^2}}{\mathrm{\Delta }I^2_\eta }}}=\left[2+{\displaystyle \frac{1}{2}}\left(\eta \overline{n}+{\displaystyle \frac{1}{\eta \overline{n}}}\right)\right]^{1/2},`$ (4.20)
which is minimum for $`\overline{n}=\eta ^1`$.

#### 4.2.2 Real Field

For single mode radiation the electric field is proportional to a quadrature $`X=(a+a^{})/2`$, which is just traced by homodyne detection at fixed zero-phase with respect to the local oscillator. The tomographic estimator is given by $`w_\eta _\eta [X](x,\phi )=2x\mathrm{cos}\phi `$, independently on $`\eta `$, whereas the squared estimator $`_\eta ^2[X]`$ can be written as

$`w_\eta ^2={\displaystyle \frac{1}{4}}\left[_\eta [a^2](x,\phi )+_\eta [a^2](x,\phi )\right]+_\eta [a^{}a](x,\phi )+{\displaystyle \frac{1}{2\eta }}+{\displaystyle \frac{\mathrm{cos}(2\phi )}{2}}.`$ (4.21)
Then one has

$`\overline{\mathrm{\Delta }w_\eta ^2}`$ $`=`$ $`{\displaystyle \frac{1}{4}}\left[a^2+a^2\right]+\overline{n}+{\displaystyle \frac{1}{2\eta }}{\displaystyle \frac{1}{4}}a+a^{}^2`$ (4.22)
$`=`$ $`\mathrm{\Delta }X^2+{\displaystyle \frac{1}{2}}\overline{n}+{\displaystyle \frac{2\eta }{4\eta }},`$
where $`\mathrm{\Delta }X^2`$ represents the intrinsic quadrature fluctuations. The tomographic noise in Eq. (4.22) can be compared with the rms variance of direct homodyne detection (see Sec. 2.3)

$`\mathrm{\Delta }X^2_\eta =\mathrm{\Delta }X^2+{\displaystyle \frac{1\eta }{4\eta }}.`$ (4.23)
Then the added noise reads

$`N[X]={\displaystyle \frac{\overline{n}}{2}}+{\displaystyle \frac{1}{4\eta }}.`$ (4.24)
For coherent states $`\mathrm{\Delta }X^2=1/4`$, and one has the noise-ratio

$`\delta x_\eta =\sqrt{{\displaystyle \frac{\overline{\mathrm{\Delta }w_\eta ^2}}{\mathrm{\Delta }X^2_\eta }}}=\sqrt{2\eta \overline{n}+2}.`$ (4.25)
#### 4.2.3 Field amplitude

The detection of the complex field amplitude of a single-mode light-beam is represented by the generalized measurement of the annihilation operator $`a`$. The tomographic estimator for $`a`$ is given by the complex function $`w_\eta _\eta [a](x,\phi )=2x\mathrm{exp}(i\phi )`$, and the precision of the measurement is evaluated as in Eq. (4.12). From Eq. (4.14) one obtains

$`w_\eta ^2_\eta ^2[a](x,\phi )=e^{i2\phi }\left[{\displaystyle \frac{1}{\eta }}+2_\eta [a^{}a](x,\phi )\right]={\displaystyle \frac{e^{i2\phi }}{\eta }}+_\eta [a^2](x,\phi ),`$ (4.26)
and

$`\left|w_\eta \right|^2\left|_\eta [a](x,\phi )\right|^2={\displaystyle \frac{1}{\eta }}\left[1+2\eta _\eta [a^{}a](x,\phi )\right],`$ (4.27)
and hence

$`\overline{\mathrm{\Delta }w_\eta ^2}={\displaystyle \frac{1}{2}}\left[{\displaystyle \frac{1}{\eta }}+2\overline{n}\left|a\right|^2\pm \left|a^2a^2\right|\right].`$ (4.28)
The optimal measurement of the complex field $`a`$ is obtained through heterodyne detection. As notice in Sec. 2.4, the probability distribution is given by the generalized Wigner function $`W_s(\alpha ,\alpha ^{})`$, with $`s=1\frac{2}{\eta }`$. Using Eq. (2.71) the precision of the measurement is easily evaluated as follows

$`\mathrm{\Delta }a^2_\eta `$ $`=`$ $`{\displaystyle \frac{1}{2}}\left[\overline{|\alpha |^2}|\overline{\alpha }|^2\pm \left|\overline{\alpha ^2}\overline{\alpha }^2\right|\right]`$ (4.29)
$`=`$ $`{\displaystyle \frac{1}{2}}\left[\overline{n}+{\displaystyle \frac{1}{\eta }}|a|^2\pm \left|a^2a^2\right|\right].`$
The noise added by quantum tomography then reads

$`N[a]={\displaystyle \frac{1}{2}}\overline{n},`$ (4.30)
which is independent on quantum efficiency. For a coherent state we have

$`\overline{\mathrm{\Delta }w_\eta ^2}={\displaystyle \frac{1}{2}}\left[\overline{n}+{\displaystyle \frac{1}{\eta }}\right],\mathrm{\Delta }a^2_\eta ={\displaystyle \frac{1}{2\eta }},`$ (4.31)
and the noise ratio then writes

$`\delta a_\eta =\sqrt{{\displaystyle \frac{\overline{\mathrm{\Delta }w_\eta ^2}}{\mathrm{\Delta }a^2_\eta }}}=\sqrt{1+\eta \overline{n}}.`$ (4.32)
#### 4.2.4 Phase

The canonical description of the quantum optical phase is given by the probability operator measure

$`d\mu (\phi )={\displaystyle \frac{d\phi }{2\pi }}{\displaystyle \underset{n,m=0}{\overset{\mathrm{}}{}}}\mathrm{exp}[i(mn)\phi ]|nm|.`$ (4.33)
However, no feasible setup is known that achieves the optimal measurement (4.33). For this reason, here we consider the heterodyne measurement of the phase, and compare it with the phase of the tomographic estimator for the corresponding field operator $`a`$, i.e. $`w_\eta =\mathrm{arg}(2xe^{i\phi })`$. Notice that the phase $`w_\eta `$ does not coincide with the local oscillator phase $`\phi `$, because $`x`$ has varying sign. The probability distribution of $`w_\eta `$ can be obtained by the following identity

$`{\displaystyle _0^\pi }{\displaystyle \frac{d\phi }{\pi }}{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑xp_\eta (x,\phi )=1={\displaystyle _\pi ^\pi }{\displaystyle \frac{dw_\eta }{\pi }}{\displaystyle _0^{\mathrm{}}}𝑑xp_\eta (x,w_\eta ),`$ (4.34)
which implies

$`p_\eta (w_\eta )={\displaystyle \frac{1}{\pi }}{\displaystyle _0^{\mathrm{}}}𝑑xp_\eta (x,w_\eta ).`$ (4.35)
The precision in the tomographic phase measurement is given by the rms variance $`\overline{\mathrm{\Delta }w_\eta ^2}`$ of the probability (4.35). In the case of a coherent state with positive amplitude $`|\beta ||\beta |`$, Eq. (4.35) gives

$`p_\eta (w_\eta )={\displaystyle \frac{1}{2\pi }}\left[1+\text{Erf}\left({\displaystyle \frac{\sqrt{2}|\beta |\mathrm{cos}w_\eta }{\sqrt{\eta }}}\right)\right],`$ (4.36)
which approaches a ”boxed” distribution in $`[\pi /2,\pi /2]`$ for large intensity $`|\beta |1`$. We compare the tomographic phase measurement with heterodyne detection, namely the phase of the direct-detected complex field $`a`$. The outcome probability distribution is the marginal distribution of the generalized Wigner function $`W_s(\alpha ,\alpha ^{})`$ ($`s=1\frac{2}{\eta }`$) integrated over the radius

$`p_\eta (\phi )={\displaystyle _0^{\mathrm{}}}\rho 𝑑\rho W_s(\rho e^{i\phi },\rho e^{i\phi }),`$ (4.37)
whereas the precision in the phase measurement is given by its rms variance $`\overline{\mathrm{\Delta }\phi _\eta ^2}`$. We are not able to give a closed formula for the added noise $`N[\phi ]=\overline{\mathrm{\Delta }w_\eta ^2}\overline{\mathrm{\Delta }\phi _\eta ^2}`$. However, for high excited coherent states $`|\beta ||\beta |`$ (zero mean phase) one has $`\overline{\mathrm{\Delta }w_\eta ^2}=\pi ^2/12`$ and $`\overline{\mathrm{\Delta }\phi _\eta ^2}=(2\eta \overline{n})^1`$. The asymptotic noise-ratio is thus given by

$`\delta \phi _\eta =\sqrt{{\displaystyle \frac{\overline{\mathrm{\Delta }y_\eta ^2}}{\overline{\mathrm{\Delta }\phi _\eta ^2}}}}=\pi \sqrt{{\displaystyle \frac{\eta \overline{n}}{6}}},\overline{n}1.`$ (4.38)
A comparison for low excited coherent states can be performed numerically. The noise ratio $`\delta \phi _\eta `$ (expressed in dB) is shown in Fig. 4.1 for some values of the quantum efficiency $`\eta `$.

It is apparent that the tomographic determination of the phase is more noisy than heterodyning also in this low-intensity regime.

In Table 4.2 a synthesis of the results of this section is reported. We have considered the ratio between the tomographic and the direct-measurement noise. This is an increasing function of the mean photon number $`\overline{n}`$, however scaled by the quantum efficiency $`\eta `$. Therefore homodyne tomography turns out to be a very robust detection scheme for low quantum efficiency.

In Fig. 4.2 the coherent-state noise ratio (in dB) for all the considered quantities are plotted for unit quantum efficiency versus $`\overline{n}`$.

In conclusion, homodyne tomography adds larger noise for highly excited states, however, it is not too noisy in the quantum regime of low $`\overline{n}`$. It is then very useful in this regime, where currently available photodetectors suffer most limitations. Indeed, it has been adopted in recent experiments of photodetection .

### 4.3 Comparison between homodyne tomography and heterodyning

We have seen that homodyne tomography allows to measure any field observable $`ff(a,a^{})`$ having normal ordered expansion $`ff^{(N)}(a,a^{})=_{nm=0}^{\mathrm{}}f_{nm}^{(N)}a^{}{}_{}{}^{n}a_{}^{m}`$ and bounded integral in Eq. (4.9). On the other hand, as shown in Sec. 2.4, heterodyne detection allows to measure field observables that admit anti-normal ordered expansion $`ff^{(A)}(a,a^{})=_{nm=0}^{\mathrm{}}f_{nm}^{(A)}a^ma^{}^n`$, in which case the expectation value is obtained through the heterodyne average

$`f={\displaystyle _{}}{\displaystyle \frac{d^2\alpha }{\pi }}f^{(A)}(\alpha ,\alpha ^{})\alpha |\rho |\alpha .`$ (4.39)
As shown in Sec. 2.4, for $`\eta =1`$ the heterodyne probability is just the $`Q`$-function $`Q(\alpha ,\alpha ^{})=\frac{1}{\pi }\alpha |\rho |\alpha `$, whereas for $`\eta <1`$ it is Gaussian convoluted with rms $`(1\eta )/\eta `$, thus giving the Wigner function $`W_s(\alpha ,\alpha ^{})`$, with $`s=1\frac{2}{\eta }`$.

Indeed, the problem of measurability of the observable $`f`$ through heterodyne detection is not trivial, since one needs the admissibility of anti-normal ordered expansion and the convergence of the integral in Eq. (4.39). We refer the reader to Refs.  for more details and to Refs.  for analysis of quantum state estimates based on heterodyne detection.

The additional noise in homodyning the complex field $`a`$ has been evaluated in Eq. (4.30), where we found that homodyning is always more noisy than heterodyning. On the other hand, for other field observables it may happen that homodyne tomography is less noisy than heterodyne detection. For example, the added noise in homodyning the intensity $`a^{}a`$ with respect to direct detection has been evaluated in Eq. (4.19). Analogously, one can easily evaluate the added noise $`N_{het}[n]`$ when heterodyning the photon number $`n=a^{}a`$. According to Eq. (2.71), the random variable corresponding to the photon number for heterodyne detection with quantum efficiency $`\eta `$ is $`\nu (\alpha )=|\alpha |^2\frac{1}{\eta }`$. From the relation

$`\overline{|\alpha |^4}=a^2a^2+4{\displaystyle \frac{1\eta }{\eta }}aa^{}+2\left({\displaystyle \frac{1\eta }{\eta }}\right)^2`$ (4.40)
one obtains

$`\overline{\mathrm{\Delta }\nu (\alpha )}^2=\mathrm{\Delta }n^2+\overline{n}\left({\displaystyle \frac{2}{\eta }}1\right)+{\displaystyle \frac{1}{\eta ^2}}.`$ (4.41)
Upon comparing with Eq. (4.17), one concludes that the added noise in heterodyning the photon number is given by

$`N_{het}[n]=\overline{\mathrm{\Delta }\nu ^2(z)}\mathrm{\Delta }I_\eta ^2={\displaystyle \frac{1}{\eta ^2}}(\eta \overline{n}1).`$ (4.42)
With respect to the added noise in homodyning of Eq. (4.19) one has

$`N_{het}[n]=N[n]{\displaystyle \frac{1}{2}}\left(n^2\overline{n}{\displaystyle \frac{1}{\eta ^2}}\right).`$ (4.43)
Since $`n^2\overline{n}^2`$, we can conclude that homodyning the photon number is less noisy than heterodyning it for sufficiently low mean photon number $`n<\frac{1}{2}\left(1+\sqrt{1+\frac{4}{\eta ^2}}\right)`$.

## Chapter 5 Multimode homodyne tomography

The generalization of homodyne tomography from a single-mode to a multimode field is quite obvious, the estimator of simple operator tensors $`O=O_1O_2\mathrm{}O_n`$ being just the product of the estimators of each single-mode operator $`O_1,O_1,\mathrm{},O_n`$. By linearity, one then obtains also the estimator for arbitrary multimode operators. Such a simple generalization, however, requires a separate homodyne detector for each mode, which is unfeasible when the modes of the field are not spatio-temporally separated. This is the case, for example of pulsed fields, for which a general multimode tomographic method is especially needed, also due to the problem of mode matching between the local oscillator and the detected fields (determined by their relative spatio-temporal overlap) , which produces a dramatic reduction of the overall quantum efficiency.

In this Chapter we review the general method of Ref.  for homodyning observables of a multimode electromagnetic field using a single local oscillator (LO), providing the rule to evaluate the estimator of an arbitrary multimode operator. The expectation value of the operator can then be obtained by averaging the estimator over the homodyne outcomes that are collected using a single LO whose mode randomly scans all possible linear combinations of incident modes. We will then specifically consider some observables for a two-mode field in a state corresponding to a twin-beam produced by parametric downconversion, and prove the reliability of the method on the basis of computer simulations.

Finally, we report some experimental results  obtained in the Prem Kumar’s lab at Northwestern University. Such experiment actually represents the first measurement of the joint photon-number probability distribution of the twin-beam state.

### 5.1 The general method

The Hilbert-Schmidt operator expansion in Eq. (3.34) can be generalized to any number of modes as follows

$`O`$ $`=`$ $`{\displaystyle _{}}{\displaystyle \frac{d^2z_0}{\pi }}{\displaystyle _{}}{\displaystyle \frac{d^2z_1}{\pi }}\mathrm{}{\displaystyle _{}}{\displaystyle \frac{d^2z_M}{\pi }}\text{Tr}\left\{O\mathrm{exp}\left[{\displaystyle \underset{l=0}{\overset{M}{}}}\left(z_la_l^{}+z_l^{}a_l\right)\right]\right\}`$ (5.1)
$`\times `$ $`\mathrm{exp}\left[{\displaystyle \underset{l=0}{\overset{M}{}}}\left(z_la_l^{}z_l^{}a_l\right)\right],`$
where $`a_l`$ and $`a_l^{}`$, with $`l=0,\mathrm{},M`$ and $`[a_l,a_l^{}^{}]=\delta _{ll^{}}`$, are the annihilation and creation operators of $`M+1`$ independent modes, and $`O`$ now denotes an operator over all modes. Using the following hyper-spherical parameterization for $`z_l`$

$`z_0={\displaystyle \frac{i}{2}}ku_0(\stackrel{}{\theta })e^{i\psi _0}`$ $``$ $`{\displaystyle \frac{i}{2}}ke^{i\psi _0}\mathrm{cos}\theta _1,`$ (5.2)
$`z_1={\displaystyle \frac{i}{2}}ku_1(\stackrel{}{\theta })e^{i\psi _1}`$ $``$ $`{\displaystyle \frac{i}{2}}ke^{i\psi _1}\mathrm{sin}\theta _1\mathrm{cos}\theta _2,`$
$`z_2={\displaystyle \frac{i}{2}}ku_2(\stackrel{}{\theta })e^{i\psi _2}`$ $``$ $`{\displaystyle \frac{i}{2}}ke^{i\psi _2}\mathrm{sin}\theta _1\mathrm{sin}\theta _2\mathrm{cos}\theta _3,`$
$`\mathrm{}`$
$`z_{M1}={\displaystyle \frac{i}{2}}ku_{M1}(\stackrel{}{\theta })e^{i\psi _{M1}}`$ $``$ $`{\displaystyle \frac{i}{2}}ke^{i\psi _{M1}}\mathrm{sin}\theta _1\mathrm{sin}\theta _2\mathrm{}\mathrm{sin}\theta _{M1}\mathrm{cos}\theta _M,`$
$`z_M={\displaystyle \frac{i}{2}}ku_M(\stackrel{}{\theta })e^{i\psi _M}`$ $``$ $`{\displaystyle \frac{i}{2}}ke^{i\psi _M}\mathrm{sin}\theta _1\mathrm{sin}\theta _2\mathrm{}\mathrm{sin}\theta _{M1}\mathrm{sin}\theta _M,`$
where $`k[0,\mathrm{})`$; $`\psi _l[0,2\pi ]`$ for $`l=0,1,\mathrm{},M`$; and $`\theta _l[0,\pi /2]`$ for $`l=1,2,\mathrm{},M`$, Eq. (5.1) can be rewritten as follows:

$`O={\displaystyle 𝑑\mu [\stackrel{}{\psi }]𝑑\mu [\stackrel{}{\theta }]_0^+\mathrm{}𝑑k\left(\frac{k}{2}\right)^{2M+1}\frac{1}{M!}\text{Tr}[Oe^{ikX(\stackrel{}{\theta },\stackrel{}{\psi })}]e^{ikX(\stackrel{}{\theta },\stackrel{}{\psi })}}.`$ (5.3)
Here we have used the notation

$`{\displaystyle 𝑑\mu [\stackrel{}{\psi }]}{\displaystyle \underset{l=0}{\overset{M}{}}}{\displaystyle _0^{2\pi }}{\displaystyle \frac{d\psi _l}{2\pi }},`$ (5.4)
$`{\displaystyle 𝑑\mu [\stackrel{}{\theta }]}2^MM!{\displaystyle \underset{l=1}{\overset{M}{}}}{\displaystyle _0^{\pi /2}}𝑑\theta _l\mathrm{sin}^{2(Ml)+1}\theta _l\mathrm{cos}\theta _l,`$ (5.5)
$`X(\stackrel{}{\theta },\stackrel{}{\psi })={\displaystyle \frac{1}{2}}\left[A^{}(\stackrel{}{\theta },\stackrel{}{\psi })+A(\stackrel{}{\theta },\stackrel{}{\psi })\right],`$ (5.6)
$`A(\stackrel{}{\theta },\stackrel{}{\psi })={\displaystyle \underset{l=0}{\overset{M}{}}}e^{i\psi _l}u_l(\stackrel{}{\theta })a_l.`$ (5.7)
From the parameterization in Eq. (5.1), one has $`_{l=0}^Mu_l^2(\stackrel{}{\theta })=1`$, and hence $`[A(\stackrel{}{\theta },\stackrel{}{\psi }),A^{}(\stackrel{}{\theta },\stackrel{}{\psi })]=1`$, namely $`A(\stackrel{}{\theta },\stackrel{}{\psi })`$ and $`A^{}(\stackrel{}{\theta },\stackrel{}{\psi })`$ themselves are annihilation and creation operators of a bosonic mode. By scanning all values of $`\theta _l[0,\pi /2]`$ and $`\psi _l[0,2\pi ]`$, all possible linear combinations of modes $`a_l`$ are obtained.

For the quadrature operator $`X(\stackrel{}{\theta },\stackrel{}{\psi })`$ in Eq. (5.6), one has the following identity for the moments generating function

$`e^{ikX(\stackrel{}{\theta },\stackrel{}{\psi })}=\mathrm{exp}\left({\displaystyle \frac{1\eta }{8\eta }}k^2\right){\displaystyle _{\mathrm{}}^+\mathrm{}}𝑑xe^{ikx}p_\eta (x;\stackrel{}{\theta },\stackrel{}{\psi }),`$ (5.8)
where $`p_\eta (x;\stackrel{}{\theta },\stackrel{}{\psi })`$ denotes the homodyne probability distribution of the quadrature $`X(\stackrel{}{\theta },\stackrel{}{\psi })`$ with quantum efficiency $`\eta `$. Generally, $`\eta `$ can depend on the mode itself, i.e., it is a function $`\eta =\eta (\stackrel{}{\theta },\stackrel{}{\psi })`$ of the selected mode. In the following, for simplicity, we assume $`\eta `$ to be mode independent, however. By taking the ensemble average on each side of Eq. (5.3) and using Eq. (5.8) one has

$`O={\displaystyle 𝑑\mu [\stackrel{}{\psi }]𝑑\mu [\stackrel{}{\theta }]_{\mathrm{}}^+\mathrm{}𝑑xp_\eta (x;\stackrel{}{\theta },\stackrel{}{\psi })_\eta [O](x;\stackrel{}{\theta },\stackrel{}{\psi })},`$ (5.9)
where the estimator $`_\eta [O](x;\stackrel{}{\theta },\stackrel{}{\psi })`$ has the following expression

$`_\eta [O](x;\stackrel{}{\theta },\stackrel{}{\psi })={\displaystyle \frac{\kappa ^{M+1}}{M!}}{\displaystyle _0^+\mathrm{}}𝑑te^{(1\frac{\kappa }{2})t+2i\sqrt{\kappa t}x}t^M\text{Tr}[Oe^{2i\sqrt{\kappa t}X(\stackrel{}{\theta },\stackrel{}{\psi })}],`$ (5.10)
with $`\kappa =2\eta /(2\eta 1)`$. Eqs. (5.9) and (5.10) allow to obtain the expectation value $`O`$ for any unknown state of the radiation field by averaging over the homodyne outcomes of the quadrature $`X(\stackrel{}{\theta },\stackrel{}{\psi })`$ for $`\stackrel{}{\theta }`$ and $`\stackrel{}{\psi }`$ randomly distributed according to $`d\mu [\stackrel{}{\psi }]`$ and $`d\mu [\stackrel{}{\theta }]`$. Such outcomes can be obtained by using a single LO that is prepared in the multimode coherent state $`_{l=0}^M|\gamma _l`$ with $`\gamma _l=e^{i\psi _l}u_l(\theta )K/2`$ and $`K1`$. In fact, in this case the rescaled zero-frequency photocurrent at the output of a balanced homodyne detector is given by

$`I={\displaystyle \frac{1}{K}}{\displaystyle \underset{l=0}{\overset{M}{}}}(\gamma _l^{}a_l+\gamma _la_l^{}),`$ (5.11)
which corresponds to the operator $`X(\stackrel{}{\theta },\stackrel{}{\psi })`$. In the limit of a strong LO ($`K\mathrm{}`$), all moments of the current $`I`$ correspond to the moments of $`X(\stackrel{}{\theta },\stackrel{}{\psi })`$, and the exact measurement of $`X(\stackrel{}{\theta },\stackrel{}{\psi })`$ is then realized. Notice that for modes $`a_l`$ with different frequencies, in the d.c. photocurrent in Eq. (5.11) each LO with amplitude $`\gamma _l`$ selects the mode $`a_l`$ at the same frequency (and polarization). For less-than-unity quantum efficiency, Eq. (5.8) holds.

Equation (5.10) can be applied to some observables of interest. In particular, one can estimate the matrix element $`\{n_l\}|R|\{m_l\}`$ of the multimode density operator $`R`$. This will be obtained by averaging the estimator

$`_\eta [|\{m_l\}\{n_l\}|](x;\stackrel{}{\theta },\stackrel{}{\psi })=e^{i_{l=0}^M(n_lm_l)\psi _l}{\displaystyle \frac{\kappa ^{M+1}}{M!}}`$
$`\times {\displaystyle \underset{l=0}{\overset{M}{}}}\left\{[i\sqrt{\kappa }u_l(\stackrel{}{\theta })]^{\mu _l\nu _l}\sqrt{{\displaystyle \frac{\nu _l!}{\mu _l!}}}\right\}`$
$`\times {\displaystyle _0^+\mathrm{}}dte^{t+2i\sqrt{\kappa t}x}t^{M+_{l=0}^M(\mu _l\nu _l)/2}{\displaystyle \underset{l=0}{\overset{M}{}}}L_{\nu _l}^{\mu _l\nu _l}[\kappa u_l^2(\stackrel{}{\theta })t],`$ (5.12)
where $`\mu _l=\text{max}(m_l,n_l)`$, $`\nu _l=\text{min}(m_l,n_l)`$, and $`L_n^\alpha (z)`$ denotes the generalized Laguerre polynomial. For diagonal matrix elements, Eq. (5.12) simplifies to

$`_\eta [|\{n_l\}\{n_l\}|](x;\stackrel{}{\theta },\stackrel{}{\psi })={\displaystyle \frac{\kappa ^{M+1}}{M!}}{\displaystyle _0^+\mathrm{}}𝑑te^{t+2i\sqrt{\kappa t}x}t^M{\displaystyle \underset{l=0}{\overset{M}{}}}L_{n_l}[\kappa u_l^2(\stackrel{}{\theta })t]`$ (5.13)
with $`L_n(z)`$ denoting the customary Laguerre polynomial in $`z`$. Using the following identity

$`L_n^{\alpha _0+\alpha _1+\mathrm{}+\alpha _M+M}(x_0+x_1+\mathrm{}+x_M)`$
$`={\displaystyle \underset{i_0+i_1+\mathrm{}+i_M=n}{}}L_{i_0}^{\alpha _0}(x_0)L_{i_1}^{\alpha _1}(x_1)\mathrm{}L_{i_M}^{\alpha _M}(x_M),`$ (5.14)
from Eq. (5.13) one can easily derive the estimator of the probability distribution of the total number of photons $`N=_{l=0}^Ma_l^{}a_l`$

$`_\eta [|nn|](x;\stackrel{}{\theta },\stackrel{}{\psi })`$ $`=`$ $`{\displaystyle \frac{\kappa ^{M+1}}{M!}}{\displaystyle _0^+\mathrm{}}𝑑te^{t+2i\sqrt{\kappa t}x}t^ML_n^M[\kappa t],`$ (5.15)
where $`|n`$ denotes the eigenvector of $`N`$ with eigenvalue $`n`$. Notice that the estimator in Eq. (5.13) does not depend on the phases $`\psi _l`$; only the knowledge of the angles $`\theta _l`$ is needed. For the estimator in Eq. (5.15), even the angles $`\theta _l`$ can be unknown.

Now we specialize to the case of only two modes $`a`$ and $`b`$ (i.e., M=1 and $`\stackrel{}{\theta }`$ is a scalar $`\theta `$). The joint photon-number probability distribution is obtained by averaging

$`_\eta [|n,mn,m|](x;\theta ,\psi _0,\psi _1)=`$
$`\kappa ^2{\displaystyle _0^+\mathrm{}}𝑑te^{t+2i\sqrt{\kappa t}x}tL_n(\kappa t\mathrm{cos}^2\theta )L_m(\kappa t\mathrm{sin}^2\theta ).`$ (5.16)
The estimator (5.15) of the probability distribution of the total number of photons can be written as

$`_\eta [|nn|](x;\theta ,\psi _0,\psi _1)=\kappa ^2{\displaystyle _0^+\mathrm{}}𝑑te^{t+2i\sqrt{\kappa t}x}tL_n^1[\kappa t].`$ (5.17)
For the total number of photons one can also derive the estimator of the moment generating function, using the generating function for the Laguerre polynomials . One obtains

$`_\eta [z^{a^{}a+b^{}b}](x;\theta ,\psi _0,\psi _1)={\displaystyle \frac{1}{(z+\frac{1z}{\kappa })^2}}\mathrm{\Phi }(2,{\displaystyle \frac{1}{2}};{\displaystyle \frac{1z}{z+\frac{1z}{\kappa }}}x^2).`$ (5.18)
For the first two moments one obtains the simple expressions

$`_\eta [a^{}a+b^{}b](x;\theta ,\psi _0,\psi _1)=4x^2+{\displaystyle \frac{2}{\kappa }}2,`$ (5.19)
$`_\eta [(a^{}a+b^{}b)^2](x;\theta ,\psi _0,\psi _1)=8x^4+\left({\displaystyle \frac{24}{\gamma }}20\right)x^2+{\displaystyle \frac{6}{\gamma ^2}}{\displaystyle \frac{10}{\gamma }}+4.`$
It is worth noting that analogous estimators of the photon-number difference between the two modes are singular and one needs a cutoff procedure, similar to the one used in Ref.  for recovering the correlation between the modes by means of the customary two-mode tomography. In fact, in order to extract information pertaining to a single mode only one needs a delta-function at $`\theta =0`$ for mode $`a`$, or $`\theta =\pi /2`$ for mode $`b`$, and, in this case, one could better use the standard one-mode tomography by setting the LO to the proper mode of interest.

Finally, we note that for two-mode tomography the estimators can be averaged by the integral

$`O`$ $`=`$ $`{\displaystyle _0^{2\pi }}{\displaystyle \frac{d\psi _0}{2\pi }}{\displaystyle _0^{2\pi }}{\displaystyle \frac{d\psi _1}{2\pi }}{\displaystyle _1^1}{\displaystyle \frac{d(\mathrm{cos}2\theta )}{2}}{\displaystyle _{\mathrm{}}^+\mathrm{}}𝑑xp_\eta (x;\theta ,\psi _0,\psi _1)`$ (5.20)
$`\times `$ $`_\eta [O](x;\theta ,\psi _0,\psi _1)`$
over the random parameters $`\mathrm{cos}(2\theta ),\psi _0`$, and $`\psi _1`$. For example, in the case of two radiation modes having the same frequency but orthogonal polarizations, $`\theta `$ represents a random rotation of the polarizations, whereas $`\psi _0`$ and $`\psi _1`$ denote the relative phases between the LO and the two modes, respectively.

#### 5.1.1 Numerical results for two-mode fields

In this section we report some Monte-Carlo simulations from Ref.  to judge the experimental working conditions for performing the single-LO tomography on two-mode fields. We focus our attention on the twin-beam state, usually generated by spontaneous parametric downconversion, namely

$`|\mathrm{\Psi }=S(\chi )|0_a|0_b=\sqrt{1|\xi |^2}{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}\xi ^n|n_a|n_b,`$ (5.21)
where $`S(\chi )=\mathrm{exp}(\chi a^{}b^{}\chi ^{}ab)`$ and $`\xi =e^{i\mathrm{arg}\chi }\text{tanh}|\chi |`$. The parameter $`\xi `$ is related to the average number of photons per beam $`\overline{n}=|\xi |^2/(1|\xi |^2)`$. For the simulations we need to derive the homodyne probability distribution $`p(x;\theta ,\psi _0,\psi _1)`$ which is given by

$`p(x;\theta ,\psi _0,\psi _1)`$ $`=`$ $`\text{Tr}[U^{}|x_a{}_{a}{}^{}x|1_bU|\mathrm{\Psi }\mathrm{\Psi }|]`$
$`=`$ $`{}_{a}{}^{}0|{}_{b}{}^{}0|S^{}(\chi )U^{}[|x{}_{a}{}^{}{}_{a}{}^{}x|1_b]US(\chi )|0_{a}^{}|0_{b}^{},`$
where $`|x_a`$ is the eigenvector of the quadrature $`x=\frac{1}{2}(a^{}+a)`$ with eigenvalue $`x`$ and $`U`$ is the unitary operator achieving the mode transformation

$`U^{}\left({\displaystyle \genfrac{}{}{0pt}{}{a}{b}}\right)U=\left(\begin{array}{cc}e^{i\psi _0}\mathrm{cos}\theta \hfill & e^{i\psi _1}\mathrm{sin}\theta \hfill \\ e^{i\psi _1}\mathrm{sin}\theta \hfill & e^{i\psi _0}\mathrm{cos}\theta \hfill \end{array}\right)\left({\displaystyle \genfrac{}{}{0pt}{}{a}{b}}\right).`$ (5.25)
In the case of two radiation modes having the same frequency but orthogonal polarizations—the case of Type II phase-matched parametric amplifier—Eq. (5.1.1) gives the theoretical probability of outcome $`x`$ for the homodyne measurement at a polarization angle $`\theta `$ with respect to the polarization of the $`a`$ mode, and with $`\psi _0`$ and $`\psi _1`$ denoting the relative phases between the LO and the two modes, respectively. By using the Dirac-$`\delta `$ representation of the $`X`$-quadrature projector

$`|xx|={\displaystyle _{\mathrm{}}^+\mathrm{}}{\displaystyle \frac{d\lambda }{2\pi }}\mathrm{exp}[i\lambda (Xx)],`$ (5.26)
Eq. (5.1.1) can be rewritten as follows

$`p(x;\theta ,\psi _0,\psi _1)={\displaystyle _{\mathrm{}}^+\mathrm{}}{\displaystyle \frac{d\lambda }{2\pi }}{}_{a}{}^{}0|{}_{b}{}^{}0|S^{}(\chi )U^{}e^{i\lambda (X_ax)}US(\chi )|0_{a}^{}|0_{b}^{}`$
$`={\displaystyle _{\mathrm{}}^+\mathrm{}}{\displaystyle \frac{d\lambda }{2\pi }}e^{i\lambda x}{}_{a}{}^{}0|{}_{b}{}^{}0|\mathrm{exp}\{i{\displaystyle \frac{\lambda }{2}}[(e^{i\psi _0}\mu \mathrm{cos}\theta +e^{i\psi _1}\nu ^{}\mathrm{sin}\theta )a`$
$`+(e^{i\psi _0}\nu ^{}\mathrm{cos}\theta +e^{i\psi _1}\mu \mathrm{sin}\theta )b+\text{H.c.}]{\displaystyle \frac{}{}}\}|0_a|0_b,`$ (5.27)
where we have used Eq. (5.25) and the transformation

$`S^{}(\chi )\left({\displaystyle \genfrac{}{}{0pt}{}{a}{b^{}}}\right)S(\chi )=\left(\begin{array}{cc}\mu \hfill & \nu \hfill \\ \nu ^{}\hfill & \mu \hfill \end{array}\right)\left({\displaystyle \genfrac{}{}{0pt}{}{a}{b^{}}}\right)`$ (5.30)
with $`\mu =\text{cosh}|\chi |`$ and $`\nu =e^{i\mathrm{arg}\chi }\text{sinh}|\chi |`$. Upon defining

$`KC=e^{i\psi _0}\mu \mathrm{cos}\theta +e^{i\psi _1}\nu ^{}\mathrm{sin}\theta ,`$
$`KD=e^{i\psi _0}\nu ^{}\mathrm{cos}\theta +e^{i\psi _1}\mu \mathrm{sin}\theta ,`$ (5.31)
where $`K`$ and $`C,D`$, with $`|C|^2+|D|^2=1`$ one has

$`K^2=\mu ^2+|\nu |^2+2\mu |\nu |\mathrm{sin}2\theta \mathrm{cos}(\psi _0+\psi _1\mathrm{arg}\nu ).`$ (5.32)
Now, since the unitary transformation

$`\left(\begin{array}{cc}C\hfill & D\hfill \\ D^{}\hfill & C^{}\hfill \end{array}\right)\left({\displaystyle \genfrac{}{}{0pt}{}{a}{b}}\right)\left({\displaystyle \genfrac{}{}{0pt}{}{a}{b}}\right)`$ (5.35)
has no effect on the vacuum state, Eq. (5.27) leads to the following Gaussian distribution

$`p(x;\theta ,\psi _0,\psi _1)`$
$`={\displaystyle _{\mathrm{}}^+\mathrm{}}{\displaystyle \frac{d\lambda }{2\pi }}e^{i\lambda x}{}_{a}{}^{}0|{}_{b}{}^{}0|\mathrm{exp}\left\{iK{\displaystyle \frac{\lambda }{2}}[(Ca+Db)+\text{H.c.}]\right\}|0_{a}^{}|0_{b}^{}`$
$`={\displaystyle _{\mathrm{}}^+\mathrm{}}{\displaystyle \frac{d\lambda }{2\pi }}e^{i\lambda x}{}_{a}{}^{}0|\mathrm{exp}\left[iK{\displaystyle \frac{\lambda }{2}}(a+a^{})\right]|0_{a}^{}={\displaystyle \frac{1}{K}}|{}_{a}{}^{}0|x/K_{a}^{}|^2`$
$`={\displaystyle \frac{1}{\sqrt{2\pi \mathrm{\Delta }^2(\theta ,\psi _0,\psi _1)}}}\mathrm{exp}\left({\displaystyle \frac{x^2}{2\mathrm{\Delta }^2(\theta ,\psi _0,\psi _1)}}\right),`$ (5.36)
where the variance $`\mathrm{\Delta }^2(\theta ,\psi _0,\psi _1)`$ is given by

$`\mathrm{\Delta }^2(\theta ,\psi _0,\psi _1)={\displaystyle \frac{K^2}{4}}={\displaystyle \frac{1+|\xi |^2+2|\xi |\mathrm{sin}2\theta \mathrm{cos}(\psi _0+\psi _1\mathrm{arg}\xi )}{4(1|\xi |^2)}}.`$ (5.37)
Taking into account the Gaussian convolution that results from less-than-unity quantum efficiency, the variance just increases as

$`\mathrm{\Delta }^2(\theta ,\psi _0,\psi _1)\mathrm{\Delta }_\eta ^2(\theta ,\psi _0,\psi _1)=\mathrm{\Delta }^2(\theta ,\psi _0,\psi _1)+{\displaystyle \frac{1\eta }{4\eta }}.`$ (5.38)
Notice that the probability distribution in Eq. (5.36) corresponds to a squeezed vacuum for $`\theta =\frac{\pi }{4}`$ and $`\psi _0+\psi _1\mathrm{arg}\xi =0`$ or $`\pi `$.

We study the tomographic measurement of the joint photon-number probability distribution and the probability distribution for the total number of photons with use of the estimators in Eqs. (5.16) and (5.17), respectively. Moreover, using the estimator in Eq. (5.12) we reconstruct the matrix elements

$`C_{n,m}{}_{a}{}^{}m|{}_{b}{}^{}m|\mathrm{\Psi }\mathrm{\Psi }|n_a|n_{b}^{},`$ (5.39)
which reveal the coherence of the twin-beam state. Theoretically one should have

$`C_{n,m}=(1|\xi |^2)\xi ^m\xi ^n.`$ (5.40)
The estimators have been numerically evaluated by applying the Gauss method for calculating the integral in Eq. (5.12), which results in a fast and sufficiently precise algorithm with use of just 150 evaluation points.

In Fig. 5.1 a Monte-Carlo simulation of the joint photon-number probability distribution is reported. The simulated values compare very well with the theoretical ones. In Ref.  a careful analysis of the statistical errors has been done for various twin-beam states by constructing histograms of deviations of the results from different simulated experiments from the theoretical ones. In comparison to the customary two-LO tomography of Ref. , where for $`\eta =1`$ the statistical errors saturate for increasingly large $`n`$ and $`m`$, here we have statistical errors that are slowly increasing versus $`n`$ and $`m`$. This is due to the fact that the range of the estimators in Eq. (5.16) increases versus $`n`$ and $`m`$. Overall we find that for any given quantum efficiency the statistical errors are generally slightly larger than those obtained with the two-LO method. The convenience of using a single LO then comes with its own price tag.

By using the estimator in Eq. (5.17) the probability distribution for the total number of photons $`N`$ of the twin beams has been also constructed (Fig. 5.2). Notice the dramatic increase of error bars versus N and for smaller $`\eta `$.

Finally, in Fig. 5.3 we report the results of the tomographic measurement of $`C_{n,m}`$ defined in Eq. (5.39). Because the reconstructed $`C_{n,m}`$ is close to the theoretically expected value in Eq. (5.40), these reveal the purity of the twin beams, which cannot be inferred from the thermal diagonal distribution of Fig. 5.1.

The first experimental results of a measurement of the joint photon-number probability distribution for a two-mode quantum state created by a nondegenerate optical parametric amplifier has been presented in Ref. . In this experiment, however, the twin beams are detected separately by two balanced-homodyne detectors. A schematic of the experimental setup is reported in Fig. 5.4, and some experimental results are reported in Fig. 5.5. As expected for parametric fluorescence, the experiment has shown a measured joint photon-number probability distribution that exhibited up to 1.9 dB of quantum correlation between the two modes, with thermal marginal distributions.

## Chapter 6 Applications to quantum measurements

In this chapter we review a number of applications of quantum tomography related to some fundamental tests in quantum mechanics.

First, we report the proposal of Ref.  for testing the nonclassicality of quantum states by means of an operational criterion based on a set of quantities that can be measured experimentally with some given level of confidence, even in the presence of loss, noise, and less-than-unity quantum efficiency.

Second, we report the experiment proposed in Ref.  for testing quantum state reduction. The state-reduction rule is tested using optical homodyne tomography by directly measuring the fidelity between the theoretically-expected reduced state and the experimental state.

Finally, we review some experimental results obtained at the Quantum Optics Lab of the University of Naples  about the reconstruction of coherent signals, together with application to the estimation of the losses introduced by simple optical components.

### 6.1 Measuring the nonclassicality of a quantum state

The concept of nonclassical states of light has drawn much attention in quantum optics . The customary definition of nonclassicality is given in terms of the $`P`$-function presented in Sec. 2.1: a nonclassical state does not admit a regular positive $`P`$-function representation, namely, it cannot be written as a statistical mixture of coherent states. Such states produce effects that have no classical analogue. These kinds of states are of fundamental relevance not only for the demonstration of the inadequacy of classical description, but also for applications, e.g., in the realms of information transmission and interferometric measurements .

We are interested in testing the nonclassicality of a quantum state by means of a set of quantities that can be measured experimentally with some given level of confidence, even in the presence of loss, noise, and less-than-unity quantum efficiency. The positivity of the $`P`$-function itself cannot be adopted as a test, since there is no viable method to measure it. As proved in Sec. 4.1 only the generalized Wigner functions of order $`s<1\eta ^1`$ can be measured, $`\eta `$ being the quantum efficiency of homodyne detection. Hence, through this technique, all functions from $`s=1`$ to $`s=0`$ cannot be recovered, i.e., we cannot obtain the $`P`$-function and all its smoothed convolutions up to the customary Wigner function. For the same reason, the nonclassicality parameter proposed by Lee , namely, the maximum $`s`$-parameter that provides a positive distribution, cannot be experimentally measured.

Among the many manifestations of nonclassical effects, one finds squeezing, antibunching, even-odd oscillations in the photon-number probability, and negativity of the Wigner function . Any of these features alone, however, does not represent the univocal criterion we are looking for. Neither squeezing nor antibunching provides a necessary condition for nonclassicality . The negativity of the Wigner function, which is well exhibited by the Fock states and the Schrödinger-cat-like states, is absent for the squeezed states. As for the oscillations in the photon-number probability, some even-odd oscillations can be simply obtained by using a statistical mixture of coherent states.

Many authors  have adopted the non-positivity of the phase-averaged $`P`$-function $`F(I)=\frac{1}{2\pi }_0^{2\pi }𝑑\varphi P(I^{1/2}e^{i\varphi })`$ as the definition for a nonclassical state, since $`F(I)<0`$ invalidates Mandel’s semiclassical formula  of photon counting, i.e., it does not allow a classical description in terms of a stochastic intensity. Of course, some states can exhibit a “weak” nonclassicality , namely, a positive $`F(I)`$, but with a non-positive $`P`$-function (a relevant example being a coherent state undergoing Kerr-type self-phase modulation). However, from the point of view of the detection theory, such “weak” nonclassical states still admit a classical description in terms of positive intensity probability $`F(I)>0`$. For this reason, we adopt non-positivity of $`F(I)`$ as the definition of nonclassicality.

#### 6.1.1 Single-mode nonclassicality

The authors of Refs.  have pointed out some relations between $`F(I)`$ and generalized moments of the photon distribution, which, in turn, can be used to test the nonclassicality. The problem is reduced to an infinite set of inequalities that provide both necessary and sufficient conditions for nonclassicality . In terms of the photon-number probability $`p(n)=n|\rho |n`$ of the state with density matrix $`\rho `$, the simplest sufficient condition involves the following three-point relation

$`B(n)(n+2)p(n)p(n+2)(n+1)[p(n+1)]^2<0.`$ (6.1)
Higher-order sufficient conditions involve five-, seven-, …, $`(2k+1)`$-point relations, always for adjacent values of $`n`$. It is sufficient that just one of these inequalities is satisfied in order to assure the negativity of $`F(I)`$. Notice that for a coherent state $`B(n)=0`$ identically for all $`n`$.

In the following we show that quantum tomography can be used as a powerful tool for performing the nonclassicality test in Eq. (6.1). For less-than-unity quantum efficiency ($`\eta <1`$), we rely on the concept of a “noisy state” $`\rho _\eta `$, wherein the effect of quantum efficiency is ascribed to the quantum state itself rather than to the detector. In this model, the effect of quantum efficiency is treated in a Schrödinger-like picture, with the state evolving from $`\rho `$ to $`\rho _\eta `$, and with $`\eta `$ playing the role of a time parameter. Such lossy evolution is described by the master equation

$`_t\rho (t)={\displaystyle \frac{\mathrm{\Gamma }}{2}}\left[2a\rho (t)a^{}a^{}a\rho (t)\rho (t)a^{}a\right],`$ (6.2)
wherein $`\rho (t)\rho _\eta `$ with $`t=\mathrm{ln}\eta /\mathrm{\Gamma }`$.

For the nonclassicality test, reconstruction in terms of the noisy state has many advantages. In fact, for non-unit quantum efficiency $`\eta <1`$ the tomographic method introduces errors for $`p(n)`$ which are increasingly large versus $`n`$, with the additional limitation that quantum efficiency must be greater than the minimum value $`\eta =0.5`$. On the other hand, the reconstruction of the noisy-state probabilities $`p_\eta (n)=n|\rho _\eta |n`$ does not suffer such limitations, and even though all quantum features are certainly diminished in the noisy-state description, nevertheless the effect of non-unity quantum efficiency does not change the sign of the $`P`$-function, but only rescales it as follows:

$`P(z)P_\eta (z)={\displaystyle \frac{1}{\eta }}P(z/\eta ^{1/2}).`$ (6.3)
Hence, the inequality (6.1) still represents a sufficient condition for nonclassicality when the probabilities $`p(n)=n|\rho |n`$ are replaced with $`p_\eta (n)=n|\rho _\eta |n`$, the latter being given by a Bernoulli convolution, as shown in Eq. (2.27). When referred to the noisy-state probabilities $`p_\eta (n)`$, the inequality in Eq. (6.1) keeps its form and simply rewrites as follows

$`B_\eta (n)(n+2)p_\eta (n)p_\eta (n+2)(n+1)[p_\eta (n+1)]^2<0.`$ (6.4)
The quantities $`B(n)`$ and $`B_\eta (n)`$ are nonlinear in the density matrix. Then, they cannot be measured by averaging a suitable estimator over the homodyne data. Hence, in the evaluation of $`B(n)`$ one has to reconstruct the photon-number probabilities $`p(n)`$, using the estimator $`_\eta [|nn|](x,\phi )`$ in Eq. (3.44). The noisy-state probabilities $`p_\eta (n)`$ are obtained by using the same estimator for $`\eta =1`$, namely without recovering the convolution effect of non-unit quantum efficiency. Notice that the estimator does not depend on the phase of the quadrature. Hence, the knowledge of the phase of the local oscillator in the homodyne detector is not needed for the tomographic reconstruction, and it can be left fluctuating in a real experiment.

Regarding the estimation of statistical errors, they are generally obtained by dividing the set of homodyne data into blocks, as shown in Sec. 3.3.1. However, in the present case, the nonlinear dependence on the photon number probability introduces a systematic error that is vanishingly small for increasingly larger sets of data. Therefore, the estimated value of $`B(n)`$ is obtained from the full set of data, instead of averaging the mean value of the different statistical blocks.

In Figs. 6.16.2 some numerical results from Ref.  are reported, which are obtained by a Monte-Carlo simulation of a quantum tomography experiment. The nonclassicality criterion is tested either on a Schrödinger-cat state $`|\psi (\alpha )(|\alpha +|\alpha )`$ or on a squeezed state $`|\alpha ,rD(\alpha )S(r)|0`$, wherein $`|\alpha `$, $`D(\alpha )`$, and $`S(r)`$ denote a coherent state with amplitude $`\alpha `$, the displacement operator $`D(\alpha )=e^{\alpha a^{}\alpha ^{}a}`$, and the squeezing operator $`S(r)=e^{r(a^2a^2)/2}`$, respectively. Fig. 6.1 shows tomographically-obtained values of $`B(n)`$, with the respective error bars superimposed, along with the theoretical values for a Schrödinger-cat state and for a phase-squeezed state ($`r>0`$). For the same set of states the results for $`B_\eta (n)`$ obtained by tomographic reconstruction of the noisy state are reported in Fig. 6.2. Let us compare the statistical errors that affect the of $`B(n)`$ and $`B_\eta (n)`$ on the original and the noisy states, respectively. In the first case the error increases with $`n`$, whereas in the second it remains nearly constant, albeit with less marked oscillations in $`B_\eta (n)`$ than those in $`B(n)`$.

The nonclassicality of the states here analyzed is experimentally verifiable, as $`B_\eta (0)<0`$ by more than five standard deviations. In contrast, for coherent states one obtains small statistical fluctuations around zero for all $`n`$. Finally, we remark that the simpler test of checking for antibunching or oscillations in the photon-number probability in the case of the phase-squeezed state (left of Figs. 6.1 and 6.2) would not reveal the nonclassical features of such a state.

#### 6.1.2 Two-mode nonclassicality

In Ref.  it is also shown how quantum homodyne tomography can also be employed to test the nonclassicality of two-mode states. For a two-mode state nonclassicality is defined in terms of non-positivity of the following phase-averaged two-mode $`P`$-function :

$`F(I_1,I_2,\varphi )={\displaystyle \frac{1}{2\pi }}{\displaystyle _0^{2\pi }}𝑑\varphi _1P(I_1^{1/2}e^{i\varphi _1},I_2^{1/2}e^{i(\varphi _1+\varphi )}).`$ (6.5)
In Ref.  it is also proved that a sufficient condition for nonclassicality is

$`C=(n_1n_2)^2(n_1n_2)^2n_1+n_2<0,`$ (6.6)
where $`n_1`$ and $`n_2`$ are the photon-number operators of the two modes.

A tomographic test of the inequality in Eq. (6.6) can be performed by averaging the estimators for the involved operators using Table 4.1. Again, the value $`\eta =1`$ can be used to reconstruct the ensemble averages of the noisy state $`\rho _\eta `$. As an example, we consider the twin-beam state of Eq. (5.21). The theoretical value of $`C`$ is given by $`C=2|\xi |^2/(1|\xi |^2)<0`$. With regard to the effect of quantum efficiency $`\eta <1`$, the same argument still holds as for the single-mode case: one can evaluate $`C_\eta `$ for the twin beams degraded by the effect of loss, and use $`\eta =1`$ in the estimators. In this case, the theoretical value of $`C_\eta `$ is simply rescaled, namely

$`C_\eta =2\eta ^2|\xi |^2/(1|\xi |^2).`$ (6.7)
In Fig. 6.3 we report $`C_\eta `$ vs. $`1\eta `$, with $`\eta `$ ranging from 1 to 0.3 in steps of 0.05, for the twin beam in Eq. (5.21) with $`|\xi |^2=0.5`$, corresponding to a total average photon number $`n_1+n_2=2`$. The values of $`C_\eta `$ result from a Monte-Carlo simulation of a homodyne tomography experiment with a sample of $`4\times 10^5`$ data. The nonclassicality test in terms of the noisy state gives values of $`C_\eta `$ that are increasingly near the classically positive region for decreasing quantum efficiency $`\eta `$. However, the statistical error remains constant and is sufficiently small to allow recognition of the nonclassicality of the twin beams up to $`\eta =0.3`$.

We conclude that quantum homodyne tomography allows one to perform nonclassicality tests for single- and two-mode radiation states, even when the quantum efficiency of homodyne detection is rather low. The method involves reconstruction of the photon-number probability or of some suitable function of the number operators pertaining to the noisy state, namely, the state degraded by the less-than-unity quantum efficiency. The noisy-state reconstruction is affected by the statistical errors; however, they are sufficiently small that the nonclassicality of the state can be tested even for low values of $`\eta `$. For the cases considered here, we have shown that the nonclassicality of the states can be proved (deviation from classicality by many error bars) with $`10^5`$$`10^7`$ homodyne data. Moreover, since the knowledge of the phase of the local oscillator in the homodyne detector is not needed for the tomographic reconstruction, it can be left fluctuating in a real experiment.

### 6.2 Test of state reduction

In quantum mechanics the state reduction (SR) is still a very discussed rule. The so–called “projection postulate” was introduced by von Neumann  to explain the results from the Compton-Simons experiment, and it was generalized by Lüders  for measurements of observables with degenerate spectrum. The consistency of the derivation of the SR rule and its validity for generic measurements have been analyzed with some criticism . In a very general context, the SR rule was derived in a physically consistent way from the Schrödinger equation for the composite system of object and measuring apparatus . An experiment for testing quantum SR is therefore a very interesting matter. Such a test in general is not equivalent to a test of the repeatability hypothesis since the latter holds only for measurements of observables that are described by self-adjoint operators. For example, joint measurements like the Arthurs-Kelly  are not repeatable, as the reduced states are coherent states, which are not orthogonal.

Quantum optics offers a possibility of testing the SR, because several observables can be chosen to perform different measurements on a fixed system. For instance, one can decide to perform either homodyne or heterodyne, or photon-number detection. This is a unique opportunity; in contrast, in particle physics the measurements are mostly quasi-classical and restricted to only a few observables. In addition, optical homodyne tomography allows a precise determination of the quantum system after the SR.

A scheme for testing the SR could be based on tomographic measurements of the radiation density matrix after nondemolition measurements. However, such a scheme would reduce the number of observables that are available for the test. Instead, one can take advantage of the correlations between the twin beams of Eq. (5.21) produced by a non-degenerate optical parametric amplifier (NOPA), in which case one can test the SR even for demolitive-type measurements. Indeed, if a measurement is performed on one of the twin beams, the SR can be tested by homodyne tomography on the other beam. This is precisely the scheme for an experimental test of SR proposed in Ref. , which is reviewed in the following.

The scheme for the SR test is given in Fig. 6.4. Different kinds of measurements can be performed on beam 1, even though here the SR only for heterodyne detection and photon-number detection will be considered.

For a system described by a density operator $`\rho `$, the probability $`p(\lambda )d\lambda `$ that the outcome of a quantum measurement of an observable is in the interval $`[\lambda ,\lambda +d\lambda )`$ is given by Born’s rule $`p(\lambda )d\lambda =\text{Tr}[\rho \mathrm{\Pi }_\lambda d\lambda ]`$, where $`\mathrm{\Pi }_\lambda `$ is the POVM pertaining to the measurement that satisfies $`\mathrm{\Pi }_\lambda 0`$ and $`𝑑\lambda \mathrm{\Pi }_\lambda =I`$. For an exact measurement of an observable, which is described by a self-adjoint operator, $`\mathrm{\Pi }_\lambda `$ is just the projector over the eigenvector corresponding to the outcome $`\lambda `$. In the case of the photon number $`a^{}a`$ the spectrum is discrete and the POVM is $`\mathrm{\Pi }_m=|mm|`$ for integer eigenvalue $`m`$. For the Arthurs-Kelly joint measurement of the position and momentum (corresponding to a joint measurement of two conjugated quadratures of the field) we have the coherent-state POVM $`\mathrm{\Pi }_\alpha =\pi ^1|\alpha \alpha |`$.

When on beam 1 we perform a measurement described by $`\mathrm{\Pi }_\lambda `$, the reduced normalized state of beam 2 is

$$\rho (\lambda )=\frac{\text{Tr}_1[|\xi \xi |(\mathrm{\Pi }_\lambda 1)]}{\text{Tr}_{1,2}[|\xi \xi |(\mathrm{\Pi }_\lambda 1)]}=\frac{\mathrm{\Xi }\mathrm{\Pi }_\lambda ^\tau \mathrm{\Xi }^{}}{p(\lambda )},$$
(6.8)
where $`O^\tau `$ denotes the transposed operator (on a fixed basis), $`\mathrm{\Xi }=(1|\xi |^2)^{1/2}\xi ^{a^{}a}`$, and $`p(\lambda )=\text{Tr}_{1,2}[\mathrm{\Xi }\mathrm{\Pi }_\lambda ^\tau \mathrm{\Xi }^{}]`$ is the probability density of the measurement outcome $`\lambda `$. In the limit of infinite gain $`|\xi |1`$ one has $`\rho (\lambda )\mathrm{\Pi }_\lambda ^\tau `$. For example, for heterodyne detection with outcome $`\alpha `$, we have $`\rho (\alpha )=|\alpha ^{}\alpha ^{}|`$.

If the readout detector on beam 1 has quantum efficiency $`\eta _r`$, Eq. (6.8) is replaced with

$$\rho ^{\eta _r}(\lambda )=\frac{\mathrm{\Xi }(\mathrm{\Pi }_\lambda ^{\eta _r})^\tau \mathrm{\Xi }^{}}{p^{\eta _r}(\lambda )},$$
(6.9)
where $`p^{\eta _r}(\lambda )=\text{Tr}_{1,2}[\mathrm{\Xi }(\mathrm{\Pi }_\lambda ^{\eta _r})^\tau \mathrm{\Xi }^{}]`$, and $`\mathrm{\Pi }_\lambda ^{\eta _r}`$ is the POVM for measurement with quantum efficiency $`\eta _r`$. As shown in Sec. 2.4, for heterodyne detection one has the Gaussian convolution

$$\mathrm{\Pi }_\alpha ^{\eta _r}=\frac{1}{\pi }_{}\frac{d^2z}{\pi \mathrm{\Delta }_{\eta _r}^2}e^{\frac{|z\alpha |^2}{\mathrm{\Delta }_r^2}}|zz|,$$
(6.10)
with $`\mathrm{\Delta }_{\eta _r}^2=(1\eta _r)/\eta _r`$. For direct photodetection $`\mathrm{\Pi }_m=|mm|`$ is replaced with the Bernoulli convolution

$$\mathrm{\Pi }_m^{\eta _r}=\underset{j=m}{\overset{\mathrm{}}{}}\left(\begin{array}{c}j\\ m\end{array}\right)\eta _r^m(1\eta _r)^{jm}|jj|.$$
(6.11)
The experimental test proposed here consists of performing conditional homodyne tomography on beam 2, given the outcome $`\lambda `$ of the measurement on beam 1. We can directly measure the “fidelity of the test”

$$F(\lambda )=\text{Tr}[\rho ^{\eta _r}(\lambda )\rho _{\mathrm{meas}}(\lambda )],$$
(6.12)
where $`\rho ^{\eta _r}(\lambda )`$ is the theoretical state in Eq. (6.9), and $`\rho _{\mathrm{meas}}(\lambda )`$ is the experimentally measured state on beam 2. Notice that we use the term “fidelity” even if $`F(\lambda )`$ is a proper fidelity when at least one of the two states is pure, which occurs in the limit of unit quantum efficiency $`\eta _r`$. In the following we evaluate the theoretical value of $`F(\lambda )`$ and compare it with the tomographic measured value.

The fidelity (6.12) can be directly measured by homodyne tomography using the estimator for the operator $`\rho ^{\eta _r}(\lambda )`$, namely

$$F(\lambda )=_0^\pi \frac{\text{d}\phi }{\pi }_{\mathrm{}}^+\mathrm{}𝑑xp_{\eta _h}(x,\phi ;\lambda )_{\eta _h}[\rho ^{\eta _r}(\lambda )](x,\phi ),$$
(6.13)
where $`p_{\eta _h}(x,\phi ;\lambda )`$ is the conditional homodyne probability distribution for outcome $`\lambda `$ at the readout detector.

For heterodyne detection on beam 1 with outcome $`\alpha `$, the reduced state on beam 2 is given by the displaced thermal state

$$\rho ^{\eta _r}(\alpha )=\eta _\xi D(\gamma )(1\eta _\xi )^{a^{}a}D^{}(\gamma ),$$
(6.14)
where

$$\eta _\xi =1+(\eta _r1)|\xi |^2,\gamma =\frac{\xi \eta _r}{\eta _\xi }\alpha ^{}.$$
(6.15)
The estimator in Eq. (6.13) is given by

$`_{\eta _h}[\rho ^{\eta _r}(\alpha )](x,\phi )={\displaystyle \frac{2\eta _h\eta _\xi }{2\eta _h\eta _\xi }}\mathrm{\Phi }(1,{\displaystyle \frac{1}{2}};{\displaystyle \frac{2\eta _h\eta _\xi }{2\eta _h\eta _\xi }}(x\gamma _\phi )^2),`$ (6.16)
where $`\gamma _\phi =\text{Re}(\gamma e^{i\phi })`$, and $`\mathrm{\Phi }(a,b;z)`$ denotes the customary confluent hypergeometric function. The estimator in Eq. (6.16) is bounded for $`\eta _h>\frac{1}{2}\eta _\xi `$, then one needs to have

$$\eta _h>\frac{1}{2}\left[1|\xi |^2(1\eta _r)\right].$$
(6.17)
As one can see from Eq. (6.17), for $`\eta _h>0.5`$ the fidelity can be measured for any value of $`\eta _r`$ and any gain parameter $`\xi `$ of the NOPA. We recall that the condition $`\eta _h>0.5`$ is required for the measurement of the density matrix. However, in this direct measurement of the fidelity, the reconstruction of the density matrix is bypassed, and we see from Eq. (6.17) that the bound $`\eta _h=0.5`$ can be lowered.

The measured fidelity $`F(\alpha )`$ in Eq. (6.13) with $`\rho ^{\eta _r}(\alpha )`$ as given in Eq. (6.14) must be compared with the theoretical value

$$F_{\mathrm{th}}=\eta _\xi /(2\eta _\xi ),$$
(6.18)
that is independent of $`\alpha `$.

For direct photodetection on beam 1 with outcome $`n`$, the reduced state on beam 2 is given by

$`\rho ^{\eta _r}(n)=\eta _\xi \left({\displaystyle \frac{\eta _\xi }{1\eta _\xi }}\right)^n\left(\begin{array}{c}a^{}a\\ n\end{array}\right)(1\eta _\xi )^{a^{}a}.`$ (6.21)
The estimator for the fidelity measurement is

$`_{\eta _h}[\rho ^{\eta _r}(n)](x,\phi )={\displaystyle \frac{(\eta _\xi _z)^n}{n!}}|_{z=0}{\displaystyle \frac{2\eta _h\eta _\xi }{2\eta _h\eta _\xi +z}}\mathrm{\Phi }(1,{\displaystyle \frac{1}{2}};{\displaystyle \frac{2\eta _h(\eta _\xi z)}{2\eta _h\eta _\xi +z}}x^2).`$ (6.22)
We see that the same bound of Eq. (6.17) holds. In this case the measured fidelity $`F(n)`$ must be compared with the theoretical value

$$F_{\mathrm{th}}(n)=\eta _\xi ^{2+2n}F(2n+1,2n+1;1;(1\eta _\xi )^2),$$
(6.23)
where $`F(a,b;c;z)`$ denotes the customary hypergeometric function.

Several simulations have been reported in Ref.  for both heterodyne and photodetection on beam 1. In the former case the quadrature probability distribution has been simulated, pertaining to the reduced state (6.14) on beam 2, and averaged the estimators in Eq. (6.16). In the latter case the reduced state (6.21) and the estimators in Eq. (6.22) have been used.

Numerical results for the fidelity was thus obtained for different values of the quantum efficiencies $`\eta _r`$ and $`\eta _h`$, and of the NOPA gain parameter $`\xi `$. A decisive test can be performed with samples of just a few thousand measurements. The statistical error in the measurement was found rather insensitive to both quantum efficiencies and NOPA gain.

### 6.3 Tomography of coherent signals and applications

Quantum homodyne tomography has been proved useful in various experimental situations, such as for measuring the photon statistics of a semiconductor laser , for determining the density matrix of a squeezed vacuum  and the joint photon-number probability distribution of a twin beam created by a non-degenerate optical parametric amplifier , and, finally, for reconstructing the quantum states of spatial modes with an array detector . In this section we review some experimental results about homodyne tomography with coherent states, with application to the estimation of the loss introduced by simple optical components .

The experiment has been performed in the Quantum Optics Lab of the University of Naples, and its the schematic is presented in Fig. 6.5. The principal radiation source is provided by a monolithic Nd:YAG laser ($``$ 50 mW at 1064 nm; Lightwave, model 142). The laser has a linewidth of less than 10 kHz/ms with a frequency jitter of less than 300 kHz/s, while its intensity spectrum is shot–noise limited above 2.5 MHz.

The laser emits a linearly polarized beam in a TEM00 mode, which is split in two parts by a beam splitter. One part provides the strong local oscillator for the homodyne detector. The other part, typically less than 200 $`\mu `$W, is the homodyne signal. The optical paths traveled by the local oscillator and the signal beams are carefully adjusted to obtain a visibility typically above $`75\%`$ measured at one of the homodyne output port. The signal beam is modulated, by means of a phase electro–optic modulator (EOM, Linos Photonics PM0202), at 4MHz, and a halfwave plate (HWP<sub>2</sub>, HWP<sub>3</sub>) is mounted in each path to carefully match the polarization state at the homodyne input.

The detector is composed by a 50$`÷`$50 beam splitter (BS), two amplified photodiodes (PD1, PD2), and a power combiner. The difference photocurrent is demodulated at 4MHz by means of an electrical mixer. In this way the detection occurs outside any technical noise and, more important, in a spectral region where the laser does not carry excess noise.

The phase modulation added to the signal beam moves a certain number of photons, proportional to the square of the modulation depth, from the carrier optical frequency $`\omega `$ to the side bands at $`\omega \pm \mathrm{\Omega }`$ so generating two weak coherent states with engineered average photon number at frequencies $`\omega \pm \mathrm{\Omega }`$. The sum sideband modes is then detected as a controlled perturbation attached to the signal beam. The demodulated current is acquired by a digital oscilloscope (Tektronix TDS 520D) with 8 bit resolution and record length of 250000 points per run. The acquisition is triggered by a triangular shaped waveform applied to the PZT mounted on the local oscillator path. The piezo ramp is adjusted to obtain a $`2\pi `$ phase variation between the local oscillator and the signal beam in an acquisition window.

The homodyne data to be used for tomographic reconstruction of the state have been calibrated according to the noise of the vacuum state. This is obtained by acquiring a set of data leaving the signal beam undisturbed while scanning the local oscillator phase. It is important to note that in case of the vacuum state no role is played by the visibility at the homodyne beam–splitter.

The tomographic samples consist of $`N`$ homodyne data $`\{x_j,\phi _j\}_{j=1,\mathrm{},N}`$ with phases $`\phi _j`$ equally spaced with respect to the local oscillator. Since the piezo ramp is active during the whole acquisition time, we have a single value $`x_j`$ for any phase $`\phi _j`$. From calibrated data we first reconstruct the quantum state of the homodyne signal. According to the experimental setup, we expect a coherent signal with nominal amplitude that can be adjusted by varying the modulation depth of the optical mixer. However, since we do not compensate for the quantum efficiency of photodiodes in the homodyne detector ($`\eta 90\%`$) we expect to reveal coherent signals with reduced amplitude. In addition, the amplitude is further reduced by the non-maximum visibility (ranging from $`75\%`$ to $`85\%`$) at the homodyne beam–splitter.

In Fig. 6.6 we report a typical reconstruction, together with the reconstruction of the vacuum state used for calibration. For both states, we report the raw data, the photon number distribution $`\rho _{nn}`$, and a contour plot of the Wigner function. The matrix elements are obtained by sampling the corresponding estimators in Eq. (3.44), whereas the confidence intervals for diagonal elements are given by $`\delta \rho _{nn}=\mathrm{\Delta }\rho /\sqrt{N}`$, $`\mathrm{\Delta }\rho `$ being the rms deviation of the estimator over data. For off-diagonal elements the confidence intervals are evaluated for the real and imaginary part separately.

In order to see the quantum state as a whole, we also report the reconstruction of the Wigner function of the field, which can be expressed in terms of the matrix elements as the discrete Fourier transform

$`W(\alpha ,\alpha ^{})=\text{Re}{\displaystyle \underset{d=0}{\overset{\mathrm{}}{}}}e^{id\phi }{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}\mathrm{\Lambda }(n,d;|\alpha |)\rho _{n,n+d}`$ (6.24)
where $`\phi =\mathrm{arg}\alpha `$, and

$`\mathrm{\Lambda }(n,d;|\alpha |)=()^n2(2\delta _{d0})|2\alpha |^d\sqrt{{\displaystyle \frac{n!}{(n+d)!}}}e^{2|\alpha |^2}L_n^d(|2\alpha |^2),`$ (6.25)
$`L_n^d(x)`$ denoting the Laguerre polynomials. Of course, the series in Eq. (6.24) should be truncated at some point, and therefore the Wigner function can be reconstructed only at some finite resolution.

Once the coherence of the signal has been established we may use homodyne tomography to estimate the loss imposed by a passive optical component like an optical filter. The procedure may be outlined as follows. We first estimate the initial mean photon number $`\overline{n}_0=|\alpha _0|^2`$ of the signal beam, and then the same quantity inserting an optical neutral density filter in the signal path. If $`\mathrm{\Gamma }`$ is the loss parameter, then the coherent amplitude is reduced to $`\alpha _\mathrm{\Gamma }=\alpha _0e^\mathrm{\Gamma }`$, and the intensity to $`\overline{n}_\mathrm{\Gamma }=\overline{n}_0e^{2\mathrm{\Gamma }}`$.

The estimation of the mean photon number can be performed adaptively on data, using the general method presented in Sec. 3.4.2. One takes the average of the estimator

$`[a^{}a](x,\phi )=2x^2{\displaystyle \frac{1}{2}}+\mu e^{i2\phi }+\mu ^{}e^{i2\phi },`$ (6.26)
where $`\mu `$ is a parameter to be determined in order to minimize fluctuations. As proved in Ref.  one has $`\mu =1/2a^2`$, which itself can be obtained from homodyne data. In practice, one uses the data sample twice: first to evaluate $`\mu `$, then to obtain the estimate for the mean photon number.

In Fig. 6.7 the tomographic determinations of $`\overline{n}_\mathrm{\Gamma }`$ are compared with the expected values for three set of experiments, corresponding to three different initial amplitudes. The expected values are given by $`\overline{n}_\mathrm{\Gamma }=\overline{n}_0e^{2\mathrm{\Gamma }}𝒱`$, where $`\mathrm{\Gamma }`$ is the value obtained by comparing the signal dc currents $`I_0`$ and $`I_\mathrm{\Gamma }`$ at the homodyne photodiodes and $`𝒱=𝒱_\mathrm{\Gamma }/𝒱_0`$ is the relative visibility. The solid line in Fig. 6.7 denotes these values. The line is not continuous due to variations of visibility. It is apparent from the plot that the estimation is reliable in the whole range of values we could explore. It is worth noting that the estimation is absolute, i.e. it does not depend on the knowledge of the initial amplitude, and it is robust, since it can be performed independently on the quantum efficiency of the homodyne detector.

One may notice that the estimation of loss can be pursued also by measuring an appropriate observable, typically the intensity of the light beam with and without the filter. However, this is a concrete possibility only for high amplitude signals, whereas losses on weak coherent states cannot be properly characterized neither by direct photocounting using photodiodes (due to the low quantum efficiency and large fluctuations), nor by avalanche photodetectors (due to the impossibility of discriminating among the number of photons). On the contrary, homodyne tomography provides the mean intensity (actually the whole photon distribution) independently on the signal level, thus allowing a precise characterization also in the quantum regime. Indeed, in Ref.  adaptive tomographic determination of the mean photon number has been extensively applied to (numerically simulated) homodyne data for coherent states of various amplitudes. The analysis has shown that the determination is reliable also for small samples and that precision is not much affected by the intensity of the signal.

## Chapter 7 Tomography of a quantum device

If we want to determine experimentally the operation of a quantum device, we need, by definition, quantum tomography. In fact, the characterization of the device operation could be done by running a basis of possible known inputs, and determining the corresponding outputs by quantum tomography. In quantum mechanics the inputs are density operators, and the role of the transfer matrix is played by the so-called quantum operation of the device, here denoted by $``$. Thus the output state $`\rho _{out}`$ (a part from a possible normalization) is given by the quantum operation applied to the input state as follows

$$\rho _{out}=(\rho _{in}).$$
(7.1)
Since the set of states $`\rho `$ actually belongs to a space of operators, this means that if we want to characterize $``$ completely, we need to run a complete orthogonal basis of quantum states $`|n`$ ($`n=0,1,2,\mathrm{}`$), along with their linear combinations $`\frac{1}{\sqrt{2}}(|n^{}+i^k|n^{\prime \prime })`$, with $`k=0,1,2,3`$ and $`i`$ denoting the imaginary unit. However, the availability of such a set of states in the lab is, by itself, a very hard technological problem. For example, for an optical device, the states $`|n`$ are those with a precise number $`n`$ of photons, and, a part from very small $`n`$—say at most n=2—they have never been achieved in the lab, whereas preparing their superpositions remains a dream for experimentalists, especially if $`n1`$ (a kind of Schrődinger kitten states).

The idea of achieving the quantum operation of a device by scanning the inputs and making tomography of the corresponding output is the basis of the early methods proposed in Refs. . Due to the mentioned problems in the availability of input states, both methods have limited application. The method of Ref.  has been designed for NMR quantum processing, whereas the method of Ref.  was conceived for determining the Liouvillian of a phase-insensitive amplifier, namely for a case in which the quantum operation has no off-diagonal matrix elements, to evaluate which one needs the superpositions $`\frac{1}{\sqrt{2}}(|n^{}+i^k|n^{\prime \prime })`$ with $`k=0,1,2,3`$ mentioned above. The problem of availability of input spates and their superpositions was partially solved by the method of Ref. , where it was suggested to use randomly drawn coherent states to estimate the quantum operation of an optical device via a maximum likelihood approach. This method, however, cannot be used for quantum systems different from the em radiation—such as finite dimensional systems, i.e. qubits—due to the peculiarity of coherent states. The solution to the problem came with the recent method of Ref. , where the problem of the availability of input states was solved by using a single bipartite entangled input, which is equivalent to run all possible input states in a kind of “quantum parallel” fashion (bipartite entangled states are nowadays easily available in most quantum systems of interest). The method is also very simple and effective, and its experimental feasibility (for single-photon polarization-encoded qubits) has been already demonstrated in a recent experiment performed in the Francesco De Martini laboratory in Roma La Sapienza . In the next sections we will review the general method and report some computer simulated results from Ref. .

### 7.1 The method

As already mentioned, the most description of a general state-transformation in quantum mechanics is given in terms of the so-called quantum operation. The state transformation due to the quantum operation $``$ is given as follows

$$\rho \frac{(\rho )}{\text{Tr}\left((\rho )\right)}.$$
(7.2)
The transformation occurs with probability given by $`p=\text{Tr}[(\rho )]1`$. The quantum operation $``$ is a linear, trace-decreasing completely positive (CP) map. We remind that a map is completely positive if it preserves positivity generally when applied locally to an entangled state. In other words, upon denoting by $``$ the identical map on the Hilbert space $`𝒦`$ of a second quantum system, the extended map $``$ on $`𝒦`$ is positive for any extension $`𝒦`$. Typically, the CP map is written using a Kraus decomposition  as follows

$`(\rho )={\displaystyle \underset{n}{}}K_n\rho K_n^{},`$ (7.3)
where the operators $`K_n`$ satisfy

$`{\displaystyle \underset{n}{}}K_n^{}K_nI.`$ (7.4)
The transformation (7.3) occurs with generally non-unit probability $`\text{Tr}[(\rho )]1`$, and the probability is unit independently on $`\rho `$ when $``$ is trace-preserving, i.e. when we have the equal sign in Eq. (7.4). The particular case of unitary transformations corresponds to having just one term $`K_1=U`$ in the sum (7.3), with $`U`$ unitary. However, one can consider also non-unitary operations with one term only, namely

$`(\rho )=A\rho A^{},`$ (7.5)
where $`A`$ is a contraction, i. e. $`A1`$. Such operations leave pure states as pure, and describes, for example, the state reduction from a measurement apparatus for a particular fixed outcome that occurs with probability $`\text{Tr}[\rho A^{}A]1`$.

In the following we will use the notation for bipartite pure states introduced in Eq. (2.60), and we will denote by $`O^\tau `$ and $`O^{}`$ the transposed and the conjugate operator of $`O`$ with respect to some pre-chosen orthonormal basis.

The basic idea of the method in Ref.  is the following. An unknown quantum operation $``$ can be determined experimentally through quantum tomography, by exploiting the following one-to-one correspondence $`R_{}`$ between quantum operations $``$ and positive operators $`R_{}`$ on two copies of the Hilbert space $``$

$$R_{}=(|II|),(\rho )=\text{Tr}_2[I\rho ^\tau R_{}].$$
(7.6)
Notice that the vector $`|I`$ represents a (unnormalized) maximally entangled state. If we consider a bipartite input state $`|\psi `$ and operate with $``$ only on one Hilbert space as in Fig. 7.1, the output state is given by

$$R(\psi )(|\psi \psi |).$$
(7.7)
For invertible $`\psi `$ the two matrices $`R(I)R_{}`$ and $`R(\psi )`$ are related as follows

$$R(I)=(I\psi _{}^{1}{}_{}{}^{\tau }R(\psi )(I\psi ^1).$$
(7.8)
Hence, the (four-index) quantum operation matrix $`R_{}`$ can be obtained by estimating via quantum tomography the following ensemble averages

$$i,j|R(I)|l,k=\text{Tr}\left[R(\psi )\left(|li|\psi ^1|kj|\psi ^1\right)\right].$$
(7.9)
Then one simply has to perform a quantum tomographic estimation, by measuring jointly two observables $`X_\lambda `$ and $`X_\lambda ^{}`$ from two quorums $`\{X_\lambda \}`$ and $`\{X_\lambda ^{}\}`$ for the two entangled quantum systems.

### 7.2 An example in the optical domain

In Ref.  it is shown that the proposed method for quantum tomography of a device can be actually performed using joint homodyne tomography on a twin-beam from downconversion of vacuum, with an experimental setup similar to that used in the experiment in Ref. . The feasibility analysis considers, as an example, the experimental determination of the quantum operation corresponding to the unitary displacement operator $`D(z)=e^{za^{}z^{}a}`$. The pertaining matrix $`R(I)`$ is given by

$`R(I)=|D(z)D(z)|,`$ (7.10)
which is the (unnormalizable) eigenstate of the operator $`ab^{}`$ with eigenvalue $`z`$, as shown in Sec. 2.4. As an input bipartite state, one uses the twin-beam from parametric downconversion of Eq. (5.21), which is clearly invertible, since

$`\psi =\sqrt{1|\xi |^2}\xi ^{a^{}a},\psi ^1={\displaystyle \frac{1}{\sqrt{1|\xi |^2)}}}\xi ^{a^{}a}.`$ (7.11)
The experimental apparatus is the same as in the experiment of Ref. , where the twin-beam is provided by a non-degenerate optical parametric amplifier (a KTP crystal) pumped by the second harmonic of a Q-switched mode-locked Nd:YAG laser, which produces a 100-MHz train of 120-ps duration pulses at 1064 nm. The orthogonally polarized twin beams emitted by the KTP crystal (one of which is displaced of $`D(z)`$ by a nearly transparent beam splitter with a strong local oscillator) are separately detected by two balanced homodyne detectors that use two independent local oscillators derived from the same laser. This provides the joint tomography of quadratures $`X_\varphi ^{}X_{\varphi ^{\prime \prime }}`$ needed for the reconstruction. The only experimental problem which still need to be addressed (even though is practically solvable) with respect to the original experiment of Ref.  is the control of the quadrature phases $`\varphi ^{}`$ and $`\varphi ^{\prime \prime }`$ with respect to the LO, which in the original experiment were random.

In Fig. 7.2 the results of a simulated experiment are reported, for displacement parameter $`z=1`$, and for some typical values of the quantum efficiency $`\eta `$ at homodyne detectors and of the total average photon number $`\overline{n}`$ of the twin beam. The diagonal elements $`A_{nn}=n|D(z)|n=[n|n|R_{D(z)}|nn]^{1/2}`$ are plotted for the displacement operator with $`z=1`$. The reconstructed values are shown by thin solid line on an extended abscissa range, with their respective error bars in gray shade, and compared to the theoretical probability (thick solid line). A good reconstruction of the matrix can be achieved in the given range with $`\overline{n}1`$, quantum efficiency as low as $`\eta =0.7`$, and $`10^6÷10^7`$ data. The number of data can be decreased by a factor $`1001000`$ using the tomographic max-likelihood techniques of Ref. , however at the expense of the complexity of the algorithm. Improving quantum efficiency and increasing the amplifier gain (toward a maximally entangled state) have the effect of making statistical errors smaller and more uniform versus the photon labels $`n`$ and $`m`$ of the matrix $`A_{nm}`$.

It is worth emphasizing that the quantum tomographic method of Ref.  for measuring the matrix of a quantum operation can be much improved by means of a max-likelihood strategy aimed at the estimation of some unknown parameters of the quantum operation. In this case, instead of obtaining the matrix elements of $`R(I)`$ from the ensemble averages in Eq. (7.9), one parametrizes $`R(I)`$ in terms of unknown quantities to be experimentally determined, and the likelihood is maximized for the set of experimental data at various randomly selected (tensor) quorum elements, keeping the same fixed bipartite input state. This method is especially useful for a very precise experimental comparison between the characteristics of a given device (e.g. the gain and loss of an active fiber) with those of a quantum standard reference.

## Chapter 8 Maximum-likelihood method in quantum estimation

Quantum estimation of states, observables and parameters is, from very basic principles, matter of statistical inference from a population sampling, and the most comprehensive quantum estimation procedure is quantum tomography. As we have shown in Chapter 3, the expectation value of an operator is obtained by averaging an estimator over the experimental data of a “quorum” of observables. The method is very general and efficient, however, in the averaging procedure, we have fluctuations which result in relatively large statistical errors.

Another relevant strategy, the maximum-likelihood (ML) method, can be used for measuring unknown parameters of transformation on a given state , or for measuring the matrix elements of the density operator itself . The ML strategy  is an entirely different approach to quantum state measurement compared to the standard quantum-tomographic techniques. The ML procedure consists in finding the quantum state, or the value of the parameters, that are most likely to generate the observed data. This idea can be quantified and implemented using the concept of the likelihood functional.

As regards state estimation, the ML method estimates the quantum state as a whole. Such a procedure incorporates a priori knowledge about relations between elements of the density matrix. This guarantees positivity and normalization of matrix, with the result of a substantial reduction of statistical errors. Regarding the estimation of specific parameters, we notice that in many cases the resulting estimators are efficient, unbiased and consistent, thus providing a statistically reliable determination.

As we will show, by using the ML method only small samples of data are required for a precise determination. However, we want to emphasize that such method is not always the optimal solution of the tomographic problem, since it suffers from some major limitations. Besides being biased due to the Hilbert space truncation—even though the bias can be very small if, from other methods, we know where to truncate—it cannot be generalized to the estimation of any ensemble average, but just of a set of parameters from which the density matrix depends. In addition, for increasing number of parameters the method has exponential complexity.

In the following we will review the ML methods proposed in Refs.  and , by deriving the likelihood functional, and applying the ML method to the quantum state reconstruction, with examples for both radiation and spin systems, and, finally, considering the ML estimation for the relevant class of Gaussian states in quantum optics.

### 8.1 Maximum likelihood principle

Here we briefly review the theory of the maximum-likelihood (ML) estimation of a single parameter. The generalization to several parameters, as for example the elements of the density matrix, is straightforward. The only point that should be carefully analyzed is the parameterization of the multidimensional quantity to be estimated. In the next section the specific case of the density matrix will be discussed.

Let $`p(x|\lambda )`$ the probability density of a random variable $`x`$, conditioned to the value of the parameter $`\lambda `$. The form of $`p`$ is known, but the true value of $`\lambda `$ is unknown, and will be estimated from the result of a measurement of $`x`$. Let $`x_1,x_2,\mathrm{},x_N`$ be a random sample of size $`N`$. The joint probability density of the independent random variable $`x_1,x_2,\mathrm{},x_N`$ (the global probability of the sample) is given by

$`(x_1,x_2,\mathrm{},x_N|\lambda )=\mathrm{\Pi }_{k=1}^Np(x_k|\lambda ),`$ (8.1)
and is called the likelihood function of the given data sample (hereafter we will suppress the dependence of $``$ on the data). The maximum-likelihood estimator (MLE) of the parameter $`\lambda `$ is defined as the quantity $`\lambda _{ml}\lambda _{ml}(\{x_k\})`$ that maximizes $`(\lambda )`$ for variations of $`\lambda `$, namely $`\lambda _{ml}`$ is given by the solution of the equations

$`{\displaystyle \frac{(\lambda )}{\lambda }}=0;{\displaystyle \frac{^2(\lambda )}{\lambda ^2}}<0.`$ (8.2)
The first equation is equivalent to $`L/\lambda =0`$ where

$`L(\lambda )=\mathrm{log}(\lambda )={\displaystyle \underset{k=1}{\overset{N}{}}}\mathrm{log}p(x_k|\lambda )`$ (8.3)
is the so-called log-likelihood function.

In order to obtain a measure for the confidence interval in the determination of $`\lambda _{ml}`$ we consider the variance

$`\sigma _\lambda ^2={\displaystyle \left[\underset{k}{}dx_kp(x_k|\lambda )\right]\left[\lambda _{ml}(\{x_k\})\lambda \right]^2}.`$ (8.4)
In terms of the Fisher information

$`F={\displaystyle 𝑑x\left[\frac{p(x|\lambda )}{\lambda }\right]^2\frac{1}{p(x|\lambda )}},`$ (8.5)
it is easy to prove that

$`\sigma _\lambda ^2{\displaystyle \frac{1}{NF}},`$ (8.6)
where $`N`$ is the number of measurements. The inequality in Eq. (8.6) is known as the Cramér-Rao bound  on the precision of the ML estimation. Notice that this bound holds for any functional form of the probability distribution $`p(x|\lambda )`$, provided that the Fisher information exists $`\lambda `$ and $`_\lambda p(x|\lambda )`$ exists $`x`$. When an experiment has ”good statistics” (i.e. for a large enough data sample) the Cramér-Rao bound is saturated.

### 8.2 ML quantum state estimation

In this section we review the method of the maximum likelihood estimation of the quantum state of Ref. , focusing attention to the cases of homodyne and spin tomography.

We consider an experiment consisting of $`N`$ measurements performed on identically prepared copies of a given quantum system. Each measurement is described by a positive operator-valued measure (POVM). The outcome of the $`i`$th measurement corresponds to the realization of a specific element of the POVM used in the corresponding run, and we denote this element by $`\mathrm{\Pi }_i`$. The likelihood is here a functional of the density matrix $`(\rho )`$ and is given by the product

$$(\rho )=\underset{i=1}{\overset{N}{}}\text{Tr}(\rho \mathrm{\Pi }_i),$$
(8.7)
which represents the probability of the observed data. The unknown element of the above expression, which we want to infer from data, is the density matrix describing the measured ensemble. The estimation strategy of the ML technique is to maximize the likelihood functional over the set of the density matrices. Several properties of the likelihood functional are easily found, if we restrict ourselves to finite dimensional Hilbert spaces. In this case, it can be easily proved that $`(\rho )`$ is a concave function defined on a convex and closed set of density matrices. Therefore, its maximum is achieved either on a single isolated point, or on a convex subset of density matrices. In the latter case, the experimental data are insufficient to provide a unique estimate for the density matrix using the ML strategy. On the other hand, the existence of a single maximum allows us to assign unambiguously the ML estimate for the density matrix.

The ML estimation of the quantum state, despite its elegant general formulation, results in a highly nontrivial constrained optimization problem, even if we resort to purely numerical means. The main difficulty lies in the appropriate parameterization of the set of all density matrices. The parameter space should be of the minimum dimension in order to preserve the maximum of the likelihood function as a single isolated point. Additionally, the expression of quantum expectation values in terms of this parameterization should enable fast evaluation of the likelihood function, as this step is performed many times in the course of numerical maximization.

For such purpose one introduces  a parameterization of the set of density matrices which provides an efficient algorithm for maximization of the likelihood function. We represent the density matrix in the form

$$\rho =T^{}T,$$
(8.8)
which automatically guarantees that $`\rho `$ is positive and Hermitian. The remaining condition of unit trace $`\text{Tr}\rho =1`$ will be taken into account using the method of Lagrange multipliers. In order to achieve the minimal parameterization, we assume that $`T`$ is a complex lower triangular matrix, with real elements on the diagonal. This form of $`T`$ is motivated by the Cholesky decomposition known in numerical analysis  for arbitrary non negative Hermitian matrix. For an $`M`$-dimensional Hilbert space, the number of real parameters in the matrix $`T`$ is $`M+2M(M1)/2=M^2`$, which equals the number of independent real parameters for a Hermitian matrix. This confirms that such parameterization is minimal, up to the unit trace condition.

In numerical calculations, it is convenient to replace the likelihood functional by its natural logarithm, which of course does not change the location of the maximum. Thus the log-likelihood function subjected to numerical maximization is given by

$$L(T)=\underset{i=1}{\overset{N}{}}\mathrm{ln}\text{Tr}(T^{}T\mathrm{\Pi }_i)\lambda \text{Tr}(T^{}T),$$
(8.9)
where $`\lambda `$ is a Lagrange multiplier accounting for normalization of $`\rho `$. Writing $`\rho `$ in terms of its eigenvectors $`|\psi _\mu `$ as $`\rho =_\mu y_\mu ^2|\psi _\mu \psi _\mu |`$, with real $`y_\mu `$, the maximum likelihood condition $`L/y_\nu =0`$ reads

$`\lambda y_\nu ={\displaystyle \underset{i=1}{\overset{N}{}}}[y_\nu \psi _\nu |\mathrm{\Pi }_i|\psi _\nu /\text{Tr}(\rho \mathrm{\Pi }_i)],`$ (8.10)
which, after multiplication by $`y_\nu `$ and summation over $`\nu `$, yields $`\lambda =N`$. The Lagrange multiplier then equals the total number of measurements $`N`$.

This formulation of the maximization problem allows one to apply standard numerical procedures for searching the maximum over the $`M^2`$ real parameters of the matrix $`T`$. The examples presented below use the downhill simplex method .

The first example is the ML estimation of a single-mode radiation field. The experimental apparatus used in this technique is the homodyne detector. According to Sec. 2.4 the homodyne measurement is described by the positive operator-valued measure

$$(x;\phi )=\sqrt{\frac{2\eta }{\pi (1\eta )}}\mathrm{exp}\left[\frac{2\eta }{1\eta }(X_\phi x)^2\right],$$
(8.11)
where $`\eta `$ is the detector efficiency, and $`X_\phi =(a^{}e^{i\phi }+ae^{i\phi })/2`$ is the quadrature operator at phase $`\phi `$.

After $`N`$ measurements, we obtain a set of pairs $`(x_i;\phi _i)`$, where $`i=1,\mathrm{},N`$. The log-likelihood functional is given by Eq. (8.9) with $`\mathrm{\Pi }_i(x_i;\phi _i)`$. Of course, for a light mode it is necessary to truncate the Hilbert space to a finite dimensional basis. We shall assume that the highest Fock state has $`M1`$ photons, i.e. that the dimension of the truncated Hilbert space is $`M`$. For the expectation $`\text{Tr}[T^{}T(x;\phi )]`$ it is necessary to use an expression which is explicitly positive, in order to protect the algorithm against occurrence of small negative numerical arguments of the logarithm function. A simple derivation yields

$`\text{Tr}[T^{}T(x;\phi )]=\sqrt{\eta }{\displaystyle \underset{k=0}{\overset{M1}{}}}{\displaystyle \underset{j=0}{\overset{k}{}}}\left|{\displaystyle \underset{n=0}{\overset{kj}{}}}k|T|n+jB_{n+j,n}n|\sqrt{\eta }xe^{in\phi }\right|^2,`$ (8.12)
where

$`B_{n+j,n}=\left[\left({\displaystyle \genfrac{}{}{0pt}{}{n+j}{n}}\right)\eta ^n(1\eta )^j\right]^{1/2},`$ (8.13)
and

$`n|x=\left({\displaystyle \frac{2}{\pi }}\right)^{1/4}{\displaystyle \frac{1}{\sqrt{2^nn!}}}H_n(\sqrt{2}x)\mathrm{exp}(x^2)`$ (8.14)
are the eigenstates of the harmonic oscillator in the position representation—$`H_n(x)`$ being the $`n`$th Hermite polynomial.

The ML technique can be applied to reconstruct the density matrix in the Fock basis from Monte Carlo simulated homodyne statistics. Fig. 8.1 depicts the matrix elements of the density operator as obtained for a coherent state and a squeezed vacuum, respectively. Remarkably, only 50000 homodyne data have been used for quantum efficiency $`\eta =80\%`$. We recall that in quantum homodyne tomography the statistical errors are known to grow rapidly with decreasing efficiency $`\eta `$ of the detector . In contrast, the elements of the density matrix reconstructed using the ML approach remain bounded, as the whole matrix must satisfy positivity and normalization constraints. This results in much smaller statistical errors. As a comparison one could see that the same precision of the reconstructions in Fig. 8.1 could be achieved using $`10^7`$$`10^8`$ data samples with conventional quantum tomography. On the other hand, in order to find numerically the ML estimate we need to set a priori the cut-off parameter for the photon number, and its value is limited by increasing computation time.

Another relevant example is the reconstruction of the quantum state of two-mode field using single-LO homodyning of Chapter 5. Here, the full joint density matrix can be measured by scanning the quadratures of all possible linear combinations of modes. For two modes the measured quadrature operator is given by

$`X(\theta ,\psi _0,\psi _1)={\displaystyle \frac{1}{2}}(ae^{i\psi _0}\mathrm{cos}\theta +be^{i\psi _1}\mathrm{sin}\theta +\text{h.c.}),`$ (8.15)
where $`(\theta ,\psi _0,\psi _1)S^2\times [0,2\pi ]`$, $`S^2`$ being the Poincaré sphere and one phase ranging between $`0`$ and $`2\pi `$. In each run these parameters are chosen randomly. The POVM describing the measurement is given by the right-hand side of Eq. (8.11), with $`X_\phi `$ replaced by $`X(\theta ,\psi _0,\psi _1)`$. An experiment for the two orthogonal states $`|\mathrm{\Psi }_1=(|00+|11)/\sqrt{2}`$ and $`|\mathrm{\Psi }_2=(|01+|10)/\sqrt{2}`$ has been simulated, in order to reconstruct the density matrix in the two-mode Fock basis using the ML technique. The results are reported in Fig. 8.2.

The ML procedure can also be applied for reconstructing the density matrix of spin systems. For example, let us consider $`N`$ repeated preparations of a pair of spin-1/2 particles. The particles are shared by two parties. In each run, the parties select randomly and independently from each other a direction along which they perform a spin measurement. The obtained result is described by the joint projection operator (spin coherent states ) $`_i=|\mathrm{\Omega }_i^A,\mathrm{\Omega }_i^B\mathrm{\Omega }_i^A,\mathrm{\Omega }_i^B|`$, where $`\mathrm{\Omega }_i^A`$ and $`\mathrm{\Omega }_i^B`$ are the vectors on the Bloch sphere corresponding to the outcomes of the $`i`$th run, and the indices $`A`$ and $`B`$ refer to the two particles. As in the previous examples, it is convenient to use an expression for the quantum expectation value $`\text{Tr}(T^{}T_i`$) which is explicitly positive. The suitable form is

$$\text{Tr}(T^{}T_i)=\underset{\mu }{}|\mu |T|\mathrm{\Omega }_i^A,\mathrm{\Omega }_i^B|^2,$$
(8.16)
where $`|\mu `$ is an orthonormal basis in the Hilbert space of the two particles. The result of a simulated experiment with only 500 data for the reconstruction of the density matrix of the singlet state is shown in Fig. 8.3.

Summarizing, the ML technique can be used to estimate the density matrix of a quantum system. With respect to conventional quantum tomography this method has the great advantage of needing much smaller experimental samples, making experiments with low data rates feasible, however with a truncation of the Hilbert space dimension. We have shown that the method is general and the algorithm has solid methodological background, its reliability being confirmed in a number of Monte Carlo simulations. However, for increasing dimension of Hilbert spaces the method has exponential complexity.

### 8.3 Gaussian-state estimation

In this Section we review the ML determination method of Ref.  for the parameters of Gaussian states. Such states represent the wide class of coherent, squeezed and thermal states, all of them being characterized by a Gaussian Wigner function. Apart from an irrelevant phase, we consider Wigner functions of the form

$`W(x,y)={\displaystyle \frac{2\mathrm{\Delta }^2}{\pi }}\mathrm{exp}\left\{2\mathrm{\Delta }^2\left[e^{2r}(x\text{Re}\mu )^2+e^{2r}(y\text{Im}\mu )^2\right]\right\},`$ (8.17)
and the ML technique with homodyne detection is applied to estimate the four real parameters $`\mathrm{\Delta },r,\text{Re}\mu `$ and $`\text{Im}\mu `$. The four parameters provide the number of thermal, squeezing and coherent-signal photons in the quantum state as follows

$`n_{th}={\displaystyle \frac{1}{2}}\left({\displaystyle \frac{1}{\mathrm{\Delta }^2}}1\right),`$
$`n_{sq}=\mathrm{sinh}^2r,`$
$`n_{coh}=|\mu |^2.`$ (8.18)
The density matrix $`\rho `$ corresponding to the Wigner function in Eq. (8.17) writes

$`\rho =D(\mu )S(r){\displaystyle \frac{1}{n_{th}+1}}\left({\displaystyle \frac{n_{th}}{n_{th}+1}}\right)^{a^{}a}S^{}(r)D^{}(\mu ),`$ (8.19)
where $`S(r)=\mathrm{exp}[r(a^2a^2)/2]`$ and $`D(\mu )=\mathrm{exp}(\mu a^{}\mu ^{}a)`$ denote the squeezing and displacement operators, respectively.

The theoretical homodyne probability distribution at phase $`\phi `$ with respect to the local oscillator can be evaluated using Eq. (2.10), and is given by the Gaussian

$`p(x,\phi )`$ $`=`$ $`\sqrt{{\displaystyle \frac{2\mathrm{\Delta }^2}{\pi (e^{2r}\mathrm{cos}^2\phi +e^{2r}\mathrm{sin}^2\phi )}}}`$ (8.20)
$`\times `$ $`\mathrm{exp}\left\{{\displaystyle \frac{2\mathrm{\Delta }^2}{e^{2r}\mathrm{cos}^2\phi +e^{2r}\mathrm{sin}^2\phi }}\left[x\text{Re}(\mu e^{i\phi })\right]^2\right\}.`$
The log-likelihood function (8.3) for a set of $`N`$ homodyne outcomes $`x_i`$ at random phase $`\phi _i`$ then writes as follows

$`L`$ $`=`$ $`{\displaystyle \underset{i=1}{\overset{N}{}}}{\displaystyle \frac{1}{2}}\mathrm{log}{\displaystyle \frac{2\mathrm{\Delta }^2}{\pi (e^{2r}\mathrm{cos}^2\phi _i+e^{2r}\mathrm{sin}^2\phi _i)}}`$ (8.21)
$``$ $`{\displaystyle \frac{2\mathrm{\Delta }^2}{e^{2r}\mathrm{cos}^2\phi _i+e^{2r}\mathrm{sin}^2\phi _i}}\left[x_i\text{Re}(\mu e^{i\phi _i})\right]^2.`$
The ML estimators $`\mathrm{\Delta }_{ml},r_{ml},\text{Re}\mu _{ml}`$ and $`\text{Im}\mu _{ml}`$ are found upon maximizing Eq. (8.21) versus $`\mathrm{\Delta },r,\text{Re}\mu `$ and $`\text{Im}\mu `$.

In order to evaluate globally the state reconstruction, one considers the normalized overlap $`𝒪`$ between the theoretical and the estimated state

$`𝒪={\displaystyle \frac{\text{Tr}[\rho \rho _{ml}]}{\sqrt{\text{Tr}[\rho ^2]\text{Tr}[\rho _{ml}^2]}}}.`$ (8.22)
Notice that $`𝒪=1`$ iff $`\rho =\rho _{ml}`$. Through Monte-Carlo simulations, one always finds a value around unity, typically with statistical fluctuations over the third digit, for number of data samples $`N=50000`$, quantum efficiency at homodyne detectors $`\eta =80\%`$, and state parameters with the following ranges: $`n_{th}<3`$, $`n_{coh}<5`$, and $`n_{sq}<3`$. Also with such a small number of data samples, the quality of the state reconstruction is so good that other physical quantities that are theoretically evaluated from the experimental values of $`\mathrm{\Delta }_{ml},r_{ml},\text{Re}\mu _{ml}`$ and $`\text{Im}\mu _{ml}`$ are inferred very precisely. For example, in Ref.  the photon number probability of a squeezed thermal state has been evaluated, which is given by the integral

$`n|\rho |n={\displaystyle _0^{2\pi }}{\displaystyle \frac{d\varphi }{2\pi }}{\displaystyle \frac{[C(\varphi ,n_{th},r)1]^n}{C(\varphi ,n_{th},r)^{n+1}}},`$ (8.23)
with $`C(\varphi ,n_{th},r)=(n_{th}+\frac{1}{2})(e^{2r}\mathrm{sin}^2\varphi +e^{2r}\mathrm{cos}^2\varphi )+\frac{1}{2}`$. The comparison of the theoretical and the experimental results for a state with $`n_{th}=0.1`$ and $`n_{sq}=3`$ is reported in Fig. 8.4. The statistical error of the reconstructed number probability affects the third decimal digit, and is not visible on the scale of the plot.

The estimation of parameters of Gaussian Wigner functions through the ML method allows one to estimate the parameters in quadratic Hamiltonians of the generic form

$`H=\alpha a+\alpha ^{}a^{}+\phi a^{}a+{\displaystyle \frac{1}{2}}\xi a^2+{\displaystyle \frac{1}{2}}\xi ^{}a^2.`$ (8.24)
In fact, the unitary evolution operator $`U=e^{iHt}`$ preserves the Gaussian form of an input state with Gaussian Wigner function. In other words, one can use a known Gaussian state to probe and characterize an optical device described by a Hamiltonian as in Eq. (8.24). Assuming $`t=1`$ without loss of generality, the Heisenberg evolution of the radiation mode $`a`$ is given by

$`U^{}aU=\gamma a+\delta a^{}+\mu ,`$ (8.25)
with

$`\gamma =\mathrm{cos}(\sqrt{\phi ^2|\xi |^2})i{\displaystyle \frac{\phi }{\sqrt{\phi ^2|\xi |^2}}}\mathrm{sin}(\sqrt{\phi ^2|\xi |^2}),`$ (8.26)
$`\delta =i{\displaystyle \frac{\xi ^{}}{\sqrt{\phi ^2|\xi |^2}}}\mathrm{sin}(\sqrt{\phi ^2|\xi |^2}),`$
$`\mu ={\displaystyle \frac{\phi \alpha ^{}\xi ^{}\alpha }{\phi ^2|\xi |^2}}(\mathrm{cos}(\sqrt{\phi ^2|\xi |^2})1)i{\displaystyle \frac{\alpha ^{}}{\sqrt{\phi ^2|\xi |^2}}}\mathrm{sin}(\sqrt{\phi ^2|\xi |^2}).`$
For an input state $`\rho `$ with known Wigner function $`W_\rho (\beta ,\beta ^{})`$, the corresponding output Wigner function writes

$`W_{U\rho U^{}}(\beta ,\beta ^{})=W_\rho [(\beta \mu )\gamma ^{}(\beta ^{}\mu ^{})\delta ,(\beta ^{}\mu ^{})\gamma (\beta \mu )\delta ^{}].`$ (8.27)
Hence, by estimating the parameters $`\gamma ,\delta ,\mu `$ and inverting Eqs. (8.26), one obtains the ML values for $`\alpha ,\phi `$, and $`\xi `$ of the Hamiltonian in Eq. (8.24). The present example can be used in practical applications for the estimation of the gain of a phase-sensitive amplifier or equivalently to estimate a squeezing parameter.

## Chapter 9 Classical imaging by quantum tomography

As we showed in Chapter 2, the development of quantum tomography has its origin in the inadequacy of classical imaging procedures to face the quantum problem of Wigner function reconstruction. In this chapter we briefly illustrate how to go back to classical imaging and profitably use quantum tomography as a tool for image reconstruction and compression: this is the method of Fictitious Photons Tomography of Ref. .

The problem of tomographic imaging is to recover a mass distribution $`m(x,y)`$ in a $`2`$-d slab from a finite collection of one dimensional projections. The situation is schematically sketched in Fig. 9.1a where $`m(x,y)`$ describes two circular holes in a uniform background. The tomographic machine, say an X-ray equipment, collects many stripe photos of the sample from various directions $`\theta `$, and then numerically performs a mathematical transformation—the so called inverse Radon transform —in order to reconstruct $`m(x,y)`$ from its radial profiles at different $`\theta `$’s. The problem which is of interest for us is when the radial profiles are not well-resolved digitalized functions, but actually represent the density distribution of random points, as if in our X-ray machine the beam is so weak that radial photos are just the collection of many small spots, each from a single X-ray photon (this situation is sketched in Fig. 9.1b). It is obvious that this case can be reduced to the previous one by counting all points falling in a predetermined $`1`$-d mesh, and giving radial profiles in form of histograms (this is what actually happens in a real machine, using arrays of photodetectors). However, we want to use the whole available information from each ”event”—i.e. the exact $`1`$-d location of each spot—in a way which is independent on any predetermined mesh. In practice, this situation occurs when the signal is so weak and the machine resolution is so high (i.e. the mesh-step is so tiny) that only zero or one photon at most can be collected in each channel. As we will see, this low-signal/high-resolution case naturally brings the imaging problem into the domain of quantum tomography. Images are identified with Wigner functions, such to obtain a description in terms of density matrices. These are still trace-class matrices (corresponding to ”normalizable” images), but are no longer positive definite, because an “image” generally is not a genuine Wigner function and violates the Heisenberg relations on the complex plane (the phase space of a single mode of radiation). Hence, such density matrices are unphysical: they are just a mathematical tool for imaging. This is the reason why this method has been named Fictitious Photons Tomography . As we will see in the following, the image resolution improves by increasing the rank of the density matrix, and in this way the present method also provides a new algorithm for image compression, which is suited to angular image scanning.

### 9.1 From classical to quantum imaging

We adopt the complex notation, with $`\alpha =x+iy`$ representing a point in the image plane. In this way $`\alpha `$ and $`\alpha ^{}`$ are considered as independent variables, and the $`2`$-d image—here denoted by the same symbol $`W(\alpha ,\alpha ^{})`$ used for the Wigner function—is just a generic real function of the point in the plane. In the most general situation $`W(\alpha ,\alpha ^{})`$ is defined on the whole complex plane, where it is normalized to some finite constant, and it is bounded from both below and above, with range representing the darkness nuance. For X-ray tomography $`W(\alpha ,\alpha ^{})`$ roughly represents the absorption coefficient as a function of the point $`\alpha `$. We consider a linear absorption regime, i.e. the image extension is negligible with respect to the radiation absorption length in the medium. At the same time we neglect any diffraction effect.

As shown in Sec. 3.2 the customary imaging technique is based on the inverse Radon transform. A tomography of a two dimensional image $`W(\alpha ,\alpha ^{})`$ is a collection of one dimensional projections $`p(x,\theta )`$ at different values of the observation angle $`\theta `$. We rewrite here the definition of the Radon transform of $`W(\alpha ,\alpha ^{})`$

$`p(x,\theta )={\displaystyle _{\mathrm{}}^+\mathrm{}}{\displaystyle \frac{dy}{\pi }}W((x+iy)e^{i\theta },(xiy)e^{i\theta }).`$ (9.1)
In Eq. (9.1) $`x`$ is the current coordinate along the direction orthogonal to the projection and $`y`$ is the coordinate along the projection direction. The situation is depicted in Fig. 9.1 where $`W(\alpha ,\alpha ^{})`$ is plotted along with its $`p(x,\theta )`$ profiles for $`\theta =0,\pi /2`$ for a couple of identical circular holes that are symmetrically disposed with respect to the origin.

The reconstruction of the image $`W(\alpha ,\alpha ^{})`$ from its projections $`p(x,\theta )`$—also called “back projection”—is given by the inverse Radon transform, which, following the derivation in Sec. 3.2, leads to the filtering procedure

$`W(\alpha ,\alpha ^{})={\displaystyle _0^\pi }{\displaystyle \frac{d\theta }{2\pi }}𝒫{\displaystyle _{\mathrm{}}^+\mathrm{}}𝑑x{\displaystyle \frac{p(x,\theta )/x}{x\alpha _\theta }},`$ (9.2)
where $`𝒫`$ denotes the Cauchy principal value and $`\alpha _\theta =\text{Re}(\alpha e^{i\theta })`$. Eq. (9.2) is commonly used in conventional tomographic imaging (see, for example, Ref. ).

Let us now critically consider the above procedure in the case of very weak signals, namely when $`p(x,\theta )`$ just represents the probability distribution of random X-ray spots on a fine-mesh multichannel: this situation is sketched in Fig. 9.1b. From Eq. (9.2) one can recover $`W(\alpha ,\alpha ^{})`$ only when the analytical form of $`p(x,\theta )`$ is known. But the experimental outcomes of each projection actually are random data distributed according to $`p(x,\theta )`$, whereas in order to recover $`W(\alpha ,\alpha ^{})`$ from Eq. (9.2) one has to evaluate the first order derivatives of $`p(x,\theta )`$. The need of the analytical form for projections $`p(x,\theta )`$ requires a filtering procedure on data, usually obtained by “splining” data in order to use Eq. (9.2).

The above procedure leads to approximate image reconstructions, and the choice of any kind of smoothing parameter unavoidably affects in a systematic way the statistics of errors. In the following we show how quantum tomography can be used for conventional imaging in presence of weak signals, providing both ideally controlled resolution and reliable error statistics.

The basic formula we will use is the expansion of the Wigner function in the number representation of Eqs. (6.24) and (6.25). In practice, the Hilbert space has to be truncated at some finite dimension $`d_{}`$, and this sets the resolution for the reconstruction of $`W(\alpha ,\alpha ^{})`$. However, as we will show, this resolution can be chosen at will, independently on the number of experimental data.

As previously noticed, in general an image does not correspond to a Wigner function of a physical state, due to the fact that the Heisenberg relations unavoidably produce only smooth Wigner functions, whereas a conventional image can have very sharp edges. However, if one allows the density matrix to be no longer positive definite (but still trace class), a correspondence with images is obtained, which holds in general. In this way every image is stored into a trace-class matrix $`\rho _{n,m}`$ via quantum tomography, and a convenient truncation of the matrix dimension $`d_{}`$ can be chosen.

The connection between images and matrices is the main point of this approach: the information needed to reconstruct the image is stored in a $`d_{}\times d_{}`$ matrix. For suitably chosen dimension $`d_{}`$ the present method can also provide a procedure for image compression. Notice that the correspondence between images and trace-class matrices retains some symmetries of the image, which manifest as algebraic properties of the matrix $`\rho _{n,m}`$. For example an isotropic image (like a uniform circle centered at the origin) is stored in a diagonal matrix. Other symmetries are given in Tab. 9.1.

The truncated Hilbert space dimension $`d_{}`$ sets the imaging resolution. The kind of resolution can be understood by studying the behavior of the kernels $`[|n+dn|](x,\theta )`$ of Eq. (3.44), which are averaged over the experimental data in order to obtain the matrix elements $`\rho _{n,n+d}`$. Outside a region that is almost independent of $`n`$ and $`d`$, all functions $`[|n+dn|](x,\theta )`$ decrease exponentially, whereas inside this region they oscillate with a number of oscillations linearly increasing with $`2n+d`$. This behavior produces the effects illustrated in Fig. 9.2, where we report the tomographic reconstruction of the font “a” for increasing dimension $`d_{}`$. The plot is obtained by numerically integrating the kernel functions from given analytic transmission profiles $`p(x,\theta )`$. As we see from Fig. 9.2 both the radial and the angular resolutions improve versus $`d_{}`$, making the details of the image sharper and sharper already form a relatively small truncation $`d_{}=48`$.

A quantitative measure of the precision of the tomographic reconstruction can be given in terms of the distance $`D`$ between the true and the reconstructed image, which, in turn, coincides with the Hilbert distance $`D`$ between the corresponding density matrices. One has

$`D`$ $`=`$ $`{\displaystyle d^2\alpha |\mathrm{\Delta }W(\alpha ,\alpha ^{})|^2}=\text{Tr}(\mathrm{\Delta }\rho )^2`$ (9.3)
$`=`$ $`{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}\mathrm{\Delta }\rho _{n,n}^2+2{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}{\displaystyle \underset{\lambda =1}{\overset{\mathrm{}}{}}}\left|\mathrm{\Delta }\rho _{n,n+\lambda }^2\right|^2,`$
where $`\mathrm{\Delta }[\mathrm{}]=[\mathrm{}]_{true}[\mathrm{}]_{reconstructed}`$. The convergence of $`D`$ versus $`d_{}`$ is given in Fig. 9.3 for a solid circle of unit radius centered at the origin. In this case the obtained density matrix has only diagonal elements, according to Tab. 9.1. These are given by

$`\rho _{n,n}=2{\displaystyle \underset{\nu =0}{\overset{n}{}}}(2)^\nu \left({\displaystyle \genfrac{}{}{0pt}{}{n}{\nu }}\right)\mathrm{\Phi }(1\nu ,2,2R^2),`$ (9.4)
where $`\mathrm{\Phi }(\alpha ,\beta ,z)`$ denotes the confluent hypergeometric function of argument $`z`$ and parameters $`\alpha `$ and $`\beta `$.

Insofar we have analyzed the method only on the basis of given analytic profiles $`p(x,\theta )`$. As already said, however, the method is particularly advantageous in the weak-signal/high-resolution situation, where the imaging can be achieved directly from averaging the kernel functions on data. In this case the procedure allows to exploit the whole available experimental resolution, whereas the image resolution is set at will. In Fig. 9.4 we report a Monte Carlo simulation of an experimental tomographic reconstruction of the font “a” for increasing number of data. All plots are obtained at the maximum available dimension $`d_{}=48`$, and using $`F=100`$ scanning phases. The situation occurring for small numbers of data is given in the first plot, where the the highly resolved image still exhibits the natural statistical fluctuations due to the limited number of data. For larger sample the image appears sharper from the random background, and it is clearly recognizable for a number of data equal to $`10^6`$. The method is efficient also from the computational point of view, as the time needed for image reconstruction is quadratic in the number of elements of the density matrix, and linear in the number of experimental data. Needless to say, imaging by quantum homodyne tomography is at the very early stages and further investigation is in order.