1. Introduction

The types of photonic states typically encountered in quantum information science (QIS)—including quantum networking, computing, and sensing—can be divided into two broad categories, depending on the fundamental unit of information employed. In qubit (or more generally, qu$d$it) encoding, logical states are defined by single photons which occupy $d$ electromagnetic modes; in this case, photonic QIS protocols are typically defined for a specific numbers of particles [1,2]. In qumode encoding, however, the electromagnetic modes themselves assume fundamental priority, with information encoded in, e.g., field variables such as the quadratures or the number of photons [3–5]. Now, because qubit-encoded optical states are typically available in a probabilistic fashion only—due to spontaneous generation processes [6], post-selection by design [1,2], or simply from loss—the qumode, or continuous-variable (CV), formalism is arguably the more general of the two viewpoints, for it accounts for variable numbers of photons and subsumes discrete-variable (DV) qubit/qudit encoding as a special case.

This generality yields a significant increase in complexity. For example, whereas a single qubit is fully described within a two-dimensional Hilbert space, an infinite-dimensional Hilbert space is required to specify a single qumode. Although this dimensionality can be truncated in practice with reasonable assumptions about maximum photon number, fully general quantum state tomography (QST) of a qumode [7] remains a considerably more computationally intensive affair than QST of a qubit [8]. Based on a collection of homodyne or heterodyne measurements with a local oscillator, numerical techniques for qumode QST have progressed from linear inversion with the Radon transform [9,10] to maximum likelihood estimation (MLE) [11,12], compressive sensing [13,14], and, very recently, neural networks [15].

In any method, one desires an estimate of the underlying quantum state which is as close as possible to the unknown ground truth. Bayesian techniques prove uniquely situated to reach this end [16–18]. Given a generic dataset ${\boldsymbol {\mathcal {D}}}$ and to-be-determined quantum state $\rho$, Bayes’ theorem expresses the possible states in terms of a probability distribution $\mathrm {Pr}(\rho |{\boldsymbol {\mathcal {D}}}) = \mathrm {Pr}({\boldsymbol {\mathcal {D}}}|\rho )\mathrm {Pr}(\rho )/\mathrm {Pr}({\boldsymbol {\mathcal {D}}})$, which (i) provides automatic uncertainty quantification for any collection of measurements or outcomes and (ii) returns the optimal estimate of any quantity of interest: the Bayesian mean estimator formed as the mean over $\mathrm {Pr}(\rho |{\boldsymbol {\mathcal {D}}})$ minimizes the squared error averaged over all states and outcomes [16]. Unfortunately, the high-dimensional integrals involved with full Bayesian inference only exacerbate the computational challenges of qumode QST. And while Markov chain Monte Carlo (MCMC) algorithms have facilitated Bayesian tomography of a variety of discrete qubit/qudit systems [19–26] and continuous qumode states that are Gaussian [27], no method for Bayesian inference of arbitrary CV states has been proposed.

In this work, we introduce and demonstrate a complete workflow for Bayesian QST of generic CV states. Capable of handling both homodyne and heterodyne measurements, our method combines physical models developed in previous CV estimation approaches with advanced MCMC algorithms [28] for full Bayesian tomography with state-of-the-art computational efficiency. After detailing our method, we demonstrate inference on several experimental datasets, including newly measured coherent and thermal state examples and legacy cat state data from Ref. [29]. As the first Bayesian method for CV QST that does not assume Gaussianity in the model, our workflow closes the existing gap between CV and DV systems in Bayesian QST and provides new opportunities for informationally efficient tomography in state-of-the-art non-Gaussian photonic QIS.

2. Model

Figure 1 depicts the envisioned tomographic scenario. An arbitrary single-mode quantum state with density matrix $\rho$ first undergoes loss described by the efficiency parameter $\eta \in [0,1]$. Then it is either sent to a homodyne setup (top), where it is mixed with a local oscillator aligned to phase $\theta$ and measured with a balanced detector, or to a heterodyne setup (bottom), which divides the signal with a 50/50 beamsplitter and performs simultaneous homodyne detection of the outputs, where the local oscillator is aligned to $\theta$ or $\theta + \pi /2$. Because loss commutes with either homodyne or heterodyne detection [30], non-unit detector efficiency can be absorbed into $\eta$ and the detection itself viewed as ideal; in the process, we do invoke the assumption that all detectors possess equal efficiencies—the desired design condition in practice.

 

Fig. 1. Bayesian QST workflow. An arbitrary single-mode quantum state with density matrix $\rho$ first undergoes loss, then is sent either to a homodyne setup (top), or to a heterodyne setup (bottom). Results ${\boldsymbol {\mathcal {D}}}$ are then fed into an inference algorithm based on Bayes’ theorem, which returns the posterior distribution over possible density matrices, given the observations, $\mathrm {Pr}(\rho |{\boldsymbol {\mathcal {D}}})$.

Download Full Size | PDF

While the ground truth state exists in an infinite-dimensional space, truncation is required for numerical tractability. To this end, we choose to express $\rho$ in the Fock basis $|{n}\rangle$, so that such truncation can be realized directly by imposing a photon-number cutoff $n_c$, i.e.,

This state is then mixed with a local oscillator (LO) with relative phase $\theta$ and detected. In the homodyne case, the quadrature measurements $x$ follow the probability density function [7]

In the tomographic context, these measurements are repeated $K$ times, leaving dataset ${\boldsymbol {\mathcal {D}}} =\{(\theta _k,x_k)\}$ (homodyne) or ${\boldsymbol {\mathcal {D}}} =\{(\theta _k,x_k,p_k)\}$ (heterodyne) where $k=1,2,\ldots,K$. Each phase value $\theta _k$ is allowed to vary but is assumed known from independent calibration. The likelihood $L_{\boldsymbol {\mathcal {D}}}(\rho ) \propto \mathrm {Pr}({\boldsymbol {\mathcal {D}}}|\rho )$ follows as the product of individual outcomes:

Up to this point, no aspects of this model are unique to Bayesian inference, but instead are common to any estimation procedure built on the likelihood. Indeed, MLE proceeds precisely by finding the physical density matrix $\rho$ which maximizes $L_{\boldsymbol {\mathcal {D}}}(\rho )$ above [Eq. (6) or Eq. (7)]. In the Bayesian paradigm, however, the likelihood is accompanied by a prior distribution as well, often chosen to be as uniform as possible to give reasonable weight to all possible states in the Hilbert space. Examples of fiducial density matrix distributions, studied extensively in DV encoding, include Hilbert–Schmidt and Bures, the former equal to the distribution of $D\times D$ density matrices found by tracing out random pure $D^{2}$-dimensional states [35], and the latter noteworthy as the only monotone distribution that is both Fisher and Fubini–Study adjusted [36].

To our knowledge, no comparable analyses of density matrix distributions exist for the CV case, leaving wide open the question of fiducial prior selection. Since the photon cutoff $n_c$ [cf. Equation (1)] should exceed the maximum photon number of all states of interest, it certainly would seem reasonable to consider a prior which preferentially weights lower photon numbers over higher ones—in contrast to uniform DV priors that place all logical basis states on equal footing. Nonetheless, given the absence of specific alternatives in the literature, we elect to implement a Bures prior distribution in the truncated $(D=n_c+1)$-dimensional space for our inference examples here. In applying uniform weight to all photon numbers up to $n_c$, Bures can be viewed as a highly cautious distribution that will permit low-uncertainty estimates only when justified by the data collected.

To construct density matrices drawn from the Bures distribution, we first define a $D^{2}$-dimensional vector of complex parameters $\mathbf {z}=(z_1,\ldots,z_{2D^{2}})$; the first $D^{2}$ populate a $D\times D$ matrix $G$, while the second $D^{2}$ parameters are used to form another complex matrix that is converted to a $D\times D$ unitary $U$ through the Mezzadri algorithm [37]. The density matrix constructed according to

The posterior density can therefore be written as

Obtaining the samples $\{\mathbf {z}^{(r)}\}$ represents the core computational bottleneck of Bayesian methods and has motivated decades of research in MCMC techniques [16,17]. Recently, we applied a particularly efficient MCMC method, known “preconditioned Crank–Nicolson” (pCN) [28], to Bayesian QST [23]. Designed to eliminate the step size/acceptance rate tradeoff intrinsic to random walk MCMC, pCN was found in our tomographic examples to obtain significant speedups compared to alternative algorithms [23]. This pCN-based Bayesian QST workflow has since been applied to several problems with a Bures prior distribution [25,26]. We adopt the same sampling procedure here, with the only substantive difference being the likelihood computed: starting with a parameter vector $\mathbf {z}^{(r)}$, we first compute $\rho$ according to Eq. (8); then the lossy $\tilde {\rho }$ follows via Eq. (3); and finally, the likelihood is evaluated according to Eqs. (4, 6) for the homodyne case, or Eqs. (5, 7) for heterodyne.

3. Experimental demonstration

3.1 Setup

Our setup for producing experimental tomographic datasets is summarized in Fig. 2. To generate a local oscillator (LO), we use a continuous-wave (CW) laser at 1548.8 nm (Pure Photonics PPCL550), then send it through a phase modulator (EOSpace PM-0K5-10-PFA-PFA-DC), which is driven with a voltage ramp (8.3 Vpp, 5 kHz) from a function generator (Stanford Research Systems DS345) producing a full $2\pi$ phase swing, then to the LO input of the heterodyne detector (Optoplex RX-KC0100C821AC), which is a 90$^{\circ }$ optical hybrid with outputs connected to a pair of amplified balanced detectors. The two electrical outputs of the detector are 50 $\Omega$ terminated into a digital sampling oscilloscope (Keysight MSOX4104A) that we use to acquire a total of $K=7998$ measurement points $(\theta _k,x_k,p_k)$, where each quadrature value is the average of four consecutive 2.5 GHz samples on our oscilloscope, repeated at an interval of 100 ns, for all examples.

 

Fig. 2. Heterodyne tomography setup for coherent and thermal states. For coherent states, the continuous-wave (CW) laser is split and used to make attenuated coherent states and a local oscillator (LO) for the tomography. For thermal states, an erbium-doped fiber amplifier (EDFA) is filtered and attenuated to create thermal states of low average photon number. The CW laser is used as the LO in the tomography. Definitions: VA = variable attenuator, VBS = variable beamsplitter, PC = polarization controller, Atten = fixed attenuator(s), POL = polarizer, PM = phase modulator, FG = function generator, and PF = programmable filter.

Download Full Size | PDF

To measure coherent test states (Fig. 2), we split about 1% of the CW laser power before the phase modulator using a variable fiber beamsplitter (Newport F-CPL-1550-N-FA) and send it through a polarization controller. It is subsequently passed through a variable optical delay line (General Photonics VDL-001-35-33-FC/APC-SS) and delay fiber to approximately match the path length to the LO, then through several fixed attenuators, and finally to the signal input of the heterodyne detector. Signal and LO path lengths are matched to within about 1 ns, which is much shorter than the approximately 0.1 ms coherence time of the laser. From the raw measurements, we calculate $g^{2}(0)\approx 1$ using the analysis method in Ref. [30], consistent with the production of high-quality coherent states.

To produce thermal states (inset in Fig. 2), we use the amplified spontaneous emission from an erbium-doped fiber amplifier (Pritel SCG-40), filter it to 1548.8$\pm$0.4 nm (Finisar Waveshaper 1000A), and send it to the signal input of the heterodyne detector. The analysis method of Ref. [30] returns $g^{2}(0)\approx 2$ for this case, as expected for thermal states. See Appendix A for descriptions of our calibration and data processing procedures. Finally, all Bayesian inference results below consist of $R=1024$ total samples, which were themselves obtained by thinning an MCMC chain of total length $RT$, where $T$ was selected empirically for convergence in the mean and standard deviation of the inferred fidelity.

3.2 Heterodyne inference

We test our Bayesian method on six separate states of varying type and amplitude; the insets of Fig. 5 plot the raw quadrature data in phase space $(x_k,p_k)$, with color signifying the respective LO phase $\theta _k$. The quadrature data of Fig. 5(a) show clear variation with LO phase, as expected for coherent states, while the data of Fig. 5(b) appear rotationally symmetric and remain centered at the origin, characteristic of thermal states. In order to provide a benchmark to compare against subsequent Bayesian results, we compute an expected coherent state parameter $\alpha _0$ or thermal state mean photon number $\mu _0$ by correcting for the LO phase and performing averages over the raw data: $\alpha _0=\langle {x_k \cos \theta _k - p_k\sin \theta _k}\rangle + \mathrm {i}\langle {x_k\sin \theta _k + p_k\cos \theta _k}\rangle$ or $\mu _0 = \langle {x_k^{2} + p_k^{2}}\rangle - 1$. These values are then used to define the expected ground truth states for fidelity calculations. In defining the states in this fashion, we are effectively taking $\eta =1$ in the model, absorbing any loss into the state itself, which has little impact in these cases since coherent and thermal states preserve their statistics under loss [30].

