Regimes of classical simulability for noisy Gaussian boson sampling

TL;DR

This paper establishes a classical simulation threshold for noisy Gaussian Boson Sampling, showing how noise and loss affect quantum advantage and proposing input squeezing as a mitigation strategy.

Contribution

It formalizes a sufficient condition for classical simulability of noisy GBS and analyzes how noise impacts quantum advantage in linear-optical architectures.

Findings

Noisy GBS becomes classically simulable under certain noise conditions.

Quantum advantage diminishes with exponential photon loss in deep circuits.

Increasing input squeezing can help GBS evade classical simulation.

Abstract

As a promising candidate for exhibiting quantum computational supremacy, Gaussian Boson Sampling (GBS) is designed to exploit the ease of experimental preparation of Gaussian states. However, sufficiently large and inevitable experimental noise might render GBS classically simulable. In this work, we formalize this intuition by establishing a sufficient condition for approximate polynomial-time classical simulation of noisy GBS --- in the form of an inequality between the input squeezing parameter, the overall transmission rate and the quality of photon detectors. Our result serves as a non-classicality test that must be passed by any quantum computationalsupremacy demonstration based on GBS. We show that, for most linear-optical architectures, where photon loss increases exponentially with the circuit depth, noisy GBS loses its quantum advantage in the asymptotic limit. Our results thus delineate intermediate-sized regimes where GBS devices might considerably outperform classical computers for modest noise levels. Finally, we find that increasing the amount of input squeezing is helpful to evade our classical simulation algorithm, which suggests a potential route to mitigate photon loss.

Equations

$$\text{sech}\left[\frac{1}{2}\Theta \left(\ln \left(\frac{1-2q_D}{\eta e^{-2r}+1-\eta }\right)\right)\right]>e^{-\frac{{ϵ}^2}{4K}}.$$

$$\eta <{\eta }_{\infty }+\frac{1-2q_D}{1-e^{-2r}}\frac{ϵ}{\sqrt{2K}}+\omega \left(\frac{1}{\sqrt{K}}\right).$$

$$W_{\rho }^{(t)}(\mathbf{x})=\frac{\exp \left[-\frac{1}{2}{\mathbf{x}}^T({\mathbf{V}}_{\rho }-t{\mathbf{I}}_2{)}^{-1}\mathbf{x}\right]}{2\pi \sqrt{\det ({\mathbf{V}}_{\rho }-t{\mathbf{I}}_2)}}.$$

$${\rho }_{\text{out}}=\mathcal{U}\left({\sigma }^{\otimes K}\otimes |0⟩⟨0{|}^{\otimes (M-K)}\right){\mathcal{U}}^{†},$$

$$\tilde{\rho }=\mathcal{U}\left({\tau }^K\otimes |0⟩⟨0{|}^{M-K}\right){\mathcal{U}}^{†}$$

$$D_T(P,\tilde{P})\le \sqrt{\frac{2}{\alpha }{D}_{\alpha }(P\Vert \tilde{P})}\le \sqrt{\frac{2}{\alpha }{D}_{\alpha }(\rho \Vert \tilde{\rho })}.$$

$$\frac{2}{\alpha }{D}_{\alpha }(\sigma \Vert \tau )<\frac{{ϵ}^2}{K}.$$

$$-\ln \left[F_{\max }(\eta ,q_D)\right]<\frac{{ϵ}^2}{4K},$$

$$s_{\sigma }$$

$$n_{\sigma }$$

$$F$$

$$\Delta$$

$$\Lambda$$

$$F_{\max }(\eta ,q_D)=\text{sech}\left[\frac{1}{2}\Theta \left(\ln \left(\frac{1-2q_D}{\eta e^{-2r}+1-\eta }\right)\right)\right].$$

$$⟨P_{\text{collision}}{⟩}_{\mathcal{U}}\le \frac{8N^2}{M},$$

$$F(S)=\left(\genfrac{}{}{0ex}{}{\frac{K}{2}+S-1}{S}\right){\text{sech}}^K(r){\tanh }^S(r).$$

$$⟨P_{\text{collision}}{⟩}_{\mathcal{U}}$$

$$M=O(K^2).$$

$$\text{sech}(r)\sum _{S=0}^{\infty }(\tanh r{)}^S|2S,2S⟩.$$

$$S(s,\varphi )=\exp \left[\frac{1}{2}s{e}^{i\varphi }{a}^{†2}-\frac{1}{2}s{e}^{-i\varphi }{a}^2\right].$$

$$\left(\begin{array}{cc}\cosh r+\cos \varphi \sinh r& -\sin \varphi \sinh r\\ -\sin \varphi \sinh r& \cosh r-\cos \varphi \sinh r\end{array}\right).$$

$${\mathbf{V}}_{\rho }={\mathbf{V}}_{STS}(s_{\rho },n_{\rho },{\varphi }_{\rho })=(2n_{\rho }+1)\left(\begin{array}{cc}\cosh 2s_{\rho }-\cos {\varphi }_{\rho }\sinh 2s_{\rho }& -\sin {\varphi }_{\rho }\sinh 2s_{\rho }\\ -\sin {\varphi }_{\rho }\sinh 2s_{\rho }& \cosh s_{\rho }+\cos {\varphi }_{\rho }\sinh 2s_{\rho }\end{array}\right).$$

$$s_{\rho }\le \frac{1}{2}\ln \frac{2n_{\rho }+1}{t}.$$

$$s_{\sigma }$$

$$n_{\sigma }$$

$$-\ln \left[F_{\max }(\eta ,q_D)\right]\le \frac{{ϵ}^2}{4K},$$

$$F_{\max }(\eta ,q_D):=\underset{t\in [1-2q_D,1]}{\max }\underset{\tau \in {\mathcal{C}}_G^{(t)}}{\max }F(\sigma ,\tau ).$$

$$F(\sigma ,\tau )=\frac{1}{\sqrt{\Delta +\Lambda }-\sqrt{\Lambda }},$$

$$\Delta$$

$$\Lambda$$

