Investigating the optimality of ancilla-assisted linear optical Bell measurements
IInvestigating the optimality of ancilla-assisted linear optical Bell measurements
Andrea Olivo
1, 2, ∗ and Frédéric Grosshans † Inria, Paris, France Laboratoire Aimé Cotton, CNRS, Université Paris-Sud,ENS Cachan, Université Paris-Saclay, 91405 Orsay Cedex, France (Dated: September 5, 2018)In the last decade Grice [1] and Ewert and van Loock [2] found linear optical networks achievingnear-unit efficiency unambiguous Bell state discrimination, when fed with increasingly complexancillary states. However, except for the vacuum ancilla case [3], the optimality of these schemes isunknown. Here, the optimality of these networks is investigated through analytical and numericalmeans. We show an analytical upper bound to the success probability for interferometers thatpreserve the polarization of the input photons, saturated by both Grice’s and Evert-van Loock’sstrategies. Furthermore, such an upper bound links the complexity of their ancilla states with thescaling of their performance. We also show a computer-aided approach to the optimization of suchmeasurement schemes for generic interferometers, by simulating an optical network supplied withvarious kinds of ancillary input states. We numerically confirms the optimality of known smallschemes. We use both methods to investigate other ancilla states, some of them never studied before.
I. INTRODUCTION
Due to its experimental and theoretical simplicity, linearquantum optics has proved to be a promising route for theearly implementation of important quantum communica-tion protocols [4]—including quantum teleportation [5–7],dense coding [8, 9] and entanglement swapping [10, 11].An essential step in these protocols is the
Bell measure-ment (BM), a projective measurement onto a basis oftwo-qubit maximally entangled states, the Bell states. Inorder to enable the common scenario where losses can betolerated but not errors, in the following we will be con-cerned with unambiguous
BM, i.e. a measurement whoseoutcome is never wrong, but is sometimes inconclusive.Lütkenhaus et al. [12] showed long ago the impossibilityof a linear optical perfect Bell measurement for dual-railphotonic qubits. In a following result, Calsamglia andLütkenhaus bound the success probability P succ of theno-ancilla case to 50% [3], thus proving the optimalityof the already known Braunstein–Mann scheme [4]. Re-cently, Grice [1], followed by Ewert and van Loock [2]showed how to overcome this bound by supplying thenetwork with ancillary states. They showed how to at-tain a success probability of P succ = 3 / with reasonableancillary states, and how to increase this probability tovalues arbitrarily close to 1, by using increasingly complexand entangled ancillæ. Other ways around the 50% bar-rier include feed-forward techniques from linear opticalquantum computation [13, 14], squeezing operations [15],Kerr nonlinearities [16], entangled coherent states, hy-brid entanglement [6, 17] or encoded qubits [18]. Whilesome of the above techniques may in principle be usedto realize a perfect BM, each one has its own disadvan-tages and present different experimental challenges in ∗ [email protected] † [email protected] their implementation. For linear optics the difficultiesare concentrated in the preparation of complex ancillarystates; this is partly compensated by the simplicity ofinterferometers.Our main motivation in this work is to investigate theoptimality of ancilla-assisted linear optical schemes: whatis the highest possible P succ for a given ancilla? Which arethe “simplest” states to achieve a given value for P succ ?Very recently, and independently from us, Smith andKaplan [19] tackled a similar problem, optimizing themutual information of a Bell measurement using singlephoton ancillæ.In Section II, we first present an analytical upperbound to the success probability of ancilla-assisted BM forpolarization-preserving interferometers. Specifically, thisbound is saturated by the Grice scheme [1]. In Section III,we expose a linear optical network optimizer based onsymbolic and numerical computations, built in order tomaximize the success probability of unambiguous BM fora given ancillary state and to argue about its optimality.This program is available as supplementary material tothis article [20]. In Section IV, we discuss the resultsobtained with the two approaches above with differentkinds of ancillæ and compare them to previously knownresults [1–4, 19]. To our knowledge, some ancillæ studiedhere were never employed for this task. We also discusssome of the hidden symmetries of the problem at hand,some of which we exploit in order to lessen the amountof computation needed. Finally, we conclude this work inSection V. II. ANALYTICAL UPPER BOUND FORPOLARIZATION-PRESERVINGINTERFEROMETERS
A desirable goal is to find an unambiguous Bell mea-surement using linear optics with the best possible successprobability P succ and the simplest possible ancillary state. a r X i v : . [ qu a n t - ph ] S e p In this section, we present an analytical upper bound forthe success probability P succ of such unambiguous discrim-ination schemes using ancillary states and polarizationpreserving linear optics. This restriction to polarizationpreserving interferometers— i.e. networks described by ablock-diagonal unitary U = diag( U h , U v ) , with U h (resp. U v ) acting on the horizontally (resp. vertically) polarizedmodes—is not motivated by experimental realities but isan artificial consequence of the proof technique. Specifi-cally, this restriction will not be enforced in the computeraided optimization of Section III. However, the previouslyknown schemes [1, 2, 4] are all polarization independent, i.e. polarization preserving with U h = U v , so this restric-tion still allows to achieve non trivial discrimination. Thecommon polarization independence of these schemes wasmotivated by the symmetry of the Bell states set. Weconjecture that the optimal interferometer presents thesame symmetry, and is polarization independent, exceptfor an initial preprocessing step of the ancillæ, as for thesingle-photon schemes of [2].We present this polarization preserving upper boundfor a generic ancilla in Section II A, and bound it itselfby a simple photon number dependent expression in Sec-tion II B. A. Generic upper bound
The ability of an interferometer to perform a BM is itsability to discriminate between the four Bell states whenthey are supplied as its input, each with equal probability / . We consider, w.l.o.g. , a pure k -photon additionalresource state | Υ (cid:105) = (cid:80) kλ =0 υ λ | Υ , λ (cid:105) , where | Υ , λ (cid:105) is apure state with λ photons polarized horizontally and k − λ photons polarized vertically. Since our interferometer ispolarization preserving, the total number of horizontally(or vertically) polarized photons is unchanged, and thisnumber can be easily deduced from the observed detectionevent. Therefore, the measurement would be the same ifone had performed a projective measurement of of thesepolarized photon numbers on the global input state, beforefeeding the projected state into the interferometer.Let us rewrite the input state according to this projec-tion: (cid:12)(cid:12) Ψ ± (cid:11) | Υ (cid:105) = k +1 (cid:88) λ =1 υ λ − √ | HV (cid:105) ± | V H (cid:105) ) | Υ , λ − (cid:105) (cid:12)(cid:12) Φ ± (cid:11) | Υ (cid:105) = k +2 (cid:88) λ =0 υ λ − √ | HH (cid:105) | Υ , λ − (cid:105) ± υ λ √ | V V (cid:105) | Υ , λ (cid:105) , (1)where we have set υ λ = 0 for λ < and λ > k , in orderto include the edge cases in the formula. Each term ofthe above sums corresponds to a term with λ horizontallypolarized photons.To obtain an upper bound, we assume the measurementto be perfect after the projection onto the state with λ horizontally polarized photons and ( k + 2) − λ verticallypolarized ones, forgetting about the linear optics restric-tion. One can easily see in the above equations that, foreach λ , the term corresponding to the states | Ψ + (cid:105) , | Ψ − (cid:105) ,and | Φ ± (cid:105) are in three orthogonal subspaces, with the onlypossible remaining ambiguity being between the | Φ + (cid:105) and | Φ − (cid:105) states. Let us now look at the distinguishability ofthose states.For an arbitrary ancillary state, the unambiguous dis-crimination has to be performed between the followingunnormalized states: (cid:12)(cid:12) Λ ± (cid:11) = υ λ − √ | HH (cid:105) | Υ , λ − (cid:105) ± υ λ √ | V V (cid:105) | Υ , λ (cid:105) . Unambiguous discrimination of such a pair of pure statesis possible with an optimal success probability [21–23] of (cid:107) Λ (cid:107) − |(cid:104) Λ + | Λ − (cid:105)| , with (cid:107) Λ (cid:107) = (cid:104) Λ + | Λ + (cid:105) = (cid:104) Λ − | Λ − (cid:105) . Thisoptimal success probability is an upper bound to what isachievable with linear optics and photon counting, leadingto P succ ,λ ≤ (cid:16) | υ λ − | , | υ λ | (cid:17) . The total success probability—assuming an equal inputprobability of / for each Bell state, and perfect discrim-ination of | Ψ ± (cid:105) —is then P succ ≤ + k +2 (cid:88) λ =0 min (cid:16) | υ λ − | , | υ λ | (cid:17) = + (cid:88) λ even min (cid:16) | υ λ − | , | υ λ | (cid:17) + (cid:88) λ odd min (cid:16) | υ λ − | , | υ λ | (cid:17) . Since the minimum is taken every two value of λ , evenand odd indices have to be considered separately. It isthen easy to notice that, among even (resp. odd) valuesof λ , all values of | υ λ | are counted once except the localmaxima, which are omitted, and the local minima, whichare counted twice, leading to P succ ≤ − (cid:88) λ even | υ λ | loc max | υ λ | + (cid:88) λ even | υ λ | loc min | υ λ | − (cid:88) λ odd | υ λ | loc max | υ λ | + (cid:88) λ odd | υ λ | loc min | υ λ | . (2)This bound can then be simplified to the following looserbound on the failure probability P fail ≥ (cid:18) max λ even (cid:16) | υ λ | (cid:17) + max λ odd (cid:16) | υ λ | (cid:17)(cid:19) , (3)the latter being equivalent to (2) iff | υ λ | has a singlelocal maximum over even values of λ and a single localmaximum over odd values of λ .In Section IV, we will compute this bound for differentancillary states | Υ (cid:105) . But for now, we can already use theabove equations to compute a bound depending only on | Υ (cid:105) ’s photon number k . B. Photon-number based upper bound
Let us now work out a bound that is independent ofthe specific form of | Υ (cid:105) . If the number of photons k inthe ancilla is odd, there are k +12 possible odd values for λ for which υ λ (cid:54) = 0 , and as many even values. Eq. (2) thenleads to the bound P fail k odd ≥ k + 1 , which can only be achieved by ancillary states where alleven values of λ are equiprobable, and all odd values of λ are equiprobable. Similarly, if k is even, there are k possible odd values for λ such that υ λ (cid:54) = 0 , and k + 1 even ones. This leads to the bound P fail k even ≥ k + 2 , which can only be achieved by ancillary states where allvalues of λ for which υ λ (cid:54) = 0 are even and equiprobable.Both limits above can be expressed by the formula P fail ≥ (cid:100) k + 1 (cid:101) even , (4)where (cid:100)·(cid:101) even is the smallest even integer greater or equalto its argument. For the trivial case k = 0 , this result isa special case of the Calsamiglia–Lütkenhaus theorem [3].For k = 1 , we find that a single extra photon does not help,at least with a polarization preserving interferometer.A non-trivial example of state achieving the limit foreven k is the { N +1 − } -photon ancillary state | Υ (cid:105) G · · ·| Υ N (cid:105) G defined by Grice in [1]. Note that, except for the2-photon state | Υ (cid:105) G = | Φ + (cid:105) , the Grice scheme needsto use GHZ-like states [24] of up to N dual-rail qubits.The | Υ n (cid:105) EvL states defined Ewert and van Loock in [2]can be written in the same form of the | Υ n (cid:105) G statesof [1] with respect to the distribution of horizontallypolarized photons; however, at variance with them, inorder to attain the same P succ two copies of each one arerequired. The photon-number dependent upper bound (4)is therefore not tight for them, even if the generic bound(3) is. However, the schemes by Ewert and van Loockstart by independently interfering each of the output of aninitial 50:50 beamsplitter with half of the ancillary state.This additional restriction changes the photon dependentbound, which is then achieved by the Ewert–van Loockschemes. III. LINEAR OPTICAL NETWORKOPTIMIZER
We present here our linear optical network optimizer,a program looking for the optical network maximizing P succ for a given ancillary state. We first restate ourproblem in terms of second quantization in Section III Abefore detailing our approach, which can be divided into two parts. As explained in Section III B, for each ancillawe want to analyze, and for each input Bell state, wegenerate a symbolic expression for the probability ampli-tudes of all output events in terms of U . Those functions,along with their gradient with respect to the U entries,are heuristically optimized in order to reduce the numberof operations needed. The optimized functions are thentranslated into a low-level language and compiled. Then,as exposed in Section III C, a constrained numerical op-timization using a nonlinear method is performed. Dueto the heavy non-smooth character of P succ , we constructa meaningful figure of merit, function of the previouslyobtained probability amplitudes.While this can seem at first glance an overkill brute-force approach, both steps present important symmetriesthat we can exploit, gaining up to two orders of magnitudein computation time in some cases, and reducing functioncomplexity. A. Polynomial representation of the network
We can represent the input (output) state to a n -modeslinear optical network through a polynomial in the input(output) modes creation operators [25] | ψ in (cid:105) = P in ( a † , . . . , a † n ) | (cid:105) , | ψ out (cid:105) = P out ( c † , . . . , c † n ) | (cid:105) . The effect of the network can be represented by a unitarytransformation U := ( u ij ) connecting the input and theoutput modes: a † i = n (cid:88) j =1 u ij c † j . (5)Notice that this only implements a subset of all the trans-formations of the modes allowed by quantum mechanics,namely photon interferences.At variance with the previous section, we representdual rail encoding with distinct spatial modes instead oforthogonal polarization modes [14]. The Bell states arerepresented by: (cid:12)(cid:12) Φ ± (cid:11) = 1 √ (cid:16) a † a † ± a † a † (cid:17) | (cid:105) , (6) (cid:12)(cid:12) Ψ ± (cid:11) = 1 √ (cid:16) a † a † ± a † a † (cid:17) | (cid:105) . (7)Let B β ( a † , . . . , a † ) , β = 1 , . . . , be the polynomial as-sociated with the above states; the generic input to thenetwork is then represented by the polynomial ( P in ) β = B β ( a † , . . . , a † ) Q ( a † , . . . , a † n ) , (8)where Q describes the ancillary state | Υ (cid:105) .At the output of the network, an array of polarizingbeamsplitters followed by photon number resolving de-tectors (PNRD) carries out a measurement in the com-putational basis. A detection event is described by thenumber of photons detected in each output mode, e.g. is the event in which three photons are detected, onein the first mode and two in the third mode. This eventcorresponds to the output state | (cid:105) out = c † c † √ | (cid:105) = c † c † √ | (cid:105) , (9)and its probability can therefore be computed form the thesquared norm of the corresponding monomial coefficientin the polynomial (8).A Bell measurement is unambiguous when at least oneoutput event occurs with nonzero probability for only oneof the four input Bell states (plus ancilla); events meetingthis condition are said to be discriminating . By writing p eβ for the square modulus of the amplitude associatedwith the detection event e when the input is the Bell state β , the total probability of successful discrimination is: P succ = 14 (cid:88) e, ˜ β ∈ S p e ˜ β with S = { e, ˜ β : ∀ β (cid:54) = ˜ β, p eβ = 0 } . (10) B. Symbolic computation
The purpose of the method presented here is to providethe optimum-finding algorithm described in the next sub-section with a fast, optimized function returning all thedetection event probabilities, along with their gradientswith respect to the entries of U , from which a figure ofmerit f ( U ) will be constructed. Working on this sepa-rately, instead of directly using sums of permanents ofsubmatrices of U [26], enables us to carefully analyze theproblem and implement some analytical shortcuts thatwill ultimately speed up the search for optima.We use SymPy [27], an open-source symbolic computa-tion library for Python. The main procedure is as follows:given an input polynomial in the form (8), we implementthe transformation (5) by direct substitution. We then re-group the resulting multivariate polynomial in c † , . . . , c † n and we extract the coefficient of each monomial; theycorrespond to the amplitudes of all possible detectionevents, and they are complex functions of the entries of U . A straightforward combinatorial argument shows thatthere are N = (cid:18) n + k + 1 k + 2 (cid:19) (11)possible detection events for each input state, where k is the number of photons in the ancillary state—making k + 2 the total number of photons entering the network.In the following, we will always assume n ≥ k + 4 .As a simple example of what the algorithm does, con-sider a network with n = 4 modes, with the Bell state | Φ + (cid:105) = B | (cid:105) as input state and no ancillary state ( k = 0 ),with B = 1 √ a † a † + a † a † ) . Let U be the unitary matrix representing the network,the output polynomial is obtained upon performing thesubstitution in eq. (5): P out = 1 √ (cid:16) (cid:88) j u j c † j (cid:17)(cid:16) (cid:88) j u j c † j (cid:17) + 1 √ (cid:16) (cid:88) j u j c † j (cid:17)(cid:16) (cid:88) j u j c † j (cid:17) . (12)After expanding all the products, we obtain a polynomialin c † , . . . , c † with N = 10 terms of degree k + 2 = 2 . Thecoefficient of the monomial c † c † is, for example, the ampli-tude of the detection event 1010. For an arbitrary event,with k i photons in mode i , a bosonic correction factor (cid:112)(cid:81) i k i ! has to be applied, as in eq. (9). Expanding (12),all the output event amplitudes read: −→ u u + u u −→ u u + u u −→ u u + u u −→ u u + u u −→ ( u u + u u + u u + u u ) / √ −→ ( u u + u u + u u + u u ) / √ −→ ( u u + u u + u u + u u ) / √ −→ ( u u + u u + u u + u u ) / √ −→ ( u u + u u + u u + u u ) / √ −→ ( u u + u u + u u + u u ) / √ . The probabilities p eβ are then obtained by taking thesquare moduli of those amplitudes. All the above stepsare automated in our program, by exploiting polynomialmanipulation routines contained in SymPy; we just haveto ‘manually’ feed as input the ancillary polynomial Q .The numerical optimization routine also needs the Ja-cobian of the figure of merit, and the above approachallows us to compute it symbolically in order to speed upthe optimization. This computation is done through thesymbolic derivatives of p eβ with respect to the real andthe imaginary part of u ij . The p eβ are the square modu-lus of a complex holomorphic function (more precisely, apolynomial) and, for any such function g ( u , . . . , u nn ) ,we have ∂ | g | ∂ Re { u ij } = 2 Re (cid:26) g ∂g ∗ ∂u ij (cid:27) and an analogous expression for the derivative with re-spect to the imaginary part. This allows for a simplesymbolic computation of the relevant derivatives. If theoptimization algorithm were not provided with a Jacobianfunction, it would have had to estimate it by evaluatingthe figure of merit at nearby points. Using the finitedifferences method, this amounts to at least n addi-tional evaluations of f ( U ) per iteration—to be comparedwith evaluating n symbolic derivatives, each of which is astrictly simpler function ( i.e. a polynomial with less terms)than the corresponding probability amplitude. Further-more the Jacobian would only be calculated to some fixedaccuracy, which adds noise to the optimization algorithmand worsen its convergence.As the output of this computation, we need N expres-sions for the probabilities, each one with n expressions forthe derivatives. But there is no need to actually computethe symbolic form of all the N events for each ancilla.Two main simplifications help to drastically reduce thenumber of symbolic computations:1. The detection events divide into equivalence classesunder permutation of the output modes, correspond-ing to a permutation of the columns of U . It is easyto check that, of all the events in the above exam-ple, it is only necessary to obtain 2000 and 1100.For example, 1010 can be obtained from 1100 byswapping the second and the third column of U ,at a negligible computational cost. The number ofindependent events does not depend on n and isequal to P k +2 , the number of integer partitions [28]of the total number of photons. This also reducesthe number of gradients to calculate, from n to atmost n ( k + 2) .2. As it is clear from eqs. (6) and (7) , knowing anevent’s amplitudes for a single Bell state allows oneto compute its amplitudes for all of them, just byswapping two rows of U and/or changing their sign.This reduces the number of functions to computeby a further factor of .With these expedients, we reduce the problem to thecomputation of just P k +2 probability functions along withtheir gradients. For example, when n = 8 , k = 2 as inthe first scheme of Grice’s paper (using a | Φ + (cid:105) ancillarystate), we end up decreasing the number of function tocompute from N = 1320 to only P = 5 , each with atmost 32 derivatives instead of 64.As SymPy is written in pure Python, we can acceleratethis section of the program using PyPy [29], an imple-mentation of the Python interpreter with a Just-In-Timecompiler. For large networks ( n ≥ , k ≥ ) we obtaina tenfold speedup over plain Python, at the cost of in-creased memory usage . Even with all these optimizationsin place, however, the problem still scales exponentially .For comparison, we obtain the functions for the first Griceiteration (“one extra Bell pair” in Table II) in about 6 I.e. the number of positive integer sums equal to k + 2 . As an example, first iteration of the Ewert–van Loock scheme( | (cid:105) ⊗ = | Υ (cid:105) ⊗ EvL ) takes 7 minutes and 30 seconds on our laptop(see Table I) using the standard Python interpreter and about150 MB of RAM. Using PyPy the time is cut down to 45 seconds,with a memory consumption of 250 MB. We nonetheless get an exponential advantage over the naiveapproach, as we can see from the asymptotic formula for P k due TABLE I. Specifications of the two computers used in this ar-ticle. The frequency is the nominal frequency of the processor.The cluster is the gmpcs-206 branch of the computing centerMésoLUM of the LUMAT research federation [32], and thespecifications refer to a single node.Name Processor seconds using 100 MB of RAM; the second iteration ofGrice’s scheme—the largest calculation we managed tocomplete, with n = 16 and k = 6 —took instead 10 days ofsingle-core CPU time on the cluster described in Table I,requiring a large portion of the 256 GB of RAM at ourdisposal. The resulting (already heavily optimized) prob-ability functions for this case consist of a grand total ofabout . million addition and multiplication operations,and the Jacobian of about million.We use Theano [31], a numerical computation libraryfor Python, as our backend in order to translate eachfunction into C and compile it to a fast numerical version. C. Numerical optimization
The second part of the procedure takes care of findingthe optimal value of P succ in eq. (10) for each input ancil-lary state. Naively, we could straightforwardly calculate P succ ( U ) from the numerical evaluation of the probabil-ity functions obtained in the previous section; however, P succ ( U ) = 0 almost everywhere in the domain U ( n ) , andit is neither continuous nor differentiable in the region ofinterest, where P succ ( U ) (cid:54) = 0 . This obviously makes mostoptimization methods highly ineffective.As a workaround, we devise a figure of merit as acontinuous alternative to P succ ( U ) . We thus search forlocal minima of: f ( U ) = (cid:88) e (cid:16) (cid:88) β p eβ − α p eα (cid:17) , (13)of which the addends of the outer sum are equal to theones of −P succ ( U ) , when the latter happen to be nonzero.Ideally, we want the figure of merit to closely mimic the to Hardy and Ramanujan [30]: P k ∼ k →∞ k √ (cid:40) π (cid:114) k (cid:41) , to be compared with the binomial coefficient in eq. (11), lowerbounded by k +2 . behavior of P succ around the optima. In particular itwould be useful to have the following holding for eachpair of locally optimal unitaries U and U : f ( U ) < f ( U ) iff P succ ( U ) > P succ ( U ) . Unfortunately, the mutual relation between f ( U ) and P succ ( U ) is not simple. While it seems reasonable toconjecture the set of local minima of f ( U ) to include theset of P succ ( U ) ’s maxima, we find that, for some U and U : f ( U ) < f ( U ) while P succ ( U ) < P succ ( U ) . (14)This forces us to conduct a search among all local op-tima of f ( U ) , rather than taking advantage of globaloptimization methods like simulated annealing.The optimization itself is implemented using SciPy’s Sequential Least-Square Programming (SLSQP) optimiza-tion method [33, 34]. As this method supports equalityconstraints, we represent U through n real variablesdescribing its entries’ real and imaginary parts, alongwith n equations enforcing orthonormality of the rows.In order to improve convergence speed, the method isallowed some leniency on the constraints, in that theyonly have to be satisfied at the local optimum. In orderto ensure uniform sampling of the starting points of eachoptimization, we choose them by picking a random ma-trix from the Haar measure on the unitary group U ( n ) .This is accomplished using the QR decomposition methoddescribed in [35].We compared the performance of the constrainedmethod above to the Broyden-Fletcher-Goldfarb-Shanno(BFGS) algorithm [36], a popular unconstrained quasi-Newton method. While the latter uses no constraintsand has to work with independent variables , the relationbetween those and the entries of U is non-trivial andthis therefore hinders our ability to input the analyticalform of the gradient. Thus the performances are com-parable with SLSQP, if not worse in some cases, even ifthe number of variables is cut by half. Furthermore, theconvergence accuracy and precision seems unaffected bythe choice of one method over the other. IV. RESULTS
Below, we work out the analytical upper bound forpolarization-preserving interferometers in Section II forthe states we used, and we compare them to the numericalresults we obtained for generic interferometers.Some of the optimization results for different input an-cillæ | Υ (cid:105) = Q ( a † , . . . , a † n ) | (cid:105) are summarized in Table II.For each one of them, we collected the local optima from We encoded the n unconstrained degrees of freedom of U ( n ) into an Hermitian matrix H , using U ( n ) = e iH . about ten thousand successful iterations of the optimiza-tion algorithm. In the table the maximum value achievedfor each input is shown; it can be noted that, for the casesalready known in the literature, we find the same maximaldiscrimination probability. Furthermore, we were able towork with other types of ancillæ. A. Vacuum extra modes
By virtue of the Calsamiglia–Lütkenhaus theorem [3],the analytical upper bound of P succ ≤ / is known tohold for general interferometers, equipped with an ancillaconsisting of an unlimited number of extra modes in thevacuum state. We can work out different version of thisbound (for photon-number- and polarization-preservingmeasurements) following the reasoning laid out in Sec-tion II, looking at the distinguishability of the states ineq. (1) after a projection onto the basis of the number-of-horizontally-polarized-photons operator. If the ancillarystate is the vacuum, or any other state with a fixed num-ber λ of horizontally polarized photons, υ λ is the onlyvalue of λ for which υ λ (cid:54) = 0 . This leads to just two distinctterms: | HH (cid:105) | Υ (cid:105) for both (cid:12)(cid:12) Φ + (cid:11) and (cid:12)(cid:12) Φ − (cid:11) , (15) ± | V V (cid:105) | Υ (cid:105) respectively for (cid:12)(cid:12) Φ + (cid:11) , (cid:12)(cid:12) Φ − (cid:11) . (16)The term (15) is identical for both inputs, while the (16)differs by a global phase. Thus, these terms are notdistinguishable at all: in this case the ancillary statecannot help to discriminate between the two differentinputs, and P succ ≤ / .Indeed, without extra photons the maximum discrim-ination probability that we find through our numericaloptimization is / , for any value of n we tried. Wequickly achieve this maximum on our laptop (see Table I),and we collect a thousand successful iterations in a matterof minutes for different values of n ≤ . The n = 14 case still takes less than an hour on the laptop, and a fewminutes of the cluster. B. Extra Bell pairs
One of the simplest linear optical network schemes withancilla achieving more-than- Bell state discriminationprobability is arguably the first iteration of Grice’s strat-egy [1]. He shows that adding one extra | Φ + (cid:105) as ancillahelps cutting the degeneracy of | Φ ± (cid:105) by half, achieving P succ = 3 / . However, the states used by Grice to increaseits success probability past / become more complex ateach iteration, since each additional ancillary state | Υ N (cid:105) G is a N -photon GHZ state. It would be experimentallymuch simpler to use multiple Bell pairs | Φ + (cid:105) ⊗ k/ as a k -photon ancilla, which motivate our research of schemesusing such resources. TABLE II. Summary of known analytical and numerical results for different ancillæ. Q is the input polynomial, n the numberof modes and k the number of photons. P numsucc is the optimum obtained through our optical network optimizer; the fractiongiven is exact up to our numerical precision (9 decimals). P anasucc is the best known explicit analytical result. P uppsucc and P uppsucc ( k ) are our analytical upper bounds for polarization-preserving networks (section II), for the ancilla and for arbitrary ancillæ withsame k . They are in bold font when matching the best known result. State
Q n k P numsucc P anasucc P uppsucc P uppsucc ( k ) Vacuum Ancilla | (cid:105) a [3] k/ Extra Bell Pairs (cid:12)(cid:12) Φ + (cid:11) ⊗ k/ a † a † + a † a † ) ... ( a † k +1 a † k +3 + a † k +2 a † k +4 )2 k/ k + 4 even — b − ( k/ (cid:98) k/ (cid:99) ) k/ c (cid:39) − √ πk k + 1 k + 2 (cid:12)(cid:12) Φ + (cid:11) = | Υ (cid:105) G ( a † a † + a † a † ) √ (cid:12)(cid:12) Φ + (cid:11) ⊗ a † a † + a † a † )( a † a † + a † a † )2
12 4 3/4 d (cid:12)(cid:12) Φ + (cid:11) ⊗ a † a † + a † a † ) ... ( a † a † + a † a † ) √
16 6 e d k Extra Photons | (cid:105) ⊗ k a † . . . a † k +4 k + 4 even — b − ( k/ (cid:98) k/ (cid:99) ) k/ c (cid:39) − √ πk k + 1 k + 2 | (cid:105) ⊗ k a † . . . a † k +4 k + 4 odd — b same as above, for k − | (cid:105) a † d d | (cid:105) ⊗ a † a † c ( ) f | (cid:105) ⊗ a † a † a † d c ( ) f | (cid:105) ⊗ = | Υ (cid:105) ⊗ a † a † a † a † c | (cid:105) ⊗ a † . . . a †
10 6 3/4 d c | (cid:105) ⊗ a † . . . a †
12 8 e c (25/32) f | (cid:105) ⊗ a † . . . a †
16 12 e c (13/16) f Grice Schemes [1] (first iteration is (cid:12)(cid:12) Φ + (cid:11) = | Υ (cid:105) G above) | Υ (cid:105) G · · · | Υ N (cid:105) G Straightforward, but long expression k + 4 2 N +1 − — k + 1 k + 2 k + 1 k + 2 k + 1 k + 2 | Υ (cid:105) G | Υ (cid:105) G ( a † a † + a † a † )( a † a † a † a † + a † a † a † a † )2
16 6 9/16 gh (first iteration is | (cid:105) ⊗ = | Υ (cid:105) ⊗ above) ( | Υ (cid:105) EvL · · ·| Υ N (cid:105) EvL ) ⊗ Straightforward, but long expression k + 4 2 N +2 − — k + 2 k + 4 k + 2 k + 4 k + 1 k + 2 (cid:18) k + 2 k + 4 (cid:19) f GHZ states | GHZ k (cid:105) a † ··· a † k +3 + a † ··· a † k +4 √ k + 4 k — 3/4 i c − (cid:100) k +1 (cid:101) even | GHZ (cid:105) a † a † a † + a † a † a † √
10 3 3/4 3/4 i c | GHZ (cid:105) = | Υ (cid:105) G a † a † a † a † + a † a † a † a † √
12 4 3/4 3/4 i c W State | W (cid:105) a † a † a † + a † a † a † + a † a † a † √ h . h c (3/4) j a Also holds for polarization non-preserving interferometers. b No generic scheme is known. c Polarization-preserving bound obtained after rotating the polarization of some or all modes by π . d The best known interferometer correspond to a smaller ancilla, together with ignoring extra modes. e Computation out of reach for our program. f For networks which start by interfering the two bell states on a 50:50 beamsplitter, analyzing each half separately. g Computation at the borderline of our computing capacity: this result is the best of just 12 optimizations over the course of three weeks. h Numerical result worse than the best known analytical scheme. i Success probability achieved by measuring all photons of the ancilla and using the remaining in a “one extra Bell pair” scheme. j Obtained through a more complex transformation of the input, exposed in the main text.
We start by working out the polarization-preservingbound of Section II. We restrict ourselves to a productof k/ states of the form | Υ (cid:105) = √ ( | H (cid:105) + | V (cid:105) ) , where | H (cid:105) (resp. | V (cid:105) ) is any state of two horizontally (resp.vertically) polarized photons. Bell pairs are a special caseof the latter, as well as the | Υ (cid:105) EvL defined in [2]. Wehave | Υ (cid:105) ⊗ k/ = 2 − k/ ( | H (cid:105) + | V (cid:105) ) ⊗ k/ = 2 − k/ k/ (cid:88) λ =0 (cid:115)(cid:18) k/ λ (cid:19) | Υ , λ (cid:105) , where | Υ , λ (cid:105) is the uniform superposition of the termswith λ horizontally polarized photons present in theexpansion of | Υ (cid:105) ⊗ k/ . Equation (2) then only have evennonzero terms, and this leads to P fail | Φ + (cid:105) ⊗ k/ ≥ − k/ − (cid:18) k/ (cid:98) k/ (cid:99) (cid:19) (17) ≥ − k/ − √ πk k/ e − / k = 1 √ πk e − / k , (18)where we have supposed k to be a multiple of and applieda second-order version of Stirling’s approximation [37] √ πn (cid:16) n e (cid:17) n ≤ n ! ≤ √ πn (cid:16) n e (cid:17) n e / n When k is even, but not a multiple of , the inequality (18)is invalid, but we still have P fail | Φ + (cid:105) ⊗ k/ ≥ √ πk (cid:18) O (cid:18) k (cid:19)(cid:19) . (19)The / √ k scaling of P fail allowed by the above bound isworse than the /k scaling achieved by Grice schemes.Nevertheless, it does not rule out strategies approachingsuccess probabilities arbitrarily close to 1 by using asinputs much simpler states, i.e. k/ Bell pairs.Our numerical search, which is not restricted topolarization-preserving schemes, converges in just abouta minute to the P succ = 3 / scheme on our laptop, using300 MB of RAM. Unfortunately, the use of two extra Bellpairs shows no improvement over a single extra Bell pair,and its optimization uses significantly more resources:about 4 hours with 20 parallel threads on the cluster(see Table I), each using 3 GB of RAM, for the collec-tion of a thousand optimizations. For three extra-bellpairs, the polarization-preserving bound of eq. (17) gives P succ ≤ /
16 = . for the latter, allowing in principlefor a scheme beyond / . However, this dimensionality isbarely out of reach for our program, even using the clus-ter. For comparison with a similar-sized case, the seconditeration of Grice’s strategy (the symbolic function’s sheersize of which we discussed at the end of Section III B) takes about 48 hours on the cluster for each starting pointto converge to a local optimum. We collected just 12optimizations, obtaining P succ = 9 / ; unfortunately thisresult is well below the known Grice’s 7/8 scheme. C. Extra single photons
The possibility of improving the discrimination proba-bility through the use of unentangled extra single photonsis of great experimental interest, especially with the recentdevelopment of high-efficiency single photon sources withnear ideal indistinguishability [38]: such ancillary stateswould be among the simplest types of input states for areal-world implementation of linear optical Bell measure-ments.Ewert and van Loock explore the use of pairs of sin-gle photon per auxiliary dual-rail mode [2, Section D ofsupp. mat.] as substitutes of their ancillary state | Υ (cid:105) EvL .While the initial transformation they apply to the inputphotons is polarization dependent, we can still use theformalism of our polarization-preserving upper bound ofSection II, restricted to the case in which each photonenters the network polarized along the ± π direction. Inthe horizontal-vertical basis, this ancilla is described bythe Hong-Ou-Mandel state [39] | Υ (cid:105) EvL = √ ( | (cid:105) + | (cid:105) ) in each mode pair. The υ λ coefficients in the case of the k -photon state (with k even) | Υ (cid:105) ⊗ k/ EvL correspond to theones of k/ Bell states | Φ + (cid:105) ⊗ k/ . We can therefore applythe same reasoning laid out in Section IV B, obtainingthe bound in eq. (17).With this restriction in place, we get for k ≥ a slightlytighter lower bound to P fail , compared to the photon-number based bound (4). For example, with 4 singlephotons the latter gives P fail ≥ / , while eq. (17) gives P fail ≥ − (cid:0) (cid:1) = 1 / . In fact, the bound is saturated bythe 4-single-photon variant of the first iteration of Ewert-van Loock strategy. They also consider the 12-single-photon state | (cid:105) ⊗ → | Υ (cid:105) ⊗ . In this case, a direct ap-plication of equation (17) leads to P fail ≥ − (cid:0) (cid:1) = 5 / ,which is indeed smaller than the actual / failure ratefound by the authors. However a look into the detailedsymmetry of the strategy, as per the same reasoning usedat the end of Section II B, leads to the better (but morerestricted) bound P fail ≥ − (cid:0) (cid:1) = 3 / , which is closer,but still below / .The aforementioned 4-photon scheme discriminates | Ψ + (cid:105) and | Ψ − (cid:105) with certainty, and | Φ ± (cid:105) only half ofthe times. Our numerical algorithm indeed finds this (1 , , , ) scheme when initialized with a 4 single photonancilla, and does not manage to improve its probabil-ity of success—an evidence of its optimality even in thepolarization-dependent case. Furthermore, we find otherschemes achieving the same total success probability witha different discrimination pattern among the Bell states: (1 , , , ) .Interestingly, with just two extra photons, we find two FIG. 1. The first “half” Ewert–van Loock single-photonsscheme. It performs a Bell measurement with P succ = 5 / on the state | β (cid:105) using two unentangled extra photons | (cid:105) ,two polarization-independent beamsplitters, four polarizingbeamsplitters, two phase shifters and six photocounters. λ / λ / | (cid:105)| (cid:105) | β (cid:105) schemes, (1 , , , ) and (1 , , , ) , achieving a discrimi-nation probability of P succ = 5 / . . The first canbe easily described as half of the 4-photon Ewert–vanLoock scheme [2] mentioned above and is described inFigure 1. It is especially relevant experimentally, sinceit is the simplest scheme achieving a success rate above / . It was independently obtained by Ewert and vanLoock [40]. By “halving” in the same way the Ewert–vanLoock 12-photon scheme that uses single photonsand achieves a probability of / , we find a similar“intermediate” scheme with 4+4=8 extra photons, achiev-ing P succ = 49 / . Unfortunately, even using the cluster,numerically testing this scheme ( n = 12 , k = 8 ) proved tobe unfeasible.We notice (as in [19]) that using an odd number k + 1 of single photons in the ancilla does not improve thediscrimination probability over the case with k photons.This is in line with the analytically-derived behavior forpolarization-preserving interferometers of Section II B.Very recently, Smith and Kaplan [19] tackled a simi-lar problem, numerically optimizing linear optical Bellmeasurements with single photons ancillæ. Their mea-surement were allowed to be ambiguous, and the chosenfigure of merit was the classical mutual information be-tween state preparation and measurement. Remarkably,despite this difference, we find corresponding results forancillæ up to five single-photons; with six photons, theyfind a slight improvement of their mutual information,but the corresponding measurement is ambiguous [41].Even if with six photons (Table II) we could not findany scheme beyond P succ = 3 / —we collected more than10 000 optimizations—the polarization-preserving boundallows for a scheme with P succ ≤ / ; an improvementover / is therefore not excluded. D. GHZ and W states
We also checked the possible use of multipartite en-tangled states and ancillæ. A three-photon GHZ stateancilla, | GHZ (cid:105) = 1 √ (cid:16) | (cid:105) + | (cid:105) (cid:17) , does not seem to help with respect to a simple Bell pair,as we still attain / discrimination probability as op-timum. So does a GHZ state, at the expense of morecomputational power; we wrongly expected the latter tobe useful, given its use (along with a Bell pair) in thesecond iteration of Grice’s scheme [1]. The analyticalpolarization-preserving bound predicts P succ ≤ / forall | GHZ k (cid:105) when k ≥ . However, for odd k , the rotationof the polarization of a single photon by an angle of ± π raises this bound to 3/4. This value can be achieved bya trivial network applying a simple π rotation on k − spatial modes of the ancilla, which leaves the remainingtwo photons in the | Φ ± (cid:105) states, which can be used asdescribed above to achieve a 3/4 success probability.Another interesting state to investigate is the three-photon W state, | W (cid:105) = 1 √ (cid:16) | (cid:105) + | (cid:105) + | (cid:105) (cid:17) . (20)Like GHZ , it is a genuinely 3-party entangled state but,unlike all other states studied above, it is not a graphstate, not even a stabilizer state. Its specific symmetryis likely the source of the interesting results we find (endof Table II). Having the same number of horizontallypolarized photons in each term, this state is as usefulas the vacuum for polarization-preserving interferome-ters, as showed in Section IV A. However, the rotationof the polarization of two photons by π gives the higherbound P succ ≤ / , and further manipulation (see below)raises this bound to P succ ≤ / . The best optimum wefind numerically, when we use a network with no extravacuum modes ( n = 10 ), is P succ = 5 / , significantlylower than the / achieved with a simpler two-photonBell pair. This optimum is extremely rare (once in morethan 20 000 optimizations), and we observe the figureof merit in this case to suffer heavily from the problemdescribed in eq. (14) about the relationship between f ( U ) and P succ ( U ) .However, in this case we could find a better scheme bymanipulating the state “by hand”. By measuring the lasttwo spatial modes we can apply a transformation such thatthe remaining modes can be, depending on the result ofthe measurement, either in the state √ ( | H, (cid:105) − | , V (cid:105) ) The phase of the Bell pair is determined by the parity of the mea-surement of the k − photons, and its effect is simply exchangesthe photon patterns for the detection of (cid:12)(cid:12) Φ + (cid:11) and (cid:12)(cid:12) Φ − (cid:11) . | Φ + (cid:105) . Applying to these modes the same unitary of theone-Bell-pair P succ = 3 / Grice strategy gives a scheme for | W (cid:105) with P succ = 7 / . While our optimization programcorrectly identifies this scheme as a local optimum whenput in as starting point, an added Gaussian noise ofaverage magnitude well below the requested convergenceaccuracy is sufficient for the optimization to diverge fromit. This numerical fragility may be the reason why wecould not find this optimum through the optimization.Applying the analytical bound to such transformed ancillagives us P succ ≤ / . Interestingly, adding at least twovacuum modes ( n ≥ ) allows the program to reach thebetter discrimination probability of 0.5785508(2). Still,this is slightly below the manually-found 7/12. V. CONCLUSION
In this work we have investigated the optimal successprobability of a linear optical Bell measurement assistedby different kinds of input ancillary states | Υ (cid:105) . In Sec-tion II, we showed how to obtain an upper bound fromthe input photon polarization distribution in | Υ (cid:105) , whenthe network is restricted to polarization-preserving inter-ferometers; we noticed that the bound is tight for somepublished schemes. With the aim of exploring the pa-rameter space of generic interferometers, we developed inSection III a linear optical network simulator, capable ofevolving a generic input state through the network andcomputing the analytical expression of the probabilities ofeach detection event in the output. We then conducted anumerical search for the optimal value of P succ in for fixed | Υ (cid:105) , and we discussed how to reduce the overall computa-tional cost by exploiting some symmetries of the problemat hand. We presented the results of both analytical andestimated numerical bounds in Section IV, and we recallthem in Table II.Through both the analytical study and the numericaloptimization we find evidences (but no proofs) for theoptimality of known small schemes. Some of them seemachievable experimentally in the short term, as they re-quire as ancilla either a small number of photons or anadditional Bell pair. While restricted to the polarization-preserving case, the photon-number based analytical up-per bound, saturated by Grice’s schemes, is evidence fortheir optimality if resources are measured in terms of thenumber of extra ancillary photons. In this setting, we have also shown that employing many copies of a Bell pairleads to a different (and worse) scaling than using Grice’sstates, giving interesting insights into the reason why thebig GHZ-like states that appear in the schemes of [1, 2]are needed. Of course, eq. (17) being only a bound, moreresearch is needed to investigate its tightness, and whethernear unity success can indeed be achieved.As pointed out in the paper, some interesting caseslie beyond the computational capabilities at our disposal.While there is still room for improvement, e.g. by furtheroptimization of the code and/or by employing more CPUtime, our numerical approach is at least as hard as comput-ing permanents of k × k submatrices of a unitary matrix.As proved by Valiant [42] and more recently pointed outby Arkhipov and Aaronson [43] in the context of linearoptics, this task pertains to the complexity class P -hardand is not believed to be solvable in polynomial time ona classical computer. However the symmetry of the Bellstates and the unambiguity constraints, which enforcea structure on the matrix entries—by imposing manynull probabilities—may enable significant speedups (evenexponential ones), even if the overall scaling could stayexponential. Recent works in [44, Appendix B] and [45,Appendix D] suggest optimized algorithms for computingthe permanent of matrices with repeated columns/rows;they may help to improve our computation.We conclude by noting that our simulator might beuseful in exploring the power of linear optics in solvingother types of problems. Due to the flexibility of Pythonand of the separation between symbolic computation andnumerical optimization, the program only requires minormodifications in order to be adapted to new tasks. As amatter of fact, it has already been used during discussionswith Chabaud et al. in order to gain insight on the effectof Hadamard networks, helping in the design of linearoptical swap-test [46]. ACKNOWLEDGMENTS
We warmly thank the Quantum Information team of theLIP6 for their hospitality, and especially Ulysse Chabaudfor stimulating discussions. We acknowledge the use ofthe computing center MésoLUM of the LUMAT researchfederation (FR LUMAT 2764). AO acknowledges finan-cial support from ANR project ANR-16-CE39-0001 andthe Erasmus+ Traineeship Programme of the EuropeanUnion. [1] W. P. Grice, Physical Review A , 042331 (2011).[2] F. Ewert and P. van Loock, Physical Review Letters ,140403 (2014).[3] J. Calsamiglia and N. Lütkenhaus, Applied Physics B ,67 (2001).[4] S. L. Braunstein and A. Mann, Physical Review A ,R1727 (1995). [5] C. H. Bennett, G. Brassard, C. Crépeau, R. Jozsa,A. Peres, and W. K. Wootters, Physical Review Let-ters , 1895 (1993).[6] S.-W. Lee and H. Jeong, in proceedings of the First Inter-national Conference on Entangled Coherent State and ItsApplication to Quantum Information Science (TamagawaUniversity, Tokyo, 2012) pp. 41–46. [7] D. Bouwmeester, J.-W. Pan, K. Mattle, M. Eibl, H. We-infurter, and A. Zeilinger, Nature , 575 (1997).[8] C. H. Bennett and S. J. Wiesner, Physical Review Letters , 2881 (1992).[9] K. Mattle, H. Weinfurter, P. G. Kwiat, and A. Zeilinger,Physical Review Letters , 4656 (1996).[10] M. Żukowski, A. Zeilinger, M. A. Horne, and A. K. Ekert,Physical Review Letters , 4287 (1993).[11] J.-W. Pan, D. Bouwmeester, H. Weinfurter, andA. Zeilinger, Physical Review Letters , 3891 (1998).[12] N. Lütkenhaus, J. Calsamiglia, and K.-A. Suominen,Physical Review A , 3295 (1999).[13] E. Knill, R. Laflamme, and G. J. Milburn, Nature ,46 (2001).[14] P. Kok, W. J. Munro, K. Nemoto, T. C. Ralph, J. P.Dowling, and G. J. Milburn, Reviews of Modern Physics , 135 (2007).[15] H. A. Zaidi and P. van Loock, Physical Review Letters , 260501 (2013).[16] D. Vitali, M. Fortunato, and P. Tombesi, Physical ReviewLetters , 445 (2000).[17] S.-W. Lee and H. Jeong, Physical Review A , 022326(2013).[18] S.-W. Lee, K. Park, T. C. Ralph, and H. Jeong, PhysicalReview Letters , 113603 (2015).[19] J. A. Smith and L. Kaplan, arXiv:1802.10527 [quant-ph](2018).[20] A. Olivo and F. Grosshans, arxiv.org/src/1806.01243/anc(supplementary information of this paper, 2018).[21] I. Ivanovic, Physics Letters A , 257 (1987).[22] D. Dieks, Physics Letters A , 303 (1988).[23] A. Peres, Physics Letters A , 19 (1988).[24] D. M. Greenberger, M. A. Horne, and A. Zeilinger, in Bell’s Theorem, Quantum Theory and Conceptions ofthe Universe (Springer Netherlands, Dordrecht, 1989) pp.69–72.[25] V. Fock, Zeitschrift für Physik , 622 (1932); englishtranslation in L. D. Faddeev, L. A. Khalfin, and I. V.Komarov, V.A. Fock–selected works : quantum mechanicsand quantum field theory (Chapman & Hall/CRC, 2004)p. 567.[26] S. Scheel, in
Quantum Information Processing (Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim, FRG, 2005)pp. 382–392. [27] A. Meurer, C. P. Smith, M. Paprocki, O. Čertík, S. B.Kirpichev, M. Rocklin, A. Kumar, S. Ivanov, J. K. Moore,S. Singh, T. Rathnayake, S. Vig, B. E. Granger, R. P.Muller, F. Bonazzi, H. Gupta, S. Vats, F. Johansson,F. Pedregosa, M. J. Curry, A. R. Terrel, S. Roučka, A. Sa-boo, I. Fernando, S. Kulal, R. Cimrman, and A. Scopatz,PeerJ Computer Science , 27 (2017).[28] M. Abramowitz, I. A. Stegun, and R. H. Romer, AmericanJournal of Physics , 958 (1988).[29] S. Pedroni, “Goals and Architecture Overview — PyPydocumentation,” .[30] G. H. Hardy and S. Ramanujan, Proceedings of the Lon-don Mathematical Society s2-17 A software package for sequential quadraticprogramming , Forschungsbericht. Deutsche Forschungs-und Versuchsanstalt für Luft- und Raumfahrt, DFVLR(DFVLR, Braunschweig, 1988).[35] M. Ozols, “How to generate a random unitary matrix,”(2009).[36] J. Nocedal and S. J. Wright,
Numerical optimization (Springer, 2006) p. 664.[37] H. Robbins, The American Mathematical Monthly , 26(1955).[38] P. Senellart, G. Solomon, and A. White, Nature Nan-otechnology , 1026 (2017).[39] C. K. Hong, Z. Y. Ou, and L. Mandel, Physical ReviewLetters , 2044 (1987).[40] P. van Loock, in presentation at “Secure Communicationvia Quantum Channels” (Bielefeld, Germany, 2017).[41] J. A. Smith, private communication (2017).[42] L. Valiant, Theoretical Computer Science , 189 (1979).[43] S. Aaronson and A. Arkhipov, Theory of Computing ,143 (2013).[44] M. C. Tichy, Entanglement and interference of identicalparticles , Ph.D. thesis, Freiburg University (2011).[45] V. S. Shchesnovich, International Journal of QuantumInformation11