We took these data and produced inferred density matrices with our Bayesian QST workflow, using an LO power of 12 mW, an MCMC thinning value of $T=2^{14}$, and taking our photon cutoff at $n_c=20$; for the expected states examined, this selection ensures a truncation error $\epsilon = 1-\sum _{n=0}^{n_c} \langle {n|\rho _g|n}\rangle$ of at most 0.003 (for the $\mu _0=3.03$ thermal state); the truncation error for all coherent states is less than $10^{-4}$. In Fig. 3, we show an example of how the elements of the Bayesian-mean estimate $\rho _B$ progress as more measurements $K$ are included in the likelihood, for the attenuated laser dataset with $\alpha _0 = -2.78 -0.54\mathrm {i}$ ($|\alpha _0|^{2}=7.97$). At a single measurement $K=1$, virtually no information has been gained about the ground truth state, and $\rho _B$ averages to the maximally mixed state in the 21-dimensional space, as needed for the Bures prior. For $K=400$, the result clearly resembles a coherent state, with a residual diagonal spike from the prior, and by $K=1600$ the estimate matches the expected ground truth with fidelity $F=0.86\pm 0.01$, improving to a fidelity of $F=0.958\pm 0.004$ at $K=7998$, which we expect would continue to increase with even more measurements.

 

Fig. 3. Scaling with the number of measurements for experimental heterodyne tomography. Retrieved Bayesian mean density matrices and fidelities for the case $|\alpha _0|^{2}=7.97$, for several values of the number of measurements $K$.

Download Full Size | PDF

The scaling of fidelity with the number of measurements for all three attenuated laser examples is plotted in Fig. 4(a), with the final results (for $K=7998$) in Fig. 5(a), along with the complete raw datasets and parameters for the expected ideal states as insets. Color denotes phase, either for the LO setting $\theta _k$ or for the estimated density matrix element $\rho _{mn}$. Nonzero phase values for the $\rho _{mn}$ elements are a consequence of the measured state’s angular position in phase space with respect to the LO; although fluctuating randomly due to environmental perturbations, it effectively remains fixed within the 0.8 ms duration of data collection for each measurement set. Final fidelities approach unity in all cases: $F=0.958\pm 0.003$ ($\alpha _0 = 1.14-0.45\mathrm {i}$, $|\alpha _0|^{2}=1.51$), $F=0.956\pm 0.004$ ($\alpha _0 = 0.24-1.76\mathrm {i}$, $|\alpha _0|^{2}=3.14$), and $F=0.958\pm 0.004$ ($\alpha _0 = -2.78-0.54\mathrm {i}$, $|\alpha _0|^{2}=7.97$).

 

Fig. 4. Fidelity scaling with heterodyne detection. (a) Fidelity of inferred states with respect to expected ground truth for an attenuated laser (coherent state). (b) Inference fidelity for amplified spontaneous emission (thermal state).

Download Full Size | PDF

 

Fig. 5. Bayesian QST with heterodyne detection. (a) Bayesian mean density matrices estimated from all 7998 measurements of each attenuated laser dataset. (b) Bayesian mean density matrices using all 7998 measurement results for each tested thermal state. Insets show the raw data points $(x_k,p_k)$ in measurement phase space, with color indicating the LO phase $\theta _k$. The values $|\alpha _0|^{2}$ ($\mu _0$) denote the mean photon number of the expected ideal coherent (thermal) state.

Download Full Size | PDF

The results for amplified spontaneous emission appear in Fig. 4(b) and Fig. 5(b). The diagonal elements of the final estimates mirror the exponential decay of ideal thermal states, although the appreciable off-diagonal entries reveal deviations that are quantified by reduced fidelities compared to the coherent state tests: $F=0.90\pm 0.03$ ($\mu _0 = 0.68$), $F=0.82\pm 0.04$ ($\mu _0 = 1.49$), and $F=0.77\pm 0.03$ ($\mu _0 = 3.03$). While we expect the thermal states to reach comparably high fidelities with additional measurements, we have noticed a similar reduction in thermal state tomographic efficiency in simulation (see Appendix B).

3.3 Comparing homodyne and heterodyne

To test the ability of our method for inference on both homodyne and heterodyne measurement schemes, as well as highlight the impact of $\eta$ selection, we can obtain homodyne data from the setup in Fig. 2 (for a 6 mW homodyne LO) by considering only the $x$-quadrature results; the traced-out $p$ measurement can then just be viewed as an additional 3 dB loss on the state to be characterized. Therefore, heterodyne inference with $\eta =1$ and homodyne inference with $\eta =0.5$ should match the same ground truth in our setup. Nonetheless, since this arrangement penalizes homodyne with additional loss, we also obtain a third dataset in which we double both the LO power (from 6 mW to 12 mW) and the signal power to counteract the internal beamsplitter loss in the $90^{\circ }$ optical hybrid. For coherent and thermal states, where loss only reduces average photon number [30], the traced-out homodyne data gathered under this “double-power” condition reflects a state of identical statistics and average photon number to the other two cases, by taking $\eta =1$ in the likelihood. And because the vacuum noise far exceeds the electronics floor at both LO powers (see Appendix A)—so that detector noise can be neglected—we can directly compare tomography for three separate conditions in which the states to be inferred are matched, thus minimizing the impact of any state-dependent effects on the results. (For squeezed states or non-Gaussian states affected by loss in nontrivial ways, we acknowledge this comparison would not be valid.)

