## Abstract

Boson sampling can provide strong evidence that the computational power of a quantum computer outperforms a classical one via currently feasible linear optics experiments. However, how to identify an actual boson sampling device against any classical computing imposters is an ambiguous problem due to the computational complexity class in which boson sampling lies. The certification protocol based on bosonic bunching fails to rule out the so-called mean-field sampling. We propose a certification scheme to distinguish the boson sampling from the mean-field sampling for any random scattering matrices chosen via the Harr measure. We numerically analyze our scheme and the influence of the imperfect input states caused by non-simultaneous arrival photons.

© 2016 Optical Society of America

## 1. INTRODUCTION

A major goal in quantum information is the implementation of a genuine and universal quantum computer that can solve computational tasks beyond the power of a classical computer. This goal, however, faces huge challenges because current technology cannot maintain coherence long enough to perform a full computation potentially involving a large number of qubits. Several quantum algorithms, such as Shor’s factoring algorithm [1,2], can compute certain problems for which there are no known classical algorithms with polynomial time complexity. However, the near-term prospects for Shor’s factoring algorithm to overturn the central role of the extended Church–Turing thesis [3] in classical computer science are low due to the difficultly of its experimental realization [4].

Boson sampling provides a different path toward showing the increased power of quantum computation [5]. The complexity of boson sampling is believed to exceed that of Shor’s factoring algorithm. If a classical algorithm existed to sample the BosonSampling distribution, even in an approximate sense, it would lead the polynomial hierarchy to collapse to the third level, something that is believed not to occur. Experimental boson sampling with up to six photons has been reported, demonstrating the theoretical prediction of the output distribution [6–10]. Temporal encoding has been used to extend the number of scattering modes to ten, with the potential to be extended much further [11,12]. In addition, research on sampling from closely rated distributions has extended our knowledge about the nature of boson sampling [13–15].

While the hardness of boson sampling is its key feature, the hardness of boson sampling requires that the certification of its output is not easy. This is because boson sampling is in a class of problems that is larger than the class NP (which Shor’s algorithm belongs to). This leaves room for a classically efficient sampling algorithm to fake the boson sampling distribution in a way that is not efficiently detectable. This has motivated the recent study of different strategies to certify the boson sampling distribution against alternative distributions that attempt to deceive a verifier. The validation test of Aaronson and Arkhipov discriminates the boson sampling distribution against the uniform one within a constant number of samples [16–18]. Another protocol based on the tendency for bosons to bunch together in interferometers rules out the classical sampling distribution, which uses distinguishable photons as input to the sampler [19,20]. A so-called mean-field sampler [21] invalidates these strategies by replicating the characteristic bosonic bunching behavior, and this type of sampling has been experimentally performed [22]. Although this certification scheme was proposed to rule out plausible physical models by observing if some forbidden events occur in the output modes, its effective range is strictly limited to the symmetric situation in which we need to artificially construct a boson sampler with a cyclic-symmetry initial state and an interaction described by a Fourier matrix [21].

In this paper, we describe a certification scheme that can rule out mean-field sampling devices from the true boson sampling without changing the internal working of the device and uses the same distribution of the input transformation matrices as required for true boson sampling. The model we use for a potential boson sampling device is as a black box whose detailed internal functioning is unknown. The black box takes as an input a Harr-random unitary transformation matrix and an input Fock state described by a sequence of zeros and ones. The device outputs Fock basis measurement results, which are sequences of integers. In our certification protocol, we inject a test input state consisting of single photons into every mode. The mean-field sampling devices produce uniform distributions for this test state, and true boson sampling does not. It is then possible to distinguish a boson sampler from a mean-field sampler via analyzing the variations in the output probability distribution. We also analyze the noise characteristics for partially distinguishable particles in our certification scheme.

## 2. SAMPLING MODELS

#### A. Boson Sampling

An ideal boson sampler scatters $n$ indistinguishable photons by an $N$-port interferometer of lossless, linear optical elements, i.e., phase shifters and beam splitters, and counts the number of photons exiting from each output mode. The input state is a Fock state with definite numbers of photons occupying each mode, which is expressed by

