Efficient simulation scheme for a class of quantum optics experiments with non-negative Wigner representation

There have been a variety of approaches to the problem of characterizing what is non-classical about quantum theory. In continuous variable quantum theory, and quantum optics in particular, the most frequently considered notions of quantumness are phrased in terms of quasi-probability distributions, such as the Wigner function and the (Glauber–Sudarshan) P-function. There is a strong tradition of considering negativity of these quasi-probability functions as some indicator of the non-classicality of a quantum state [1–4]. There are other approaches to identifying signatures of quantum theory, but with the rise of quantum information, the exponential speedup of some quantum algorithms over the best known classical algorithms have increasingly attracted attention as an important signature of quantum theory. This approach is relevant in the quantum optics setting due to the Knill–Laflamme–Milburn (KLM) model for universal quantum computation with quantum optical elements [5] and especially poignant in light of the recent work by Aaronson and Arkhipov wherein a simple non-universal linear optical system is shown to be able to perform computational tasks believed to be hard for classical computers [6]. Whether quantum computational speedups (and the boundaries between computational complexity classes) are reflected in the traditional measures of non-classicality based on negative Wigner functions was an open question answered in the affirmative, for finite-dimensional systems of qudits, by the recent work of Veitch et al [7, 8]. It is therefore natural to suspect that negative quasi-probability is intimately linked to quantum computational speedups also in the case of continuous variable quantum information processing.

Continuous variable quantum information theory provides a potentially powerful alternative to the usual discrete formalism and many of the seminal results in discrete variable quantum computation have analogues in the continuous variable setting. Perhaps the most important example is the 'continuous variable Gottesman–Knill theorem', which states that a computation restricted to the subset of quantum theory containing only Gaussian states and operations is classically efficiently simulatable [9, 10]. More concretely, unitary Gaussian quantum information is defined to be the following set of operators (see, e.g. [11]): n mode Gaussian input state; quadratic Hamiltonians; and measurements with (or without) post-selection onto Gaussian states. Bartlett et al [9, 10] have shown explicitly that there exists a classical algorithm that reproduces the output probabilities of the measurement results and executes in time that scales polynomially with the number of modes. This shows that some non-Gaussian resources are necessary to obtain an exponential speed-up with quantum optical experiments, but leaves open the question of whether they are sufficient.

Recently, Veitch et al [7, 8] have established an important connection between negative quasi-probability and quantum speed-up. Their work considers access to Clifford operations on qudits (the case of qubits is not covered by their proof) and measurements of stabilizer states and develops an efficient classical simulation of this model via a discrete analogue of the Wigner function [12]. Furthermore, they show that the onset of negativity in the discrete Wigner function can be used to identify a necessary condition for a mixed quantum state to enable an exponential speed-up through quantum computation. Their result implies, somewhat surprisingly, that there exist a class of bound universal states outside of the convex set of stabilizer states that can still be efficiently simulated and therefore also can not serve as a resource for an exponential speed-up with quantum computation.

There is a tight mathematical correspondence between the discrete and continuous Wigner representations, the Clifford/stabilizer model for qudits [13] and the Gaussian model for quantum optics considered by Bartlett et al [9, 10]. It is therefore natural to ask whether the restriction to Gaussian states in the model of Bartlett et al can be relaxed to allow more general class of initial states that have non-negative Wigner representation while still permitting an efficient classical simulation.

This work affirms an answer in the positive by showing that a large class of quantum states with positive Wigner representation exists outside the convex hull of the n-mode Gaussian states that can be efficiently simulated using a classical computer, given restrictions to quadratic Hamiltonians and Gaussian measurement. We refer to these states as non-Gaussian bound states in analogy to the discrete variable case, although no protocol has yet been proposed for distillation of non-Gaussian states using Gaussian resources, to the best of our knowledge. Our simulation result shows, in particular, that linear optical quantum devices are essentially no more computationally powerful than classical computers under such restrictions. We show this by providing an explicit classical simulation algorithm that can be used to simulate sampling the output probability distributions of the evolved initial states. As a practical application we apply our results to determine a threshold on the computational power of single-photon-added-thermal states (SPATS) [14–16] for variable efficiencies. In this sense, our work serves as both a conceptual and practical generalization of the continuous variable Gottesman–Knill theorem to a broader class of input states. This complements prior work that shows that Gaussian cluster states with finite squeezing cannot be used to perform universal quantum computation using existing proposals for measurement based quantum computing, even when local non-Gaussian measurements are permitted [17, 18]. This is despite the fact that such states can be highly entangled [17].

This paper is outlined as follows. We begin with a brief review of the Wigner function formalism and Gaussian quantum mechanics in section 2. We then provide our simulation protocol for states with positive Wigner representation in section 3. In section 4, we discuss positivity of the Wigner function as a necessary condition for quantum computation. We illustrate the bound state region via the recently studied class of limited-efficiency SPATS (LESPATS) and show that quantum efficiencies of 50% are a necessary threshold for computational speed-up. Finally, section 5 contains our conclusion and further discussion about our findings.

Wigner functions provide a natural quantum analogue of the classical phase space distribution of a dynamical system. We provide below a brief review of the properties of Wigner functions. For simplicity, we focus our attention on Wigner functions for a single particle (or equivalently a single mode) in one dimension. The generalization to higher dimensions and more particles is straightforward [19].

The Wigner representation of a state ρ is defined to be [19]

$$\begin{equation} W_{\rho}(q,p)=(2\pi)^{-1}\int_{-\infty}^{\infty}\left\langle q-y/2 \right|\rho\left|q+y/2\right\rangle\mathrm{e}^{\mathrm{i}py/\hbar}\mathrm{d}y, \end{equation} \tag{ 1 }$$

where |q〉 is a position eigenstate. The Wigner function is both positive and negative in general. However, it otherwise has many of the same properties as a classical probability density on phase space. For these reasons, the Wigner function is often referred to as a quasi-probability function. Intuitively, if we could find a bona fide joint probability distribution of non-commuting observables, then there would be no difference between quantum and classical theories. It is not surprising, then, that negativity is necessary in all possible quasi-probability representations of a quantum state [20].

The time-evolution of the Wigner function for a Hamiltonian of the form H = p²/2m + V (q) is given by [19, 21]

$$\begin{equation} \fl \frac{\partial W_{\rho}(q,p)}{\partial t} =\{H,W_\rho\} + \sum_{\ell=1}^\infty\frac{1}{(2\ell+1)!}\left(-\frac {\mathrm{i}}{2}\right)^{2\ell}\frac{\partial^{2\ell+1}V(q)}{\partial q^{2\ell+1}}\frac{\partial^{2\ell+1}W_{\rho}(q,p)}{\partial {p}^{2\ell+1}}, \end{equation} \tag{ 2 }$$

where {·,·} is the Poisson bracket, which governs classical Hamiltonian equations of motions. This result is important because it states that the time-evolution of W_ρ(q,p) is given by Liouville's equation, plus a quantum correction. The quantum correction is zero for the case of the quadratic Hamiltonian:

$$\begin{equation} \frac{\partial W_{\rho}(q,p)}{\partial t} = \{H,W_\rho\}. \end{equation} \tag{ 3 }$$

Hence the evolution equation agrees precisely with the classical predictions. This observation will be vital for our simulation algorithm because the Hamiltonians permitted by linear optics are quadratic (harmonic oscillators), which (along with our non-negativity assumption) implies that the evolution of W_ρ(q,p) can be simulated using an ensemble of classical trajectories.

This discussion above implies that a Wigner function that is initially classical, meaning that it is non-negative and hence interpretable as a probability density function (pdf), will remain classical under the action of a quadratic Hamiltonian. In this context is therefore useful to determine the conditions under which a Wigner function is non-negative as this gives a practically relevant boundary between quantum and classical states. Hudson's theorem [22] was the first attempt to characterize the positive Wigner functions and it was later generalized to the following [23]. Let ψ be a pure quantum state of n oscillators (modes). Then its Wigner function is positive if and only if

$$\begin{equation} \psi(\skew4\vec{Q})=\mathrm{e}^{-\frac{1}{2}(\skew4\vec{Q}\cdot A\skew4\vec{Q}\,+\,B\cdot\skew4\vec{Q}\,+\,c)}, \end{equation} \tag{ 4 }$$

where A is an n × n Hermitian matrix, B is an n-dimensional complex vector and c is a normalization constant. In quantum optics terminology, these are either coherent states or squeezed states. That is, plugging these states into the definition of the Wigner function yields multivariate Gaussian distributions in phase space. Convex combinations of these states (incoherent mixtures of them) also have positive Wigner function since the mapping is linear. Early on, these were incorrectly conjectured to be the only such mixed states with positive Wigner function. The question of mixed states was given a full treatment in [24] and latter in [25]. Both references independently found that a theorem in classical probability attributed to Bochner [26] and generalization thereof can be used to characterize both the valid Wigner functions and the subset of positive ones. What was shown is that there exist a large class of states with positive Wigner function that are not convex combinations of Gaussian states. So far, these states have received little attention. In section 4, we show that such states are more than a mathematical curiosity; they arise naturally in quantum optics.

Gaussian measurements are also easily modeled in the Wigner representation. Recall that for non-negative states the Wigner function picture allows us to represent the system as a probability density over underlying physical states in phase space, u_f. Gaussian measurements in this picture are also modeled as probability densities for outcomes k, conditioned on the value of the underlying physical state u_f. Specifically, consider the case of measurement ${\cal M}$ of a Gaussian state G with covariance matrix $V_{{\cal M}}$. The POVM representation of this measurement is [27]

