Simulating Gaussian boson sampling quantum computers

In standard GBS, $N$ single-mode squeezed vacuum states are sent into an $M$-mode linear photonic network. If $N\neq M$, the remaining $N-M$ ports are vacuum states. If each input mode is independent, the entire input state is defined via the density matrix

$$\hat{\rho}^{(\text{in})}=\prod_{j}^{M}\left|r_{j}\right\rangle\left\langle r_{j}\right|.$$

(1)

Here, $\left|r_{j}\right\rangle=\hat{S}(r_{j})\left|0\right\rangle$ is the squeezed vacuum state and $\hat{S}(r_{j})=\exp(r_{j}(\hat{a}_{j}^{\dagger(\text{in})})^{2}/2-r_{j}(\hat{a}_{j}^{(\text{in})})^{2}/2)$ is the squeezing operator [43, 44, 45]. We assume the squeezing phase is zero and $\hat{a}^{\dagger(\text{in})}$ $\left(\hat{a}_{j}^{(\text{in})}\right)$ is the input creation (annihilation) operator for the $j$-th mode. The vacuum state $\left|0\right\rangle$ corresponds to letting $r_{j}=0$ in the squeezing operator.

The ideal GBS assumes the inputs are pure squeezed states [46]. These non-classical states are characterized by their quadrature variance and for normal ordering are defined as [47, 48]

	$\displaystyle\Delta_{x_{j}}^{2}$	$\displaystyle=2\left(n_{j}+m_{j}\right)=e^{2r_{j}}-1$
	$\displaystyle\Delta_{y_{j}}^{2}$	$\displaystyle=2\left(n_{j}-m_{j}\right)=e^{-2r_{j}}-1.$		(2)

Here, $\Delta_{x_{j}}^{2}=\left\langle(\Delta\hat{x}_{j})^{2}\right\rangle=\left\langle\hat{x}_{j}^{2}\right\rangle-\left\langle\hat{x}_{j}\right\rangle^{2}$ is the simplified variance notation for the input quadrature operators $\hat{x}_{j}=\hat{a}_{j}^{(\text{in})}+\hat{a}_{j}^{\dagger(\text{in})}$ and $\hat{y}_{j}=-i(\hat{a}_{j}^{(\text{in})}-\hat{a}_{j}^{\dagger(\text{in})})$, which obey the commutation relation $[\hat{x}_{j},\hat{y}_{k}]=2i\delta_{jk}$.

Squeezed states, while still satisfying the Heisenberg uncertainty principle $\Delta_{x_{j}}^{2}\Delta_{y_{j}}^{2}=1$, have altered variances such that one quadrature has variance below the vacuum noise level, $\Delta_{x_{j}}^{2}<1$. By the Heisenberg uncertainty principle, this causes the variance in the corresponding quadrature to become amplified well above the vacuum limit, $\Delta_{y_{j}}^{2}>1$ [43, 49, 50]. Other squeezing orientations are available and remain non-classical so long as one variance is below the vacuum limit.

For pure squeezed states, the number of photons in each mode is $n_{j}=\left\langle\hat{a}_{j}^{\dagger(\text{in})}\hat{a}_{j}^{(\text{in})}\right\rangle=\sinh^{2}(r_{j})$, whilst the coherence per mode is $m_{j}=\left\langle(\hat{a}_{j}^{(\text{in})})^{2}\right\rangle=\cosh(r_{j})\sinh(r_{j})=\sqrt{n_{j}\left(n_{j}+1\right)}$.

Linear photonic networks made from a series of polarizing beamsplitters, phase shifters and mirrors act as large scale interferometers, and crucially conserve the Gaussian nature of the input state.

In the ideal lossless case, as seen in Fig.(1), photonic networks are defined by an $M\times M$ Haar random unitary matrix $\boldsymbol{U}$. The term Haar corresponds to the Haar measure, which in general is a uniform probability measure on the elements of a unitary matrix [51].

For lossless networks, input modes are converted to outputs following

$$\hat{a}_{i}^{(\text{out})}=\sum_{j=1}^{M}U_{ij}\hat{a}_{j}^{(\text{in})},$$

(3)

where $\hat{a}_{i}^{(\text{out})}$ is the output annihilation operator for the $i$-th mode. Therefore, each output mode is a linear combination of all input modes.

Practical applications will always include some form of losses in the network, for example photon absorption. This causes the matrix to become non-unitary, in which case its denoted by the transmission matrix $\boldsymbol{T}$, and the input-to-output mode conversion now contains additional loss terms [52, 53]:

$$\hat{a}_{i}^{(\text{out})}=\sum_{j=1}^{N}T_{ij}\hat{a}_{j}^{(\text{in})}+\sum_{j=1}^{M}B_{ij}\hat{b}_{j}^{(\text{in})}.$$

(4)

