Sufficient Conditions for Efficient Classical Simulation of Quantum Optics

We provide general sufficient conditions for the efficient classical simulation of quantum-optics experiments that involve inputting states to a quantum process and making measurements at the output. The first condition is based on the negativity of phase-space quasiprobability distributions (PQDs) of the output state of the process and the output measurements; the second one is based on the negativity of PQDs of the input states, the output measurements, and the transition function associated with the process. We show that these conditions provide useful practical tools for investigating the effects of imperfections in implementations of boson sampling. In particular, we apply our formalism to boson-sampling experiments that use single-photon or spontaneous-parametric-down-conversion sources and on-off photodetectors. Considering simple models for loss and noise, we show that above some threshold for the probability of random counts in the photodetectors, these boson-sampling experiments are classically simulatable. We identify mode mismatching as the major source of error contributing to random counts and suggest that this is the chief challenge for implementations of boson sampling of interesting size.


I. INTRODUCTION
It is generally believed that quantum computers can perform certain tasks faster than their classical counterparts. Identifying the resource that enables this speedup is of particular interest in quantum information science. Attempts to identify the elusive quantum feature are generally back-door attacks, studying not what is essential for speedup, but rather what is lacking in quantum circuits that can be efficiently simulated classically. A promising candidate resource comes from the result that, in general, there is no quantum speedup for circuits whose initial states and operations have nonnegative Wigner functions [1][2][3]. This suggests that negativity of the Wigner function [4], which can also be viewed as quantum interference [5], is a necessary resource for quantum speedup.
Of particular interest are quantum-optical models of computation that seem achievable in the near future. There has been considerable interest in boson sampling [6] as an intermediate model of quantum computation that, despite its simple physical implementation, is believed to be classically hard. In this model, N single photons are injected into N ports of an M -port (M N ), lossless, passive linear-optical network (LON). 1 The remaining M − N input ports receive vacuum states. Using on-off photodetectors at the output of the network, one samples from the output photon-counting probability distribution. This output distribution, in general, is propor-tional to the squared modulus of the permanent of a complex matrix. Computing permanents is a difficult problem, known to be #P hard in complexity theory [7,8]. In their original proposal, Aaronson and Arkhipov (AA) provided strong evidence that sampling from the output distribution cannot be simulated efficiently classically [6], and this has come to be known as the boson-sampling problem.
Subsequent studies showed that there are other product input states for which sampling from the output probability distribution is classically hard [9][10][11][12]. This gives rise to questions about the classes of input states and measurements for which sampling the output distribution is classically intractable. Given the well-developed theory of phase-space quasiprobability distributions (PQDs) for bosonic states and measurements, an inevitable question asks whether negativity of such PQDs is required for classical intractability of the sampling problem. In addition, a question of both fundamental and practical importance concerns the effect of imperfections on the classical intractability of sampling problems. There have been various investigations of the effect on boson sampling of imperfections in the LON [13][14][15][16] and of mode mismatching [17].
The present paper makes two contributions to this discussion. The first, developed in Sec. II, is to formulate two sufficient conditions for efficient classical simulation of generic quantum-optics experiments: M bosonic modes prepared in an arbitrary bosonic state undergo an M -mode (tracepreserving) quantum process; one generates samples by making a measurement on the M output modes (see Fig. 1). The first condition is based on expressing the probability distribution of the measurement outcomes in terms of a PQD for the output state of the process and PQDs of the elements of the Positive-Operator-Valued Measure (POVM) that describe the output measurement. If the PQD of the output state can be efficiently computed and if for some operator orderings all the PQDs are nonnegative, then efficient classical simulation of sampling is possible. Our second condition generalizes a previous no-go theorem [2], which was given in terms of the Wigner function, and it is particularly useful when one cannot efficiently compute the PQD of the output state. For this condition, we derive a relation that decomposes the output probability distribution into three parts: a PQD for the phase-space complex amplitudes of the input state; a transition function associated with the quantum process, which is a conditional PQD for the output complex amplitudes of the process given the input complex amplitudes; and PQDs of the measurement POVM elements. If specific operator orderings exist such that all these PQDs-input, transition, and output-are efficiently describable and nonnegative, sampling from the output probability distribution can be efficiently simulated classically. These conditions show that negativity is a necessary resource for a generic quantum-optics experiment not to be efficiently simulatable; the result includes boson sampling as the special case where the quantum process is a LON. We emphasize that efficient classical simulation might still be possible using other methods even if our conditions are not satisfied.
Our second contribution, developed in Sec. III, is to apply the results of the first one to investigate the effects of imperfections on implementations of boson sampling that use single-photon states [6] or two-mode squeezed-vacuum states [9] as inputs and photodetection at the output. The imperfections we consider are loss and mode mismatching at the input to and within the LON and subunity efficiency and random counts in the photodetectors. Considering simple models for these errors, we find necessary and sufficient conditions for the relevant PQDs to be nonnegative, and thus for such boson-sampling implementations to be efficiently simulated classically using these methods. These conditions say that an experiment can be classically simulated when the probability of random counts per photodetector exceeds some threshold in the experiment.
The various sources of error we consider are not completely independent. The random counts at the photodetectors include both intrinsic dark counts and counts due to mode-mismatched photons, i.e., nonoverlapping parts of photon wave packets. These mode-mismatched photons are lost to the interference that gives rise to the desired output photocount distribution. They are part of the losses within the apparatus, but they can make their way through the LON to the photodetectors and be counted within the spatiotemporal windows of the detectors. They thus contribute essentially random counts within the photodetectors.
As we discuss in Secs. III B and III C, our conditions for classical simulatability are not a challenge for situations with practical losses and high-quality photodetectors, if the only source of random counts is the intrinsic dark counts in the detectors. The chief challenge for boson sampling, we believe, comes when a substantial number of mode-mismatched photons reach the detectors and are counted as random counts. A good, but not exact rule of thumb is that our methods can clas- The input state ρin is processed through an M -mode quantum process described by quantum operation E, producing the output state ρout = E(ρin); an output probability distribution p(n) = Tr ρout Πn is sampled by measuring a POVM {Πn}.
sically simulate a boson-sampling experiment when the number of mode-mismatched photons reaching the photodetectors exceeds the number of mode-matched photons. The analysis in Secs. III B and III C suggests that this is an important challenge for implementations of boson sampling of interesting size, i.e. when the size of the system is sufficiently large to represent a challenge for a classical computer to sample.
The paper concludes with a discussion of outstanding issues in Sec. IV.

