Tomographic reconstruction of the Wigner function on the Bloch sphere

We present a filtered backprojection algorithm for reconstructing the Wigner function of a system of large angular momentum j from Stern-Gerlach-type measurements. Our method is advantageous over the full determination of the density matrix in that it is insensitive to experimental fluctuations in j, and allows for a natural elimination of high-frequency noise in the Wigner function by taking into account the experimental uncertainties in the determination of j, its projection m, and the quantization axis orientation. No data binning and no arbitrary smoothing parameters are necessary in this reconstruction. Using recently published data [Riedel et al., Nature 464:1170 (2010)] we reconstruct the Wigner function of a spin-squeezed state of a Bose-Einstein condensate of about 1250 atoms, demonstrating that measurements along quantization axes lying in a single plane are sufficient for performing this tomographic reconstruction. Our method does not guarantee positivity of the reconstructed density matrix in the presence of experimental noise, which is a general limitation of backprojection algorithms.


Introduction
The reconstruction of the quantum-mechanical state of a system from measurements is an important topic of the emerging field of quantum technology [1]. Through partial or full state reconstruction we can estimate entanglement properties of multipartite quantum systems, and judge their usefulness for further experimental progress in fields such as quantum metrology [2,3,4,5,6,7,8,9], quantum simulation [10], and quantum computation [11,12,13,14,15].
Particularly in quantum metrology, experiments often involve large numbers of particles, and single-particle resolution is unavailable in both control and measurement. Because of this limitation, standard methods for reconstructing the quantum-mechanical density matrix [16,17,13] cannot be applied. For instance, and centrally to this work, in a Bose-Einstein condensate consisting of N atoms, with each atom representing a pseudo-spin-1/2 subsystem, the total spin length j = N/2 can take on very large values and the known reconstruction procedures become problematic. In a single Stern-Gerlach measurement on the atomic ensemble we measure the numbers of up and down spins N ↑ and N ↓ , in terms of which the total spin is j = (N ↑ + N ↓ )/2 and the projection quantum number is m = (N ↑ − N ↓ )/2. Since it is very difficult to determine the populations N ↑ and N ↓ with atomic accuracy [18,19], the density matrix, which requires knowledge of j, becomes impossible to reconstruct arXiv:1101.4131v2 [quant-ph] 10 Mar 2011 in full. Further, reconstructing the (2j + 1) 2 degrees of freedom of the density matrix [16,17,20] requires at least as many uncorrelated measurements, and therefore the experimental uncertainty in m will hinder this full determination. In the absence of reliable data, there will be significant uncertainty and noise throughout the density matrix in its Dicke representation ρ mm = jm|ρ|jm , which severely limits its usefulness. We need a method for calculating those components ofρ which are significant even in the presence of noise and for very large values of j, and a way of determining which components must remain unknown.
The Wigner function [21] is ideal for such a controlled reconstruction. It is a realvalued function on a sphere of radius j(j + 1), represented in terms of orthonormal Laplace spherical harmonics as [22] W (ϑ, ϕ) = 2j k=0 k q=−k where ϑ is the polar angle measured from the +z axis, and ϕ is the azimuthal angle around the z axis. While this sphere is commonly called a generalized Bloch sphere [4], its surface actually represents a two-dimensional phase space instead of a Hilbert space as for the original Bloch sphere. This Wigner function contains the same information as the density matrix for any spin-j system. While the marginals of the better-known Wigner function in planar space [21,23,24,25] are real-space or momentum-space probability distributions, the marginals of the spherical Wigner function are the projection quantum number distributions along all quantization axes [see (6) below]; further, the expectation value of the angular momentum vector is proportional to the "center of mass" of the Wigner function, × π 0 sin ϑdϑ 2π 0 dϕ{sin ϑ cos ϕ, sin ϑ sin ϕ, cos ϑ}W (ϑ, ϕ). Most importantly, the Wigner function allows us to differentiate between more significant components ρ kq (with smaller values of k) and more noise-prone components (with larger values of k) in a natural way. Further, if only components with k 2j are reconstructed, then accurate knowledge of j is not necessary. As detailed in section 2, the transformation from j-space (the Dicke representation ρ mm of the density matrix) to k-space (the spherical harmonic decomposition ρ kq of the Wigner function) proceeds though coupling coefficients which, at low k, are smooth in both j and m; this significantly reduces the impact of uncertainties in the experimental determination of (j, m).
Methods for reconstructing planar Wigner functions by inverse Radon transform are well established in the context of nonlinear optics [24,25]. In the past they have also been applied to tomographic data on large-spin quantum systems, locally approximating the Bloch sphere by a tangental plane and neglecting its curvature [6]. While this approximation is valid for spin states which are very localized on the Bloch sphere and do not wrap around it, future experimental progress is expected to produce ever more delocalized states (e.g., Schrödinger-cat states) whose properties are strongly influenced by the spherical shape of the Bloch sphere. Previous work on the reconstruction of the Wigner function on the full Bloch sphere has used the Husimi-Q distribution as input [26], which is the convolution of the system's Wigner function with that of a coherent state (see section 3). This convolution washes out features of the Wigner function that are smaller than a coherent state. Since the principal characteristic of spin-squeezed states is that their Wigner function possesses a peak width smaller than that of a coherent state, such a deconvolution-based reconstruction approach is ill suited for studying spin-squeezed states, which is the goal of much current research in atomic physics [2,3,4,5,6,7,8,9]. We therefore require a new method for reconstructing the complex quantum-mechanical states of large-spin systems from experimental data in the absence of simplifying circumstances, such as strong phase-space localization and/or lack of spin squeezing. This paper is organized as follows. In section 2 we present a novel filtered backprojection algorithm for reconstructing the Wigner function from experimental Stern-Gerlach data. Section 3 specializes this algorithm to data acquired with quantization axes lying in a single plane. Finally, section 4 applies the latter algorithm to a data set acquired in our group [6]. In what follows, a "single Stern-Gerlach measurement" describes a single determination of the projection quantum number m of a quantum system along a certain quantization axis. In our case this corresponds to a single run of state preparation and population determination of a two-component Bose-Einstein condensate, yielding a single tuple (j n , m n ). The equivalent for the original experiment [27] is sending a single silver atom through the experimental apparatus, and determining its deflection by the magnetic field gradient. On the other hand, a "Stern-Gerlach experiment" is a series of many single Stern-Gerlach measurements with fixed quantization axis, sufficient to determine the probability distribution {p −j , p −j+1 , . . . , p j } while j is presumed fixed.