$$F_{\max }(\eta ,q_D)=\text{sech}\left[\frac{1}{2}\Theta \left(\ln \left(\frac{1-2q_D}{\eta e^{-2r}+1-\eta }\right)\right)\right].$$

$$D_{\alpha }({\rho }_1\Vert {\rho }_2)=\frac{1}{\alpha -1}\ln \mathrm{Tr}\left\{{\left({\rho }_1^{\frac{1-\alpha }{2\alpha }}{\rho }_2{\rho }_1^{\frac{1-\alpha }{2\alpha }}\right)}^{\alpha }\right\},$$

$$D_{\alpha }({\rho }_1\Vert {\rho }_2)\ge D_{\alpha }(\mathcal{E}({\rho }_1)\Vert \mathcal{E}({\rho }_2)),$$

Full text

Quantum computational supremacy  Harrow and Montanaro (2017); Preskill (2018)—when a quantum device performs a computational task beyond the capabilities of classical computers—is a long-anticipated milestone towards practical quantum computation. A recent publication has claimed to have achieved such milestone on a quantum processor with 53 superconducting qubits Arute et al. (2019). However, subsequent work managed to significantly speed up the classical simulation of this multi-qubit device Pednault et al. (2019). This suggests that quantum computational supremacy, instead of being a one-shot experimental proof, will be the result of a long-term competition between continuously improved quantum devices and faster classical simulations.

Since near-term quantum devices do not benefit from error correction, in a realistic scenario noise increases with the size of the experiment and potentially negates its quantum advantage. Numerous works have studied the effects of noise on Boson Sampling (BS), including partial photon distinguishability Rohde (2015); Shchesnovich (2015); Renema et al. (2018a), fabrication imperfections  Kalai and Kindler (2014); Leverrier and García-Patrón (2015); Arkhipov (2015), losses Rahimi-Keshari et al. (2016); Oszmaniec and Brod (2018); García-Patrón et al. (2019); Renema et al. (2018b); Aaronson and Brod (2016), and detector dark counts Rahimi-Keshari et al. (2016). In contrast with BS, no rigorous analysis has been carried out for noisy Gaussian Boson Sampling (GBS), a variant of BS which finds applications in specific computational problems such as dense subgraph searching Arrazola and Bromley (2018); Banchi et al. (2019), perfect matching counting Brádler et al. (2018a), graph isomorphism Brádler et al. (2018b) and the simulation of vibrational spectra Huh et al. (2015). In addition, GBS is related to non-Gaussian probabilistic state engineering Sabapathy et al. (2019); Su et al. (2019); Gagatsos and Guha (2019).

Concretely, GBS Hamilton et al. (2017) describes a computational task in which photon statistics is directly measured from a Gaussian state. An arbitrary GBS instance can be implemented by: (i) deterministic preparation of $K$ single-mode squeezed vacuum states; (ii) interference over an $M$-mode interferometer Reck et al. (1994); Clements et al. (2016); (iii) sampling of output statistics by photon number resolving detectors. The deterministic sources, together with its high generation probability and sampling rate, render GBS a promising alternative to other proposals.

In this work, we take the first step towards answering the following questions: what is the minimum squeezing needed to demonstrate quantum advantage with GBS? How much photon loss would render the device classically simulable? What level of dark counts and inefficiencies can be tolerated in the detectors? These questions are crucial for a deeper understanding of the origin of quantum complexity in GBS, as well as for the design of any experiment aiming at a demonstration of quantum computational supremacy.

Main result and ramifications — Our main result states that a noisy GBS device (see Fig. 1(a) for a schematic) can be classically efficiently simulated up to error $\epsilon$ if the following condition is satisfied:

[TABLE]

Here the photon source is characterized by the squeezing parameter $r$, and $\eta$ is the overall transmission rate. The quality of the photon detectors is quantified by $q_{D}:=p_{D}/\eta_{D}$, where $p_{D}$ is the dark count rate and $\eta_{D}$ is the quantum efficiency. Finally, $\Theta$ is the ramp function $\Theta(x):=\max(x,0)$.

We briefly summarize our classical algorithm as follows:

If inequality (1) holds, proceed; otherwise exit the algorithm. 2. 2.

Calculate the classical Gaussian state which is closest to the output state of the noisy device. We derive an efficient analytical formula for this procedure. 3. 3.

Sample from the classical Gaussian state according to the method given in Ref. Rahimi-Keshari et al. (2016), which can be done efficiently.

See Sec. VI of the Supplemental Material sup for a detailed description of our algorithm, which includes Refs. Arkhipov and Kuperberg (2012); Bremner et al. (2010); Lund et al. (2014); Aaronson ; Brod (2015); Marian et al. (2002); Spedalieri et al. (2012); Müller-Lennert et al. (2013); Wilde et al. (2014); Frank and Lieb (2013); Beigi (2013); Van Erven and Harremos (2014); Seshadreesan et al. (2018); Wilde et al. (2017); Killoran et al. (2019); Qir .

A finite-sized experiment will count as strong evidence towards quantum computational supremacy when the best known classical algorithm to simulate it, running on a classical supercomputer, takes a reasonably large amount of time (for further discussion, see Bremner et al. (2016); Neville et al. (2017); Boixo et al. (2018)). Therefore, inequality (1), for reasonably small $\epsilon$, defines a natural and useful test that any alleged quantum advantage demonstration with noisy GBS must pass.

