Generation and sampling of quantum states of light in a silicon chip

TL;DR

This paper demonstrates the on-chip generation and processing of up to eight-photon quantum states in silicon photonics, implementing various boson sampling protocols and benchmarking a molecular vibronic spectra algorithm, advancing scalable quantum photonic computing.

Contribution

It introduces integrated silicon photonic circuits capable of generating and processing multiple quantum states with different sampling protocols, enabling scalable quantum algorithms.

Findings

Achieved on-chip generation of up to eight-photon quantum states.

Implemented multiple boson sampling protocols in a single silicon chip.

Benchmarked a quantum algorithm for molecular vibronic spectra.

Abstract

Implementing large instances of quantum algorithms requires the processing of many quantum information carriers in a hardware platform that supports the integration of different components. While established semiconductor fabrication processes can integrate many photonic components, the generation and algorithmic processing of many photons has been a bottleneck in integrated photonics. Here we report the on-chip generation and processing of quantum states of light with up to eight photons in quantum sampling algorithms. Switching between different optical pumping regimes, we implement the Scattershot, Gaussian and standard boson sampling protocols in the same silicon chip, which integrates linear and nonlinear photonic circuitry. We use these results to benchmark a quantum algorithm for calculating molecular vibronic spectra. Our techniques can be readily scaled for the on-chip implementation of specialised quantum algorithms with tens of photons, pointing the way to efficiency advantages over conventional computers.

Figures

Tables

Table 1. Table S1: Estimated count rates near-term experiments using current technologies. The rates are calculated considering a silicon device where integrated sources deliver quantum states to a 100-modes interferometer ( m = k = 100 𝑚 𝑘 100 m=k=100 ). This configuration requires in total 200 detectors for SBS (100 for the signal modes and 100 for the idler modes), and 100 detectors for GBS. Note that in SBS n 𝑛 n signal photons detection correspond to 2 ​ n 2 𝑛 2n -fold events including the n 𝑛 n idler photons. The calculation, via Eq.( S40 ), is performed using the average values for coupling, filtering, on-chip and off-chip transmission losses, and detection efficiency as characterised in our experiment.

Number of signal photons	Event Rate (Hz)
	Spiral Sources		Spiral Sources & Int.Det.		Ring Sources		Ring Sources & Int.Det.
	SBS	GBS	SBS	GBS	SBS	GBS	SBS	GBS
10	$0.3$	$6\times {10}^3$	$2\times {10}^3$	$5\times {10}^5$	65	$6\times {10}^4$	$4\times {10}^5$	$5\times {10}^6$
16	-	1	$8\times {10}^{-3}$	$1\times {10}^3$	$2\times {10}^{-3}$	$4\times {10}^2$	$3\times {10}^3$	$4\times {10}^5$
20	-	$2\times {10}^{-3}$	-	$15$	-	$7$	30	$5\times {10}^4$
32	-	-	-	$4\times {10}^{-6}$	-	$1\times {10}^{-5}$	-	15
40	-	-	-	-	-	-	-	$2\times {10}^{-2}$

Equations

$$p(\bar{\mathcal{M}}|D)$$

$$=\frac{1}{1+\frac{p(D|{\mathcal{M}}_1)p_0({\mathcal{M}}_1)}{p(D|\bar{\mathcal{M}})p_0(\bar{\mathcal{M}})}}.$$

$$\frac{1}{1+\frac{p(D|{\mathcal{M}}_{k_0})}{p(D|\bar{\mathcal{M}})}}\ge p(\bar{\mathcal{M}}|D)$$

$$\ge \frac{1}{1+m\frac{p(D|{\mathcal{M}}_{k_{\text{max}}})}{p(D|\bar{\mathcal{M}})}},$$