Results for both an attenuated laser and amplified spontaneous emission appear in Fig. 6, where the average photon number in all cases is $\sim$1.5. A thinning of $T=2^{14}$ and cutoff of $n_c=20$ are again used. Unsurprisingly, the lowest fidelities are obtained by the homodyne case with $\eta =0.5$: compared to heterodyne, its data are a subset that ignores the information from the entire $p$ quadrature; compared to the homodyne case with $\eta =1$, the amount of quadrature information is the same, but embedded within a higher proportion of vacuum noise. Yet the two cases with $\eta =1$ prove more subtle, showing nearly identical efficiency on the coherent state example, but an appreciable edge in favor of homodyne for the thermal state. The relative merits of homodyne and heterodyne have been explored in several contexts, including tomography of Gaussian states [39,40] and the channel capacities of Gaussian receivers [34]. Yet although our examples here are also Gaussian, Gaussianity is not assumed in the Bayesian prior, which considers all physically allowed states. Therefore it is unclear how previous findings which assume Gaussianity a priori would apply to our tests. In any regard, Fig. 6 serves to reveal the ability of our Bayesian approach to handle multiple measurement configurations and efficiency parameters, returning a justifiable state estimate under the prior distribution and available data—whatever the data may be.

 

Fig. 6. Comparing heterodyne (HET) and homodyne (HOM) tomography at multiple efficiencies $\eta$. (a) Attenuated laser results. (b) Amplified spontaneous emission. All examples have a mean photon number of $\sim$1.5. Homodyne datasets are obtained by considering only the $x$ quadrature data in Fig. 2, with signal and LO power either equal to the heterodyne case (for HOM, $\eta =0.5$) or double (for HOM, $\eta =1$).

Download Full Size | PDF

4. Cat state inference

Our Bayesian QST workflow is specifically designed for tomography of arbitrary CV states, including non-Gaussian states of relevance to QIS. While the experimental production of such states is beyond the current capabilities of our setup in Fig. 2, in this section we specifi cally test our approach on previous tomographic data from the non-Gaussian experiment described by Gerrits et al. [29], in which cat states were produced via photon subtraction of squeezed vacuum. In that work, MLE was employed on the homodyne data to recover the state produced, followed by parametric bootstrap resampling to estimate error bars. We note that Fisher information could also be used for uncertainty quantification with MLE-based QST [41]; in any case, the MLE and Bayesian findings are expected to agree perfectly in the asymptotic limit of infinite observations. Yet no clear relationship exists in the finite-measurement regime explored here, making it constructive to compare, so in this section the predictions of our Bayesian method to those of MLE on identical sets of quadrature samples.

In Fig. 7, we show the results from our Bayesian QST workflow (with $T=2^{13}$ and $n_c=10$) applied to data in Ref. [29] which used a transition edge sensor (TES) to implement three-photon subtraction and produce an odd cat state of the form $|{C_{-}(\alpha )}\rangle \propto |{\alpha }\rangle - |{-\alpha }\rangle$. To aid visual analysis, we plot the Wigner function of the Bayesian mean $\rho _B$, rather than the Fock basis elements as before; Fig. 7 can thus be compared directly against the MLE Wigner function in Fig. 4 of Ref. [29]. Adopting the same optimization approach and error bar percentiles as in Gerrits et al., we calculate that the nearest ideal cat state has $|\alpha |=1.64^{+0.09}_{-0.08}$, and the fidelity between that cat state and the state reconstructed from the data is $F=0.49^{+0.09}_{-0.09}$, in good agreement with the previous MLE findings $|\alpha |=1.76^{+0.02}_{-0.19}$ and $F=0.59^{+0.04}_{-0.14}$. The relatively large error bars result from the limited number of quadrature measurement values ($K=1087$), caused by the low probability of subtracting three photons. Bayesian-estimated Wigner functions for all other datasets from Ref. [29], and a comparison of our findings with those from MLE, can be found in Appendix C. For the purposes of the present investigation, however, the key point of our results is to highlight the successful application of Bayesian inference in reconstructing a non-Gaussian ground truth quantum state—a first for Bayesian QST of CV systems.

 

Fig. 7. Wigner function of the Bayesian mean density matrix $\rho _B$ estimated from homodyne data on a three-photon-subtracted squeezed vacuum state from Ref. [29]. The nearest ideal cat state has $|\alpha |=1.64^{+0.09}_{-0.08}$ with fidelity $F=0.49^{+0.09}_{-0.09}$.

Download Full Size | PDF

5. Discussion

In our experimental examples so far, we have observed appreciably lower fidelities for the inference of thermal states compared to coherent states (Figs. 4–6), which is corroborated in simulation as well [see Fig. 9(b) in Appendix B]. We hypothesize that this effect derives from the mixedness of thermal states: mixed states require additional parameters for their specification compared to pure states, suggesting in general the need for more measurements to reach a given error level. Thus, Bayesian inference should attain higher fidelity with respect to a ground truth pure state than to a thermal one at the same number of measurements, although further investigation into this effect would be useful.

Moving forward, as with any MCMC method, computational efficiency represents the strongest impediment toward the continued scaling to higher dimensions and multiple quantum modes. For reference, the individual chains behind the results in Figs. 4–7 required up to a maximum of $\sim$40 hours to complete on our desktop computer (for $n_c=20$, $K=7998$, and $RT=2^{24}$ total steps)—a number which will increase both with photon cutoff $n_c$ and number of measurements $K$. Speedup through parallelization is generally hampered by the intrinsically serial nature of Markov chains, but recent developments in parallel approaches [42–45] may present a path to estimation of vastly more complex quantum systems than hitherto possible in Bayesian QST.

Unexpectedly, the development and testing of our new method brought us into contact with several fascinating general problems in CV quantum state characterization that extend well beyond the Bayesian focus here. The fact that the relative tomographic efficiency of heterodyne and homodyne measurement varies with the ground truth signal (cf. Figure 6) aligns roughly with previous observations on Gaussian states [34,39]. Yet no systematic study comparing the efficiency of homodyne and heterodyne measurements has been performed for arbitrary CV states. From a theoretical perspective, both homodyne and heterodyne detection are informationally complete, in that they measure the state in complete representations—namely, the Wigner [7] and Husimi-$Q$ [46] functions, respectively. But the number of measurements required for a given level of accuracy need not be the same. In the future, therefore, it would be interesting to formally investigate the efficiency of homodyne and heterodyne detection on arbitrary quantum states, to elucidate what features of some ground truth state $\rho _g$ determine which method should be selected.