To illustrate this, we test inequality (1) on several recent small-scale GBS experiments Clements et al. (2018); Paesani et al. (2019); Zhong et al. (2019). For instance, in Ref. Paesani et al. (2019), $K=4$ squeezed vacuum states with $r\approx 0.1$ were input into a 12-mode random-walk circuit with overall transmission $\eta=0.088$ and detector efficiency $\eta_{D}=0.78$. For typical superconducting nanowire single photon detectors, the dark count rate is $p_{D}\approx 10^{-4}$ Hadfield et al. (2005). Using these numbers, the experiment can be efficiently simulated by our algorithm with error $2\%$. A simple python implementation takes $2$ms to output a sample on a laptop . In the most recent GBS demonstration, $K=6$ squeezed vacuum states with an average $r\approx 0.38$ are coupled into a 12-mode interferometer implemented with bulk optics in free space Zhong et al. (2019). With overall transmission $\eta=0.89$ 111The transmission rate of the interferometer is 0.99, which can be found in in Ref. Zhong et al. (2018). We assume that the efficiency of coupling the output light to the detector is $0.9$, which is reported from a recent BS experiment Wang et al. (2019). This gives an overall transmission around $0.89$, $\eta_{D}=0.82$ and assuming $p_{D}=10^{-4}$, we find that inequality (1) has no solution for any $\epsilon\in[0,1]$, which makes this experiment pass our classicality test.

We can also study the behavior of inequality (1) in the asymptotic limit. Assuming that $q_{D}$ and $r$ are constants as $K$ increases, in the asymptotic limit of large $K$, we have the following condition for efficient classical simulability

[TABLE]

Here $\eta_{\infty}:=q_{D}(1+\coth r)$, and $\eta<\eta_{\infty}$ is the corresponding condition that follows for exact sampling ($\epsilon=0$) obtained previously in Ref. Rahimi-Keshari et al. (2016). See Fig. 2(a) and (b) for plots of this inequality for various squeezing levels, transmission rates, and detector qualities.

From our result cast into the form (2), we observe the following: (i) by including some error tolerance in the classical simulation, our result improves on previous algorithms of Ref. Rahimi-Keshari et al. (2016) and gives tighter bounds on noise parameters for finite-size experiments. (ii) our result allows more general types of noise, including a combination of photon loss and dark counts, compared to previous analyses for standard BS García-Patrón et al. (2019); Oszmaniec and Brod (2018); (iii) increasing the input squeezing reduces the upper bound on the transmission rate, which underpins the notion of squeezing as a non-classical resource.

An important scenario which allows comparison with previous BS results is that of pure photon loss by assuming perfect detectors, $q_{D}=0$. In this case, for exact sampling $(\epsilon=0)$, inequality (1) gives trivial result $\eta<0$—a major drawback of the methods in Ref. Rahimi-Keshari et al. (2016). This can be understood by observing that a squeezed state, under the effect of pure loss, never becomes exactly a classical state (though it approximates one arbitrarily well).

In Sec.II of the Supplemental Material sup , we argue that exact sampling of GBS with arbitrary loss is computationally hard. The seemingly counter-intuitive result suggests that the notion of exact simulability is often too stringent to be considered in the error analysis of quantum devices, as even a real-world experiment is incapable of sampling from the exact theoretical distribution. Therefore, below we analyze the effect of pure loss by considering the more relevant approximate simulation.

In a realistic setup, the input squeezing (and the energy) per mode does not scale with $K$. Multiplying both sides of inequality (2) by the average photon number $\bar{N}=K\sinh^{2}r$, we find that the inequality is satisfied if $\eta\bar{N}<O(\epsilon\sqrt{\bar{N}})$. Thus, lossy GBS can be efficiently simulated when the average number of surviving photons is less than $O(\sqrt{\bar{N}})$, which matches the scaling obtained for standard BS Oszmaniec and Brod (2018); García-Patrón et al. (2019).

In most linear-optical architectures, photon loss is defined by unit depth of the circuit, leading to an exponential decrease of transmission with the depth. Our results imply that, asymptotically, GBS implemented on these platforms is rendered efficiently simulable by losses if the depth is linear in the number of modes. It is easy to see that condition (2) is always satisfied for these circuits when they have super-logarithmic depth, which also holds for BS García-Patrón et al. (2019); Oszmaniec and Brod (2018). In the other extreme, for planar circuits (i.e., with only nearest-neighbor beam splitters) of logarithmic depth we can construct a tensor network simulation that runs in quasi-polynomial time. This is analogous to a similar algorithm for BS García-Patrón et al. (2019). The only difference is that we need to introduce an additional cutoff on the Hilbert space, for large photon numbers, that does not degrade too much the Gaussian state nor slows down the simulation.

From this asymptotic analysis we conclude that lossy GBS inevitably loses its quantum advantage for sufficiently large sizes, at least within most current architectures. Nonetheless, a large running time advantage might be still be possible for GBS devices with intermediate sizes and modest (but non-negligible) noise rates. A similar situation has been discussed in the context of qubits Gao and Duan (2018).

Gaussian states and classicality— Our results crucially follow from the fact that classical Gaussian states can be sampled efficiently Rahimi-Keshari et al. (2016). Here we briefly review relevant definitions and a few well-known facts.

Consider a single-mode bosonic system and let $\hat{{\bm{x}}}:=(\hat{q},\hat{p})$ denote the row vector of its quadrature operators. They satisfy the commutation relations ($\hbar=2$) $\left[\hat{x}_{j},\hat{x}_{k}\right]=2i\bm{\Omega}_{jk}~{},$ where ${\bm{\Omega}}=\left(\begin{smallmatrix}0&-1\\ 1&0\end{smallmatrix}\right)$. A state $\rho$ is Gaussian if it is fully characterized by its mean vector $\bar{{\bm{x}}}_{\rho}=\operatorname{Tr}\left\{\rho\hat{{\bm{x}}}\right\}$ and its covariance matrix $\bm{V}_{\rho}^{j,k}=\frac{1}{2}\operatorname{Tr}{}\left\{\rho\left[\hat{x}_{j}-\bar{x}_{\rho}^{j},\hat{x}_{k}-\bar{x}_{\rho}^{k}\right]_{+}\right\}~{},$ where $\left[\cdot,\cdot\right]_{+}$ is the anti-commutator Weedbrook et al. (2012). All distance measures we consider are minimized when two Gaussian states are displaced along the same direction by the same amount Seshadreesan et al. (2018), so we set $\bar{{\bm{x}}}=0$ from now on.

The ${(t)}$-ordered phase-space quasi-probability distribution [$({t})$-PQD] of a single mode Gaussian state is given by

[TABLE]

