Optimal time-resolved photon number distribution reconstruction of a cavity field by maximum likelihood

We present a method for reconstructing the average evolution of the photon number distribution of a field decaying in a high-Q cavity. It applies an iterative maximum likelihood state reconstruction algorithm to the diagonal elements of the field density operator. It is based on quantum non-demolition measurements carried out with atoms crossing the cavity one by one. A small set of successively detected atoms defines a positive operator valued measure (POVM). The reconstruction is performed by applying this POVM to a large ensemble of field realizations. An optimal POVM based on the detection of a minimal number of atoms is shown to be sufficient to ensure an unambiguous convergence of the reconstruction. The cavity crossing time of this minimal number of atoms must be much shorter than the lifetime of the largest photon number present in the field. We apply the method to monitor the evolution of number states prepared by quantum feedback in a recent experiment. The method could also be useful in circuit QED experiments.


Introduction
Reconstruction of quantum states [1] is an important issue in quantum information science. Reconstruction of quantum oscillator states is particularly useful in quantum optics [2][3][4][5], trapped ion physics [6,7], cavity [8] and circuit [9] quantum electrodynamics. It is a powerful tool for studying decoherence [8,10] of systems at the boundary between the quantum and the classical worlds. For states evolving rapidly in time such as non-classical fields trapped in cavities, the reconstruction must be based on information acquired on an ensemble of realizations, each of them providing, in a short time interval, only a limited amount of information. Adequate estimation procedures must be implemented in order to minimize both the time window in each run and the number of realizations required for a faithful state reconstruction.
In many cavity and circuit QED experiments, in which the field is stored in a high-Q resonator (damping time T c ) a partial reconstruction of the state restricted to the photon number distribution provides all the relevant information. This is the case for experiments preparing and studying the evolution of photon number states (Fock states) in cavities [9,11,12]. In these experiments, the time window for data acquisition in each run must be much shorter than the decoherence time T c /n m of the Fock state corresponding to the maximum number of photons n m involved in the measured field. We have recently used a photon number distribution reconstruction method based on the maximum likelihood (MaxLik) procedure [13,14] for process tomography of Fock state relaxation [15] and for the assessment of the fidelity of a quantum feedback procedure preparing on demand and stabilizing Fock states in a cavity [11,12]. A microwave field stored in a superconducting cavity is dispersively probed by nonresonant circular Rydberg atoms. The detection of the final state of each atom provides quantum non-demolition (QND) information on the photon number [16].
In this paper we describe in detail this efficient reconstruction method whose principle was only briefly outlined in [15]. We show that detecting in each realization a minimum number of atoms in a time short compared to T c /n m is sufficient to reconstruct faithfully the field evolution. In section 2, we recall the basics of our experimental setup and of the QND photon number measurements. In section 3, we present the reconstruction algorithm, and discuss in section 4 the required minimal atom number per realization. In section 5, we confirm the results of section 4 by numerical simulation of an ideal experiment. Realistic experimental parameters are taken into account in the numerical simulations of section 6. In section 7 we apply the method to the actual data of the feedback experiment presented in [11] and present some concluding remarks in section 8. Figure 1 is a sketch of the experimental setup described in detail in [11]. A microwave field is confined in a high-finesse Fabry-Pérot cavity C made of two superconducting niobium mirrors, resonant at ω c /2π = 51 GHz. When cooled down to 0.8 K, it has a very long field energy lifetime (T c = 65 ms in [11]), with a low average number of thermal photons n th = 0.05.