Here, $\hat{b}_{j}^{(\text{in})}$ is the annihilation operator for the $j$-th loss mode whilst $\boldsymbol{B}$ is a random noise matrix. The loss matrix conserves the total unitary of the system as $\boldsymbol{T}\boldsymbol{T}^{\dagger}+\boldsymbol{B}\boldsymbol{B}^{\dagger}=\boldsymbol{I}$, where $\boldsymbol{I}$ is the identity matrix.

Although linear networks are conceptually simple systems, they introduce computational complexity from the exponential number of possible interference pathways photons can take in the network. These are generated from the beamsplitters. As an example, the input ports of a 50/50 beamsplitter are defined as superpositions of the output modes:

	$\displaystyle\hat{a}_{1}^{\dagger(\text{in})}$	$\displaystyle=\frac{1}{\sqrt{2}}\left(\hat{a}_{3}^{\dagger(\text{out})}+\hat{a}_{4}^{\dagger(\text{out})}\right)$
	$\displaystyle\hat{a}_{2}^{\dagger(\text{in})}$	$\displaystyle=\frac{1}{\sqrt{2}}\left(\hat{a}_{3}^{\dagger(\text{out})}-\hat{a}_{4}^{\dagger(\text{out})}\right).$		(5)

Classical states such as coherent and thermal states input into large scale linear networks are readily simulated on a classical computer [34]. This leads to the question, why are non-classical states such as Fock and squeezed states computationally hard to simulate on large scale photonic networks?

In short, non-classical states input into a linear network generate large amounts of quantum correlations, which are non-trivial to simulate classically if combined with photo-detection. Continuing with the beamsplitter example above, if input ports $1$ and $2$ now contain indistinguishable Fock states, as is the case with standard boson sampling, the input state can be written as

$$\left|1,1\right\rangle_{1,2}=\frac{1}{\sqrt{2}}\left(\left|2,0\right\rangle_{3,4}-\left|0,2\right\rangle_{3,4}\right).$$

(6)

The absence of a $\left|1,1\right\rangle_{3,4}=\hat{a}_{3}^{\dagger(\text{out})}\hat{a}_{4}^{\dagger(\text{out})}\left|0,0\right\rangle_{3,4}$ term means the output photons are bunched, that is, they are entangled and always arrive in pairs [54]. This phenomena is known as the Hong-Ou-Mandel (HOM) interference effect and is named after the authors who first observed this phenomena [55].

If input ports $1,2$ have input squeezed vacuum states, one generates two-mode Einstein-Podolsky-Rosen (EPR) entanglement in the quadrature phase amplitudes [56, 57]. More details, including a summary of the derivation, can be found in [58, 31].

In GBS experiments, photon count patterns, denoted by the vector $\boldsymbol{c}$, are generated by measuring the output state $\hat{\rho}^{(\text{out})}$ . The corresponding distribution changes depending on the type of detector used.

The original proposal for GBS [7] involved normally ordered photon-number-resolving (PNR) detectors that can distinguish between incoming photons and implement the projection operator [54, 59]

$$\hat{p}_{j}(c_{j})=\frac{1}{c_{j}!}:(\hat{n}^{\prime}_{j})^{c_{j}}e^{-\hat{n}^{\prime}_{j}}:.$$

(7)

Here, $:\dots:$ denotes normal ordering, $c_{j}=0,1,2,\dots$ is the number of counts in the $j$-th detector and $\hat{n}^{\prime}_{j}=\hat{a}_{j}^{\dagger(\text{out})}\hat{a}_{j}^{(\text{out})}$ is the output photon number. Practically, a PNR detector typically saturates for a maximum number of counts $c_{\text{max}}$, the value of which varies depending on the type of PNR detector implemented.

When the mean number of output photons per mode satisfies $\left\langle\hat{n}^{\prime}_{j}\right\rangle\ll 1$, PNR detectors are equivalent to threshold, or click, detectors which saturate for more than one photon at a detector. Threshold detectors are described by the projection operator [59]

$$\hat{\pi}_{j}(c_{j})=:e^{-\hat{n}^{\prime}_{j}}\left(e^{\hat{n}^{\prime}_{j}}-1\right)^{c_{j}}:,$$

(8)

where $c_{j}=0,1$ regardless of the actual number of photons arriving at the detector.

In situations where $\left\langle\hat{n}^{\prime}_{j}\right\rangle\gg 1$, threshold detectors become ineffective since they cannot accurate discriminate between the different numbers of incoming photons. However, by sending photons through a demultiplexer (demux) before arriving at the detector, the output pulses of light become ’diluted’ such that the mean photon number is small enough to be accurately counted using threshold detectors [60, 61]. In this approach the demux acts as a type of secondary optical network of size $1\times M_{S}$, where $M_{S}$ is the number of secondary modes. For output photons to become sufficiently dilute, one must satisfy $M_{S}\gg\left\langle\hat{n}^{\prime}_{j}\right\rangle$ [60, 61, 62]. When this is the case, the detectors can be thought of pseudo-PNR (PPNR), because the threshold detectors are accurately measuring photon counts, thus becoming equivalent to PNR detectors [62].