$$\begin{equation} {\cal M}(V_{{\cal M}})=\{|{G(\boldsymbol{k},V_{{\cal M}}}\rangle\langle{G(\boldsymbol{k},V_{{\cal M}})}|:\ \boldsymbol{k}\in\mathbb{R}^{2n}\}, \end{equation} \tag{ 5 }$$

which selects a Gaussian state with mean k and covariance $V_{{\cal M}}$ from all possible Gaussian states with mean k and covariance $V_{{\cal M}}$. In the Wigner function picture the representation of this measurement is

$$\begin{equation} M_{{\cal M}(V_{{\cal M}})}(\boldsymbol{k}|\boldsymbol{u}_f)=\frac{1}{{\cal X}}\exp(-(\boldsymbol{k}-\boldsymbol{u}_f)^{\mathrm{T}}V_{{\cal M}}^{-1}(\boldsymbol{k}-\boldsymbol{u}_f)), \end{equation} \tag{ 6 }$$

where ${\cal X}$ is the normalization constant. We introduce the notation in (6) to emphasize the difference between the representations of measurements and states. The interpretation of this equation is that if the system is actually at the point u_f the effect of a measurement will be to produce an outcome k according to the probability density $M_{{\cal M}(V_{{\cal M}})}(\boldsymbol {k}|\boldsymbol {u}_f)$. Using this equation and the law of total probability we can find the probability density of measurement outcomes k for the measurement ${\cal M}(V_{{\cal M}})$ on a quantum state with Wigner representation W_ρ(k):

$$\begin{equation} p({\boldsymbol{k}}|{\cal M}(V_{{\cal M}}),\rho)=\int_{\boldsymbol{u}}W_\rho(\boldsymbol{u})M_{{\cal M}(V_{{\cal M}})}({\boldsymbol{k}}|{\boldsymbol{u}})\mathrm{d}{\boldsymbol{u}}. \end{equation} \tag{ 7 }$$

Of course, this agrees with the probability assigned by the Born rule

$$\begin{equation} p(\boldsymbol{k}|{\cal M}(V_{{\cal M}}),\rho)=\mathrm{Tr}(|{G(\boldsymbol{k},V_{{\cal M}}}\rangle \langle{G(\boldsymbol{k},V_{{\cal M}})}|\rho). \end{equation} \tag{ 8 }$$

The simulation algorithm that we propose uses none of the special properties of Gaussian measurements other than the fact that they have a positive Wigner representation and that we can efficiently draw samples from a Gaussian distribution. This means that our results will apply to any measurements that satisfy these properties. We focus our attention on Gaussian measurements rather than these more general measurements because of their simplicity and physical relevance.

At first glance it may seem that a simulation algorithm for linear optics may be difficult owing solely to exponential size of the Hilbert space dimension that is generated by the evolution. We overcome this problem by exploiting the fact that our Hamiltonians are quadratic in p and q, which implies that the evolution of the Wigner function follows the Liouville equation as shown in equation (3). Since the Liouville equation preserves non-negativity and probability mass, the Wigner function will remain a classical distribution throughout the evolution. This allows us to model evolution of the Wigner function using an ensemble of classical trajectories, each of which can be efficiently simulated. The resulting trajectories can be used to efficiently draw samples from the final distribution of measurement outcomes prescribed by the Born rule without needing to know the final quantum state.

It is critical to understand that we are not simulating the evolution of the quantum state, rather we are simulating measurement outcomes that occur at the end of a linear-optical protocol; this is exactly analogous to the difference between knowing a probability distribution and being able to sample from it. In particular, the ability to efficiently draw samples does not imply the ability to efficiently learn the underlying distribution because the dimension of the probability distribution on n modes is exponentially large.

We show in this section that this simulation strategy can be used to efficiently sample from the output of the following class of quantum algorithms:

Algorithm Class 1. Family of efficiently simulatable quantum algorithms. Input. Number of modes n, an initial n mode quantum state $\rho =\rho _1\otimes \ldots \otimes \rho _{n}$ where each ρj has positive Wigner representation Wρj(u), a description of a linear optical transformation UT,x which is parameterized by $T\in \mathbb {R}^{2n\times 2n}$ and $\boldsymbol {x}\in \mathbb {R}^{2n}$. Output. A string of measurement outcomes k sampled according to the probability density p(Kquant = k) determined by the Born rule.
1.	Apply the linear optical transformation UT,x.
2.	Perform the separable Gaussian measurement ${\cal M}(V_{{\cal M}})={\cal M}(V_{{\cal M},1})\otimes {\cal M}(V_{{\cal M},2})\otimes \ldots$ $\otimes {\cal M}(V_{{\cal M},n})$ where we follow the naming convention of equation (5) and the tensor product is understood to mean that the POVM elements of ${\cal M}(V_{{\cal M}})={\cal M}(V_{{\cal M},1})\otimes \cdots \otimes {\cal M}(V_{{\cal M},n})$ are tensor product combinations of the POVM elements of ${\cal M}(V_{{\cal M}})$ in the obvious way.
3.	Return the measurement outcome $\boldsymbol {k}=(\boldsymbol {k}_{1},\boldsymbol {k}_{2},\ldots ,\boldsymbol {k}_{n})\in \mathbb {R}^{2n}$ corresponding to the mean of a Gaussian POVM element.

Here we conceive of any quantum algorithm in this class as a way of sampling outcome strings k distributed according to the probability densities given by the Born rule. We label the corresponding random variable K_quant. Here we are not simulating the evolution of the Wigner distribution, which would be equivalent to simulating the full quantum state. Rather, we show that there is a corresponding classical algorithm that produces outcome strings k with (very nearly) the same distribution those from algorithm class 1.

Using intuition similar to that in [7], we note that quantum algorithms in class 1, can be simulated using the following classical algorithm, provided access to classical resources with infinite numerical precision and a blackbox function that can be used to draw samples from W_ρj(u) for j = 1,...,n. We will later provide an algorithm that does not require infinite precision, but we provide the infinite precision algorithm first because it conveys the necessary intuition without focusing on the technical issues that arise when discretizing the distributions.

Algorithm 1. Infinite precision classical simulation algorithm for algorithms in class 1. Input. As algorithms in class 1, except ρ is not provided. Output. A string of measurement outcomes k sampled according to the probability density p(K(2)class = k).
1.	Sample $\boldsymbol {u}\in \mathbb {R}^{2n}$ according to the distribution Wρ(u) = Wρ1(u1)···Wρn(un) by drawing a sample independently from each mode using the blackbox function.
2.	Apply the affine transformation $\tilde {\boldsymbol {u}}=T\boldsymbol {u}+\boldsymbol {x}$ corresponding to the linear optical transformation UT,x to the sampled phase space point u. This transformation is an affine mapping due to Louiville's theorem.
3.	Return the outcome string $\boldsymbol {k}=(\boldsymbol {k}_{1},\boldsymbol {k}_{2},\ldots ,\boldsymbol {k}_{n})\in \mathbb {R}^{2n}$ from the distribution $M_{{\cal M}(V_{{\cal M}})}(\boldsymbol {k}|\tilde {\boldsymbol {u}})=M_{{\cal M}(V_{{\cal M},1})}(\boldsymbol {k}_{1}|\tilde {\boldsymbol {u}}_{1})\cdots M_{{\cal M}(V_{{\cal M},n})}(\boldsymbol {k}_{n}|\tilde {\boldsymbol {u}}_{n})$, where $M_{{\cal M}(V_{{\cal M},j})}(\boldsymbol {k}_{i}|\tilde {\boldsymbol {u}}_{i})$ is given as in equation (6),
	$\begin{eqnarray*} &&\displaystyle M_{{\cal M}(V_{{\cal M},j})}(\boldsymbol{k}_{i}|\tilde{\boldsymbol{u}}_{i})=\frac{1}{{\cal X}_j}\exp(-(\boldsymbol{k}_{i}-\tilde{\boldsymbol{u}}_{i})^{\mathrm{T}}V_{{\cal M},j}^{-1}(\boldsymbol{k}_{i}-\tilde{\boldsymbol{u}}_{i})). \end{eqnarray*}$

The intuition behind this class of algorithms is to use the classical phase space model afforded to us by the non-negative Wigner functions and quadratic evolutions in order to turn the quantum problem into one that can be efficiently simulated by a classical computer. In this context we can think of the quantum system as actually being definitely at some point $\boldsymbol {u}\in \mathbb {R}^{2n}$ that is unknown to the user of the algorithm. The point then moves under a fully deterministic classical evolution and measurement on each register amounts to picking a point k from a normal distribution centered at the location of the system. Each classical algorithm samples outcomes k according the probability density

$$\begin{eqnarray*} &&p(\boldsymbol{K}_{\mathrm{class}}^{(2)}=\boldsymbol{k})=\int_{\boldsymbol{u}}W_\rho(\boldsymbol{u})M_{{\cal M}(V_{{\cal M}})}(\boldsymbol{k}|\boldsymbol{u})\,\mathrm{d}\boldsymbol{u}, \end{eqnarray*}$$

which agrees with the density given by the Born rule. Thus the outcomes K⁽²⁾_class from the classical simulator are distributed in exactly the same way as K_quant, which are outcomes drawn from the actual quantum system.

Unfortunately, algorithms similar to 1 cannot be executed precisely on digital computers and instead would require an analogue computer (often referred to as a real computer). If physical, such computers would have unrealistic computational powers such as being able to solve NP-complete problems in polynomial time [28] and would also violate the holographic principle [29]. For these reasons, we need to discretize 1 in order to assess the cost of simulating linear optics on realistic classical computers. The major technical difference between the continuous variable case and the discrete case considered in [7] involves showing that finite-precision errors can be made negligible with efficient overhead costs given a set of reasonable assumptions about the input states, dynamics and measurements.

Since infinite precision is requried in the continuous variable setting, in order to specify a quantum state we must make some finite precision truncation. To this end we shall assume access to a family of oracles ${\cal W}_{\rho _{j},\eta }(l,m)$ that takes integers l,m and returns a value satisfying $|{\cal W}_{\rho _{j},\eta }(l,m)-W_{\rho _{j}}(\boldsymbol {\mu }_{\rho _j}+(l,m)\delta )|<\eta$, for discretization length δ and discretization error η > 0. That is, each oracle queries the Wigner function at points on a grid centered at the mean of the distribution. This is a weak assumption as it does not require us to even know the state we are simulating. Using this resource the algorithm can be written as:

Algorithm 2. Finite-precision classical simulation algorithm for quantum algorithms in class 1. Input. As algorithm 1, but also require δ, a discretization length for the input, |2|, a bound for the numerical error involved in applying the affine transformation, Γ, a discretization length for the output, $|{\cal A}|$, the area truncated square region of phase space that the simulator considers and μρj, the mean of the Wigner distribution of the quantum state on each mode j. We require $\sqrt {|{\cal A}|}$ to be an odd integer multiple of δ and Γ to be an odd integer multiple of δ. Output. A string of measurement outcomes k sampled according to $\mathrm {Pr}(\boldsymbol {K}_{\mathrm {class}}= \boldsymbol {k})\equiv \mathrm {Pr}_{\mathrm {sim}}(\boldsymbol {k})$.
1.		For each j = 1,...,n execute (a) through (d). (This step approximates sampling a point from phase space.)
	(a)	For each integer $l,m \in \left[-\frac {\sqrt {|{\cal A}|}}{\delta },\frac {\sqrt {|{\cal A}|}}{\delta }\right]$ set $\mathrm {Pr}_{\mathrm {sim},\rho_{j}}(l,m)={\cal W}_{\rho _{j},\eta }(l,m)\cdot \delta ^{2}.$
		(The phase space is truncated to a size $|{\cal A}|$ and discretized into boxes of size δ. This step sets a pdf over the centers of the boxes.)
	(b)	For each (l,m) set $\mathrm {Pr}_{\mathrm {sim},\rho_{j}}(l,m)=\mathrm {Pr}_{\mathrm {sim},\rho_{j}}(l,m)/\sum _{l,m}\mathrm {Pr}_{\mathrm {{sim}},\rho _{j}}(l,m).$ (This step normalizes the pdf.)
	(c)	Draw a sample (l,m) from the pdf $\mathrm {Pr}_{\mathrm {sim},\rho_{j}}(l,m).$
	(d)	Set uj = μj + (l,m)δ
2.		Set ${\tilde {\boldsymbol {u}}}=T\boldsymbol {u}+\boldsymbol {x}$ using enough digits of precision such that the numerical error is at most |2|, where $\boldsymbol {u}\equiv \boldsymbol {u}_{1}\oplus \ldots \oplus \boldsymbol {u}_{n}$ and ${\tilde {\boldsymbol {u}}}\equiv {\tilde {\boldsymbol {u}}}_{1}\oplus \ldots \oplus {\tilde {\boldsymbol {u}}}_{n}$. (This step corresponds to updating the sampled state according to the linear optical transformation.)
3.		For each j = 1,...,n execute (a) through (e). (This step is to simulate drawing a measurement outcome from the Gaussian measurement distribution centered at ${\tilde {\boldsymbol {u}}}$.)
	(a)	For each integer $l,m \in \left[-\frac {\sqrt {|{\cal A}|}}{\delta },\frac {\sqrt {|{\cal A}|}}{\delta }\right]$ set $\mathrm {Pr}_{\mathrm {sim},{\cal M}(V_{{\cal M},j})}(l,m):=\exp (-\delta (l,m)^{\mathrm {T}}V_{{\cal M},j}^{-1}$ $\delta (l,m))\cdot \delta ^{2}.$
		(The outcome space is truncated to a size $|{\cal A}|$ and discretized into boxes of size δ. This sets a pdf over the centers of the boxes.)
	(b)	Set $\mathrm {Pr}_{\mathrm {sim},{\cal M}(V_{{\cal M},j})}(l,m)=\mathrm {Pr}_{\mathrm {sim},{\cal M}(V_{{\cal M},j})}(l,m)/\sum _{l,m}\mathrm {Pr}_{\mathrm {sim},{\cal M}(V_{{\cal M},j})}(l,m).$ (This step normalizes the pdf.)
	(c)	Draw a sample (l,m) from the pdf $\mathrm {Pr}_{\mathrm {sim},{\cal M}(V_{{\cal M},j})}(l,m)$.
	(d)	Find integers (r,s) such that |δ(l,m) − Γ(r,s))|∞ ⩽ Γ/2.
		(This just amounts to rebinning the outcome distribution into hypercubes of sidelength Γ; this introduces no errors if Γ is an integer multiple of δ.)
	(e)	Set measurement outcome $\boldsymbol {k}_{j}=\tilde {\boldsymbol {u}}_j+\Gamma (r,s)$.
4.		Return measurement outcome $\boldsymbol {k}\equiv \boldsymbol {k}_{1}\oplus \ldots \oplus \boldsymbol {k}_{n}$.

Our simulation protocol can necessarily only sample from a discrete distribution so we must introduce some notion of how a discrete distribution can be close to the continuous probability density given by the Born rule. The most natural way to do this is to discretize the outcome distribution into boxes of side length Γ according to,

Definition 1. Let ${\cal N}_{\Gamma }(\boldsymbol {k})\subset \mathbb {R}^{2n}$ be a hypercube in outcome space with side length Γ centered at the point k. We define the Γ discretization of the quantum outcome distribution to be $\mathrm {Pr}(\boldsymbol {K}_{\mathrm {quant},{\Gamma }}=\boldsymbol {k})\equiv \mathrm {Pr}_{\mathrm {quant},{\cal M}(V_{\cal M})}(\boldsymbol {k})\equiv \int _{{\cal N}_{\Gamma }(\boldsymbol {k})}p(\boldsymbol {K}_{\mathrm {quant}}=\tilde {\boldsymbol {k}})\,\mathrm {d}\tilde {\boldsymbol {k}}$.

We can now fix a Γ according to our operational requirements for the simulation and ask how well a simulation protocol samples from this distribution. Notice this is an unavoidable consequence of trying to approximate a continuous quantity with a discrete system. With this in hand we can give a precise classical simulation protocol by discretizing our native algorithm, which results in the family of protocols described in algorithm 2.

It now easy to see that both the cost of the simulation and the error in our sampling will be a function of the discretization parameters δ(n,) and $|{\cal A}(n,\epsilon )|$. If we can pick these parameters such that for fixed error our simulation scheme scales as poly(n) then the simulation is efficient. Concretely,

Definition 2. Let $\mathrm {Pr}_{\mathrm {sim},{\cal M}(V_{\cal M})}(\boldsymbol {k})$ and $\mathrm {Pr}_{\mathrm {quant} {\cal M}(V_{\cal M})}(\boldsymbol {k})$ be the simulated and actual probabilities of obtaining a measured value that is inside a hypercube of volume Γ²ⁿ centered at a point in phase space $\boldsymbol {k}\in \mathbb {R}^{2n}$ for the separable Gaussian measurement ${\cal M}(V_{{\cal M}})$. A simulation algorithm is efficient if for inputs n,Γ and there exists a choice of $|{\cal A}|$ and δ such that:

• 1.  
  The one-norm distance between the Γ-discretized quantum distribution and the simulator distribution is at most ,

  $$\begin{eqnarray*} &&|\mathrm{Pr}_{\mathrm{quant},{\cal M}(V_{\cal M})}(\boldsymbol{k})-\mathrm{Pr}_{\mathrm{sim},{\cal M}(V_{\cal M})}(\boldsymbol{k})|_{1}\leqslant\epsilon, \end{eqnarray*}$$

  where we take $\mathrm {Pr}_{\mathrm {sim}, {\cal M}(V_{\cal M})}(\boldsymbol {k})=0$ whenever k is outside the domain of $\mathrm {Pr}_{\mathrm {sim},{\cal M}(V_{\cal M})}(\boldsymbol {k})$ defined by algorithm 2.
• 2.  
  The computational complexity of the simulation scales as poly(n/).

We now can show that the simulation of sufficiently smooth separable positive Wigner functions under linear optical operations and Gaussian measurements is efficient. This result is formally stated in the following theorem and proof is given in the appendix.

Theorem 1. The output of algorithm 2 satisfies $|\mathrm {Pr}_{\mathrm {quant},{\cal M}(V_{\cal M})}(\boldsymbol {k})-\mathrm {Pr}_{\mathrm {sim},{\cal M}(V_{\cal M})}(\boldsymbol {k})|_{1}\leqslant \epsilon$ for input $(\{\boldsymbol {\mu }_{\rho _j}\},\{V_{\rho _j}\},\{V_{{\cal M},j}\},T,x,\delta ,|\boldsymbol {\epsilon }_2|,\Gamma ,|{\cal A}|)$ if

• 1.  
  n, $\max _j\{|\boldsymbol {\mu }_{\rho _j}|,\|V_{\rho _j},\|V_{{\cal M},j}\|\},\|T\|$ and ||x|| are bounded,
• 2.  
  there exist finite β, Λ such that $|\nabla W_{\rho }(\boldsymbol {u})|\leqslant n\beta /|{\cal A}|^{n}$ and $\left |\nabla _{\boldsymbol {k}}M_{{\cal M}(V_{{\cal M}})}(\boldsymbol {k}|{\boldsymbol {u}})\right |\leqslant n\Lambda /|{\cal A}|^{n}$ for $\boldsymbol {k}\in \mathbb {R}^{2n},\ \boldsymbol {u}\in \mathbb {R}^{2n}$,
• 3.  
  $\delta \leqslant \min \left \{\frac {\epsilon }{16\left [\left (1+\|T\|\right )\Lambda +\beta \right ]n\sqrt {2n}},\Gamma \right \}$,
• 4.  
   < 1, $|\boldsymbol {\epsilon }_2|\leqslant \|T\|\delta \sqrt {\frac {n}{2}}$,
• 5.  
  $|{\cal A}|\geqslant 16n\max _{i,j} ( [V_{\rho _i} ]_{11}+ [V_{\rho _i} ]_{22}+ [V_{{\cal M},j} ]_{11}+ [V_{{\cal M},j} ]_{22} )/\epsilon$,
• 6.  
  The finite precision error from each oracle ${\cal W}_{\rho _j,\eta }$ satisfies $\eta \leqslant \frac {\epsilon }{8n|{\cal A}|^{n}}$.

Furthermore, if we assign unit cost to evaluations of ${\cal W}_{\rho _j,\eta }$ and Gaussian functions unit cost then the computational complexity of the algorithm is $O (n^{5} (\max _i\|V_{\rho _i}\|+\max _j\|V_{{\cal M},j}\| ) [\Lambda ^2\|T\|^2+\beta ^2 ]/\epsilon ^3 ),$ which implies efficiency.

The key insight of this theorem is that the assumption of Gaussian preparations made in the continuous variable Gottesman–Knill theorem can be relaxed [9, 10], and that a much wider class of quantum dynamics can be efficiently simulated than previously thought. Indeed, although we stated the algorithm only for product state inputs and product measurements we can now see that this restriction was unnecessary. In fact, our simulation scheme works for any positively represented input and measurement as long as it is possible to efficiently sample from the corresponding distribution. The product assumption is a sufficient but not necessary condition for this efficient sampling. We also note that the algorithm requires us to know the mean and covariance matrices of the distributions, which might be hard to compute analytically. However, since we already require efficient sampling we may appeal to Monte Carlo estimation protocols to determine these quantities within acceptable error tolerances. This extension of the continuous variable Gottesman–Knill theorem places much stronger limitations on the input states that can be used for continuous variable quantum computation and underscores the significance of negativity in the Wigner function as a resource for quantum computation (in analogy to recent results for discrete systems). In particular, we will show that theorem 1 places a minimum efficiency for a class of photonic thermal states beyond which the states cannot be used as a resource for linear optical quantum computation with Gaussian measurements.

The debate over the 'correct' definition of classicality for quantum states of light has been a long and, at times, fierce one. One of the most common notions of classicality is whether a state can be represented as a convex combination of Gaussian states. Here we have shown that, in the context of computational power, such a condition is superseded by the condition of positivity of the Wigner function. In this section we give a concrete example of an interesting class of states which are not Gaussian but which have positive Wigner representation and thereby admit an efficient classical simulation.

We consider the experimentally accessible class of states called SPATS [14–16]:

$$\begin{eqnarray*} &&\rho(\overline{n})=\frac{1}{\overline{n}(\overline{n}+1)}\sum_{n=0}^{\infty}\left(\frac{\overline{n}}{\overline{n}+1}\right)^{n}n|n\rangle\!\langle n|, \end{eqnarray*}$$

where $\overline {n}$ is the mean photon number—given by, in terms of temperature T, the Planck distribution $\overline {n}=1/(\exp(1/T)-1)$. It is known that all states in this class are outside the convex hull of Gaussian states and have negative Wigner function for finite temperatures.

Under experimentally realistic conditions we must consider states subjected to losses. In general, losses can be modeled as an interaction with a vacuum mode at a beam-splitter with transmittance η, also called the quantum efficiency. The loss rate is then 1 − η. The Wigner function of the SPATS after this channel, which we call LESPATS is [30]

$$\begin{eqnarray*} &&\fl W_{\rho(\overline{n},\eta)}(q,p)=\frac{2}{\pi}\frac{1+2\eta[\overline{n}+2(\overline{n}+1)(q^{2}+p^{2})-2\overline{n}\eta-1]}{(1+2\overline{n}\eta)^{3}}\exp\kern-1pt\left(-\frac{2(q^{2}+p^{2})}{1+2\overline{n}\eta}\right). \end{eqnarray*}$$

Note that the most negative value of the LESPATS is at (q,p) = (0,0) for all η and $\overline {n}$. Thus, we consider the quantity

$$\begin{eqnarray*} &&W_{\rho(\overline{n},\eta)}(0,0)=\frac{2}{\pi}\frac{1+2\eta(\overline{n}-2\overline{n}\eta-1)}{(1+2\overline{n}\eta)^{3}}. \end{eqnarray*}$$

By inspection, we can see that for efficiencies of η ⩽ 0.5, the Wigner function of the LESPATS is positive $W_{\rho (\overline {n},\eta )}(q,p)>0$. Thus, efficiencies of η > 0.5 are necessary for quantum computational speed-up with SPATS.

Note however that for η ⩽ 0.5, the LESPATS are not inside the convex hull of Gaussian states. To see this most clearly, we require a different phase space distribution. The Glauber P-function (see e.g. [31]) is defined as

$$\begin{eqnarray*} &&\rho=\int\!\!\int P_{\rho}(q,p)|\alpha\rangle\langle\alpha|\ \mathrm{d}q\,\mathrm{d}p, \end{eqnarray*}$$

where $\left |\alpha \right \rangle$ are the coherent states, which are vacuum states that have been displaced in phase space (symmetric Gaussian states). Note that if P_ρ is a probability distribution then ρ is in the convex hull of coherent states. The P-function of the LESPATS is [14]

$$\begin{eqnarray*} &&P_{\rho(\overline{n},\eta)}(q,p)=\frac{1}{\pi\overline{n}^{3}\eta}\left[(\overline{n}+1)\frac{q^{2}+p^{2}}{\eta}-\overline{n}\right]\exp\kern-1pt\left(-\frac{q^{2}+p^{2}}{\overline{n}\eta}\right). \end{eqnarray*}$$

Note that this function is negative for all allowed values of η and $\overline {n}$. Thus, the LESPATS are always outsides the convex hull of Gaussian states but are bound universal states [7] for η ⩽ 0.5.

To illustrate this, we compare the negativity of the Wigner function with the distance to the convex hull of Gaussian states. The distance we use is based on fidelity, which can be computed using phase space distributions as

$$\begin{eqnarray*} &&\fl F(\rho(\overline{n},\eta),|0\rangle \langle0\rangle)=\mathrm{Tr}[\rho(\overline{n},\eta)|0\rangle \langle0\rangle]=\pi\int\!\!\int P_{\rho(\overline{n},\eta)}(q,p)Q_{|0\rangle \langle0\rangle}(q,p)\,\mathrm{d}q\,\mathrm{d}p, \end{eqnarray*}$$

where the Q-function

$$\begin{eqnarray*} &&Q_{|0\rangle \langle0\rangle}(q,p)=\frac{1}{\pi}\exp(-(q^{2}+p^{2})) \end{eqnarray*}$$

is dual to the P-function⁵. Using a standard table of Gaussian integrals we find

$$\begin{eqnarray*} &&F(\rho(\overline{n},\eta),|0\rangle \langle0\rangle)=\frac{1-\eta}{(1+\overline{n}\eta)^{2}}. \end{eqnarray*}$$

This effect is demonstrated in figure 1. Note that, for any state, a quantum efficiency of η ⩽ 0.5 is sufficient to ensure membership of the convex hull of states that have positive Wigner representation. This effect is mirrored in the discrete case [7], where depolarizing noise of 50% is sufficient to ensure membership of the polytope of states with positive discrete Wigner function, when starting from any qudit state.

Zoom In Zoom Out Reset image size

We have developed an efficient classical simulation scheme for Gaussian quantum computations utilizing separable initial preparations with positive Wigner function. Because there exist states with positive Wigner function that lie outside the convex hull of Gaussian states, we have identified a large class of bound states: states that cannot be prepared using Gaussian operations, yet do not permit universal quantum computation. We illustrated this class using the example of SPATS, showing that quantum efficiencies of 50% are necessary for quantum computation.

The negativity of the Wigner function has long been used as a qualitative indicator of quantumness. By showing that negative Wigner representation has a clear operational signficance as a necessary resource for comptuation our work lends credence to efforts to extend negativity to a quantitative notion of quantumness. In terms of the Wigner function, the volume of the negative parts of the represented quantum state has been suggested as the appropriate measure of quantumness [33]. The distance (in some preferred norm on the space of Hermitian operators) to the convex subset of positive Wigner functions was suggested to quantify quantumness in [34].

Bartlett and Sanders [35] nicely summarized what was known at the time about efficient simulation of continuous variable quantum computation. The table presented there is reproduced below in table 1 with some more recent results. The field began with Lloyd and Braunstein's observation that nonlinear optical processes are sufficient for universal continuous variable quantum computation. Later, it was shown for discrete variable encodings that linear optics is sufficient provided photon counting measurements are available [5, 36]. The continuous variable analogue of the measurement-based model shows that preparation of single photon state preparation is also sufficient [37].

More recently, the result of Aaronson and Arkhipov [6] shows that preparing and measuring single photon states (without post-selection) is equivalent to a sampling problem that is thought to be hard classically—but it still manages to (probably) not be universal for quantum computation. It is possible that the Aaronson and Arkhipov model may be intermediately between classically efficiently simulatable and universal for quantum computation. Another suspected model of this type is the 'one-clean-qubit' model of Knill and Laflamme [38]. The key point for this latter model is that uses highly mixed states. Mixed states have not been given full consideration for continuous variable quantum computation. Here we have shown, via the Wigner phase space formalism and independent of purity, negative representation is necessary for universal quantum computation. Moreover, any computation that uses states possessing a positive Wigner function is classically efficiently simulatable. It would be quite interesting if this condition turned out also to be sufficient as this would provide a sharp boundary between quantum and classical systems with regard to their computational power.

We thank Earl Campbell and Christian Weedbrook for helpful comments. The authors acknowledge financial support from the Government of Canada through NSERC, CIFAR, USARO-DTO. After completion of this work, we were made aware of [39], who obtain a similar result for a different (namely, local) set of dynamical transformations.

Proof of theorem 1. Our goal is to show that with this choice of discretization parameters δ and $|{\cal A}|$ the simulation algorithms outlined in algorithm 2 require $O\!\left (\mathrm {poly}(n/\epsilon )\right )$ resources and result in error at most . Since we bin the data at the end of the protocol into hypercubes of volume Γ²ⁿ ⩾ δ²ⁿ and Γ is an integer multiple of δ, no error is introduced by first binning the outcomes into hypercubes of volume δ²ⁿ because every hypercube of volume δ²ⁿ is contained in precisely one hypercube of volume Γ²ⁿ. We therefore may take the quantum distribution to be binned into hypercubes of side length δ (i.e. Γ = δ) without loss of generality.

We start by analyzing the error. Following scheme outlined above we can decompose this into four broad parts:

• 1.  
  The error introduced by the use of finite precision output of ${\cal W}_{\rho _j,\eta }$.
• 2.  
  The error introduced by truncating the sampling distribution over phase space.
• 3.  
  The error introduced by discretizing this truncated distribution.
• 4.  
  The error introduced by truncating the outcome distribution.
• 5.  
  The error introduced by discretizing the outcome distribution.

Denoting the region in phase space that the initial states are confined to as ${\cal A}_{\rho }={\cal A}_{\rho _1}\otimes \cdots \otimes {\cal A}_{\rho _n}$ and the region that the observations are confined to as ${\cal A}_{{\cal M}(V_{{\cal M}})}={\cal A}_{{\cal M}(V_{{\cal M},1})}\otimes \cdots \otimes {\cal A}_{{\cal M}(V_{{\cal M},n})}$ (where $|{\cal A}_{\rho _i}|=|{\cal A}_{{\cal M}(V_{{\cal M},i})}|=|{\cal A}|$), we can use the triangle inequality to express these errors as

$$\begin{eqnarray} &&\fl |\mathrm{Pr}_{\mathrm{sim},{\cal M}(V_{\cal M})}(\boldsymbol{k})-\mathrm{Pr}_{\mathrm{quant},{\cal M}(V_{\cal M})}(\boldsymbol{k})|_{1} \leqslant~\left|\mathrm{Pr}_{\mathrm{sim},{\cal M}(V_{\cal M})}(\boldsymbol{k})-\mathrm{Pr}_{\mathrm{trunc-quant},{\cal M}(V_{\cal M})}(\boldsymbol{k})\right|_1\nonumber\\ &&+\mathrm{Pr}_{\mathrm{quant}}(\boldsymbol{u}\notin {\cal A}_{\rho})+\mathrm{Pr}_{\mathrm{quant}}(\boldsymbol{k}\notin {\cal A}_{{\cal M}(V_{{\cal M}})}), \end{eqnarray} \noindent \tag{ A.1 }$$

which says that the total error is at most the distance between the truncated distributions plus the probability that a point is sampled, or measured, outside the truncated region. In other words, the one–norm error introduced by truncating is, even in the most pathological case conceivable, the sum of the probabilities of sampling an initial trajectory outside of ${\cal A}_\rho$ and measuring an outcome outside of ${\cal A}_{{\cal M}(V_{{\cal M},i})}$.

It then follows that $|\mathrm {Pr}_{\mathrm {sim},{\cal M}(V_{\cal M})}(\boldsymbol {k})-\mathrm {Pr}_{\mathrm {quant},{\cal M}(V_{\cal M})}(\boldsymbol {k})|_{1} \leqslant \epsilon /2$ if

$$\begin{equation} \mathrm{Pr}_{\mathrm{quant}}(\boldsymbol{u}\notin{\cal A}_{\rho})+\mathrm{Pr}_{\mathrm{quant}}(\boldsymbol{k}\notin{\cal A}_{{\cal M}(V_{{\cal M}})})\leqslant \epsilon/4, \end{equation} \tag{ A.2 }$$

$$\begin{equation} \left|\mathrm{Pr}_{\mathrm{sim},{\cal M}(V_{\cal M})}(\boldsymbol{k})-\mathrm{Pr}_{\mathrm{trunc-quant},{\cal M}(V_{\cal M})}(\boldsymbol{k})\right|_1\leqslant \epsilon/4. \end{equation} \noindent \tag{ A.3 }$$

We will first demonstrate that (A.2) is an immediate consequence of assumption 5. We then will show that (A.3) is satisfied if assumption 3 holds.

We begin by bounding the truncation error. For the ith mode,

$$\begin{eqnarray*} \mathrm{Pr}_{\mathrm{quant}}\left(\boldsymbol{u}_{i}\notin{\cal A}_{\rho_i}\right) & = \mathrm{Pr}\left(|q_{i}-\mu_{Q}|>\frac{\sqrt{|{\cal A}|}}{2},|p_{i}-\mu_{P}|>\frac{\sqrt{|{\cal A}|}}{2}\right)\\ & \leqslant \mathrm{Pr}\left(|q_{i}-\mu_{Q}|>\frac{\sqrt{|{\cal A}|}}{2}\right)+\mathrm{Pr}\left(|p_{i}-\mu_{P}|>\frac{\sqrt{|{\cal A}|}}{2}\right). \end{eqnarray*} \noindent$$

Our upper bounds for both of these probabilities are established using Chebyshev's inequality. Chebyshev's inequality states that, for a probability distribution P with mean μ and standard deviation σ that P(|x − μ| ⩾ kσ) ⩽ (σ/k)². In our case, the two standard deviations are $\sqrt {[V_{\rho _i}]_{11}}$ and $\sqrt {[V_{\rho _i}]_{22}}$ which gives us

$$\begin{equation} \mathrm{Pr}_{\mathrm{quant}}\left(\boldsymbol{u}_{i}\notin{\cal A}_{\rho_i}\right)\leqslant \frac{4([V_{\rho_i}]_{11}+[V_{\rho_i}]_{22})}{|{\cal A}|}. \end{equation} \noindent \tag{ A.4 }$$

It immediately follows from the independence of the n distributions that compose W_ρ that

$$\begin{eqnarray*} &&\mathrm{Pr}_{\mathrm{quant}}(\boldsymbol{u}\notin{\cal A}_{\rho})\leqslant4n\max_{i=1,\ldots,n}\left(\frac{\left[V_{\rho_i}\right]_{11}+\left[V_{\rho_i}\right]_{22}}{|{\cal A}|}\right). \end{eqnarray*} \noindent$$

By identical reasoning,

$$\begin{eqnarray*} &&\mathrm{Pr}_{\mathrm{quant}}(\boldsymbol{k}\notin{\cal A}_{{\cal M}(V_{{\cal M}})})\leqslant4n\max_{j=1,\ldots,n}\left(\frac{[V_{{\cal M},j}]_{11}+[V_{{\cal M},j}]_{22}}{|{\cal A}|}\right). \end{eqnarray*} \noindent$$

Thus we have

$$\begin{eqnarray*} &&\fl \mathrm{Pr}_{\mathrm{quant}}(\boldsymbol{k}\notin{\cal A}_{\rho})+\mathrm{Pr}_{\mathrm{quant}}(\boldsymbol{u}\notin{\cal A}_{{\cal M}(V_{{\cal M}})})\leqslant4n\max_{i,j}\left(\frac{\left[V_{\rho_i}\right]_{11}+\left[V_{\rho_i}\right]_{22}+[V_{{\cal M},j}]_{11}+[V_{{\cal M},j}]_{22}}{|{\cal A}|}\right) \end{eqnarray*} \noindent$$

so choosing $|{\cal A}|\geqslant 8n\max _{i,j}\left (\frac {\left [V_{\rho _i}\right ]_{11}+\left [V_{\rho _i}\right ]_{22}+[V_{{\cal M},j}]_{11}+[V_{{\cal M},j}]_{22}}{\epsilon /2}\right )$ guarantees

$$\begin{equation} \mathrm{Pr}_{\mathrm{quant}}(\boldsymbol{u}\notin{\cal A}_{\rho})+\mathrm{Pr}_{\mathrm{quant}}(\boldsymbol{k}\notin{\cal A}_{{\cal M}(V_{{\cal M}})})\leqslant\frac{\epsilon}{4}. \end{equation} \noindent \tag{ A.5 }$$

We must now bound the discretization error on the truncated distributions. It is natural to break this error up into three pieces corresponding to the first three steps of the algorithm. First, there is the error introduced by discretizing the initial sampling distribution over phase space. Next there is the numerical error introduced in implementing the affine transformation $\tilde {\boldsymbol {u}}=T\boldsymbol {u}+\boldsymbol {x}$. Finally, there is the error introduced by discretizing the outcome distribution. The remainder of the proof is devoted to bounding these three errors and thereby bounding the total error using (A.1).

Step 1 of the simulation algorithm breaks up ${\cal A}_{\rho }$ into hypercubes of volume δ²ⁿ and samples from the set of centers of these hypercubes according to

$$\begin{eqnarray} &&\fl \mathrm{Pr}_{\mathrm{sim},\rho}(\boldsymbol{u}=\boldsymbol{\mu}+(l_{1},\ldots,l_{n},m_{1},\ldots,m_{n})\delta)= \mathrm{Pr}_{\mathrm{sim},\rho_{1}}(l_{1},m_{1})\cdots\mathrm{Pr}_{\mathrm{sim},\rho_{n}}(l_{n},m_{n})\nonumber\\ &&={\cal W}_{\rho_{1},\eta}(l_{1},m_{1})\cdots{\cal W}_{\rho_{n},\eta}(l_{n},m_{n})\cdot\delta^{2n}\nonumber\\ &&=W_{\rho_{1}}(\boldsymbol{u}_{1})\cdots W_{\rho_{n}}(\boldsymbol{u}_{n})\cdot\delta^{2n}+\boldsymbol{\epsilon}_1, \end{eqnarray} \noindent \tag{ A.6 }$$

where ₁ is the numerical error introduced by using η > 0. The probability weighting assigned to a hypercube center u is only approximately equivalent to the probability mass contained in the hypercube; we must bound the error introduced by this non-equivalence. Define ${\cal N}_{\delta }(\boldsymbol {u})\subset \mathbb {R}^{2n}$ to be the hypercube with side length δ centered at the point u. The for a fixed u in the domain of $\mathrm {Pr}_{\mathrm {sim},\rho}(\boldsymbol {u})$ we have that every point $\boldsymbol {w}\in {\cal N}_\delta (\boldsymbol {u})$ is also in the region of truncation ${\cal A}$ hence W^trunc_ρj(w) = W_ρj(w) for all such w. We then use this simplifying observation, the mean value theorem and the triangle inequality to find

$$\begin{eqnarray} &&\fl \left|\mathrm{Pr}_{\mathrm{sim},\rho}(\boldsymbol{u})-\int_{{\cal N}_{\delta}(\boldsymbol{u})}W_{\rho}(\boldsymbol{v})\,\mathrm{d}\boldsymbol{v}\right| \leqslant \left(\max_{\boldsymbol{v}\in{\cal N}_{\delta}(\boldsymbol{u})}|\nabla W_{\rho}(\boldsymbol{v})|\right)\delta^{2n}\max_{\boldsymbol{v}\in{\cal N}_{\delta}(\boldsymbol{u})}|\boldsymbol{v}-\boldsymbol{u}|+|\boldsymbol{\epsilon}_1|\nonumber\\ &&\leqslant \frac{\delta^{2n+1}\beta n\sqrt{n/2}}{|{\cal A}|^{n}}+|\boldsymbol{\epsilon}_1|, \end{eqnarray} \tag{ A.7 }$$

where $\max _{\boldsymbol {v}\in {\cal N}_{\delta }(\boldsymbol {u})}|\boldsymbol {v}-\boldsymbol {u}|\leqslant \sqrt {2n}\frac {\delta }{2}$ from Pythagoras' theorem and $\left (\max _{\boldsymbol {v}\in {\cal N}_{\delta }(\boldsymbol {u})}|\nabla W_{\rho }(\boldsymbol {v})|\right )\leqslant n\beta /|{\cal A}|^n$ by assumption 2 of theorem 1.

In the second step of the simulation algorithm we move the sampled point $\boldsymbol {u}\in \mathbb {R}^{2n}$ to the point $\tilde {\boldsymbol {u}}=T\boldsymbol {u}+\boldsymbol {x}+\boldsymbol {\epsilon }_{2}$, simulating the evolution due to the linear optical transformation U_T,x. The numerical error ₂ depends only on numerical precision which can be made exponentially small using a linear amount of memory. The other source of error in the simulation of the dynamics is due to error in the initial conditions caused by sampling the point $\boldsymbol {u}\in \mathbb {R}^{2n}$ as opposed to the point $\boldsymbol {v}\in \mathbb {R}^{2n}$ that would have been sampled if δ = 0. The triangle inequality then implies that

$$\begin{equation} \fl |\tilde{\boldsymbol{u}}-\tilde{\boldsymbol{v}}|:=\boldsymbol|\left(T\boldsymbol{u}+\boldsymbol{x}+|\boldsymbol{\epsilon}_{2}|\right)-\left(T\boldsymbol{v}+\boldsymbol{x}\right)|\leqslant|\boldsymbol{\epsilon}_2|+\|T\|\max_{\boldsymbol{v}\in{\cal N}_{\delta}(\boldsymbol{u})}|\boldsymbol{v}-\boldsymbol{u}| = |\boldsymbol{\epsilon}_2|+\|T\|\delta\sqrt{\frac{n}{2}}. \end{equation} \noindent \tag{ A.8 }$$

Since ₂ can be made exponentially small using a polynomial number of computational steps, we can choose $|\boldsymbol {\epsilon }_2|\leqslant \|T\|\delta \sqrt {\frac {n}{2}}$ without affecting the efficiency of the algorithm. Therefore, by making such a choice, the total error in the simulated dynamics is at most

$$\begin{equation} |\tilde{\boldsymbol{u}}-\tilde{\boldsymbol{v}}|\leqslant \|T\|\delta\sqrt{{2n}}. \end{equation} \noindent \tag{ A.9 }$$

The final step of the simulation algorithm samples a measurement outcome on each mode. To do so we break up the truncated outcome space ${\cal A}_{{\cal M}(V_{{\cal M}})}$ into hypercubes of volume δ²ⁿ and sample from the set of centers of these hypercubes according to

$$\begin{eqnarray} &&\fl \mathrm{Pr}_{\mathrm{sim},{\cal M}(V_{\cal M})}(\boldsymbol{k}=\tilde{\boldsymbol{\mu}}+(l_{1},\ldots,l_{n},m_{1},\ldots,m_{n})\delta|\tilde{\boldsymbol{u}}) = \mathrm{Pr}_{\mathrm{sim},{\cal M}(V_{{\cal M},1})}(l_{1},m_{1})\cdots\mathrm{Pr}_{\mathrm{sim},{\cal M}(V_{{\cal M},n})}(l_{n},m_{n})\nonumber\\ &&=\frac{1}{\zeta}\exp\kern-1pt\left(-\sum_{i=1}^{n}\delta(l_{i},m_{i})^{\mathrm{T}}V_{{\cal M},i}^{-1}\delta(l_{i},m_{i}))\right)\cdot\delta^{2n}\nonumber\\ &&=\frac{1}{\zeta}\exp\kern-1pt\left(-\sum_{i=1}^{n}(\boldsymbol{k}_{i}-\tilde{\boldsymbol{u}}_{i})^{\mathrm{T}}V_{{\cal M}(V_{{\cal M},i})}^{-1}(\boldsymbol{k}_{i}-\tilde{\boldsymbol{u}}_{i})\right)\cdot\delta^{2n}, \end{eqnarray} \noindent \tag{ A.10 }$$