Another open question in quantum theory that surfaced in our analysis centers on distributions of random density matrices for CV quantum states. The Bures prior we have selected enjoys several commending theoretical properties, as the lone monotone metric which is both Fisher and Fubini–Study adjusted—i.e., it reduces to desired distributions in the classical and pure-state limits [36]. However, as a DV-based distribution, Bures depends on the specific Hilbert space dimension which, in the CV tomographic case, is set somewhat arbitrarily by the photon cutoff $n_c$. Accordingly, it would prove theoretically satisfying to develop a fiducial CV density matrix distribution with a less abrupt stepwise $n_c$ dependence. This objective would seem to bear profitable overlap with other advanced CV tomographic techniques, such as the neural network approach of Ref. [15]; there the authors utilized a restricted Boltzmann machine ansatz, acknowledging the open question of whether practically interesting quantum states may be eliminated by this construction. Thus, we suspect that a more detailed understanding of random CV quantum state distributions could provide valuable insight into this problem as well.

Nevertheless, we emphasize that the existence of such interesting open problems in no way diminishes the immediate practical value of Bayesian inference for CV QST. Indeed, one of the most attractive features of Bayes’ theorem is its ability to incorporate any assumptions an experimenter wishes to bring to a problem, via an explicit prior distribution, and then determine the impact of both these assumptions and subsequent observations for estimation of any quantity of interest. Consequently, there is no need for a prior to align with some archetypal distributions; instead, it need only reflect whatever prior knowledge is imposed by the particular experimenter. And regardless of how this prior emerges, the Bayesian paradigm ensures that all underlying assumptions are enumerated and exposed for discussion and analysis—a situation of immense value as the field of QST continues to probe new regimes of state complexity.

6. Conclusion

In this work, we have introduced a workflow for Bayesian QST of generic CV states. This general method is capable of handling homodyne and heterodyne measurements and uses advanced MCMC algorithms for full Bayesian tomography with state-of-the-art computational efficiency. We experimentally prepared multiple types of Gaussian states, showing the inferred density matrices obtained from Bayesian estimation. Moreover, we used legacy measurement data of cat states to verify our workflow’s ability to reconstruct non-Gaussian states as well. This powerful computational analysis tool could be helpful in the precise characterization of up-and-coming CV quantum computing hardware [47], CV quantum networking [48], and CV quantum sensors [49], where it is crucial to accurately characterize the states produced.

Appendix A: experimental procedures

For polarization alignment, we found it necessary to minimize the cross-talk between the fast- and slow-axes of the polarization-maintaining-fiber at the couplings before and after the phase modulator, using polarizers and polarization controllers ; if the LO polarization is not aligned perfectly to the axis of the phase modulator, it is possible for the phase modulator to rotate the polarization of the LO, which undesirably makes the shot noise oscillate with time. To calibrate the shot noise, for a given LO power, we record two million samples (spanning 0.8 ms at our 2.5 GHz sampling rate) for each output of the heterodyne detector. At a spacing of 250 samples (100 ns) of each output channel, adjusted by the known electronic delay mismatch (separately calibrated with bright pulses), we average groups of four adjacent samples and save the average values as single quadrature points. We perform this procedure for the whole data file, which produces 7998 measurements, from which we calculate the mean, $SN_{\text {m}}$, and shot-noise variance, $SN_{\text {var}}$, for each detector channel. We calibrate the shot noise for a range of relevant optical powers [see Fig. 8(a)]. We find a linear operating range from 1–15 mW, as measured at the input to the heterodyne detector, and primarily operate at $\sim$12 mW to be far from the electronics noise. We also measured the electronics noise and plot the ratio of the shot-noise variance and the electronics noise variance in Fig. 8(b), which at our primary operating point is about 17.

 

Fig. 8. Shot-noise calibration. Channel 1 data (fit) correspond to the blue circle (line). Similarly, channel 2 data (fit) correspond to the orange square (line). (a) Shot-noise variance versus LO power. From the fit: R$^{2}_\text {CH1}$=0.9998 and R$^{2}_\text {CH2}$=0.9997. (b) Ratio of shot-noise variance to electronics noise variance versus LO power. We normally operate at about 12 mW to balance linearity and electronics noise.

Download Full Size | PDF

To calibrate the phase modulator’s applied voltage for a total phase swing of $2\pi$, we send in a bright coherent input signal and adjust the amplitude of a 5 kHz sawtooth voltage signal to achieve a continuous sinewave output on the oscilloscope. Any peak-to-peak amplitude not corresponding to a multiple of $2\pi$ will produce a discontinuous jump at the falling edge of each sawtooth, and when exactly $2\pi$ is achieved, a single sinewave period will be traversed during a single ramp. A square-wave synchronization signal from the function generator is also sent to the oscilloscope and saved to provide a marker for the beginning and end of each ramp. Assuming linear modulator operation, we then know the phase change between each sample.

When there is an input signal, we convert our collected data from voltage into dimensionless quadrature units, using the shot-noise calibration. As before, we average adjacent groups of four oscilloscope samples, spaced every 250 samples, to provide individual voltage points $V_{\text {avg}}[k]$. For our total 0.8 ms record length, we therefore complete a full $2\pi$ phase sweep four times. We then normalize these values such that a vacuum input has mean zero ($\langle {x}\rangle _\text {vac}=\langle {p}\rangle _\text {vac}=0$) and variance 1/2 ($\Delta x^{2}_\text {vac} = \Delta p^{2}_\text {vac} = 1/2$). Thus a single pair of quadratures $(x_k,p_k)$ is given by

(11)$$x_k = (V_{X,\text{avg}}[k]-SN_{\text{X,m}})\sqrt{\frac{1/2}{SN_{\text{X,var}}}}$$

(12)$$p_k = (V_{P,\text{avg}}[k]-SN_{\text{P,m}})\sqrt{\frac{1/2}{SN_{\text{P,var}}}}$$

for all $k\in \{1,2,\ldots,K\}$. These normalized points can then be input into our tomographic method.

Appendix B: simulations

For a further test on our workflow, we produced simulated data for several representative classes of states and perform Bayesian inference as a function of measurement number. As ground truth $\rho _g$ ($\rho _g=|{\psi _g}\rangle \langle {\psi _g}|$ if pure), we consider coherent states

(13)$$|{\psi_g}\rangle = e^{-|\alpha_0|^{2}/2} \sum_n \frac{\alpha_0^{n}}{\sqrt{n!}}|{n}\rangle,$$

thermal states

(14)$$\rho_g = \sum_n \frac{\mu^{n}}{(1+\mu)^{n+1}} |{n}\rangle\langle{n}|,$$

