Boson Sampling with Gaussian input states: toward efficient scaling and certification

TL;DR

This paper proposes a practical approach to scale Boson Sampling experiments by integrating continuous-variable quantum information, temporal encoding, and advanced detection techniques, aiming to demonstrate quantum advantage more efficiently.

Contribution

It introduces a novel combination of methods including temporal loops and scattershot Boson Sampling to enhance scalability and certification in Boson Sampling experiments.

Findings

Proposes a scalable experimental setup for Boson Sampling.

Details assumptions for computational hardness in the new configuration.

Moves towards efficient certification of quantum advantage.

Abstract

A universal quantum computer of large scale is not available yet, however, intermediate models of quantum computation would still permit demonstrations of a quantum computational advantage over classical computing and could challenge the Extended Church-Turing Thesis. One of these models based on single photons interacting via linear optics is called Boson Sampling. Although Boson Sampling was demonstrated and the threshold to claim quantum computational advantage was achieved, the question of how to scale up Boson Sampling experiments remains. To make progress with this problem, here we present a practically achievable pathway to scale Boson Sampling experiments by combining continuous-variable quantum information and temporal encoding. We propose the combination of switchable dual-homodyne and single-photon detections, the temporal loop technique, and scattershot-based Boson Sampling. We detail the required assumptions for concluding computational hardness for this configuration. Furthermore, this particular combination of techniques moves towards an efficient scaling and certification of Boson Sampling, all in a single experimental setup.

Equations

$$P(n)=\underset{\text{scattershot Boson Sampling}}{\underbrace{\left(\genfrac{}{}{0ex}{}{m}{n}\right)\underset{\text{standard Boson Sampling}}{\underbrace{{\chi }^{2n}(1-{\chi }^2{)}^m}}}},$$

$$\chi =\sqrt{\frac{n}{m+n}}$$

$$s_i=(x_{1,i},p_{1,i},x_{2,i},p_{2,i},\dots ,x_{m,i},p_{m,i}{)}^T$$

$${\xi }_i=\frac{1}{C^2}{s}_i{s}_i^T,$$

$$⟨{\xi }_i⟩=\frac{1}{2C^2}(\Sigma +I)$$

$$\bar{\xi }=\frac{1}{K}\sum _{i=1}^K{\xi }_i$$

$$\frac{1}{2C^2}(\Sigma +I)\ge \frac{1}{2C^2}I$$

$$\displaystyle\text{Pr}\Big{\{}\frac{(1-\eta)}{2C^{2}}(\Sigma+I)\leq\bar{\xi}\leq\frac{(1+\eta)}{2C^{2}}(\Sigma+I)\Big{\}}$$

$$\ge 1-8m{e}^{-K{\eta }^2/(8C^2\ln 2)}.$$

$$\Sigma \ge bI.$$

$$\displaystyle\text{Pr}\Big{\{}\frac{(1-\eta)\Sigma-\eta I}{C^{2}}\leq 2\bar{\xi}-\frac{I}{C^{2}}\leq\frac{(1+\eta)\Sigma+\eta I}{C^{2}}\Big{\}}$$

$$\ge 1-8m{e}^{-K{\eta }^2/(8C^2\ln 2)}.$$

$$\displaystyle\text{Pr}\Big{\{}\frac{[1-\eta(1+b^{-1})]\Sigma}{C^{2}}\leq 2\bar{\xi}-\frac{I}{C^{2}}\leq\frac{[1+\eta(1+b^{-1})]\Sigma}{C^{2}}\Big{\}}$$

$$\ge 1-8m{e}^{-K{\eta }^2/(8C^2\ln 2)}.$$

$$\displaystyle\text{Pr}\Big{\{}\frac{(1-\eta)}{2C^{2}}(\Sigma+I)\leq\bar{\xi}\leq\frac{(1+\eta)}{2C^{2}}(\Sigma+I)\Big{\}}$$

$$\ge 1-8m{e}^{-K{\eta }^2/(8C^2\ln 2)}$$

$$\displaystyle\text{Pr}\Big{\{}\frac{(1-\eta)\Sigma-\eta I}{C^{2}}\leq 2\bar{\xi}-\frac{I}{C^{2}}\leq\frac{(1+\eta)\Sigma+\eta I}{C^{2}}\Big{\}}$$

$$\ge 1-8m{e}^{-K{\eta }^2/(8C^2\ln 2)}$$