$$C_{i+1}=\{\begin{array}{ll}{C}_i+1,& \text{if }{R}_{\vec{j},\vec{k}}>(n/m{)}^m\\ C_i-1,& \text{if }{R}_{\vec{j},\vec{k}}\le (n/m{)}^m.\end{array}$$

$$C_{i+1}=\{\begin{array}{ll}{C}_i+2,& \text{if }{L}_{\vec{j},\vec{k}}\ge a_2\\ C_i+1,& \text{if }1/a_1\le L_{\vec{j},\vec{k}}<a_2\\ C_i,& \text{if }a1<L_{\vec{j},\vec{k}}<1/a_1\\ C_i-1,& \text{if }1/a_2<L_{\vec{j},\vec{k}}\le a_1\\ C_i-2,& \text{if }{L}_{\vec{j},\vec{k}}\le 1/a_2.\end{array}$$

$$Q(\vec{\alpha })=\frac{1}{{\pi }^m\sqrt{\text{det}{\sigma }_{\text{Q}}}}\exp \left[-\frac{1}{2}{\vec{\mathbf{\alpha }}}^{†}{\sigma }_{\text{Q}}^{-1}\vec{\mathbf{\alpha }}\right],$$

$${\sigma }_{\text{Q}}=\frac{1}{2}\left[\begin{array}{cc}U& 0\\ 0& U^{\ast }\end{array}\right]S{S}^{†}\left[\begin{array}{cc}{U}^{†}& 0\\ 0& U^{\text{T}}\end{array}\right]+I_{2m}/2,$$

$$S=\left[\begin{array}{cc}{\oplus }_{j=1}^m\cosh {\xi }_j& {\oplus }_{j=1}^m\sinh {\xi }_j\\ {\oplus }_{j=1}^m\sinh {\xi }_j& {\oplus }_{j=1}^m\cosh {\xi }_j\end{array}\right]$$

$$p_{\text{gbs}}(\vec{k})={\pi }^m{\int }_{{\mathbb{C}}^m}\prod _{i=1}^{m}d{\alpha }_{i}Q(\vec{\alpha })P_{|\vec{k}⟩⟨\vec{k}|}(\vec{\alpha }),$$

$$P_{|\vec{k}⟩⟨\vec{k}|}(\vec{\alpha })=\prod _{i=1}^m{e}^{|{\alpha }_i{|}^2}{\left(\frac{{\partial }^2}{\partial {\alpha }_i\partial {{\alpha }_i}^{\ast }}\right)}^{k_i}\delta ({\alpha }_i)\delta ({\alpha }_i^{\ast })$$

$$p_{\text{gbs}}(\vec{k})=\frac{|\text{Haf }{B}_{\vec{k}}{|}^2}{k_1!\cdots k_m!\sqrt{\text{det}{\sigma }_{\text{Q}}}}.$$

$$\text{Haf}X=\sum _{\mu \in {\text{C}}_{2K}}\prod _{j=1}^k{X}_{\mu (2j-1),\mu (2j)},$$

$$\text{Perm}X=\sum _{\mu \in {\text{P}}_k}\prod _{j=1}^k{X}_{1,\mu (j)},$$

$$p_{\text{sbs}}(\vec{k})=\frac{|\text{Perm}{U}_{\vec{k},\vec{j}}{|}^2}{k_1!\cdots k_M!}.$$

$$|{\alpha }_i⟩=e^{-\frac{|\alpha {|}^2}{2}}\sum _{n_i=0}^{\infty }\frac{{\alpha }^{n_i}}{\sqrt{n_i!}}|n_i⟩.$$

$${\beta }_k=\sum _{j=1}^M{U}_{kj}{\alpha }_j.$$

$$p_{\text{cs}}(\vec{k})=\prod _{i=1}^m\frac{e^{-|{\alpha }_i{|}^2}|{\alpha }_i{|}^{2n_i}}{n_i!},$$

$${\rho }_i^{\text{th}}=(1-{\tau }_i)\sum _{n_i=0}^{\infty }{\tau }_i^{n_i}|n_i⟩⟨n_i|,$$

$$p_{\mathrm{th}}(\vec{k})=\frac{1}{k_1!\cdots k_m!}\frac{\mathrm{Perm}{A}_{\vec{k}}}{{\prod }_{i=1}^m(1+⟨n_i⟩)},$$

$$A=UD{U}^{†},$$

$$D=\text{diag}\left\{{\tau }_1,...,{\tau }_m\right\}$$

$${\tilde{\sigma }}_{\text{Q}}^{(j)}=\frac{1}{2}\left[\begin{array}{cc}U& 0\\ 0& U^{\ast }\end{array}\right]{\tilde{S}}_j{\tilde{S}}_j^{†}\left[\begin{array}{cc}{U}^{†}& 0\\ 0& U^{\text{T}}\end{array}\right]+I_{2m}/2$$

$${\tilde{S}}_j=\left[\begin{array}{cc}{\oplus }_{i=1}^m(1+{\delta }_{j,i}\cosh {\xi }_j)& {\oplus }_{i=1}^m{\delta }_{j,i}\sinh {\xi }_j\\ {\oplus }_{i=1}^m{\delta }_{j,i}\sinh {\xi }_j& {\oplus }_{i=1}^m(1+{\delta }_{j,i}\cosh {\xi }_j)\end{array}\right].$$

$$p_{\text{dss}}^{(4)}(\vec{k})=$$

$$p_{\text{tms}}(\vec{x})=\sum _{\begin{array}{c}\vec{h},\vec{k}:\\ \vec{x}=\vec{h}+\vec{k},\\ {\sum }_i{h}_i={\sum }_i{k}_i=n/2\end{array}}{p}_{\text{gbs}}(\vec{h},\vec{k}),$$

$$U_{\text{ps}}=\left[\begin{array}{cc}1& 0\\ 0& i\end{array}\right],U_{\text{bs}}=\left[\begin{array}{cc}1& 1\\ -1& 1\end{array}\right]/\sqrt{2}.$$

$${\sigma }^{(2k)}=\frac{1}{2}\left[\begin{array}{cc}{U}_{\text{bs}}& 0\\ 0& U_{\text{bs}}^{\ast }\end{array}\right]S_{2k}{S}_{2k}^{†}\left[\begin{array}{cc}{U}_{\text{bs}}^{†}& 0\\ 0& U_{\text{bs}}^{\text{T}}\end{array}\right]$$

$$S_{2k}=\left[\begin{array}{cc}{\oplus }_{j=1}^{2k}\cosh {\xi }_j^{\prime }& {\oplus }_{j=1}^{2k}\sinh {\xi }_j^{\prime }\\ {\oplus }_{j=1}^{2k}\sinh {\xi }_j^{\prime }& {\oplus }_{j=1}^{2k}\cosh {\xi }_j^{\prime }\end{array}\right],U_{\text{bs}}={\left[\begin{array}{cc}\sqrt{\eta }& \sqrt{1-\eta }\\ -\sqrt{1-\eta }& \sqrt{\eta }\end{array}\right]}^{\oplus 2k}.$$

$$p_{\text{sms}}^{(k,\xi )}(2n)=\left(\genfrac{}{}{0ex}{}{k/2+n-1}{n}\right){\text{sech}}^k(\xi ){\tanh }^{2n}(\xi ).$$

$$p_{\text{tms}}^{(k,\xi )}(n)=\left(\genfrac{}{}{0ex}{}{k}{n}\right){\text{sech}}^{2k}(\xi ){\tanh }^{2n}(\xi ).$$

$$p_{\times }^{(k,\xi )}(2n)$$

$$={\text{sech}}^{5k}(\xi ){\tanh }^{2n}(\xi )\sum _{m=1}^n\left(\genfrac{}{}{0ex}{}{k/2+n-m-1}{n-m}\right)\left(\genfrac{}{}{0ex}{}{2k}{2m}\right){\tanh }^{2m}(\xi ).$$

Full text

††preprint: APS/123-QED

Devices that address customised problems with quantum algorithms are expected to demonstrate an efficiency advantage over conventional computers. Boson sampling is a specific model of quantum computing that is suited to the platform of photonics Aaronson and Arkhipov (2011); Broome et al. (2013); Spring et al. (2013); Crespi et al. (2013); Tillmann et al. (2013) and has been mapped to the calculation of molecular vibronic spectra Huh et al. (2015), simulation of spin Hamiltonians González Olivares et al. (2016), simulation of molecular quantum dynamics Sparrow et al. (2018), and the enhancement of classical optimisation heuristics Arrazola and Bromley (2018). Implementing such applications at a size that challenges conventional computers Neville et al. (2017) demands the integration and high fidelity operation of a large number of different components, including circuitry Carolan et al. (2015); Harris et al. (2017); Wang et al. (2018a), detectors Pernice et al. (2012), filters Harris et al. (2014), and photon sources Silverstone et al. (2014); Spring et al. (2017). The low efficiency of individual spontaneous photon sources has motivated the adoption of deterministic solid-state photon sources Somaschi et al. (2016); Loredo et al. (2017); Wang et al. (2017, 2018b). However, the low-loss integration of solid-state sources into photonic circuitry is an on-going challenge.

Creative approaches to realise boson sampling with high numbers of photons from spontaneous sources have seen the design of variant models. In principle, the simultaneous optical pumping of a number $k$ of spontaneous sources that exceeds the number $n$ of desired photons, boosts the overall rate of photon-pair production combinatorially. In Scattershot boson sampling (SBS) Lund et al. (2014); Bentivegna et al. (2015); Zhong et al. (2018), one photon from each pair heralds the location of its partner, such that a Fock state of $n$ photons is prepared over a random subset of modes. In Gaussian boson sampling (GBS) Hamilton et al. (2017), $k$ single mode optical nonlinearities are coherently pumped to produce $k$ modes of squeezed (vacuum) light, before linear optical processing and $n$-photon detection at the output.

The complex photonic circuitry required to scale these approaches can be addressed with integrated photonics. Here, by pumping 4 integrated spontaneous four wave mixing (SFWM) sources with either a single-colour or a two-colour laser, we select between the Fock state required for SBS and the squeezed state required for GBS. Both states of light are routed on-chip to the same linear optical circuit, which implements a random unitary operation over 12 waveguides. In the GBS mode of operation, we benchmark this class of device for calculation of vibronic spectra. In the limit of the SBS mode of operation, with $n=k=4$, we implement standard boson sampling with 4 heralded photons, generating and processing 8 photons on-chip. Our analysis shows that larger versions of our silicon photonic chip, that exploit the combinatorial boost in photon rate available through the SBS and GBS protocols, open up the regime of efficiency advantages over conventional computers.

The silicon circuitry and configuration of these experiments can be understood with reference to Fig. 1. Four SFWM spiral sources are coherently pumped by on-chip splitting of the near-$1550\text{ nm}$ pump laser via multi-mode interference (MMI) near 50:50 beam-splitters; pump light is then removed by asymmetric Mach-Zehnder interferometers (AMZIs). In the dual-wavelength pumping (GBS) scheme (where photons are generated at the same signal wavelength) two spectral regions of the temporally compressed (spectrally broadened) pump are selected and recombined using wavelength-division multiplexers (WDMs), before injection to the chip. In the single-wavelength pumping (SBS) scheme, (where signal and idler photons are generated at different wavelengths) idler photons are separated using a second layer of AMZIs and out-coupled to a fibre array to herald the presence of the signal photons. In both regimes, signal photons are routed to the four central modes of a continuous random walk, implemented over 12 evanescently coupled waveguides, then out-coupled to the fibre array. Ultra low-loss out-coupling is implemented by aluminium assisted apodized grating couplers Ding et al. (2014, 2013). An array of 16 superconducting nano-wire single-photon detectors (SNSPDs) with approximately $78\%$ efficiency detect the 4 heralding modes and the 12 modes from the random walk. (See Supplementary Materials for details).

The SBS protocol is designed to tackle the inefficiency that results from implementing standard boson sampling with spontaneous photon sources Lund et al. (2014). In the original model of boson sampling Aaronson and Arkhipov (2011), a linear optical circuit is configured to implement a random unitary operation over $m$ modes. A number $n<m$ of input ports are each populated with a single photon, such that $n$ photons undergo an $n\times m$ random operation, before detection with photon counters. Because the probability amplitude for each $n$-photon transition is equal to the permanent of the corresponding transfer matrix, which is in general intractable to classical computation, an ideal experiment would efficiently produce samples from an essentially classically forbidden probability distribution. However, the rate at which the required $n$-photon state is produced with spontaneous sources decreases exponentially with $n$. As shown in Fig. 1c, the SBS protocol addresses this inefficiency by pumping $k$ sources that each produce weak two-mode squeezed light (with $n<k\leq m$). One (idler) photon from each pair is directly detected to identify the input port of its partner (signal) photon. In comparison to the original scheme, the probability to generate $n$-photons is boosted by a factor $\propto{k\choose{n}}$  Lund et al. (2014); Bentivegna et al. (2015).

We implemented the SBS protocol in our silicon chip, with $k=4$ sources, using the single wavelength pumping scheme shown in Fig. 1b. The experimentally measured distribution for the $n=3$ case (6-photon events) shown in Fig. 2a, has a mean fidelity Carolan et al. (2015) with the theoretical distribution of $92\%$. Figure 2b shows the measured difference in photon pair generation rates between the SBS and standard boson sampling (where only $n$ sources are pumped) protocols for $n=1\text{ to }3$, and for $n=4$ where the two regimes converge. As predicted, the $n$-photon pair generation rate, is enhanced by a factor of approximately $\binom{4}{n}$. In the SBS regime, the two-pair and three-pair photon rates were measured at $5.8\text{ kHz}$ and $4\text{ Hz}$, respectively, while the 4-pair (8-photon) rate was approximately 4 events per hour. The average purity of the signal photons, estimated via unheralded second-order correlation measurements Christ et al. (2011), was calculated at $86\%$. (See Supplementary Materials for details).

While fidelity comparisons with theoretical distributions are not a scalable method of validating boson sampling, an efficient alternative is to compare measurements against those predicted by a classically tractable and plausible distribution. As shown in Fig. 2c, we used Bayesian model comparison Bentivegna et al. ; Spagnolo et al. (2014); Carolan et al. (2014); Bentivegna et al. (2015); Carolan et al. (2015) to validate our $n=4$ (8-photon) statistics against those predicted by distinguishable photons (the $n=3$ case is inset). The dynamically updated confidence that samples are drawn from the distribution of indistinguishable particles versus distinguishable particles reaches higher than 99.9%. (See Supplementary Materials for details).

In contrast to SBS, the GBS protocol does not project the input state onto a single Fock state. Rather, the input in GBS is an ensemble of single mode squeezed states, as shown in Fig. 1d, which further increases the $n$-photon detection probability as compared to SBS Hamilton et al. (2017). After processing with linear optical circuitry the probability for a particular pattern at single photon detectors is given by a function known as the Hafnian of the relevant transfer matrix. Similar to the permanent of a matrix, the Hafnian is computationally hard to calculate  Caianiello (1973); Hamilton et al. (2017); Kruse et al. (2018). (See Supplementary Materials for details).

We implemented the GBS protocol using the two-colour pumping scheme described in Fig. 1b, which generates weak single-mode squeezed light at each source. While this spectral selection reduces the pump power below that of the single-colour configuration, we observed statistics with up to 4 signal photons at $1.1$ Hz rate. The experimentally measured distribution for $n=4$ photons, with one photon per detector (full anti-bunching), shown in Fig. 3a, has a mean statistical fidelity of $87\%$ with the ideal theoretical distribution. In validating GBS, a wide range of alternative models is available where output distributions arising from general Gaussian input states are classically tractable. As illustrated in Fig. 3b, we focus on models that are plausible in our experiment. We validate the ideal input state of 4 single-mode squeezed states against input with: thermal states, resulting from excessive loss; coherent states, from unfiltered pump light; distinguishable single-mode squeeze light, due to spectral mismatch; and two-mode squeezed states, from spurious photons generated at different wavelengths. For each test, a confidence $>99.9\%$ is reached after $\approx 120$ samples for an ideal model using single mode squeezed states rather than the alternative models. (See Supplementary Materials for details).

An application of GBS is the calculation of molecular vibronic spectra Huh et al. (2015); Clements et al. (2017); Shen et al. (2018), where programmable linear optical circuitry, together with squeezed and displaced light, can approximate probabilities for the vibrational transition (Frank-Condon factors) between ground and excited states of a given molecule. To investigate the performance of silicon photonic chips for this application, we mapped our random walk circuit to a synthetic molecule and considered the difference $\mathcal{C}=F_{Q}-F_{C}$ between the fidelity of the reconstructed Franck-Condon (FC) profile $F_{Q}$ and the optimal fidelity obtainable from a classical strategy $F_{C}$ Clements et al. (2017). While the FC-profile directly reconstructed from GBS measurements has a fidelity $>99\%$ to the ideal FC-profile, the low level of squeezing means that classically-tractable vacuum contributions dominate and the improvement over a classical strategy is small, ($\mathcal{C}=0.4\%$).

Data post-processing, based on characterised losses and detection efficiencies, makes it possible to investigate how the fidelity of the FC-profile depends on the amount of simulated squeezing, as shown in Fig. 4(inset). For moderate levels of simulated squeezing, an improvement of up to $\mathcal{C}=9\%$, corresponding to $F_{Q}=86\%$ (shown in Fig. 4), is obtained from post-processing. For high values of simulated squeezing, due to the increased contribution of higher order photon terms, the approach does not provide any further advantage. (See Supplementary Materials for details).

Compared to SBS with 4 heralded photons, we observed much higher rates for 4-photon GBS using a pump with lower power, and requiring fewer detectors. This type of resource saving for GBS over SBS increases as the size of demonstration increases, though GBS requires approximately twice as many photons as SBS to demonstrate efficiency advantages over classical algorithms. Increasing the number and purity of integrated photon sources Caspani et al. (2017) while decreasing the loss in photonic circuitry to reduce noise and distinguishability among photons, leads to similar fidelities when considering the scaling of both the GBS and SBS protocols (e.g. loss acts to add noise via thermal components to the initial SBS-two-mode or GBS-single-mode squeezed state resource states). The maximum number of photons we generated and processed in these experiments is 8, which is double the largest reported to date in integrated photonics. Based on these results and further analyses (see Supplementary Materials), we estimate that arrays of several hundred integrated detectors, possible with current fabrication technologies Schuck et al. (2013), would allow fully-integrated experiments with tens of photons. In the context of calculating molecular transition probabilities, our class of photonic chip will be useful for reconstructing Franck-Condon profiles from photonic quantum sampling algorithms, if relevant instances requiring $>20$ photon events are discovered. More generally, for problems that can be mapped to variant models of boson sampling, such as interrogating the vibrational dynamics of molecules Sparrow et al. (2018), our results and analysis show that efficiency advantages over conventional computers are a realistic prospect with the platform of integrated photonics.

Appendix A Supplementary Materials

Appendix B Experimental setup details

A detailed schematic of the experimental setup and images of the integrated device and the packaging are shown in Fig. S1. Single photons were generated and guided in silicon waveguides with a cross-section of $450\text{ nm}\times 250\text{ nm}$, using spontaneous four-wave mixing (SFWM). In the $\chi^{(3)}$ SFWM process weak squeezed states are generated as pump light is transmitted through four $1.4\text{\ cm}$ long spirally-shaped waveguides Silverstone et al. (2014), see Fig. 1. In this process photons from the pump are annihilated while correlated pairs of photons are generated. Each of the 8 asymmetric Mach-Zehnder interferometer (AMZI) on-chip filters consist of two multi-mode interferometers (MMIs), acting as near-50% beam splitters, and a thermal phase shifter on the shorter arm of the Mach-Zehnder interferometer. The interferometer is implemented via a continuous random walk obtained via evanescently coupling 12 waveguides, with the four sources injecting photons in the four central modes. The distance between near-neighbour waveguides in the random walk is 450 nm and the coupling region has a length of $110\mu\text{m}$, ensuring coupling between all 12 waveguides. In order to fine tune the AMZI filters, the chip was wire-bounded to a PCB to individually address each thermal phase shifter via an electronic voltage driver with 12 bits resolution (Qontrol™).

To pump the sources, laser pulses were generated via a tunable laser (PriTel™Ultrafast Optical Clock) with a repetition rate of 500 MHz, emitting 2 nm bandwidth transform limited pulses at near-telecom wavelengths, and amplified via an Erbium-doped fiber amplifier (EDFA). In the single-wavelength pumping scheme, the emitted pulses (central wavelength of 1542.9 nm) are filtered off-chip via a wavelength-division multiplexer (WDM) ($>100$ dB extinction, 1.6 nm bandwidth) to remove spurious tails in the pump spectrum (see Fig. S2a), and injected into the chip via a single mode fiber. Photons are emitted via non-degenerate SFWM at $1549.3$ nm (signal wavelength) and $1536.6$ nm (idler wavelength). In the dual-pump regime, the pulses, with a central wavelength of $1549.3$ nm are temporally compressed via a fiber optical pulse compressor (PriTel™Femtopulse${}^{\text{\textregistered}}$), achieving a 10 nm wide spectrum from which two slices (at 1552.5 nm and 1546.1 nm, each with 1.6 nm bandwidth) were selected via the WDM (see Fig. S2b). A piezo-controlled optical delay line is used to ensure temporal overlap between the two obtained pulses of different wavelengths, which are then multiplexed together into a single mode fiber via a second WDM and injected into the chip. Photons are emitted on-chip via degenerate SFWM at the signal wavelength (1549.3 nm). In both pumping schemes, fiber polarization controllers were used to maximise the coupling. The input power is monitored via the use of a $99:1$ fiber beam-splitter and a photo-diode just before injecting the pump light into the chip.

Photons are collected out of the chip via a 16-channel fiber array, filtered via off-chip single channel fiber filters ($>100$ dB extinction, 0.8 nm bandwidth, 0.3 dB insertion loss on average) to remove spurious pump photons and enhance the photon indistinguishability, and finally sent to an array of 16 SNSPDs (Photon Spot™, $\sim 100$ Hz dark counts, $\sim 80\%$ efficiency) for detection. Fast counting logic (Swabian Instruments™GmbH), supporting more than 40 million events per second, was used to collect the single photon detections and process them on the fly.

Appendix C Chip design, fabrication, and components characterisation

Optical losses, in particular insertion loss of integrated optical components and coupling loss, play a critical role for on-chip experiments. On our device, in order to decrease the coupling loss, a fully-etched apodized grating coupler using a photonic crystal is designed Ding et al. (2013), and aluminium (Al) mirrors are used below the grating couplers via a flip-bonding process Ding et al. (2013). The fabrication process starts from a commercial silicon-on-insulator (SOI) wafer with top silicon thickness of 250 nm and buried oxide layer of 3 $\mu$m, and proceeds as follows. Firstly, 1.6 $\mu$m thick SiO2 is deposited by the plasma-enhanced chemical vapour deposition (PECVD) process on the SOI wafer. After that, the Al mirror is deposited by electron-beam (ebeam) evaporator, and followed by another thin layer of SiO2 deposition with thickness of 1 $\mu$m. Following that, the wafer is flip-bonded with Benzocyclobutene (BCB) bonding process to another silicon carrier wafer. The substrate and buried oxide (BOX) layers of the original SOI wafer are removed by fast dry-etching and buffered hydrofluoric (BHF) acid chemical-etching, successively, resulting in the final Al-introduced SOI wafer. After that, the silicon photonic circuit is fabricated by standard ebeam lithography (EBL) followed by Inductively Coupled Plasma (ICP) etching and ebeam resist stripping. In order to simplify the fabrication process, grating couplers are designed as fully-etched so that they can be fabricated with the rest fully-etched silicon nanowires and other photonic components simultaneously. After the photonic part is fabricated, 1.3 $\mu$m thick SiO2 is deposited by PECVD, followed by polishing process to planarize the surface with approximately 300 nm sacrifice, resulting in a final top SiO2-cladding layer of 1 $\mu$m. The micro-heaters are patterned afterwards by another EBL process followed by 100 nm titanium (Ti) deposition and liftoff process. The conducting wires and electrode pads are obtained by standard Ultraviolet (UV) lithography followed by Au/Ti deposition and lift-off process.

Our fabrication platform provide a propagation loss of $\sim$2 dB/cm measured by the cut-back method for the fully-etched silicon waveguide with geometry of 450 nm$\times$250 nm. Figure S3 presents the characterisations of the integrated components. The Al-mirror assisted grating couplers exhibit coupling efficiency of -1.1 dB at the wavelengths used, with 1 dB coupling bandwidth of 40 nm, as shown in Fig. S3a. The 2$\times$2 MMI structures are used to implement AMZI filters with a titanium micro-heaters applied on one arm as phase shifters in the AMZIs tunable filters. In this situation, applying a voltage to the Ti micro-heater results in a power dissipated and heating the optical waveguide underneath, changing the refractive index in the waveguide and inducing a phase shift. As shown in Fig. S3b, 3  V heating voltage results in a transmission shift of more than one free-spectral range (FSR). Moreover, the insertion loss presented in the inset of Fig. S3b is less than 0.1 dB, indicating less than 0.05 dB insertion loss for each 2$\times$2 MMI.

Inside the AMZI filters, the phase shifters allow to finely tune the filtered wavelengths. In order to achieve an accurate tuning, a precise characterisation of the heaters is required. The behaviour of the resistance for all heaters are presented in Fig. S4a, indicating a consistent behaviour of all heaters. A non-Ohmic behaviour is observed, as is usual for Ti-based phase-shifters Santagati et al. (2018); Paesani et al. (2017); Wang et al. (2018a), and characterised in order to calculate the voltage required for a given heating power. The phase-shift as a function of applied heating power is finally characterized for all the heaters, showing a consistent heating efficiency of all heaters with capability of $\sim$3$\pi$ for an applied heating power of approximately 18 mW.

In order to implement 4 $\times$12 random walk, a 12 $\times$12 waveguides evanescent coupler is designed. The widths of the 12 waveguides is designed to be 450 nm with coupling gap of 180 nm. The coupling length was chosen to be of $110\mu\text{m}$ in order to ensure a coupling between all 12 waveguides, as was preliminarily determined by three-dimensional finite-difference time-domain (3D FDTD) simulations of the circuit, as reported in Fig. S5a. The 3D FDTD simulations had only the scope to establish the length of the coupling region for the design, and were thus performed considering no phase fluctuations between the modes. In practice, random phase fluctuations naturally arise in silicon circuits due to nano-scale fabrication imperfections, which insert randomness in the unitary evolution of the random walk. Due to such effects it has been shown that the distribution of the obtained transform matrices converges to the Haar random distribution exponentially fast with the interaction length Banchi et al. (2017), making random walks a promising approach for the implementation of low-loss-interferometers for boson sampling Muraleedharan et al. (2018). The transfer matrix of the implemented $4\times 12$ continuous random walk, which is the rectangular $4\times 12$ submatrix given by the 4 central rows of the total $12\times 12$ unitary matrix, is shown in Fig. S5b-c, which was preliminarily characterised using standard methods based on classical and two-photon interference Laing and O’Brien (2012); Peruzzo et al. (2011); Carolan et al. (2014, 2015).

Appendix D System efficiency

Optical losses in the experimental set-up and in the integrated device were preliminarily measured. The average channel efficiency, that is the total loss experienced by each photon from generation to detection, measured via the coincidence to singles ratio Silverstone et al. (2014) and averaged over all 16 channels used, was measured to be $-11.5$ dB. Characterisation of the 16 SNSPDs presented an average efficiency of $0.78$ (-1.0 dB). Mean transmission efficiencies through the off-chip filters and through the fiber connections to the detectors, averaged over the 16 channels used, were measured to be $0.87$ (-0.57 dB) and $0.94$ (-0.26 dB) respectively. Measured efficiency for the grating couplers is $0.77$ (-1.1 dB) at the wavelengths used, and the transmission for the AMZI on-chip filters is $0.98$ (-0.1 dB). Propagation loss through straight waveguides was estimated via cut-back measurements to be approximately 2 dB $\text{cm}^{-1}$, indicating a transmission efficiency of $0.995$ for the 12-mode random walk implemented. Due to bending losses, the propagation loss in the spiral waveguides is higher, with average losses of $-7.1$ dB measured in the 1.4 cm long spirals.

Appendix E Degenerate and non-degenerate SFWM experimental details

The performance of the individual sources was directly characterised by using the control grating couplers after the first array of AMZIs (outputs 1 in Fig S1b), configuring the filters accordingly to collect the photons before the interferometer. Non-degenerate and degenerate SFWM effects were used for the generation of the input states required for SBS and GBS. The pumping schemes used to excite the processes are reported in Fig. 1b and Fig. S1a. For both regimes, the measured two-fold coincidences count rates from each single source are reported in Fig. S6a for different input pump powers (the input pump power is measured off-chip, as shown in S1a). To obtain these curves, for each source the two-fold coincidences are analysed using the same two channels and detectors, so that channel efficiencies are the same for all cases and differences in the curves are only due to the different efficiencies between the sources. The efficiencies of the two channels used, measured as $\eta_{i}=C_{i}/CC$, where $C_{i}$ is the singles count rates on the associated channel and $CC$ are the coincidence counts, are $\eta_{1}=-9.7$ dB and $\eta_{2}=-8.9$ dB.

In the non-degenerate SFWM regime (blue-green curves in Fig. S6a), the coincidence count rates are given in term of the source efficiency (i.e. two-photon emission probability) $\epsilon=\tanh(\xi)^{2}/\cosh^{2}(\xi)$ as $CC=\epsilon R\eta_{1}\eta_{2}$ (see for example Ref. Silverstone et al. (2014, 2015)), where $\xi$ is the squeezing parameter and $R$ is the repetition rate of the laser, from which $\epsilon$ can be calculated. In SBS measurement conditions (9 dBm input pump power) we obtain for the four non-degenerate SFWM sources the two-mode squeezing parameters $\vec{\xi}=\{0.25,0.21,0.18,0.17\}$.

In the scheme used for the dual-pump regime, the pump is initially broadened in spectrum through a pulse compressor and the two different wavelengths are selected via WDM filters. This is done to avoid synchronisation and phase locking issues that would arise when using two different pump lasers directly. While this significantly simplifies the experimental setup, it introduces some issues that have to be taken into account. First, a temporal delay between the two different frequencies is introduced when passing through the WDM filters, which is compensated via a delay line (see Fig. 1b and Fig. S1a). The photon count rates dependence on the temporal delay between the pump pulses with different wavelengths is shown in Fig. S6c. A peak of coincidences is observed when the pulses are overlapped, that is when photons are emitted via degenerate SFWM. A low level of noise can observed when the two pumps are not overlapped, with a signal-to-noise ratio of approximately 30, which is mainly due to multi-pair emission of non-degenerate SFWM photons from the two different pumps. Secondly, by selecting only two frequencies on a broad pump spectrum, we significantly reduce the power injected into the chip, which in turn reduces the photon rates in the degenerate SFWM when compared to the non-degenerate SFWM case. In Fig. S6a the photon pair rates (two-fold coincidences count rates) for each single degenerate SFWM source are reported (red-yellow curves) for different input pump powers. These were measured as in the non-degenerate SFWM case previously mentioned, with a only difference that a 50:50 off-chip fiber beam splitter was used to probabilistically separate the two degenerate photons (which succeeds with a one-half probability). Due to this probabilistic splitting, in the degenerate SFWM regime the coincidence count rates are given in term of the source efficiency $CC=\epsilon R\eta_{1}\eta_{2}/2$ (again, see for example Ref. Silverstone et al. (2014, 2015)). As we are now generating weak single-mode squeezing, the source efficiency (i.e. two-photon emission probability) is given by $\epsilon=\tanh(\xi)^{2}/2\cosh(\xi)$. In the GBS measurement conditions (1 dBm total input pump power, measured off-chip) we obtain for the four degenerate SFWM sources the single-mode squeezing parameters $\vec{\xi}=\{0.11,0.09,0.07,0.07\}$. These squeezing values have been used, together with the characterised unitary of the interferometer, to construct the matrix $B=T[\oplus_{i}\tanh(\xi_{i})]T^{t}$ for GBS. In GBS measurement conditions (1 dBm input power), four-photon events are observed at approximately a $5\text{ Hz}$ rate.

While the count rate obtained using the implemented dual-wavelength pumping scheme for degenerate SFWM are greatly reduced by the low pump power injected, we remark that such scheme offers plenty of room for improvement. For example, the use of two separate pulsed lasers, each tuned at one of the two wavelength desired, would certainly allow much higher pump powers in the dual-wavelength regime, at the cost of requiring precise synchronisation between the two lasers.

Appendix F Photon purity characterisation

In SBS high spectral purity of each signal photon emitted, heralded by the associated idler photon, is important to achieve a good quality quantum interference between the signal photons. While the good fidelities of the output distribution in the boson sampling experiments and the validation tests indicate good photon purities, in the non-degenerate SFWM case, where two-mode squeezed states are generated, a quantitative estimate of the photon purity can be obtained through unheralded second order correlation ($g_{2}$) measurements Christ et al. (2011). To perform these measurements, photons are collected at the output of each source (through the ports labelled as output1 in Fig. S1b), the idler and signal photons are separated through an off-chip filter and second order correlation measurements are performed on the signal photon while the idler is not measured. Following Ref. Christ et al. (2011), the photon purity $\mathcal{P}$ can be estimated from the obtained second order correlation $g_{2}(0)$ via $\mathcal{P}=g_{2}(0)-1$. In Fig. S6b the obtained purities for the four different sources are reported for various values of the input pump power. A decrement of the purity for increasing pump powers is observed, plausibly due to non-linear noises in the sources and multi-pair emission errors, so that a trade-off between photon purity and generation rate has to be chosen. In measurement conditions for SBS via non-degenerate SFWM (9 dBm input pump power) the average purity is approximately $86\%$.

Appendix G Bayesian model comparison for boson sampling validation.

Boson sampling validation consists in providing supporting evidence to show that the measured samples are obtained from a true boson sampler, dictated by quantum interference, rather than from a classically tractable noisy device. In order to have a consistent verification protocol throughout the work, we have performed the validation via Bayesian model comparison techniques adapted from previous works Bentivegna et al. ; Carolan et al. (2014); Spagnolo et al. (2014); Carolan et al. (2015), whose generality allowed to validate all tested models using the same approach. The procedure is as follows. Suppose that, given a set of data $D=\{x_{i}\}_{i=0}^{N}$ constisting of $N$ output samples $x_{i}$ from the interferometer, we want to verify if the data is more likely to arise from an ideal model $\bar{\mathcal{M}}$ or from a test model $\mathcal{M}_{1}$. For example, $\bar{\mathcal{M}}$ can represent an ideal implementation of SBS or GBS, while $\mathcal{M}_{1}$ a classical implementation where all photons are distinguishable. Given the measured set of samples $D$, Bayes’ rule offers an immediate way to estimate the confidence in model $\bar{\mathcal{M}}$, that is probability for model $\bar{\mathcal{M}}$ to represent the underlying experiment, which is given by:

[TABLE]

Assuming statistical independence between different events, the probability $p(D|\bar{\mathcal{M}})$ is given by $p(D|\bar{\mathcal{M}})=\prod_{i=1}^{N}p(x_{i}|\bar{\mathcal{M}})$, where $p(x_{i}|\bar{\mathcal{M}})$ is the probability of obtaining the measured outcome $x_{i}$ according to the ideal model (dictated by the permanent function for SBS and by the Hafnian in GBS), and similarly for $p(D|\mathcal{M}_{1})$.

Note that, in the case where the detection is restricted to collision-free events, the probabilities $p(x_{i}|\bar{\mathcal{M}})$ have to be normalised to take it into account, which corresponds to dividing by the sum of the probabilities over all collision-free events. The reason for this is that the probability of a sample being collision-free differs between the different models; the normalisation constants are thus in general different.

The probabilities $p_{0}(\bar{\mathcal{M}})=p_{0}(\mathcal{M}_{1})=0.5$ represent a prior distribution on the true model, which we here assume to be uniform to avoid any bias. The confidence $p(\mathcal{M}_{1}|D)$ for the adversary model $\mathcal{M}_{1}$ can be calculated in the same way, and both confidences are plotted in Fig. 2c for $\bar{\mathcal{M}}=\text{SBS}$ and in Fig. 3b for $\bar{\mathcal{M}}=\text{GBS}$ as a function of the number of experimental samples, dynamically updating the confidence as new samples were collected. Note that, has the protocol requires to estimate probabilities $p(x_{i}|\bar{\mathcal{M}})$, which is not efficient for classical machines, the approach is not scalable. However, it is general and can be used to validate against any model which allows to calculate probabilities $p(x_{i}|\mathcal{M}_{1})$.

One could also ask if it would be more convenient to obtain a simultaneous verification from a larger set of $m$ test models $\{\bar{\mathcal{M}},\mathcal{M}_{1},\ldots,\mathcal{M}_{m}\}$, instead of performing multiple two-models comparisons $\{\bar{\mathcal{M}},\mathcal{M}_{k}\}$. In such case, considering a uniform prior distribution between all models, Bayes’ rule gives for the confidence in $\bar{\mathcal{M}}$

[TABLE]

with $k_{\text{max}}=\operatorname*{arg\,max}_{k}p(D|\mathcal{M}_{k})$ and $k_{0}\in\{1,\ldots,m\}$ an arbitrary index value. From the equation above a link between the multiple comparisons with single models and the simultaneous comparison with multiple models can be observed. In particular, it can be easily verified that the confidence in $\bar{\mathcal{M}}$ in the simultaneous validation converges to one if and only if also the confidences in all the multiple two-model comparison validations converge to one.

Appendix H Scattershot boson sampling

In the standard approach to boson sampling Aaronson and Arkhipov (2011), $n$-photon Fock states are generated by $n$ single photon sources and injected into the interferometer. However, if each source has a probability $\epsilon<1$ to generate a photon, then the probability that all $n$ sources fire is given by $\epsilon^{n}$, which decreases exponentially with $n$ and prevents any exponential quantum advantage. The SBS approach has been proposed to avoid these limitations and to increase the complexity of photonic experiments with realistic non-ideal photon-pair sources, such as those based on spontaneous parametric down-conversion or SFWM Lund et al. (2014); Bentivegna et al. (2015). As shown in Fig. 1c, in the SBS scenario, $k>n$ two-mode squeezing sources are simultaneously pumped, and for each source a possible detection of an idler photon on the idler mode heralds the presence of a signal photon in the associated mode. In this way a $n$ photon state is generated whenever a random, but heralded, set of exactly $n$ between the $k$ sources fire, which can happen with an enhanced probability $\binom{k}{n}\epsilon^{n}(1-\epsilon)^{k}$, where now $\epsilon=\tanh(\xi)^{2}$ with $\xi$ the squeezing parameter. In the limit $k\gg n$, this probability represents an exponential combinatorial speed-up in the generation rate with respect to the standard approach, indicating that scalability can in principle be achieved even in presence of non-deterministic photon sources Lund et al. (2014). As the input state is an heralded $n$ photons Fock state, output probabilities are measured in the same way as for standard boson sampling, via the squared permanent of the submatrix associated to the input/output pattern Aaronson and Arkhipov (2011); Lund et al. (2014).

Appendix I Additional validation tests for scattershot boson sampling

In Fig. 2c validation against a distinguishable sampler are reported for the SBS 6-fold and 8-fold data, using Bayesian model comparison techniques. The results give supporting evidence for the correct functioning of the experiment. To give further supporting evidence, we also performed validation protocols based on the likelihood ratio test and the row-norm estimator test, which are standards in previous experimental works Spagnolo et al. (2014); Bentivegna et al. (2015); Wang et al. (2017); Loredo et al. (2017); Wang et al. (2018b).

The row-norm estimator test, originally proposed by Aaronson and Arkiphov Aaronson and Arkhipov (2014), is designed to validate the outputs of a boson sampler from those obtained from a trivial uniform distribution. In contrast to the other validation protocols, this test can be efficiently performed (i.e. is scalable) as it does not require the calculation of ideal output probabilities. Considering a boson sampling experiment where $m$ output modes and $n$ photons are used, the row-norm estimator test proceeds as follows. The counter $C$ is initialised to $C=0$ and is iteratively updated each time an input/output sample $\{\vec{j},\vec{k}\}$ (here, $\vec{j}=\{j_{1},\dots,j_{m}\}$ and $\vec{k}:=\{k_{1},\dots,k_{m}\}$) is obtained from the experiment via calculating the row-norm estimator $R_{\vec{j},\vec{k}}:=\prod_{p}\sum_{q}|T_{\vec{j},\vec{k}}|^{2}$, which is efficiently computable Aaronson and Arkhipov (2014), and using the updating rules

[TABLE]

After the data is collected, the protocol concludes that the outcomes are drawn by a boson sampler if the final value of the counter is $C>0$, otherwise it concludes that the outcomes are drawn from a uniform distribution. The results of the row norm test using 4-fold and 6-fold data are shown in Fig. S7a.

The likelihood ratio test is an alternative approach to distinguish an ideal boson sampler from a device which operates using distinguishable photons Spagnolo et al. (2014); Bentivegna et al. (2015), and works as follows. For an ideal sampler and a distinguishable sampler the probability to observe an output state $\vec{k}$ given the input state $\vec{j}$ is given by $p_{\mathrm{ind}}(\vec{k}|\vec{j})=|\text{Perm}(T_{\vec{j},\vec{k}})|^{2}$ and $p_{\mathrm{dist}}(\vec{k}|\vec{j})=\text{Perm}(|T_{\vec{j},\vec{k}}|^{2})$, respectively. In the likelihood ratio test, a counter $C$ is initialised to $C=0$ and is iteratively updated each time an input/output sample $\{\vec{j},\vec{k}\}$ is obtained from the experiment via calculating the estimator $L_{\vec{j},\vec{k}}=p_{\mathrm{ind}}(\vec{k}|\vec{j})/p_{\mathrm{dist}}(\vec{k}|\vec{j})$ and using the updating rules

[TABLE]

After the data is collected, the protocol concludes that the outcomes are drawn by a boson sampler if the final value of the counter is $C>0$, otherwise it concludes that the outcomes are given by sampling distinguishable photons. In our tests, for the control parameters we used $a_{1}=0.75$ and $a_{2}=2$, with validation results from 4-fold and 6-fold experimental data reported in Fig. S7b.

Appendix J Gaussian boson sampling

The task of Gaussian boson sampling consists in generating samples from the photon-counting probability distribution of linear-optically evolved vacuum squeezed states. That is, we here deal with the probability $p_{\text{gbs}}(\vec{k})$ to detect a pattern of single photons $\vec{k}=\{k_{1},\dots,k_{m}\}$ $(\sum_{i=1}^{m}k_{i}=n)$ at the output of a linear-optical circuit injected with $m$ vacuum squeezed states $\otimes_{i=1}^{m}\ket{\xi_{i}}$. In what follows, we shall characterize each state $\ket{\xi_{i}}$ in terms of its squeezing parameter $0\leq\xi_{j}<1$, while the linear-optical circuit in terms of the unitary matrix $U$ that transforms the input-mode operators $a_{j}^{\dagger}$ into output mode operators $b_{k}^{\dagger}$, $b_{k}^{\dagger}=\sum_{j}U_{kj}a^{{\dagger}}_{j}$.

In the context of Gaussian boson sampling, it is useful to describe the linear-optically evolved state $\otimes_{i=1}^{m}\ket{\xi_{i}}$ via its Husimi Q-function Leonhardt (2010):

[TABLE]

where $\vec{\boldsymbol{\alpha}}=(\alpha_{1},\dots,\alpha_{m},\alpha_{1}^{*},\dots,\alpha_{m}^{*})^{\text{T}}$, $\sigma_{\text{Q}}$ is the $2m\times 2m$ Q-covariance matrix of the evolved state,

[TABLE]

and ${I}_{2m}$ is the $2m\times 2m$ identity matrix. Next, using the formulas Leonhardt (2010)

[TABLE]

and integrating Eq. (S3) by parts using the generalized Faá di Bruno’s formula for partial derivatives Hardy (2006), we obtain the following expression Hamilton et al. (2017); Kruse et al. (2018),

[TABLE]

Here, $B_{\vec{k}}$ is a $n\times n$ matrix obtained from the $m\times m$ matrix $C={U}\cdot\text{diag}[\tanh\xi_{1},\dots,\tanh\xi_{m}]\cdot{U}^{\text{T}}$ by repeating $k_{i}$ times its $i$th column and row. In turn, the Hafnian of a $2k\times 2k$ matrix $X$ is defined as Caianiello (1973)

[TABLE]

where C2k is the set of canonical permutations on $2k$ elements, obeying $\mu(2j-1)<\mu(2j)$ and $\mu(2j)<\mu[2(j+1)]$, $\forall j$. In fact, the Hafnian can be seen as a generalization of another matrix function, the permanent, which governs the statistics of indistinguishable bosons and is at the heart of the complexity of standard boson sampling models,

[TABLE]

where summation runs over all permutations of the numbers $\{1,2,...,k\}$. The computation of both matrix functions, permanents and Hafnians, is hard for a classical computer. Consequently, since the photon counting probability distribution $p_{\text{gbs}}(\vec{k})$ is given in terms of Hafnians, under certain complexity-theoretic assumption, the task of Gaussian boson sampling is intractable for a classical computer Hamilton et al. (2017); Kruse et al. (2018).

Interestingly, scattershot boson sampling can be seen as a special case of Gaussian boson sampling. Namely, one can interfere pairs of adjacent vacuum squeezed states (input to the Gaussian boson sampler) on a beam splitter and use one of the arms of emerging two-mode squeezed states for heralding, while sending the other arms into a linear optical circuit $U$. This setting, being a special case of Gaussian boson sampling, results in a scattershot scenario Chakhmakhchyan and Cerf (2017); Lund et al. (2014). The corresponding photon counting probability distribution $p_{\text{sbs}}({\bf n})$ is now given in terms of permanents of the unitary matrix $U$,

[TABLE]

Here, $U_{\vec{k},\vec{j}}$ is obtained from $U$ by repeating $k_{i}$ times its $i$th row and $j_{i}$ times its $i$th column, where the set $\vec{k}=\{k_{1},\dots,k_{m}\}$ corresponds to the heralded pattern of single photons input to the circuit $U$.

We also remark that the linear-optical transformation that we deal with in our experiments has 12 output but 4 input modes and it is therefore given by a $4\times 12$ unitary transfer matrix $T$ (rather than a $12\times 12$ square matrix $U$). Such a scenario is equivalent to a standard setting with 12 input and output modes described by a $12\times 12$ square matrix $U$, just that all the eight extra input modes are injected with the vacuum. Effectively, this means that the corresponding rows of $U$ are deleted and we always deal with its $4\times 12$ submatrix $T$ obtained selecting the 4 central rows, which represent the transfer matrix of interest in our experiment.

Appendix K Test models for validating Gaussian boson sampling

Bayesian approaches allow to validate experimental data by comparing an ideal experimental implementation of GBS against a general test model for which output probabilities are computable. In this section we describe in more details the test models used and how output probabilities can be computed and used in the validation tests.

Coherent states: We first consider a set of $m$ input coherent states $\otimes_{i=1}^{m}\ket{\alpha_{i}}$, each of which has the following representation in the Fock basis,

[TABLE]

A linear optical circuit $U$ transforms a tensor product of coherent states $\otimes_{i=1}^{m}\ket{\alpha_{i}}$ into another tensor product of coherent states $\otimes_{i=1}^{m}\ket{\beta_{i}}$ with amplitudes

[TABLE]

In other words, coherent states remain separable while evolved through a linear optical circuit. Thus, the joint probability $p_{\mathrm{cs}}(\vec{k})$ of detecting a pattern of $n$ single photons $\vec{k}=\{k_{1},\dots,k_{m}\}$, with an $m$-mode coherent state $\otimes_{i=1}^{m}\ket{\alpha_{i}}$ at the input, admits a simple product form

[TABLE]

which is a product of Poisson distributions. Consequently, both sampling from it and computing its elements can be done efficiently on a classical computer. Although this problem is trivial from a computational point of view, boson sampling from coherent states comprises a physically relevant test model for validating our GBS experiments. In the model implemented the input coherent states amplitudes were considered to be uniform between all modes, as is physically plausible.

Thermal states: Our next test model deals with sampling from linear-optically evolved thermal states. That is, we consider an $M$-mode input thermal state $\otimes_{i=1}^{m}\rho_{i}^{\text{th}}$. In turn, each state $\rho_{i}^{\text{th}}$ can be expressed as an incoherent mixture of Fock states,

[TABLE]

where $\tau_{i}=\langle{n_{i}}\rangle(\langle{n_{i}}\rangle+1)$ and $\langle{n_{i}}\rangle$ is the average photon number of the state. The probability to detect a set $\vec{k}=\{k_{1},\dots,k_{m}\}$ of $n$ single photons at the output of a linear-optical circuit $U$ injected with an $m$-mode thermal state $\otimes_{i=1}^{m}\rho_{i}^{\text{th}}$ then reads

[TABLE]

where

[TABLE]

and ${A}_{\vec{k}}$ is obtained from $A$ by repeating $k_{i}$ times its $i$th column and row.

Interestingly, although the photon-counting probability distribution for input thermal states is given in terms of matrix permanents (namely, of Hermitian positive semi-definite matrix permanents, since $\tau_{i}\geq 0$, $\forall i$) sampling from the probability distribution in Eq. (S13) is classically tractable Rahimi-Keshari et al. (2015). Even more, computing its elements can also be done efficiently for a restricted set of linear-optical transformations $U$ and input states $\rho_{i}^{\text{th}}$ Chakhmakhchyan et al. (2017).

Distinguishable squeezed states: The treatment of GBS with distinguishable squeezed states is analogous to that of standard boson sampling with distinguishable single photons. Namely, since distinguishable squeezed states do not interfere with each other, the corresponding GBS experiment with $k$ distinguishable squeezed states is equivalent to a set of $k$ experiments where a single squeezed state $\ket{\tilde{\xi}_{j}}:=\ket{0,\dots,0,\xi_{j},0,\dots,0}$ evolves according to the $m$-mode linear-optical transformation $U$ ($j=1,\dots,k$). In turn, the photon-counting statistics at its output is obtained by accumulating photon detection events from these experiments. To find the corresponding probability distribution, we first write down the Q-covariance matrix of the evolved states $\ket{\tilde{\xi}_{i}}$:

[TABLE]

where

[TABLE]

Next, in order to find the probability $p_{\text{dss}}(\vec{k})$ of detecting a pattern $\vec{k}=\{k_{1},\dots,k_{m}\}$ of $n$ single photons at the output of $k$ distinguishable squeezed states distributed among $m$ modes, one has to consider all the possible detection events upon the evolved states $\ket{\tilde{\xi}_{j}}$, $\forall j$, that yield the given pattern $\vec{k}$. Since the corresponding expression is rather bulky for arbitrary $n$, $k$ and $m$, we present here the expression relevant to our experimental setting only. That is, we assume that four distinguishable squeezed states are injected in the $a$th, $b$th, $c$th and $d$th mode of a $12$-mode linear-optical circuit. The probability $p_{\text{dss}}^{(4)}(\vec{k})$ of detecting $n=4$ single photons in the $q$th, $r$th, $s$th and $t$th modes (no more than two photons per mode) at its output then reads

[TABLE]

Here, $G$ is the set of all six distinct pairs of input mode numbers, i.e., $G=[\{a,b\},\{a,c\},\{a,d\},\{b,c\},\{b,d\},\{c,d\}]$, $l$ is the list of mode numbers where photons were detected, $l=\{q,r,s,t\}$, $\tilde{{C}}_{i}^{\prime}$ is obtained from the $12\times 12$ matrix $\tilde{C}_{i}={U}\cdot\text{diag}[0,\dots,\tanh\xi_{i},\dots,0]\cdot{U}^{\text{T}}$ by keeping its $q$th, $r$th, $s$th and $t$th rows and columns and $p^{(j)}_{ij}=\left|\tilde{{C}}^{(j)}_{ij}\right|^{2}$. If two photons were detected at a given mode, the corresponding mode number is repeated twice in the list $l$. For instance, if two photons were detected in the $r$th mode, the list $l$ reads $l=\{q,r,r,t\}$. Note also that the first line in Eq. (S18) corresponds to events where a pair of squeezed states produced a detection of two photons, while the second line corresponds to events where a single squeezed state produced a detection of four photons.

We also remark that although the computation of each element of the photon-counting probability distribution $p_{\text{dss}}(\vec{k})$ is hard, sampling from it is classically tractable, analogous to the case of standard boson sampling with distinguishable single-photons (see, e.g., Ref. Aaronson and Arkhipov (2014)).

Two-mode squeezed states: To treat the case where the input states are two-mode squeezed states instead of single-mode ones, we can make use of the analogy between the photonic circuits represented in Fig. S8a and Fig. S8b. Similar connections between sampling protocols using TMS and SMS have already been highlighted in previous works, for example in Ref. Kruse et al. (2018) it was used to relate SBS and GBS. The idea is to note that $m$ two-mode squeezers can be obtained from $2m$ single-mode squeezers combined pairwise into phase shifters and beam-splitters, as shown in Fig. S8b. As photons in different modes of the two-mode squeezers would not interfere in the interferometer, we can separate all the $m$ top output modes of the beam-splitters and the $m$ bottom output modes, and send them into separate but equal interferometers ${T}$. Suppose that single-photon detection at the outputs of both interferometers, obtaining patterns $\vec{h}$ and $\vec{k}$ in the top and bottom $m$ output modes, respectively (note that $\sum_{i}h_{i}=\sum_{i}k_{i}=n/2$ must hold, where $n$ is the total number of photons). This would correspond to a detection patter $\vec{x}=\vec{k}+\vec{h}$ in the original TMS sampler. The output probabilities of the TMS sampler $p_{\text{tms}}(\vec{x})$ can then be calculated from the probabilities of the analogue SMS Gaussian boson sampler $p_{\text{tms}}(\vec{h},\vec{k})$ as

[TABLE]

where $p_{\text{gbs}}(\vec{h},\vec{k})$ is the calculated probability of obtaining patterns $\vec{h}$ and $\vec{k}$ in the top and bottom $m$ output modes of the scheme in Fig. S8b. The probability $p_{\text{gbs}}(\vec{h},\vec{k})$ is obtained using the total transfer matrix $T_{\text{gbs}}=(T\oplus T)\cdot T_{\text{bs}}\cdot T_{\text{ps}}$ which includes the matrices associated to the array of phase shifters $T_{\text{ps}}=U_{\text{ps}}^{\oplus m}$ and beam-splitters $T_{\text{bs}}=U_{\text{bs}}^{\oplus m}$, where

[TABLE]

We remark that, although sampling from TMS input states deviates from an ideal implementation of GBS, it still represents a task which is hard classically, as it can be mapped into an analogous GBS problem. Therefore, rather then a validation on the computational complexity of the sampler, the test performed should here be interpreted as supporting evidence for the correct experimental implementation of the protocol.

Uniform distribution: As for standard boson sampling, the output probabilities for a uniform sampler are simply given by $p_{\text{uni}}(\vec{k})=1/n_{\text{conf}}$ where $n_{\text{conf}}$ is the total number of possible output patterns such that $\sum_{i=1}^{m}k_{i}=n$. For completeness, validation results for GBS data against a uniform sampler are reported in Fig. S9.

Appendix L Gaussian boson sampling with pseudo photon number resolving detection

To study GBS in a non-collision-free regime (i.e. including cases where more than one photon per output mode is possible), we additionally implemented GBS in a configuration where pseudo number resolving photodection was performed, where up to 2 photons could be probabilistically resolved in the output modes. This can be achieved with a slight modification at the detection apparatus. In particular, inserting 50:50 fiber beam splitters on the output fibers and detecting the outputs with the SNSPDs, as pictured in Fig. S10a. In this way, if two photons are present in a single output mode, they can be probabilistically split at the beam splitter with a probability 1/2, and the simultaneous detection of both the associated detectors identifies the presence of two photons on the mode. The cases where a single photon is on a mode are not affected, as it will be detected by either one of the two SNSPDs associated to his mode.

Results for the validation tests performed in the non-collision-free regime are reported in Fig. S10b. As for the collision-free case reported in the main text, the validation are performed against hypothetical models of the experiment where thermal states, coherent states, distinguishable SMS and TMS states are used at the input instead of ideal SMS states, using Bayesian model comparison techniques. In all tests, we obtain high confidence in the ideal model after 125 events.

Appendix M Noises in Gaussian boson sampling with integrated SFWM sources

We briefly discuss here some noise effects that may arise when performing GBS in an integrated implementation based on SFWM waveguide or ring-resonator sources Silverstone et al. (2014, 2015).

The first type of noise we study is the one inserted from losses before the circuit. For example, losses due to cross-two-photon absorption, in addition to transmission losses, are known to limit the heralding efficiency of integrated non-linear sources Husko et al. (2013). While in the standard approach to boson sampling, where single-photon sources are used and the inputs are Fock photon number states, such losses can be simply accounted into a reduced efficiency of the sources Loredo et al. (2017); Wang et al. (2017, 2018b), in GBS the situation is different: losses change the input Gaussian state. In fact, under the action of losses a single-mode squeezed state transforms towards a thermal state, inserting noises in the output distribution. Note that such effect is present also in SBS implementations based on parametric non-linear sources Spagnolo et al. (2014); Carolan et al. (2014); Bentivegna et al. (2015); Carolan et al. (2015), where it is typically referred to as multi-pair emission noise. The GBS formalism allows to describe and study it analytically.

Considering a GBS experiment where $n$ photons are emitted by $k$ SMS sources and injected into the interferometer, the presence of losses before the interferometer can be described as shown in Fig S11a. Losses are modelled by combining each of the $k$ sources with an ancillary vacuum mode through a beam-splitter with transmission $0\leq\eta\leq 1$, and tracing out the (undetected) ancillary modes before sending the signal ones into the interferometer. Considering for simplicity uniform squeezing and losses through the modes, this procedure can be described analytically in the following way. First, we label the odd modes as the signal ones, i.e. with SMS sources on them, and the even ones as the auxiliary modes initialized in the vacuum state, so that the $2k$ squeezing parameters are given by $\vec{\xi}^{\prime}=\{\xi_{1},0,\xi_{2},0,\ldots,\xi_{k-1},0,\xi_{k},0\}$. Following Eq. (S16) and Eq. (S17), the covariance matrix of the state after the action of the losses, modelled by the beam-splitters with transmission $\eta$, is given by:

[TABLE]

where

[TABLE]

The Gaussian state injected in the interferometer is obtained by tracing out the ancillary modes. Its covariance matrix $\sigma$ is therefore simply given by deleting the even rows and columns from $\sigma^{(2k)}$ in Eq. (S21). The output probabilities can finally be calculated via the Hafnian of the matrices obtained as described in Ref. Hamilton et al. (2017); Kruse et al. (2018) and in previous supplementary sections.

In Fig. S11b-c the mean statistical fidelities of the output 4-photon and 6-photon distributions for a 12 modes interferometer are reported as a function of the losses, for various values of squeezing and number of sources used. For a systems of this size, the protocol appears to have a good resilience to this type of noise. As expected, the noise becomes more significant increasing the system size (number of photons and sources) and the squeezing.

The second type of noise we discuss is the noise due two spurious emission of photons via non-degenerate SFWM. This is relevant in our scheme based on SFWM as, using a dual-wavelength pump scheme in order to obtain SMS via degenerate SFWM, non-degenerate SFWM emission can be obtained at each sources by both pumps at different wavelengths. Effectively, each source comprises a single-mode squeezer and also two two-mode squeezers, one per pump. The presence of TMS induces noise in the system. For example, consider the case when we observe two photons in the signal modes (were only the signal wavelength is selected). In the ideal case these would be generated from the first term of a SMS state, but there is also a chance that these two photons could incorrectly originate from two-pairs of non-degenerate pairs emitted via two-mode squeezing, where only the two photons in the signal wavelengths are observed while the other two photons are discarded. Such event would then give an erroneous sample.

The signal-to-noise ratio (SNR) between the correct events originated by ideal SMS and the noisy events where at least a photon pair is originated via spurious TMS can be quantified. For the ideal case, the probability of generating $n$ pairs of photons ($2n$ total signal photons) from $k$ SMS sources with squeezing $\xi$ is given by the negative binomial distribution (see for example Ref. Hamilton et al. (2017); Kruse et al. (2018)):

[TABLE]

The probability of obtaining $n$ signal photons ($n$ pairs) via $k$ TMS sources with squeezing $\xi$ is instead given by (see for example Ref. Hamilton et al. (2017); Kruse et al. (2018)):

[TABLE]

Suppose now performing GBS with $k$ sources based on SFWM, that is $k$ SMS sources and $2k$ TMS sources with same squeezing $\xi$, and suppose we observe a $2n$ photons event. We consider it a correct event if and only if all photons originate from SMS, while all other cases are considered as noise. The probability of observing a noisy event is then given by:

[TABLE]

The SNR is then finally given by:

[TABLE]

The SNR as a function of photons and squeezing is plotted in Fig. S12, where we take $k=2n$ (which saturates the complexity for GBS Hamilton et al. (2017)). It can be observed that for low levels of squeezing the SNR decreases rapidly with the number of photons. However, for higher squeezing the $\cosh^{4k}(\xi)$ factor in Eq. (S28) becomes dominant and a good SNR is recovered.

To reduce this type of noise, different methods can be adopted. First, one may consider to detect the idler photons emitted in the TMS case, instead of discarding them, via filtering them before the interferometer and detecting the idler modes, in the same way as in the scattershot approach. If any idler photon is detected in addition to the photons in the signal modes, then the event is discarded. In this way the SNR can be in principle arbitrarily high (though in practice it would be limited by the loss). One other option would be to consider an hybrid scheme between GBS and SBS, where detecting the idler photons would herald the injection of a non-Gaussian, but known, state into the interferometer, while if no idler is detected then the event is treated in the standard GBS approach. Such scheme would be more suitable for integrated approaches based on SFWM, increasing the count rate by not erasing all the unideal events, and is likely to be of the same complexity of GBS.

Appendix N Franck-Condon profiles

The analogy between photons distributed among optical modes and molecular phonons among vibrational modes comprises a platform for photonic quantum simulations of complex molecular dynamics Huh et al. (2015); Clements et al. (2017); Sparrow et al. (2018). As we detail in this section, based on such an analogy, a setting similar to Gaussian boson sampling can be used, e.g., to simulate an arbitrary molecular Franck-Condon (FC) profile Huh et al. (2015). The latter corresponds to the molecular vibrational transition profile within the harmonic approximation and the assumption of a coordinate-independent electronic transition moment.

According to the Franck-Condon principle Franck and Dymond (1926); Condon (1928), the probability of a vibrational transition from an initial state $|\vec{j}\rangle$, which we assume to be the ground state $|\vec{j}\rangle=\ket{0,\dots,0}$, to a final state $|\vec{k}\rangle=\ket{k_{1},\dots,k_{m}}$ is given as

[TABLE]

where $m$ is the number of vibrational modes and $\mathcal{U}_{\text{Dok}}$ is the Doktorov transformation Doktorov et al. (1977). The latter consists of linear transformations that conserve the number of excitations, mode displacement and single-mode squeezing,

[TABLE]

Here, $\mathcal{D}(\alpha_{i})=\exp\left(\alpha_{i}{a}_{i}+\alpha_{i}^{*}{a}_{i}^{\dagger}\right)$ and $\mathcal{S}(\xi_{i})=\exp\left[\xi_{i}/2\left({a}_{i}^{2}+{a}_{i}^{\dagger^{2}}\right)\right]$ denote, respectively, the displacement and squeezing operators acting on the $i$th vibrational mode ($\alpha_{i}$ and $\xi_{i}$ are displacement and squeezing parameters), $\mathcal{U}_{\mathrm{L}}$ and $\mathcal{U}_{\mathrm{R}}$ are linear transformations and ${a}_{i}^{\dagger}$ $({a}_{i})$ are the bosonic creation (annihilation) operators.

The displacement and squeezing parameters, as well as linear transformations $\mathcal{U}_{\mathrm{L}}$ and $\mathcal{U}_{\mathrm{R}}$ are defined via the initial and final molecular configurations. Namely, these parameters can be found from the following expressions

[TABLE]

Here, $\omega_{i}$ and $\omega^{\prime}_{i}$ are the harmonic angular frequencies of the initial and final electronic states. In turn, $U_{\text{D}}$ and $d_{i}$ are, respectively, the Duschinsky rotation matrix and displacements along the normal coordinates $q_{i}$ (associated with mode operator ${a}_{i}$), which are responsible for molecular structural changes induced by an electronic transition. These structural changes or, in other words, normal coordinate transformations can be written as Duschinsky (1937)

[TABLE]

where $U_{\text{D}}^{[ij]}$ denote matrix elements of $U_{\text{D}}$.

In turn, squeezing parameters $\xi_{i}$ and linear transformations $\mathcal{U}_{\mathrm{L}}$ and $\mathcal{U}_{\mathrm{R}}$ are obtained from the singular value decomposition of the matrix $J$,

[TABLE]

That is, $\mathcal{U}_{\mathrm{L}}$ and $\mathcal{U}_{\mathrm{R}}^{\dagger}$ are the state-space representations of $U_{\mathrm{L}}$ and $U_{\mathrm{R}}^{\text{T}}$ (T denotes matrix transposition), while squeezing parameters $\xi_{i}$ are the natural logarithms of the diagonal elements of the matrix $\Sigma$, i.e., of the singular eigenvalues of $J$.

Having defined all the quantities appearing in Eq. (S29), we can now interpret $p_{\text{FC}}(\vec{k})$ (also known as an FC factor) as the probability of detecting an excitation pattern $\vec{k}$ at the output of a linearly evolved squeezed coherent state $\ket{\Psi}:=\left[\otimes_{j=1}^{m}\mathcal{S}(\xi_{j})\right]\mathcal{U}_{\mathrm{R}}^{\dagger}\left[\otimes_{i=1}^{m}\mathcal{D}(\alpha_{i})\right]\ket{0,\dots,0}$, i.e.,

[TABLE]

Finally, the molecular FC profile at a given frequency $\omega$ reads

[TABLE]

where $\delta_{x,y}$ is the Kronecker delta. Importantly, $\text{FCP}(\omega)$ can be approximated by generating samples $\vec{k}$ according to the probability distribution $p_{\text{FC}}(\vec{k})$. On the other hand, knowing that FC factors $p_{\text{FC}}(\vec{k})$ are the probabilities of detecting excitation patterns $\vec{k}$ at the output of linearly evolved squeezed coherent states, we conclude that FC profiles can be simulated in a photonic experiment. That is, one can approximate $\text{FCP}(\omega)$ by generating photonic displaced squeezed states, evolving them through a linear-optical network and performing boson sampling at its output Huh et al. (2015). Such a simulation procedure is, therefore, similar to Gaussian boson sampling. However, as opposed to Gaussian boson sampling, photonic simulations of FC profiles require data corresponding to the detection of various photon numbers $\sum_{i=1}^{m}k_{i}=n$, i.e., $n$ is not fixed here. Moreover, in general, simulations of FC profiles necessitate displaced squeezed states rather than vacuum squeezed states as per Gaussian boson sampling (Note, however, that molecular structural changes with no displacement do exist. For instance, photonic FC profile simulations of the tropolone molecule, C7H6O2, require squeezed vacuum states only Clements et al. (2017)).

As explained in the main text, the described procedure for simulating FC profiles via a photonic quantum device can be reversed to provide a benchmarking tool. Namely, for a given photonic experiment, an FC profile corresponding to a synthetic molecule can be constructed and simulated. The merit of such a procedure is to understand the fidelities that one could expect to find in a similar setting engineered to investigate a particular molecule. To achieve this, given the description of our experiment, we deduce the corresponding molecular parameters from Eqs. (S29)-(S34) (this procedure yields not a unique but a set of synthetic molecules). Since our experiment relies on squeezed vacuum states, $\ket{\Psi}$ appearing in Eq. (S35) is a squeezed vacuum state as well, $\ket{\Psi}=\otimes_{i=1}^{m}\mathcal{S}(\xi_{i})\ket{0,\dots,0}$. Thus, we are able to simulate squeezing molecular structural transformations only, i.e., $\alpha_{i}=0$, and, consequently, $d_{i}=0$, $\forall i$, in Eq. (S33). In this configuration, the output probabilities can be described using the GBS formalism [Eq. (J)-(S5)], so that the FC factors $p_{\text{FC}}(\vec{k})$ are given by

[TABLE]

where $B_{\vec{k}}$ is obtained from ${B}={\mathcal{U}_{\text{L}}}\cdot\text{diag}[\tanh\xi_{1},\dots,\tanh\xi_{m}]\cdot\mathcal{U}^{\text{T}}_{\text{L}}$ by repeating $k_{i}$ times its $i$th column and row, and $\sigma_{\text{Q}}$ is given in Eq. (J).

To quantify the enhancement of a quantum experiment over classical approximation strategies we adopt the approach developed in Ref. Clements et al. (2017). The idea is to quantify this improvement in terms of the fidelity of experimentally reconstructed profiles with respect to the fidelity obtainable by the best classical approximation. The latter represents FC profiles with highest fidelity achievable from classical experiments, i.e., experiments employing classical states only (that is, states with a non-negative Glauber-Sudarshan function Leonhardt (2010)). As the fidelity between states is invariant under unitary evolutions, finding the optimal classical strategy corresponds to finding the closest classical state to initial displaced squeezed states Clements et al. (2017). In turn, the single-mode classical state closest to a single-mode displaced Gaussian state is a coherent state with the same level of displacement Marian et al. (2002). Consequently, in absence of displacement, vacuum is the optimal classical state. To quantify the quantum enhancement obtained in the experiment we then consider the difference $\mathcal{C}=F_{\text{Q}}-F_{\text{C}}$ between the fidelity $F_{\text{Q}}$ of the FC profile reconstructed from the experiment to the theoretical FC profile, and the fidelity obtained using the optimal classical strategy $F_{\text{C}}$.

As a first test, we reconstruct the FC profile for a synthetic molecule directly associated to the device. Collected experimental data using a pseudo number resolving detection scheme were used. In this case, reversing the protocol and using it as a benchmarking tool corresponds to using the transfer matrix of the device as $\mathcal{U}_{\text{L}}$ and the characterised squeezing values as $\xi_{i}$. As the choice of frequencies $\omega^{\prime}_{i}$ does not affect the fidelity of an FC profile (only its shape), which is the parameter of interest here, we choose them randomly and with arbitrary unit of measure. This defines the synthetic molecule used for benchmarking. In this case, although the fidelity of the profile is very high ($F_{\text{Q}}>99\%$), the quantum enhancement is only $\mathcal{C}=0.4\%$, i.e. while there is still an improvement, it is actually very low. This is due to the fact that the squeezing in the sources is very low, which means that the FC profile is dominated by classical vacuum contributions.

To obtain higher quantum enhancements, molecules involving higher squeezing values have to be involved. Below we describe a prescription to simulate such molecules even when the actual squeezing in the device is lower. This may be useful in practical implementations of Franck-Condon quantum simulations.

For simplicity, we focus on the case where no displacement is used, although the approach can be easily generalised to the case with non-zero displacement. Consider the situation where we want to reconstruct the FC profile $p_{\text{FC}}(\vec{k})$ for a molecule with corresponding squeezing parameters $\xi_{i}$ and unitary operation $\mathcal{U}_{\text{L}}$. However, in our imperfect experiment we have access to a limited amount of squeezing only, $\bar{\xi}_{\text{max}}\leq\max_{i}\xi_{i}$, so that for each source we can only implement a squeezing $\bar{\xi}_{i}\leq\bar{\xi}_{\text{max}}$. The transition amplitudes can be then estimated in the imperfect device by tuning the squeezing parameters in the sources so that $\tanh(\bar{\xi}_{i})=\gamma\tanh(\xi_{i})$, with $0<\gamma\leq 1$ a real constant rescaling and $\text{diag}[\tanh(\bar{\xi}_{1}),\dots,\tanh(\bar{\xi}_{m})]=\gamma\text{\ diag}[\tanh(\xi_{1}),\dots,\tanh(\xi_{m})]$. In this way, considering also uniform losses $\eta$ in the device, from Eq. (S37), we have a simple relation between the ideal FC factors $p_{\text{FC}}(\vec{k})$ and the ones reconstructed in the imperfect device $\bar{p}_{\text{FC}}(\vec{k})$, given by

[TABLE]

with $n=\sum_{i=1}^{m}k_{i}$ the total number of photons in the output configuration $\vec{k}$, and $\mathcal{N}$ is a normalization constant. As $\eta$ and $\gamma$ are known from the characterisation of the device, Eq. (S38) can be inverted to reconstruct the FC profile via

[TABLE]

Note that this post-processing procedure does not increase the computational complexity of the experiment, in the sense that, if from the experiment we are able to collect only up to $n$-photon events, then we are truncating the FC profile in Eq. (S36) at $n$-photon contributions. This implies that an FC profile reconstructed with this approach is faithful only if the ideal contributions involving more than $n$ photons are negligible, which, for large values of $n$, is in general true only for $\eta$ and $\gamma$ close to unity. This can be also observed in Fig. 4 (inset), where the quantum enhancement initially increases with $\gamma$, but starts to decrease as higher order terms become dominant. Therefore, the approach is clearly not scalable, i.e. it will not provide much help in future large-scale implementations. Nonetheless, it can provide significant advantages in near-term quantum devices, as shown in its application to our experiment.

We implemented the protocol described above to reconstruct FC profiles of synthetic molecules requiring higher squeezing using the post-processing procedure. The synthetic molecule to be simulated is obtained as in the first test above, with the only difference that now the associated squeezing values are $\xi_{i}=\bar{\xi}_{i}/\gamma$, where $\bar{\xi}_{i}$ are the squeezing parameters in the device, preliminarily characterised as described in previous sections. In Fig. 4 (inset) the improvement over optimal classical strategies is reported for different values of simulated squeezing (reported as $\max_{i}\xi_{i}$), that is for different values of $\gamma$, both for the case where the post-processing is performed (solid black line) and when raw data is instead applied (dashed black line). It can be noted that the procedure provides a significant improvement, with the quantum enhancement initially increasing and then decreasing when the FC terms, due to more than 4-photon contributions, start to be dominant. The maximum improvement $\mathcal{C}=9\%$ is obtained for $\gamma\approx 5$, for which the FC profile (plotted in Fig. 4) has a fidelity of $86\%$ with the ideal one. In all the cases, applying the post-processing procedure provides a significant advantage as compared to raw data usage.

Appendix O Estimation of photon number scaling with current silicon photonics technology

Finally, we analyse the potential of the integrated quantum photonics approach as a platform for multi-photon boson sampling protocols in terms of estimated count rates when scaling up the component integration with current technologies. By using as relevant parameters for each component the values measured in our experimental set-up, we aim to provide an approximate evaluation of what event rates can be expected from current photonic technologies in near-term devices. The rates for events with $n$ signal photons, emitted from $k$ integrated sources and injected in a $m$ mode interferometer, for SBS and GBS are respectively

[TABLE]

Here $R_{0}$ represents the repetition rate of the pump, $\xi$ the squeezing at the sources, $\eta_{\text{det}}$ the detection efficiency, $\eta_{\text{ch}}$ the transmission of the optical channels between the interferometer and the detectors (which includes chip-to-fiber coupling loss for off-chip detection). The factor $\eta_{\text{u}}^{m}$ represents the total losses in the interferometer, where $\eta_{\text{u}}$ represents the losses in each coupling operation (e.g. evanescently coupled waveguides in a universal scheme Reck et al. (1994); Carolan et al. (2015); Clements et al. (2016)), and the number of single operations each photon undergoes in the $m$ mode interferometer is assumed to be $m$, as in the universal scheme proposed in Ref. Clements et al. (2016).

In the analysis below we consider a source in each input mode of the interferometer ($m=k$), and used the system efficiencies characterised in our silicon photonic device and experimental set-up: $\eta_{\text{det}}=80\%$, $\eta_{\text{ch}}=64\%$, $\eta_{\text{u}}=99.95\%$, $R_{0}=500\text{ MHz}$. We consider two standard types of integrated sources: a low efficiency one, i.e. the spiral sources used in our device, with a squeezing parameter $\xi=0.17$, and more efficient sources based on integrated ring resonator cavities (see, for example, Ref. Silverstone et al. (2015)), with $\xi=0.31$. We also consider the case where on-chip detection is performed, where no losses associated to off-chip coupling and fiber transmission are present ($\eta_{\text{ch}}=1$).

In Table S1 we report the event rate estimates for these cases with a system size corresponding to a 100-modes interferometer ($k=m=100$), which is realistic on near-term devices thanks to the fabrication scalability of silicon photonics. In particular, an array of $100$ detectors to detect the signal output modes (200 in total for SBS considering the idler modes) would be required, which is possible with current technologies. For example, arrays of $>200$ individually addressable on-chip high-efficiency SNSPDs have already been demonstrated (see e.g. Ref. Schuck et al. (2013)). The rates suggest that experiments with tens of interfering photons should be available with further scaling of current components in silicon photonics, without any further major technological breakthrough required.

To investigate the limits of the approach, namely scaling up the number of components with current silicon quantum photonics technology, we have tested the event rates consider larger-scale circuits with up to 1000-modes interferometers and 1000 detectors on the signal modes. Again, we assume the number of sources and signal modes and detectors to be the same, and consider a technology including ring sources and integrated detectors with efficiency as above. As can be observed in Fig. S13a-b, while increasing the number of sources provides an initial improvement to the event rates, it actually turns to be detrimental if the circuit size becomes too large. The reason for this is that losses in the interferometer, which scale as $\eta_{\text{u}}^{m}=\eta_{\text{u}}^{k}$ become dominant for values of $k$ that are too large, suppressing the combinatorial enhancement provided by the scattershot or Gaussian boson sampling approaches. A trade-off thus has to be adopted to obtain an optimal event rate in SBS and GBS, which is reported in Fig. S13d (also shown as black lines in Fig. S13a-b). The optimal event rates as a function of photon number are reported in Fig. S13c for the different approaches. As expected, the presence of losses implies an exponential decrease in the event rate. However, experiments with a large number of photons are still possible before the event rate becomes impractical, which we identify as experiments with an event rate lower than a threshold of 1 event/week (shaded areas in Fig. S13). For GBS such threshold is reached with $\approx 70$ signal photons, while for SBS with $\approx 48$ signal photons. These values are expected to be at the limit of what is tractable for classical supercomputers Neville et al. (2017).

For comparison, we also perform the analysis for another approach which has been recently investigated, that is standard boson sampling using a time-demultiplexed high-efficiency quantum dot source to deliver multi-photon Fock states into the circuit Wang et al. (2017). For the relevant parameters we use values similar to the ones proposed in Ref. Wang et al. (2017): photon generation probability from the quantum dot $p_{\text{qd}}=65\%$ and repetition rate $R_{0,\text{qd}}=76\text{ MHz}$. While we optimistically consider no losses in delivering the photons to the integrated circuit (achievable for example by integrating the quantum dot), the losses in the photon demultiplexing scheme, given by $\eta_{\text{demux}}=\lceil\eta_{\text{switch}}\log_{2}n\rceil^{n}$, where $\eta_{\text{switch}}$ is the value for each switch in the scheme. We consider low-loss switches based on Pockels cells with $\eta_{\text{switch}}=99.5\%$, as proposed in Ref. Wang et al. (2017). The scaling of the event rate in this scenario is reported as a dashed line in Fig. S13c, where, for fair comparison, the size of the interferometer is considered to be the same as for the GBS case. We estimate a significant decrease of the event rates with this approach compared to SBS or GBS.

To further increase the complexity, e.g. to target experiments with hundreds of photons, the capability to scale up the number of current optical components in silicon quantum photonics is not enough: technological progress is required. A first improvement would be to develop materials with lower transmission losses and more efficient sources of squeezed light. Promising integrated photonic platforms are being developed in this direction Dutt et al. (2015). Ultimately, the introduction of rudimentary error-correction techniques in boson sampling, tackling the effect of losses, is likely to be required for applications far beyond current capabilities.