II. SIMULATION OF GENERIC SAMPLING PROBLEMS
A. Generic quantum-optics sampling problem We consider the generic quantum-optics sampling problem depicted in Fig. 1: M input bosonic modes, with overall density operator ρ in , traverse a quantum process described by a (trace-preserving) quantum operation E, leading to the output state ρ out = E(ρ in ); at the output, one measures the M modes, thus sampling from a probability distribution where POVM elements Π n , with n denoting the joint outcome, characterize the measurement. The POVM satisfies a completeness relation, n Π n = I, with I denoting the identity operator in the M -mode Fock space. The quantum operation E is a linear, trace-preserving, completely positive map from the set of all density operators associated with a quantum system to itself. In general, such a quantum operation describes the system dynamics associated with a joint unitary transformation on the system and an environment and arises formally from tracing out the environment [18]. The question is whether sampling from the output probability distribution (2.1) can be efficiently simulated on a classical computer. If such classical simulation is possible, then using Stockmeyer's approximate counting algorithm [19], one can approximate the output probability to within a multiplicative error in BPP NP , which is contained in the third level of the polynomial hierarchy;p(n) approximates p(n) to within a multiplicative factor g if p(n)/g ≤p(n) ≤ p(n)g.
Ideal boson sampling is a special case of this general problem, for which the input state is a multimode Fock state with N single photons, the quantum process is a lossless LON described by an M × M unitary matrix U with M N , and photon-counting measurements are made on each output mode. The output probabilities, in general, are proportional to the squared modulus of permanents of complex matrices, which are, in the likely event of single-photon detections, submatrices of U [20]. Computing permanents is a difficult problem, known to be #P hard [7,8], and in a class that contains the entire polynomial hierarchy of complexity classes [21].
The key observation by Aaronson and Arkhipov (AA) was that multiplicative approximation of the squared modulus of permanents of real matrices is also a #P-hard problem, and it is likely this is the case for general complex matrices [6]. If boson sampling were classically simulatable, one could use Stockmeyer's approximate counting algorithm to approximate the output probability to within a multiplicative error, and this would lead to the collapse of the polynomial hierarchy to the third level, which is thought to be unlikely [6]. Given two plausible conjectures, AA further showed that it is likely that classical simulation of sampling from probability distributions that closely approximate the ideal output probability, known as approximate boson sampling, is also hard.
The approximate sampling problem is of more practical interest than exact sampling because in an experiment the input quantum state, quantum process, and output measurement are only characterized approximately, so one does not sample from the exact probability distribution (2.1). Moreover, one might not be able to distinguish efficiently sampling from two probability distributions that are close to one another. Hence, in practice, an interesting sampling problem is the one for which approximate sampling is hard. In this paper, we do not consider this form of sampling; instead we focus on simulating exact sampling from output distributions that arise in the presence of errors and imperfections. Even though we do not consider approximate sampling explicitly, our simulation methods for exact sampling do lead to a sufficient condition for approximate sampling to be classically simulatable.
We motivate our methods for classical simulation by considering a simple special case of the generic sampling problem. Suppose the multimode input state ρ in has a nonnegative Glauber-Sudarshan P function [22,23] P (α|ρ in ), i.e., here |α is an M -mode coherent state with phase-space complex amplitudes α = (α 1 , α 2 , . . . , α M ). 2 Such states, as mixtures of coherent states, are often called classical states. Suppose further that the quantum process transforms multimode coherent states to classical states; such processes are known as classical processes [24]. Then the output state ρ out is classical as well and has a nonnegative P function. Using the linearity of quantum processes over density operators, the output state can be expressed as where P E (β|α) is the P function of the state E |α α| . Hence, the output probability distribution (2.1) is given by where the Husimi Q functions [25] of the POVM elements, Q Π (n|β) = β| Π n |β /π M , are always nonnegative and satisfy π M n Q Π (n|β) = 1. As all the PQDs in the expression (2.4) are nonnegative, sampling from the output photon-counting probability distribution can be efficiently simulated on a classical computer, provided the PQDs can be efficiently generated: α is chosen from the input probability distribution P (α|ρ in ); given α, β is chosen according to the probability distribution P E (β|α); finally, given β, n is chosen according to the measurement distribution π M Q Π (n|β). Applying this procedure to input thermal states for LONs and using Stockmeyer's approximate counting algorithm, it was shown that permanents of Hermitian positive-semidefinite matrices can be approximated to within a multiplicative constant in BPP NP [26].
In the next subsection, we generalize this approach to nonclassical states and processes. The key idea, taken over from classical states and processes, is to use phase-space complex amplitudes and associated PQDs to describe the input state, the transition through the quantum process, and the output measurements. The generalization is to use operator orderings more general than the normal ordering of the P function and the antinormal ordering of the Q function, thus expanding the possibilities for finding nonnegative PQDs.

