Quantum polarization tomography of bright squeezed light

We reconstruct the polarization sector of a bright polarization squeezed beam starting from a complete set of Stokes measurements. Given the symmetry that underlies the polarization structure of quantum fields, we use the unique SU(2) Wigner distribution to represent states. In the limit of localized and bright states, the Wigner function can be approximated by an inverse three-dimensional Radon transform. We compare this direct reconstruction with the results of a maximum likelihood estimation, finding an excellent agreement.


Introduction
Polarization of light is a robust characteristic that can be efficiently manipulated using modest equipment without introducing more than marginal losses. It is thus not surprising that this is often the variable of choice to encode quantum information, as one can convince oneself by looking at some recent cutting-edge experiments, including quantum key distribution [1], quantum dense coding [2], quantum teleportation [3], rotationally invariant states [4], phase super-resolution [5], and weak measurements [6].
In the discrete-variable regime of single, or few, photons, one is mostly interested into two-mode states, which for all practical purposes can be regarded as a spin system [7,8]. As a result, the polarization state can be determined from correlation functions of different orders [9][10][11][12][13][14][15][16][17][18][19]. Given the small dimensionality of the Hilbert space involved, the state reconstruction can be readily performed.
In the continuous-variable case, polarization properties are exploited for an expedient generation, manipulation, and measurement of nonclassical light. Polarization squeezing [20][21][22][23], which has been observed in numerous experiments [24][25][26][27][28], is perhaps the most tantalizing illustration. Full Stokes polarimetry [29] is the method employed by the majority of the practitioners in this area.
However, the reconstruction in this limit is a touchy business and requires special care. The origin of the problem can be traced back to the fact that the characterization of the polarization state by the whole density operator is superfluous, because it contains much more than polarization information. This redundancy can easily be handled for low number of photons, but becomes a real hurdle for highly excited states. An adequate solution has been proposed recently: it suffices with a subset of the density matrix that has been called the "polarization sector" [30,31] or the "polarization density operator" [32]. Its knowledge allows for a complete specification of the state on the Poincaré sphere (actually on a set of nested spheres that can be appropriately called the Poincaré space). The technique was devised by Karassiov and coworkers [33][34][35] and implemented experimentally in [36].
In this paper we present a comprehensive treatment of polarization tomography. As with any reliable quantum tomographical scheme, we need to supply three key ingredients [37]: the availability of a tomographically complete measurement, a suitable representation of the quantum states, and a robust algorithm for inverting the experimental data. In this respect, we use a standard Stokes scheme that implements the first ingredient in a very simple way; for the second, we resort to the well-known SU(2) Wigner distribution [38][39][40][41][42][43][44][45], and finally, we prove that the inversion of the data in terms of that Wigner function is an inverse three-dimensional (3D) Radon transform.
To support the experimental feasibility of our scheme, we carry out the full tomography of a bright polarization-squeezed state generated in a Kerr medium [46]. The reconstruction is accomplished in three-different ways: by the direct inversion of the Radon transform, by a maximum-likelihood estimation and, finally, by a Gaussian approximation. The final results are compared and the sources of uncertainty are analyzed.