Wigner function reconstruction by filtered backprojection
The density matrixρ of a system of total angular momentum j (assumed fixed here; this condition will be relaxed in section 2.1) is usually expressed in one of the two forms with the transformation coefficients (in the following simply termed Clebsch-Gordan coefficients) [22] t jmm kq = (−1) j−m−q j, m; j, −m |k, q , nonzero only if q = m − m . Both forms contain the same information and are completely interchangeable. While form (2a) is more common, form (2b) allows expressing the Wigner function on the Bloch sphere (1). Since our goal is the reconstruction of the Wigner function from experimental data, we focus on form (2b), in particular its low-k components. In order to determine the unknown quantum-mechanical state of a system of total spin j, it is necessary that many instances of this state can be generated experimentally [1], on which destructive measurements are performed. Further, projective Stern-Gerlach measurements must be performed along many different quantization axis orientations (ϑ n , ϕ n ).
For the correctness of the following reconstruction method it is crucial that these measured quantization axes are distributed as evenly as possible over the hemisphere of orientations. Since this requirement may be difficult to fulfill experimentally, we assign weights c n to the individual measurements in order for the weighted measurement density to approximate a homogeneous distribution of quantization axes as best possible. Notice that these weights are independent of the outcomes m n of the Stern-Gerlach measurements. In the ideal case of homogeneously distributed quantization axis orientations (for example through the vertices of a geodesic hemisphere), all these weights are chosen equal and the data are used most efficiently.
In this way, the results from M single Stern-Gerlach measurements along various quantization axes orientations are assembled into a data set of tuples (ϑ n , ϕ n , c n , m n ) with n = 1 . . . M and M n=1 c n = 1. Our filtered backprojection algorithm for reconstructing the Wigner function coefficients is then given by with D j m m (α, β, γ) = jm |e −iαĴz e −iβĴy e −iγĴz |jm a Wigner rotation matrix [28]; . This is formally equivalent to the filtered backprojection algorithm used for planar inverse Radon transforms [29], with the factor 2k + 1 representing the "filter", and the summand representing the backprojection. Our algorithm has all of the typical properties of planar inverse Radon transforms by filtered backprojection: no data binning is required, and there are no ad hoc parameters to be chosen or optimized. Further, as the backprojection algorithm is a direct sum and does not include an inversion step (such as a straight inversion of the Radon transform would require), the impact of experimental noise is bounded in the result. It is this last property which makes backprojection algorithms fast and reliable in practical applications such as X-ray computed tomography [29].
Our specific backprojection (4) can be interpreted in an intuitive way. The measured values of m n in the coordinate frame attached to the quantization axis (ϑ n , ϕ n ) are distributed according to the diagonal elements ρ mnmn and are converted from j-space into k-space via the Clebsch-Gordan coefficients t jmnmn kq with q = 0 (see section Appendix A for a numerical procedure). They are then rotated into the lab frame through the rotation matrices D k qq (ϕ n , ϑ n , χ n ) with the value of χ n irrelevant (set to zero) since q = 0.
In the following, we demonstrate that this algorithm (4) works in the limit of infinite data. If all quantization axis orientations have been used with equal frequency, and infinitely many measurements have been performed along each quantization axis, the sum over measurements M n=1 c n can be replaced by a normalized integral 1 2π π/2 0 sin ϑdϑ 2π 0 dϕ over the hemisphere of axis orientations (by symmetry the other hemisphere yields an identical result) and a sum over the measurement outcomes m, where the Stern-Gerlach probability distribution along a quantization axis (ϑ, ϕ) is given by the diagonal elements ρ mm of (2a) in the rotated frame, Using the orthogonality relations of Clebsch-Gordan coefficients, and spherical harmonics, it is easy to show that indeed ρ (fbp) kq = ρ kq , proving the validity of the reconstruction method in the limit of infinitely many homogeneously distributed Stern-Gerlach experiments.
In the more experimentally relevant case of a finite data set, the literature on the two-dimensional inverse Radon transform by filtered backprojection [29] indicates that excellent results can still be recovered, albeit with aliasing artifacts present to some degree. As a rough estimate, if Stern-Gerlach experiments are performed only along certain quantization axes spaced by an average angle ∆η, then the reconstructed partial waves of the Wigner function become unreliable for k k max = π/∆η. Further, if the number M of measurements is much less than the number of degrees of freedom (k max + 1) 2 , then the reconstructed coefficients ρ kq will be dominated by noise, in particular at large k. Both of these effects are mitigated in section 2.2 for the present reconstruction scheme.