B. Sufficient conditions for efficient classical simulation
To formulate our condition for efficient classical simulation of the generic sampling problem, we use the s-ordered phasespace quasiprobability distributions [(s)-PQDs] of a Hermitian operator ρ, which are defined by [27,28] is the corresponding s-ordered characteristic function, with being the M -mode displacement operator, a = (a 1 , . . . , a M ) the row vector of modal annihilation operators, and a † = (a † 1 , . . . , a † M ) T the column vector of creation operators. The diagonal matrix s = diag(s 1 , s 2 , . . . , s M ) has the various ordering parameters on the diagonal. Equation (2.5) is a Fourier transform, which can be inverted using Because ρ is Hermitian, the characteristic function satisfies Φ (2.10) These are usually called PQDs when ρ is a density operator, but we generalize the terminology to any Hermitian operator so we can apply it to POVM elements. The outcome probabilities (2.1) can be expressed in terms of the PQDs of the output state and the POVM as [27] These measurement (−s)-PQDs are normalized according to for any values of β and s, as one sees by applying Tr[D(ξ)] = π M δ 2M (ξ) to Eq. (2.12).
First condition. We now present a sufficient condition for efficient classical simulation of the sampling problem. If there exist values of s such that the PQDs in Eq. (2.11) are nonnegative, they can be used to simulate sampling from p(n) in two steps: (i) The vector of complex amplitudes β is chosen from the probability distribution W (s) (β|ρ out ).
(ii) For the given β, the outcome n is chosen from the prob- This condition is particularly useful if the (s)-PQD of the output state can be efficiently computed, as it can be, for example, for Gaussian input states and Gaussian processes [29]. For cases where efficient computation of the output (s)-PQD is not possible, we now derive a general expression that relates the (s)-PQD of the output state ρ out to the (t)-PQD of the input state ρ in . This derivation introduces the transfer function, which transfers complex amplitudes from input to output of the quantum process and which depends on both the input and output operator orderings.
An M -mode input state can be expanded in terms of displacement operators, (2.14) In this expression, we can replace D † (ξ) = D(−ξ) with D(ξ) if we wish, because ρ in is Hermitian. Linearity of quantum processes implies that from which we obtain the (s)-ordered characteristic function of the output state, Using the Fourier transform (2.5) and its inverse (2.9), we can obtain the relation between the input and output PQDs, where the transition function associated with the quantum process is defined by The quantity Tr E D † (ξ) D(ζ) gives the "matrix elements" of the quantum process E in the displacement-operator basis. We can use the antinormally ordered form of the displacement operator, combined with the coherent-state resolution of the identity, I = d 2M γ |γ γ| /π M , to obtain This allows us to convert E D † (ξ) into the action of the quantum process on coherent states:  .17), we can assemble the ingredients for a classical simulation of sampling from the output distribution of the quantum circuit shown in Fig. 1, Second condition. We can carry out a classical simulation, using the following procedure, if there exist values of t and s such that the PQDs of the input, the transition function, and the measurement are all nonnegative and efficiently describable: (i) The vector of complex amplitudes α is chosen from the input probability distribution W (t) (α|ρ in ).
(ii) For the given α, the vector β is chosen from the transition probability distribution T (s,t) E (β|α).
(iii) For the given β, the outcome n is chosen from the output That the three probability distributions are efficiently describable must be judged on a case-by-case basis. For the input, this is generally achieved by assuming that the input state ρ in is a product state of the M modes or, perhaps, a product of blocks of modes of fixed size; likewise, for the output, the measurements are generally product measurements of the M modes or products of measurements on blocks of fixed size. The complexity of the transition function depends on the quantum process; for the LONs used in boson sampling, the transition function is Gaussian and can be generated trivially from the matrix that describes the LON, as we show in Sec. III.
This second condition includes the previous results as special cases. For classical states and classical processes, we can choose s = t = I M , and Eq. (2.22) reduces to Eq. (2.4). In addition, for s = t = 0, we have that when the transition function is nonnegative and the input quantum state and output POVM elements have nonnegative Wigner functions, sampling from the output distribution can be simulated classically [2].
A procedure for determining if there are input and output orderings that give nonnegative quasidistributions is the following. The (t)-PQD W (t) (α|ρ in ) of the input state is nonnegative for ordering parameters t ≤t, wheret ≥ −I M , and the (−s)-PQD W (−s) Π (n|β) of the POVM elements is nonnegative for s ≥s, i.e., s k ≥s k , for all k, wheres ≤ I M . The necessary and sufficient condition for there also to be a nonnegative transition function is that T (s,t) (β|α) be nonnegative and no more singular than a δ function. This is a necessary and sufficient condition for our second condition to yield a classical simulation.
A crucial point, in general and for what follows, is that simulations using our first condition provide tighter bounds for classical simulation than the second condition because the first condition, unlike the second, does not require anything about the input PQD, in particular, that it be nonnegative. To make the difference clear, suppose the input state is a nonclassical Gaussian state, and the process is also a nonclassical process, but one that cancels the nonclassicality of the input state in such a way as to make the output state classical; an example is provided by input squeezed states and antisqueezing operations. In this situation, for any measurement at the output, the experiment is classically simulatable according to the first condition. For the second condition, however, the input PQD and the transition function are only nonnegative for a limited range of ordering parameters t ≤t < I M , and only for certain measurements does the second condition hold.

