Optical homodyne tomography with polynomial series expansion
aa r X i v : . [ qu a n t - ph ] N ov Optical homodyne tomography with polynomial series expansion
Hugo Benichi ∗ and Akira Furusawa Department of Applied Physics, The University of Tokyo, Tokyo, Japan (Dated: October 23, 2018)We present and demonstrate a method for optical homodyne tomography based on the inverse Radon trans-form. Different from the usual filtered back-projection algorithm, this method uses an appropriate polynomialseries to expand the Wigner function and the marginal distribution and discretize Fourier space. We show thatthis technique solves most technical difficulties encountered with kernel deconvolution based methods and re-constructs overall better and smoother Wigner functions. We also give estimators of the reconstruction errorsfor both methods and show improvement in noise handling properties and resilience to statistical errors.
I. INTRODUCTION
In quantum mechanics it is not possible to directly observea quantum state | ψ i . In order to obtain full knowledge about | ψ i it is necessary to accumulate measurement statistics ofobservables, such as position ˆ x or momentum ˆ p , on many dif-ferent bases. In quantum optics, this statistical measurementcan be achieved by angle resolved homodyne measurementof the operator ˆ x θ = ˆ x cos θ + ˆ p sin θ to acquire statistics ofthe squared modulus of the wave function |h x θ | ψ i| . Insteadof the quantum state | ψ i , one is rather usually interested inreconstructing the more general density matrix ˆ ρ of the sys-tem. Fully equivalent to ˆ ρ , it is also possible to reconstructthe Wigner function W ( q, p ) from |h x θ | ψ i| . However, thereconstruction of ˆ ρ or W ( q, p ) is not immediate and requiresthe reconstruction of the complex phase of the quantum sys-tem from the many angle resolved measurements. With themeasurement of |h x θ | ψ i| , these two operations together arereferred to as quantum homodyne tomography or optical ho-modyne tomography [1].While some tomography algorithms reconstruct the for-mer density matrix, others rather reconstruct the latter Wignerfunction. Independently, tomography algorithms can beroughly classified into two species. Historically the first to beproposed and used for optical homodyne tomography, linearmethods exploit and inverse the linear relationship betweenthe experimentally measurable quantity |h x θ | ψ i| on one handand ˆ ρ or W ( q, p ) on the other hand. Among them, the filteredback-projection algorithm [1, 2] based on the inverse Radontransform [3] is the most commonly used. Similar in nature,there also exist methods based on quantum state sampling ofindividual components of the density matrix ˆ ρ with samplefunctions [4, 5]. The linear methods, however, suffer in gen-eral from technical difficulties associated with the numericaldeconvolution necessary to perform the linear inversion of theRadon transform (see Sec. II for details). In addition, theyusually do not guarantee the physicality of the reconstructedstate, the positivity of ˆ ρ . Finally they perform weakly againststatistical noise and show numerical instabilities for higherfrequency components and fine details of the reconstructedobjects. Variational methods, such as the maximum entropy[6] and maximum likelihood [7] algorithms, were latter ap- ∗ [email protected] plied to optical homodyne tomography to address these prob-lems. These methods can be designed to enforce the physi-cality of the reconstructed state and are usually more resilientto statistical errors. Since the reconstructed states are not de-fined constructively, an approximation procedure, typically it-erative, is used to achieve the reconstruction in practice [8].Notice that in theory it is actually possible to bypass thesenumerical reconstructions and directly observe the Wignerfunction W ( q, p ) with repeated measures of the parity oper-ator ˆ P = e iπ ˆ n where ˆ n is the number operator [9]. This mea-surement technique uses the link between the Wigner functionvalue at point ( q, p ) and the expectation value of ˆ P for the dis-placed density matrix ˆ ρW ( q, p ) = 2 π tr h ˆ D ( − α )ˆ ρ ˆ D ( α ) e iπ ˆ n i , (1.1)where ˆ D is the displacement operator and α = ( q + ip ) / √ . Aclose tomography technique has been experimentally demon-strated in coupled systems of atoms and light [10]. Unfortu-nately, a parity detector is a highly non-linear detector whichcan only be partially implemented for light beams with time-multiplexing and single photon detectors. Therefore with cur-rent state-of-the-art technologies in quantum optics, it is notpossible to rely on count statistics alone for quantum state to-mography and one has to use optical homodyne tomographybased on Gaussian measurements.While the linear methods look inferior to the variationalmethods, most of their associated problems are only technicalin nature and can in principle be solved. In this paper we showthat is it possible to use a linear reconstruction algorithm withbetter resilience to noise and better physical properties overallthan the usual filtered back-projection method. The successof this approach lies in a systematic expansion of both theWigner function W ( q, p ) and the marginal distribution p ( x, θ ) in polar coordinates. This circular harmonic expansion tech-nique has been applied in the past to other problems wherethe Radon transform plays a role in tomography [11, 12], andhere we adapt it to the quantum framework of optical homo-dyne tomography. In Sec. II we first review the basics of theinverse Radon transform and the usual filtered back-projectionalgorithms for optical homodyne tomography. In Sec. III weintroduce the expansion method: we first conduct a spectralanalysis of the angular components of p ( x, θ ) and W ( q, p ) ;from this analysis we argue that a polynomial approximationis an efficient way to expand the radial components. In Sec.IV we give details about the implementation of the algorithmand also provide an estimator of the reconstruction errors. Us-ing our estimator we study the performances relatively to thefiltered back-projection algorithm on simulated and experi-mental data sets. We complete this comparison with numer-ical studies of the distance between target and reconstructedquantum states. II. FILTERED BACK-PROJECTION
In 1917, Radon introduces the integral transform R of two-dimensional functions integrated along straight lines and pro-vides the formula for the inverse transform R − [3]. Todaythe Radon and inverse Radon transforms are ubiquitous in to-mography and find applications in many different area of sci-ence. The Radon transform is as well applicable to opticalhomodyne tomography. First we recall the definition of theobservable operator ˆ x θ of an homodyne measurement, ˆ x θ = ˆ U † θ ˆ x ˆ U θ = ˆ x cos θ + ˆ p sin θ, (2.1)where ˆ U θ is the rotation operator in phase space, or phase-shifting operator. The marginal distribution of the homodynecurrent p ( x, θ ) is then distributed according to the squaredmodulus of the wave function p ( x, θ ) = |h x θ | ψ i| = h x | ˆ U θ | ψ ih ψ | ˆ U † θ | x i , (2.2)where | x θ i is the eigenvector of ˆ x θ . The Radon transform R links the Wigner function W ( q, p ) of the quantum state | ψ i and p ( x, θ ) the marginal distribution of the homodyne currentwith a projection of W ( q, p ) on a particular angle of observa-tion θ [13] p ( x, θ ) = R ( W )= Z Z R W ( q, p ) δ ( x − q cos θ − p sin θ ) dqdp = Z + ∞−∞ W ( x cos θ − p sin θ, x sin θ + p cos θ ) dp. (2.3)In his original paper, Radon mathematically inverses his trans-form with the back-projection B of the derivative of theHilbert transform H of p ( x, θ ) W ( q, p ) = 12 π B (cid:18) ∂∂y H ( p ( x, θ ))( y ) . (cid:19) , (2.4)where the back-projection operator B of a function f ( x, θ ) isthe function F ( q, p ) defined by F ( q, p ) = Z π f ( q cos θ + p sin θ, θ ) dθ. (2.5)Expanding Eq. (2.4) we obtain the inversion formula W ( q, p ) = − P π Z π Z + ∞−∞ p ( x, θ )( q cos θ + p sin θ − x ) dxdθ, (2.6) where P is the principal-value operator. Although exact,this expression is nevertheless unusable with experimen-tal data as the algebraic expression of p ( x, θ ) is unknown.However, the projection-slice theorem or Fourier slice theo- p ( x , θ ) p ( k , θ ) W ( q , p ) W ( u , v ) ∼ ∼
1D Fouriertransform 2D FouriertransformRadon transformprojection slice theorem
Figure 1. Different transforms for different paths from p to W . rem [14] gives another reverse path from p ( x, θ ) to W ( q, p ) to work around the difficulties of the principal-value opera-tor (see Fig.1). If ˜ p ( k, θ ) and ˜ W ( u, v ) are, respectively, theone-dimensional and two-dimensional Fourier transforms of p ( x, θ ) and W ( q, p ) , then the projection-slice theorem statesthat ˜ p ( k, θ ) = ˜ W ( k cos θ, k sin θ ) . (2.7)Simply computing the Fourier transform ˜ p ( k, θ ) from themeasured data would seem like the most efficient way to ob-tain W ( q, p ) after a second inverse Fourier transform, butEq. (2.7) shows that it is necessary to interpolate ˜ W ( u, v ) in Fourier space, which leads to significant numerical difficul-ties [15]. To avoid this interpolation Eq. (2.7) can be used toreplace ˜ W ( u, v ) in the inverse Fourier transform of W ( q, p ) to obtain the inversion formula, W ( q, p ) = 12 π Z π Z + ∞−∞ p ( x, θ ) K ( q cos θ + p sin θ − x ) dxdθ. (2.8)Here, the marginal distribution is convoluted with an integra-tion kernel K ( x ) and then back-projected into phase space,where K ( x ) is defined as the inverse Fourier transform of | k | K ( x ) = 12 π Z + ∞−∞ | k | e ikx dk. (2.9)To use Eq. (2.8) in practice it is necessary to regularize K ( x ) and replace it with some numerical approximation. This ispossible with the use of a window function g ( k ) such that theintegral, π Z + ∞−∞ g ( k ) | k | e ikx dk, (2.10)converges. The most common way to regularize Eq. (2.9) is tochoose g ( k ) = [ − k c , + k c ] ( k ) and introduce a hard frequencycutoff parameter k c so that K ( x ) ≈ πx (cos( k c x ) + k c x sin( k c x ) − . (2.11) -3-2-1 0 1 2 3 4 5 -4 -2 0 2 4 K ( x ) x k c = 3 k c = 5 k c = 7 Figure 2. Regularized integration kernel K ( x ) for different values of k c . In practice, the choice of k c affects how much high fre-quency components of the Wigner function will get recon-structed. If k c is set too low the convolution in Eq .(2.8) willfilter out the fine physical details of the Wigner function. If k c is set too high, the convolution will introduce unphysical highfrequency noise from the statistical errors in the measurementof p ( x, θ ) . Figure 2 shows the integration kernel K ( x ) for dif-ferent high frequency sensitivities. Choosing the right value of k c is a trade off between these two regimes. From Eqs. (2.7)and (2.8) it is also possible to insert other filter functions atdifferent steps of the inversion to obtain modified algorithmswith enhanced and more selective noise filtering properties. Inany case the numerical implementation of Eq. (2.8) will relyon deconvolution of the marginal distribution, an operationvery sensitive to statistical noise. III. HARMONIC SERIES EXPANSION
To numerically perform optical homodyne tomography, it isnecessary at some point to apply an approximation procedurefrom the infinite dimensional space which features the un-known physical state to a finite dimensional space used to de-scribe the reconstructed state. In the filtered back-projectionalgorithm, the discretization is achieved by direct evaluationof W ( x i , p i ) on the set of points { ( x i , p i ) } i chosen to probethe phase space. Rather than this point-by-point reconstruc-tion, a discretization of another space should help to solve thenumerical issues encountered in Sec. II. Since we are deal-ing with objects behaving like probability distributions, thestatistical moments of p ( x, θ ) and W ( q, p ) might be a solu-tion to the problem. In Ref.[16], Ourjoumtsev et al. describessuch a technique where they parametrize the Wigner functionof a photon subtracted squeezed vacuum with the second andfourth moments of the marginal distribution p ( x, θ ) . Gener-alizing this approach for any quantum state to higher ordermoments requires the use of the moment generating function (cid:10) e λx (cid:11) , where h x i is the expectation value of x with regards to p ( x, θ ) . Superior to the moment generating function the char-acteristic function (cid:10) e iλx (cid:11) only needs the mean and varianceto be defined to exist. This and the projection-slice theoremof Eq. (2.7) hint that Fourier space is a good candidate for anefficient discretization.We decompose our discretization procedure in two steps: (1) an angular harmonic decomposition with Fourier series;(2) a polynomial series expansion of the radial components.We express W ( q, p ) in radial coordinates ( r, φ ) and notice that W ( r, φ + 2 π ) = W ( r, φ ) . Therefore we write the radial partof W ( r, φ ) in terms of a Fourier series and we define the set ofradial functions, or angular harmonic components { w n ( r ) } n by w n ( r ) = 12 π Z + π − π W ( r, φ ) e − inφ dφ, (3.1)which allows us to write W ( r, φ ) , W ( r, φ ) = ∞ X n = −∞ w n ( r ) e inφ , (3.2)with the symmetry relation w n ( r ) = w ∗− n ( r ) . The 2D Fouriertransform ˜ W ( u, v ) of W ( q, p ) is written in radial coordinates, ˜ W ( k, θ ) = Z + ∞ Z + π − π W ( r, φ ) e − irk cos( θ − φ ) rdrdφ, (3.3)with the change of variables ( u, v ) → ( k, θ ) . ˜ W ( u, v ) is re-lated to the Weyl function χ ( u, v ) = tr (ˆ ρe − iv ˆ q + iu ˆ p ) by asimple π/ rotation, ˜ W ( u, v ) = χ ( − v, u ) , (3.4) ˜ W ( k, θ ) = χ ( k, θ + π . (3.5)We can easily write ˜ W in polar coordinates in terms of theangular harmonic components w n ( r ) of W ( r, φ ) , ˜ W ( k, θ ) = ∞ X n = −∞ Z + ∞ w n ( r ) rdr × Z + π − π e − ikr cos( θ − φ )+ inφ dφ. (3.6)With a Jacobi-Anger expansion of e iz cos φ using Bessel func-tions J n , e iz cos φ = ∞ X n = −∞ i n J n ( z ) e inφ , (3.7)it is possible to conduct the angular integration in Eq. (3.6) toobtain the expression, ˜ W ( k, θ ) = 2 π ∞ X n = −∞ ( − i ) n e inθ Z ∞ w n ( r ) J n ( kr ) rdr. (3.8)Notice that R ∞ w n ( r ) J n ( kr ) rdr is the n th order Hankeltransform of w n ( r ) .In the same fashion, since p ( x, θ + 2 π ) = p ( r, θ ) we de-compose the marginal distribution as p θ ( x ) = ∞ X n = −∞ c n ( x ) e inθ , (3.9)with the sets of radial functions c n ( x ) defined by c n ( x ) = 12 π Z + π − π p ( x, θ ) e − inθ dθ. (3.10)Using the projection-slice theorem of Eq. (2.7) and the or-thogonality of e inθ on [ − π, + π ] we are able to write for everyangular harmonic order n , i n π Z + ∞−∞ c n ( x ) e − ikx dx = Z ∞ w n ( r ) J n ( kr ) rdr. (3.11)We have obtained a relation between, on one side the Fouriertransform of the angular harmonics of p ( x, θ ) , and on theother side, the Hankel transform of the angular harmonics of W ( r, φ ) . If we inverse the Hankel transform with the orthog-onality relation, or closure relation of Bessel functions, Z ∞ kdkJ n ( kr ) J n ( kr ′ ) = 1 r δ ( r − r ′ ) , (3.12)we finally obtain w n ( r ) = i n π Z ∞ J n ( kr ) kdk Z + ∞−∞ c n ( x ) e − ikx dx. (3.13)At that point it would be natural to convey some radial de-composition of w n ( r ) and c n ( x ) . However, there is no simpleway to achieve this. Looking at Eq. (3.13), we notice that theFourier transform of kJ n ( k ) , or at least J n ( k ) , should be in-volved in the process. The latter is written in terms of theChebysheff’s polynomials of the first kind T n Z + ∞−∞ J n ( k ) e − ikx dk = 2( − i ) n √ − x T n ( x ) [ − , +1] ( x ) . (3.14)Equation (3.14) hints at the use of the polynomial series toachieve this radial decomposition. It is safe to assume for ap-plications that the Wigner function will only take nonzero val-ues from the origin up to a certain limit L ≥ r . Since we arecarrying the decomposition in polar coordinates what we arelooking after is a polynomial family which is orthogonal on adisk of radius L . There are of course infinitely many such fam-ilies but one which proves to be particularly adequate to thetask is the set of Zernike polynomials Z ns ( r, ϕ ) = R ns ( r ) e inϕ originally introduced for the study of optical aberrations inlenses and other circular optical systems [17]. The polynomi-als are defined for s ≥ | n | ≥ and s − | n | even. While theangular part gives straightforward orthogonality and fits withour previous approach using Fourier series, the radial compo-nents R ± ns defined for t = | n | ≥ by R ± ns ( r ) = ( s − t ) / X k =0 ( − k ( s − k )! k ! (cid:0) s + t − k (cid:1) ! (cid:0) s − t − k (cid:1) ! r s − k , (3.15)are orthogonal on [0 , with respect to the weight function r for all positive and negative orders n , Z R ns ( r ) R ns ′ ( r ) rdr = 12( s + 1) δ s ′ s . (3.16) Furthermore it turns out that the Radon transform of Zernikepolynomials happens to have the simple expression, R (cid:0) R ns ( r ) e inφ (cid:1) = 2 s + 1 p − x U s ( x ) e inθ , (3.17)where U s ( x ) are the Chebysheff’s polynomials of the sec-ond kind [18, 19] (see also the last paragraph of this sectionfor a proof). The critical aspect for tomography lies in thefact that U s ( x ) is again an orthogonal polynomial family on [ − , with respect to the weight function √ − x . In otherwords by finding a family of orthogonal polynomials whoseRadon transform element by element is yet another family oforthogonal polynomials, we have in some sense diagonalizedthe Radon transform. The inverse Radon transform can alsobe exactly calculated and any technical difficulties associatedwith kernel functions or regularization immediately vanish.With the use of Eq. (3.16) we are eventually able to expandthe angular harmonic functions w n ( r ) on the n th order radialpolynomials R ns ( r ) , w n ( r ) = ∞ X s =0 w sn R ns ( r ) . (3.18)Given that R ns ( r ) is non zero only when s ≥ | n | ≥ and s −| n | is even, we introduce the change of variable s → | n | +2 m ,re-index the sequence w sn and rewrite Eq. (3.18) w n ( r ) = ∞ X m =0 w mn R n | n | +2 m ( r ) . (3.19)Putting Eqs. (3.19) and (3.2) together we obtain the completeexpansion of W ( r, φ ) inside the unit disk D (0 , , W ( r, φ ) = ∞ X n = −∞ ∞ X m =0 w mn R | n || n | +2 m ( r ) e inφ . (3.20)Notice from Eq. (3.15) that R + ns ( r ) = R − ns ( r ) which justi-fies the use of R | n || n | +2 m although w mn are in general complexconstants. Applying the relation (3.17) on Eq. (3.20), p ( x, θ ) is also written in terms of the coefficients w mn as p ( x, θ ) = ∞ X n = −∞ ∞ X m =0 w mn | n | + 2 m + 1 p − x U | n | +2 m ( x ) e inθ . (3.21)To justify the use of Zernike polynomials and prove Eq.(3.17), the relation, Z R nm ( r ) J n ( rk ) rdr = ( − ( m − n ) / J m +1 ( k ) k , (3.22)between Zernike polynomials and Bessel functions [17] is es-sential. If we recall Eq. (3.11), replace w n ( r ) by its expansionon R ns ( r ) in Eq. (3.18) and cut the integration from + ∞ tounity, we obtain ∞ X m =0 w mn ( − m J | n | +2 m +1 ( k ) k = i | n | π Z + ∞−∞ c n ( x ) e − ikx dx. (3.23)To finally obtain the complete inversion of R and the expan-sion of c n ( x ) as in Eq. (3.21), we only need to inverse theFourier transform in Eq. (3.23) from the rhs to the lhs and usethe Fourier transform of J s ( k ) /k , Z + ∞−∞ J s +1 ( k ) k e ikx dk = 2 i s s + 1 U s ( x ) p − x [ − , +1] ( x ) , (3.24)to obtain c n ( x ) = ∞ X m =0 w mn | n | + 2 m + 1 U | n | +2 m ( x ) p − x [ − , +1] ( x ) . (3.25)Notice that Eqs. (3.22) and (3.24) close the link between U s ( x ) and R nm ( r ) , the first two families of orthogonal func-tions used in the analysis, and the Bessel functions J n ( k ) or-thogonal with respect to the weight function /k , Z ∞ J s ( k ) J t ( k ) dkk = 12 s δ ks (3.26)In summary by identifying three families of orthogonal func-tions related together by the Radon transform R and theFourier transform F , we have been able to find an expansionof the Wigner function W ( q, p ) that allows to greatly simplifythe technical difficulties of tomography with inverse Radontransform. IV. RECONSTRUCTION ALGORITHMA. The algorithm
The algorithm works in four steps: (1) choosing the size L of the reconstruction disk, (2) evaluating the coefficients w mn ,(3) choosing the cutoffs N and M of the angular and radialseries, and (4) calculating W ( r, φ ) . Step 1 is necessary forthe orthogonal relations given in Sec. II on [0 , and [ − , +1] to hold. In practice we have to normalize the marginal dis-tribution p ( x, θ ) → p ( x/L, θ ) /L and the Wigner function W ( r, φ ) → W ( r/L, φ ) /L . Step 2 is easily conducted by in-verting the relation (3.21) with the orthogonal Chebysheff’spolynomials U | n | +2 m ( x ) , w mn = | n | + 2 m + 12 π Z + π − π dθe − inθ × Z +1 − dx p ( x/L, θ ) L U | n | +2 m ( x ) . (4.1)The recurrence relation, U s +1 ( x ) = 2 xU s ( x ) − U s − ( x ) , (4.2)allows one to efficiently calculate U s ( x ) for any s and any x given U ( x ) = 1 and U ( x ) = 2 x . After obtaining the coef-ficients w mn and choosing cutoff orders N and M , the Wignerfunction W ( r, φ ) is then approximated by the partial sums, W ′ ( r, φ ) = N X n = − N M X m =0 w mn R | n || n | +2 m (cid:16) rL (cid:17) e inφ /L, (4.3) -4-2024-4 -2 0 2 4-0.2-0.15-0.1-0.0500.050.10.15(h) qp -4-2024-4 -2 0 2 4-0.2-0.15-0.1-0.0500.050.10.15(g) qp -4-2024-4 -2 0 2 4-0.2-0.15-0.1-0.0500.050.10.15(f) qp -4-2024-4 -2 0 2 4-0.2-0.15-0.1-0.0500.050.10.15(b) qp -4-2024-4 -2 0 2 4-0.2-0.15-0.1-0.0500.050.10.15(c) qp -4-2024-4 -2 0 2 4-0.2-0.15-0.1-0.0500.050.10.15(d) qp -4-2024-4 -2 0 2 4-0.2-0.15-0.1-0.0500.050.10.15(e) qp -4-2024-4 -2 0 2 4-0.2-0.15-0.1-0.0500.050.10.15(a) qp Figure 3. Comparison between polynomial series tomography (leftpanels: N = 8 , M = 30 ) and filtered back-projection tomography(right panels: k c = 8 . ) for the state ρ = 0 . | ih | + 0 . | ih | .(a) J = 5 × ; (b) J = 20 × ; (c) J = 80 × ; (d) J =320 × ; (e) J = 5 × ; (f) J = 20 × ; (g) J = 80 × ;(h) J = 320 × . All data sets have been synthetically generatedwith rejection sampling. Using the symmetry relation w m − n = ( w mn ) ∗ , we keep the realpart of Eq. (4.3) and simplify the sum on n to W ′ ( r, φ ) = M X m =0 N X n =0 R nn +2 m (cid:16) rL (cid:17) /L × ( a mn cos( nφ ) + b mn sin( nφ )) . (4.4)where we have defined w mn = ( a mn + ib mn ) / for n ≥ and w m = a m . Figures 3 and 4 show examples of reconstructedWigner functions for a mixture of | i and | i , and a thermalstate respectively. In comparison to filtered back-projectiontomography, polynomial series tomography converges fasterwith fewer numbers of experimental points J . The recon-structed Wigner functions also show less visible artifacts andare overall smoother. To evaluate efficiently R mn ( r ) we noticethat R nn ( r ) = r | n | and then use the recurrence relation [20], R nn +2( m +1) ( r ) = n + 2( m + 1)( m + 1)( n + m + 1) × n(cid:18) ( n + 2 m + 1) r − ( n + m ) n + 2 m − ( m + 1) n + 2( m + 1) (cid:19) R nn +2 m ( r ) − m n + mn + 2 m R nn +2( m − ( r ) (cid:27) . (4.5)In contrast to setting the value of k c , the values of N and M have a real physical meaning. This is a major advantage -4-2024-4 -2 0 2 400.020.040.0 qp -4-2024-4 -2 0 2 4-0.0200.020.040.060.080.10.120.14(e) qp -4-2024-4 -2 0 2 400.020.040.060.080.10.120.14(b) qp -4-2024-4 -2 0 2 4-0.0200.020.040.060.080.10.120.14(f) qp -4-2024-4 -2 0 2 400.020.040.060.080.10.120.14(c) qp -4-2024-4 -2 0 2 4-0.0200.020.040.060.080.10.120.14(g) qp -4-2024-4 -2 0 2 4-0.0200.020.040.060.080.10.120.14(h) qp -4-2024-4 -2 0 2 400.020.040.060.080.10.120.14(d) qp Figure 4. Comparison between polynomial series tomography (leftpanels: N = 8 , M = 30 ) and filtered back-projection tomography(right panels: k c = 8 . ) for a thermal state of mean photon number h ˆ n i = 1 . (a) J = 5 × ; (b) J = 20 × ; (c) J = 80 × ; (d) J = 320 × ; (e) J = 5 × ; (f) J = 20 × ; (g) J = 80 × ;(h) J = 320 × . All data sets have been synthetically generatedwith rejection sampling. -4-2024-4 -2 0 2 4-0.2-0.15-0.1-0.0500.050.10.15( a) qp -4-2024-4 -2 0 2 4-0.2-0.15-0.1-0.0500.050.10.15(b) qp -4-2024-4 -2 0 2 4-0.2-0.15-0.1-0.0500.050.10.15(c) qp -4-2024-4 -2 0 2 4-0.2-0.15-0.1-0.0500.050.10.15(d) qp -4-2024-4 -2 0 2 4-0.2-0.15-0.1-0.0500.050.10.15(e) qp -4-2024-4 -2 0 2 4-0.2-0.15-0.1-0.0500.050.10.15(f) qp Figure 5. Effect of increased radial resolution on the stabilityof tomography of an experimentally measured photon subtractedsqueezed vacuum (same data as in Ref. [21]). For all panels J = 1 × . (a) Polynomial series tomography, N = 8 , M = 20 ;(b) M = 30 ; (c) M = 40 ; (d) filtered back-projection tomography, k c = 7 ; (e) k c = 9 ; (f) k c = 11 . of this method compared to the usual filtered back-projectionalgorithm. M will decide what will be the highest polynomialorder of the radial features of W . Therefore it is equivalent Figure 6. Effect of N and M on the convergence of polynomialseries tomography. Same experimental data as in Fig. 5. (a) Circularcut at constant r and effect of N for M = 32 , J = 2 × . (b)Radial cut at constant φ and effect of M for N = 10 , J = 2 × . to choosing the maximum photon number of the density ma-trix diagonal elements. N will set the resolution of the an-gular features of W , which decides how many off-diagonalcomponents of the density matrix will be reconstructed. Fur-thermore it is easy to change N and M after computing thecoefficients w mn . Figure 5 shows the effect of increasing M on the precision of polynomial series tomography. In compar-ison to filtered back-projection tomography when increasingthe kernel sensitivity k c , increasing the radial resolution M does not produce artifacts in the Wigner function. Figure 6further shows the effect of increasing N and M on the pre-cision of the tomography reconstruction of experimental data.While the angular components show quick convergence, theradial components require higher M values to be faithfullyreconstructed. Figure 7 illustrates the advantage of polyno-mial series tomography in radial resolution for quantum stateswith a higher number of photons. Both M and k c where set atvalues high enough to recover the original Schroedinger’s catstate negativity at the origin of phase space. While the back-filtered projection shows numerical uinstability when k c is sethigh, the Wigner function reconstructed by polynomial seriestomography is smoother at the equivalent resolution.Finally the value of R nn +2 m in r = 0 will be non-zero onlyfor n = 0 , therefore we have the useful formula to evaluatethe Wigner function at the origin of phase space, W ′ (0 ,
0) = M X m =0 ( − m a m /L, (4.6)which is similar to the formulation of W (0 , using the diag-onal elements of the density matrix. B. Unbiased error estimator
To quantitatively compare our algorithm with the usualback-filtered tomography algorithm we give a consistentmethod to estimate the reconstruction error and obtain confi-dence intervals when calculating the value of W ( q, p ) . If W ′ and W ′′ are the reconstructed value of W ( q, p ) with Eqs. (4.4)and (2.8) respectively, we call σ W ′ and σ W ′′ the variance of -4 -2 0 2 4 -4 -2 0 2 4-0. q p -4 -2 0 2 4 -4 -2 0 2 4-0.3-0.2-0.100.10.2(a) q p -4 -2 0 2 4 -4 -2 0 2 4-0.3-0.2-0.100.10.2(b) q p Figure 7. Effect of increased radial resolution on the stability of to-mography of a Schroedinger’s cat states with h ˆ n i = 3 . For all panels J = 4 × . (a) Original Wigner function; (b) polynomial seriestomography, N = 8 , M = 46 ; (c) filtered back-projection tomogra-phy, k c = 11 . the reconstruction errors assuming they are distributed accord-ing to a Gaussian for both algorithms. We also assume thatthere are no systematic errors but only statistical errors. Let’sassume an optical homodyne measurement set consists of J experimental points { ( x j , θ j ) } j independently and identicallydistributed according to the underlying marginal distribution p ( x, θ ) . To begin with we give an estimator of σ W ′′ ( q, p ) for the usual filtered back-projection method using formula(2.8). To calculate the value of W at point ( q, p ) , p ( x, θ ) willbe replaced either by a binned histogram made from the dataset { ( x j , θ j ) } j , or by a sum of delta functions approximating p ( x, θ ) p ( x, θ ) = 1 J X j δ ( x − x j ) × δ ( θ − θ j ) . (4.7)In the latter case, the swap of p ( x, θ ) for expression (4.7) inEq. (2.8) leads to W ′′ ( q, p ) = 12 πJ J X j =1 K ( q cos θ i + p sin θ i − x i ) . (4.8) Since p ( x, θ ) is a valid probability distribution W ′′ ( q, p ) isnothing else than h K ( q cos θ + p sin θ − x ) i the expectationvalue of the kernel function. Therefore Eq. (4.8) can be re-garded as a Monte Carlo integral where the expectation valueof the kernel function is calculated by randomly sampling K according to the distribution p ( x, θ ) . In other words, the op-tical homodyne tomography with filtered back-projection isin effect an analogical Monte Carlo integration where the ho-modyne measurement plays the part of the random numbergenerator. In that familiar case the statistical properties of thereconstruction error are well known. First of all we are as-sured of the unbiased convergence of the sum in Eq. (4.8).The central limit theorem also states that the error will indeedconverge to a Gaussian distribution of zero mean and whosestandard deviation σ W ′′ ( q, p ) for J experimental points is σ W ′′ ( q, p ) = σ K / √ J − , (4.9)which exhibits a / √ J rate of convergence, and where σ K = p h K i − h K i / π . By using the approximations, h K i ≈ J J X j =0 K ( q cos θ j + p sin θ j − x j ) , (4.10) h K i ≈ J J X j =0 K ( q cos θ j + p sin θ j − x j ) , (4.11)we can actually estimate σ K in a straightforward way easy toinclude in the implementation of Eq. (4.8).The same analysis for the coefficients { w mn } yields the re-construction sum, w mn = | n | + 2 m + 12 π J X j =1 U | n | +2 m ( x j /L ) e − inθ j /L. (4.12)As previously errors are Gaussian distributed for every coeffi-cient w mn with a / √ J rate of convergence. If a quantity Y iscalculated through the measure of the variables { y i } i ≤ I withthe formula, Y = f ( y , . . . , y I ) , (4.13)then the variance σ Y of Y can be approximated by σ Y = I X i =1 ( ∂ y i f ) σ y i + 2 X j>i ( ∂ y i f ) (cid:0) ∂ y j f (cid:1) σ y i y j , (4.14)where σ xy = h xy i − h x ih y i . Using Eq. (4.4) we can ap-ply this formula to estimate the variance σ W ′ anywhere inphase space, but because of its simple formulation thanks toEq. (4.6), we will only study it at the origin (0 , : σ W ′ (0 ,
0) = 1( J − L M X m =0 σ a m + 2 M X k>m ( − m + k σ a m a k ! . (4.15) Figure 8. Estimation of σ W (0 , with filtered back-projection to-mography (plain line) and polynomial series tomography (dottedlines). (a) ρ = 0 . | ih | + 0 . | ih | ; (b) thermal state with h ˆ n i = 1 ;(c) photon subtracted squeezed vacuum (same data as in Fig.5). Notice that in this case the variance estimator formula of Eq.(4.14) is not an approximation anymore due to the linear com-bination nature of Eqs. (4.4) or (4.6). We can compute an esti-mate of σ a m when computing the coefficients w mn in the sameway we did with Eqs. (4.10) and (4.11). Figure 8 shows es-timation of the reconstruction errors for different states usingEq. (4.9) and (4.15). We have found that the value of k c hasvery little influence on σ W ′′ at the center of phase space. On the contrary M has a strong influence on σ W ′ (0 , . However,as was shown in Figs.3, 4 and 5, far from the origin the poly-nomial series tomography algorithm shows less uncertainties. Figure 9. Effect of M on the convergence of W ′ (0 , and the mag-nitude of σ W ′ (0 , . (a) Thermal state with h ˆ n i = 1 , rejection sam-pling. (b) Experimental photon subtracted squeezed vacuum state(same data as in Fig.5). We also assumed the convergence error due to finite trun-cation N and M of the expansion to be smaller than the sta-tistical error itself. This can be checked in the algorithm byiteratively calculating σ W ′ (0 , for increasing values of M and stop when the magnitude of the M th and last coefficient w M is less than σ W ′ (0 , (see Fig. 9). This technique can berepeated independently for every point of phase space ( q, p ) ,and different values of N and M can even be used for differentpoints of phase space. C. Monte Carlo error estimation
Independently from the estimators of the previous para-graph, we also use Monte Carlo simulations to generate manysynthetic data sets and evaluate the reconstruction errors. Thismethod is easily applied if we know precisely which state | ψ i is under investigation. For example, we can choose a knowndensity matrix or Wigner function and calculate the associatedmarginal distribution p ( x, θ ) . From this marginal distributionwe generate K synthetic data sets of J points { ( x j , θ j ) } ( k ) j using, for example, rejection sampling. With the algorithmof our choice we repeat the tomography reconstruction andcalculate a set of K Wigner function { W ( k ) } k . Finally for a Figure 10. Comparison between Monte Carlo simulation and directestimation of σ W (0 , . Black curves are the estimation of σ W (0 , with Monte-Carlo simulation using K data sets. Dashed curves arethe direct estimation of σ W (0 , using Eqs. (4.9) and (4.15) forthe K th data set. (a) Data sets of J = 10 points generated usingrejection sampling for the state . | ih | + 0 . | ih | . (b) Data setsof J = 10 points generated with bootstrapping resampling fromthe same experimental data as Fig. 5. (i) Filtered back-projectiontomography with k c = 7 ; (ii) polynomial series tomography with M = 10 ; (iii) M = 20 ; (iv) M = 30 ; (v) M = 40 . given point of phase space ( x , p ) , we calculate the averagevalue ¯ W of the set { W ( k ) } k : ¯ W = 1 K K X k =1 W ( k ) ( x , p ) , (4.16)and obtain an estimate of the error σ ¯ W at point ( x , p ) by σ W = 1 K K X k =1 (cid:16) W ( k ) ( x , p ) − ¯ W (cid:17) . (4.17)Since it is a Monte Carlo based simulation, every quantityshows again a / √ K convergence rate.With experimental data, we can sample p ( x, θ ) only onceand therefore we need a technique to generate the syntheticdata sets after the experimental measurement. Resampling isthe easiest approach and here we estimate the reconstructionerror of experimental data sets with the bootstrapping resam-pling method [22]. The results of both techniques are illus-trated in Fig. 10 and overall there is a good agreement be-tween the estimated values of Monte Carlo simulations andthe predicted value of σ W (0 , using Eq. (4.9) or (4.15). D. Distance to a target state
To conclude this comparative study of polynomial seriesexpansion and filtered back-projection-based tomography, wenumerically estimate in this final paragraph the distance be-tween some original target quantum state and reconstructedstates using both algorithms. For this purpose we will considerone distance for the Wigner function and one distance for thedensity matrix. We use the L Euclidian distance d L ( ., . ) for Figure 11. Estimation of the distance between the target thermal stateof mean photon number h ˆ n i = 1 and reconstructed quantum statesaveraged over 1000 samples of J data points for different tomogra-phy settings. (a) L distance h d L ( W target , W tomo ) i . (b) Frobeniusdistance h d F (ˆ ρ target , ˆ ρ tomo ) i .Figure 12. Estimation of the distance between the target state . | ih | + 0 . | ih | and reconstructed quantum states averagedover 1000 samples of J data points for different tomography set-tings. (a) L distance h d L ( W target , W tomo ) i . (b) Frobenius distance h d F (ˆ ρ target , ˆ ρ tomo ) i .Figure 13. Estimation of the distance between the target oddSchroedinger’s cat state ∝ | α i − | − α i with h ˆ n i = 3 andreconstructed quantum states averaged over 1000 samples of J data points for different tomography settings. (a) L distance h d L ( W target , W tomo ) i . (b) Frobenius distance h d F (ˆ ρ target , ˆ ρ tomo ) i . the Wigner function defined by d L ( W A , W B ) = (cid:18)Z Z dxdp | W A ( x, p ) − W B ( x, p ) | (cid:19) / , (4.18)0and with the Frobenius norm k . k F defined by k A k F = p tr ( A ∗ A ) = X i,j | A ij | / , (4.19)we define a distance d F ( ., . ) for density matrix as d F (ˆ ρ A , ˆ ρ B ) = k ˆ ρ A − ˆ ρ B k . (4.20)First we choose a target state and derive its exact Wigner func-tion W target and density matrix ˆ ρ target . We then evaluate the dis-tances from the target state according to Eqs. (4.18) and (4.20)using as before Monte Carlo sampling techniques. Rather thanaveraging a reconstructed state over many simulated data sets,we average the distance computed over many reconstructedstates and estimate the numbers: h d L ( W target , W tomo ) i and h d F (ˆ ρ target , ˆ ρ tomo ) i . (4.21)Numerical simulation results are shown in Figs. 11-13 for,respectively, a thermal state with h ˆ n i = 1 , a mixture of vac-uum and one-photon state . | ih | + 0 . | ih | , and an oddSchroedinger’s cat state with h ˆ n i = 3 . In agreement withthe previous results on tomography uncertainties, we observethat polynomial series expansion tomography performs bet-ter than filtered back-projection for these two first cases. Inthe case of the Schroedinger’s cat state ∝ | α i − | − α i , bothdistances behave differently for higher J and tend to reach aprecision limit which depends on the tomography algorithmand settings. Although the exact cause of this saturation isunknown, we believe it is due to the significantly more com-plex structure of the Schroedinger’s cat state. According toour simulations, it seems to depend only on the radial and an-gular precision settings, more precisely on parameters M , N ,and k c . In this case again, polynomial series expansion provesto reach a higher precision level than filtered back-projectionfor a relevant range of tomography settings. To conclude thisparagraph, it is interesting to notice that in the case of the d L ( ., . ) distance there is an intrinsic limitation on the preci-sion of polynomial series expansion tomography due to thecircular geometry of the reconstruction space [23]. This couldbe the reason for the saturation phenomenon visible in Fig.13. V. CONCLUSION
We have shown and demonstrated a technique for opticalhomodyne tomography based on polynomial series expansion of the Wigner function. In Sec.II we have given the basisof the usual filtered back-projection algorithm and explainedthe main reason for its weak performances against statisti-cal noise. We have also introduced the projection-slice the-orem and the relation between phase space, Fourier space andthe marginal distribution. In Sec.III we have shown that itis possible to link three families of orthogonal functions be-tween these three spaces to decompose p ( x, θ ) the marginaldistribution, W ( q, p ) the Wigner function, and their Fouriertransforms. We have shown that the Radon transform pre-serves the orthogonality of these families and therefore takesan especially simple form in this case. In Sec.IV we haveexplained and applied to experimental and simulated data themost straightforward implementation of that technique with adirect linear estimation of the coefficients of the polynomialseries expansion. We have also provided estimators of the re-construction errors and shown that it performs better than fil-tered back-projection tomography with respect to reconstruc-tion artifacts and statistical errors. More precisely, polynomialseries tomography is superior with fewer experimental datapoints and when higher radial resolution is needed for higherphoton number states. These results are confirmed when look-ing at the distance between a chosen target state and states re-constructed with both tomography techniques. Furthermorethis technique exploits the projection slice theorem directlyand therefore is faster than convolution based filtered back-projection. Finally we remark that it is in principle possible touse the maximum likelihood technique to find the set of coef-ficients w nm that maximizes the probability of measuring theexperimentally measured data set. ACKNOWLEDGMENTS
This work was partly supported by the Strategic Informa-tion and Communications R& D Promotion (SCOPE) pro-gram of the Ministry of Internal Affairs and Communicationsof Japan, Project for Developing Innovation Systems, Grants-in-Aid for Scientific Research, Global Center of Excellence,Advanced Photon Science Alliance, and Funding Program forWorld-Leading Innovative R&D on Science and Technology(FIRST) commissioned by the Ministry of Education, Cul-ture, Sports, Science and Technology of Japan, and ASCR-JSPS, the Academy of Sciences of the Czech Republic andthe Japanese Society for the Promotion of Science. [1] D. T. Smithey, M. Beck, M. G. Raymer, A. Faridani, Phys. Rev.Lett. , 1244 (1993).[2] K. Vogel, H. Risken, Phys. Rev. A , 2847 (1989).[3] J. Radon, Berichte der Sachsischen Akadamie der Wissenschaft , 262 (1917)[J.Radon (by P. C. Parks), IEEE Transactions on medical imaging MI-5, 170(1986)].[4] G. M. D’Ariano, U. Leonhardt, H. Paul, Phys. Rev. A , 1801(1995).[5] U. Leonhardt, M. G. Raymer, Phys. Rev. Lett. , 1985 (1996).[6] G. Drobny, V. Buzek, Phys. Rev. A , 053410 (2002). [7] Z. Hradil, Phys. Rev. A , R1561 (1997).[8] A. I. Lvovsky, J. Opt. B: Quantum Semiclass. Opt. , S55(2004).[9] H. Moya-Cessa, P. L. Knight, Phys. Rev. A , 2479 (1993).[10] S. Del´eglise, I.Dotsenko, C. Sayrin, J. Bernu, M. Brune, J.-M.Raymond, S. Haroche, Nature , 510-514 (2008).[11] W. G. Hawkins, H. H. Barrett, SIAM J. Numer. Anal. , 873(1986).[12] N. C. Rouze, V. C. Soon, G. D. Hutchins, Pattern RecognitionLett. , 636-642 (2006).[13] U. Leonhardt, Measuring the Quantum State of Light (Cam-bridge University Press, Cambridge, 1997).[14] R. N. Bracewell, Aust. J. Phys. , 198 (1956).[15] H. Stark, J. W. Woods, I. Paul, R. Hingorani, IEEE Trans, Biomed. Eng. , 496-505 (1981).[16] A. Ourjoumtsev, R. Tualle-Brouri, J. Laurat, and P. Grangier,Science , 83 (2006).[17] M. Born, E. Wolf, Principles of Optics (Cambridge UniversityPress, Cambridge, 1999).[18] S. R. Deans,
The Radon Transform and Some of Its Applications (Dover Publications, New York, 2007).[19] E. Zeitler, Optik , 396-415 (1974).[20] A. Prata, W. V. T. Rusch, Appl. Opt. , 749-754 (1989).[21] N. Lee, H. Benichi, Y. Takeno, S. Takeda, J. Webb, E. Hunting-ton, and A. Furusawa, Science , 330 (2011).[22] B. Efron, R. J. Tibshirani, An Introduction to the Bootstrap (Chapman & Hall/CRC, New York, 1994).[23] M. Pawlak, S. X. Liao, IEEE Trans. Info. Theory48