Accounting for fluctuations in the total angular momentum j
We recall that for systems composed of many spin-1/2 components, such as twocomponent Bose-Einstein condensates, the total angular momentum j = (N ↑ + N ↓ )/2 often varies between single Stern-Gerlach measurements, as each such measurement requires the preparation of a new condensate. Instead of constructing a separate Wigner function for each occurring value of j, we notice that for k 2j the Clebsch-Gordan coefficients t jmm k0 depend smoothly on the total angular momentum j. This allows us to reconstruct the low-resolution part of the Wigner function even if j varies slightly between single Stern-Gerlach measurements. To this end we include the measured values of j in the data tuples, extending them to (ϑ n , ϕ n , c n , j n , m n ); the filtered backprojection formula is modified to Again we refer to Appendix A for a numerical method to evaluate this expression. The same smoothness of the Clebsch-Gordan coefficients at low k is used in section 2.2 to treat measurement uncertainties in both j n and m n in a perturbative manner in (9). This is fundamentally different from a direct tomographic reconstruction of the Dicke matrix elements ρ mm , where such uncertainties introduce large but correlated errors throughout the density matrix and make such a perturbative treatment impossible.

Measurement uncertainties and high-k damping
It is natural to assume that M uncorrelated experimental measurements can only serve to reconstruct M coefficients ρ kq , suggesting an upper limit k max ≈ √ M (assuming again a homogeneous distribution of quantization axis orientations). For larger values of k the angular power spectrum [30] C (fbp) k tends to acquire large fluctuations because of insufficient experimental data (see figure 3 for an example). However, simply cutting the reconstruction off at k max is unsatisfactory because it disregards that some useful information is still present in these high-k partial waves. A more natural cutoff is introduced through the kdependent sensitivity to experimental uncertainties. Assuming experimental variances of N 2 N /2 (with no covariance, jm = j m ) yield a leading order damping of the Clebsch-Gordan coefficients The rotation matrix elements are damped similarly: if the pointing direction of the quantization axis Ω = (ϑ, ϕ) has an uncertainty of σ Ω 1 {in terms of the expectation value of the angle η ΩΩ between the ideal axis orientation Ω and its true experimental value Ω we define σ 2 Ω = sin 2 η ΩΩ = 1 − [cos ϑ cos ϑ + sin ϑ sin ϑ cos(ϕ − ϕ )] 2 }, then for large k we find the rotation matrix elements to be damped as If σ N and σ Ω are equal for all measurements, the linearity of (9) yields a simple smoothing ρ kq → ρ kq e −αk(k+1) with α = . In this way, these two damping formulas (11,12) cut off the reconstruction at large k in a natural and smooth way.