The experimental setup
The microwave field is probed by a succession of rubidium atom samples, prepared in box B into the circular Rydberg states g or e with principal quantum numbers 50 and 51, respectively. Each atom undergoes π/2 microwave pulses mixing e and g in the low-Q cavities R 1 and R 2 fed by the classical source S, which constitute a Ramsey interferometer. Atomic samples are separated by a time interval T a = 82 µs. They are excited by 2 µs laser pulses and have a well-defined velocity v = 250 m s −1 so that the position of each sample is known within 1 mm at any time from its preparation in B to its detection (in the e-g basis) by field ionization in D. The number of Rydberg atoms in a single sample obeys a Poisson law with an average ≈0.8. This number is low enough so that the probability of having more than two atoms in a sample is negligible. The atomic detection efficiency is 35%.
The Rydberg atoms interact with the field in C with a coupling strength measured by the vacuum Rabi frequency, 0 /2π = 48 kHz. The g → e frequency is ω a = ω c + δ with δ/2π ≈ 250 kHz. The atom-field interaction is thus dispersive (δ 0 ). The atomic coherence undergoes in C a light-induced phase shift linear in the photon number n: ϕ(n) = (n + 1/2)ϕ 0 . The phase shift per photon ϕ 0 is where t eff = √ π/2 w 0 /v is the effective interaction time for an atom crossing the Gaussian profile of the cavity mode with a waist w 0 . The measurement of ϕ(n) by Ramsey interferometry provides information about the photon number distribution of the cavity field [17]. From now on, we will consider that ϕ 0 = 2π/(n m + 1), which is optimal for discriminating photon numbers ranging from 0 to n m .
Assuming the field to be in the n-photon Fock state |n , the probability to detect an atom in state j ( j = 0, 1 for e and g, respectively) is ideally given by where the Ramsey interferometer phase φ can be tuned by varying the relative phase between the two π/2-pulses. We introduce the photon number operatorn and we define for any function f (n) the operator f (n) = n f (n) |n n|. The positive operator valued measure (POVM) associated with the atomic detection in state j iŝ For a field with a photon number distribution P(n) = n|ρ|n (ρ is the field density operator), the probability π( j|P, φ) for detecting the atom in state j is More generally, detecting successively N D atoms in states J = { j k } 1 k N D , with the respective Ramsey phase settings = {φ k } 1 k N D , corresponds to the measurement of the POVMÊ Note that the single-atom POVMs π( j k |n, φ k ) commute since they are all diagonal in the photon number basis. The probability of finding the detection event J with phase setting is then π(J |P, ) = Tr ρÊ ,J .
After this detection, the field density operatorρ is changed tô For a density operator that is diagonal in the Fock states basis, this equation reduces tô leading for the photon number distribution to where Z is a normalization factor. We recognize here a direct application of Bayes law expressing how the photon number distribution is updated once the N D atoms have been detected.
We have shown in [16] that this measurement projects the field into a Fock state when N D is of the order of n 2 m . These atomic detections achieve thus an ideal projective QND measurement of the photon number. The result of this projection is random. By performing many such realizations, we could reconstruct the initial photon number distribution P(n). In the conditions of [16], where n m = 7, this progressive collapse of the field state takes a time of the order of 20 ms, which is not negligible as compared to T c /n m . This direct determination of P(n) leads thus to a skewed probability distribution for large n, reducing the fidelity of the state estimation at a given time. In order to improve this fidelity, we used instead in [11,12,15] a MaxLik reconstruction method using many realizations of incomplete projections with N D n 2 m atoms detected in a time interval much shorter than T c /n m .