The computational complexity of GBS arises in the probabilities of count patterns $\boldsymbol{c}$. These probabilities are determined using a straightforward extension of the projection operators for both PNR and click detectors. In the case of PNR detectors, the operator for a specific photon count pattern is

$$\hat{P}(\boldsymbol{c})=\bigotimes_{j=1}^{M}\hat{p}_{j}(c_{j}),$$

(9)

the expectation value of which corresponds to the $\#P$-hard matrix Hafnian [7, 33]

$$\left\langle\hat{P}(\boldsymbol{c})\right\rangle\approx|\text{Haf}(\boldsymbol{D}_{S})|^{2}.$$

(10)

Here, $\boldsymbol{D}_{S}$ is the sub-matrix of $\boldsymbol{D}=\boldsymbol{T}\left(\bigoplus_{j=1}^{M}\tanh(r_{j})\right)\boldsymbol{T}^{T}$, corresponding to the channels with detected photons, where the superscript $T$ denotes the transpose of the lossy experimental transmission matrix.

The Hafnian is a more general function than the matrix permanent, although both are related following [7, 33]

$$\text{perm}(\boldsymbol{D})=\text{Haf}\left(\begin{array}[]{cc}0&\boldsymbol{D}\\ \boldsymbol{D}^{T}&0\end{array}\right).$$

(11)

Therefore, one of the advantages of GBS with PNR detectors is that one can compute the permanent of large matrices from the Hafnian instead of performing Fock state boson sampling.

For threshold detectors, the operator for a specific binary count pattern $\boldsymbol{c}$ is

$$\hat{\Pi}(\boldsymbol{c})=\bigotimes_{j=1}^{M}\hat{\pi}_{j}(c_{j}),$$

(12)

whose expectation value corresponds to the Torontonian [8]

$$\left\langle\hat{\Pi}(\boldsymbol{c})\right\rangle=\text{Tor}(\boldsymbol{O}_{S}),$$

(13)

where $\boldsymbol{O}_{S}$ is a matrix constructed from the covariance matrix of the output state. The Torontonian has also been shown to be $\#P$-hard to compute, and is directly related to the Hafnian [8].

Only a few short years after the original proposal of GBS with threshold detectors [8], the first large scale experimental network, called Jiuzhang 1.0, was implemented by Zhong et al [15]. In this experiment, $N=50$ single-mode squeezed states were sent into a lossy $M=100$ mode linear network which generated over 50 million binary count patterns in $\approx 200\text{s}$.

The experimentalists estimated it would take Fugaku, the largest supercomputer at the time, $600$ million years to generate the same number of samples using the fastest exact classical algorithm at the time [20]. Exact samplers are algorithms designed to replicate an experiment by directly sampling the Torontonian or Hafnian using a classical supercomputer. When experiments claim quantum advantage, they typically do so by comparing the run-time of the experiment with such exact samplers.

The algorithm implemented in [20] uses a chain rule sampler that samples each mode sequentially by computing the conditional probability of observing photons in the $k$-th mode given the photon probabilities of the previous $k-1$-th modes. The computational complexity of this algorithm scales as $\mathcal{O}(N_{D}^{3}2^{N_{D}})$ where $N_{D}=\sum_{j=1}^{M}c_{j}$ is the total number of detected photon counts. Due to this scaling with detected counts and system size, more recent experiments aim to beat these exact samplers by increasing the number of observed counts. The mean number of counts per pattern in Jiuzhang 1.0 was $\approx 41$, and the largest number of counts in a single pattern was $76$.

Shortly after this initial experiment, in order to reduce the probability of an exact sampler replicating the experiment, Zhong et al [16] implemented an even larger GBS network, named Jiuzhang 2.0, which increased the number of modes to $M=144$ and performed multiple experiments for different input laser powers and waists. This expectedly produced an increase in the number of observed clicks, with a mean click rate of $\approx 68$ and a maximum observed click number of $N_{D}=113$.

At the same time, in classical computing, the increased probability of multiple photons arriving at the same detector was exploited by Bulmer et al [22] to improve the scaling of the chain rule sampler, achieving $\mathcal{O}(N_{D}^{3}2^{N_{D}/2})$ for GBS with both PNR and click detectors. Despite this apparently modest improvement, it was estimated that generating the same number of samples as Jiuzhang 1.0 on Fugaku would now only take $\sim 73$ days [22]. Due to the substantial speed-up over previous exact algorithms [20, 21], it was predicted that experiments using either PNR or click detectors needed $N_{D}\gtrsim 100$ to surpass exact samplers [22].