III. EFFICIENT CLASSICAL SIMULATIONS OF IMPLEMENTATIONS OF BOSON SAMPLING
A. General considerations for passive linear optical networks We now consider the case where the quantum process is a lossy M -mode LON. In this case the quantum process takes coherent states to coherent states according to [30] E LON |γ γ| = |γL γL| , loss (environment) modes that are initially in vacuum and that carry away photons lost within the LON; the larger LON that includes the loss modes is described by a unitary operatorŨ, which transforms annihilation operators according tõ where a 0 is the row vector of annihilation operators for the M loss modes andŨ is the unitary matrix that describes the complex-amplitude transformation within the larger LON. The larger LON takes overall coherent states to overall coherent states according tõ The quantum operation (3.1) follows from tracing out the loss modes: What the model teaches is that L is a submatrix of the larger unitary matrixŨ and thus satisfies L † L = I M −P † P ≤ I M .
In an experiment, the transfer matrix L of any LON can be efficiently characterized by inputting coherent states [30].
Thus the transition function (2.18) becomes T (s,t) The transition function is well behaved and nonnegative, and has the final (normalized) Gaussian form, if and only if i.e., Σ is positive (semidefinite). Note that if we choose the same ordering at input and output, i.e., s = t = sI M , then To apply our second method for generating an efficient classical simulation of sampling, we should apply the procedure outlined at the end of Sec. II B. Suppose the input state has nonnegative (t)-PQD W (t) (α|ρ in ) for t ≤t, and the output measurement has nonnegative (−s)-PQD W (−s) Π (n|β) for s ≥s. Then the necessary and sufficient condition for our second method to yield an efficient classical simulation of sampling from the output probability distribution is that (3.12) Two special cases deserve attention. For a lossless LON, the transfer matrix L = U is unitary, and the condition (3.12) becomess ≤ U †t U . In the case of identical measurements on all the output modes, the POVM elements become a product Π n = M k=1 Π n k , where {Π n k } is the POVM for the measurement on output mode k, and the (−s)-PQD of the measurements is also a product, W (n k |β k ). In this situation, the optimal output ordering parameters are the same for all M modes, i.e.,s = sI M . Thus, for a lossless LON with identical product measurements, the condition (3.12) simplifies tosI M ≤ U †t U , which is equivalent tosI M ≤t.
In the next two subsections we apply our conditions for classical simulations to two schemes for boson sampling in the presence of errors. Before turning to that, however, we digress briefly to note that since a LON is a classical process, we can provide a classical simulation for all classical input states, since we can choose t = s = I M , i.e., the P function for the input state and the (always nonnegative) Q function for the measurements; this leads to the δ transition function of Eq. (3.11). This is the motivating case considered at the end of Sec. II A. A particular example is provided by inputting coherent states to an LON and performing any measurements at the output.
The flip side of classical input states is classical measurements, such as heterodyne measurements, for which the P functions of the POVM elements are nonnegative, allowing us to choose s = −I M ; in this situation, we can choose t = −I M , i.e., the Q function for the input state, and have a nonnegative transition function according to Eq. (3.10). Hence, in the case of classical measurements, efficient classical simulation is possible for any input state.
In the symmetric case, where both the input state and the POVM elements have nonnegative Wigner functions, we can choose t = s = 0, and given Eq. (3.10), the transition function is always nonnegative. An example is Gaussian input states to a LON and Gaussian measurements at the output [31].