For ${t}=-1,0$ and $1$, we obtain the Husimi, Wigner and Glauber-Sudarshan functions, respectively. Equation (56) only holds when ${\bm{V}}_{\rho}-t{\bm{I}}_{2}$ is positive definite. Since the covariance matrix ${\bm{V}}_{\rho}$ is positive definite, there always exists $t^{*}_{\rho}\in[0,1]$ such that: for $t<t^{*}_{\rho}$ the $(t)$-PQD is a Gaussian function; for $t=t^{*}_{\rho}$, the $(t)$-PQD has $\delta$-function singularities; and for $t>t^{*}_{\rho}$, the $(t)$-PQD does not exist. A Gaussian state is called classical if its $P$ function (i.e., for $t=1$) is well-behaved. We slightly generalize this notion of classicality and say that $\rho$ is $(t)$-classical if

${\bm{V}}_{\rho}-t{\bm{I}}_{2}$ is positive definite, or equivalently, if its $(t)$-PQD is non-singular (i.e., if $t\leq t^{*}_{\rho}$). Finally, we denote $\mathcal{C}_{G}^{(t)}$ as the set of $(t)$-classical Gaussian states.

Noise model—The input squeezed light is usually noisy due to photon loss in its preparation Vernon et al. (2018). We model such source imperfection as an ideal squeezed vacuum state followed by a lossy channel with transmission $\eta_{S}$. An $M$-mode lossy interferometer is mathematically described by a sub-unitary transformation matrix ${\bm{A}}$. By assuming uniform loss, which is usually a good approximation for integrated platforms, the interferometer can be simplified to $M$ single-mode lossy channels, each with transmission $\eta_{I}$, followed by an ideal unitary transformation ${\bm{U}}$ (see García-Patrón et al. (2019); Oszmaniec and Brod (2018) for a rigorous treatment). That is, we can write ${\bm{A}}=\eta_{I}\bm{U}$. As a final simplification, we combine the source and interferometer losses into a single loss channel with transmission $\eta:=\eta_{S}\eta_{I}$ and encapsulate it into a mixed input Gaussian state $\sigma$ (see Fig. 1).

The current proofs for hardness of (Gaussian) Boson Sampling require a no-collision regime (sup, , Sec. I), where the probability of detecting more than one photon per mode is negligible. In this regime we can replace photon counting detectors by threshold detectors without loss of generality Rahimi-Keshari et al. (2016); Quesada et al. (2018). Each threshold detector has sub-unity quantum efficiency denoted by $\eta_{D}$ and registers random dark counts with probability $p_{D}$. Following Ref. Rahimi-Keshari et al. (2016); Barnett et al. (1998) we describe them by the POVM elements: $\Pi_{0}=(1-p_{D})\sum_{n=0}^{\infty}(1-\eta_{D})^{n}|n\rangle\langle n|$ and $\Pi_{1}=I-\Pi_{0}$. We define $q_{D}=p_{D}/\eta_{D}$ as the figure of merit that characterizes the detectors.

With this noise model, the probability distribution at the output reads $P({\bm{n}})=\operatorname{Tr}\left\{\rho_{\text{out}}\Pi_{{\bm{n}}}\right\}$, with ${\bm{n}}=(n_{1},n_{2},\ldots,n_{M})$ and $n_{i}\in\{0,1\}$, where

[TABLE]

and $\mathcal{U}$ is the Hilbert space unitary associated to the unitary matrix $\bm{U}$ acting in mode space. Proving the conditions for efficient classical sampling from these distribution, namely inequality (1), is the subject of the remainder of this paper.

Approximate simulation of noisy GBS—In Ref. Rahimi-Keshari et al. (2016), the authors studied the problem of exact sampling from an $M$-mode quantum state of the form

[TABLE]

using threshold detectors of quality $q_{D}$. They proved this task can be simulated exactly with polynomial running time if $\tau$ is a $(t)$-classical state, for $t\in[1-2q_{D},1]$ (See (sup, , Sec. VI) for a brief review of this result). In this work, we restrict our consideration to $(t)$-classical Gaussian states for a given $q_{D}$.

It is natural to expect that, when the mixed input state $\sigma$ is close to some $(t)$-classical Gaussian state $\tau$, the corresponding noisy GBS experiment can be efficiently simulated with small error. Since any such state $\tau$ leads to an efficient classical simulation, we minimize the distance to $\sigma$ over all possible choices of $\tau$. It is this intuition that eventually leads to our main result given by condition (1).

We start our formal derivation by connecting the distance between the two output distributions to the distance between the two output states. The former is usually measured by the total variation distance, $D_{T}(\tilde{P},P):=\frac{1}{2}\sum_{{\bm{n}}}\left|\tilde{P}({\bm{n}})-P({\bm{n}})\right|$. For the distance between states we choose the sandwiched Rényi relative entropy due to its generality, denoted by $D_{\alpha}(\sigma\|\tau)$, for $\alpha\in[\frac{1}{2},1)$ Müller-Lennert et al. (2013); Wilde et al. (2014). The two measures can be connected by the following inequality chain:

[TABLE]

We first used a generalized Pinkser’s inequality and, in a slight abuse of notation, $D_{\alpha}(P\|\tilde{P})$ also denotes the Rényi divergence between two distributions Van Erven and Harremos (2014). The second inequality follows from the data-processing inequality under quantum measurements Müller-Lennert et al. (2013).

The last quantity can be greatly simplified due to the similar structures of Eqs. (4) and (5). Since the sandwiched Rényi relative entropy is invariant under the action of the unitary transformation and is additive under tensor products Müller-Lennert et al. (2013), we have that $D_{\alpha}(\rho\|\tilde{\rho})=KD_{\alpha}(\sigma\|\tau)$. Substituting this into Eq. (6) we see that $D_{T}(P,\tilde{P})\leq\epsilon$ whenever

[TABLE]