Polarization structure of quantum fields
A satisfactory description of the polarization structure of quantum fields and the corresponding observables that specify this structure is of paramount importance for our purposes.
We restrict our attention to the case of a monochromatic plane wave (the formalism can be easily extended to more involved multimode wavefronts [47,48]), which we assume to propagate in the z direction, so its electric field lies in the xy plane. Under these conditions, we are dealing with a two-mode field that can be fully characterized by two complex amplitude operators. They are denoted byâ H andâ V , where the subscripts H and V indicate horizontally and vertically polarized modes, respectively. The commutation relations of these operators are The description is greatly simplified if we use the Schwinger representation [49,50] together with the total number These operators coincide, up to a factor 1/2, with the Stokes operators [51], whose average values are precisely the classical Stokes parameters [52]. Using equation (2.1), one immediately notices thatĴ = (Ĵ 1 ,Ĵ 2 ,Ĵ 3 ) satisfy the commutation relations of the su(2) algebra where ε k m is the Levi-Civita fully antisymmetric tensor. This noncommutability precludes the simultaneous exact measurement of the physical quantities they represent. Among other consequences, this implies that no field state (apart from the vacuum) can have sharp nonfluctuating values of all the operatorsĴ simultaneously. This is expressed by the uncertainty relation where the variances are given by In other words, the electric vector of a monochromatic quantum field never traces out a definite ellipse.
In classical optics, the total intensity is a well-defined quantity and the Poincaré sphere appears thus as a smooth surface with radius equal to that intensity. In quantum optics we haveĴ and, as fluctuations in the number of photons are unavoidable (leaving aside photon-number states), we are forced to talk of a three-dimensional Poincaré space (with axis J 1 , J 2 and J 3 ) that can be envisioned as foliated in a set of nested spheres with radii proportional to the different photon numbers that contribute significantly to the state. The Hilbert space H of these two-mode fields has a convenient orthonormal basis in the form of Fock states for both polarization modes, namely |n H , n V . However, since each subspace with a fixed number of photons N must be handled separately. In other words, in the previous onion-like picture of the Poincaré space, each shell has to be addressed independently. This can be emphasized if instead of the Fock basis, we employ the relabeling According to (2.6), we have that J = N/2 and this basis can be also seen as the common eigenstates of {Ĵ 2 ,Ĵ 3 }. In this way, for each fixed J (i.e., fixed number of photons N), m runs from −J to J and these states span a (2J + 1)-dimensional subspace whereinĴ acts in the usual way (in unitsh = 1) It is clear from all this previous discussion that the moments of any energy-preserving observable (such asĴ) do not depend on the coherences between different subspaces. The only accessible information from any state described by the density matrixρ is thus its polarization sector, which is specified by the block-diagonal form whereρ (J) is the reduced density matrix in the J subspace. Anyρ and its associated blockdiagonal formρ pol cannot be distinguished in polarization measurements (and, accordingly, we drop henceforth the subscript pol). This is consistent with the fact that polarization and intensity are, in principle, independent concepts: in classical optics the form of the ellipse described by the electric field (polarization) does not depend on its size (intensity).

Polarization squeezing and the dark plane
The variances of the angular-momentum operators (2.2) are not independent, for they are constrained by (2.11) It is always possible to find pairs of maximally conjugate operators for this uncertainty relation. This is equivalent to establishing a basis in which only one of the operators (2.2) has a nonzero expectation value, say Ĵ k = Ĵ = 0 and Ĵ m = 0. The only nontrivial Heisenberg inequality reads thus (2.12) Polarization squeezing can then be sensibly defined as [20][21][22][23]: 13) The choice of the conjugate operators {Ĵ k ,Ĵ } is by not means unique: there exists an infinite set {Ĵ ⊥ (θ ),Ĵ ⊥ (θ + π/2)} that are perpendicular to the state classical excitationĴ m , for which Ĵ ⊥ (θ ) = 0 for all θ . All these pairs exist in the J k -J plane, which is called the "dark plane" because it is the plane of zero mean intensity. We can express a genericĴ ⊥ (θ ) aŝ J ⊥ (θ ) =Ĵ k cos θ +Ĵ sin θ , θ being an angle defined relative toĴ k . Condition (2.13) is then equivalent to whereĴ ⊥ (θ sq ) is the squeezed parameter andĴ ⊥ (θ sq + π/2) the antisqueezed parameter.
In the experiments presented in this paper, a focal role will be played by the example in which the horizontal and vertical modes have the same amplitude but are phase shifted by π/2: 2, α being a real number. This light is circularly polarized and fulfills Ĵ 1 = Ĵ 3 = 0, Ĵ 2 = α 2 . It is advantageous to work in the circular polarization basis, whose right (+) and left (−) amplitudes are given in terms of the linear ones bŷ In this manner, â + = α and â − = 0. The operators in the J 1 -J 3 plane correspond to the quadrature operators of the dark left-polarized mode. In fact, expressing the fluctuations ofĴ in terms of the noise of the circularly polarized modes δâ ± and assuming | δâ ± | α we find [53] √ 2 are the rotated quadratures for the ith amplitude. On the other hand, we have that 17) and the intensity exhibits no dependence on the dark mode. In consequence, the condition (2.14) can be recast for this example as 18) that is, polarization squeezing is equivalent to vacuum squeezing in the orthogonal polarization mode.
In the dark-plane measurements, the beam is divided equally between two photodetectors. Such measurements are then identical to balanced homodyne detection: the classical excitation is a local oscillator for the orthogonally polarized dark mode. The phase between these modes is varied by rotating the measurement through the dark plane, allowing a full characterization of the noise properties. This is a unique feature of polarization measurements and has been used in many experiments [54][55][56][57][58].