The dimensionality of a unitary transformation $\widehat{\mathbb{S}}$ between the input and output state vectors will increase exponentially as $n$ and $N$ grow. However, when restricting the transformations to those of linear optical elements, the transformation can be represented by an $N\times N$ complex unitary matrix $\mathrm{\Lambda}$, which describes the scattering of photons. One can write this transformation down in terms of the Heisenberg evolution of the creation operators as

The hardness of boson sampling is shown by using the boson sampling process to make approximations of the permanent of a matrix. The proof depends on some plausible conjectures about approximating the permanents of i.i.d. Gaussian matrices [5]. To be close to the Gaussian matrices, it is desirable to choose $\mathrm{\Lambda}$ randomly according to the Harr measure on $N\times N$ unitary matrices and for $N$ to scale at least as ${n}^{2}$. Physically, it requires that the scattering process is arbitrarily tunable to perform the corresponding transformation $\mathrm{\Lambda}$. With lossless optical elements and efficient photon detectors, a physically constructed boson sampler would accurately generate samples from this probability distribution. The existence of a classical algorithm that simulates such a quantum device in polynomial time would imply a collapse in the polynomial hierarchy [5]. The collapse of the polynomial hierarchy has quite a technical definition, but essentially, a collapse of the hierarchy is akin to a statements of the form $P=NP$, but depending on exactly which level the hierarchy collapses, the weaker the statement is about the relationship between $P$ and $NP$. Nevertheless, it is believed that the hierarchy does not collapse at any level, and therefore, the boson sampling problem is hard for any classical algorithm.

#### B. Classical Sampling

One way to construct a classical analog of boson sampling is to take as an input distinguishable photons instead of indistinguishable photons. Given an input state ${S}_{D}=({s}_{1},{s}_{2},\dots ,{s}_{N})$, the probability to detect the output configuration ${T}_{D}=({t}_{1},{t}_{2},\dots ,{t}_{N})$ is

#### C. Mean-Field Sampling

The mean-field sampler devised by Tichy *et al.* [21] is an efficiently computable and physically plausible model of a boson sampler imposter. The interference of Bose–Einstein condensates motivates the physical implementation of the mean-field sampler, which provides an efficient approximation for the probability distribution of the boson sampling. The input state ${S}_{Mf}$ is a set of $n$ identical single-photon states

On the side of output ports, each photon would be detected on the $k$th mode with the probability

Injecting $n$ single-photon states as per Eq. (5) into the scattering setup shot by shot, we can detect the output event ${T}_{Mf}=({t}_{1},{t}_{2},\dots ,{t}_{N})$, whose mode arrangement is $\mathrm{k\u20d7}=({k}_{1},{k}_{2},\dots ,{k}_{n})$, with the probability

In our case, we supposed the relative phases ${\theta}_{p}$ with the same subscript $p$ among these single-photon states were identical but would be randomly regenerated for the next run of sampling. Repeatedly running the sampler time after time, we will get the average probability by integrating ${\theta}_{p}$ as

#### D. Coherent-State Sampling

As the name implies, the coherent-state sampler takes the input as a product of coherent states ${S}_{Cs}=|{e}^{i{\chi}_{1}}|{\alpha}_{1}|,{e}^{i{\chi}_{2}}|{\alpha}_{2}|,\dots ,{e}^{i{\chi}_{N}}|{\alpha}_{N}|\u27e9$, where $|{\alpha}_{j}|$ is the amplitude of the state, and ${\chi}_{j}$ is the relative phase randomly chosen from 0 to $2\pi $. To implement the coherent-state sampler, one only needs a source of coherent light from a laser which is split into $N$ modes, which results in the required value of ${\alpha}_{j}$. This is, generally speaking, easier than preparing many indistinguishable single-photon states. After propagating through the linear scattering unitary, the transformation of coherent states is

The state of the output ports is still a product of the coherent states. After renormalization, the final state is

The probability to detect the output Fock state ${T}_{Cs}=({t}_{1},{t}_{2},\dots ,{t}_{N})$ is

After averaging over the random phase ${\chi}_{j}$, the probability becomes