$$\displaystyle\text{Pr}\Big{\{}\frac{[1-\eta(1+b^{-1})]\Sigma}{C^{2}}\leq 2\bar{\xi}-\frac{I}{C^{2}}\leq\frac{[1+\eta(1+b^{-1})]\Sigma}{C^{2}}\Big{\}}$$

$$\ge 1-8m{e}^{-K{\eta }^2/(8C^2\ln 2)}$$

$$s_i\sim {\mathcal{N}}_{4m}\left(0,(\Sigma +I)/2\right)$$

$${\xi }_i=\left(\frac{s_i}{C}\right)\left(\frac{s_i^T}{C}\right)=\frac{1}{C^2}{s}_i{s}_i^T$$

$$\bar{\xi }=\frac{1}{K}\sum _{i=1}^K{\xi }_i$$

$$\bar{z}=[{\bar{\xi }}_1,{\bar{\xi }}_2,\cdots ,{\bar{\xi }}_{\text{runs}}]$$

$$\frac{[1-\eta (1+b^{-1})]\Sigma }{C^2}\le 2\bar{\xi }-\frac{I}{C^2}$$

$$⟺2\bar{\xi }-\frac{I}{C^2}-\frac{[1-\eta (1+b^{-1})]\Sigma }{C^2}\ge 0$$

$$2\bar{\xi }-\frac{I}{C^2}\le \frac{[1+\eta (1+b^{-1})]\Sigma }{C^2}$$

$$⟺\frac{[1+\eta (1+b^{-1})]\Sigma }{C^2}-2\bar{\xi }+\frac{I}{C^2}\ge 0.$$

