Photonic implementation of boson sampling: a review

Abstract. Boson sampling is a computational problem that has recently been proposed as a candidate to obtain an unequivocal quantum computational advantage. The problem consists in sampling from the output distribution of indistinguishable bosons in a linear interferometer. There is strong evidence that such an experiment is hard to classically simulate, but it is naturally solved by dedicated photonic quantum hardware, comprising single photons, linear evolution, and photodetection. This prospect has stimulated much effort resulting in the experimental implementation of progressively larger devices. We review recent advances in photonic boson sampling, describing both the technological improvements achieved and the future challenges. We also discuss recent proposals and implementations of variants of the original problem, theoretical issues occurring when imperfections are considered, and advances in the development of suitable techniques for validation of boson sampling experiments. We conclude by discussing the future application of photonic boson sampling devices beyond the original theoretical scope.


Introduction
In recent years, much work has been directed to the development of suitable technologies for scalable quantum computation. This effort is motivated by the promise of quantum algorithmic speed-up, with a notable early example provided by Shor's algorithm for integer factoring. 1 Despite the tremendous advances in quantum technologies [2][3][4][5][6][7] reported in the last few years, the implementation of a large-scale universal quantum computer is still far from our current capabilities. Hence, intermediate milestones need to be identified in the longterm effort toward harnessing the computational potential of quantum systems. A first fundamental step in this direction would be the achievement of the regime of quantum computational supremacy, 8 i.e., the experimental demonstration of a quantum device capable of performing a computational task unambiguously faster than present-day classical computers.
Within this framework, Aaronson and Arkhipov (AA) 9 formulated a well-defined computational problem that they called boson sampling. This problem consists in sampling from the output distribution of n indistinguishable bosons that interfere during the evolution through a Haar-random-chosen linear network. In Ref. 9, strong evidence was provided that boson sampling is an intractable problem for classical computers, as it is related to the evaluation of permanents of matrices with complex entries (a problem known, in computational complexity terms, to belong to the #P-hard class 10 ). On the other hand, boson sampling can be tackled with a dedicated quantum hardware, which, despite not being universal for quantum computation, is capable of implementing the required dynamics. To this end, one of the most suitable platforms is provided by photonic systems, as the necessary elements (sources, linear evolution, and detection) are available with present technology. 2 Given the lack of error-correction-an issue shared by all current proposals for quantum computational supremacy 8 -it is an open question whether current technology is capable of scaling boson sampling to arbitrarily large sizes while maintaining a quantum advantage, which has also motivated research in more robust variants of the model. [11][12][13][14][15][16][17] Following these theoretical developments, a strong experimental effort was then initiated to realize progressively larger instances of boson sampling experiments. This is leading to a race aimed at reaching the quantum advantage 18 regime in a photonic platform, namely the condition where the experiment is outperforming a classical computer. In parallel, theoretical work was required to define suitable methods to verify that a given device is sampling from the correct distribution. 19 This is indeed a relevant issue since the very complexity of the problem forbids the application of usual verification methods, and classical simulations become progressively intractable for the increasing system size.
In this paper, we discuss recent advances in the field of photonic boson sampling, with particular attention to experimental implementations and validation methods. The remainder of this paper is organized as follows. In Sec. 2, we first provide a theoretical overview on the problem, discussing the computational model, classical simulation algorithms, and conditions that turn the system amenable to efficient classical simulation. Then in Sec. 3, we discuss variants of the original task that have been proposed to improve the efficiency of the quantum simulator, without affecting the problem's complexity. In Sec. 4, we review the experimental implementations reported so far, discussing the employed photonic platforms. In Sec. 5, we provide an overview of the validation techniques for boson sampling that have been proposed and tested experimentally. In Sec. 6, we then discuss the scalability of photonic platforms toward achieving the quantum supremacy regime. Finally, we provide an outlook in Sec. 7, where we also discuss recently highlighted applications of boson sampling in contexts different from the original proposal.

Boson Sampling: Theoretical Overview
In this section, we review the theoretical underpinnings of the boson sampling problem. In Sec. 2.1, we define the problem and discuss its computational complexity. In Sec. 2.2, we review the known classical simulation algorithms, and in Sec. 2.3 we discuss some regimes, in which the classical simulation of boson sampling is known to become tractable.

Model
Consider a set of m modes with associated creation operators a † i , for i ¼ 1; …; m, satisfying bosonic commutation relations. A Fock state of n photons in these modes can be written as where s i are the non-negative integers that count the number of photons in each mode and P s i ¼ n. When s i ≤ 1 for all i, the state is a no-collision state, meaning that there is no individual mode containing two or more photons.
Consider now an m-mode linear-optical transformation (i.e., an interferometer) described by U ∈ SUðmÞ. Its action is fully determined by the evolution of the creation operators: U ij a † j : (2) It follows 9,20,21,22 that the transition probability between an input state jSi ¼ js 1 s 2 …s m i and an output state jTi ¼ jt 1 t 2 …t m i can be written as where U S;T is an n × n submatrix of U constructed by taking t i copies of the i'th row of U and s j copies of its j'th column, 9 and is the permanent 10 of matrix B, and S n is the symmetric group.
Since the permanent is, in general, hard to compute, 10 Eq. (3) underpins the complexity of a particular class of linear optical experiments. In complexity theory terms, the permanent is #P-hard, 10 and the best known classical algorithm for computing it, due to Ryser, takes Oðn2 n Þ steps for an n × n matrix. The connection to complexity theory was noted by Troyansky and Tishby 21 who attempted to leverage Eq. (3) to build a quantum-mechanical algorithm to compute permanents. They realized that their algorithm was not efficient, as exponentially many experimental samples would be needed to produce a sufficiently good approximation. AA explored this connection further by shifting to sampling problems, where the computational task is to produce a sample from some probability distribution sufficiently close to the ideal one. 9 We now outline the argument of AA, emphasizing the requirements it raises for experimental demonstrations of boson sampling. Further developments that attempted to soften the initial requirements are reviewed in Sec. 3.
The standard boson sampling setup 9 consists of the following ingredients (see Fig. 1).
(i) Preparation of a n-photon, m-mode input state, with each mode containing either zero or one photon. Without loss of generality, we can choose an input state with a single photon in each of the first n modes. Having more photons per input mode means choosing repeated columns in the submatrix of Eq. (3), which is likely to make the permanent easier to compute (in the situation where all photons are input in the same mode, the computation becomes trivial). (ii) An m-mode interferometer described by an mdimensional Haar-random unitary operator U, where m ¼ Oðn 2 Þ. An arbitrary interferometer can be built as a circuit with Oðm 2 Þ two-mode elements and depth OðmÞ. 23,24 Since only some input modes are occupied, it is possible to reduce this further to a circuit of OðmnÞ elements and depth Oðn log mÞ. 9 The randomness of U has two main purposes. The first is so U does not have any special structure that a classical simulation algorithm could exploit. The second is that, for Haar-random unitaries in Oðn 2 Þ modes, the outcomes are dominated by no-collision events due to the bosonic birthday paradox. 9,25 (iii) Photon detectors on every output. Given (ii), it suffices for them to be bucket detectors, i.e., to only distinguish between vacuum and nonvacuum states.
Let D be the probability distribution predicted by quantum mechanics for the no-collision outcomes of such an experiment, and let k · k represent the total variation distance (TVD). 9 The main result of AA states (up to two additional complexitytheoretic conjectures) that, if a classical algorithm existed that could produce a sample from any distribution D 0 such that in time polyðn; 1∕ϵÞ, the polynomial hierarchy (PH) would collapse to the third level. Since PH is a tower of complexity classes that is strongly believed to be infinite, the result of AA provides strong evidence against the possibility of an efficient classical simulation of boson sampling. The proposal of AA uses an approximate notion of simulation, as seen in Eq. (5). Thus one could attempt a classical simulation using an algorithm for approximating the permanent, rather than Ryser's exact algorithm. The problem with this approach is that it can be shown that the permanent would also be #P-hard to approximate on the worst case 9 to the required precision. One classical algorithm to approximate the permanent is due to Gurvits. 26 It takes time Oðn 2 ∕ϵ 2 Þ to compute the permanents of the n × n matrices in a boson sampling distribution to precision ϵ. This is not sufficient to give a simulation of boson sampling because there we typically have exponentially many permanents that are exponentially small, and so approximating them to a nontrivial precision via this algorithm would take exponential time.
In very broad strokes, the argument of AA is as follows. If U is a Haar-random matrix and m ¼ Oðn 2 Þ, it can be proven that n × n submatrices of U look approximately Gaussian. Since Gaussian matrices do not seem to possess any special structure that a classical algorithm could exploit, it is reasonable to conjecture that, with high probability, the permanent of a Gaussian matrix is among the hardest ones to compute (this is the permanent-of-Gaussians conjecture). However, if a classical algorithm could sample from a distribution ϵ-close in TVD to the ideal one, it would have to be approximating most probabilities very well. Using an algorithm attributable to Stockmeyer, 27 it would then be possible to solve a #P-hard problem inside BPP NP , which is a complexity class that resides inside the third level of the PH. Using a theorem due to Toda, 28 this would imply that the entire hierarchy would have to collapse to its third level. This result, thus, establishes that boson sampling is hard to simulate, even in an approximate sense, up to three conjectures that: (i) PH is infinite, (ii) permanents of Gaussian matrices are #P-hard to approximate, and (iii) permanents of Gaussian matrices are not too concentrated around zero. The third conjecture is more technical but is necessary to relate two versions of the problem of approximating permanents of Gaussian matrices: to within additive and multiplicative errors.

Simulation Algorithms
The AA argument provides evidence that any boson sampling simulation on a classical computer must take a time which grows exponentially with the number of photons n. It is important to find the best classical algorithms for boson sampling simulation. This will help to establish the goals for experimentally demonstrating unequivocal quantum computational advantage. From a more practical point of view, simulation algorithms are required for analyzing current, small-scale implementations. Optimal simulation algorithms also help clarify how inevitable experimental imperfections affect the computational power of the model; this includes partial photon distinguishability and photon loss. In this section, we review known classical simulation algorithms for boson sampling.
The brute-force approach for simulating boson sampling would involve calculating the probabilities associated with all input-output possibilities (in case we are using the variableinput, scattershot approach to boson sampling, see Sec. 3.1). Doing this for n photons in m ¼ n 2 modes with Ryser's algorithm, and for all n 2 n input-output combinations, would take Oðn2 n Þ n 2 n time steps. As regards space (memory) usage, it is not necessary to store all calculated permanents. 29 This bruteforce approach is computationally intensive in the extreme; we now review recently proposed algorithms that have considerably improved on this. In Ref. 18, two different new algorithms were proposed. The first is based on the well-known technique of rejection sampling: we sample from a simple (e.g., uniform) distribution over all events, then test whether the drawn event could be an output by computing its probability (which, in our case, is given in terms of a matrix permanent). This results in exact sampling, at a reasonable computational cost if the largest event probability is known, and not too large. The implementation of this algorithm for exact simulation of boson sampling requires an estimate of the largest event probability. If this maximum is underestimated, the simulation will be only approximate, and if this maximum is too large, it will result in a computational bottleneck due to the large number of rejected samples.
The second algorithm of Ref. 18 is a variant of the Markov chain Monte Carlo method known as metropolised independence sampling. The key idea involves constructing a Markov chain whose stationary distribution is the boson sampling distribution. To improve the performance, it is necessary to choose an initial, efficiently computable distribution as close as possible to the target distribution; Neville et al. 18 proposed to use the distribution corresponding to distinguishable photons as the initial distribution. The algorithm also requires a burn-in phase, in which the Markov chain approaches the stationary case. Using validation methods we review in Sec. 5, Neville et al. 18 showed that the algorithm seems to perform well in simulating a 20-photon boson sampling experiment on a standard laptop while enabling one to simulate a 30-photon, 900-mode apparatus on a local server.
Clifford and Clifford 30 proposed a breakthrough boson sampling simulation algorithm. Their result applies methods of hierarchical sampling commonly used in Bayesian computation. 31,32 The algorithm samples from the exact boson sampling distribution and requires only O½n2 n þ polyðm; nÞ steps to output each sample, where n is the number of photons and m is the number of modes. The time it takes to output each sample corresponds roughly to the time of computing two permanents of a general n × n matrix. This algorithm has the same complexity also for the scattershot variation of boson sampling (see Sec. 3.1). Using recent feasibility studies on the calculation of permanents of large matrices, 33 it is possible to estimate 18 that boson sampling problems will only become reasonably hard for classical computers for experiments with around n ¼ 50 photons.

Classically Simulable Limits of Boson Sampling
One important way to investigate the limitations of boson sampling is to search for physically realistic circumstances that might alter its computational capabilities, for example, by making it efficiently simulable. In this section, we review the known results regarding the complexity of boson sampling under realistic types of errors.
The most physically relevant experimental imperfection in boson sampling is photon loss. Formally, loss has been described in two alternative models. One model, more common in the quantum optics literature, replaces losses in linear-optical elements by beam splitters that have some probability of routing photons into unmeasured modes. 34 The second model considers the modification of Eq. (3) when a specific number of photons is lost, i.e., by averaging the final probability over all possible input photons that might have been lost. 35, 36 Brod and Oszmaniec 36 discuss in detail in which regimes these two models are computationally equivalent.
On the positive side, Aaronson and Brod 35 showed that boson sampling retains its computational complexity if the number of lost photons is constant. This is done by replacing the permanent with a different matrix function in Eq. (3) and showing that in this limit that function is as hard to compute as the permanent.
On the other hand, upper bounds on allowed photon losses are also known. If only Oðlog nÞ photons are left, Aaronson and Brod 35 claimed that boson sampling becomes efficiently classically simulable. Although this was claimed without proof, it follows immediately from the algorithm of Clifford and Clifford. 30 One can also provide an efficient classical simulation by showing that the lossy boson sampling state is sufficiently close to some well-known state that is itself classically simulable, such as thermal 37 or particle-separable states. 36 This approach was used to show that boson sampling becomes classically simulable when less than Oð ffiffiffi n p Þ photons remain.
In a very recent development, Renema et al. 38 proposed an algorithm that gives an efficient classical simulation when a constant fraction of photons remain. Although the previous results all regard lossy boson sampling in some asymptotic limit, the latter would be directly applicable to finite-size experiments. It is also possible to combine the effects of losses with other imperfections, such as done by Rahimi-Keshari et al. 34 Using phase-space methods (i.e., based on positive quasiprobability distributions), they showed that boson sampling becomes classically simulable under a combination of (sufficiently high) constant rates of losses and dark counts.
Another relevant experimental imperfection that has been extensively considered in the boson sampling literature is partial photon distinguishability, e.g., due to some internal unmeasured degree of freedom such as frequency, polarization, and arrival time. Already in their original paper, AA pointed out that if all photons are completely distinguishable boson sampling becomes classically simulable. Interestingly, the transition probabilities for distinguishable photons are also given by permanents, but in this case they are permanents of matrices with only positive entries and can be approximated classically in polynomial time. 39 Intermediate regimes of distinguishability have been considered by several authors. [40][41][42][43] In particular, the formalism of Tichy 40 has been used to describe an algorithm for simulating boson sampling when the photons have high-pairwise distinguishability. 44 From the other side, it is also known that if partial pairwise distinguishability decreases sufficiently fast as the number of photons increases, 45 the complexity of boson sampling is unaffected by it.
Beyond losses and partial distinguishability, boson sampling has also been investigated under the effects of other sources of noise, such as fabrication imperfections 46,47 and Gaussian noise in the output data. 48 Finally, very recently, a tensor network approach for the simulation of noisy quantum circuits has been reported in Ref. 49. More specifically, such a simulation approach leads to the consequence that most circuits with constant levels of noise can be simulated classically with polynomial resources. Although this approach has been tailored for quantum circuits, the application of this method can be potentially envisaged to other quantum many-body systems including boson sampling. modes; the use of different input quantum states and detectors different from single-photon bucket detectors. In this section, we review some variants of the boson sampling problem (see Fig. 2).

Scattershot Boson Sampling
The initial implementations of boson sampling used probabilistic sources based on spontaneous parametric down-conversion (SPDC). One of the two generated photons is detected, heralding the other photon. Those SPDC single-photon sources suffer from two drawbacks. The first is the probabilistic nature of the source; the second is the requirement of low pump intensity, so as to avoid higher-order nonlinear effects, which can generate more than two photons per pulse.
Scattershot boson sampling 11 is a variant of the boson sampling problem, which is particularly suitable to implementation using SPDC sources. Simply put, we connect k (k > n) such single photon, SPDC sources to different input ports of the linear interferometer [see Fig. 2(a)]. The expected number of single photons per pulse is increased by a factor k n , which for k ≫ n represents an exponential improvement in generation rate with respect to the fixed input version of boson sampling. This scattershot boson sampling set-up naturally solves a version of the boson sampling problem, in which the inputs are chosen uniformly at random, and we are required to sample from the output distribution of this interferometer. Note that this variant does not require an increase in the pump laser power, as the laser can pump the SPDC sources sequentially, with little loss to down-converted photons.
One variant of scattershot boson sampling is the so-called driven boson sampling of Barkhofen et al. 13 The key idea is to add k beam splitter layers prior to the m-mode Haar-random interferometer [see Fig. 2(b)]. These beam splitters can be used to couple up to k × m heralded, SPDC single-photon sources to the main interferometer, increasing the n-photon event rate by a factor of k with respect to the original scattershot proposal with m sources.

Gaussian Boson Sampling
The scattershot boson sampling variant described in Sec. 3.1 is already an example of a variant using different input states-in that case, Gaussian states that, upon measurement of a single photon, herald single photons. Hamilton et al. 14,50 analyzed the advantages and complexity of using Gaussian states directly as input states in a linear interferometer.
Hamilton et al.'s proposal of Gaussian boson sampling 14,50 relies on inputs of single-mode squeezed states, linear interferometers, and single-photon detectors at the output [see Fig. 2(c)]. They showed that the output statistics is given by the hafnian 20 of a matrix that includes information about both the linear interferometer and the Gaussian input state. Like the permanent, calculating the hafnian is a problem in the #P complexity class. This made it possible for them to define a variant of the boson sampling problem called Gaussian boson sampling. By relying on some conjectures about hafnians, which mirror the conjectures about the permanent used in the usual boson sampling problem, they showed that also this problem is intractable for classical computers. Classical benchmarking of Gaussian boson sampling was also performed with the Titan supercomputer. 51 In that paper, it was shown that classical simulation of a 800-mode Gaussian boson sampling experiment where 20 output detectors over 800 modes provide a click can be performed in ∼2 h.
Gaussian boson sampling uses all the squeezed photons, rather than using half of them to herald the other half. This advantage, however, is offset by the fact that hafnians are slightly easier to calculate. 50,52 Gaussian boson sampling allows for the use of higher pump power, which in scattershot boson sampling had to be kept low to avoid multiple-photon generation. This variant also inherits the exponential advantage in event rate of the scattershot variant. In comparison with scattershot boson sampling, Gaussian boson sampling events form a much smaller sampling space, as variable inputs are not required. This results in an experimental advantage over the scattershot variant if one wants to completely characterize the probabilities associated with each event. For larger implementations, however, this is not a tenable goal; we will discuss the problem of boson sampling validation in Sec. 5.
By considering the time-reversed version of Gaussian boson sampling, we see that a variant of boson sampling with singlephoton inputs and Gaussian basis measurements is exactly as hard to simulate as Gaussian boson sampling. This variant has been recently analyzed in detail. 53,54 Furthermore, a variant of Gaussian boson sampling with threshold detectors has been also investigated. 55 Some possible computational applications of Gaussian boson sampling, which we will revisit in Sec. 7, include finding dense subgraphs 56 and perfect matching of graphs, 57 improving the performance of stochastic optimization algorithms, 58 and simulation of molecular dynamics. 59,60 Gaussian boson sampling has also been linked to the dynamical Casimir effect in multimode waveguide resonators. 61

Other Variants
Besides scattershot and Gaussian boson sampling, other variants of the task have been theoretically proposed and investigated. In Refs. 16 and 17, the sampling complexity has been analyzed for photon-added and photon-subtracted states, and different detection schemes [see Fig. 2(d)]. Depending on the specific configuration, the complexity is maintained for all system parameters or shows a transition toward a classically simulable regime when using a large number of photons.
Finally, the complexity of multiboson correlation sampling with photons having no spectral overlap has been analyzed in Ref. 15. In this case, the output photons after evolution through a random linear network are measured by a timeand polarization-resolved correlation detection apparatus [see Fig. 2(e)]. In this case, the complexity is provided by interference of all n! paths that are indistinguishable to the detectors. This variant of the original problem has also provided an insight on the role of the measurement type in the particle indistinguishabilty regime. 43 Finally, a scattershot version of multiboson correlation sampling has been proposed in Ref. 62, enabling enhanced events rate by sampling also from time and frequency modes of the input photons.

Experimental Boson Sampling
To perform a standard boson sampling experiment with photons, one needs the following three building blocks: single-photon sources, unitary transformations (representing possible linear interferometers), and single-photon detectors. In Fig. 3, we summarize the technical solutions adopted in the experimental demonstrations of boson sampling so far. A list of the experiments reported in the literature is provided in Table 1, where the main features of the three building blocks are also detailed.
While diverse in the experimental details, the first demonstrations of boson sampling [63][64][65][66][67]69,76 share the main features of their experimental platforms. In particular, in all these experiments, several identical photons are generated by SPDC in nonlinear crystals [ Fig. 3(a) In order to produce more than two identical photons, the pump intensity of the SPDC process is increased so that the probability of simultaneous generation of two or more pairs becomes significant. If the different pairs of photons are collected on distinct spatial or polarization modes, the photons can easily be directed to separate input ports of the linear optical circuit. However, emission of multiple pairs is exponentially less efficient than the generation of a single pair. In fact, in the mentioned experiments, at most four identical photons were generated on different modes, pumping one 63,65,66,69,76 or two 64,67 distinct nonlinear crystals. In addition, it should be noted that in this process the emission of two pairs on the same two modes (an outcome that is not desired) has the same probability of the emission of two pairs on four distinct modes. The correct event is identified by employing one of the four output photons as trigger while only the remaining three are injected in the interferometer, and postselection is operated on the overall detection of four photons. This is a relevant aspect since multiple emissions of photons on the same mode must be avoided in boson sampling experiments. In this case, transition amplitudes are related to matrix permanents with repeated rows that can be easier to compute. As a result, all these first experiments report at most three-photon genuine-boson-sampling instances. In Ref. 67, the output distribution of four photons on different modes is reconstructed a posteriori from a probabilistic input.
The implementation of the unitary transformation was done with a guided-wave architecture, which enables one to cascade interferometers with compactness and stability impossible to achieve by a conventional bulk-optics approach. While boson sampling experiments have been reported also making use of fiber splitters 63,70 or arrays of parallel and continuously coupled waveguides 67 [Fig. 3(g)], the most common layout employs an integrated network of directional couplers (the integrated version of the bulk beam splitter) and phase shifters. Indeed, by arranging such discrete elements in a triangular 23,66,69 or rectangular 24 layout, with judiciously chosen values of the reflectivities and phase shifts, it is possible to realize any arbitrary linear transformation of the optical modes. In Ref. 77, it was also shown how to construct an operational approach to directly implement Haar-random matrices from appropriate distributions for the optical components. One should note that the capability of producing arbitrary transformations is definitely needed to extract and implement Haar-random unitaries, and thus it is essential for the demonstration of a real quantum advantage in a boson sampling experiment.
Up to now, the realization of arbitrarily chosen transformations by the mentioned triangular-network approach has been experimentally reported in femtosecond laser written 66 or silica-on-silicon 69 waveguide circuits. The femtosecond laser writing technology 78 gives the unique possibility to realize three-dimensional (3-D) waveguide circuits; 79 the third dimension was exploited in Ref. 66 to deform the waveguide paths out of the plane and thus independently adjust phases and reflectivities within the interferometer. In particular, such an innovative design enabled the realization of an arbitrarily chosen five-mode transformation [ Fig. 3(e)]. More recently, a six-mode interferometer with full reconfiguration capabilities was demonstrated in Ref. 69, using state-of-art silica-on-silicon lithographic technology. The authors developed an integrated interferometer provided with 30 electronically controlled thermo-optic phase shifters, which allow the dynamical tuning of the circuit to obtain any six-mode linear interferometer. In that device, the individual splitting elements in the triangular network are replaced by Mach-Zehnder interferometers: in this way, it is possible to control the equivalent directional coupler's reflectivity by acting on the internal phase of the Mach-Zehnder.
As a matter of fact, all the experiments discussed above suffer from a few common and very relevant limitations, which hamper an effective scaling of photonic boson sampling up to the level of a provable quantum advantage over current classical computers. As a first point, the SPDC generation scheme presented above hardly allows the scaling of the number of input photons beyond the small numbers already demonstrated. One possible experimental route to achieve multiphoton input states with higher efficiency while still relying on nonlinear frequency conversion, is the scattershot approach [ Fig. 3(b)], whose theory was discussed in Sec. 3.1. An experimental demonstration of this method was given in Ref. 68. There, six sources of down-converted photon pairs (implemented using three distinct BBO crystals) operate in parallel in order to inject three-photon states into integrated interferometers with up to 13 modes. One of the sources feeds a photon pair to a fixed couple of inputs. The other five sources produce random, but heralded, single photons that are coupled to other input ports of the interferometer. Detection of three photons at the output, in coincidence with one heralding event, postselects the occurrence of a proper input state with three photons in separate modes. This technique enabled a 4.5-times enhanced generation rate of three-photon states with respect to the multipair SPDC scheme adopted previously. In a recent paper, 74 the scattershot approach has been scaled up to 12 sources, located in six nonlinear crystals. Each of the crystals can emit probabilistically up to two heralded photons with horizontal and/or vertical polarization, multiplexed in the same spatial mode. The sources were then connected to the input ports of an interferometer consisting of six spatial modes. By considering the polarization degree of freedom, this experimental scheme implemented a 12-mode transformation composed by two 6 × 6 unitary-independent blocks. This led to the observation of three, four, and fivephoton events thanks to the increased signal rate provided by the scattershot approach. The measured event rates were, respectively, 3.9 kHz, 44 Hz, and 0.3 Hz. When trying to scale up the number of sources, a bulk implementation of SPDC will eventually reach a limit size, whereas a more promising platform for larger experiments can be provided by the adoption of arrays of integrated sources 80-83 and multiplexing. 84 Demonstrations of arrayed integrated sources have indeed been reported using different technologies and approaches, including spontaneous-four-wave-mixing (SFWM) in silicon photonics 75,83 and glass waveguides, 80 or SPDC in nonlinear waveguides. 81,82 Very recently, scattershot and Gaussian boson sampling in a silicon chip was reported online, 75 in a device that includes four waveguide-integrated photon-pair sources, based on nondegenerate SFWM. These sources consist in spiralling waveguides (to increase the nonlinear interaction region), where two-mode vacuum squeezing is produced via a single-wavelength laser pump [ Fig. 3(c)]. In the chip, the sources are directly connected to an array of 12 evanescently coupled waveguides, enabling detection of 8 photons in the chip corresponding to four-photon boson sampling experiments. The measured rates were four events/hour for the scattershot regime, for an overall number of ∼20 recorded events and 1.1 Hz for the Gaussian scenario.
A different strategy involves the use of quantum-dot sources to generate trains of identical single photons 70-73 with high efficiency. Such quantum-dot sources are based on single-defect states in the semiconductor lattice, localized within a resonant microcavity [ Fig. 3(d)]. A laser pulse excites the defect state, which, upon spontaneous decay, emits a single photon. These defect states are discrete and well-defined in energy: consecutive excitation pulses thus result in the emission of highly indistinguishable single photons. The use of mode-locked oscillators as pump lasers enables single-photon trains at tens of megahertz repetition rates from the quantum dot. To efficiently route the photons into different spatial optical modes, i.e., distinct inputs of the unitary circuit, active time-to-space demultiplexing [71][72][73] is preferably adopted (passive demultiplexing has also been used, 70 but it is not efficient when scaling up the number of photons). With this kind of source, boson sampling experiments with up to five detected photons have been demonstrated, 72,73 with five-photon events rates of 4 Hz 72 and 0.78 kHz, 73 respectively. In those works, active demultiplexing was performed by a chain of Pockels cells. However, the switching frequency of those bulk-optics devices did not exceed the megahertz. Miniaturized guided-wave electro-optic switches, 85 fabricated in nonlinear substrates, might enable scaling up the operating frequencies to match the pulse frequency of mode-locked oscillators. Indeed, it is still a challenge to produce such devices with reasonably low insertion losses. As a further aspect to note, semiconductor quantum dot sources typically emit at near-infrared wavelengths around 900 nm, close to the sensitivity limits of silicon detectors. Superconducting nanowire single-photon detectors (SNSPDs), see Fig. 3(j) may thus be a convenient choice in this case. 71,73 A second important problem of the first experiments of photonic boson sampling is the relevant insertion loss of the circuit implementing the linear interferometer: in fact, the typical transmission of an integrated-optics device with reasonable complexity and 10 to 20 modes is in the order of 30%. Losses are a very relevant problem as the number n of photons increases. Indeed, the success rate of the boson sampling experiment scales as the circuit losses to the power of n. Two innovative experimental implementations have emerged recently that may allow one to substantially reduce the loss of the interferometric apparatus.
Wang et al. 72,73 proposed and demonstrated boson sampling using multimode interferometers operating in free-space that consist in micro-optical assemblies. These are built of several fused-quartz trapezoids, covered with differently reflecting optical coatings, stacked together [see Fig. 3(f)] to produce a device with the size of a few centimetres. Collimated light beams, injected at the sides, experience multiple reflections and interference inside the structure, each interface between two trapezoids corresponding to a row of beam splitters. Such structures have been reported with geometries that mimic both the triangular 72 and the rectangular 73 interferometric layouts that would allow one, in principle, to implement arbitrary linear-optical transformations. However, while the achieved experimental results prove the optimal transmission (over 99%) and high phase-stability of these structures, it is not clear whether they could be adopted to realize linear-optical transformations that are truly Haar-random, as the free parameters of the circuit appear to be quite reduced. In fact, the reflectivity of one interface defines the characteristics of an entire row of beam splitters, and phases also cannot be controlled punctually within a given trapezoid.
A fundamentally different approach was proposed theoretically in Ref. 86 and recently reported experimentally in Ref. 71, and relies on the use of temporal modes instead of spatial ones. The available modes are a train of time bins separated by an interval τ. The unitary transformation is implemented by means of two different low-loss fiber loops where the photons can be routed by means of rapidly switching modulators [see Fig. 3(h)]. The smaller loop, introducing a τ delay on the injected photons, enables one to interfere adjacent time bins. The larger loop, instead, produces a delay much larger than τ and feeds back the circuit with the same photon packets, thus allowing for successive interference steps. It is easily shown that with this scheme one can implement the same beam-splitter networks that, in space, allow for arbitrary linear interferometers. 86 Such an apparatus was used, in conjunction with a quantum dot source, to realize a boson sampling experiment with four photons interfering in eight modes. 71 Notably, in this approach, the photonpulse train generated by the quantum dot source does not need to be spatially demultiplexed.
An experimental implementation of the multiboson correlation sampling scheme 15 with photons generated from three atomic-ensemble quantum memories has been reported 87 very recently, showing that measurement over inner modes can be a useful resource in multiphoton networks. 88

Validation of Boson Sampling
By its very nature, boson sampling 9 is a sampling task, whose solution might not be efficiently verifiable (contrarily to problems in the complexity class NP, such as factoring). It is nevertheless crucial to identify suitable methodologies able to tackle this fundamental aspect. Indeed, any experimental instance approaching the regime of quantum advantage must be necessarily accompanied by appropriate evidence that the employed device is correctly performing the sampling process, rather than drawing data from other distributions that can be efficiently classically sampled that merely resemble the correct one. The first experimental instances of boson sampling [63][64][65][66] have employed a direct comparison between the measured sample and the expected distribution, evaluated classically with a numerical approach, as a means to verify the device operation. However, such an approach cannot be scaled to larger instances since calculating the expected distribution is exponentially hard, and an exponentially large sample has to be collected in order to obtain a meaningful estimate of the output probabilities.
While full certification is believed to be out of reach, a bottom-up approach to assess the correct functioning of a boson sampling machine is currently under investigation. The idea behind such schemes is to identify suitable classically computable mock-up distributions that represent plausible error models for the device operation. Appropriate methods should then be applied to exclude that the data were sampled from such mock-up distributions. Finally, progressively refined models representing more stringent tests can be included in the toolbox.
Starting from Ref. 89, which raised the fundamental question of boson sampling verification, some initial approaches to validation were developed (a schematic view is shown in Fig. 4). The very first example of mock-up distribution is the uniform one, for which an efficient method requiring a very small sample size has been proposed 19 and tested experimentally. 67,76 The validation technique is based on the evaluation of an efficiently computable quantity, which is obtained from the elements of the unitary matrix and is only partially correlated to the permanents expressing the input-output probabilities. A less trivial mock-up distribution is the one obtained by sampling an m-mode interferometer, described by the same unitary transformation U, with distinguishable particles as inputs. In the latter case, no quantum interference occurs and an efficient classical sampling algorithm can be employed to draw data samples. 9 To discriminate between this distinguishable-photon scenario and a boson sampling machine, a method based on likelihood ratio tests has been developed and experimentally tested 76 with up to 3 photons in a 13-mode integrated interferometer, also considering a Bayesian variant 90 and its extension to the analysis of scattershot boson sampling data. 68 More recently, likelihood ratio tests have been employed in more complex boson sampling experiments [70][71][72][73] and even used to assess a method for the classical simulation of boson sampling based on Markov chains. 18 The upside of likelihood ratio tests can be found in the fast convergence rate in terms of sample size, whereas the main downside is that they rely on the evaluation of matrix permanents, and thus they are not computationally efficient. A different approach to discriminate against distinguishable particles can be found in exploiting collective properties of multiphoton interference such as bunching or clouding, 67 the latter being the tendency of bosons to exit on nearby output modes in quantum walks unitaries. By combining the observation of such property with Bayesian hypothesis testing, it was possible to discriminate multiphoton interference with respect to distinguishable particles for 3 photons in an integrated quantum walk unitary of 21 modes.
As a matter of fact, approaches employing bunching or clouding as a signature of the correct behavior of a boson sampling machine can be tricked by a classically computable configuration called mean-field sampler. 91 The latter is a model originating from independent, noninterfering single-photon states with random phases and is able to mimic some (but not all) features of genuine multiphoton interference. To help discriminate against this model, zero transmission laws 92 based on peculiar transformations can be employed. More specifically, for certain unitary matrices and input states possessing specific symmetries, it is possible to efficiently predict (without relying on permanent calculations) that a significant subset of inputoutput combinations will never be observed when genuine multiphoton interference is present. Hence, by implementing such transformations and by directly counting how many events are recorded in the forbidden combinations, it is possible to verify whether the input state is a mean-field sampler. These zerotransmission laws have been proposed 91,92 for Fourier interferometers, and tested experimentally in integrated devices by employing a universal design 69 or an efficient 3-D layout, 93 leading to the observation of zero transmission law with 2-and 3-photons in devices with up to 8 modes. Theoretical 94,95 and experimental 96 papers have extended this approach to so-called Sylvester interferometers, which were analyzed using a general mathematical framework. 97,98 Sylvester interferometers have been also suggested and tested as being optimal transformations for certain number of photons and modes, 99 in terms of minimizing the achievable error rate as a function of the sample size in a general test of genuine multiphoton interference. Further studies focusing on specific interferometer designs have shown the possibility of constructing a genuine indistinguishability witness in a similar spirit to the detection of genuine multipartite entanglement. 100 Moving a step further, other validation methods that can be adopted for any linear evolution and do not require the implementation of purpose-built interferometers have been devised. A first class of validation tests based on pattern recognition techniques 101 was developed for both boson sampling validation and its scattershot variant. Such an approach constructs a compatibility test 102 that compares two different samples and determines whether they have been both generated by indistinguishable (distinguishable) particles or not. This is obtained using machine learning techniques for the identification of internal patterns within the measured samples. This approach has been tested in 3-photon, 13-mode experiments. A different methodology relying on statistical signatures of many-body quantum interference has been recently proposed 103 (with an associated perspective article 104 ). This method is based on advanced and mostly analytic tools from statistical physics and random matrix theory and is efficient in the number of photons n and modes m. This approach has been implemented experimentally 105 in a three-photon, seven-mode experiment, by including machine learning techniques to optimize the identification whether the sample has been generated by indistinguishable or distinguishable photons. Further methods have been proposed recently, 106,107,108 but they are still lacking of experimental implementation.
The latter two methods we have discussed in some detail are actually representatives of the two main categories of validation schemes: (i) approaches based on computational criteria and (ii) physically motivated approaches. In the first case, the validator employs computational techniques to find patterns or peculiar features in data samples. Such techniques can be general-purpose and do not rely on knowledge of the system under investigation (multiphoton evolution in multimode interferometers). Those approaches directly look at the output data produced by the quantum device and search for common properties without making any specific assumption on the system under investigation. In the second case, the validation methods exploit the physical phenomenon behind the boson sampling task, and thus design techniques that rely on the knowledge that output data have been produced by multiphoton interference. These physically motivated approaches can also employ tools inspired by other physical systems, such as many-body physics.
Note that the two distinct methodologies described above (computationally or physically based) can also be combined and benefit from each other. On the one side, computationally oriented approaches can provide new insights in the physical problem, thus highlighting hidden features otherwise not detectable. On the other hand, the performances of physically based approaches can be boosted by adding computational modules in the data identification stage.

Scalability Toward Quantum Supremacy
It is important to note that many of the requirements for boson sampling may be relaxed if future theoretical work manages to provide evidence for hardness of simulation of simpler experiments. For example, the requirements of using m ¼ Oðn 2 Þ modes and that interferometers be drawn from the uniform, Haar ensemble, were assumptions made for technical convenience in the original work by AA, 9 but it is an open question whether they are strictly necessary.
One major scalability issue that seems to follow from the requirement of Haar-random interferometers is related to losses. Arbitrary linear-optical networks have worst-case depth that is linear (in the number of modes). 23,24,109 If each beam splitter has a fixed loss factor, this implies an exponentially small per-photon probability of survival in the interferometer, which would imply that boson sampling is not scalable. One important avenue of theoretical work is to investigate the hardness of boson sampling with other interferometer ensembles that hopefully can be built with shorter depth (similar to the idea of t designs in quantum computing 110 ) and thus circumvent losses. Another theoretical improvement that would help in this direction would be proving the hardness of boson sampling with fewer than Oðn 2 Þ modes [hopefully only OðnÞ modes]. A caveat is that this might introduce the need for number-resolving detectors. Interestingly, in the context of exact boson sampling it is known that interferometers of depth 4 and Oð2nÞ modes are already hard to simulate, 12 although it remains an important open question whether a similar result holds also for approximate simulation.
The other main scalability issue to reach the quantum advantage regime with photonic boson sampling devices is represented by photon sources. The main requirements toward a scalable implementation are high generation probability, high level of photon indistinguishability, and low probability of generating higher order photon number terms. Current approaches employ parametric down-conversion, spontaneous four wave mixing, and solid-state quantum-dot sources.
The first category of photon generation modules suffers from an intrinsic small generation probability. Indeed, parametric down-conversion is a probabilistic process where the generation efficiency g must be kept low to avoid significant contributions from higher order terms with more than one photon per mode. This feature becomes an issue when using parametric downconversion sources for fixed input boson sampling since the probability of generating the correct n-photon input is exponentially decreasing as g n . Nevertheless, some of the variants of the original task, such as scattershot and Gaussian boson sampling, employ a careful design leading to a significant enhancement in the generation rate with parametric down-conversion sources. Scaling up to a larger number of photons requires reaching high values of heralding efficiencies and photon indistinguishability. Recent studies have shown the possibility of tailoring the source characteristics to obtain significant improvements in those parameters. 74 Further perspectives for scaling up the number of photons can be envisaged with multiplexing, 84 within an integrated platform. 80 The second approach employing solid-state sources has been recently introduced in the boson sampling context. Such a platform presents significantly higher efficiency and generation rate 73 than parametric down-conversion, maintaining a negligible probability of generating multiple photons at the same time. Current experiments employ single quantum-dot sources combined with time-to-spatial multiplexing, 85 thus distributing a train of n single photons generated at subsequent time bins in n different spatial modes. Such an approach requires implementing multiplexing with a low level of losses for scaling up to large photon number. An alternative approach would require generating highly indistinguishable photons from different quantum-dot sources.
The problem of devising photon sources of high efficiency and indistinguishability is also strongly related to the transition between a classically simulable regime and the hard-to-simulate regime. More specifically, recent studies 38,44 have shown that partial photon indistinguishability p limits the hardness of a classical simulation. Indeed, given a certain value of p, it has been shown that an n-photon system presents a complexity that corresponds to a lower number of photons k < n. Furthermore, such threshold only depends on the photon indistinguishability p and is actually independent from n, meaning that it is not possible to overcome this issue only by increasing n. Thus strong effort should be dedicated to building sources with a high level of photon indistinguishability.

Outlook
In this review, we have discussed recent advances in the implementation of photonic boson sampling. Starting from the original proposal by AA 9 and from the first instances reported in Refs. 63-66, significant progress has been reported on both theory and experiments. Several theoretical studies have focused on different aspects of the problem, including the definition of the limits for classical simulations, the role of experimental imperfections in the complexity of classical simulation, as well as proposals of new variants of the original problem. In parallel, technological advances have led to the implementation of progressively larger experimental instances. Indeed, tools such as quantum dot single-photon sources or superconducting nanowire single-photon detector have enabled the possibility to generate and measure samples with larger number of photons, opening the route to further advances. Additionally, progress in the implementation of integrated multimode interferometers allowed the inclusion of reconfigurable elements in the structures, in order to realize arbitrary linear operations.
A fundamental issue toward achieving the regime of quantum advantage is the capability to identify whether a given set of data has been sampled from the correct distribution. This is particularly important for future, larger-scale implementations where the classical simulation will become unfeasible. The first results focusing on identifying specific errors (such as sampling from a distribution with distinguishable particles) via different techniques have been developed and tested, but additional work in this direction is required.
Finally, the first proposals of application of boson sampling devices outside its original scope have been reported. For instance, it has been shown that Gaussian boson sampling represents a tool for photonic simulation of vibronic spectra in molecular dynamics. 59,60,75 Promising perspectives have been also pointed out in the context of hybrid quantum computation: in this approach, specific quantum routines are embedded in a more complex computation task performed with a classical approach, providing an overall quantum boost. Some examples have very recently been proposed for the Gaussian boson sampling variant, which is related to the sampling of hafnians of the evolution matrix. More specifically, two relevant scenarios have been identified. On one side, Arrazola et al. 58 have introduced an NP-hard optimization problem (Max-Haf). A Gaussian boson sampling device can be inserted within the computation performed by stochastic algorithms (such as random search or simulated annealing) to improve the algorithm performance.
The key ingredient can be found in replacing the uniform sampling part of the algorithm with a quantum device that performs sampling from a probability distribution proportional to matrix hafnians, thus leading to a faster convergence to the output of the computation. In parallel, the connection between hafnians and graph theory 111 has led to a few proposals aiming to exploit Gaussian boson sampling in the analysis of complex graphs. For instance, in Ref. 57, it has been shown that the estimation of the number of perfect matchings in undirected graphs can be performed more efficiently by encoding the adjacency matrix of a graph in a Gaussian state, allowing one to improve the sampling success rate. Furthermore, stochastic algorithms that tackle an NP-hard problem, such as the densest k-subgraph problem, can have its performance enhanced when sampling is performed via a Gaussian boson sampling device. 56 Finally, a very recent work 112 has shown the connection between such computational model and graph isomorphism. These ideas can represent a starting point for the application of quantum boson sampling devices also to problems not envisaged in the original theoretical proposal. Fabio Sciarrino received his PhD in 2004 with a thesis in experimental quantum optics. He is a full professor and head of the Quantum Information Lab in the Department of Physics of Sapienza Università di Roma. Since 2013, he has been a fellow of the Sapienza School for Advanced Studies. His main field of research is quantum information and quantum optics, with works on quantum teleportation, optimal quantum machines, fundamental tests, quantum communication, and orbital angular momentum.