where ζ is a normalizing constant and k_i are components of k. As in step one we may use the mean value theorem to bound the error introduced by sampling this way rather than according to the true probability mass over each hypercube. Using ${\cal N}_{\delta }(\boldsymbol {k})\subset \mathbb {R}^{2n}$ to be the hypercube with side length δ centered at the point k we find from the mean value theorem and the assumptions of theorem 1 that the error in the probability enclosed in a single hypercube centered at k is at most

$$\begin{eqnarray} &&\fl \left|\mathrm{Pr}_{\mathrm{sim},{\cal M}(V_{{\cal M}})}(\boldsymbol{k}|\tilde{\boldsymbol{u}})-\int_{{\cal N}_{\delta}(\boldsymbol{k})}M_{{\cal M}(V_{{\cal M}})}(\boldsymbol{\kappa}|\tilde{\boldsymbol{u}})\,\mathrm{d}\boldsymbol{\kappa}\right| \leqslant \max_{\boldsymbol{\kappa}\in{\cal A}_{{\cal M}(V_{{\cal M}})}} \left|\nabla_{\boldsymbol{\kappa}} M_{{\cal M}(V_{{\cal M}})}(\boldsymbol{\kappa}|\tilde{\boldsymbol{u}})\delta^{2n}\right|\max_{\boldsymbol{\kappa}\in{\cal N}_{\delta}(k)}|\boldsymbol{\kappa}-\boldsymbol{k}|\nonumber\\ &&\leqslant \frac{\Lambda\delta^{2n+1}n\sqrt{2n}}{|{\cal A}|^{n}}. \end{eqnarray} \tag{ A.11 }$$