The reconstruction process
The MaxLik reconstruction method seeks, in general, a density operatorρ maximizing the joint likelihood of the results of all measurements carried out on M identical copies of the system. We are here interested in the diagonal elements ofρ only, but for the sake of generality we first describe the reconstruction method in terms of the full density operator. For each copy of the field, one performs N D -atom detections whose outcomes are represented by a sequence J of N D atomic detection results. Let us call M the number of times the phase setting has been used and M ,J the number of times J has been obtained for this setting. We denote by f ,J = M ,J /M the corresponding frequencies.
If all M ,J are very large, the frequencies f ,J are equal to the probabilities p ,J = (M φ /M)Tr(ρÊ ,J ) = Tr(ρÊ ,J ), where we have defined the POVM normalized by ,JÊ ,J = 1. Ideally, provided that the set of measurements is large enough, one could retrieveρ by inverting the linear equations In practice, the M ,J numbers are finite, leading to statistical fluctuations of the f ,J and this simple reconstruction procedure may lead to unphysical density operators with negative eigenvalues. The MaxLik method [13] circumvents this difficulty by searching for the density operatorρ maximizing the likelihood function Using the results of [18,19], an iterative algorithm leading to the density matrix maximizing L has been proposed in [13]. It introduces the nonlinear function ofρ, which obviously reduces to the unity operator ifρ satisfies equation (12). The searched density These expressions suggest to obtainρ f as the fixed point of the iterations [14] ρ i+1 =R(ρ i )ρ i (15) or, within a normalization, starting from an initial guess, for instance a flat photon number distribution. The iterative expression (16) ensures the positivity of the density operator at each step and is thus preferred in the general case. For diagonal density operators, which will only interest us from now on, the simpler iteration (15) is well behaved (i.e. yields at each iteration positive P i (n)s). Moreover, it increases the likelihood at each step. The density operatorρ i converges towardρ f if it is the unique fixed point of (15). We determine in the next section the minimum number N D of detected atoms in each realization, ensuring the uniqueness of the fixed point. Before this discussion, let us note that the iteration (15) has a simple physical interpretation. The first step of the iterative procedure replaces our initial guessρ 0 bŷ where we recognize the projected density operator updated by the measurement , J (see equation (9)) averaged over a large ensemble of realizations (weights f ,J ). If N D were large enough to ensure the complete collapse of the photon number distribution in each realization, ρ 1 would obviously be the searched distribution and the procedure would converge in a single step [16]. For smaller N D values,ρ 1 , which includes information provided by many measurements, is a better guess of the actual state thanρ 0 . Resuming then the iteration defined by equation (15) withρ 1 as a new initial guess leads to an estimationρ 2 better than ρ 1 and so on.

The minimum number of detected atoms per realization
Let us now examine how many atoms have to be detected in each realization in order to ensure the uniqueness of the MaxLik iteration fixed point toward which the state estimate converges. For this purpose, it is convenient to introduce a simple geometrical interpretation of the detection probabilities π(J |P, ). We associate with each photon number n from 0 to n m a unit vector n in an O x y-plane at an angle ϕ(n) with the O x-axis. Similarly, for each detection in state j with a Ramsey interferometer phase φ, we define a unit vector d at an angle −φ + jπ with the O x-axis. For a single atom detected per realization (N D = 1), we then obtain from equation (2) π( j|P, φ) = with P = n P(n)n. We recognize here the dipole of n m + 1 charges distributed on a circle at the positions defined by the n vectors as represented in figure 2. The total charge n P(n) of this distribution is 1.
For two atoms per realization (N D = 2) corresponding to the detection vectors d 1 and d 2 , we similarly obtain with Q (1) = P, D (1) = d 1 + d 2 , and where the rank 2 tensors Q (2) and D (2) are defined by their components and their scalar product by The tensor Q (2) is the quadrupole of the charge distribution. While the dipole is sufficient to determine the one-atom probabilities, dipole and quadrupole are required for two-atom joint probabilities.
The quadrupole of the photon number distribution is thus explicitly determined in terms of simple linear combinations of two-atom joint probabilities. These equations generalize to an arbitrary rank : with = {(π/2)δ p,y , (π/2)δ q,y , . . . , (π/2)δ t,y } and J = { j p , j q , . . . , j t } such as j p + j q + · · · + j t is even. Let us now establish the main result of this section. A photon number distribution over a range of n m + 1 values requires at least the measurement (with a minimum of two Ramsey interferometer phases) of N min D = I [(n m + 1)/2] atoms (I denotes the integer part). We first note that all multipoles up to rank of a charge distribution localized on a circle whose total charge is unity (figure 2) are determined by 2 real parameters. For = 1, the two parameters are the two coordinates of the dipole. Going to = 2 adds the three components Q (2) x x , Q (2) yy and Q (2) x y which are constrained by the identity Q (2) x x + Q (2) yy = 1 resulting from the unit length of the n vectors. The total number of parameters is thus 4. Let us assume that the multipole expansion up to rank − 1 is determined by 2 − 2 parameters. The Q ( ) multipole has + 1 components constrained by the − 1 relations Q ( ) x x pqr ... + Q ( ) yypqr ... = Q ( −2) pqr ... which again result from the normalization of the n vectors. It can be shown that there are no other independent relations linking the Q ( ) components to the Q (s) ones with s < . The number of independent parameters determining all multipole components up to rank is thus 2 − 2 + 2 = 2 . The information contained in the measurements realized with N D atoms thus results in 2N D constraints for the n m + 1 charges whose sum is unity, which demonstrates that the minimum number of detected atoms should be N min D = I [(n m + 1)/2]. The above argument shows that the MaxLik strategy converges if we have perfect knowledge of the N min D -atom joint probabilities. We examine in the next section the rate of convergence and the influence of a finite number of realizations leading to statistical noise on the measured probabilities.