To reduce the probability of multiple photons arriving at a single detector, Deng et al [18] added a $1\times 8$ demux to the output modes of Jiuzhang 2.0, dubbing this upgraded network Jiuzhang 3.0. The demux is made of multiple fiber loop beamsplitters that separate photons into four temporal bins and two spatial path bins and aims to ensure $\left\langle\hat{n}^{\prime}_{j}\right\rangle\apprle 1$, increasing the likelihood the threshold detectors are accurately counting the oncoming photons. This simple addition generated patterns with almost double the largest number of clicks obtained in Jiuzhang 2.0, with one experiment observing a maximum of $N_{D}=255$ with a mean click number of more than $100$.

The linear networks implemented in the three experiments above are all static, that is, once the networks have been fabricated, one cannot reprogram any of the internal optics. This has the advantage of circumventing the exponential accumulation of loss that arises from increased depth [17, 63], although one sacrifices programmability.

To this end, Madsen et al [17] implemented a programmable GBS called Borealis using three fiber loop variable beamsplitters, each containing a programmable phase-shifter. In this experiment, $N$ squeezed state pulses are generated at a rate of $6\text{MHz}$ and are sent into each variable beamsplitter sequentially, allowing pulses in different time-bins to interfere. Photons output from the final fibre loop are then sent into a $1\times 16$ demultiplexer and then counted by $16$ PNR detectors, which are superconducting transition edge sensors (TES).

The demultiplexer is required for two reasons. The first is to partially ’dilute’ the number of photons in each output pulse for the PNR detectors used [17]. The second is to ensure each TES has reached its baseline operating temperature before another pulse is detected, ensuring an accurate determination of photon numbers.

The two largest networks using this setup sent $N=216$ and $N=288$ input pulses into the network, detecting a mean photon number of $\approx 125$ and $\approx 143$, respectively. However, this network was more susceptible to losses than the previous systems [15, 16, 18].

There are a number of “exact” classical simulation methods that generate equivalent photon-counts, where we use quotes for “exact”, because even these methods have round-off errors due to to finite arithmetic. These algorithms are all exponentially slow [20, 21, 22]. It is known that the non-classicality of the squeezed input states of GBS is essential to providing quantum advantage in the GBS experiments, since classical states which have a positive Glauber-Sudarshan P-representation [72, 73] are known to have classically simulable photon counts. This will be treated in Section (III.2.2). However, the non-Gaussianity arising from the photon-detection measurement process is also important.

This is because the measurement set-up where quadrature-phase-amplitude measurements are made on a system described by Gaussian states can be modeled by a positive Wigner function, and hence becomes classically simulable [34, 35]. Gaussian states are defined as those with a Gaussian Wigner function. We note that squeezed states, while non-classical, are Gaussian states. Examples of the classical simulation of entanglement for systems with a positive Wigner function are well known [74, 75, 76, 77, 78, 79].

The role of the measurement set-up in GBS is clarified by the work of Chabaud, Ferrini, Grosshans and Markham [65] and Chabaud and Walschaers [66]. These authors illustrate a connection between quantum advantage and the non-Gaussianity of both the input states and measurement set-up. In [65], conditions sufficient for the strong simulation of Gaussian quantum circuits with non-Gaussian input states are derived. Non-Gaussian states are those for which the Gaussian factorization of correlation functions is not applicable [80].

By demonstrating a mapping between bosonic quantum computation and continuous-variable sampling computation, where the measurement comprises a double quadrature detection, Chabaud and Walschaers [66] adapt classical algorithms derived in [66] to derive a general algorithm for the strong simulation of bosonic quantum computation, which includes Gaussian boson sampling. They prove that the complexity of this algorithm scales with a measure of non-Gaussianity, the stellar rank of both the input state and measurement set-up. This enables a quantification of non-Gaussianity, including from the nature of the measurement, as a resource for achieving quantum advantage in bosonic computation.

Originally developed by Wigner for symmetrically ordered operators [89], phase-space representations establish a mapping between quantum operators of different orderings to probability distributions defined on the classical phase-space variables of position and momentum [90, 54], which are more commonly rewritten in terms of the complex coherent state amplitude vectors $\boldsymbol{\alpha}$, $\boldsymbol{\alpha}^{*}$.

Moyal [91] later showed how one can use these methods can be used to compute the dynamics of operators. Due to this, phase-space representations are frequently used to compute the operator master equation [92, 90], which for some representations, corresponds to the second-order Fokker-Planck equation (FPE), which is commonly used in statistical mechanics.

The FPE in turn, can be mapped to a stochastic differential equation (SDE) that introduces randomly fluctuating terms and, for some real or complex vector $\boldsymbol{a}$, takes the general form [93]

$$\frac{d}{dt}\boldsymbol{a}=\boldsymbol{A}(\boldsymbol{a})+\boldsymbol{B}(\boldsymbol{a})\xi(t),$$

(16)

