{'fraction_non_alphanumeric': 0.04320565893653792, 'fraction_numerical': 0.014468871829910373, 'mean_word_length': 4.308726264134841, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Tomography of a number-resolving detector by reconstruction of an atomic many-body quantum state

The high-fidelity analysis of many-body quantum states of indistinguishable atoms requires the accurate counting of atoms. Here we report the tomographic reconstruction of an atom-number-resolving detector. The tomography is performed with an ultracold rubidium ensemble that is prepared in a coherent spin state by driving a Rabi coupling between the two hyperfine clock levels. The coupling is followed by counting the occupation number in one level. We characterize the fidelity of our detector and show that a negative-valued Wigner function is associated with it. Our results offer an exciting perspective for the high-fidelity reconstruction of entangled states and can be applied for a future demonstration of Heisenberg-limited atom interferometry.

High-fidelity preparation, manipulation, and detection of quantum states of many indistinguishable atoms have been greatly improved during the last decades.
Recent experiments creating entangled atomic quantum states in ensembles of indistinguishable atoms have reported counting noise that ranges from 3 atoms at a total number of a 600 atoms [25], to 1.6 atoms at 3000 atoms [26], to 10 atoms at 10 4 [27,28], to better than 17 atoms at 10 5 [14], and to 50 atoms at 5 × 10 5 [11]. An improvement of the counting noise below the value of the single atom, where the quantization of the atomic signal becomes apparent, promises a quantitative and qualitative improvement. For example, such a counting resolution would allow for the direct detection of Bell correlations between two separated atomic ensembles [29] and the observation of parity signals in Hong-Ou-Mandel-like interference experiments with many-particle states [30,31]. In metrology, a number-resolving counting can be applied to demonstrate a Heisenberg-limited resolution in atom interferometry [32,33]. A number-resolving counting has been obtained in a cavity-based detection [34], but is so far restricted to a discrimination between 0 and 1, and the scaling   [35], and was applied to extract local correlations. A number-resolved counting for up to 1000 atoms has been demonstrated in a millimeter-sized magnetooptical trap (mMOT) [36,37], but was so far not applied to the detection of many-body quantum states. The fine calibration of quantum measurement devices generally requires quantum detection tomography (QDT) techniques [38,39]. The goal of QDT is to provide a set of positive-operator-valued measures (POVM) that fully characterize the detector, beyond the assumption of projective measurements. So far, QDT has been mainly investigated for optical photocounting and homodyne detection [40][41][42][43][44] and very recently also applied to characterize qubit readout for pairs of trapped ions [45] and quantum computing machines [46]. Although QDT is necessary for quantum state preparation, control, reconstruction, and error mitigation, its potential has not yet been leveraged for the characterization of neutral-atom quantum systems. In this case, the technique is specifically promising, because of the small inherent detection loss compared to optical systems and the large achievable atom numbers compared to ion systems.
In this Letter, we apply a MOT-based number counting [47] to analyze the dynamics of a many-body spin state. We generate an atomic Bose-Einstein condensate (BEC) in one atomic clock level, apply a microwave (MW) coupling pulse of variable duration on the atomic clock transition, and count the number of atoms in the other, initially empty level. We are able to follow the time evolution of the coherent spin state with a clear resolution of the number quantization. By applying a stochastic matrix approach to the recorded histograms, we obtain a set of nonclassical POVM operators that fully characterize the detection process. The expected Poissonian distributions can be reproduced with a statistics-limited fidelity of up to 99%. We predict that the single-level detection is capable of operating an interferometric measurement with up to 7.8 dB squeezing-enhanced phase sensitivity gain, if the same total number can be provided with smaller fluctuations. This sensitivity gain can be further enhanced by increasing the mean atom number. Our detection capability promises a novel quality for the analysis of entangled quantum states and demonstration of Heisenberg-limited interferometry.
In our experiments, we generate a BEC of 10 5 87 Rb atoms in a crossed-beam optical dipole trap with a preparation time of 3.3 s. The details of the BEC production can be found in Ref. [48]. We prepare the BEC in the hyperfine level |F, m F = |2, 2 and reduce the number of atoms to enter the regime of our number-resolved counting. The reduction is realized by transferring 34 atoms to the level |1, 1 , on average, and a subsequent optical removal of the residual atoms in the F = 2 manifold. A further MW pulse transfers the remaining atoms to the level |2, 0 . A final resonant light push on the F = 1 manifold terminates our state preparation with 34 and 0 atoms in the clock states |2, 0 and |1, 0 , respectively. The total number of particles fluctuates by 6.4 atoms, which is dominated by projection noise (5.8 atoms).
We apply a resonant MW pulse on the clock transition with a variable duration ranging from t = 2.5 µs to 56 µs. The many-body state in the pseudo-spin-1/2 system can thus be represented by a coherent spin state (CSS), with maximal total spin, but variable rotation angle θ. The analysis of the CSS is based on counting the number of atoms in the level |1, 0 . To this end, the atoms in the level |2, 0 are removed and the remaining atoms are counted by fluorescence detection in the mMOT. The quality of the state analysis thus depends on the efficiency of the removal procedure. Therefore, the detection process starts with a strong reduction of the atomic density by switching off one of the two dipole trap laser beams. A σ + -polarized light push at a magnetic field of 6.7 G quickly pumps all F = 2 atoms into a closed cycling transition reducing the probability to fall into a nonresonant state. The resonant atoms are accelerated and leave the trap, while the probability of unwanted collisions is reduced by the low density. The removal of atoms in the level |2, 0 has a finite extinction ratio of 42.4 dB, resulting in an unwanted, Poisson-distributed remainder of 0.27 atoms maximally. These atoms are produced by two processes: (i) They escape the removal process to the level F = 1 because of imperfect optical pumping. (ii) They are captured from the background gas, which is temporally increased after the operation of the two-dimensional magneto-optical trap. We detect the remaining atoms in the mMOT setup, consisting of a magnetooptical trap with millimeter-sized illumination beams [47,48]. The optical dipole trap is switched off to start an equilibration phase in the mMOT of 50 ms, during which the magnetic fields settle and the atoms are compressed and cooled. Subsequently, the main atom counting signal is obtained by collecting fluorescence light for 65 ms with a charge-coupleddevice (CCD) camera. Finally, a second image without atoms is recorded for background subtraction. The spin preparation and detection processes require a total of 1.8 s. After nine measurement runs, the system is halted for 60 s to avoid a slow increase of the mMOT capture rate from the background gas. Fig. 1 (a) shows a sketch of the experimental setup including the mMOT and ODT beams and the high-numericalaperture detection objective. Fig. 1 (b) shows a time trace of 100 consecutive number measurements in |1, 0 shown for three different MW pulse lengths. The measured number of atoms accumulate at integer numbers, enabling a number assignment fidelity ranging from 99.7% at 1 atoms to 99.0% at 15 atoms [48]. Fig. 2 (a) shows the mean number of the transferred atoms as a function of the microwave pulse duration. The mean atom number follows a sinusoidal Rabi oscillation with a Rabi frequency Ω = 2π × 8.2 kHz (see below). Figs. 2 (b)-(f) present the exemplary histograms, which can be associated to rotation angles θ = Ωt. Without rotation (b), the distribution shows the detection of recaptured atoms, which can be treated as statistical dark counts in the detection system. For finite rotations (c)-(f), the distributions shift to higher atom number and increased width.
Under the assumptions that the microwave generates a ho- mogeneous coupling to the cloud and that the level |1, 0 is initially completely empty -that are both very well fulfilled in our case -we can employ the recorded data for QDT. We associate the two clock levels |1, 0 and |2, 0 with the letters a and b, respectively, to simplify the notation. We model the detection by expressing the probability of a measurement result n (the number of atoms in clock level a) as in terms of a stochastic matrix (Markov mapping) V with non-negative elements V n,m ≥ 0, which satisfies the normalization property n V n,m = 1 for all m. Physically, the quantity V n,m can be interpreted as the probability to measure n atoms if m atoms reach the detector. We use the ideal probability P id (m|t) = Tr[|m m|ρ(t)], whereρ = +∞ N=0 ρ N |N b |0 a 0| a N| b is the generic atomic state before starting the dynamics,Û(t) = exp − iΩ R t(â †b +âb † )/2 describes the Rabi coupling andρ(t) =Û(t)ρÛ(t) † . The assumption that the initial state is diagonal is well justified experimentally. The matrix V provides a full characterization of the detection process, including finite resolutions and biases. It should be noticed that Eq. (1) can be rewritten as P V (n|t) = Tr ρ(t)Π n , in terms of a POVM set {Π n }, wherê V being positive semi-definite guarantees thatΠ n ≥ 0, while the condition n V n,m = 1 for all m guarantees the completeness relation nΠn = 1.
Our QDT protocol consists of finding the coefficients ρ N , V n,m and Ω R that minimize a cost function C = j d 2 H (t j ). This is given by the sum over all times t j of the squared statistical distance d 2 H (t j ) [49] between the probability distribution P V (n|t j ) and the experimental histogram P exp (n|t j ) [e.g. Fig. 2 The constrained minimization of Eq. (3) is performed with a gradient descent algorithm [50]. We emphasize that a reliable characterization of the POVM set requires a substantial overlap between probability distributions at different times in order to avoid overfitting [50]. For our experimental parameters, see Fig. 2, this is guaranteed for 0 ≤ n 20.
In Fig. 2 we compare the experimental histograms (bars) with the probabilities derived from the QDT (orange circles), namely Eq. (1) with V, ρ N and Ω R calculated using the minimization algorithm. The agreement is excellent, as the obtained probability distribution P V (n|t j ) achieves a very high fidelity with P exp (n|t j ), for all t j (notice that the iterative optimization algorithm is stopped when C = 0.01, which is a value close to saturation [50]). The histograms are consistent with that calculated with a binomial distribution (green squares). The latter assume a Gaussian distribution of the total number of atoms with measured mean and standard deviation, ideal Rabi transfer and the convolution with a binomial distribution with a mean number of 0.27 atoms to account for the unwanted detection of background atoms. The mean number of atoms as a function of time, n P V (n, t)n for the reconstructed V, ρ N and Ω R interpolates well the detec-tion events as shown in Fig. 2(a). It should be noticed that the Rabi frequency extracted from the tomographic reconstruction Ω R = 8.2 ± 0.2 kHz agrees precisely with the result of a sinusoidal fit to the data.
In Fig. 3 we show the results of our joint detection and state reconstruction. Figure 3(a) shows the elements V n,m . For most values of m, the weights V n,m concentrate around the diagonal n = m, where they reach their maximal value. In other words, if m atoms reach the detector, the most probable event is to detect n = m. The probability of such detection events is quantified below. For instance, in panel (b) we show V n,m as a function of n and for the specific values m = 0, 5 and 10: the histograms are cuts of the plot of panel (a). For m 20, the reconstructed V n,m spreads away from the diagonal. Here, the QDT becomes uncertain because the recorded probability distributions do not overlap sufficiently and the optimization method is affected by overfitting of the data [51]. To recognize the overfitting effect, we have performed a 'learning test', see Ref. [50], in which the experimental histogram at time t j is compared with the reconstructed P V (n|t j ), where the coef-ficients Ω R ,ρ, and V are calculated from the minimization algorithm using all the experimental data except those at time t j . In this case, we observe a fidelity between P V (n|t j ) and P exp (n|t j ) above 99% for times t j up to 18.48 µs [50]. In the future, a QDT at larger atom numbers can be obtained by taking more histograms with larger statistics. In panel (c) we show the reconstructed elements ρ N as a function of the number of particles. As we see, the reconstructed diagonal state has approximately a Gaussian shape with mean N = 35.4 and root mean square error ∆N = 6.4 ≈ N 1/2 .
For m 20, as V n,m is strongly peaked around n = m. It is thus convenient to calculate P(n − m) = m V n−m,m P m , giving the reconstructed probability that n − m particles are detected if m particles hit the detector. Here, P m = j P id (m|t j ) is the overall probability that m particles hit the detector, for the considered measurement times t j and takes into account the most likely detection events and it is almost negligible for m 20. For a noiseless detector P(n − m) is a delta peak at n = m, regardless the P m distribution. In the case of our noisy detector, P(n − m) is still strongly peaked at n = m, with an overall probability of about 70%, see Fig. 3 (d). The slight asymmetry of the distribution V n−m reflects the unwanted recapture of atoms described above, which biases V n−m to positive values of n − m. By calculating the variance of the P(n − m) distribution for n ≤ m (thus not affected by the atom recapture) we can extract a detection sensitivity σ = 0.4 ± 0.02. This counting uncertainty is larger than the uncertainty obtained from Fig. 1, because the finite number of measurements additionally deteriorates the QDT.
As shown in Eq. (2), accessing the matrix V allows us to characterize the POVM elementsΠ n . For instance, in Fig. 3(d), we plot the Wigner distribution of the reconstructed POVM operatorsΠ n , W n (x, p) = m V n,m W(x, p; m), where W(x, p; m) is the Wigner function of the Fock state |m [52]. For n = 0 the Wigner function is positive, as expected, corresponding to the detection of vacuum. On the contrary, for n ≥ 1, the Wigner functions W n (q, p) have negative values, indicating the absence of a classical analogue of these operators. To emphasise the fundamental quantum nature of our detection we notice that a POVM with negative Wigner function is necessary to prove Bell's non-locality with Gaussian states (which have positive Wigner distributions) [40,53].
Measuring the atom number in a single level, as done in our experiment, still allows to surpass the standard quantum limit of phase sensitivity, for instance provided that (i) each fixed-N state is spin squeezed and (ii) the distribution of the total atom number ρ N has sufficiently low fluctuations. Using the results of our QDT and the atom number distribution, we can predict an optimal gain of about 1.3 over the standard quantum limit, when squeezing the relative atom number distribution [50]. A higher gain, up to 7.8 dB is possible when also reducing ∆N, for our detection andN ≈ 36. The gain further increases when increasingN.
In summary, we have employed a number-resolving detector to analyze the dynamics of a coherent spin state derived from an atomic BEC. We have characterized the detection process by the simultaneous reconstruction of the diagonal quantum state and the detector's POVM operators. The latter are characterized by negative Wigner functions, thus unveiling the inherent quantum nature of the detector. In the future, the presented detector and the developed QDT techniques will be directly extended to entangled many-body states, promising the detection of entanglement with unprecedented fidelity in the regime of up to 100 atoms.
Tomography of a number-resolving detector by reconstruction of an atomic many-body quantum state (Supplemental Material) In the Supplemental Material we provide more details about the optimization algorithm for the quantum detection tomography as well as a detailed analysis about the possibility to reach sub-shot noise sensitivities with the noisy detector and single-level counting.

Optimization algorithm
The minimization of the cost function C = j d 2 H (t j ), with d 2 H (t j ) given by Eq. (3) of the main text, is done using a gradient descent algorithm. The gradient descend technique is a widely used method to find the minimum of a cost function with respect to a large number of free parameters. It is typically used in neural networks and variational quantum algorithm. Starting from an initial guess, the parameters are updated iteratively by following the direction of the negative gradient of C. Notice that, in place of d 2 H (t j ) in the cost function, we could have used the the Kullback-Leibler divergence between P V (n|t j ) and P exp (n|t j ), obtaining similar results.
In our case, the minimization is divided in different steps. The numerical optimization provides the stochastic matrix V, the coefficients ρ N and the Rabi frequency Ω R . First, we run a number of iterations to minimize V, followed by a minimization of ρ N , and finally, of Ω R . The procedure is then repeated iteratively until C = 0.01, which is set as a convenient cutoff. This cutoff is very close to saturation, meaning that increasing the number of optimization steps of the algorithm, C can reach values slightly smaller than 0.01 at the price of an rapidly increasing computation time.
The minimization algorithm does not require any input parameter. It relies on assuming the initial state in the diagonal form (see main text) and that the transformation follows a unitary evolution. Both these assumptions are well justified experimentally. The number of particles domain is cut to the maximum detected value, n max (equal to 59 in our case), such that V is a square matrix of dimension (n max + 1) 2 . We impose the following constraints at each step of the algorithm: i) V ≥ 0; ii) n V n,m = 1 for all m; iii) ρ N ≥ 0 for all N; and iv) N ρ N = 1. As initial conditions we take a random and normalized (n max + 1) × (n max + 1) matrix V, a Gaussian distribution for ρ N and an initial guess for Ω R . The center of the Gaussian number distribution and the Rabi frequency are extracted from a fit of the mean. number of atoms N(t) as a function of time t, according to the function N(t) = N sin 2 (πΩ R t) + C off , with N, Ω R and the offset C off as fitting parameters. The width of the Gaussian, ∆N is extracted from a fit of N 2 (t) as a function of time according to N 2 (t) = 4 k=0 b k sin k (πΩ R t), where b k are fitting parameters. We identify b 4 = N 2 . Notice that for these initial conditions the initial cumulative fidelity is very small, F ≈ 0.5. We have verified that the output of the gradient descent algorithm does not depend on the choice of the initial random matrix V.
Error bars on quantities of interest are obtained from a resampling of the experimental probability distribution. Specifically, given the reconstructed V, ρ N , and Ω R , we generate about 20 new sets of n max resampled data for each time t j ( j = 1, ..., 8), mimicking re-acquisitions of experimental data sets. For each data set we run the gradient descent optimization. Overall, this procedure provides a set of 20 matrices V, ρ N , and Ω R , from which we extract error bars.