We complete the error analysis by bounding the distance between the simulator distribution and truncated, discretized quantum distribution. Our aim is to rewrite this error in terms of quantities we have already bounded. Let ${\cal D}_{\boldsymbol {u}}$ be the domain of $\mathrm {Pr}_{\mathrm {sim},{\cal M}(V_{\cal M})}(\boldsymbol {u})$; i.e. the discrete set of points in ${\cal A}_\rho$ that can be sampled by the simulator. The difference between the probabilities of a fixed outcome k occurring for both the simulation and the truncated quantum distribution

$$\begin{eqnarray} &&\fl \left|\mathrm{Pr}_{\mathrm{sim},{\cal M}(V_{\cal M})}(\boldsymbol{k})-\mathrm{Pr}_{\mathrm{trunc-quant},{\cal M}(V_{\cal M})}(\boldsymbol{k})\right|\nonumber\\ &&\fl\quad =\left|\sum_{\boldsymbol{u}\in{\cal D}_{\boldsymbol{u}}}\mathrm{Pr}_{\mathrm{sim},{\cal M}(V_{\cal M})}(\boldsymbol{k}|\tilde{\boldsymbol{u}}+\boldsymbol{\epsilon}_{2})\mathrm{Pr}_{\mathrm{sim},\rho}(\boldsymbol{u})-\int_{{\cal N}_{\delta}(\boldsymbol{k})}\int_{{\cal A}_{\rho}}M_{{\cal M}(V_{{\cal M},i})}^{\mathrm{trunc}}(\boldsymbol{\kappa}|\tilde{\boldsymbol{v}})W_{\rho}(\boldsymbol{v})\,\mathrm{d}\boldsymbol{v}\mathrm{d}\boldsymbol{\kappa}\right|\nonumber \\ &&\fl\quad =\left|\sum_{\boldsymbol{u}\in{\cal D}_{\boldsymbol{u}}}\mathrm{Pr}_{\mathrm{sim},{\cal M}(V_{\cal M})}(\boldsymbol{k}|\tilde{\boldsymbol{u}}+\boldsymbol{\epsilon}_{2})\mathrm{Pr}_{\mathrm{sim},\rho}(\boldsymbol{u})-\sum_{\boldsymbol{u}\in{\cal D}_{\boldsymbol{u}}}\int_{{\cal N}_{\delta}(\boldsymbol{k})}\int_{{\cal N}_{\delta}(\boldsymbol{u})}M_{{\cal M}(V_{{\cal M},i})}^{\mathrm{trunc}}(\boldsymbol{\kappa}|\tilde{\boldsymbol{v}})W_{\rho} (\boldsymbol{v})\,\mathrm{d}\boldsymbol{v}\mathrm{d}\boldsymbol{\kappa}\right|\nonumber\\ &&\fl \quad \leqslant\sum_{\boldsymbol{u}\in{\cal D}_{\boldsymbol{u}}}\mathrm{Pr}_{\mathrm{sim},{\cal M}(V_{\cal M})}(\boldsymbol{k}|\tilde{\boldsymbol{u}}+\boldsymbol{\epsilon}_{2})\left|\mathrm{Pr}_{\mathrm{sim},\rho}(\boldsymbol{u})-\int_{{\cal N}_{\delta}(\boldsymbol{u})}W_{\rho}(\boldsymbol{v})\,\mathrm{d}\boldsymbol{v}\right|\nonumber \\ &&\fl \qquad+\sum_{\boldsymbol{u}\in{\cal D}_{\boldsymbol{u}}}\int_{{\cal N}_{\delta}(\boldsymbol{u})}\left|\mathrm{Pr}_{\mathrm{sim},{\cal M}(V_{\cal M})}(\boldsymbol{k}|\tilde{\boldsymbol{u}}+\boldsymbol{\epsilon}_{2})-\int_{{\cal N}_{\delta}(\boldsymbol{k})}M_{{\cal M}(V_{{\cal M},i})}^{\mathrm{trunc}}(\boldsymbol{\kappa}|\tilde{\boldsymbol{v}})\mathrm{d}\boldsymbol{\kappa}\right|W_{\rho}(\boldsymbol{v})\,\mathrm{d}\boldsymbol{v}. \end{eqnarray} \tag{ A.12 }$$