where $\boldsymbol{A}$ is a vector function and $\boldsymbol{B}$ is a matrix, both of which are typically known, while $\xi(t)$ is a real Gaussian noise term with $\left\langle\xi\right\rangle=0$ and $\left\langle\xi_{i}(t)\xi_{j}(t^{\prime})\right\rangle=\delta(t-t^{\prime})\delta_{ij}$. These randomly fluctuating terms, defined as the derivative of a Wiener process [93], play an analogous role to quantum and thermal fluctuations for applications in quantum mechanics.

Each solution of an SDE is a stochastic trajectories in phase-space and physically corresponds to a single experiment. Therefore, averages over an ensemble of trajectories $E_{S}$ corresponds to a mean value from multiple experimental shots.

Such dynamical methods have been successfully applied to a variety of quantum optics systems [94, 95, 96, 97] (also see Ref.([71]) for a brief review), including the CIM [26, 27]. However, as is clear from subsection II.2, linear networks are not modeled dynamically and hence cannot produce an SDE.

Instead, phase-space representations model linear networks as stochastic difference equations (SDFE) where the $M$-th order SDFE takes the general form [98]

$$\boldsymbol{a}_{n+1}=\sum_{m=0}^{M}\boldsymbol{A}_{m}(n-m,\boldsymbol{a}_{n-m})+\boldsymbol{B}(n,\boldsymbol{a}_{n})\xi_{n},$$

(17)

where $\boldsymbol{a}_{n}$, $\boldsymbol{A}_{m}$, $\boldsymbol{B}$ and $\xi_{n}$are discrete analogs to their continues variable definitions in Eq.(16). This becomes clearer when comparing with Eq.(4).

Due to the randomly fluctuating term, SDEs don’t have analytical solutions and must be computed numerically using a variety of methods [99]. As is also the case with numerical computations of non-stochastic differential equations, SDEs are approximated as difference equations for numerical computation. Although one usually has practical issues such small step size limits for SDEs [99, 98], SDEs and SDFEs are two sides of the same coin.

Due to the use of PNR and threshold detectors, we focus on simulating linear networks using the normally ordered Glauber-Sudarshan diagonal P-representation and generalized P-representations. Non-normally ordered Wigner and Q-function methods are possible also [30, 31, 32], but these have too large a sampling error to be useful for simulations of photon counting.

For GBS with threshold detectors, the number of possible binary patterns obtained from an experiment is $\approx 2^{M}$ with each pattern having a probability of roughly $\left\langle\hat{\Pi}(\boldsymbol{c})\right\rangle\approx 2^{-M}$. Samples from large scale networks are exponentially sparse, requiring binning to obtain meaningful statistics. It is therefore necessary to define a grouped count observable that corresponds to an experimentally measurable quantity.

The most general observable of interest is a $d$-dimensional grouped count probability (GCP), defined as [30]

$$\mathcal{G}_{\boldsymbol{S}}^{(n)}(\boldsymbol{m})=\left\langle\prod_{j=1}^{d}\left[\sum_{\sum c_{i}=m_{j}}\hat{\Pi}_{S_{j}}(\boldsymbol{c})\right]\right\rangle,$$

(24)

which is the probability of observing $\boldsymbol{m}=(m_{1},\dots,m_{d})$ grouped counts of output modes $\boldsymbol{M}=(M_{1},M_{2},\dots)$ contained within disjoint subsets $\boldsymbol{S}=(S_{1},S_{2},\dots)$. Here, similar to Glauber’s original definition from intensity correlations [83], $n=\sum_{j=1}^{d}M_{j}\leq M$ is the GCP correlation order.

Each grouped count $m_{j}$ is obtained by summing over binary patterns $\boldsymbol{c}$ for modes contained within a subset $S_{j}$. For example, a $d=1$ dimensional GCP, typically called total counts, generates grouped counts as $m_{j}=\sum_{i}^{M}c_{i}$, whilst $d>1$ gives $\boldsymbol{m}=(m_{1}=\sum_{i=1}^{M/d}c_{i},\dots,m_{d}=\sum_{i=M/d+1}^{M}c_{i})$. This definition also includes more traditional low-order marginal count probabilities and moments.

Although a variety of observables can be generated using GCPs, such as the marginal probabilities commonly used by spoofing algorithms [23, 24], multi-dimensional GCPs are particularly useful for comparisons with experiment. The increased dimension causes the number of grouped counts to grow, e.g. for $d=2$ $\boldsymbol{m}=(m_{1},m_{2})$, which in turn increases the number of count bins (data points) available for comparisons, providing a fine grained comparison of results. This also causes effects of higher order correlations to become more statistically significant in the data [32].

One the most useful applications arises by randomly permuting the modes within each subset $S_{j}$, which results in a different grouped count for each permutation. The number of possible permutation tests scales as [32]