squeezed vacuum

(15)$$|{\psi_g}\rangle = \frac{1}{\sqrt{\cosh r}} \sum_n \sqrt{\frac{(2n)!}{2^{n} n!}} \left(\frac{1}{2}\tanh r \right)^{n} |{2n}\rangle,$$

and Fock states

(16)$$|{\psi_g}\rangle=|{n_0}\rangle.$$

Selecting an LO phase $\theta _k$ uniformly at random from a $2\pi$ interval, and taking $\eta =1$ for these tests, we computed the corresponding probability density [Eq. (6) or (7)] over a predefined grid of points (we chose a resolution $\Delta x = \Delta p = 0.07$) and then drew a single result $x_k$ or $(x_k,p_k)$ according to a multinomial distribution, repeating this procedure $K=8000$ times.

We simulated both homodyne and heterodyne tomography for four ground truth states from each class, with parameters chosen for a mean photon number $\langle {n}\rangle \in \{0,1,2,3\}$. With our chosen cutoff $n_c=10$, the truncation error $\epsilon = 1-\sum _{n=0}^{n_c} \langle {n|\rho _g|n}\rangle$ is less than 0.07 for all cases; however, since the same $n_c$ is used in both data generation and Bayesian estimation, the only impact of this error is to modify the ground truth state in a known way, and thus not the inference accuracy. For each state and measurement subset of size $K\in \{1,400,800,1600,3200,6400,8000\}$, we perform MCMC with a thinning of $T=2^{11}$ and obtain a set of samples $\{\rho ^{(r)}\}$ from which we compute the mean and standard deviation of sample fidelity

(17)$$F^{(r)}=\left(\textrm{Tr} \sqrt{\sqrt{\rho_g} \rho^{(r)} \sqrt{\rho_g}} \right)^{2}$$

for plotting.

In Fig. 9, we show the progression of the fidelity as more measurements—i.e., repeated preparations of the input $\rho _g$—are included in the likelihood. In all cases, the fidelity increases monotonically with $K$ as expected, although the rate varies with both ground truth state and measurement type. The thermal states show the slowest convergence of all state families, consistent with experimental findings in Figs. 4 and 6. The Fock states (d) experience the most striking difference between homodyne and heterodyne measurement schemes, with the former converging in fidelity significantly more rapidly; the other states (a–c) appear less sensitive to the chosen measurement type , although the thermal states still evince noticeably lower tomographic efficiency for heterodyne.

 

Fig. 9. Inference results from simulated datasets. (a) Coherent states. (b) Thermal states. (c) Squeezed vacuum. (d) Fock states. A photon number cutoff of $n_c=10$ is utilized for both simulation and inference.

Download Full Size | PDF

Appendix C: complete cat state results

In Fig. 10, we show the results from our Bayesian QST workflow for all four datasets in Ref. [29], which used avalanche photodiodes (APDs) or a TES to implement one-, two-, and three-photon subtraction on squeezed vacuum to herald coherent-state superpositions (cat states), according to the ideal form

(18)$$|{C_\pm (\alpha)}\rangle = \frac{1}{\sqrt{2(1\pm e^{{-}2|\alpha|^{2}})}}\left(|{\alpha}\rangle \pm |{-\alpha}\rangle \right),$$

where $|{\pm \alpha }\rangle$ are as defined in Eq. (13), and $\pm$ corresponds to the even (odd) cat state. The Wigner function for a given density matrix $\rho$ is defined under our $\hbar =1$ convention as

(19)$$W(x,p) = \sum_m\sum_n \rho_{mn} \frac{({-}1)^{n}}{\pi} \sqrt{\frac{2^{m} n!}{2^{n} m!}} e^{-(x^{2}+p^{2})}[-(x-\mathrm{i} p)]^{m-n} L_n^{m-n}(2(x^{2} + p^{2})),$$

where $L_n^{m-n}(\cdot )$ denotes the generalized Laguerre function.

 

Fig. 10. Wigner functions of Bayesian mean density matrices for all four experimental datasets from Ref. [29]. Labels denote the type of detector used (APD or TES) and the number of photons subtracted (1–3).

Download Full Size | PDF

Following the procedure of Gerrits et al., to each MCMC sample $\rho ^{(r)}$ we assign a cat state amplitude $\alpha ^{(r)}$ as that with maximum overlap:

(20)$$\alpha^{(r)} = \underset{\alpha}{\arg\max} \langle {C_\pm(\alpha)|\rho^{(r)}|C_\pm(\alpha)}\rangle,$$

so that the sample fidelity follows as

(21)$$F^{(r)} = \langle {C_\pm(\alpha^{(r)})|\rho^{(r)}|C_\pm(\alpha^{(r)})}\rangle.$$

We then report $|\alpha |$ and $F$ estimates as $B^{U-B}_{L-B}$, where $B$ is the Bayesian mean, $U$ the 84th percentile, and $L$ the 16th percentile, again following Ref. [29]. Wigner functions of the Bayesian mean density matrices $\rho _B$ for all cases appear in Fig. 10. Our procedure assumes a system efficiency $\eta =0.853$ and cutoff $n_c=10$, for which the truncation error for all expected cat states is less than $7\times 10^{-6}$; $R=1024$ samples are again used for calculations, with a thinning value of $T=2^{13}$ for all cases except APD1, which required $T=2^{15}$ (due to the larger number of quadrature samples $K$).

Estimates for $|\alpha |$ and $F$ with uncertainties are plotted in Fig. 11, along with the MLE findings of Gerrits et al. for comparison. The variation in error box length follows from the diversity in dataset size. APD1 has 324 510 quadrature samples, TES2 has 24 790, APD2 has 41 223, and TES3 has 1087. Both methods show very similar uncertainties for each state, highlighting the congruity between Bayesian uncertainty quantification and parametric bootstrap resampling in these examples. Agreement in the absolute values is also generally good, though the estimated ranges do not overlap in all cases. Given the use of experimental data and differences between Bayesian inference and MLE, there is no fundamental reason to expect perfect agreement. Nevertheless, we do note that the $\alpha$ landscape for optimization via Eq. (20) is quite flat. For example, in Fig. 11(c) we show the fidelity between an ideal even cat state of variable complex amplitude $\alpha$ and the Bayesian mean ($F= \langle {C_+(\alpha )|\rho _B|C_+(\alpha )}\rangle$) for the APD2 result $\rho _B$. Despite the appreciable separation between the Bayesian and MLE estimates for $|\alpha |$, both cat state estimates (shown as dots in the contour plot) possess very similar overlap with $\rho _B$ ($F=0.56$ when rounded), indicating closer similarity between the Bayesian and MLE estimates than the uncertainty intervals in Fig. 11(a) might suggest. This points to the general difficulty in assigning a unique cat state to a particular density matrix $\rho ^{(r)}$.

 