Assembling the Wigner function
Inserting the resulting coefficients (9) into the form of the Wigner function (1) we find the tomographically reconstructed Wigner function where the contributions can be simplified to As is to be expected in spherical symmetry, the contribution of an individual Stern-Gerlach measurement (see figure 1) depends only on the relative angle cos η ΩΩn = cos ϑ cos ϑ n + sin ϑ sin ϑ n cos(ϕ − ϕ n ) between the quantization axis orientation Ω n = (ϑ n , ϕ n ) of the measurement and the point Ω = (ϑ, ϕ) on the Bloch sphere ( figure 2).
Similarly to technical implementations of the planar inverse Radon transform [29], the Wigner function is thus assembled from additive contributions due to the individual 0 π π/2 π/4 3π/4 η  Stern-Gerlach measurements (see figure 2 for an example). The constructive or destructive interference of these contributions is what yields the reconstructed Wigner function (see figure 4). The spatial resolutions of the Ξ jm (cos η) ultimately determine the spatial resolution of the reconstructed Wigner function: if the Wigner function is composed predominantly of contributions with m n ≈ ±j n its angular resolution is limited by that of a coherent state, ∆η 1/ j ; if on the other hand the majority of contributions has m n ≈ 0 the resolution can be significantly higher, ∆η 1/ j . We make use of this observation in sections 3 and 4, where a spin-squeezed state is reconstructed and the increased spatial resolution is critically important.

Positivity of the density matrix
It is well known that only positive semi-definite density matrices represent valid quantum-mechanical states of a system [1]. Unfortunately, the filtered backprojection method (9) does not assure that the reconstructedρ is positive semi-definite when used with a finite and noisy data set. For the purpose of displaying the Wigner function graphically, this is of no concern (see figure 4); however, when the tomographically reconstructed coefficients ρ (fbp) kq are used in quantitative calculations (see section 4) positivity can be crucial. This is a similar problem as the requirement for a positive absorption density in medical computed tomography (CT) imaging [29]. It is also present in many quantum-state reconstruction schemes, and has been discussed extensively in the quantum tomography literature [1].
We do not offer a solution for assuring the positivity of the reconstructed density matrix. Here we merely point out that in other reconstruction schemes, such as maximum-likelihood estimates [31], the ansatzρ =T †T forces the density matrixρ to be positive semi-definite; but a direct tomographic reconstruction ofT similar to (9) is currently lacking.