Using the probabilities in Eqs. (3), (4), (8), and (12), we have simulated boson sampling, classical sampling, mean-field sampling, and coherent-state sampling in the cases of two and three photons averaging over 100 Harr-random unitary matrices, as depicted in Fig. 1. The counting statistics of boson sampling is similar to mean-field sampling, while classical sampling obviously behaves higher on coincidence outcomes but lower on those with multiply populated modes. This phenomenon implies that mean-field sampling and coherent-state sampling abide by the bunching tendency of boson sampling. In the simulation, we set the input state of the coherent-state sampler to ${\alpha}_{j}=1$ on those modes that would have a single photon input and ${\alpha}_{j}=0$ for those modes that would have been a vacuum. As the coherent state has an indefinite number of photons, there are many events output from the coherent-state sampler that are not exactly two-photon or three-photon events in the respective simulations. The coherent-state sampler therefore yields many events containing more or less photons in total, and these events are not shown in Fig. 1. When compared with boson sampling, the counting statistics of the events shown for the coherent-state sampler is consequently lower.

## 3. CERTIFICATION PROTOCOL

#### A. Certification via the $\mathsf{N}$ Photon Single-Photon Input State

Since the certification based on bosonic bunching does not distinguish between boson sampling and the mean-field sampling, Tichy *et al.* [21] proposed a scheme using a symmetric Fourier matrix and a particular initial state with cyclic symmetry. In their scheme, some forbidden outcomes in the boson sampling will occur in the mean-field sampling. However, their scheme has the obvious disadvantage that it must change the sampling method and deviates from the standard boson sampling construction.

Suppose there is a sampling algorithm where the output can be measured, and the input is a linear unitary interaction of single photons and vacuum states. It is assumed that no details of the internal workings of the implementation of the algorithm are known. Is it possible to identify whether this sampler is a fully quantum boson sampler or a possibly classical mean-field sampler? In this paper, we present numerical evidence that inputting a test state that is constructed from single photons can solve this problem. The test state of the boson sampler is ${S}_{T}=(1,1,\dots ,1)$, whose mode arrangement is $\mathrm{j\u20d7}=(1,2,\dots ,N)$. A correspondingly complex test state for the mean-field sampler is $N$ runs using a single-photon state prepared as $|\psi \u27e9=\frac{1}{\sqrt{N}}\sum _{l=1}^{N}{e}^{i{\theta}_{l}}{\widehat{a}}_{l}^{\u2020}|0\u27e9$. In this case, the sum term in Eq. (9) becomes $\sum _{l=1}^{N}{|{\mathrm{\Lambda}}_{{k}_{q},l}|}^{2}$, which is equal to 1 as $\mathrm{\Lambda}$ is unitary. Thus, given the test state as the input state of a mean-field sampler, the outcome ${T}_{Mf}$ will be detected with the probability

From this expression, we can see that the probability distribution is not a function of the unitary matrix $\mathrm{\Lambda}$, but only a function of the number of photons $n$ and the output configuration ${T}_{Mf}$. Since the number of photons is equal to the number of input modes, the probability distribution is determined only by the dimension of the unitary matrix. The invariance of this distribution with respect to changes in the unitary matrix suggests that there exists an efficient method for discrimination between the mean-field sampler and the fully quantum boson sampler [17]. Here we study numerically a scenario with small numbers of photons where the entire probability distribution can be feasibly studied.

In Fig. 2, we randomly selected four $4\times 4$ unitary matrices by the Harr measure and numerically simulated boson sampling and mean-field sampling with the test input state ${S}_{T}=(1,1,1,1)$. The number of samples simulated for all situations is 10,000. The counting statistics of the mean-field sampling present as relatively uniform when compared with boson sampling in agreement with the theory above. The sample size here is not large enough for the effect of the randomly varying relative phase ${\theta}_{l}$ to completely vanish, so the counting statistics of mean-field sampling do not completely satisfy the probability given in Eq. (14). Fortunately, the sampling statistic of the mean-field sampler approximates the uniform distribution in the case of an arbitrary unitary matrix, while the sampling statistics of the boson sampler will vary for different matrices. For a particular matrix, the difference between the maximum event and the minimum event in boson sampling statistics is larger than those in the relatively uniform mean-field sampling. By using this particular test state as an input, the distributions of boson sampling and mean-field sampling can be distinguished without needing to modify the sampler, even without any prior knowledge about the linear scattering matrix.