Fig. 11. Comparison of analysis techniques on experimental cat state data. (a) Estimated cat state amplitude $|\alpha |$. (b) Estimated fidelity $F$. The dots represent either the Bayesian mean calculated here or the MLE point estimate from Ref. [29], and the boxes enclose the 16th–84th percentiles. (c) Fidelity between an even cat state $|{C_+(\alpha )}\rangle$ and the Bayesian mean $\rho _B$ for the APD2 dataset. Points show the corresponding estimates from (a).

Download Full Size | PDF

Funding

U.S. Department of Energy (DE-AC05-00OR22725); Advanced Scientific Computing Research (ERKJ353, ERKJ355).

Acknowledgments

We are grateful to T. Gerrits for providing the data from Ref. [29] for use in our analyses. We thank B. T. Kirby, S. Guha, C. N. Gagatsos, and A. J. Pizzimenti for valuable discussions. This work was performed at Oak Ridge National Laboratory, operated by UT-Battelle for the U.S. Department of Energy under contract no. DE-AC05-00OR22725.

Disclosures

The authors declare no conflicts of interest.

Data availability

Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.

References

1. E. Knill, R. Laflamme, and G. J. Milburn, “A scheme for efficient quantum computation with linear optics,” Nature 409(6816), 46–52 (2001). [CrossRef]  

2. P. Kok, W. J. Munro, K. Nemoto, T. C. Ralph, J. P. Dowling, and G. J. Milburn, “Linear optical quantum computing with photonic qubits,” Rev. Mod. Phys. 79(1), 135–174 (2007). [CrossRef]  

3. S. L. Braunstein and P. van Loock, “Quantum information with continuous variables,” Rev. Mod. Phys. 77(2), 513–577 (2005). [CrossRef]  

4. C. Weedbrook, S. Pirandola, R. García-Patrón, N. J. Cerf, T. C. Ralph, J. H. Shapiro, and S. Lloyd, “Gaussian quantum information,” Rev. Mod. Phys. 84(2), 621–669 (2012). [CrossRef]  

5. O. Pfister, “Continuous-variable quantum computing in the quantum optical frequency comb,” J. Phys. B: At., Mol. Opt. Phys. 53(1), 012001 (2020). [CrossRef]  

6. Y. Shih, “Entangled biphoton source - property and preparation,” Rep. Prog. Phys. 66(6), 1009–1044 (2003). [CrossRef]  

7. A. I. Lvovsky and M. G. Raymer, “Continuous-variable optical quantum-state tomography,” Rev. Mod. Phys. 81(1), 299–332 (2009). [CrossRef]  

8. D. F. V. James, P. G. Kwiat, W. J. Munro, and A. G. White, “Measurement of qubits,” Phys. Rev. A 64(5), 052312 (2001). [CrossRef]  

9. D. T. Smithey, M. Beck, M. G. Raymer, and A. Faridani, “Measurement of the Wigner distribution and the density matrix of a light mode using optical homodyne tomography: Application to squeezed states and the vacuum,” Phys. Rev. Lett. 70(9), 1244–1247 (1993). [CrossRef]  

10. D. T. Smithey, M. Beck, J. Cooper, and M. G. Raymer, “Measurement of number-phase uncertainty relations of optical fields,” Phys. Rev. A 48(4), 3159–3167 (1993). [CrossRef]  

11. K. Banaszek, G. M. D’Ariano, M. G. A. Paris, and M. F. Sacchi, “Maximum-likelihood estimation of the density matrix,” Phys. Rev. A 61(1), 010304 (1999). [CrossRef]  

12. A. I. Lvovsky, “Iterative maximum-likelihood reconstruction in quantum homodyne tomography,” J. Opt. B: Quantum Semiclassical Opt. 6(6), S556–S559 (2004). [CrossRef]  

13. A. Smith, C. A. Riofrío, B. E. Anderson, H. Sosa-Martinez, I. H. Deutsch, and P. S. Jessen, “Quantum state tomography by continuous measurement and compressed sensing,” Phys. Rev. A 87(3), 030102 (2013). [CrossRef]  

14. A. Kyrillidis, A. Kalev, D. Park, S. Bhojanapalli, C. Caramanis, and S. Sanghavi, “Provable compressed sensing quantum state tomography via non-convex methods,” npj Quantum Inf. 4(1), 36 (2018). [CrossRef]  

15. E. S. Tiunov, V. V. Tiunova (Vyborova), A. E. Ulanov, A. I. Lvovsky, and A. K. Fedorov, “Experimental quantum homodyne tomography via machine learning,” Optica 7(5), 448–454 (2020). [CrossRef]  

16. C. P. Robert and G. Casella, Monte Carlo Statistical Methods (Springer, 1999).

17. D. J. C. MacKay, Information Theory, Inference, and Learning Algorithms (Cambridge University Press, 2003).

18. R. Blume-Kohout, “Optimal, reliable estimation of quantum states,” New J. Phys. 12(4), 043034 (2010). [CrossRef]  

19. Y.-L. Seah, J. Shang, H. K. Ng, D. J. Nott, and B.-G. Englert, “Monte Carlo sampling from the quantum state space. II,” New J. Phys. 17(4), 043018 (2015). [CrossRef]  

20. C. Granade, J. Combes, and D. G. Cory, “Practical Bayesian tomography,” New J. Phys. 18(3), 033024 (2016). [CrossRef]  

21. B. P. Williams and P. Lougovski, “Quantum state estimation when qubits are lost: a no-data-left-behind approach,” New J. Phys. 19(4), 043003 (2017). [CrossRef]  

22. T. T. Mai and P. Alquier, “Pseudo-Bayesian quantum tomography with rank-adaptation,” J. Stat. Plan. Inference 184, 62–76 (2017). [CrossRef]  