Polarization Wigner function
The structure discussed so far highlights that SU(2) is the symmetry group for polarization. To provide an appropriate phase-space description of states, we take advantage of the pioneering work of Stratonovich [38] and Berezin [39], who introduced quasi-probability distribution functions on the sphere satisfying all the proper requirements. This construction was later generalized by others [40][41][42][43][44][45] and has proved to be very useful in visualizing properties of spinlike systems [59][60][61][62][63] To gain physical insights into this approach, let us start by representing the density matrix with respect to the polarization basis. Instead of using directly the states {|J, m }, it is more convenient to write such an expansion aŝ where the irreducible tensor operators T

(J)
Kq are [64] T (J) with C Jm Jm,Kq being the Clebsch-Gordan coefficients that couple a spin J and a spin K (0 ≤ K ≤ 2J) to a total spin J. These tensor operators have the right transformation properties under rotations and they indeed constitute the most suitable orthonormal basis Although at first sight they might look a little bit intimidating, they are nothing but the multipoles used in atomic physics [65]: one can check that and similarly theT Kq can be related to the Kth power of the generators (2.2). Accordingly, the expansion coefficients are known as state multipoles. The Wigner function associated with the state (2.19) is and Y Kq (θ , φ ) are the spherical harmonics. This kernel is unitary and satisfies the normalization conditions The integral extends over the unit sphere S 2 and dΩ is the invariant measure therein, namely, dΩ = sin θ dθ dφ .
From (2.24) and the properties of the irreducible tensors, one can immediately express the Wigner function in the very suggestive form which clearly shows that determining this Wigner function is tantamount to the knowledge of all the state multipoles. (i.e., all the moments of the Stokes parameters).
One can obtain the marginal of W (J) (θ , φ ) once summed over all the values of J where the factor has been introduced to ensure the proper normalization.
In the above-mentioned example of a strong circularly polarized state, we can consider that the sphere can locally be replaced by its tangent plane since J α 2 . Using simple geometrical relations between the coordinates (θ , φ ) and the Cartesian coordinates (q, p) in that tangent plane, we get 29) which confirms that the dark plane is equivalent to the standard phase space for continuous variables.
In the limit of large photon numbers the representation (2.24) is not very useful. In such a case, a remarkably effective approximation is given by [66] where n = (cos θ sin φ , sin θ sin φ , cos θ ) is the unitary vector in the direction (θ , φ ).

Tomograms and tomographic inversion
A general polarimetric apparatus consists of a half-wave plate, with axis at angle α, followed by a quarter-wave plate at angle β . For fixed values of the angles (α, β ) of the wave plates, the selected direction in the Stokes space is The polarization transformations performed by the wave plates can be represented byĴ 2 , which generates rotations about the direction of propagation, andĴ 3 , which generates phase shifts between the modes. Their joint action is given by the operator which describes displacements over the sphere. After that, a polarizing beam splitter projects onto the basis |J, m . In physical terms, the wave plates transform the input polarization allowing the measurement of different Stokes parameters by the projection onto the basis |J, m . This can be modeled byΠ (2.34) The experimental histograms recorded for each direction (θ , φ ) correspond to the tomographic probabilities (2.35) The reconstruction in each (2J + 1)-dimensional invariant subspace can now be carried out exactly since it is essentially equivalent to a spin J [67][68][69][70]. One can proceed in a variety of ways, but perhaps the simplest one is to look for an integral representation of the tomograms (2.35); as soon as one realizes that the tomograms read as that is, they appear as the Fourier transform of the characteristic function for the observablê J · n. After some direct manipulations, we find that where dn indicates integration over the unit sphere and the kernel K (x) is where χ J (ω ) is the character of the (2J + 1)-dimensional representation of the SU(2) group [64] and ω is given by cos ω 2 = n · n sin ω 2 . (2.41) For J 1, m can be taken as a continuous variable. Replacing the sum by an integral, integrating by parts and taking into account that for localized states n · n 1, the Wigner function simplifies to where we have included J as an argument to stress that it must be treated as continuous. The reconstruction in this limit turns out to be equivalent to an inverse Radon transform of the measured tomograms.