We could optimize this bound over $\alpha\in\left[\tfrac{1}{2},1\right)$ — since $D_{\alpha}$ is non-decreasing over $\alpha$ Müller-Lennert et al. (2013), we expect an optimal $\alpha^{*}$ exists. For simplicity we only consider the case $\alpha=\tfrac{1}{2}$, which admits analytical calculations; a numerical optimization over $\alpha$ can be found in (sup, , Sec. V). In this special case, we have $D_{\frac{1}{2}}(\sigma,\tau)=-\ln F(\sigma,\tau)$, where $F(\sigma,\tau):=\operatorname{Tr}\left\{\sqrt{\sqrt{\sigma}\tau\sqrt{\sigma}}\right\}^{2}$ is the quantum fidelity. As we mentioned earlier, we will have a tighter bound if $\tau$ is further optimized over $C^{(t)}_{G}$ with $t\in[1-2q_{D},1]$. We thus have

[TABLE]

with $F_{\max}(\eta,q_{D}):=\max_{\tau\in\mathcal{C}_{G}^{(t)},t\in[1-2q_{D},1]}F(\sigma,\tau)$.

The argument above shows that a noisy GBS can be efficiently simulated with error no more than $\epsilon$ if the above inequality holds; in other words, we have reduced our multi-mode problem into a single-mode optimization.

Optimization of the fidelity—To perform the optimization, it is convenient to use the squeezed thermal states (STS) parameterization Marian and Marian (1993). Any single-mode Gaussian state $\rho$ can be written as $\rho=S(s_{\rho},\phi_{\rho})\rho_{T}^{n_{\rho}}S^{\dagger}(s_{\rho},\phi_{\rho})$. Here $\rho_{T}^{n_{\rho}}$ is a thermal state with average photon number ${n_{\rho}}$, and $S(s_{\rho},\phi)=\exp\left[\frac{1}{2}s_{\rho}e^{i\phi_{\rho}}\hat{a}^{\dagger 2}-\frac{1}{2}s_{\rho}e^{-i\phi_{\rho}}\hat{a}^{2}\right]$ is the squeezing operator with squeezing parameter $s_{\rho}$ and phase rotation $\phi_{\rho}$. For the lossy squeezed state $\sigma$ appearing in our noise model (Fig. 1), according to Ref. sup , Sec. II, we have $\phi_{\sigma}=0$ and

[TABLE]

A straightforward calculation using an analytical expression of the quantum fidelity between two single-mode Gaussian states  Holevo (1975) leads to

[TABLE]

Since $\tau$ is a $(t)$-classical Gaussian state, it is straightforward to show that this implies $s_{\tau}\leq\frac{1}{2}\ln\frac{2n_{\tau}+1}{t}$ (sup, , Sec. III). After solving this constrained optimization problem, we obtain that the maximum fidelity is $\text{sech}\left[\Theta\left(s_{\sigma}-\frac{1}{2}\ln\frac{2n_{\sigma}+1}{t}\right)\right]$  (sup, , Sec. IV). Finally, we optimize this quantity again over $t\in[1-2q_{D},1]$. Clearly, the maximum value is achieved at $t=1-2q_{D}$. Using Eqs. (9)-(10) we finally arrive at the following solution:

[TABLE]

Plugging this into Eq. (8) we attain our main result for classical simulability of noisy GBS, expressed in Eq. (1). This inequality draws a boundary in the parameter space $(\eta,q_{D},K,\epsilon)$. In Fig. 2 (a) and (b) we plot them in the $\eta-K$ plane for several values of the squeezing parameter $r$ and detector quality $q_{D}$. For noise parameters below the corresponding solid line, GBS can be efficiently simulated with error $\epsilon\leq 0.01$. We observe that the region where quantum computational supremacy has not been ruled out expands as we increase the quality of detectors or the amount of input squeezing; And the latter underpins the notion of squeezing as a non-classical resource. For bounds obtained by using general sandwiched Rényi relative entropies, our numerical simulation results shown in Fig. 2(c) and (d) suggest that the quantum fidelity (corresponding to $\alpha=1/2$) is optimal, while quantum relative entropy (corresponding to $\alpha\rightarrow 1$) gives the worst bound. For more details see the Supplemental Material sup .

Conclusion.—In this work, we established a non-classicality condition that any quantum computational supremacy demonstration based on GBS must satisfy, connecting important experimental parameters, such as the squeezing in the photon source, the overall photon transmission rate, and the dark count rates in photon detection. Our result shows that lossy GBS can be efficiently simulated when the average number of surviving photons is quadratically related to the number of input photons, which matches the scaling obtained for BS García-Patrón et al. (2019); Oszmaniec and Brod (2018).

Motivated by recent works for BS Renema et al. (2018b, a); Moylett et al. (2019), exploiting the properties of the combinatorial quantity associated to the statistics of bosons is an important direction of future research, which could lead to similar improvement on classical simulation of noisy GBS. Another interesting topic for future research is whether non-uniformity of photon loss can be leveraged to improve our simulation Brod and Oszmaniec (2019).

I Non-collision regime for GBS

For standard boson sampling with $N$ single-photon inputs and $M$ modes, it was shown Aaronson and Arkhipov (2011); Arkhipov and Kuperberg (2012) that the probability of detecting collision events is bounded as follows:

[TABLE]

where the average is over Haar-random unitaries. The proof of Eq. (13) relies on the fact that Haar-random unitaries map any $N$-photon, $M$-mode state onto the maximally mixed state (its density matrix is given by the identity on the corresponding Hilbert space).

For GBS, the input state has indefinite photon number. Specifically, the probability of generating $S$ photon pairs is given by Hamilton et al. (2017)

[TABLE]

Therefore, from Eq. (13), the probability of detecting collision events at the output of GBS satisfies

[TABLE]

$F(S)$ in Eq. (14) is a negative binomial (or Pascal) distribution. In the large $K$ limit, $F(S)$ converges to a Gaussian distribution with mean value $\frac{K}{2}\sinh^{2}r$ and variance $\frac{K}{2}\sinh^{2}r\cosh^{2}r$, which gives $\langle S^{2}\rangle_{F}=O(K^{2})$. Therefore, we also expect no-collision outputs in GBS to dominate whenever

[TABLE]