23. J. M. Lukens, K. J. H. Law, A. Jasra, and P. Lougovski, “A practical and efficient approach for Bayesian quantum state estimation,” New J. Phys. 22(6), 063038 (2020). [CrossRef]  

24. J. M. Lukens, K. J. H. Law, and R. S. Bennink, “A Bayesian analysis of classical shadows,” npj Quantum Inf. 7(1), 113 (2021). [CrossRef]  

25. S. Lohani, J. M. Lukens, D. E. Jones, T. A. Searles, R. T. Glasser, and B. T. Kirby, “Improving application performance with biased distributions of quantum states,” Phys. Rev. Res. 3(4), 043145 (2021). [CrossRef]  

26. H.-H. Lu, K. V. Myilswamy, R. S. Bennink, S. Seshadri, M. S. Alshaykh, J. Liu, T. J. Kippenberg, D. E. Leaird, A. M. Weiner, and J. M. Lukens, “Bayesian tomography of high-dimensional on-chip biphoton frequency combs with randomized measurements,” arXiv:2108.04124 (2021).

27. J. DiGuglielmo, C. Messenger, J. Fiurášek, B. Hage, A. Samblowski, T. Schmidt, and R. Schnabel, “Markov chain Monte Carlo estimation of quantum states,” Phys. Rev. A 79(3), 032114 (2009). [CrossRef]  

28. S. L. Cotter, G. O. Roberts, A. M. Stuart, and D. White, “MCMC methods for functions: Modifying old algorithms to make them faster,” Statist. Sci. 28(3), 424–446 (2013). [CrossRef]  

29. T. Gerrits, S. Glancy, T. S. Clement, B. Calkins, A. E. Lita, A. J. Miller, A. L. Migdall, S. W. Nam, R. P. Mirin, and E. Knill, “Generation of optical coherent-state superpositions by number-resolved photon subtraction from the squeezed vacuum,” Phys. Rev. A 82(3), 031802 (2010). [CrossRef]  

30. B. Qi, P. Lougovski, and B. P. Williams, “Characterizing photon number statistics using conjugate optical homodyne detection,” Opt. Express 28(2), 2276–2290 (2020). [CrossRef]  

31. R. A. Campos, B. E. A. Saleh, and M. C. Teich, “Quantum-mechanical lossless beam splitter: SU(2) symmetry and photon statistics,” Phys. Rev. A 40(3), 1371–1384 (1989). [CrossRef]  

32. J. Appel, D. Hoffman, E. Figueroa, and A. I. Lvovsky, “Electronic noise in optical homodyne tomography,” Phys. Rev. A 75(3), 035802 (2007). [CrossRef]  

33. T. Richter, “Determination of photon statistics and density matrix from double homodyne detection measurements,” J. Mod. Opt. 45(8), 1735–1749 (1998). [CrossRef]  

34. M. Takeoka and S. Guha, “Capacity of optical communication in loss and noise with general quantum Gaussian receivers,” Phys. Rev. A 89(4), 042309 (2014). [CrossRef]  

35. K. Życzkowski and H.-J. Sommers, “Induced measures in the space of mixed quantum states,” J. Phys. A: Math. Gen. 34(35), 7111–7125 (2001). [CrossRef]  

36. H.-J. Sommers and K. Życzkowski, “Bures volume of the set of mixed quantum states,” J. Phys. A: Math. Gen. 36(39), 10083–10100 (2003). [CrossRef]  

37. F. Mezzadri, “How to generate random matrices from the classical compact groups,” Not. Am. Math. Soc. 54, 592–604 (2007). [CrossRef]  

38. V. A. Osipov, H.-J. Sommers, and K. Życzkowski, “Random Bures mixed states and the distribution of their purity,” J. Phys. A: Math. Theor. 43(5), 055302 (2010). [CrossRef]  

39. C. R. Müller, C. Peuntinger, T. Dirmeier, I. Khan, U. Vogl, C. Marquardt, G. Leuchs, L. L. Sánchez-Soto, Y. S. Teo, Z. Hradil, and J. Řeháček, “Evading vacuum noise: Wigner projections or Husimi samples?” Phys. Rev. Lett. 117(7), 070801 (2016). [CrossRef]  

40. J. Řeháček, Y. S. Teo, Z. Hradil, and S. Wallentowitz, “Surmounting intrinsic quantum-measurement uncertainties in Gaussian-state tomography with quadrature squeezing,” Sci. Rep. 5(1), 12289 (2015). [CrossRef]  

41. J. Řeháček, D. Mogilevtsev, and Z. Hradil, “Tomography for quantum diagnostics,” New J. Phys. 10(4), 043022 (2008). [CrossRef]  

42. D. N. VanDerwerken and S. C. Schmidler, “Parallel Markov chain Monte Carlo,” arXiv:1312.7479 (2013).

43. B. Calderhead, “A general construction for parallelizing Metropolis–Hastings algorithms,” Proc. Natl. Acad. Sci. 111(49), 17408–17413 (2014). [CrossRef]  

44. P. E. Jacob, J. O’Leary, and Y. F. Atchadé, “Unbiased Markov chain Monte Carlo methods with couplings,” J. R. Stat. Soc. B 82(3), 543–600 (2020). [CrossRef]  

45. J. Heng, A. Jasra, K. J. H. Law, and A. Tarakanov, “On unbiased estimation for discretized models,” arXiv:2102.12230 (2021).

46. U. Leonhardt and H. Paul, “Phase measurement and Q function,” Phys. Rev. A 47(4), R2460–R2463 (1993). [CrossRef]  

47. M. V. Larsen, X. Guo, C. R. Breum, J. S. Neergaard-Nielsen, and U. L. Andersen, “Deterministic multi-mode gates on a scalable photonic quantum computing platform,” Nat. Phys. 17(9), 1018–1023 (2021). [CrossRef]  

48. X. Guo, C. R. Breum, J. Borregaard, S. Izumi, M. V. Larsen, T. Gehring, M. Christandl, J. S. Neergaard-Nielsen, and U. L. Andersen, “Distributed quantum sensing in a continuous-variable entangled network,” Nat. Phys. 16(3), 281–284 (2020). [CrossRef]  

49. B. J. Lawrie, P. D. Lett, A. M. Marino, and R. C. Pooser, “Quantum sensing with squeezed light,” ACS Photonics 6(6), 1307–1318 (2019). [CrossRef]