In fact, when returning to the case where the input state does not have all modes occupied by a single photon, the most likely event still occurs in the boson sampling distribution and not the mean-field distribution when a large number of samples are taken, as shown in Fig. 3. One can firmly believe that the highest sampled event is from the boson sampling distribution. We conjecture the reason is that random phase of input states for the mean-field sampler limits the maximum value in probability distribution. However, this method may need a number of samples that depends exponentially on the number of input photons, and so this may be an inefficient protocol to distinguish the boson sampling and mean-field distributions.

Our verification protocol can be compared and contrasted with those described in the introduction [17,21]. One can see that this protocol is different from these previous ones but has some related structures. The protocol in [17] tests if the output distribution is uniform in the 0 and 1 subspace of detections or if it is a true boson sampling distribution proportional to matrix permanents. This protocol is efficient and stays completely within the boson sampling model of $n$ photons into $N>{n}^{2}$ modes. It can, however, be efficiently deceived with the classical sampling distribution. The protocol of [21] deviates from the standard boson sampling model but, given some assumptions about the nature of the linear interaction, is an efficient test of a necessary condition to perform boson sampling. In particular, it has the ability to perform non-classical interference. Our protocol also deviates from the standard boson sampling, but in a much simpler way than that in [21]. We add more photons to the input state so that all input modes contain single photons. In this case, one is now in the regime of [17], where to discriminate the mean-field sampler and boson sampling, one must distinguish uniformity from non-uniformity. When applying our criterion to a distribution that does scatter into $N>{n}^{2}$ modes, the distinction between the maximal and minimally sampled events still exists, but as stated above, one may need to gather an exponentially scaling number of samples to distinguish these distributions. Depending on the cost of single photons relative to the modes, it may be desirable to work in the small photon number regime and take an exponential number of samples to perform our test instead of using the $N=n$ case. It also must be remembered that no single test can efficiently verify the device’s proper functioning as a boson sampler, as this is believed to not be possible [5].

#### B. Robustness for the Partially Distinguishable Boson Sampling

Boson sampling in practice is expected to be influenced by experimental inaccuracies in the inputting, scattering, and detecting of single photons. We now consider the case of imprecision in preparing input states. In particular, we consider input states to the boson sampler constructed from partially distinguishable photon states. In this analysis, the degree of freedom used to induce distinguishability is the mutual delay of injected photons. The theory we use to describe this effect considers orthogonal spatio-temporal modes. By performing a Gram–Schmidt orthogonalization on the spatio-temporal degrees of freedom of the input state, the single particle in each mode can be represented as

This equation can be qualitatively understood as mode 1 defining a basis into which all other modes are compared, so the first mode is the “pure” mode. The second mode then has some common overlap with the first mode and hence can be divided up into two components that overlap with the first mode and some other orthogonal components that result in the sum over the modes only being required up to two. The third mode in which a photon is created will then have some components that overlap with the first and second modes and then an orthogonal component, which means there will be three terms in this sum, and so on in this manner until all photons have been expanded. The terms that have identical subscripts ${t}_{{k}_{r}}$ describe an input state of temporally indistinguishable photons, which relates to a standard boson sampling. The term with the different subscript ${t}_{{k}_{r}}$ corresponds to a classical sampling with distinguishable particles. The complex case is terms with some identical and some different ${t}_{{k}_{r}}$, which relates to the partially distinguishable boson sampling. Adjusting the coefficient ${C}_{r,{k}_{r}}$, the input state gives rise to a sampling in the range from fully indistinguishable to fully distinguishable particles.

By injecting the test state ${S}_{T}=(1,1\dots 1)$, we still can differentiate the mean-field sampling from the partially distinguishable boson sampling. As shown in Fig. 4, we simulate the partially distinguishable boson sampling and the mean-field sampling with the test state ${S}_{T}=(1,1,1)$. For the partially distinguishable boson sampling, the fluctuations vary considerably depending on random matrices, while mean-field sampling presents a stable distribution.

## 4. CONCLUSIONS