$$\text{Pr}\Big{\{}\frac{[1-\eta(1+b^{-1})]\Sigma}{C^{2}}\leq 2\bar{\xi}-\frac{I}{C^{2}}\leq\frac{[1+\eta(1+b^{-1})]\Sigma}{C^{2}}\Big{\}}=\frac{\text{\# of True cases}(\bar{z};\eta)}{\text{runs}}$$

$$1-8m{e}^{-K{\eta }^2/(8C^2\ln 2)}$$

$$\forall i\in [K]:0\le {\xi }_i\le I.$$

$$\Sigma =O^T{\Sigma }_{\text{diag}}O,$$

Full text

I Introduction

Boson Sampling is a model of intermediate—as opposed to universal—quantum computation, initially proposed to confront the limits of classical computation compared to quantum computation AA . An efficient classical computation of the Boson Sampling protocol would support the Extended Church-Turing Thesis “which asserts that classical computers can simulate any physical process with polynomial overhead”harrow17 , i.e., polynomial time and memory requirements. But an efficient classical algorithm for Boson Sampling would also imply that the Polynomial Hierarchy of complexity classes, which is believed to have an infinite number of distinct levels, would reduce (or “collapse”) to just three levels. Consequently, a computer scientist could not simultaneously support the Extended Church-Turing Thesis and an infinite structure of the Polynomial Hierarchy. Hence, one is cornered into a position that either a fundamental change in our understanding of computational complexity is needed or quantum-enabled algorithms must be able to perform some tasks efficiently that cannot be performed efficiently on a classical computer.

Even in the absence of a full-scale quantum computer, a physically constructed Boson Sampling device could outperform a classical device sampling from the same distribution harrow17 ; lund2017 ; review_road_QCS ; brod_2019_review . The first experimental demonstration of quantum computational advantage was published in October 2019 for a sampling problem using 53 qubits in a superconducting circuit architecture google_q_supremacy . Recently, quantum computational advantage was demonstrated by two other groups: one using photons with 50 single-mode squeezed states in a 100-mode linear optical network q_advantage_2020_jwpan , and the other one using 216 photon and temporal encoding in a 3D topology Madsen2022 .

In a simplified view, an implementation of the Boson Sampling protocol can be summarized as follows: $n$ indistinguishable single photons are inputs into the ports of a $m$ modes linear optical network, represented by a unitary matrix ($U$), and at the output ports single-photon detection is performed (Fig. 1). Any alleged Boson Sampling device must give samples from this output distribution for any given $U$. The computational hardness is expressed as the scaling of runtime in terms of the number of input photons $n$. The classical hardness argument has been given when $m$ scales as $n^{5}\log n$, but it likely also holds for $m$ scaling as $n^{2}$.

Implementations of Boson Sampling have been successfully demonstrated, initially using single-photon pairs from Spontaneous Parametric Down Conversion (SPDC) and later using single photons from Quantum Dots BroomeBS ; SpringOX-BS ; tillmannVienna ; CrespiRome ; ScattershotRome ; LoredoBS ; JWP_BS_2017 ; JWPan_BS_2018_PRL ; flamini2018photonic . A more advanced implementation achieved the threshold of quantum computational advantage using 50 single-mode squeezed states as inputs for an optical network with 100 modes q_advantage_2020_jwpan . Another important milestone was the photonic implementation of programmable Boson sampling with 216 squeezed modes entangled with three-dimensional connectivity using a time-multiplexed and photon-resolving detection Madsen2022 . This work showed a runtime advantage over 50 million times as extreme as reported by earlier photonic machines.

The question is then how to keep scaling up Boson Sampling experiments. Three factors are currently contributing against quantum demonstrations of Boson Sampling: (a) better classical algorithms which move the threshold of quantum computational advantage to greater number of input single photons nevilleBristol2017 ; Clifford&Clifford ; Clifford&Clifford_2020 , e.g. the classical algorithm of Neville et al. nevilleBristol2017 solved the Boson Sampling problem with 30 photons in a standard computer efficiently; (b) difficulties on the scaling of the preparation of multiple single photons (i.e $\ket{1}^{\otimes n}$); and (c) scaling of photon losses in the linear optical network Patron_Renema ; Oszmaniec_Brod_err ; Renema_classical_simulability_arxiv . Different theoretical directions have been pursued to study different types of quantum states Lund2014 ; gaussian_BS_prl_2017 ; Chabaud2017 and measurements Lund2017measurement ; Chakhmakhchyan2017 to circumvent these problems. However, the challenge of scaling preparation and certification remains a topic of interest.

Implicit in the context of the implementation of Boson Sampling devices is that the operation of the device does indeed satisfy the theoretical requirements for an advantage over classical computational devices. The condition presented in Ref. AA was on the total variation distance for the output distribution compared with the ideal distribution. The standard view is aligned with that of cascaded quantum networks (e.g. constructions like those of Gardiner1993 ; Carmichael1993 ) where the performance of individual components are characterised and their performance as a network is extrapolated. This approach, though conceptually straightforward, suffers from exponential sensitivity on the performance of interconnects between components. Ideally, the performance of a Boson Sampling device could be characterised in situ. A number of efficient schemes have been proposed to perform these kinds of tests, such as tests against uniform distributions Aaronson2014 , particle distinguishability Shchesnovich2016 , mean-field approximations Tichy2014 , and others (see Ref. Brod2019_review for a detailed list of such tests). The drawback of all these tests to date is that they test a particular property of the output distributions, rather than verifying that they confirm the output distributions meet the theoretical requirements for a quantum advantage with Boson Sampling.

Due to the exponential sized sampling space at the output, it is a difficult problem in general to find many individual conditions required to be simultaneously satisfied by a non-ideal output distribution. For the exponentially sized sample space, an exponentially sized ensemble is required to directly estimate the total variation distance without any prior assumptions (see Hangleiter2019 for properties of distributions needed to avoid this issue).

At this point in time, Boson Sampling faces difficult challenges involved with building and characterising scaled devices. Motivated by these constraints, we propose a method to scale Boson Sampling experiments using continuous-variable quantum states and time-bin encoding. Our analysis of this proposal takes into account finite squeezing and given some reasonable assumptions hold, operational performance can be characterized efficiently. Our scheme presents an efficient way to implement Boson Sampling experiments, and also an efficient way to certify it using the same experimental setup. Here, by changing the detection type and making the assumption that the state before detection is Gaussian, we certify the total variation distance via the state covariance matrix. This is a quantity that has a polynomial dependence on the input size, and hence this statistical distance can be efficiently estimated.

It is noteworthy that the current search for applications of Boson Sampling goes beyond the scope of computational complexity. For instance, Boson Sampling has been adapted to simulate molecular vibrational spectra molecularBS_theory ; molecularBS_exp and may be used as a tool for quantum simulation aspuru2012photonic ; georgescu2014quantum . Other Boson Sampling-inspired applications are the verification of NP-complete problems verifying_NP , quantum metrology sensitivity improvements MORDOR , and a quantum cryptography protocol BS_quantum_crypt .

II Proposal

Here we propose a technique to scale Boson Sampling experiments based on continuous variables (CV) and temporal encoding. In the CV case, the information is encoded on the quantum modes of light, specifically, on the eigenstates of operators with continuum spectrum kok2010introduction . Continuous-variables quantum information has achieved impressive results. An initial report of 10,000 entangled modes in a continuous-variable cluster mode Furusawa_10000 was later upgraded to one million modes Furusawa_one_million . Some of these systems were conceived to perform measurement-based quantum computation, and here we show they can be adapted to Boson Sampling. Moreover, while some of the theoretical work for measurement-based quantum computation assume unrealistically infinite squeezing, here we require only finite squeezing. The world record for detected squeezed light is 15 db Roman_15db , while 20.5 db of squeezing would give rise to conditions for fault-tolerant quantum computation under the Gottesman-Kitaev-Preskill encoding Menicicci2014 ; GKP_encoding .

The work of Lund et al. Lund2014 (a.k.a. Scattershot Boson Sampling) demonstrated Gaussian states can be used as inputs in Boson Sampling experiments and only bounded squeezing is necessary, provided each output is projected in the number basis by single-photon detection. It is important at this point to emphasise that the specific task of Scattershot Boson Sampling requires the generation of a different distribution to standard Boson Sampling. One result presented in Ref. Lund2014 was to show that sampling from this alternative distribution is also a hard problem when employing classical computing resources. The single photon Boson Sampling distribution is contained within the Scattershot Boson Sampling distribution, but this is merely used as part of the proof for computational hardness. The requisite task is the efficient generation of samples from a Gaussian state measured in the Fock basis without any further processing.

A detailed discussion of how much squeezing is necessary for Scattershot Boson Sampling experiments can be found on Ref. Lund2014 . Interestingly, the authors Lund2014 showed that for a two-mode squeezer, like SPDC, there is a trade-off between the strength of the SPDC (linked to $\chi$) and the most likely number of photons detected represented by the variable $n$. This indicates that Scattershot Boson Sampling experiments that are performed with less photons would require higher $\chi$ levels:

[TABLE]

where $P(n)$ is the probability of detecting $n$ photons and this probability is locally maximised when

[TABLE]

and this maximum probability is lower bounded by $1/\sqrt{n}$ if $m\geq n^{2}$ and $n>1$. In this regime, $\chi$ decreases as $n$ increases and when taking $m=n^{2}$ at $n=8$ only $3$dB of squeezing is required to achieve this optimal probability. See 111In the case of Scattershot Boson Sampling Lund2014 , $m$ refers to the number of two-mode squeezers, i.e., the number of SPDC single-photon sources. In that paper, the condition $m=n^{2}$ is imposed. In the present work, we use a more general definition: $m$ is simply the number of modes in the linear optical network represented by a $m$ x $m$ unitary matrix ($U$). for clarification of the notation on the Scattershot Boson Sampling.

More recently, a new variant of Boson Sampling was proposed in which the Boson Sampling protocol is formulated in terms of the Hafnian of a matrix, a more general function belonging to the #P class, and was proven for the case of exact sampling. Moreover, it was shown in Ref. Grier2022 that under the standard Anti-Concentration and Permanent-of-Gaussians conjectures, there is no classical algorithm to efficiently sample from the ideal (Hafnian-based) distributions, even approximately, unless the polynomial hierarchy collapses. This protocol is called Gaussian Boson Sampling gaussian_BS_prl_2017 ; gaussian_BS_arxiv_supp . In this present work, we use the terminology scattershot Boson Sampling and Gaussian input Boson Sampling as synonyms. A demonstration of scattershot, Gaussian, and standard Boson Sampling using integrated optics and single photons from Spontaneous Four Wave Mixing sources can be found in Reference BS_integrated_bristol_2019 .

II.1 Scaling

Now, we will analyse the scaling of our Boson Sampling experimental proposal. Consider two pulsed-squeezed-light sources, with $\tau$ time interval between subsequent pulses, where these two states interfere in a 50:50 beamsplitter, followed by a controllable delay, where a pulse can be delayed by $N\tau$ before being released, Fig.(2). This may be a loop architecture Time-Bin_BS_motes or a quantum memory, provided it returns the given input state with enough high fidelity. Here we will use the loop architecture. The modulator should implement the desired unitary $U$ operation by interfering delayed pulses. At the end of each spatial path, there are two possible measurement schemes that can be performed. Either, the light can be sent to a single-photon detector to record the samples (output) for Boson Sampling, or the light is directed towards a homodyne detection leonhardt1997measuring ; Bachor&Ralph setup that is used to characterize the output state, including the output state from the optical network.

The use of time-bin (loop architecture) for Boson Samplings has been proposed and experimentally demonstrated Time-Bin_BS_motes ; JWP_BS_2017 . This technique converts spatial modes into temporal modes. This gives a quadratic reduction in required spatial modes and provides a significant benefit in scaling to a large number of single photons as input states for the Boson Sampling protocol. However, the original time-bin proposal was conceived for the discrete case, i.e., manifold single photons. The certification in the discrete case, as addressed in the article Aolita2015 requires the tested state to be generated $\mathcal{O}(nm^{n})$ times when verifying a process with $n$ photons in $m$ modes. So, for $m=n^{2}$, this is worse than an exponential growth. Therefore, the discrete case does not constitute a scalable approach for Boson Sampling when also considering certification. The certification of the Boson Sampling protocol is an important aspect of our proposal and is discussed in section Certification below.

A significant benefit to our approach is that, under some reasonable assumptions, the operation of the sampling device can be characterized using the sampling state itself without the need for other probe input states. To achieve this, the following assumptions are needed: (i) the output state received by the single-photon detection is the same as that received by the homodyne detection, which is achievable by movable mirrors, for example as in the procedure given by webb2006homodyne ; (ii) the two squeezed input states are Gaussian and that the modulation network changes the states but leaves the output still in a Gaussian form, a standard Gaussian optics property; (iii) the output is fully characterized by a multi-mode covariance matrix, and finally (iv) the choice of when to make a sampling run and when to make a characterization run is irrelevant. In other words, the experimental setup is assumed stable, and the output will not change over the time as one changes between the two different measurement schemes.

A Gaussian output state can be fully characterized by the mean vector (which we will assume zero) and covariance matrix. For a $m$ mode state and $n$ detected photons, the number of possible photon number detection events scales as $m^{n}$. However, to describe a Gaussian state before the detection has occurred, only the number of entries in a covariance matrix for a $m$ mode state is required and this scales as $4m^{2}$. For the case of Gaussian input Boson Sampling (a.k.a. scattershot) where there are two groups of $m$ modes and $n$ photon detections, the size of the Fock basis detection sample space is $m^{2n}$, but the full covariance matrix for the state prior to detection will require $16m^{2}$ entries.

Performing the characterization involves reconstructing the covariance matrix from the CV measurement samples. The measurements chosen must be sufficient in number to estimate all elements of the covariance matrix, including terms involving the correlations between X and P in the same mode. To avoid repeated changes to measurement settings, we propose performing this by means of dual homodyne. In a dual-homodyne arrangement, the signal mode is split at a 50:50 beamsplitter and both modes undergo a CV homodyne detection, one measured in X and the other in P. This permits a simultaneous measurement of the X and P quadratures at the cost of adding $1/2$ a unit of vacuum noise to the diagonal elements of the state covariance matrix. So, if $\Sigma$ is the state covariance matrix, then the dual homodyne modes will see Gaussian statistics with a covariance matrix of $(\Sigma+I)/2$ (under units where the variance of vacuum noise is unity), where $I$ is the identity operator. This covariance matrix can then be estimated by constructing matrix-valued samples from each sampling run. Let

[TABLE]

be a $2m$-dimensional real vector representing the $i^{\mathrm{th}}$ data sample from the dual-homodyne measurement, with the first subscript representing the mode to which the corresponding homodyne detector is attached. From this sample vector, a sample matrix can be formed from the outer product of the $s_{i}$

[TABLE]

where $C$ is a scaling factor. The reason for introducing this scaling factor will become apparent shortly. For now, one can think of it as being related to the maximum energy of the state. This sample matrix is a positive semi-definite matrix for all $i$. The expectation value for each sample $\xi_{i}$ over the incoming Gaussian distribution is then

[TABLE]

and so a sample average over $K$ samples

[TABLE]

will be an unbiased estimator for $(\Sigma+I)/(2C^{2})$.

To see how close the sample average is to the true average, we apply the operator Chernoff bound wilde_2013 ; nielsen_chuang_book (following the notation of Wilde wilde_2013 , Section 16.3). This provides a lower bound for the probability that the sample average does not deviate significantly from the expected value. Let $K$ be the number of sample matrices and $\bar{\xi}$ the sample average of $K$ samples as defined in Eq. 6. The input state covariance matrix $\Sigma$ is positive definite, and we have

[TABLE]

which is the expectation of each operator forming the sum in Eq. 6. For this situation, the operator Chernoff bound for any $0<\eta<1/2$ is given by

[TABLE]

To then bound the probability for making a multiplicative estimate of $\Sigma$, the spectrum of $\Sigma$ needs to be bounded away from zero.

Let the parameter $b$ represent the variance of the quadrature for the maximum possible squeezing for the state being estimated. This means that

[TABLE]

The Chernoff bound for the estimator $\bar{\xi}$ can be rewritten as

[TABLE]

Then using the inequality in Eq. 9, this can be written as

[TABLE]

This means the rewritten estimate $2\bar{\xi}-I/C^{2}$ gives a multiplicative estimate of the scaled covariance matrix $\Sigma/C^{2}$.

The interpretation of Equation 11 is that the chance that the finite sample estimate of the covariance matrix deviates from the true value decays exponentially in the number of covariance matrix samples $K$ and the square of deviation permitted $\eta$, but depends linearly on $m$, the number of modes. In our application of Gaussian input Boson Sampling, the value of $b$ is fixed, as too much squeezing can actually degrade performance. So for fixed $\eta$, as the number of modes $m$ increases, the number of samples $K$ required to achieve the same probability bound in the operator Chernoff bound only grows logarithmically $\mathcal{O}(\ln m)$.

The validity of the operator Chernoff bound applied to the characterization of Gaussian states can be demonstrated numerically. We draw samples from a multivariate Gaussian distribution with mean-vector zero and covariance matrix $(\Sigma+I)/2$ (and scale the vectors by $C$). The vectors are used to construct the matrix samples $\xi_{i}$. The probability on the left-hand side of Eq. 11 is found by computing the sample average and finding how often the two conditions (in curly brackets) are satisfied. More details can be found in the Supplemental Material.

Numerical simulations were performed for the cases of $m=1,4,9$ for several choices of $\Sigma$, $K$, and $b$. All the cases that were examined showed an agreement with Eq. 11. In Fig. 3, we present some of the numerical results for the specified parameters. The figures presented here consider only the case of diagonal covariance matrices. These cases show the essential features of the bound as the eignvalues of the covariance matrix contain the most critical information with respect to the Chernoff bound. Nevertheless, we present numerical simulations of non-diagonal covariance matrices, confirming this heuristic, in the Supplemental Material.

Note that a prerequisite for applying the operator Chernoff bound is that all samples $\xi_{i}$ have all of their eigenvalues in $[0,1]$ wilde_2013 . Although a smaller $C$ value (with other parameters being fixed) can provide a tighter bound, this prerequisite for applying the bound would be satisfied at an increasingly infrequent rate. For the purposes of these numerical simulations, the value of $C$ was chosen large enough such that this prerequisite is satisfied at a reasonable rate.

II.2 Certification

Finally, one would like to certify if the generated state is sufficient to perform the task at hand, that is Boson Sampling. For approximate Boson Sampling, one does not need to generate the state ideally but within some trace distance bound $\varepsilon$. Using the Fuchs-van de Graff inequality fuchs1997cryptographic , a trace distance is upper bounded by the fidelity by $\varepsilon<1-F$. A robust certification strategy is given by Aolita et al. Aolita2015 , which tests if a fidelity lower bound (or equivalently maximum trace distance) holds between a pure Gaussian target state and a potentially mixed preparation state. Here, we argue that this result can also be employed for certification. In order to perform the certification, the Gaussian covariance matrix elements need to be estimated and manipulated with knowledge of the target pure state. This produces a bound of the fidelity, which can be used to test for appropriateness of the apparatus to perform Gaussian input Boson Sampling. The samples needed to achieve a fixed fidelity bound (or fixed trace distance) is higher than the Chernoff bound and scales as $\mathcal{O}(m^{4})$ times, where $m$ is the number of modes in the state being certified. This certification process can require considerable amounts of data to scale, but the scaling with system size is polynomial, making the process more likely to be feasible. This is as opposed to the certification in the discrete variables model as addressed in the same article Aolita2015 which requires the tested state to be generated $\mathcal{O}(nm^{n})$ samples when certifying a process with $n$ photons in $m$ modes.

II.3 Sampling

After the stage of certification is finished, the movable mirrors must be removed and direct the light toward the single-photon detectors. Doing so, one is projecting the Gaussian states into a Fock basis, and thus obtaining the output of the Boson Sampling experiment, similarly as in the Gaussian input Boson Sampling (a.k.a. scattershot) Lund2014 . Our proposed method greatly simplifies the numbers of required resources for scaling Boson Sampling experiments. Here we benefit from having well certified states, with certification growing polynomially as discussed above, and from having only two squeezed light sources, and thus simplifying the preparation of input states. Our method also requires fewer detectors. For instance, if one wishes to implement a 20 input single photons Boson Sampling, then it requires a $400\times 400$ linear optical network, and therefore 400 single-photon detectors. Note that not obeying $m\gg n$, at least $m=n^{2}$, violates the mathematical assumptions upon which the approximate Boson Sampling problem is currently formulated. In our proposal, due to the time-bin implementation of the linear optical network, only 2 single-photon detectors are required.

II.4 The role of imperfections

An important consideration in the performance of any sampling device is the role that imperfections, originating from any process, have on the ability one has to make conclusions about the classical easiness or hardness of computing random samples. Crucially, the dominant imperfection using current technology is photon losses. The effects of photon losses are included within the approximate sampling requirements, i.e., how much one can deviate from the perfect sampling. Unfortunately, this does not mean that losses can be neglected, as they will in many cases give rise to exponential scaling in the total variation distance between the lossless distribution and lossy distribution. For example, a constant loss rate for each mode will induce an exponential scaling in the total variation distance as a function of the number of photons to be detected. It would seem that this implies that any level of loss would render classically hard Boson Sampling impossible, but this is not the case. The goal to demonstrate the quantum computational advantage only requires producing samples that are close to some distribution, given some $\varepsilon$ tolerance, that is hard for a classical computer to reproduce. Aaronson and Brod Aaronson_Brod_error showed that the distribution generated from losing a fixed number of photons (i.e. not scaling with the number of photons) is hard for a classical computer to sample. A lossy Boson Sampling device will be close to this distribution at some scale. Unfortunately, the total variation distance for constant loss per mode will still asymptotically scale exponentially against this distribution with a fixed number of lost photons. Oszmaniec and Brod Oszmaniec_Brod_err showed that if the total number of photons that remained after loss scaled as $\sqrt{n}$ of the number input photons $n$, then a simple distribution from sampling distinguishable bosons would satisfy the total variation distance requirement for approximate sampling. Hence an efficient classical computation could reproduce samples close to the required Boson sampling distribution, nullifying the quantum computational advantage. In a later article, Brod and Oszmaniec studied the case of nonuniform losses. brod_oszmaniec_nonuniform_losses .

Despite advances in understanding the effect of losses Arkhipov_error_analysis ; rohde_ralph_error ; Aaronson_Brod_error ; L_Patron_error ; rahimi_ralph_caves_sufficient ; Oszmaniec_Brod_err ; Patron_Renema ; Renema_classical_simulability_arxiv ; brod_oszmaniec_nonuniform_losses , there remains a gap between the necessary and sufficient criteria for the hardness of approximate sampling with lossy Boson Sampling devices, which makes this an important open question. All implementations will be subject to imperfections, including the proposed implementation presented here. However, in this present work, we are primarily concerned with the experimental implications of scaling and certifying continuous-variables Boson Sampling. We expect that future results on the hardness of lossy Boson Sampling would be able to be incorporated into the proposal we present.

III Discussion

Continuous-variables (CV) quantum information, particularly in the context of optical Gaussian states weedbrook2012gaussian , has been put forward as an alternative for quantum computation. Due to the scaling aspects discussed in the previous sections, we point out that Boson Sampling can greatly benefit from the current optical CV technology Furusawa_10000 ; Furusawa_one_million ; takeda2017universal ; spie_travelling_wave_of_light ; Roman_15db ; furusawa_review_apl_2019 . In this sense, all the building blocks for this proposal have been successfully demonstrated and have good performance for scaling. The threshold for quantum computational advantage in the CV regime is currently uncertain and the subject of further investigation, but our Boson Sampling proposal certainly does not suffer from the scaling issues of discrete Boson Sampling, in particular the issue of preparing a large number of indistinguishable single photons.

In summary, we revisited the motivation behind Boson Sampling and the experimental challenges currently faced. Photon losses and scaling of many input single photons are factors working against quantum implementations of Boson Sampling. These facts pose great challenges and a new scalable approach is desirable. Here we presented a new method to do so based on continuous variables and temporal encoding. Our method assumes finite squeezing and also provides a feasible way to perform the characterization of the input states and the certification of the Boson Sampling protocol, providing viable scaling as the system size increases. With this approach, we aim to propose a pathway to scale Boson Sampling experiments, and thus, make this intermediate mode of quantum computing more compelling for its foundational questions in Computer Science as well as its applications molecularBS_theory ; molecularBS_exp ; verifying_NP ; MORDOR ; BS_quantum_crypt .

IV ACKNOWLEDGMENTS

We thank Andrew G. White, T. C. Ralph, T. Weinhold, Nathan Walk, and Daniel Brod for helpful discussions. This work was supported by the U.S. Department of Energy QuantISED award, the Brookhaven National Laboratory LDRD grant 22-22. APL acknowledges support from BMBF (QPIC) and the Einstein Research Unit on Quantum Devices. This work was supported by the Australian Research Council Centre of Excellence for Engineered Quantum Systems (Grants No. CE170100009) and the Australian Research Council Centre of Excellence for Quantum Computation and Communication Technology (Grant No. CE170100012). This work was also supported by the CFREF award on Transformative Quantum Technologies and the NSERC of Canada.

Supplemental Material

Output State Characterization Scaling: Numerical Simulations

A Gaussian state can be fully characterized by the mean vector (which we will assume zero) and the covariance matrix. Therefore, it suffices to estimate the $4m\times 4m$ covariance matrix of the Gaussian state prior to single-photon detection.

This can be done using CV homodyne detection. To estimate the scaling with the number of dual-homodyne measurement samples needed for the sample average to be “close” to the true average, one can apply the operator Chernoff bound. It provides a lower bound for the probability that the sample average does not deviate significantly from the true average.

In the main text, we present three versions of the operator Chernoff bound.

[TABLE]

In Eq. 12, the operator Chernoff bound is written for $\frac{1}{2C^{2}}(\Sigma+I)$ (includes half a unit of vacuum noise). Equations 13 and 14 describe the operator Chernoff bound for the scaled covariance matrix $\Sigma/C^{2}$.

One can numerically demonstrate the validity of these inequalities. To remind ourselves, $K$ is the number of samples, $m$ is the number of modes, $b$ is the variance of the quadrature for the maximum possible squeezing, and $C$ is a scaling factor for the CV vector samples. The steps used in the numerical simulation are as follows.

We draw $K$ $4m$-dimensional vector samples from a multivariate Gaussian distribution with mean-vector zero and covariance matrix $(\Sigma+I)/2$.

[TABLE] 2. 2.

We calculate the matrix samples from the outerproduct of each scaled vector sample, $s_{i}/C$, with itself.

[TABLE] 3. 3.

We check that each $\xi_{i}$ has all of its eigenvalues in $[0,1]$. This condition is necessary for applying the operator Chernoff bound. If this fails, we restart from step 1. A larger scaling factor $C$ gives a higher success rate. However, a smaller $C$ value provides a tighter bound. 4. 4.

We calculate the sample average

[TABLE] 5. 5.

We perform multiple runs of steps 1-4 (typically 50 runs). We then obtain a list $\bar{z}$ of all the sample averages:

[TABLE] 6. 6.

For each $\eta\in(0,0.5)$ (step size $\Delta\eta=0.001$), we use each $\bar{\xi}_{j}\in\bar{z}$ to see if the two conditions of the operator Chernoff bound under consideration are true. For example, for Eq. 14, we would check that

[TABLE]

and

[TABLE]

This can be done by verifying whether all of the eigenvalues of the matrix on the left-hand side of Eq. 15, and all of the eigenvalues of the matrix on the left-hand side of Eq. 16 are non-negative. As such, we obtain a corresponding boolean list for $\bar{z}$. 7. 7.

The probability on the left-hand side of the operator Chernoff bound under consideration is calculated for each $\eta$ as the number of true cases divided by the total number of runs. For instance, for Eq. 14:

[TABLE] 8. 8.

The right-hand side of the operator Chernoff bound under consideration is calculated (using the analytic expression) for all $\eta\in(0,0.5)$ as

[TABLE] 9. 9.

Both sides of the operator Chernoff bound under consideration are plotted as a function of $\eta$.

Note that it is important for the scaling factor $C$ to be chosen large enough such that each $\xi_{i}$ has all of its eigenvalues between zero and one:

[TABLE]

This is a necessary condition for applying the operator Chernoff bound. One can check if the chosen scaling parameter $C$ ensures that the condition above holds with high enough probability. This can be verified by creating a histogram of the top eigenvalues of $\xi_{i}$. Although a smaller scaling $C$ (with other parameters being fixed) can provide a tighter bound, the required condition above for applying the bound would be satisfied at an increasingly infrequent rate. This can be estimated as the joint probability for all of the $K$ samples $\xi_{i}$ to have their eigenvalues in $[0,1]$.

In addition to the results in the main paper, we also performed this simulation for non-diagonal covariance matrices of the form

[TABLE]

where $O$ is an arbitrary orthogonal matrix and $\Sigma_{\text{diag}}$ is the diagonal covariance matrix originally chosen. This preserves the eigenvalues of the original covariance matrix $\Sigma_{\text{diag}}$.

In Fig. 4 and Fig. 5, we present plots of the probabilities in Eq. 14 for the diagonal and the non-diagonal cases respectively. In both cases, we see agreement with the operator Chernoff bound Eq. 14.