Learning test
Here we provide a "learning test" of our method. For this, we reconstruct Ω R ,ρ, and V using the experimental data collected at all time except at time t j . Numerical methods and initial conditions are the same as in the previous analysis. The reconstructed quantities (labelled with index j, Ω ( j) R ,ρ ( j) , and V ( j) ) are used to calculate the expected output probability P V j (n|t j ) at the time t j not seen by the minimization algorithm. In Fig. S1(a) we compare the probability P V j (n|t j ) (red circles) with the experimental histograms P exp (n|t j ) (grey bars), at different times. While at small times (e.g. at t = 0 and t = 3.64 µs), the method is able to learn the unseen distribution, this capability is lost for large times. This is due to the fact that, at large times, distributions do not overlap substantially: the algorithm simply "overfits" the experimental distributions without been able to generalize to unseen distributions. To be more quantitative, we calculate the fidelity between the experimental distribution and the reconstructed one. Results are shown in Fig. S1(b).  FIG. S1. Learning test. We apply the tomography method using all data expect at time t j in order to guess the experimental distribution at time t j . (a) The learned distribution P V j (n|t j ) (red circles) is compared to the experimental distribution (grey bars), at different times. In panel (b) we show 1 − F (t j ), where F (t j ) is the fidelity between the experimental distribution P exp (n, t j ) and P V j (n|t j ), Eq. (S1).