Boson sampling with identical bosons provides evidence that quantum computers can outperform classical computers. However, its verification is known to be difficult due to its computational complexity being outside the class $NP$. In our paper, we have proposed a certification scheme for the specific case of discriminating the boson sampler from the mean-field sampler. Without any limitation on the scattering matrices, the protocol only needs to prepare a test state and primitively analyze the output probability distribution. We have numerically simulated samples of the boson sampling and the mean-field sampling distributions when inputting the test state. For the output probability distributions, the statistics of boson sampling fluctuate more significantly than mean-field sampling. Considering the influence of imperfect input states, we have shown the certification capability of the protocol for partially distinguishable boson sampling against mean-field sampling.

In summary, we describe an experimental approach to verify the actual boson sampling against a physically plausible imposter. Compared to previous methods, our scheme retains a scattering matrix randomly chosen by the Harr measure and continues to work under the noise of partial distinguishability caused by non-simultaneously arriving photons. This work does not address the efficiency problem of our scheme. Optimizing the scheme to reduce the number of samples needed will be the direction of our future work.

## Funding

National Natural Science Foundation of China (NSFC) (11475160, 61575180); Natural Science Foundation of Shandong Province (ZR2014AQ026, ZR2014AM023); Fundamental Research Funds for the Central Universities of China (201313012); China Scholarship Council (CSC) (201306330035); National Key Laboratory of Electromagnetic Environment (201500005); Australian Research Council (ARC) (CE110001027).

## REFERENCES

**1. **P. W. Shor, “Algorithms for quantum computation: discrete logarithms and factoring,” in *Proceedings of the 35th Annual Symposium on Foundations of Computer Science* (IEEE, 1994), pp. 124–134.

**2. **P. W. Shor, “Polynomial-time algorithms for prime factorization and discrete logarithms on a quantum computer,” SIAM Rev. **41**, 303–332 (1999). [CrossRef]

**3. **M. A. Nielsen and I. L. Chuang, *Quantum Computation and Information* (Cambridge University, 2000).

**4. **E. Martin-Lopez, A. Laing, T. Lawson, R. Alvarez, X. Q. Zhou, and J. L. O’Brien, “Experimental realization of Shor’s quantum factoring algorithm using qubit recycling,” Nat. Photonics **6**, 773–776 (2012). [CrossRef]

**5. **S. Aaronson and A. Arkhipov, “The computational complexity of linear optics,” in Proceedings of the Forty-Third Annual ACM Symposium on Theory of Computing (STOC), New York, 2011, pp. 333–342.

**6. **M. A. Broome, A. Fedrizzi, S. Rahimi-Keshari, J. Dove, S. Aaronson, T. C. Ralph, and A. G. White, “Photonic boson sampling in a tunable circuit,” Science **339**, 794–798 (2013). [CrossRef]

**7. **J. B. Spring, B. J. Metcalf, P. C. Humphreys, W. S. Kolthammer, X. M. Jin, M. Barbieri, A. Datta, N. Thomas-Peter, N. K. Langford, D. Kundys, J. C. Gates, B. J. Smith, P. G. R. Smith, and I. A. Walmsley, “Boson sampling on a photonic chip,” Science **339**, 798–801 (2013). [CrossRef]

**8. **M. Tillmann, B. Dakic, R. Heilmann, S. Nolte, A. Szameit, and P. Walther, “Experimental boson sampling,” Nat. Photonics **7**, 540–544 (2013). [CrossRef]

**9. **A. Crespi, R. Osellame, R. Ramponi, D. J. Brod, E. F. Galvao, N. Spagnolo, C. Vitelli, E. Maiorino, P. Mataloni, and F. Sciarrino, “Integrated multimode interferometers with arbitrary designs for photonic boson sampling,” Nat. Photonics **7**, 545–549 (2013). [CrossRef]

**10. **J. Carolan, C. Harrold, C. Sparrow, E. Martín-López, N. J. Russell, J. W. Silverstone, P. J. Shadbolt, N. Matsuda, M. Oguma, M. Itoh, G. D. Marshall, M. G. Thompson, J. C. F. Matthews, T. Hashimoto, J. L. O’Brien, and A. Laing, “Universal linear optics,” Science **349**, 711–716 (2015). [CrossRef]

