Inefficiency of classically simulating linear optical quantum computing with Fock-state inputs
Bryan T. Gard, Jonathan P. Olson, Robert M. Cross, Moochan B. Kim, Hwang Lee, Jonathan P. Dowling
aa r X i v : . [ qu a n t - ph ] F e b Inefficiency of classically simulating linear optical quantum computing with Fock-stateinputs
Bryan T. Gard, ∗ Jonathan P. Olson, Robert M. Cross, Moochan B. Kim, Hwang Lee, and Jonathan P. Dowling
1, 3 Department of Physics & Astronomy, Hearne Institute for Theoretical Physics,Louisiana State University, 202 Nicholson Hall, Baton Rouge, Louisiana 70803, USA Department of Physics & Astronomy, University of Rochester,P.O. Box 270171, Rochester, New York 14627, USA Computational Science Research Center, No.3 HeQing Road, Beijing, China (Dated: November 2, 2018)Aaronson and Arkhipov recently used computational complexity theory to argue that classicalcomputers very likely cannot efficiently simulate linear, multimode, quantum–optical interferometerswith arbitrary Fock–state inputs [S. Aaronson and A. Arkhipov, Theory Comput., , 143 (2013)].Here we present an elementary argument that utilizes only techniques from quantum optics. Weexplicitly construct the Hilbert space for such an interferometer and show that its dimension scalesexponentially with all the physical resources. We also show in a simple example just how theSchr¨odinger and Heisenberg pictures of quantum theory, while mathematically equivalent, are notin general computationally equivalent. Finally, we conclude our argument by comparing the symme-try requirements of multi–particle bosonic to fermionic interferometers and, using simple physicalreasoning, connect the non–simulatability of the bosonic device to the complexity of computing thepermanent of a large matrix. I. INTRODUCTION
There is a history of attempts to use linear quan-tum interferometers to design a quantum computer.ˇCern´y showed that a linear interferometer could solveNP–complete problems in polynomial time but only withan exponential overhead in energy [1]. Clauser and Dowl-ing showed that a linear interferometer could factor largenumbers in polynomial time but only with exponentialoverhead in both energy and spatial dimension [2]. Cerf,Adami, and Kwiat showed how to build a programmablelinear quantum optical computer but with an exponentialoverhead in spatial dimension [3].Nonlinear optics provides a well–known route to uni-versal quantum computing [4]. We include in this nonlin-ear class the so–called “linear” optical approach to quan-tum computing [5], because this scheme contains an ef-fective Kerr nonlinearity [6].In light of these results there arose a widely held beliefthat linear interferometers alone, even with nonclassicalinput states, cannot provide a road to universal quan-tum computation and, as a corollary, that all such de-vices can be efficiently simulated classically. However,recently Aaronson and Arkhipov (AA) gave an argumentthat multimode, linear, quantum optical interferometerswith arbitrary Fock–state photon inputs very likely couldnot be simulated efficiently with a classical computer [7].Their argument, couched in the language of quantumcomputer complexity class theory, is not easy to followfor those not skilled in that art. Nevertheless, White andcollaborators, and several other groups, carried out ex- ∗ [email protected] FIG. 1. Quantum pachinko machine for numerical depth L = 3. We indicate an arbitrary bosonic dual–Fock in-put | N i | M i at the top of the interferometer and then thelattice of beam splitters ( B ), phase shifters ( ϕ ), and pho-ton–number–resolving detectors ( ∇ ). The vacuum inputmodes | i (dashed lines) and internal modes | ψ i (solid lines)are also shown. The notation is such that the superscripts la-bel the level ℓ and the subscripts label the row element fromleft to right. periments that demonstrated that the conclusion of AAholds up for small photon numbers [8–11]. Our goal hereis to understand—from a physical point of view—whysuch a device cannot be simulated classically.In their paper, AA prove that both strong and weaksimulation of such an interferometer is not efficient clas-sically. In the context of Fock-state interferometers, astrong simulation implies the direct computation of thejoint output probabilities of a system. However, one canconsider a “weak” simulation where one could efficientlyestimate the joint output probabilities to within someacceptably small margin of error. There are many exam-ples of systems for which weak simulation is efficient evenwhen strong simulation is not, such as finding the perma-nent of an n × n matrix with real, positive entries. Butas our goal is to provide the most straightforward andphysical explanation for this phenomenon, we do so onlyfor the strong case. Since many classical systems cannoteven be strong simulated, it may at first seem unsurpris-ing that this is the case. However, we note that here notonly does our system’s classical counterpart—Galton’sboard—admit an efficient strong simulation, but so doesa myriad of other quantum interferometers with non-Fockstate inputs as we will show.We then independently came to the same conclusionas AA in our recent analysis of multi–photon quantumrandom walks in a particular multi–mode interferome-ter called a quantum “pachinko” machine shown in Fig.1[12]. The dual–photon Fock state | N i | M i is inputtedinto the top of the interferometer and then the photonsare allowed to cascade downwards through the lattice ofbeam splitters ( B ) and phase shifters ( ϕ ) to arrive at anarray of photon–number–resolving detectors ( ∇ ) at thebottom. Our goal was to compute all the joint probabil-ities for, say, the q th detector received p photons whilethe r th detector received s photons, and so forth, for ar-bitrary input photon number and lattice dimension. Wefailed utterly. It is easy to see why.Working in the Schr¨odinger picture, we set out to com-pute the probability amplitudes at each detector by fol-lowing the Feynman–like paths of each photon throughthe machine, and then summing their contributions atthe end. For a machine of numerical depth L , as shownin Fig. 1, it is easy to compute that the number of suchFeynman–like paths is 2 L ( N + M ) . So for even a meagernumber of photons and levels the solution to the problemby this Schr¨odinger picture approach becomes rapidly in-tractable. For example, choosing N = M = 9 and L = 6,we have 2 ∼ = 5 × total possible paths, which isabout four orders of magnitude larger than the numberof atoms in the observable universe. We were puzzled bythis conclusion; we expected any passive linear quantumoptical interferometer to be efficiently simulatable classi-cally. With the AA result now in hand, we set out hereto investigate the issue of the complexity of our quantumpachinko machine from an intuitive physical perspective.The most mathematics and physics we shall need is ele-mentary combinatorics and quantum optics.Following Feynman, we shall explicitly construct thepachinko machine’s Hilbert state space for an arbitrarylevel L , and for arbitrary photon input number, and showthat the space’s dimension grows exponentially as a func-tion of each of the physical resources needed to build andrun the interferometer [13]. Because interference only oc-curs when the input state has been symmetrized (withrespect to interchange of mode), we compute the size ofthe symmetrized subspace and show that it too growsexponentially with the number of physical resources. Weremark that while a classical pachinko machine (or “Gal-tons board”) will also have an exponential large state space, because no interference occurs there is only aquadratic increase with L in the number of calculationsnecessary to simulate the output (corresponding to thenumber of beam splitters in the interferometer). Fromthis result we conclude that it is very likely that any clas-sical computer that tries to simulate the operation of thequantum pachinko machine will always suffer an expo-nential slowdown. We will also show that no exponentialgrowth occurs if Fock states are replaced with photoniccoherent states or squeezed states, which elucidates partof the special nature of photonic Fock states. Howeveran exponentially large Hilbert space, while necessary forclassical non–simulatability, is not sufficient. We thenfinally examine the physical symmetry requirements forbosonic versus fermionic multi–particle states and showthat in the bosonic case, in order to simulate the interfer-ometer as a physics experiment, one must compute thepermanent of a large matrix, which in turn is a problem ina computer algorithm complexity class strongly believedto be intractable on a classical or even a quantum com-puter. This concludes our elementary argument, whichinvokes only simple quantum mechanics, combinatorics,and a simplistic appeal to complexity theory. II. THE PACHINKO MACHINE MODEL
As our argument is all about counting resources, wehave carefully labeled all the components in the pachinkomachine in Fig. 1 to help us with that reckoning. Themachine has a total of L levels of physical depth d each. The input state at the top is the dual–Fock state | N i | M i , where the superscripts label the level num-ber and the subscripts the element in the row at thatlevel (from left to right). We illustrate a machine of to-tal numerical depth of L = 3. For 1 ≤ ℓ < L , we showthe vacuum input modes along the edges of the machine.The resources we are most concerned about are energy,time, spatial dimension, and number of physical elementsneeded to construct the device. All of these scale eitherlinearly or quadratically in either L or N + M . The to-tal physical depth is D = Ld and so the spatial area is A = ( √ D ) = 2 L d . Using identical photons of fre-quency ω , the energy per run is E = ( N + M ) ~ ω . Thetime it takes for the photons to arrive at the detectors is T = √ Ld/c , where c is the speed of light. In each levelthe photons encounter ℓ number of beam splitter (BS) sothe total number is B = P Lℓ =1 ℓ = L ( L + 1) /
2. Beloweach BS (with the exception of the L th level) there aretwo independently tunable phase shifters (PS) for a totalnumber of PS that is ϕ = P L − ℓ =1 ℓ = L ( L − ∇ = 2 L . The total num-ber of input modes is equal to the total number of outputmodes and is I = O = 2 L . The total number of in-ternal modes is ψ = P L − ℓ =1 ℓ = L ( L − L or N + M .The input state may be written in the Heisenbergpicture as | N i | M i = (ˆ a † ) N (ˆ a † ) M | i | i / √ N ! M !,where ˆ a † is a modal creation operator. Each BS per-forms a forward unitary mode transformation, which weillustrate with B , of the form ˆ a = ir ˆ a + t ˆ a andˆ a = t ˆ a + ir ˆ a where the reflection and transmissioncoefficients r and t are positive real numbers such that r + t = R + T = 1. The choice r = t = 1 / √ iϕ ˆ n ) onmode | ψ i , where ˆ n := ˆ a † ˆ a is the number operator, ˆ a is the annihilation operator conjugate to ˆ a † , and ϕ isa real number. Finally the 2 L detectors in the final level L are each photon number resolving [14].To argue that this machine (or any like it) cannotbe simulated classically, in general, it suffices to showthat this is so for a particular simplified example. Wenow take N and L arbitrary but M = 0 and turn offall the phase shifts and make all the BS identical bysetting ϕ ℓk = 0, t ℓk = t , and r ℓk = r for all ( k, ℓ ).We then need the backwards BS transformation on thecreation operators, which is, ˆ a † = ir ˆ a † + t ˆ a † andˆ a † = t ˆ a † + ir ˆ a † . Similar transforms apply down themachine at each level. With M = 0 the input simpli-fies to | N i | i = (ˆ a † ) N | i | i / √ N ! and now we ap-ply the first backwards BS transformation | ψ i | ψ i =( ir ˆ a † + t ˆ a † ) N | i | i / √ N ! to get the state at levelone.At every new level each ˆ a † will again bifurcate accord-ing to the BS transformations for that level, with thetotal number of bifurcations equal to the total numberof BS, and so the computation of all the terms at thefinal level involves a polynomial number of steps in L . Itis instructive to carry this process out explicitly to level L = 3 to get, | ψ i = 1 √ N ! ( irt ˆ a † − r t ˆ a † + ir ( t − r )ˆ a † − r t ˆ a † + irt ˆ a † + t ˆ a † ) N Y ℓ =1 | i ℓ , (1)where we have used a tensor product notation for thestates. If r ∼ = 0 or r ∼ = 1 the state is easily computed.Since we are seeking a regime that cannot be simulatedclassically we work with r ∼ = t ∼ = 1 / √ III. SOLUTION IN THE HEISENBERG ANDSCHR ¨ODINGER PICTURES
It is now clear from Eq.(1) what the general form ofthe solution will be. We define | ψ i L := X { n ℓ } N = P Lℓ =1 n ℓ | ψ i L { n ℓ } , | i L := L Y ℓ =1 | i Lℓ , (2) and the general solution has the form, | ψ i L = 1 √ N ! L X ℓ =1 α Lℓ ˆ a † Lℓ ! N | i L = 1 √ N ! X N = P Lℓ =1 n ℓ (cid:18) Nn , n , . . . , n L (cid:19) × Y ≤ k ≤ L ( α Lk ˆ a † Lk ) n k | i L , (3)where all the coefficients α Lℓ will be nonzero in general.Since all the operators commute, as they each operateon a different mode, we have expanded Eq. (3) usingthe multinomial theorem where the sum in the expan-sion is over all combinations of non–negative integers con-strained by N = P Lℓ =1 n l and (cid:18) Nn , n , . . . , n L (cid:19) = N ! n ! n ! . . . n L ! (4)is the multinomial coefficient [15]. The state | ψ i L ishighly entangled over the number–path degrees of free-dom. Each monomial in the expansion of Eq. (3) isunique and so the action of the set of all monomial oper-ators on the vacuum will produce a complete orthonor-mal basis set for the Hilbert space at level L , given by | ψ i L { n ℓ } := Q Lℓ =1 | n ℓ i Lℓ , where the n ℓ are subject to thesame sum constraint. Let us call the dimension of thatHilbert space dim[ H ( N, L )], which is therefore the totalnumber of such basis vectors.Taking L = 3 and N = 2, we can use Eq.(3) to com-pute the probability a particular sequence of detectorswill fire with particular photon numbers. What is theprobability detector one gets one photon, detector twoalso gets one, and all the rest get zero? This is the mod-ulus squared of the probability amplitude of the state | i | i | i | i | i | i . Setting r = t = 1 / √ α = irt = i/ (2 √
2) and α = − r t = − / (2 √ P ∼ = 0 . L and N ) tocompute the single and binary joint probabilities, thatdetector p gets n photons and detector q gets m [16].However computing arbitrary joint probabilities betweentriplets, quadruplets, etc., of detectors rapidly becomesintractable. We can provide a closed form expressionfor dim[ H ( N, L )] by realizing that it is the same as thenumber of different ways one can add up non–negativeintegers that total to fixed N . More physically this is thenumber of possible ways that N indistinguishable pho-tons may be distributed over 2 L detectors. The answeris well known in the theory of combinatorics and is:dim[ H ( N, L )] = (cid:18) N + 2 L − N (cid:19) , (5)where this is the ordinary binomial coefficient [17]. Forour example with L = 3, N = 2, Eq.(5) implies that thenumber of distinct probabilities P npqrst to be tabulatedis again 21.We first examine two “computationally simple” ex-amples. Taking N arbitrary and L = 1 we getdim[ H ( N, N + 1, which is easily seen to be thenumber of ways to distribute N photons over two de-tectors. Next taking N = 1 and L arbitrary we getdim[ H (1 , L )] = 2 L , which is the number of ways to dis-tribute a single photon over 2 L detectors. If we wereto invoke Dirac’s edict—“Each photon then interferesonly with itself.”—we would then expect that addinga second photon should only double this latter result[18]. Instead the effect of two–photon interference onthe state space can be seen immediately by computingdim[ H (2 , L )] = L (2 L + 1). That is, adding a second pho-ton causes a quadratic (as opposed to linear) jump in thesize of the Hilbert space. Dirac was wrong; photons dointerfere with each other, and that multiphoton interfer-ence directly affects the computational complexity. Allthese three cases are simulatable in polynomial time stepswith N and L , but we see a quadratic jump in dimensionas soon as we go from one to two photons. These jumpsin complexity continue for each additional photon addedand the dimension grows rapidly.We therefore next investigate a “computationally com-plex” intermediate regime by fixing N = 2 L − L and then choose an odd–numbered photon inputso that this restriction holds. Equation (5) becomesdim[ H ( N )] = (2 N )! / ( N !) . Deploying Sterling’s approx-imation for large N , in the form n ! ∼ = ( n/e ) n √ πn wehave dim[ H ( N )] ∼ = 2 N / √ πN . This is one of our pri-mary results. The Hilbert space dimension scales expo-nentially with N = 2 L −
1. Since all the physical param-eters needed to construct and run our quantum pachinkomachine scale only linear or quadratically with respectto N or L , we have an exponentially large Hilbert spaceproduced from a polynomial number of physical resources—Feynman’s necessary condition for a potential univer-sal quantum computer.Let us suppose we build onto an integrated optical cir-cuit a machine of depth L = 69 and fix N = 2 L − H (137)] = 10 , which is again on the order ofthe number of atoms in the observable universe. Follow-ing Feynman’s lead, we conclude that, due to this expo-nentially large Hilbert space, we have a sufficient condi-tion that a classical computer can not likely efficientlysimulate this device. However this is not a necessarycondition. From the Gottesman-Knill theorem we knowthat quantum circuits that access an exponentially largeHilbert space may sometimes be efficiently simulated [20].We will strengthen our argument (below) by discussingthe necessity of properly symmetrizing a multi-particlebosonic state and tie that physical observation back tothe complexity of computing the permanent of a largematrix. Let us now compare our Heisenberg picture result tothat of the Schr¨odinger picture. In the computationallycomplex regime where N = 2 L − LN = 2 N ( N +1) / ∼ = 2 N / . Taking N = 137and L = 69, as in the previous example, we get an as-tounding 2 ∼ = 4 × total paths. Dirac provedthat the Heisenberg and Schr¨odinger pictures are math-ematically equivalent, that they always give the samepredictions, but we see here that they are not alwaysnecessarily computationally equivalent [21]. Calculationsin the Heisenberg picture are often much simpler than inthe Schr¨odinger picture. The fact that the two picturesare not always computationally equivalent is implicit inthe Gottesman–Knill theorem; however, it is satisfyingto see here just how that is so in a simple optical inter-ferometer [20]. IV. SAMPLING WITH COHERENT &SQUEEZED STATE INPUTS
To contrast this exponential overhead from the re-source of bosonic Fock states, let us now carry out thesame analysis with the bosonic coherent input state in-put | β i | i , where we take the mean number of photonsto be | β | = n . In the Heisenberg picture this input be-comes ˆ D ( β ) | i | i , where ˆ D ( β ) = exp( β ˆ a † − β ∗ ˆ a )is the displacement operator [22]. Applying the BS trans-formations down to final level L we get | ψ i L = exp β L X ℓ =1 α Lℓ ˆ a † Lℓ − β ∗ L X ℓ =1 α L ∗ ℓ ˆ a Lℓ ! | i L = L Y ℓ =1 exp( βα Lℓ ˆ a † Lℓ − β ∗ α L ∗ ℓ ˆ a Lℓ ) | i L = L Y ℓ =1 | βα Lℓ i Lℓ . (6)At the output we have 2 L coherent states that have beenmodified in phase and amplitude. This is to be expected,as it is well known that linear interferometers transforma coherent state into another coherent state. Since allthe coefficients α Lℓ are computable in B = L ( L + 1) / L , independent of n . The mean number of photons ateach detector is then simply n Lℓ = | βα Lℓ | = n | α Lℓ | .A similar analysis may be carried out for bosonicsqueezed input states. Taking, for example, asingle–mode squeezed vacuum input | ξ i | i =ˆ S ( ξ ) | i | i , with the squeezing operator defined asˆ S ( ξ ) = exp { [ ξ ∗ (ˆ a ) − ξ (ˆ a † ) ] / } , we arrive at, | ψ i L = exp ξ ∗ L X ℓ =1 α ∗ ℓ ˆ a Lℓ ! − ξ L X ℓ =1 α ℓ ˆ a † Lℓ ! / | i L , (7)which does not in general decompose into a separableproduct of single–mode squeezers on each output port.Nevertheless the probability amplitudes may still be com-puted in a time polynomial in L by noting that, fromEq.(5) with N = 2, there are at most 2 L ( L + 1) termsin this exponent that must be evaluated. This resultgeneralizes to arbitrary Gaussian state inputs [23]. Theoutput of the interferometer may be then calculated onthe transformed device in polynomial steps in L .The exponential scaling comes from the bosonic Fockstructure | N i = (ˆ a † ) N | i / √ N ! and the rapid growthof the number–path entanglement in the interferome-ter. It is well known that beam splitters can gener-ate number–path entanglement from separable bosonicFock states. For example, the simplest version of theHOM effect at level one with separable input | i | i becomes | ψ i | ψ i = ( i ˆ a † + ˆ a † )(ˆ a † + i ˆ a † ) | i | i / i [ | i | i + | i | i ] / √ V. COMPARISON OF BOSONIC TOFERMIONIC FOCK STATE INPUTS
We now compare the multimode bosonic Fock state in-terferometer to the multimode fermionic interferometer.We will restrict ourselves to spin–1/2 neutral fermionssuch as neutrons that are commonly used in interferom-etry. Now the number of fermions per input mode isrestricted to zero, one, or two and we can have two onlyif they have opposite spin states to be consistent with thePauli exclusion principle. The exclusion principle is de-rived from the requirement that the total multi–particlefermionic wave function, which is the product of thespin and spatial wave functions, is antisymmetric underthe interchange of any two particle state labels. Like-wise there is a constraint on the bosonic multi–particlemulti–mode wave function that the total wave functionbe symmetric. The symmetry of the wave function mustbe enforced at each beam splitter where the particles be-come indistinguishable and the spatial part of the wavefunction experiences maximal overlap for multi–particleinterference to occur. For the sake of argument we takethe coherence length of the particles to be infinite (orat least much larger than the depth of the interferome-ter Ld ) so that enforcing the correct symmetry at eachbeam splitter requires enforcing the correct symmetryeverywhere in space.Some care must now be used in the notation. For ex-ample when we write the bosonic spatial wave function input state | i bA in | i bB in , we are assuming both bosonshave the same spin state, since clearly this state isspatially symmetric under particle interchange its spinstate must also be, so that the product of the two(total wave function) remains symmetric. To denotethis point we instead write |↑i bA in |↑i bB in to explicitlyshow the spin state. [More properly we should write ψ b ( x A in ) ψ b ( x B in ) |↑i bA in |↑i bB in but this notation is a bitcumbersome.] Thence for a 50:50 BS the HOM ef-fect for bosons in the same spin state can be written, |↑i bA in |↑i bB in BS → |↑↑i bA out | i bB out + | i bA out |↑↑i bB out , so bothbosons “stick” at the beam splitter and emerge together.This effect arises as a direct result of the fact that the spatial part of the wave function, which gives rises toan effective attraction at the BS, is symmetric. Wecould instead prepare an antisymmetric bosonic singletspin state input |↑i bA in |↓i bB in − |↓i bA in |↑i bB in , in which casethe spatial wavefunction must be also antisymmetric, ψ b ( x A in ) ψ b ( x B in ) − ψ b ( x B in ) ψ b ( x A in ) , so that the prod-uct of the two remains symmetric. In this case the parti-cles behave fermionically as far as the spatial wavefunc-tion overlap is concerned at the BS and they repel eachother in an anti–HOM effect, always exiting out separateports and never together; |↑i bA in |↓i bB in − |↓i bA in |↑i bB in BS →|↑i bA out |↓i bB out − |↓i bA out |↑i bB out [26, 27]. The reverse hap-pens for fermions.For example, the symmetric spin input state |↑i fA in |↑i fB in is allowed for fermions only if the spatialwave wave function is antisymmetric, ψ f ( x A in ) ψ f ( x B in ) − ψ f ( x B in ) ψ f ( x A in ), so that the entire wave function prod-uct remains antisymmetric. Since the spatial partgoverns the HOM effect they repel at the BS andobey an anti–HOM effect and always exit out separateports, consistent with the exclusion principle, namely |↑i fA in |↑i fB in BS → |↑i fA out |↑i fB out . However, we can makethe fermions behave spatially bosonically by preparingthem in a spin–antisymmetric singlet input state, whichthen must be symmetric in the spatial part, and so theybehave as bosons as far as the spatial overlap is con-cerned, and we recover the usual HOM effect, where nowthey always exit the same port together: |↑i fA in |↓i fB in −|↓i fA in |↑i fB in → |↑↓i fA out | i fB out − |↓↑i fA out | i fB out . There isno violation of the exclusion principle as they also alwaysexit with opposite spins. (This type of effective spatial at-traction between fermions in a spin singlet state explainswhy the ground state of the neutral hydrogen moleculeis a bound state.) It is clear then that even fermions canexperience number-path entanglement in a linear inter-ferometer, although not to the same degree as bosons.However this entanglement is still sufficient to lead to anexponential growth in the fermionic Hilbert space, as weshall now argue.Now we are ready to apply our resource counting ar-gument to the fermionic case. For fermions the com-putationally complex regime may be accessed when thenumber of input particles N is half the number of inputmodes 2 L . The dimension of the Hilbert space may becomputed as before and turns out to be, for this exam-ple, (cid:0) LN (cid:1) . This also grows exponentially as a functionof the resources choosing N = L . Following the sameSterling’s approximation argument as above we get ex-actly the same exponential formula for the Hilbert spacedimension as with bosons, namely 2 N / √ πN .So, in general, in both the fermionic and bosoniccase the Hilbert space dimension grows exponen-tially with respect to the resources: particle num-ber and mode number. However, Feynman’s argu-ments notwithstanding, an exponential growth in theHilbert space is only sufficient but not necessary toattain classical non–simulatability. For example, fromthe Gottesman–Knill theorem, we can construct a Clif-ford–algebra–based quantum computer circuit that ac-cesses an exponentially large Hilbert space but still canbe simulated efficiently classically [20]. Sometimes thereare shortcuts through Hilbert space, as we shall now ar-gue is the case here for fermions but not for bosons.In order to access these large Hilbert spaces in theinterferometer one must require that multi–particle in-terference take place at each beam splitter, where theparticles must be indistinguishable, and the spatial wavefunction overlap determines the type of particle–modeentanglement that will result. The overall bosonic wavefunction (spatial multiplied by spin) must be totally sym-metric and the overall fermionic wave function must betotally antisymmetric at each row of BS, and so theymust have these symmetries everywhere in space and par-ticularly at the input. Now if we give up on a completetabulation of the Hilbert state space at level L , due toits exponential growth, and treat the interferometer us-ing a standard quantum optical input–output formalism,there is an efficient way to take a given multi–particle,multi–mode input state at the top of the interferome-ter to the bottom of the interferometer. This method iscalled matrix transfer and is accomplished by encodingeach level of BS transformations in terms of L matricesof size (2 L ) × (2 L ) and then multiplying them together.This can be done in the order of O ( L ) steps and so it isefficient.We must now address the issue of computing the sam-pling output of the interferometer. While the one- andtwo- particle joint particle detection probabilities at thedetectors may be computed efficiently, computation ofthe higher order joint probabilities rapidly become in-tractable [16]. In order to compute the complete jointprobability distribution, we must compute the determi-nant (if the input is fermionic) or the permanent (if theinput is bosonic) of the (2 L ) × (2 L ) matrix found above.Using the method of Laplace decomposition for con-structing the determinant of a matrix, one decomposesthe large determinant into a sum over ever-smaller deter-minants, appending alternating plus and minus signs toeach in a checkerboard pattern. Constructing the perma-nent follows the same process but all the signs are nowonly plus. However, for the determinant, there is a polynomialshortcut through the exponential Hilbert space—therow–reduction method. For fermions we may always con-struct the most general input state efficiently. On theother hand, there is no known method such as row re-duction to compute the permanent of an arbitrary ma-trix efficiently. The most efficient known protocols forthe permanent computation are variants on the Laplacedecomposition and all scale exponentially with the size ofthe matrix. This problem of computing the permanent,in the lingo of computer complexity theory, is that it isin the class of “ O (2 L L ) number of steps [28]. Finally, we have reachedthe snag that undermines our ability to efficiently com-pute the output and so renders simulation of the deviceclassically intractable. VI. CONCLUSION
In conclusion, we have shown that a multi–mode linearoptical interferometer with arbitrary Fock input statesis very likely not simulatable classically. Our result isconsistent with the argument of AA. Without invokingmuch complexity theory, we have argued this by explic-itly constructing the Hilbert state space of a particularsuch interferometer and showed that the dimension growsexponentially with the size of the machine. The outputstate is highly entangled in the photon number and pathdegrees of freedom. We have also shown that simulat-ing the device has radically different computational over-heads in the Heisenberg versus the Schr¨odinger picture,illustrating just how the two pictures are not in generalcomputationally equivalent within this simple linear op-tical example. Finally we supplement our Hilbert spacedimension argument with a discussion of the symmetryrequirements of multi–particle interferometers and par-ticularly tie the simulation of the bosonic device to thecomputation of the permanent of a large matrix, whichis strongly believed to be intractable. It is unknown (butthought unlikely) if such bosonic multi–mode interfer-ometers as these are universal quantum computers, butregardless they will certainly not be fault tolerant. Aspointed out by Rohde [29], it is well known that Fockstates of high photon number are particularly sensitiveto loss [30]. They are also super–sensitive to dephas-ing as well [31]. This implies that even if such a machineturns out to be universal it would require some type of er-ror correction to run fault tolerantly. Nevertheless, suchdevices could be interesting tools for studying the rela-tionship between multi–photon interference and quantuminformation processing for small numbers of photons. Ifwe choose each of the PS and BS transformations inde-pendent of each other, we have a mechanism to programthe pachinko machine by steering the output into anyof the possible output states. Even if universality turnsout to be lacking we may very well be able to exploitthis programmability to make a special purpose quantumsimulator for certain physics problems such as frustratedspin systems [32].
ACKNOWLEDGMENTS
B.T.G would like to acknowledge support from theNPSC/NIST fellowship program. J.P.D. would like to ac- knowledge support from the NSF. This work is also sup-ported by the Intelligence Advanced Research ProjectsActivity (IARPA) via Department of Interior NationalBusiness Center contract number D12PC00527. The U.S.Government is authorized to reproduce and distributereprints for governmental purposes notwithstanding anycopyright annotation thereon. Disclaimer: The viewsand conclusions contained herein are those of the authorsand should not be interpreted as necessarily representingthe official policies or endorsements, either expressed orimplied, of IARPA, DoI/NBC, or the U.S. Government.We also would like to acknowledge interesting and use-ful discussions with S. Aaronson, S. T. Flammia, K. R.Motes, and P. P. Rohde. [1] V. ˇCern´y, Phys. Rev. A , 116 (1993).[2] J. F. Clauser and J. P. Dowling, Phys. Rev. A , 4587(1996).[3] N. J. Cerf, C. Adami, and P. G. Kwiat, Phys. Rev. A , R1477 (1998).[4] G. J. Milburn, Phys. Rev. Lett. , 2124 (1989).[5] E. Knill, R. LaFlamme, and G. J. Milburn, Nature (Lon-don) , 46 (2001).[6] G. G. Lapaire, P. Kok, J. P. Dowling, and J. E. Sipe,Phys. Rev. A , 042314 (2003).[7] S. Aaronson and A. Arkhipov, Theory Comput. , 143(2013).[8] M. A. Broome, A. Fedrizzi, S. Rahimi-Keshari, J. Dove,S. Aaronson, T. C. Ralph, and A. G. White, Science , 794 (2013).[9] A. Crespi, R. Osellame, R. Ramponi, D. J. Brod, E. F.Galvao, N. Spagnolo, C. Vitelli, E. Maiorino, P. Mat-aloni, and F. Sciarrino, Nat. Photon. , 545 (2013).[10] M. Tillmann, B. Daki´c, R. Heilmann, S. Nolte, A. Sza-meit, and P. Walther, Nat. Photon. , 540 (2013).[11] J. Spring, B. J. Metcalf, P. C. Humphreys, W. S.Kolthammer, X. Jin, M. Barbieri, A. Datta, N. Thomas-Peter, N. K. Langford, D. Kundys, et al. , Science ,798 (2013).[12] B. T. Gard, R. M. Cross, P. M. Anisimov, H. Lee, andJ. P. Dowling, J. Opt. Soc. Am. B , 1538 (2013. Theparticular triangular arrangement of pins in Fig. 1 is alsocalled a quantum quincunx [B. C. Sanders, S. D. Bartlett,B. Tregenna, et al. , Phys. Rev. A67, 042305 (2003)]).[13] R. P. Feynman, J. Opt. Soc. Am. B , 464 (1984).[14] H. Lee, U. Yurtsever, P. Kok, G. Hockney, C. Adami,S. L. Braunstein, and J. P. Dowling, J. Mod. Opt. ,1517 (2004).[15] NIST Digital Library of Mathematical Functions (Release1.0.5 of 2012–10–01), § Sec 26.4, http://dlmf.nist.gov/ (2012).[16] K. Mayer, M. C. Tichy, F. Mintert, T. Konrad, andA. Buchleitner, Phys. Rev. A , 062307 (2011).[17] A. Benjamin, A. T. Benjamin, and J. J. Quinn, Proofsthat Really Count: The Art of Combinatorial Proof (Mathematical Association of America, Oberlin, OH,,2003) p. 71.[18] P. A. M. Dirac,
The Principles of Quantum Mechanics,4th ed. (Clarendon Press, London, 1958) p. 9.[19] D. Bonneau, E. Engin, K. Ohira, N. Suzuki, H. Yoshida,N. Iizuka, M. Ezaki, C. M. Natarajan, M. G. Tanner,R. H. Hadfield, et al. , New J. Phys. , 045003 (2012).[20] D. Gottesman, in Group22: Proceedings of the XXII In-ternational Colloquium on Group Theoretical Methods inPhysics , edited by S. P. Corney, R. Delbourgo, and P. D.Jarvis (International Press,Cambridge, MA, 1999) pp.32–43.[21] P. A. M. Dirac, Proc. R. Soc. A , 621 (1927).[22] C. C. Gerry and P. L. Knight,
Introductory QuantumOptics. (Cambridge University Press, Cambridge, UK,2005).[23] S. D. Bartlett, B. C. Sanders, S. L. Braunstein, andK. Nemoto, Phys. Rev. Lett. , 097904 (2002); V. Ve-icht, N. Wiebe, and C. Ferrie, New J. Phys. , 013037(2013).[24] Z. Y. Ou, C. K. Hong, and L. Mandel, Opt. Commun. , 118 (1987).[25] C. F. Wildfeuer, A. P. Lund, and J. P. Dowling, Phys.Rev. A , 052101 (2007).[26] R. Loudon, in Coherence and Quantum Optics VI , editedby J. H. Eberly, L. Mandel, and E. Wolf (Plenum, NewYork, 1989) p. 703.[27] R. Loudon, Phys. Rev. A , 4904 (1998).[28] H. J. Ryser, Combinatorial Mathematics, The CarusMathematical Monographs (Mathematical Association ofAmerica, Oberlin, OH, 1963).[29] P. P. Rohde, arXiv:1010.4608.[30] S. D. Huver, C. F. Wildfeuer, and J. P. Dowling, Phys.Rev. A , 063828 (2008).[31] A. A. Qasimi and D. F. V. James, Opt. Lett. , 268(2009).[32] J. W. Britton, B. C. Sawyer, A. C. Keith, C. –C. J. Wang,J. K. Freericks, H. Uys, M. J. Biercuk, and J. J.Bollinger, Nature (London)484