The final inequality is found by adding and subtracting $\sum _{\boldsymbol {u}\in {\cal D}_{\boldsymbol {u}}}\mathrm {Pr}_{\mathrm {sim},{\cal M}(V_{\cal M})}(\boldsymbol {k}|\tilde {\boldsymbol {u}}+\boldsymbol {\epsilon }_{2})\int _{{\cal N}_{\delta }(\boldsymbol {u})}W_{\rho }(\boldsymbol {v})\,\mathrm {d}\boldsymbol {v}$ and applying the triangle inequality.

We need one more intermediary bound before arriving at the final result. Let $\boldsymbol {v}\in {\cal N}_{\delta }(\boldsymbol {u})$. It follows from the mean value theorem and $M_{{\cal M}(V_{{\cal M}})}^{\mathrm {trunc}}=M_{{\cal M}(V_{{\cal M}})}$ throughout the domain of integration that

$$\begin{eqnarray} &\left|\int_{{\cal N}_{\delta}(\boldsymbol{k})}M_{{\cal M}(V_{{\cal M}})}^{\mathrm{trunc}}(\boldsymbol{\kappa}|\tilde{\boldsymbol{u}}+\boldsymbol{\epsilon}_{2})\,\mathrm{d}\boldsymbol{\kappa}-\int_{{\cal N}_{\delta}(\boldsymbol{k})}M_{{\cal M}(V_{{\cal M}})}^{\mathrm{trunc}}(\boldsymbol{\kappa}|\tilde{\boldsymbol{v}})\mathrm{d}\boldsymbol{\kappa}\right| \nonumber\\ &\qquad\qquad=\left|\int_{{\cal N}_{\delta}(\boldsymbol{k})}M_{{\cal M}(V_{{\cal M}})}(\boldsymbol{\kappa}|\tilde{\boldsymbol{u}}+\boldsymbol{\epsilon}_{2})\,\mathrm{d}\boldsymbol{\kappa}-\int_{{\cal N}_{\delta}(\boldsymbol{k})}M_{{\cal M}(V_{{\cal M}})}(\boldsymbol{\kappa}|\tilde{\boldsymbol{v}})\mathrm{d}\boldsymbol{\kappa}\right| \nonumber\\ &\qquad\qquad \leqslant \int_{{\cal N}_{\delta}(\boldsymbol{k})}\max_{\boldsymbol{x}\in{\cal A}_{\rho}}\left|\nabla_{\boldsymbol{u}}M_{{\cal M}(V_{{\cal M}})}^{\mathrm{trunc}}(\boldsymbol{\kappa}|\boldsymbol{u})|_{\boldsymbol{u}=\boldsymbol{x}}\right|\cdot\left|T(\boldsymbol{u}-\boldsymbol{v})+\boldsymbol{\epsilon}_{2}\right|\mathrm{d}\boldsymbol{\kappa}\nonumber\\ &\qquad\qquad \leqslant \int_{{\cal N}_{\delta}(\boldsymbol{k})}\frac{n\Lambda}{|{\cal A}|^{n}}\cdot|T(\boldsymbol{u}-\boldsymbol{v})+\boldsymbol{\epsilon}_{2}|\mathrm{d}\boldsymbol{\kappa}\nonumber\\ &\qquad\qquad \leqslant \delta^{2n}\frac{n\Lambda}{|{\cal A}|^{n}}\left(\|T\|\delta\sqrt{\frac{n}{2}}+\|T\|\delta\sqrt{\frac{n}{2}}\right), \end{eqnarray} \tag{ A.13 }$$