Fisher information
As a further consistency check of our results, we calculate -using the reconstructed distribution P V (n|t) -the Fisher information with respect to the parameter θ = Ω R t: The quantity F V (θ)/N is shown by the solid black line in Fig. S2. The grey region are fluctuations obtained via bootstrapping. The Fisher information is easily calculated using ideal probabilities, giving F id (θ)/N = cos 2 (θ/2) and shown as dashed blue line in the figure. Equation (S2) follows the ideal behaviour for relatively large values of θ, where the influence of (small) detection noise is less relevant.

POSSIBLE OBSERVATION OF SUB SHOT NOISE SENSITIVITIES
The aim of this subsection is to evaluate the possibility to reach a sub-shot noise phase sensitivity by using i) a superposition of spin squeezed states, ii) measuring the number of particles in one atomic level, and iii) using the noisy detectors characterized in this work. We model the probe state as a superposition, of Gaussian states with different number of particles, where N(N) provides the normalization and |µ x are eigenstates of the operatorĴ x with eigenvalue µ. HereĴ x = (â †b +b †â )/2 andâ (b) are annihilation operators for the two modes a and b corresponding to two atomic hyperfine levels. The quantity s is a parameter directly related to the variance ofĴ x : 4(∆Ĵ x ) 2 = s N ρ N N. The case s = 1 approximates well, up to corrections of the order 1/N, a superposition of coherent spin states. When s < 1, each state |ψ(s, N) is spin squeezed. Following the experimental results, here we model ρ N as a Gaussian distribution with meanN and root mean-square fluctuations ∆N (we recall thatN ≈ (∆N) 2 ≈ 36 in our experiment). Given the stateρ, Eq. (S3), we study the sensitivity to phase shift estimation when applying the transformation e −iθĴ y . The sensitivity is calculated via error propagation when measuring the number of particles in one atomic level only, namely eventually taking into account the finite detection resolution. When optimized over θ, the sensitivity depends strongly on the squeezing parameter as well as the finite width of ρ N . In Fig. S3(a) we plot the sensitivity gain over the shot noise, namely G = 1/[N(∆θ) 2 ] as a function of ∆N and s. Here we consider N ≈ 36, as in the experiment, and use the reconstructed experimental detection. The white region corresponds to G ≤ 1. We see that the gain increases when decreasing ∆N and s, and can reach G ≈ 6 for our relatively small number of particles and finite detection resolution. In Fig. S3(b), we plot the gain as a function of s, for N ≈ 36 and fixed (∆N) 2 = 36, as in the experiment. The blue line is the gain obtained with ideal detection, while the red line is obtained for our noisy detector (namely obtained by using our reconstructed detector matrix V). For s = 1 (no squeezing), we find G ≈ 0.75, consistent with the analysis of the Fisher information shown in Fig. S2. In Fig. S3(c) we plot the gain as a function of N, optimized over s, and for ∆N = N. The blue line (obtained for ideal detection) follows G ∼N 0.25 . The red line (obtained for our noisy detection) tends to approach the blue line for N 1.