Monte Carlo simulations, the ideal case
We have numerically checked the convergence and the accuracy of the MaxLik procedure for different N D values. We perform Monte Carlo simulations of QND measurement sequences and feed the results of these simulations into the iterative reconstruction algorithm. An initial photon number distribution P f (n) (shown by the green histogram in figure 3(a)) with n m = 7 is arbitrarily chosen. It has n m + 1 = 8 non-zero elements spanning the photon numbers from 0 to n m . The measurement is simulated using a phase shift per photon equal to π/4 and four Ramsey phases (φ 0, π/4, π/2, 3π/4). For each atom, φ is chosen randomly among these four values. Although two phases are, in principle, enough according to section 4, the choice of four phases ensures that the number of iterations required for convergence is the same for all eight photon numbers. These simulations are performed assuming an ideal cavity with negligible field decay and a perfect Ramsey interferometer (fringe contrast equal to 1). The influence of finite decay time and limited Ramsey fringe contrast will be taken into account in the next section.
A single realization of a simulated measurement consists of N D detected atoms. The number of realizations leading to a reconstruction is M = 19 000. Thirty independent simulations are performed in order to estimate the dispersion of the final result. Figures 3(b)-(d) display the evolution of the eight P i (n) values versus the iteration rank i, averaged over the 30 simulations, for N D = 3, 4 and 6. An initial flat photon number distribution is assumed (P 0 (n) = 1/8). The dotted lines indicate the values of P f (n). For N D = 3 (figure 3(b)) the field state seems to stabilize at a first fixed point different from P f (n) after several hundreds of iterations. It then evolves after i max = 100 000 iterations into the distribution shown in figure 3(a) (red histogram). The error bars represent the variance of the 30 independent simulations. The occurrence of a 'metastable' point in the iteration and the poor fidelity at i = i max (compare the red to the green histogram in figure 3(a)) is consistent with the analysis of section 3, which states that N D = 3 is not enough for an unambiguous reconstruction of P f (n) with n m = 7. Figures 3(c) and (d) present in the same way the evolution of P i (n) for N D = 4 and 6. The P i (n)s reach a stable point after about 300 iterations for N D = 4 and about 60 iterations for N D = 6. These fixed points remain unchanged thereafter until i = i max . The final histograms (orange and yellow bars in figure 3(a)) are now very close to P f (n). The figure of merit of the reconstruction, defined as the standard deviation of P i max (n) with respect to P f (n), is (29) We find that σ = 0.01 for N D = 3 and 0.001 for N D = 4 and 6. For comparison the initial guess P 0 (n) corresponds to σ 0 = 0.025. These numerical results confirm that N D = 4 atoms are necessary and sufficient to reconstruct faithfully an arbitrary photon number distribution with n m = 7. In addition, they indicate the order of magnitude of the number of required iterations. Increasing N D above N min D does not improve the reconstruction fidelity but decreases the number of required iteration steps.