where we have used the assumption that $|\boldsymbol {\epsilon }_2|\leqslant \|T\|\delta \sqrt {n/2}$ and identical reasoning to bound $\left |T(\boldsymbol {u}-\boldsymbol {v})\right |$. In conjunction with (A.11) this implies

$$\begin{equation} \fl \left|\mathrm{Pr}_{\mathrm{sim},{\cal M}(V_{\cal M})}(\boldsymbol{k}|\tilde{\boldsymbol{u}}+\boldsymbol{\epsilon}_{2})-\int_{{\cal N}_{\delta}(\boldsymbol{k})}M_{{\cal M}(V_{{\cal M}})}^{\mathrm{trunc}}(\boldsymbol{\kappa}|\tilde{\boldsymbol{v}})\,\mathrm{d}\boldsymbol{\kappa}\right|\leqslant\left(1+\|T\|\right)\frac{n\Lambda\sqrt{2n}\delta^{2n+1}}{|{\cal A}|^{n}}+|\boldsymbol{\epsilon}_1|. \end{equation} \tag{ A.14 }$$

Since there are numerical errors in the computation of the probability distribution for measurement outcomes, it is conceivable that the algorithm assigns more than unit probability to outcomes in the region ${\cal A}$ (although this is unlikely for small values of ). Although we can claim that $\sum _{\boldsymbol {u}\in {\cal D}_{\boldsymbol {u}}}\int _{{\cal N}_{\delta }(\boldsymbol {u})}W_{\rho }(\boldsymbol {v})\,\mathrm {d}\boldsymbol {v}\leqslant 1$ we can only claim that $\sum _{\boldsymbol {u}\in {\cal D}_{\boldsymbol {u}}}\mathrm {Pr}_{\mathrm {sim},{\cal M}(V_{\cal M})}(\boldsymbol {k}|\tilde {\boldsymbol {u}}+\boldsymbol {\epsilon }_{2})\leqslant 1+\epsilon \leqslant 2$ under the assumption that  < 1. Using these facts, we substitute (A.7) and (A.14) into (A.12) to find