$$\frac{\binom{M}{M/d}}{d}=\frac{M!}{d(M/d)!(M-M/d)!},$$

(25)

where different correlations are compared in each comparison test. This leads to a very large number of distinct tests, making good use of the available experimental data.

Despite these advantages, caution must be taken when comparing very high dimensional GCPs, since the increased number of bins means that there are fewer counts per bin. This causes the experimental sampling errors to increase, reducing the significance of statistical tests [32].

A more realistic model is one that includes decoherence in the input state. Such a model was proposed in [30] and assumes the intensity of the input light is weakened by a factor of $1-\epsilon$ whilst adding $n_{j}^{\text{th}}=\epsilon n_{j}$ thermal photons per mode. The number of input photons per mode remains unchanged from the pure squeezed state definition, but the coherence of photons is now altered as $\tilde{m}=(1-\epsilon)m_{j}$. To account for possible measurement errors, a correction factor $t$ is also applied to the transmission matrix, to improve the fit to the simulated distribution.

The advantage of this decoherence model, which is a type of thermal squeezed state [104, 47], is that it allows a simple method for generating phase-space samples for any Gaussian state following

	$\displaystyle\alpha_{j}$	$\displaystyle=\frac{1}{2}\left(\Delta_{x_{j}}w_{j}+i\Delta_{y_{j}}w_{j+M}\right)$
	$\displaystyle\beta_{j}$	$\displaystyle=\frac{1}{2}\left(\Delta_{x_{j}}w_{j}-i\Delta_{y_{j}}w_{j+M}\right),$		(26)

where $\left\langle w_{j}w_{k}\right\rangle=\delta_{jk}$ are real Gaussian noises that model quantum fluctuations in each quadrature and

	$\displaystyle\Delta_{x_{j}}^{2}$	$\displaystyle=2(n_{j}+\tilde{m}_{j})$
	$\displaystyle\Delta_{y_{j}}^{2}$	$\displaystyle=2(n_{j}-\tilde{m}_{j}),$		(27)

are the thermal squeezed state quadrature variances.

As $\epsilon\rightarrow 1$ the inputs become classical, since $\Delta_{x_{j}}^{2},\,\Delta_{y_{j}}^{2}\geq 1$, but small amounts of thermalization causes the inputs to remain non-classical given, for example, with a squeezing orientation of $\Delta_{x_{j}}^{2}<0$. So long as the state is Gaussian, one can model a variety of inputs, both classical and non-classical, by simply varying $\epsilon$, where $\epsilon=0$ corresponds to a pure squeezed state and $\epsilon=1$ to a thermal state.

The input stochastic samples are also straightforwardly extended to non-normally ordered representations, as outlined in detail in [30, 31, 32].

Performing the summation over binary patterns can now be efficiently simulated in phase-space using the multi-dimensional inverse discrete Fourier transform [30]

$$\mathcal{G}_{\boldsymbol{S}}^{(n)}\left(\boldsymbol{m}\right)=\frac{1}{\prod_{j}(M_{j}+1)}\sum_{\boldsymbol{k}}\tilde{\mathcal{G}}_{M}^{(n)}(\boldsymbol{k})e^{i\sum_{j}k_{j}\theta_{j}m_{j}},$$

(28)

where the Fourier observable simulates all possible correlations generated in an experimental network and is defined as

$$\mathcal{\tilde{G}}_{\boldsymbol{S}}^{(n)}\left(\boldsymbol{k}\right)=\left\langle\prod_{j=1}^{d}\bigotimes_{i\in S_{j}}\left(\pi_{i}(0)+\pi_{i}(1)e^{-ik_{j}\theta_{j}}\right)\right\rangle_{P}.$$

(29)

Here, $\theta_{j}=2\pi/(M_{j}+1)$, $k_{j}=0,\dots,M_{j}$ and $\pi_{j}=\exp(-n^{\prime}_{j})(\exp(n^{\prime}_{j})-1)^{c_{j}}$ is the positive-P click observable, obtained from the equivalence Eq.(20), where $n^{\prime}_{j}=\alpha^{\prime}_{j}\beta^{\prime}_{j}$ is the output photon number.

This simulation method is highly scalable, with most observables taking a few minutes to simulate on a desktop computer for current experimental scales. To highlight this scalability, much larger simulations of network sizes of up to $M=16,000$ modes have also been performed [30].

If classical states are input into a GBS, the network can be efficiently simulated using the diagonal P-representation, which arises as a special case of the positive P-representation if $P(\boldsymbol{\alpha},\boldsymbol{\beta})=P(\boldsymbol{\alpha})\delta(\boldsymbol{\alpha}^{*}-\boldsymbol{\beta})$. Due to this, initial stochastic samples are still generated using Eq.(26), except now one has rotated to a classical phase-space with $\boldsymbol{\beta}=\boldsymbol{\alpha}^{*}$.

Similar to non-classical states, the input density matrix for any classical state is defined as