The setup
To validate our approach, we perform the tomography of a polarization squeezed state, generated in a polarization-maintaining optical fibre through the nonlinear Kerr effect [46]. The setup is shown in figure 1. The light source is a shot-noise limited ORIGAMI laser from Onefive GmbH emitting 220 fs pulses at a repetition rate of 80 MHz and centered at 1560 nm. The light is fed into a 13 m-long polarization-maintaining birefringent fibre (3M FS-PM-7811, 5.6 µm mode-field diameter) so that quadrature squeezed states are simultaneously and independently generated in both polarization modes. The strong birefringence of the fibre (beat length 1.67 mm) causes a delay between the emerging quadrature-squeezed pulses. We precompensate for this delay in an unbalanced Michelson-like interferometer placed before the fibre. A small part (0.1 %) of the fibre output serves as the input to a control loop to maintain the relative phase between the exiting pulses locked to π/2, so the light is circularly polarized.
The quantum state is detected with a Stokes measurement, as sketched in the previous section. The two detectors (InGaAS PIN photodiodes, custom-made by Laser Components GmbH with 98 % quantum efficiency at DC) are balanced and have a sub-shot noise resolution at a frequency range between 5 MHz and 30 MHz. Each detector has two separate outputs: DC, providing the average values of the photocurrents, and AC, providing the photocurrents amplified in radio-frequency (RF) spectral range. The RF currents of the photodetectors are mixed with an electronic local oscillator at 12 MHz, amplified (FEMTO DHPVA-100), and digitized by an analog/digital exit converter (Gage CompuScope 1610) at 10 Megasamples per second with a 16-bit resolution and 10 times oversampling.
The measurements are performed at a pulse energy of 93 pJ. In the dark plane a total squeezing of about 3.8 dB is observed. In the orthogonal quadrature, the noise was enhanced by several tens of dB due to Guided Acoustic Wave Brillouin Scattering (GAWBS) [71][72][73]. In the direction of the classical excitation, the state is expected to be shot-noise limited, since the Kerr effect only influences the phase and does not contribute to the photon number.
To perform the reconstruction, histograms of the Stokes variables are recorded for different angles (θ , φ ). This is done by rotating the wave plates with motorized stages (OWIS DMT 40-D20-HSM) and scanning one eighth of the Poincaré sphere in 8100 steps, a measurement that took over eight hours. The unmeasured octants were deduced from symmetry. For each setting of the wave plates, the photocurrent noise of both detectors was  simultaneously sampled 0.5 × 10 6 times. Noise statistics of the detectors difference current were acquired in histograms with 751 bins. Additionally, the optical intensity was recorded.
In figure 2 we show typical histograms at different angles on the Poincaré sphere. As the widths largely vary from squeezing to antisqueezing ranges, there are two plots in which the amplitude scale differs in more than one order of magnitude. The histograms labeled 1, 2 and 3 are measured in the dark plane. Tomogram 1 denotes the angle of maximum squeezing, while 3 corresponds to the antisqueezing. Tomogram 4 is at the classical mean value, where the measured noise is almost shot-noise limited. Due to the high number of samples, the measured histograms are smooth and, at the same time, the number of bins makes it possible to resolve the large dynamical range of amplitudes, so no data interpolation was needed. We also plot histograms showing the electronic noise and the shot noise.
For all these histograms we have performed a Gaussianity check, using the Kolmogorov-Smirnov and the χ 2 tests, as well as the Kullback-Leibler divergence [74]. We can conclude that all the histogramas are Gaussian within confidence levels ranging from 95 % to 98 %.