II Exact lossy GBS is hard

In this section we give evidence that exact classical simulation of a lossy GBS device cannot be efficient, unless the polynomial hierarchy collapses to its third level. The post-selection based argument we use is standard and was used to prove similar claims for many different restricted models of quantum computation (see e.g. Bremner et al. (2010); Aaronson and Arkhipov (2011)), so we only detail the parts of the argument that pertain to GBS. The construction we use is directly inspired by the scattershot boson sampling model Lund et al. (2014); Aaronson , though our purpose is different, as we are interested especially in the effect of losses and the complexity of the model.

For this proof, we assume that losses are uniform within the interferometer $\bm{U}$, and so we can move all losses to the end (this is a standard assumption that is a good approximation for e.g. integrated photonic devices, but was also shown to hold under more general conditions Oszmaniec and Brod (2018)). In contrast to the results in the main paper, here we also ignore all other sources of losses. This is an important caveat, but can be justified as losses in photon sources and detectors are effectively constant, whereas losses inside the interferometer $\bm{U}$ scale with its depth (which, for boson sampling, also typically scales with the number of photons). Therefore, photon loss within the linear optical network is the main scalability botteleneck. We leave it as an open question whether this caveat can be eliminated.

Theorem 1

If there is an efficient classical algorithm to sample from the output distribution of a lossy Gaussian boson sampling instance exactly (or up to multiplicative error), then the polynomial hierarchy collapses to its third level.

Proof. Consider the following lossy GBS setup. We prepare $2K$ identical SMSV states with identical squeezing paramater $r$. These states are input, in pairs, into 50:50 beam splitters, generating two-mode squeezed vacuum (TMSV) states of the form

[TABLE]

For each TMSV state we couple one mode directly to a number-resolving detector (these are the heralding registers H), whereas the other half are sent into the lossy interferometer $\bm{U}$ (which may also require some additional vacuum inputs). The detectors at the output of $\bm{U}$ are called the boson sampling registers R. This entire setup is shown in Fig. 3.

We now run this device, post-selecting on outcomes that satisfy two properties:

• (i)

Exactly one photon is observed in each of the heralding modes H, and

• (ii)

There are exactly $n$ photons in total in the R registers.

The two properties above guarantee that, in every event accepted by the post-selection, exactly one photon was injected into each non-vacuum input of $\bm{U}$ [due to the form of Eq. (17)], and no photons were lost within $\mathcal{U}$. Therefore, the resulting conditional probability distribution is the same as an ideal boson sampling instance with single-photon inputs and interferometer $\bm{U}$.

Now note that standard boson sampling, when augmented with the power of post-selection, can perform universal quantum computation Aaronson and Arkhipov (2011). Thus, by choosing the interferometer $\bm{U}$ properly, the same is true for the device of Fig. 3. From this it immediately follows that, by a standard argument (see e.g. Bremner et al. (2010)), there can be no efficient classical algorithm to simulate the output distribution of the lossy GBS device exactly (or up to multiplicative error), otherwise the polynomial hierarchy collapses to its third level.

It is a well-accepted complexity-theoretic conjecture that the polynomial hierarchy is infinite, and so Theorem 1 can be taken as evidence that an efficient classical algorithm which exactly simulates a lossy GBS device does not exist. Note that the theorem did not require any assumptions on the strength of either losses or squeezing, so it holds for any squeezing parameter $r>0$ and any interferometer transmissivity $\eta>0$. Interestingly, if we replace $\bm{U}$ by the construction described in Brod (2015), Theorem 1 also proves that GBS (lossy or not) is hard to simulate even if the entire linear-optical sector of Fig. 3 has only five layers of (long-range) beam splitters.

Just like previous similar results Aaronson and Arkhipov (2011); Rahimi-Keshari et al. (2016), Theorem 1 is not too relevant in a realistic scenario. The requirement of simulating the output distribution exactly (or with multiplicative error) is too strict, since a realistic device with experimental imperfections is not simulating the idealized device to that precision either. At best, Theorem 1 places bounds on how far a proposed efficient classical algorithm can be extended (for example, it shows that any exact simulation of lossy GBS based on the algorithm of Rahimi-Keshari et al. (2016) can only be efficient in the presence of dark counts).

III Mapping Squeezed lossy states to squeezed thermal states

Any zero-mean single-mode Gaussian state can be decomposed into a squeezed thermal state $\rho=S(s_{\rho},\phi_{\rho})\rho_{T}^{n_{\rho}}S^{\dagger}(s_{\rho},\phi_{\rho})$ Marian and Marian (1993). Here $\rho_{T}^{n_{\rho}}$ is a thermal state with average photon number $n_{\rho}$, and its covariance matrix is given by $(2n_{\rho}+1){\bm{I}}_{2}$. $S(r,\phi)$ is the Stoler squeezing operator,

[TABLE]

Its transformation on the quadratures $x$ and $p$ is given by a symplectic matrix

[TABLE]

The covariance matrix of $\rho$ can be written as

[TABLE]

Recall that a single-mode Gaussian state is $(t)$-classical when ${\bm{V}}_{\rho}-t{\bm{I}}_{2}$ is positive definite. From Eq. (20) we know that this happens when

[TABLE]

For a lossy squeezed state $\sigma$ with ${\bm{V}}_{\sigma}=\operatorname{diag}\left\{a_{+},a_{-}\right\}$ where $a_{\pm}=\eta e^{\pm 2r}+(1-\eta)$, we can easily solve $s_{\sigma}$ and $n_{\sigma}$ by directly comparing ${\bm{V}}_{\sigma}$ with Eq. (20), which gives

[TABLE]

IV Sufficient conditions derived by using the quantum fidelity

In the main text we derive the following bound by using the quantum fidelity,

[TABLE]

where

[TABLE]

We explicitly calculate this quantity in what follows.

We start by computing the quantum fidelity between the lossy squeezed state $\sigma$ and a $(t)$-classical Gaussian state $\tau$. The fidelity between two single-mode Gaussian states is given by  Spedalieri et al. (2012)