$$\begin{eqnarray*} &&\fl \left|\mathrm{Pr}_{\mathrm{sim},{\cal M}(V_{\cal M})}(\boldsymbol{k})-\mathrm{Pr}_{\mathrm{trunc-quant},{\cal M}(V_{\cal M})}(\boldsymbol{k})\right|\leqslant\left[\left(1+\|T\|\right)\Lambda+\beta\right]\frac{n\sqrt{2n}\delta^{2n+1}}{|{\cal A}|^{n}}+|\boldsymbol{\epsilon}_1|. \end{eqnarray*}$$

Since we have assumed that the outcomes are discretized into hypercubes of volume δ²ⁿ, the discretization produced $\frac {|{\cal A}|^{n}}{\delta ^{2n}}$ hypercubes. This implies that the 1-norm distance between the the two distributions is:

$$\begin{eqnarray*} &&\fl \sum_{\boldsymbol{k}}\left|\mathrm{Pr}_{\mathrm{sim},{\cal M}(V_{\cal M})}(\boldsymbol{k})-\mathrm{Pr}_{\mathrm{trunc-quant},{\cal M}(V_{\cal M})}(\boldsymbol{k})\right|\leqslant\left[\left(1+\|T\|\right)\Lambda+\beta\right]n\sqrt{2n}\delta+|\boldsymbol{\epsilon}_1||{\cal A}|^n/\delta^{2n}. \end{eqnarray*} \noindent$$