$$\hat{\rho}^{(\text{in})}=\int P(\boldsymbol{\alpha})\left|\boldsymbol{\alpha}\right\rangle\left\langle\boldsymbol{\alpha}\right|\text{d}^{2}\boldsymbol{\alpha},$$

(30)

where the form of the distribution $P(\boldsymbol{\alpha})$ changes depending on the state. Input amplitudes are again transformed to outputs following $\boldsymbol{\alpha}^{\prime}=\boldsymbol{T}\boldsymbol{\alpha}$ and are used to define the output state

$$\hat{\rho}^{(\text{out})}=\int P(\boldsymbol{\alpha})\left|\boldsymbol{\alpha}^{\prime}\right\rangle\left\langle\boldsymbol{\alpha}^{\prime}\right|\text{d}^{2}\boldsymbol{\alpha}.$$

(31)

Using the output coherent amplitudes, one can now efficiently generate binary patterns corresponding to any classical state input into a linear network. In order to conserve the simulated counts for each ensemble, corresponding to a single experimental shot, the $j$-th output mode of the $k$-th ensemble is independently and randomly sampled via the Bernoulli distribution [32]

$$P_{j}^{(k)}(c_{j}^{(\text{class)}})=(p_{j}^{(k)})^{c_{j}^{(\text{class)}}}(1-p_{j}^{(k)})^{1-c_{j}^{(\text{class)}}}.$$

(32)

Here, $c_{j}^{(\text{class})}=0,1$ is the classically generated bit of the $j$-th mode where the probability of $c_{j}^{(\text{class})}=1$, the ’success’ probability, is

$$p_{j}^{(k)}=(\pi_{j}(1))^{(k)}=\left(1-e^{-n^{\prime}_{j}}\right)^{(k)}.$$

(33)

This is simply the click probability of the $k$-th stochastic ensemble with an output photon number of $n^{\prime}_{j}=\left|\alpha_{j}\right|^{2}$.

For an $M$-mode network, each stochastic ensemble outputs the classical faked pattern of the form

$$(\boldsymbol{c}^{(\text{class})})^{(k)}=[P_{1}^{(k)},P_{2}^{(k)},\dots,P_{M}^{(k)}],$$

(34)

which are binned to obtain GCPs, denoted $\mathcal{G}^{(\text{class})}$, that can be compared to simulated distributions.

In the limit of large numbers of patterns, binned classical fakes are approximately equal to the simulated output distribution corresponding to the density matrix Eq. (31). An example of this can be seen in Fig. (3) where $4\times 10^{6}$ binary patterns are generated by sending $N=20$ thermal states into an $M=N$ mode lossless linear network represented by a Haar random unitary. The binned classical patterns produce a total count distribution that agrees with simulations within sampling error, because for this classical input example there is no quantum advantage.

For comparisons to be meaningful, statistical tests are needed to quantify differences between experiment and theory. To this end, Refs. [30, 31] implement a slightly altered version of a standard chi-square test, which is frequently used in random number generator validation [105, 106], and is defined as

$$\chi_{ES}^{2}=\sum_{i=1}^{k}\frac{\left(\mathcal{G}_{i,E}-\bar{\mathcal{G}}_{i,S}\right)^{2}}{\sigma_{i}^{2}}.$$

(35)

Here, $k$ is the number of valid photon count bins, which we define as a bin with more than $10$ counts, and $\mathcal{\bar{G}}_{i,S}$ is the phase-space simulated GCP ensemble mean for any input state $S$, which converges to the true theoretical GCP, $\mathcal{G}_{i,S}$, in the limit $E_{S}\rightarrow\infty$. The experimental GCP is defined as $\mathcal{G}_{i,E}$, and $\sigma_{i}^{2}=\sigma_{T,i}^{2}+\sigma_{E,i}^{2}$ is the sum of experimental, $\sigma_{E,i}^{2}$, and theoretical, $\sigma_{T,i}^{2}$, sampling errors of the $i$-th bin.

The chi-square test result follows a chi-square distribution that converges to a normal distribution when $k\rightarrow\infty$ due to the central limit theorem [107, 108]. One can use this to introduce another test that determines how many standard deviations away the result is from its mean. The aim of this test is to obtain the probability of observing the output $\chi_{ES}^{2}$ result using standard probability theory. For example, an output of $6\sigma$, where $\sigma$ is the standard deviation of the normal distribution, indicates the data has a very small probability of being observed.

To do this, Ref.[32] implemented the Z-score, or Z-statistic, test which is defined as

$$Z_{ES}=\frac{\left(\chi_{ES}^{2}/k\right)^{1/3}-\left(1-2/(9k)\right)}{\sqrt{2/(9k)}},$$

(36)