Quantization axes lying in a single plane
When the spin-j system's quantum-mechanical state is fairly localized on the Bloch sphere, not every choice of quantization axis orientation has the same potential for extracting information about the state. When the axis is close to parallel to the state, most Stern-Gerlach measurements will yield |m| ≈ j, with a limited angular resolution ∼ 1/ √ j given by the size of a coherent state on the Bloch sphere [26]. If the axis is close to perpendicular to the state, on the other hand, the distribution of measured values m represents the structure of the state's Wigner function much more accurately, with an angular resolution ∼ 1/j. This difference in scaling of the angular resolution, visible in figure 1, suggests that for large j it may be advantageous to focus on performing Stern-Gerlach measurements with quantization axes in a plane perpendicular to the quantum state, instead of covering the entire hemisphere of axis orientations. As a consequence much fewer measurements are needed, and we can get much more rapid convergence of the reconstruction in practice. But it is not a priori clear that this restriction of the quantization axes to a single plane has the potential for reconstructing the full quantum-mechanical state of the system.
As it turns out, a modification to the "filter" function in (9) results in a full reconstruction of the mirror-symmetric part of the Wigner function. Defining the coordinate system such that the state is localized near the +z axis and all quantization axes lie in the xy plane, the in-plane filtered backprojection formula is where (a) n = Γ(a + n)/Γ(a) is a Pochhammer symbol, and (a) . We again prove this reconstruction in the infinite-data limit. In the case of a homogeneous distribution of all azimuthal axis orientation angles ϕ we use the relationship 1 π π 0 dϕ[D k q 0 (ϕ, which remains true in the experimentally more relevant case of a finite number A of equally-spaced axis orientations (replacing 1 π π 0 dϕ → 1 A A−1 a=0 with ϕ = aπ/A) as long as k < A. Together with (6) and (7) we thus find that Thus in the infinite-data limit such an in-plane reconstruction exactly determines the coefficients ρ kq for which k + q is even, while giving no information on the coefficients for which k + q is odd. Since the parity of k + q is the z ↔ −z reflection parity of the spherical harmonics Y kq (ϑ, ϕ), the in-plane formula (15) reconstructs the positiveparity component W + (ϑ, ϕ) of the Wigner function W (ϑ, ϕ) = W + (ϑ, ϕ) + W − (ϑ, ϕ), with W ± (π − ϑ, ϕ) = ±W ± (ϑ, ϕ). If we know from other measurements that the state is fully localized on the "northern" Bloch hemisphere (z > 0), then the correct Wigner function is which has the decomposition ρ (fbp,P,N) kq in terms of the overlap integrals Υ q kk given in Appendix B. We conclude that the data acquired by Stern-Gerlach measurements with quantization axes lying solely within a plane are sufficient for an exact reconstruction of the Wigner function.

Measurement uncertainties and high-k damping
Measurement uncertainties can be introduced in (15) in the same way as in section 2.2. However, in an in-plane measurement series we can additionally separate out the azimuthal axis orientation uncertainty: since the rotation matrix elements D k q0 (ϕ, ϑ, 0) are proportional to e −iqϕ , a variance ϕ 2 − ϕ 2 = σ 2 ϕ leads to a damping D k q0 (ϕ,

Demonstration with experimental data
In this section we reconstruct the Wigner function from a data set describing ensembles of N = 1250(45) atoms acquired in our group [6]. In contrast to [6] we rotate the coordinate system such that all quantization axes lie in the xy plane and the state is localized around the +z axis; in this way the procedure of section 3 can be employed directly. The data set consists of three experimental runs spanning different ranges of ϕ with different angular resolutions, owing to the fact that the need for homogeneity in ϕ for the filtered backprojection algorithm (15) was not known at the time of data acquisition. We use weights c n adjusted such that the weighted density of Stern-Gerlach measurements is as close to homogeneous as possible over the range ϕ = 0 . . . π of azimuthal quantization axis orientations. As discussed in section 3 the planar arrangement of quantization axis orientations leads to a Wigner function which is peaked along both the +z and −z directions, featuring two identical copies of the quantum state. An additional Ramsey experiment [6] was used to experimentally determine the correct location of the state on the northern (z > 0) Bloch hemisphere.
High-k damping (section 3.1) is achieved with an experimental uncertainty of σ N ≈ 11 atoms (11) and with an experimental error model dominated by phase noise: (20), with phase noise amplitude σ ph = 8.2 • [6]. In figure 3 the effect of this damping is shown to be crucial for partial waves k 70.
The resulting reconstructed Wigner function is shown in figure 4. The highfrequency artifacts in the Wigner function far from the central peak are due to incomplete destructive interference of the contributions from the individual Stern-Gerlach measurements (see section 2.3). We expect that a more complete data set, including more quantization axis orientations, will lead to a smoother Wigner function at large angles ϑ.

Spin-squeezing measurement
We demonstrate the quantitative use of the reconstructed Wigner function by estimating the amount of spin squeezing in the system. Given a set of reconstructed Wigner function coefficients ρ kq we can calculate the probability distribution for the angular momentum projection quantum number onto any quantization axis orientation (ϑ, ϕ) from (6). In principle the variance V (ϑ, ϕ) = m 2 (ϑ, ϕ)−[ m (ϑ, ϕ)] 2 measures the amount of spin noise obtained experimentally. The expectation values of small integer powers of the projection quantum number m depend only on the low-k components of the Wigner function, which are particularly insensitive to experimental noise (11,12,20); in particular, In figure 5 we plot the resulting variances for quantization axes in the xy plane, and compare them to a coherent state centered on the +z axis. In the presence of imaging noise the variance of such a coherent state is given by (21) with ρ from data [6] from ρ kq from p m fits through (21). The solid line shows the results of Gaussian fits (figure 6) to the probability distributions pm( π 2 , ϕ) given in (6). As in [6] we first subtract the experimental noise (σ 2 N /2 with σ N = 11) from the calculated variances, and then divide by the variance of a coherent state, V coh = j /2 with j = 630 [see (22)]. A spin-squeezed state is characterized by negative values (in dB).
In figure 5 the experimental variance and the coherent-state variance are compared without imaging noise, i.e., the leading-order imaging noise contribution σ 2 N /2 is subtracted from the experimental variance before comparison with the variance of a noise-free coherent state. We notice that the resulting variance of the reconstructed state (dashed line in figure 5) is much larger than what was determined directly from the variances of the Stern-Gerlach data sets along the different quantization axes (open circles). We believe that this is a result of the lack of positivity of the density matrix (see section 2.4), owing to the finite and noisy data set used for its tomographic reconstruction. In fact, the probability distributions in figure 6 clearly show negative values, which strictly speaking render the reconstructed density matrix unphysical.
As an alternative extraction method we calculate the probability distribution p m (ϑ, ϕ) along a given quantization axis through (6) and fit it with a Gaussian curve (see figure 6); the variance of this fit then serves as an estimate of V (ϑ, ϕ). In such a fit the positivity of the p m is no longer a crucial ingredient. In figure 5 we show that this produces results that are very close to the variances calculated directly from Stern-Gerlach experiments along the various quantization axes. The deviations close to the squeezing maximum (ϕ s ≈ 6.7 • ) result from the fact that the reconstructed Wigner function contains contributions from all measurements, and therefore the extracted variance along a given quantization axis may be contaminated. Nonetheless the reconstructed Wigner function delivers a very concise picture of the structure of the multiparticle state, even for a data set with a non-uniform distribution over quantization axis orientations and with fluctuating values of j n . Further, with our method the variance V (ϑ, ϕ) can be calculated along any quantization axis orientation.
In practice any proof of spin squeezing will not proceed through the reconstruction of the Wigner function followed by either a fit to the projection (6) or a direct study of the projection noise (21). Instead, once the direction of squeezing ϕ s has been determined, a full Stern-Gerlach experiment will be performed along this axis in order to directly estimate the probability distribution p m (ϕ s ), as in [6] and in figure 5 (circles). In this way, problems associated with the positivity of the reconstruction (section 2.4) and with the influence of experimental data and noise from directions ϕ = ϕ s are strictly eliminated.
In future experiments providing data for the present tomographic reconstruction method, we plan to perform Stern-Gerlach measurements along many more quantization axes, but with as little as a single measurement per axis. Further, we will pay attention to cover the entire range of quantization axes uniformly [either the entire equator for (15) or the entire sphere for (9)] in each experimental run. In this way, we expect to need only minimal data preprocessing before reconstructing the Wigner function, and will be able to use the acquired data in the most efficient way by using equal weights c n = 1/M for all data points. We also expect that for such an improved data set the variance of the simple estimate given by (21) will be closer to that of the quantum-mechanical state.

Conclusions
We have presented a simple method for a tomographic reconstruction of the Wigner function of a spin-j system, applicable even to experimental settings where j is large and fluctuates between measurements. While the general procedure (9) requires Stern-Gerlach type measurements spread uniformly over all possible quantization axis orientations, a more specialized and faster procedure (15) determines the Wigner function using only a single plane of quantization axis orientations. We have shown that this latter procedure is capable of reconstructing the Wigner function of a spinsqueezed state from a recently published experimental data set [6].