Simulating the reconstruction of a time-varying P(n) under realistic experimental conditions
We now take into account in the simulations the realistic parameters of our cavity QED experiment described in section 2. We recall that our goal is the reconstruction of a decaying photon number distribution with a temporal resolution much better than T c /n m . We thus now introduce in the simulation field relaxation between atomic samples by numerical integration of the field master equation. This influences the simulated detection results. However, the reconstruction part in the simulation ignores relaxation assuming that the average measurement time for detecting N min D atoms is short enough with respect to the field relaxation time. The Monte Carlo simulation now also accounts for the Poisson distribution of the atoms in each sample. In addition, it replaces the ideal single-atom detection probability π( j|n, φ) of equation (2) by with c s = 0.76 being the limited contrast of our Ramsey interferometer [11]. This replacement modifies accordingly the expressions ofÊ φ,J and hence the MaxLik iteration operatorR(ρ). A single realization of a simulated experiment involves 500 successive samples crossing C over a 41 ms time interval. We have investigated the evolution of two initial fields: a nearly pure n = 4 Fock state (green histogram in figure 4(a)) corresponding to the field state prepared in the quantum feedback experiment of [11] and a coherent state with an average photon number of 3 (green histogram in figure 4(c)). As in a real typical experiment, we resum the sequence 4000 times. We reconstruct the photon number distribution P(n, t) by applying the iterative procedure to the subset of atoms detected in N S = 20 samples crossing C in a time interval centered at t. The corresponding time window is 1.6 ms, which corresponds to N D 4 detected atoms on the average. It defines the time resolution of the reconstruction. Note that the atoms detected before the time window at t do not influence the average photon number distribution at this time because of the QND nature of the atom-cavity interaction. Hence, we exploit information provided by all the atoms crossing C in each realization and not only that contained in a single time window. The full reconstruction process is repeated 100 times in order to obtain the statistical error bars. The time origin t = 0 corresponds to the detection of the first atomic sample used for reconstruction.
The thin lines in figures 4(b) and (d) display one of the reconstructed P(n, t)s corresponding to the initial states described by the green histograms in figures 4(a) and (c), respectively. The statistical noise in the signals results from the finite number of realizations used for the reconstruction. The thick solid lines are theoretical fits to the reconstructions. The only fit parameters are the initial photon number probabilities denoted by P ini (n) from which we calculate the theoretical values of P(n, t), knowing T c , n th and using the field density matrix master equation [17]. The fit is adjusted to the reconstructed signals up to t = 16 ms. Note that this fitting procedure directly gives the initial photon number distribution. It uses the a priori knowledge of cavity damping time and of the residual blackbody radiation for eliminating systematic errors due to the effect of field damping during the first N S atom sample measurement. The averages of 100 independent reconstructions of P ini (n) are shown by red histograms in figures 4(a) and (c) for the two initial states. The standard deviation of the reconstruction with respect to the 'real' state is σ = 0.004 for the Fock state and σ = 0.002 for the coherent one, while the initial flat distribution corresponds to σ 0 = 0.09 and 0.03, respectively. In the case of the Fock state, the reconstructed population P ini (4) = 80 ± 3% is, within the error bar, equal to the expected one (82%). In addition to yielding a good estimation of the initial state, these reconstructions faithfully reproduce the expected field evolution even for times outside the 16 ms window used for the fit. This makes us confident that our reconstruction procedure is able to reproduce with good fidelity and high time resolution the evolving photon number distribution of a relaxing field in a cavity. The field to be reconstructed is close to a Fock state (photon number distribution P f (n) shown by the green histogram in panel (a)). The red histogram in panel (a) is obtained by averaging 100 independent reconstructions of P ini (n). Error bars represent the standard deviation of these reconstructions. In (b) the thin lines represent one of the 100 reconstructed P(n, t) versus time. The solid lines present the theoretical fit from which the initial distribution P ini (n) is extracted. (c) and (d) The same as (a) and (b) for an initial coherent state with three photons on average. For the sake of clarity, only P(n = 2, t) and P(n = 4, t) are shown in the reconstructed signal.