where $\left(\chi_{ES}^{2}/k\right)^{1/3}$ is the Wilson-Hilferty (WH) transformed chi-square statistic [107], which allows a faster convergence to a normal distribution when $k\geq 10$ [107, 108], and $\mu=1-\sigma^{2}$, $\sigma^{2}=2/(9k)$ are the corresponding mean and variance of the normal distribution.

The Z-statistic allows one to determine the probability of the count patterns. A result of $Z_{ES}>6$ would indicate that experimental distributions are so far from the simulated results that patterns may be displaying non-random behavior. Due to this, and to present a unified comparison notation, we compute the Z-statistic of the chi-square results presented in [30] using Eq.(36).

Valid experimental results with correct distributions should have $\chi_{EI}^{2}/k\approx 1$, where the subscript $I$ denotes the ideal GBS distribution, with $Z_{EI}\approx 1$. For claims of quantum advantage, one must simultaneously prove that $Z_{CI}\gg 1$, where the subscript $C$ denotes binary patterns from the best classical fake, such as the diagonal-P method described above. This is the ’gold standard’ and would show that, within sampling errors, experimental samples are valid, and closer to the ideal distribution than any classical fake.

In the more realistic scenario of thermalized squeezed inputs, one may still have quantum advantage if $Z_{ET}\approx 1$ while $Z_{CT}\gg 1$. Therefore, these four observables are of the most interest for comparisons of theory and experiment, and are given below.

Throughout this section, we primarily review comparisons from Jiuzhang 1.0, for purposes of illustration [30]. A thorough comparison of all data sets obtained in Jiuzhang 2.0 is presented elsewhere [32].

We first review comparisons of total counts, which is the probability of observing $m$ clicks in any pattern and is usually one of the first observables experimentalists compare samples to. This is because in the limit of a large number of clicks, one can estimate the ideal distribution as a Gaussian distribution via the central limit theorem.

To simulate the ideal total counts distribution in phase-space using GCPs we let $n=M$ and $\boldsymbol{S}=\{1,2,\dots,M\}$, giving $\mathcal{G}_{\{1,\dots,M\}}^{(M)}(m)$. For Jiuzhang 1.0, one obtains a Z-statistic of $Z_{EI}\approx 340$ for $k=63$ valid bins. Clearly, experimental samples are far from the ideal and indicate photon count distributions are not what would be expected from an ideal GBS.

To determine whether these differences either increase or decrease when higher-order correlations become more prevalent in the simulations, the dimension of the GCP is increased to $d=2$. In this case, Jiuzhang 1.0 sees an almost doubling in Z values for comparisons with the simulated ideal distribution (see Fig.(4)), where $Z_{EI}\approx 647$ is obtained from $k=978$. [30]. The increase in the number of valid bins with GCP dimension causes the normal distribution of the WH transformed chi-square statistic to have a smaller variance. When compared to a single dimension, where $k=63$ gives $\sigma^{2}\approx 3.5\times 10^{-3}$, binning with $d>1$ causes experimental samples to now pass a more stringent test, as $k=978$ produces a normal distribution with variance $\sigma^{2}\approx 2.3\times 10^{-4}$.

Simulating the more realistic scenario of thermalized squeezed states, a thermalized component of $\epsilon=0.0932\pm 0.0005$ and a transmission matrix correction of $t=1.0235\pm 0.0005$ is used as to compare a simulated model to samples from Jiuzhang 1.0 [30]. In this case, an order of magnitude improvement in the resulting Z value is observed where $Z_{ET}\approx 17\pm 2$.

Despite this significant improvement, differences are still large enough that $Z_{ET}>1$. Similar results were also obtained for data sets with the largest recorded photon counts in Jiuzhang 2.0, that is data sets claiming quantum advantage [32]. However, this is not the case for data sets with small numbers of recorded photons which typically give $Z_{ET}\approx 1$ [32], although these experiments should be easily replicated by exact samplers. The large amount of apparent thermalization is the likely reason why simulated squashed states described Jiuzhang 1.0 samples just as well as the ideal GBS in Ref. [85].

When higher-order correlations are considered, samples from Jiuzhang 1.0 are far from the expected correlation moments of the ideal distribution. Although including the above fitting parameters improves this result with $Z_{ET}\approx 31$, the samples still deviate significantly from the theoretical thermalized distribution [30].

Unlike comparisons of total counts, samples from all data sets in Jiuzhang 2.0 satisfy $Z_{EI},Z_{ET}>1$ for simulations with $d>1$, meaning higher-order correlations also display significant departures from the expected theoretical results, even for the simpler cases with low numbers of photon counts [32]. The reasons for this are not currently known.

Phase-space GBS simulation software xQSim can be downloaded from public domain repositories at [88].

The Authors declare no competing interests.

This research was funded through grants from NTT Phi Laboratories and the Australian Research Council Discovery Program.

A.S.D wrote the original draft, created figures and performed numerical simulations. All authors contributed to reviewing and editing this manuscript.