Matrix phase-space representations in quantum optics

TL;DR

This work addresses the challenge of verifying quantum computational advantage in large-scale Gaussian boson sampling by introducing matrix phase-space representations (matrix-P) that incorporate global symmetries via coherent projection matrices. The approach unifies previous normally-ordered phase-space methods with stochastic gauges through symmetry projections, enabling a complete and positive representation. The parity-projected matrix-P framework dramatically reduces sampling errors, achieving around $10^{-3}$ relative accuracy with millions of samples for networks with up to $10^{4}$ modes, and enabling efficient verification of lossless GBS where traditional methods struggle. This method offers a scalable, accurate tool for validating large quantum photonic experiments and has potential to extend to other symmetry-constrained quantum systems.

Abstract

We introduce matrix quantum phase-space distributions. These extend the idea of a quantum phase-space representation via projections onto a density matrix of global symmetry variables. The method is applied to verification of low-loss Gaussian boson sampling (GBS) quantum computational advantage experiments with up to 10,000 modes, where classically generating photon-number counts is exponentially hard. We demonstrate improvements in sampling error by a factor of 1000 or more compared to unprojected methods, which are infeasible for such cases.

Figures (6)

• Figure 1: Matrix-P simulations of total photo-count probability (solid lines) versus the exact multi-mode squeezed state photon counting distribution for GBS (dashed lines) for squeezed states with squeezing parameter $\boldsymbol{r}=[0.5,\dots,0.5]$ and a Haar random unitary matrix of size $M=200$. Numerical probabilities are obtained by averaging over ensembles of size $E_{S}=1.2\times10^{6}$ . Sampling errors and difference errors are both $<6\times10^{-5}$, and are not visible. A size of $M=10^{6}$ reduced the errors to $\sim6\times10^{-6}$, with a maximum probability of 0.0096.
• Figure 2: Positive-P phase-space simulations of total photo-count probability (solid lines) versus the exact distribution (dashed lines) for pure state GBS. All other parameters as in Fig(1). Difference errors are $\sim0.03$, and are clearly visible. This is caused by a skewed distribution leading to sampling errors, requiring enormous sample numbers to reach full convergence, as explained in the SM.
• Figure 3: Comparisons of positive-P simulations of the total count distribution versus the exact multi-mode squeezed state photon counting distribution Eq.(\ref{['eq:Exact_multi_mode_squeezed_state']}) (solid black line) for two GBS set-ups where pure squeezed states with uniform squeezing parameter $\boldsymbol{r}=[0.5,\dots,0.5]$ are transformed by a Haar random unitary matrix $\boldsymbol{U}$ of size a) $M=N=50$ and b) $M=N=20$. Positive-P moments are obtained by averaging over sample ensembles of size $E_{S}=2.4\times10^{5}$ (solid orange line) and $E_{S}=2.4\times10^{8}$ (solid blue line). Upper and lower lines correspond to $\pm1\sigma_{T,j}$ theoretical sampling errors for the $j$-th photon count bin.
• Figure 4: The $M=N=20$ mode lossless network with $\boldsymbol{r}=[0.5,\dots,0.5]$ of Fig.(\ref{['fig:Exact_+P_PNR_comp-1']}b) is simulated using the matrix-P representation for an ensemble size of $E_{S}=1.2\times10^{6}$ (solid blue line) and compared to the exact photon counting distribution (dashed black line) Eq.(\ref{['eq:Exact_multi_mode_squeezed_state']}). Sampling and difference errors are of the order $<10^{-3}$, where the matrix-P moments converge to their exact values for all $m$.
• Figure 5: Logarithmic plot of estimated positive-P probability densities $P(m)=P(P_{m})$ versus binned stochastic trajectories of $P_{m}$ for $m=4$ (solid blue line), such that $P_{4}=\frac{1}{4!}\left(n\right)^{4}e^{-n}$, and $m=3$ (solid orange line), with $P_{3}=\frac{1}{3!}\left(n\right)^{3}e^{-n}$. Trajectories are obtained from simulations of the lossless $M=N=20$ mode GBS network with $E_{S}=2.4\times10^{8}$, whose total count distribution is presented in Fig.(\ref{['fig:Exact_+P_PNR_comp-1']}b). $N_{b}=4\times10^{3}$ bins are used to obtained density estimates, on a range $[-100,100]$ with spacing $\Delta_{b}=0.05$. Upper and lower lines correspond to $\pm1\sigma_{T,j}$, which increase as $P_{m}\rightarrow\pm\infty$ due to the exponentially small probabilities one is required to resolve. The inset is a close-up of the bin range $[-0.5,0.5]$, where the probability densities $P(3)$ overlap with the $P(4)$ probabilities. The bins correspond to the largest probability to obtain a specific trajectory value.