B. Boson sampling with single-photon sources
We turn now to investigating the effect of errors in a practical implementation of boson sampling that uses singlephoton sources and on-off photodetectors. Recall that in this model, which is the one proposed originally by Aaronson and Arkhipov [6], N single photons are injected into the first N ports of an M -port LON, with M N , and the remaining N − M ports receive the vacuum state. To avoid having more than one count at an output detector, one generally requires that the number of photons counted at the detectors be √ M [6,9]. We consider the following sources of error: impurity of the input photons and mode mismatching of these photons into the LON, losses and mode mismatching within the LON, and inefficiency and random counts in the detectors.
It is a considerable practical challenge to generate a singlephoton state. We assume that the output of the single-photon sources is a statistical mixture of vacuum and a single pho-ton, (1 − µ) |0 0| + µ |1 1|, µ ∈ [0, 1]. Note that this state is the output of a beamsplitter with transmissivity √ µ when the beamsplitter is illuminated by a pure single photon. In addition to the impurity of the input, the input photons are generally not mode-matched to the temporal, frequency, and polarization modes that interfere ideally through the LON. The nonoverlapping parts of the input photons are lost to the ideal interference that leads to the probability distribution one wants to sample at the output, so we treat them as a loss and model that loss by virtual beamsplitters with transmissivity √ η B . Taking into account both the impurity and the mode mismatching, we have that the state input into the first N ports is whereη = µη B . We return to a discussion of mode mismatching, at the input to and throughout the LON, at the end of this subsection. By using Eqs. (2.5) and (2.6) for a single mode, the (t)-PQD of the mixed input state (3.13) is given by [32] which is nonnegative for t ≤t = 1 − 2η = 1 − 2µη B . As the vacuum state (η = 0) is a classical state whose (t)-PQD is nonnegative for t ≤ 1, the overall (t)-PQD of the input is nonnegative for 15) where J N is the diagonal matrix with 1s in the first N diagonal positions and 0s otherwise. Losses within the LON are taken into account by the transfer matrix L. For a particular implementation of boson sampling, one should use the measured transfer matrix to analyze the system [30]. A good part of these losses is mode mismatching within the network, about which we say more below. For our analysis, we adopt a simple model of losses that allows us to investigate how the effect of losses scales with the size of the network. In particular, we assume that all paths through the LON suffer the same amount of loss and thus describe the network by a transfer matrix where U is the unitary transfer matrix for a lossless LON. We make this more specific in the following way. The network consists of -port optical elements, each with a uniform transmissivity √ η 0 , and has depth d. Thus each input port speaks to d output ports. We assume that the network is fully connected, so M = d ; hence, each input photon sees a loss . For the on-off photodetectors, which we assume to be identical, we use a model similar to that devised in Ref. [33] for detectors with subunity efficiency and dark counts. We think of the dark-count probability in the model more generally than in the original model, however; it is not just the probability for intrinsic dark counts in the detector, but it also includes any sort of random counts. We discuss below how modemismatched photons at the input and within the LON can propagate through the LON and contribute random counts at the photodetectors. The POVM elements associated with the on-off outcomes-zero denotes the off state, i.e., no detector click, and 1 denotes a click-are (3.17) where η D , satisfying 0 ≤ η D ≤ 1, is the detector efficiency, and 1 − p D is the probability of no random count.
The sum in Π 0 is an unnormalized thermal state, so by using the (s)-PQD of a thermal state, we can find the (−s)-PQD of Π 0 (η D , p D ) to be this is nonnegative provided s ≥ 1 − 2/η D , which is really no restriction at all. The (−s)-PQD of Π 1 (η D , p d ), given by is nonnegative provided that s ≥s = 1 − 2p D /η D . Writing this in terms of the ordering parameters for all the detectors, we have nonnegative output (s)-PQDs if (3.20) Putting Eqs. (3.15) and (3.20), plus our description of the LON, into Eq. (3.12), we find that the condition for an efficient classical simulation is that U ΣU † = (2p D /η D − η L )I M + η Lt ≥ 0; provided there is even one single-photon input, this reduces to the simple condition where η characterizes the overall loss in the experiment. Because of the simplicity of our model for losses within the LON, the condition (3.21) does not apply precisely to LONs with general transfer matrices L, but the dependence of η L on and M does indicate how the condition for simulability scales with the size of the LON. A recent study [34] shows that if a fixed number K of photons are lost, boson sampling remains classically hard, provided K is not too large. This suggests that in the presence of loss, one can inject more single photons into the LON so that on average, an interesting boson-sampling problem is realized. The mean number of photodetector counts is ηN ; if we require that the number of counts does not exceed An obvious question is why our method needs random counts for classical simulability. The answer is that in the absence of random counts, sampling from the (exact) output probability distribution cannot be efficiently simulated classically; this can be shown using Stockmeyer's approximate counting algorithm. Even for large losses, it is still possible that all the input photons get counted by the detectors at the output. As any lossy LON can be thought of as part of a larger, lossless LON, probabilities of these events are proportional to the squared modulus of permanents of complex matrices, which are submatrices of a unitary matrix for the larger LON. If sampling were classically simulatable, using Stockmeyer's approximate counting algorithm one could approximate one of these probabilities to within a multiplicative error in BPP NP [Stockmeyer's algorithm allows for the proportionality factors to be of the order 2 −poly(N ) ]; this would lead to the collapse of polynomial hierarchy to the third level because as observed by Aaronson and Arkhipov, multiplicative approximation of these probabilities is #P hard. Note, however, that this argument does not imply that boson-sampling experiments with losses and very low random counts are still practically interesting. One can expect that above some threshold for losses, a classical algorithm can efficiently generate samples from an approximate probability distribution, and in practice, this cannot be distinguished from the outcomes of the experiment.
The importance of random counts prompts us to return to the question of mode mismatching at the input to and within the LON. Mode mismatching occurs when temporal, frequency, and polarization properties of photon wave packets do not overlap ideally at the input to the LON and at the optical elements used to implement a specific LON. The nonoverlapping parts of the photon wave packets are lost to the ideal interference that leads to the probability distribution one wants to sample at the output. Mode mismatching is thus a loss mechanism and is likely to be the dominant loss mechanism within a large optical network. Without some intervention, however, nonoverlapping parts of the photon wave packets continue through the LON and are counted within the temporal and spatial windows defined by the photodetectors. These photocounts are effectively random and contribute to the random-count probability of the detectors (they might be correlated between different output modes, but it is hard to see how this correlation could be used to our advantage); indeed, they are very likely to be the dominant contribution to the random-count probability, as in high-quality detectors, the intrinsic dark-count rate is very low.
In principle, one can use active filters (mode cleaners) to remove nonoverlapping parts of photon wave packets at the input to and output from and perhaps within the LON and, hence, to turn mode mismatching into a genuine loss where the mode-mismatched parts of the photon wave packets do not contribute counts at the photodetectors. To assess how serious this problem is, suppose that mode mismatching is the dominant loss mechanism. Suppose further that a fraction f B of the photons lost at the input, numbering µ(1 − η B ), continue into the LON and on to the photodetectors and that a fraction f L of the photons lost within the LON, numbering µη B (1−η L ), continue to the detectors. Assuming these modemismatched photons are counted with efficiency η D , they contribute random-count probability With the same assumptions and same values for loss parameters as above, we now also assume that f B = 0.1, on the grounds that the input loss η B = 0.1 already reflects a major attempt to clean up the input wave-packet modes, and f L = 0.9, on the grounds that it would be quite difficult to clean up the output photons without introducing additional losses. With these assumptions, we get p D = 0.046 for M = 10, p D = 0.049 for M = 100, and p D = 0.034 for M = 1 600. Comparing these random counts with the corresponding thresholds from the previous page indicates that mode mismatching is indeed a challenge for boson-sampling experiments of interesting size; recall that this assumes that additional single photons are fed into the LON to compensate for losses in order keep the number of detected photons as close to √ M as possible. The scaling with M is such that if large LONs can be constructed without compromising the loss parameters, the situation gets better as M increases.
It is worth noting that for the range of parameters we have considered, for which N M , the condition for classical simulatability is that the number of mode-mismatched photons counted at the photodetectors, M p D , exceed the number of mode-matched photons, N η. This is a useful rule of thumb for assessing the simulatability of a boson-sampling experiment.

C. Boson sampling with SPDC sources
A major practical challenge for implementing boson sampling is reliable single-photon sources. In most quantumoptics experiments, spontaneous parametric down-conversion (SPDC) is used as a probabilistic source for preparing single photons [35][36][37][38]. If the two-mode squeezed vacuum state generated by a SPDC source has weak squeezing, photon counting on the heralding mode prepares vacuum or a single photon in the signal mode, which can then be used as one of the inputs to the M input ports of a boson-sampling LON. This scheme can be viewed as sampling from the output photoncounting probability distribution of a larger LON with 2M modes; the larger LON consists of the identity process acting on the heralding modes and the original LON acting on the signal modes. This scenario implements randomized boson sampling, in which when N photons are randomly detected in the heralding modes, N single photons are injected into the corresponding ports of the original LON. (With the loss parameters we consider here, in boson sampling with singlephoton sources, the single photons are also randomly injected into a LON, but one does not know to which input ports.) In the absence of any losses or inefficiencies, the average number of photons input to the signal-mode LON is N = M sinh 2 r, where r is the squeezing parameter, assumed to be positive without loss of generality; to achieve N = √ M M , one chooses sinh 2 r = 1/ √ M [9]. We consider the following sources of error: mode mismatching of the signal modes into the smaller, signal-mode LON, described by virtual beamsplitters with transmissivity √ η B ; losses in the signal-mode LON, described by the transfer matrix L, but no losses for the heralding modes, so that the overall transfer matrix is I M ⊕ L; and for all modes, the model for random counts and inefficiency that we introduced previously for on-off detectors. Since the input squeezed vacuum states are Gaussian and the LON is a Gaussian process, we can efficiently find the (t)-PQD of the output state and use our first condition to check whether efficient classical simulation of the sampling problem is possible. To simplify our analysis and to compare directly with our results for singlephoton inputs, we specialize to the simple model for loss in the signal-mode LON in which L = √ η L U . Given this assumption, all the signal modes suffer the same loss in the LON, so we can refer the LON losses to the input, combine them with the mode mismatching of the signal modes, and thus describe both by virtual beamsplitters with transmissivity √ η BL = √ η B η L , which act on each signal mode before it enters the signal-mode LON. The upshot is that the larger LON is fed by M copies of the two-mode state Here ρ hs is the two-mode squeezed vacuum state generated by a SPDC source, U s0 (η BL ) is the unitary operator for a beamsplitter with transmissivity √ η BL that acts on the signal mode of the SPDC and a vacuum input, and the trace is taken over the mode reflected from the beamsplitter. With the LON losses referred to the input, the larger LON is now described by the unitary transfer matrix I M ⊕ U , which corresponds to a δ-function transfer function that does not alter the negativity of the input (t)-PQD. The state ρ hs is a Gaussian state. The Wigner function (t = 0) of any Gaussian state is a Gaussian function, but if the state is nonclassical, there existst ∈ (0, 1] such that for t ≤t, the (t)-PQD is a Gaussian, for t =t, the (t)-PQD has δ-function singularities, and for t >t, the (t)-PQD is more singular than a δ function. In order to findt for ρ hs , we need to use the covariance matrix of the Gaussian (t)-PQD and find the value of t at which the covariance matrix transitions from positive to negative, i.e., the smallest eigenvalue goes to zero.
The covariance matrix of the Wigner function of ρ hs0 in Eq. (3.23) is given by where σ hs = cosh 2rI 2 sinh 2rZ 2 sinh 2rZ 2 cosh 2rI 2 (3.25) is the covariance matrix of the two-mode squeezed vacuum state with squeezing parameter r, with Z 2 = diag(1, −1) be-ing the Pauli z matrix, and is the symplectic transformation of a beamsplitter with transmissivity √ η BL [39]. The 4 × 4 top-left submatrix of σ hs0 is then the covariance matrix of the Wigner function of ρ hs , (3.27) The covariance matrix of the (t)-PQD is given by σ hs − tI 4 ; what we need to know is when, as t increases from zero, the smallest eigenvalue of this 4 × 4 matrix goes to zero. Interchanging rows and columns of σ hs − tI 4 separates it into the direct sum of two 2 × 2 matrices, which have the same eigenvalues (X 2 is the Pauli x matrix). The smaller eigenvalue goes to zero when where we have restored η BL = η B η L . For the (t)-PQD of the overall input, we havet =tI 2M . The analysis of on-off photodetection in Sec. III B shows that nonnegativity of the measurement (s)-PQDs requires s ≥s = 1 − 2p D /η D , and our first method of simulation requires that s = t ≤t, so the condition for efficient classical simulation of the SPDC scheme is thats ≤t, which gives The mean number of photons input to the signal modes is N = M sinh 2 r, meaning that sinh 2 r in the above expressions is a surrogate for N/M . The average number of counts at the photodetectors is η N = η M sinh 2 r, where η = η D η L η B gives the total loss through the system. As in our analysis of single-photon boson sampling, we choose η N = , and sinh 2 r = 0.33; the threshold for classical simulability is p D ≥ 0.060. It is notable that these thresholds are close to twice those we found under comparable conditions for single-photon boson sampling; the single-photon thresholds would be higher if the single-photon sources produced photons with no impurity (µ = 1).
As in single-photon boson sampling, SPDC boson sampling suffers from the problem of mode-mismatched photons becoming random counts in the photodetectors. The same analysis as for single-photon boson sampling yields random-count probability (3.22) with µ = 1. Again assuming f B = 0.1 and f L = 0.9, with all the other parameters the same as above, the random-count probability is p D = 0.091 for M = 10, p D = 0.096 for M = 100, and p D = 0.033 for M = 1 600. This indicates that mode mismatching is a challenge for SPDC boson-sampling experiments of interesting size. Just as for single-photon sources, this conclusion assumes that additional photons are input to compensate for losses, but again the scaling is favorable provided one can keep losses and mode mismatching under control as system size increases.

IV. CONCLUSION
In this paper we established sufficient conditions for efficient classical simulation of general quantum-optical experiments that involve a quantum state that is subjected to an M -mode quantum process and measurement at the output of the process. These conditions support the notion that negativity is a quantum resource by showing that efficient classical simulation of sampling from the output probability distribution is possible when there are (i) nonnegative output-state and output-measurement quasiprobability distributions or (ii) nonnegative input-state and output-measurement quasiprobability distributions and a nonnegative transition function associated with the quantum process.
We applied our conditions for classical simulability to two implementations of the boson-sampling problem. We considered simple models of errors and imperfections to assess the effects of mode mismatching, loss in the LON, and inefficiency and random counts of on-off photodetectors. We found that these errors have a significant impact and obtained random-count thresholds beyond which efficient classical simulation is possible. For any actual implementation of boson sampling, however, one should go beyond the simple examples given here and use our methods to model all the imperfections, noise, and errors, particularly, formulating and analyzing a detailed model of losses and mode mismatching within the particular LON, in order to determine when it is possible to do classical simulations using our methods. In the case of mode mismatching, nonoverlapping parts of photon wave packets that proceed to and are counted at the detectors are likely to be the major contribution to the random-count probability; hence, it is particularly important to assess the need for and effectiveness of active mode cleaning (so-called quantum filters) to mitigate this effect.
We caution that we do not warrant that there is no other method of efficient classical simulation when our conditions are not satisfied. Indeed, we have only considered the problem of sampling from the exact output probability distribution of measurement outcomes. A more general problem is approximate sampling, i.e., sampling from a close approximation to the exact probability distribution, in which case the question is whether sampling from the approximate distribution can be efficiently simulated classically. We have shown that in the presence of losses in boson-sampling experiments and with zero or very low random counts, the exact sampling problem cannot be simulated using our methods. Yet, under the same conditions, one might be able to simulate approximate sampling. A possible approach might be to simulate sampling from a nonnegative distribution that approximates a slightly negative quasidistribution, perhaps using techniques like those recently introduced for discrete-variable systems [40]. We leave this as a subject for future research.
Several lessons might be drawn from our work in this paper. First, in any protocol that uses probabilistic state preparation, the state preparation should be included when one searches for efficient classical simulations. If classical simulation is possible for sampling from the whole distribution, then it is also possible for sampling from a subdistribution that is chosen by postselection. This is the approach we used in our analysis of SPDC boson sampling, where we included the heralding modes explicitly in the search for a classical simulation. Second, our random-count thresholds are hard boundaries. These hard boundaries might be moved closer to the ideal problem by considering approximate sampling, as discussed above, but the point here is that such hard boundaries might not be found by considering perturbations about an ideal protocol. This might be a general property of analogue quantum protocols like boson sampling. A third lesson is that it is generally easier to devise analogue quantum protocols than it is to show that the protocol does not have an efficient classical simulation. Confronted with a new analogue quantum protocol, the responsibility of theorists and experimenters alike is to put on the classical thinking cap and to focus on whether classical simulations are possible in the presence of noise; this is essential for designing experiments that are meaningful implementations of the quantum protocol.
Our methods for classical simulation are based on the wave aspects of boson-sampling experiments, as opposed to the particle aspects. One classical analogue of boson sampling replaces the identical input bosons with classical distinguishable particles undergoing probabilistic transitions within a network; in this situation, output probabilities are given by permanents of matrices with nonnegative elements [6]. In contrast, in our methods, we deal with waves undergoing interference within a LON and try to mimic quantum mechanics by using quasiprobability distributions to translate from particle inputs and particle measurements to the complex amplitudes of interfering waves. This is a natural way to try to simulate an analogue quantum protocol like boson sampling. We close by noting that the mode-mismatched photons that make their way to the detectors, which we identify as the chief challenge for boson sampling, are effectively the distinguishable photons of a particle description.