Reconstruction of real experimental data
In [11,12], quantum feedback protocols have been used to prepare on demand photon number states. The field is probed by QND measurements carried out with atoms allowing a control computer to estimate in real time the field state and to decide which correction action must be taken in order to bring this state closer to the target. The procedure is repeated in a succession of loops. When the controller estimates that the target is reached with a threshold fidelity, it stops acting on the field. Repeating this operation 4000 times, we prepare a statistical ensemble on which we apply the MaxLik reconstruction procedure with n m = 7 as an independent estimation of the prepared field state.   [11] when the fidelity with respect to the target state n = 4 reaches the 0.8 threshold. This initial state was chosen for the simulation of section 6, so that the green histograms in figures 4(a) and 5(a) are identical. Figure 5(b) shows the evolution of P(n, t) during the 100 ms following the interruption of the feedback loop. The reconstruction method is the same as the one used in section 6 with N D 4 atoms in each time window. The thin lines describe the decay of the initial approximate four-photon Fock state after the termination of the feedback loop. The solid lines are fits realized as described in section 6. The real data are in good agreement with theoretical predictions over the total 100 ms measurement time. The fit parameters P ini (n) are displayed in figure 5(a) (red histogram). The P ini (4) probability is 0.74, slightly smaller than the one (0.82) estimated by the controller. This indicates that the controller overestimates the fidelity of the state it prepares, its task being much more complex than merely reconstructing the photon number distribution. This explains the small discrepancy between the green and the red histograms in figure 4(a). However there is great similarity between the reconstructions of the experimental and simulated data (the same statistical noise and time evolution in figures 4(b) and 5(b)). Note that in the reconstructed distribution P ini (7) 0 ( figure 5(a)). The vanishing values of P ini (0) and P ini (1) also prove that the eight and nine photon Fock states are not populated (with the ϕ 0 = 2π/(n m + 1) = π/4 setting, the atomic POVMs do not distinguish n from n + n m ). This justifies our choice of the maximum photon number n m = 7. Non-zero populations over the whole 0-n m interval may indicate an improper choice of n m . A similar time-resolved reconstruction procedure was applied to study the evolution of fields prepared in Fock states by random projective measurements [15]. The detailed analysis presented here validates the procedure that was only briefly outlined in [15].

Conclusion
We have presented here a powerful MaxLik method for the statistical reconstruction of fields evolving in a cavity. It uses a stream of non-resonant atoms as a QND dispersive probe. We have shown that a minimal sequence of N D n m /2 atoms around a given time in each realization is necessary and sufficient to reconstruct a snapshot of the field photon number distribution, provided it does not contain more than n m photons. Each realization contributes to the reconstruction of the field at all times since information can be extracted from a long atomic sequence divided into elementary N D atom bins. Due to the QND character of our measurement scheme, data coming from different bins contribute independently to the statistical reconstruction. In this way, information is acquired more rapidly than in procedures involving destructive measurements, which are limited to one measurement in a given time bin for each realization.
The presented method has been successfully tested on numerically simulated as well as on experimentally measured data. For an ideal atomic beam delivering one atom every T a , the time resolution of this method of state reconstruction is T a n m /2, which should be smaller than T c /n m the lifetime of the Fock state with the maximum photon number in the field. The method can thus reconstruct the evolution of arbitrary photon number distribution provided n m < √ 2T c /T a . If we had a deterministic source of atoms and a perfect detector, the limit in our experiment would be n m 40, this figure being limited to 10 by the randomness of the present atomic source and the limited detection efficiency. Improving on these factors and on the cavity damping time should allow us to time resolve the evolution of fields involving up to a few tens of photons.
Our method can be generalized to the full reconstruction of field density operators [8] by applying controlled field displacements before measurement of the N D -atom POVMs. The minimum number of detected atoms in a time bin is still given by N min D = I [(n m + 1)/2], where n m is now the number of Fock states populated in the displaced field. Note finally that our reconstruction method, based on non-resonant interaction with the cavity, may be of interest in the context of circuit QED with three-dimensional cavities [20]. In these experiments, an increase of the coherence time of the qubit is obtained at the expense of reduced control, which makes it difficult to implement the reconstruction method used in [9].