Experimental reconstruction
As clearly expressed in (2.42), for high photon numbers the tomography turns out to be equivalent to an inverse Radon transform of the measured histograms. In practice, this onestep 3D reconstruction is very demanding in computational resources. Therefore, we divide the process into two steps: first, a set of 2D projections is reconstructed from the recorded histograms; subsequently, the Wigner function is slice-wise generated from those projections (to which we apply a Hamming filter to smooth the noise). The symmetry of the state is used as a prior information to reduce the range of measured angles to an octant. This minimizes the systematic errors stemming from imperfections of the polarization optics.
In figure 3 (right panel) we show the final result of the inverse Radon transform for our polarization squeezed state. More concretely, we plot an isocontour surface of W (J, θ , φ ) = constant (with the constant being 1/e from the maximum) in the Poincaré space having J 1 , J 2 , and J 3 as orthogonal axes. As coordinate units we use the shot noise set by a coherent state. The ellipsoidal shape of the state is clearly visible. The center of the ellipsoid is far away from the origin, since we have 10 × 10 11 photons per measurement time (using 1.9 MHz resolution bandwidth). The antisqueezed direction of the ellipsoid is dominated by excess noise stemming largely from GAWBS, as we have already mentioned.
In the left panel of figure 3 we sketch density plots of the projections on the coordinate planes of the previous Wigner function (including the particular case of a coherent state). The contours agree with the 3.8 ± 0.3 dB squeezing that was directly measured from the variances. The projections on the planes J 1 -J 2 and J 2 -J 3 show an additional spreading of the state in the J 2 direction caused by the imperfect polarization contrast in the measurement setup that mixes some of the antisqueezing on the J 2 direction.
This Radon reconstruction requires a large set of measured data to get a reasonably accurate representation of the state. There are two main reasons for this: integrals are approximated by finite sums (in our case, we used 751 bins in 91 steps) and the kernel (2.39) is singular, so some ad hoc filtering of the raw data is needed. Acquiring such large data sets may be unwise, for it demands long measurement times. Ensuring the proper stability of the setup is thus essential and might be difficult depending on the quantum state measured.
This limitation may be circumvented by adopting a statistically-motivated method, such as the maximum likelihood (ML) [75]. In our case, the relation between the Wigner function W and the tomograms w can be written as a system of linear equations where the subscripts in W k and w j is a shorthand notation for the respectives coordinates. The coefficients c jk can be interpreted as the overlap of the jth projector with the kth volume element of the Wigner function and can be readily determined from equations (2.27) and (2.35). The most likely Wigner function is then found by minimizing the Kullback-Leibler divergence between the normalized vectors of the computed tomograms w j and recorded onesw j . Technically, this can be achieved by the iterative expectation-maximization algorithm [76][77][78] which converges monotonously to the ML estimate from any strictly positive initial vector W The significantly greater stability of the statistical inversion allows us to get reconstructions of the same quality but from far smaller data sets. This is illustrated in figure 4, where we draw a comparison between Radon and ML methods, although in the latter case using only nine different settings of the angles (θ , φ ), which amounts to reducing the measurements by two orders of magnitude. In other words, the measuring time is shortened from eight hours to less than five minutes! This result indicates that the experimental characterization of considerably more complicated quantum states with less symmetries should still be within the reach of the present measurement setup.
As we have discussed in section 2.2, the dark plane is of special interest. The theory shows that the reconstruction therein can be obtained in two different ways: either by reconstructing the dark mode directly from the histograms or by calculating projection of the 3D Wigner function along the J 2 direction. The two results are compared in figure 5 and good agreement within the experimental uncertainties is found. Since the Radon transform should provide a plausible explanation for all the measured histograms, such a comparison may serve as an independent test of the quality of the 3D tomography.
Finally, the high confidence levels of the Gaussianity tests seems to call for a Gaussian ML reconstruction. The Gaussianity is used as a prior information about the signal, which helps to reduce drastically the number of free parameters. In this case, the Wigner function is represented by the 3 × 3 covariance matrix G: and the calculated variances σ j = n j Gn j are matched (in the ML sense) to the actually measured variances [79]. Since G must be positive semidefinite, only six real parameters describe the measured system and the problem is highly overdetermined, in consequence, the Gaussian state can be obtained from a few histograms. In principle, by comparing Gaussian reconstructions based on different subsets of measured data, various imperfections of the setup, such as instabilities and biases, can be detected.
The matrix G turns out to be which once diagonalized gives the principal variances 0.43962, 1.13853, and 309.22177 (in shot-noise units). This agrees well with the standard and ML reconstructions, as can be also appreciated in figure 4. The Gaussian reconstruction was done without assuming a particular orientation or symmetry of the state with respect to the Stokes coordinates. The covariance matrix suggests that the misalignment of the principal axis is less than 0.5 degrees within the measurement errors, in accordance with the definition of angles adopted in the experiment. This Gaussian approach allows for a simple estimate of the errors: just take the pseudoinversion of the measurement matrix as a linear model and use the standard theory of error propagation. The errors to be propagated are actually the errors in the estimated variances for each tomogram, which are found from the χ 2 distribution. In addition, we can assume that the variances of different tomograms are uncorrelated.
Taking a 97.5 % confidence interval (which corresponds to three standard deviations), the principal variances can be written as 0.440 ± 0.002, 1.139 ± 0.001, 309.2 ± 0.3 . (3.5) Note that the relative errors in the two larger variances are roughly the same (∼ 0.1 %), while for the smallest variance is four times larger. This is a consequence of the experimental setup: the smallest variance is directly revealed only in one of the recorded projections used for the reconstruction.

Concluding remarks
In summary, we have presented a complete programme for the full polarization tomography of quantum states. Using the SU(2) Wigner function, we have provided an exact inversion formula in terms of the histograms of a standard Stokes measurement and derived a simplified version for very localized, high intensity states which turns out to be an inverse Radon transform. As a test of the theory, the reconstruction of an intense polarization squeezed state has been performed. A ML reconstruction algorithm has also been presented and has been compared to the direct method, thereby yielding an excellent agreement. Of course, the technique can be readily used for any other polarization state.