[TABLE]

where $\Delta=\frac{1}{4}\det({\bm{V}}_{\sigma}+{\bm{V}}_{\tau})$ and $\Lambda=\frac{1}{4}(\det{\bm{V}}_{\sigma}-1)(\det{\bm{V}}_{\tau}-1)$. To optimize Eq. (26) we make use of the STS parameterization for $\sigma$ and $\tau$. The fidelity is minimized when the squeezing axes are aligned Marian and Marian (1993), so we can set $\phi_{\tau}=0$ for simplicity. Straightforward calculations then lead to

[TABLE]

Since $\tau$ is a $(t)$-classical Gaussian state, using Eq. (4) in the main text we have $s_{\tau}\leq s_{0}(\tau):=\frac{1}{2}\ln\frac{2n_{\tau}+1}{t}$. The task is then to find the point $(s^{*}_{\tau},n_{\tau}^{*})$ that maximizes the quantum fidelity subject to that constraint. Note from Eq. (26) that the fidelity monotonically decreases with $\left|s_{\sigma}-s_{\tau}\right|$. So its optimization has two regimes.

First, when $s_{\sigma}\leq s_{0}(\sigma)$, a maximum of the fidelity is reached at $s_{\tau}^{*}=s_{\sigma}$ and $n_{\tau}^{*}=n_{\sigma}$, which gives $F(\sigma^{*}_{\tau},n^{*}_{\tau})=1$. This corresponds to the case when $\sigma\in\mathcal{C}_{G}^{(t)}$, i.e., $\sigma$ itself is a $(t)$-classical Gaussian state. This regime reproduces the previous result for exact simulation of GBS.

Second, when $s_{\sigma}\geq s_{0}(\sigma)$, the fidelity is maximized at $s_{\tau}^{*}=s_{0}(\tau)$. Substitute $s_{0}(\tau)$ into Eqs. (26)-(28), we have a function of $n_{\tau}$ to optimize. It follows that its maximum is reached at $n_{\tau}^{*}=-\frac{1}{2}+\frac{1}{2}\sqrt{1+2t\sinh(2s_{c})\exp(2s_{\sigma})}$, where $s_{c}=\frac{1}{2}\ln(2n_{\sigma}+1)$. The corresponding maximum fidelity is $F(s^{*}_{\tau},n^{*}_{\tau})=\text{sech}(s_{\sigma}-\frac{1}{2}\ln\frac{2n_{\tau}+1}{t})$. Our result agrees with previous result when $t=1$ Marian et al. (2002).

Combining both regimes we write the maximum fidelity compactly as $\text{sech}\left[\Theta\left(s_{\sigma}-\frac{1}{2}\ln\frac{2n_{\sigma}+1}{t}\right)\right]$, where $\Theta$ is the ramp function. We now need to further optimize this over $t\in[1-2q_{D},1]$. It is clear that the maximum value is achieved at $t=1-2q_{D}$. Using Eqs. (22)-(23) we finally write the maximum fidelity as

[TABLE]

Plugging this into Eq. (24) we obtain our final sufficient condition for classical simulability of noisy GBS. Notice that if we set $\epsilon=0$ we recover the corresponding condition for exact simulation given in Eq. (8) in the main text, as expected.

V Optimizing the simulability condition using sandwiched Rényi relative entropy

The additive error of our approximate classical algorithm of GBS is upper bounded by using the sandwiched Rényi relative entropy, which is formally defined as follows for two quantum states $\rho_{1}$ and $\rho_{2}$ Müller-Lennert et al. (2013); Wilde et al. (2014),

[TABLE]

for $\alpha\in(0,1)\cup(1,\infty)$. Below we list some of its properties which are useful for us later:

• •

Unitary invariance: Given any unitary $U$, $D_{\alpha}(U\rho_{1}U^{\dagger}\|U\rho_{2}U^{\dagger})=D_{\alpha}(\rho_{1}\|\rho_{2})$ for $\alpha\in[\frac{1}{2},1)\cup(1,\infty)$ Müller-Lennert et al. (2013).

• •

Additivity: $D_{\alpha}(\rho_{1}\otimes w_{1}\|\rho_{2}\otimes w_{2})=D_{\alpha}(\rho_{1}\|w_{1})+D_{\alpha}(\rho_{2}\|w_{2})$ for $\alpha\in[\frac{1}{2},1)\cup(1,\infty)$ Müller-Lennert et al. (2013).

• •

Data processing: For any completely positive trace-preserving map $\mathcal{E}$, we have

[TABLE]

for $\alpha\in[\frac{1}{2},1)\cup(1,\infty)$ Frank and Lieb (2013); Beigi (2013).

• •

Monotonicity: $D_{\alpha}(\rho_{1}\|\rho_{2})\geq D_{\alpha^{\prime}}(\rho_{1}\|\rho_{2})$ for $\infty>\alpha\geq\alpha^{\prime}>0$ Müller-Lennert et al. (2013); Beigi (2013).

Properties listed above makes the sandwiched Rényi relative entropy a valid quantum distance measure. Moreover, it reduces to well-known quantum distance measures for specific values of $\alpha$. Specifically, we have Müller-Lennert et al. (2013); Wilde et al. (2014)

[TABLE]

which correspond to the logarithm of the quantum fidelity, the quantum relative entropy and the max-relative entropy, respectively.

Another essential ingredient in our deriviation of sufficient conditions is the generalized Pinsker’s inequality Van Erven and Harremos (2014):

[TABLE]

for $\alpha\in(0,1]$ and $D_{\alpha}(P\|Q)$ is the Rényi divergence between two distributions,

[TABLE]

Notice that it is only proved for for $\alpha\in(0,1]$, which is the reason why we have to restrict our optimization over $\alpha\in[\tfrac{1}{2},1)$.

Similar to our calculation in the main text, by using the aforementioned properties, we can derive the following sufficient conditions for efficient simulation of GBS:

[TABLE]

where $D^{\min}_{\alpha}(\eta,q_{D})$ is the $\alpha$-order Rènyi relative entropy minimized over all permitted $(t)$-classical Gaussian states:

[TABLE]

Since $D_{\alpha}$ is non-decreasing over $\alpha$, the l.h.s of Eq. (37) is expected to reach it’s minimum at some $\alpha^{*}$. To optimize over $\alpha$ we first try to calculate $D^{\min}_{\alpha}(\eta,q_{D})$ for fixed $\alpha$.

To facilitate our calculation we first define $Q(\sigma\|\tau)$ by

[TABLE]

From Ref. (Seshadreesan et al., 2018, Theorem 21), we have the following expression for $Q_{\alpha}(\sigma\|\tau)$ between two single-mode Gaussian states with zero mean:

[TABLE]

where

[TABLE]

From our optimization for $\alpha=1/2$, we expect the same landscape for general $\alpha$-order sanwiched Rényi relative entropy: it is minimized at

[TABLE]

This will give us a function of $n_{\tau}$ to minimize. The expression of $D_{\alpha}$ is too complicated to analytically show that above assumption is true. However, it can be verified analytically for $\alpha=1$, when the Rényi relative entropy reduces to the quantum relative entropy Wilde et al. (2017).

For fixed $\alpha$, we obtain $D_{\alpha}^{\min}(\eta,q_{D})$ by using Eq. (46) and then numerically minimizing over $n_{\tau}$. As shown in the Fig. 2 in the main text, We find out that the l.h.s of Eq. (37) is minimized at $\alpha=\frac{1}{2}$. That is when $D_{\alpha}(\sigma,\tau)=-\ln F(\sigma,\tau)$. Therefore, the bound we calculate in the main text by using quantum fidelity is the tightest one we can get. To give an explicit example, in the figure above we plot $\frac{2}{\alpha}D_{\alpha}^{\min}$ against $\alpha$ for $\eta=0.1,r=0.11,q_{D}=0$.

VI Efficient classical algorithm

As explained in the main text, we devise our efficient approximate sampling algorithm in two steps: 1) for a given noisy Gaussian state find its closest classical state and 2) sample exactly from this classical state according to the algorithm given in Rahimi-Keshari et al. (2016). Therefore, in order to present our approximate sampling algorithm, we first review the phase-space quasi-probability distribution (PQD) method used in Rahimi-Keshari et al. (2016).

Consider an $M$-mode bosonic system and let $\hat{{\bm{x}}}:=(\hat{q}_{1},\hat{p}_{1},\dots,\hat{q}_{M},\hat{p}_{M})$ denote the row vector of its quadrature operators. They satisfy the commutation relations ($\hbar=2$)

[TABLE]

Here ${\bm{I}}_{M}$ is the $M\times M$ identity matrix. The annihilation operators are given by $\hat{a}_{j}=\tfrac{1}{2}\left(\hat{q}_{j}+i\hat{p}_{j}\right)$ and $\hat{{\bm{a}}}:=(\hat{a}_{1},\ldots,\hat{a}_{M})^{T}$.

The ${\bm{t}}$-ordered PQD [$({\bm{t}})$-PQD] of an $M$-mode Hermitian operator $\rho$ is defined as

[TABLE]

where ${\bm{x}}\in\mathbb{R}^{2M}$ are the eigenvalues of quadrature $\hat{{\bm{x}}}$ and $-{\bm{I}}_{M}\leq{\bm{t}}\leq{\bm{I}}_{M}$. Here $\Phi^{({\bm{t}})}_{\rho}({\bm{\xi}})=\operatorname{Tr}[\rho D({\bm{\xi}})]\exp({\bm{\xi}}{\bm{t}}{\bm{\xi}}^{\dagger}/2)$ is the ${\bm{t}}$-ordered characteristic function and $D({\bm{\xi}})=\exp(i\hat{{\bm{x}}}^{T}{\bm{\Omega}}{\bm{\xi}})$ is the displacement operator. We slightly abuse the notation that ${\bm{t}}$ is a diagonal matrix with various order parameters on its diagonal. For ${\bm{t}}=-{\bm{I}}_{M}$, ${\bm{t}}=0$ and ${\bm{t}}={\bm{I}}_{M}$, we obtain the Husimi, Wigner and Glauber-Sudarshan functions, respectively.

In Ref. Rahimi-Keshari et al. (2016) it was shown that the probability of detecting a photon pattern can be written as the overlap between PQDs:

[TABLE]

where $W_{\rho_{\text{out}}}^{({\bm{t}})}({\bm{x}})$ is the $({\bm{t}})$-PQD of the pre-measurement state and $W_{\Pi}^{(-{\bm{t}})}({\bm{n}}|{\bm{x}})$ is the $(-{\bm{t}})$-PQD of the measurement operator. The crucial observation in Ref. Rahimi-Keshari et al. (2016) is the following: If both PQDs are positive and can be simulated efficiently for some ${\bm{t}}$, the device as a whole can be efficiently simulated by successively sampling from the chain of distributions given in Eq. (51).

For threshold detector with quantum efficiency $\eta_{D}$ and dark count rate $q_{D}$, its POVM elements are given by

[TABLE]

It was shown that their PQDs

[TABLE]

are non-negative provided $t\geq\bar{t}=1-2p_{D}/\eta_{D}$ Rahimi-Keshari et al. (2016). Since $({\bm{t}})$-classical Gaussian states have positive Gaussian $(t{\bm{I}}_{M})$-PQDs,

[TABLE]

any $(\bar{t}{\bm{I}}_{M})$-classical Gaussian states can be exactly and efficiently sampled according to  (51). Among those state, the one which is closest to the output state of a noisy GBS device is solved in the main text and with more details in Sec. IV.

We outline our approximate sampling algorithm in Alg. 1. Since the construction of the covariance matrix of a multi-mode Gaussian state, sampling from a multi-variant Gaussian distribution and sampling from a Bernoulli distribution can be done efficiently, our algorithm can be executed in polynomial time. An implementation using the Strawberry Fields library Killoran et al. (2019) can be found at 222Github repository at https://github.com/dionysos137/GBS_classicality.