Since there are n modes and numerical error η, it is straightforward to see that |₁| ⩽ nηδ²ⁿ; therefore choosing

$$\begin{eqnarray*} &&\delta\leqslant\frac{\epsilon}{8\left[\left(1+\|T\|\right)\Lambda+\beta\right]n\sqrt{2n}} \end{eqnarray*} \noindent$$

and

$$\begin{eqnarray*} &&\eta\leqslant \frac{\epsilon}{8n|{\cal A}|^{n}} \end{eqnarray*} \noindent$$

implies

$$\begin{equation} |\mathrm{Pr}_{\mathrm{sim},{\cal M}(V_{\cal M})}(\boldsymbol{k})-\mathrm{Pr}_{\mathrm{trunc-quant},{\cal M}(V_{\cal M})}(\boldsymbol{k})|_{1}\leqslant\frac{\epsilon}{4}. \end{equation} \noindent \tag{ A.15 }$$

This sets the discretization error to be /4. We found previously that the truncation error is at most /4 under assumption 5 of theorem 1. We then combine (A.1), (A.5) and (A.15) and find:

$$\begin{eqnarray*} &&|\mathrm{Pr}_{\mathrm{sim},{\cal M}(V_{\cal M})}(\boldsymbol{k})-\mathrm{Pr}_{\mathrm{quant},{\cal M}(V_{\cal M})}(\boldsymbol{k})|_{1}\leqslant\epsilon/2 \end{eqnarray*}$$

if Γ = δ. This also implies that the result holds for Γ ⩾ δ given that Γ is an integer multiple of δ (which ensures that the hypercube of size δ²ⁿ that the measurement is assigned to is inside the hypercube of size Γ²ⁿ that it should be assigned to given Γ-discretization); therefore, it generally holds if

$$\begin{equation} \delta\leqslant\min\left\{\frac{\epsilon}{8\left[\left(1+\|T\|\right)\Lambda+\beta\right]n\sqrt{2n}},\Gamma\right\}. \end{equation} \tag{ A.16 }$$

A final issue remains: we did not consider errors that are introduced in the sampling steps in the algorithm that arise because the simulated probability distributions are not normalized. We will show that, from the assumption that  < 1, it suffices to divide by 2 in all previous calculations. To see this, let y be a discrete probability distribution and let x be a vector such that |y − x|₁ ⩽ ' for some 1 > ' ⩾ 0. It then follows from the triangle inequality that 1 + ' ⩾ |x|₁ ⩾ 1 − '. Thus

$$\begin{eqnarray} \big|\boldsymbol{y}-\boldsymbol{x}/|\boldsymbol{x}|_1\big|_1&=\big|\boldsymbol{y}|\boldsymbol{x}|_1-\boldsymbol{x}\big|_1/|\boldsymbol{x}|_1\nonumber\\ &=\big|\boldsymbol{y}|\boldsymbol{x}|_1-\boldsymbol{x}|\boldsymbol{x}|_1+\boldsymbol{x}|\boldsymbol{x}|_1-\boldsymbol{x}\big|_1/|\boldsymbol{x}|_1\nonumber\\ &\leqslant|\boldsymbol{y} -\boldsymbol{x}|_1+ ||\boldsymbol{x}|_1-1|\leqslant 2\epsilon', \end{eqnarray} \tag{ A.17 }$$

under the assumption that ' < 1. Therefore, by combining these results with (A.16) we see that the difference between the (now normalized) simulated probabilities and the quantum predictions is at most if

$$\begin{equation} \delta\leqslant\min\left\{\frac{\epsilon}{16\left[\left(1+\|T\|\right)\Lambda+\beta\right]n\sqrt{2n}},\Gamma\right\} \end{equation} \tag{ A.18 }$$

and

$$\begin{equation} |{\cal A}|\geqslant 16n\max_{i,j}\left(\left[V_{\rho_i}\right]_{11}+\left[V_{\rho_i}\right]_{22}+\left[V_{{\cal M},j}\right]_{11}+\left[V_{{\cal M},j}\right]_{22}\right)/\epsilon,\end{equation} \tag{ A.19 }$$

as claimed by theorem 1.

We complete the proof by showing that algorithm 2 is computationally efficient with this choice of $|{\cal A}|$ and δ. Again we will analyze the simulation algorithms step by step. In the first step we sample a point u_i on the phase space of each register from the distribution $\mathrm {Pr}_{\mathrm {{sim}},\rho }(\boldsymbol {u}_{i})$. This distribution has support on $\frac {|{\cal A}|}{\delta ^{2}}$ squares, and thus if we take $|{\cal A}|$ and δ to be proportional to their respective lower and upper bounds then the number of times that ${\cal W}_{\rho _{j},\eta }$ must be evaluated scales as $\Theta \left (n^{4}\left (\|V_\rho \|+\|V_{{\cal M}}\|\right )\left [\Lambda ^2\|T\|^2+\beta ^2\right ]/\epsilon ^3\right )$ (here Θ(·) is Bachmann–Landau notation meaning asymptotically bounded above and below by a constant multiplied by (·)). If we ascribe unit computational cost to every such access, then the computational complexity of this step is proportional to the number of times that ${\cal W}_{\rho _{j},\eta }$ is queried. This task must be repeated n times, and hence the total computational complexity of this step is

$$\begin{eqnarray*} &&\Theta\left(n^{5}\left(\|V_\rho\|+\|V_{{\cal M}}\|\right)\left[\Lambda^2\|T\|^2+\beta^2\right]/\epsilon^3\right). \end{eqnarray*}$$

In step two we apply the affine transformation u → Tu + x, whose cost is dominated by the cost of performing a matrix multiplication using $O(\log (\|T\|\delta \sqrt {n}))$ bits of precision. Since the matrix multiplication requires a number of arithmetic operations that scales quadratically with the matrix dimension and addition and multiplication scale at most quadratically with the number of bits of precision, the total cost of this step is $O(n^{2}\log ^2(\|T\|\delta \sqrt {n}))$, which is subdominant to the cost of the previous step and therefore does not affect the scaling.

In step three, we are confronted with the task of measuring the resultant trajectory using a separable Gaussian measurement. There is a strong duality between drawing a sample trajectory from the initial separable Wigner function and drawing a measurement outcome for the separable Gaussian measurement of the final trajectory. In particular, we discretize the space surrounding the outcome space of each of the n separable measurements into $|{\cal A}|/\delta ^2$ points. The approximate calculation of the measurement probability requires that we perform a number of operations that are proportional to the number of points. Therefore, identically to step one, the computational cost is

$$\begin{eqnarray*} &&\Theta\left(n^{5}\left(\|V_\rho\|+\|V_{{\cal M}}\|\right)\left[\Lambda^2\|T\|^2+\beta^2\right]/\epsilon^3\right), \end{eqnarray*}$$

which verifies the computational complexity claimed by theorem 1 and shows that under the assumptions of the theorem all three steps are computationally efficient, and hence the algorithm is efficient as well.   □