**11. **K. R. Motes, A. Gilchrist, J. P. Dowling, and P. P. Rohde, “Scalable boson sampling with time-bin encoding using a loop-based architecture,” Phys. Rev. Lett. **113**, 120501 (2014). [CrossRef]

**12. **Y. He, Z.-E. Su, H.-L. Huang, X. Ding, J. Qin, C. Wang, S. Unsleber, C. Chen, H. Wang, Y.-M. He, X.-L. Wang, C. Schneider, M. Kamp, S. Höfling, C.-Y. Lu, and J.-W. Pan, “Scalable boson sampling with a single-photon device,” arXiv:1603.04127 [quant-ph] (2016).

**13. **P. P. Rohde and T. C. Ralph, “Error tolerance of the boson-sampling model for linear optics quantum computing,” Phys. Rev. A **85**, 022332 (2012). [CrossRef]

**14. **A. P. Lund, A. Laing, S. Rahimi-Keshari, T. Rudolph, J. L. O’Brien, and T. C. Ralph, “Boson sampling from a Gaussian state,” Phys. Rev. Lett. **113**, 100502 (2014). [CrossRef]

**15. **M. Bentivegna, N. Spagnolo, C. Vitelli, F. Flamini, N. Viggianiello, L. Latmiral, P. Mataloni, D. J. Brod, E. F. Galvãoand, A. Crespi, R. Ramponi, R. Osellame, and F. Sciarrino, “Experimental scattershot boson sampling,” Sci. Adv. **1**, e1400255 (2015). [CrossRef]

**16. **C. Gogolin, M. Kliesch, L. Aolita, and J. Eisert, “Boson-sampling in the light of sample complexity,” arXiv:1306.3995v2 (2013).

**17. **S. Aaronson and A. Arkhipov, “BosonSampling is far from uniform,” Quantum Inf. Comput. **14**, 1383–1423 (2013).

**18. **N. Spagnolo, C. Vitelli, M. Bentivegna, D. J. Brod, A. Crespi, F. Flamini, S. Giacomini, G. Milani, R. Ramponi, P. Mataloni, R. Osellame, E. F. Galvao, and F. Sciarrino, “Experimental validation of photonic boson sampling,” Nat. Photonics **8**, 615–620 (2014). [CrossRef]

**19. **J. Carolan, D. A. MeineckeJasmin, P. J. Shadbolt, N. J. Russell, N. Ismail, K. Wörhoff, T. Rudolph, M. G. Thompson, J. L. O’Brien, J. C. F. Matthews, and A. Laing, “On the experimental verification of quantum complexity in linear optics,” Nat. Photonics **8**, 621–626 (2014). [CrossRef]

**20. **J. Carolan, J. Meinecke, P. Shadbolt, N. J. Russell, N. Ismail, K. Wörhoff, T. Rudolph, M. Thompson, J. L. O’Brien, J. Matthews, and A. Laing, “Verifying quantum complexity in linear optical experiments,” in *Conference on Lasers and Electro-Optics* (Optical Society of America, 2014), paper FM2A.7.

**21. **M. C. Tichy, K. Mayer, A. Buchleitner, and K. Mølmer, “Stringent and efficient assessment of boson-sampling devices,” Phys. Rev. Lett. **113**, 020502 (2014). [CrossRef]

**22. **A. Crespi, R. Osellame, R. Ramponi, M. Bentivegna, F. Flamini, N. Spagnolo, N. Viggianiello, L. Innocenti, P. Mataloni, and F. Sciarrino, “Suppression law of quantum states in a 3D photonic fast Fourier transform chip,” Nat. Commun. **7**, 10469 (2016). [CrossRef]

**23. **M. Jerrum, A. Sinclair, and E. Vigoda, “A polynomial-time approximation algorithm for the permanent of a matrix with non-negative entries,” in Proceedings of the Thirty-Third Annual ACM Symposium on Theory of Computing (STOC), New York, 2001, pp. 712–721.

**24. **M. C. Tichy, H. T. Lim, Y. S. Ra, F. Mintert, Y. H. Kim, and A. Buchleitner, “Four-photon indistinguishability transition,” Phys. Rev. A **83**, 062111 (2011). [CrossRef]