1.

Introduction

Quantum devices and quantum algorithms promise substantial advantages in many computational tasks. Such applications include as notable examples, quantum computing, simulation, communication, and sensing. In recent years, an ever-increasing number of advancements have been made in this field, with the first experiments^1–6 tackling the quantum advantage regime, namely, the scenario where quantum devices are capable of outperforming classical computers in specific tasks. These results have thus opened the way for the development and application of noisy intermediate-scale quantum processors⁷ for quantum-enhanced information processing.

An example of a nonuniversal quantum processor is the Gaussian boson sampling (GBS) scheme. It is a variant of the original proposal of boson sampling (BS), which is a classically hard computational problem that can be tackled through the use of dedicated quantum photonic devices. More precisely, GBS is the problem of generating samples from the photon-counting output distribution of indistinguishable Gaussian states of light after the evolution through a multimode random linear optical interferometer.^8–12 This problem is intractable for a classical computer when the input states are indistinguishable sources of single-mode squeezed states. Then, a dedicated quantum photonic device can tackle such a task more efficiently, thus corresponding to a quantum computational advantage for a problem instance of sufficient size. The GBS problem has thus drawn attention in the photonic community due to the practical chance to achieve the quantum advantage regime with the technology available today. The latest GBS instances have reached the condition where the quantum device has solved the task faster than current state-of-the-art classical strategies in several experiments,^2,4,6 including the “Borealis” machine by the company Xanadu.⁵ Further interest in the GBS scheme is motivated by the connection between the sampling process and the problem of counting perfect matchings of arbitrary graphs. In fact, an important property of GBS is that it is possible to encode an adjacency matrix of a graph in the device by proper tuning of the interferometer parameters and the squeezing values. Thus the collected samples could be used to learn several properties of the graph and be employed to solve some relevant problems in the field such as finding the dense subgraphs and max-clique, simulating vibronic spectra, and graph similarity.¹³ Recently, GBS-based algorithms to solve such well-known problems in graph theory have been formulated.^14–16 First tests on quantum devices have been then performed on small-scale integrated photonics devices¹⁷ and on larger dimensions in bulk optics through time-^18,19 and path-encoded interferometers.²⁰

In recent years, great effort has been devoted to developing more efficient classical algorithms capable of simulating the sampling process,^21–30 with the motivation of investigating the classical simulatability threshold of GBS devices. In this context, further studies had individuated sources of experimental noise that could undermine the hardness of the problem and allow for classical simulations,³¹ such as photon losses^32–37 and photon distinguishability.^38–41 Thus the benchmarking of a GBS experiment can be performed by following a validation approach, i.e., a test that discerns when samples are drawn from classical simulative models. The methods vary from Bayesian tests,^4,42,43 statistical properties of two-point correlation functions^44–46 or higher order ones,⁴ detector binning,⁴⁷ and properties of marginal probabilities,^40,48,49 to a very recent method based on the feature vector components of the graph encoded in the device.⁵⁰

In this work, we focus on testing graph theory-based methods to assess the operation of a Gaussian boson sampler in the presence of experimental noise. As a test bench system for this analysis, we employed the Gaussian boson sampler Borealis,⁵ which recently claimed quantum advantages and has been made available on the cloud on Xanadu Cloud⁵¹ and Amazon Braket.⁵² Some properties of the device have been recently investigated in Ref. 53 by remote users. The authors measured the quadrature coherence scale to find genuine signatures of the features of single-mode bosonic systems in phase space representation. Here, we analyze the capabilities of the method based on graph feature vectors⁵⁰ to exclude that Borealis samples collected in a real experiment subjected to noise can be compatible with some relevant classically simulative modes, that is, thermal, coherent, squashed, or distinguishable particles samplers. First, in Sec. 2, we provide background information on GBS. Then, in Sec. 3, we review the main features of the structure of Borealis. After that, in Sec. 4, we explain the process for interfacing with Borealis. Finally, in Sec. 5, we analyze the collected data and perform the validation of Borealis against the aforementioned alternative models.

2.

Background on Gaussian Boson Sampling and Connection with Graph Theory

GBS is a linear optics scheme to generate samples drawn from the photon-counting distribution generated by Gaussian light sources at the outputs of a multimode interferometer.^10–12 Such a sampling task is hard to simulate for certain classes of Gaussian states, such as indistinguishable sources of single-mode squeezed vacuum (SMSV) states. The hardness is preserved in noisy experimental conditions as long as the levels of losses and photon distinguishability are limited.

Consider $m$ independent sources of Gaussian states ${\rho }_i$ without displacement such that the input is $\rho ={\otimes }_{i=1}^m{\rho }_i$ with a $2m\times 2m$ covariance matrix $\sigma$. The probability distribution of detecting $n$ photons in the configuration $\overrightarrow{n}=(n_1,n_2,\dots ,n_m)$, where $n_i$ is the number of photons in the output $i$ and $\sum _{i=1}^m{n}_i=n$, is

The quantity ${\sigma }_Q$ corresponds to $\sigma +I_{2m}$, being $I_{2m}$ the $2m\times 2m$ identity matrix, and Haf to the hafnian operation of the $2n\times 2n$ submatrix $A_{\overrightarrow{n}}$. Such a submatrix ($A_{\overrightarrow{n}}$) is individuated by taking $n_i$ times the $i$’th row and the $i+m$’th column of $A$, and the hafnian is the operation that counts the number of a perfect matching in a graph represented by a symmetric adjacency matrix.⁵⁴ The whole $2m\times 2m$ matrix $A$ has the following structure:

In this scenario, the adjacency matrix of a graph can be encoded in $B$, through the Takagi–Autonne factorization, by tuning $U$ and $s_i$ values. This provides a direct link between the output of specific instances of GBS and graph theory. For this reason, GBS with SMSV states have attracted a great deal of attention not only for proving a quantum advantage in the sampling task but also for applications on graphs. Indeed, GBS devices have been suggested as possible tools to tackle problems, such as finding densest subgraphs, the max-clique,^14,18 and the graph similarity.^15–17 The opposite scenario with respect to pure SMSV inputs is a thermal sampler, in which there is no coherence among the possible number of emitted photons that result in $A$ matrix with $C\ne 0$ only. Note that the latter case corresponds to a classically simulative model, given that an efficient sampling algorithm from thermal inputs can be defined.

Many approaches have been proposed to validate GBS experiments in the quantum advantage regime.^{42,44–46,50} The goal of a validation algorithm is to exclude the samples drawn from classical simulative distributions, such as the outputs of thermal, coherent, squashed, and distinguishable Gaussian samplers. In this work, we apply a recent method introduced in Ref. 50 that makes use of graph feature vectors estimated from GBS samples. The components of such vectors are called “orbits,” which result from a coarse graining of the output configurations. More precisely, the method consists of the classification of different samplers in the space spanned by the three feature vector components identified by the probability of the orbits $O_1=[\mathrm{1,1},\mathrm{1,1},\dots ]$, $O_2=[\mathrm{2,1},\mathrm{1,1},\dots ]$, and $O_3=[\mathrm{2,2},\mathrm{1,1},\dots ]$ for a given number $n$ of detected photons. An orbit $\overrightarrow{n}$ is defined as the set of possible index permutations of $\overrightarrow{n}$. The three orbits $\{O_1,O_2,O_3\}$, employed for the validation method, collect, respectively, the output states in which the number of photons in the modes is 0 or 1, in which only one detector measures two photons, and in which two modes host two photons. Plotting the feature vectors in the space $\{O_1,O_2,O_3\}$, retrieved by summing the normalized frequencies of the output configuration corresponding to each orbit, can be used to distinguish between various types of classical samplers. Indeed, the orbit distributions for the different models are characterized by different behaviors, such as, for instance, lying on different hyperplanes. We underline that such a choice is effective in the $n\ll m$ regime, in which we expect that bunching configuration has a very small probability to occur. In general, one must choose the orbits that have the highest probability in order to be estimated from a finite sample with good accuracy.

3.

Structure of Borealis as a Time-Bin GBS

Borealis⁵ is based on time-domain multiplexing^55–57 with limited connectivity. A single squeezed-light source emits batches of $m$ time-ordered squeezed-light pulses that interfere with one another with the help of optical delay loops, programmable beam splitters (BSgate), and phase shifters (Rgate). The loops arrangement shown in Fig. 1(a) is an example of a universal time-bin interferometer able to encode any operation over the modes (see also Refs. 55 and 56). The transformation is controlled by time modulations of the splitting ratios of two BSgates [$r_m(t)$ and $r(t)$] and the two phases of the Rgates [${\varphi }_m(t)$ and $\varphi (t)$]. The two loops are concatenated and cover the time separation $\tau$ between consecutive bins and the entire time duration of the $m$ modes, ${\tau }_m=m\tau$, respectively. This architecture has been recently employed to realize a programmable and universal interferometer encoded in the time bin for GBS applications.¹⁹

Fig. 1

Time-bin encoding machines and Borealis structure. (a) Example of universal time-bin interferometers. In such an encoding, the optical modes are discrete time bins. The unitary operation over the modes is performed by two concatenated fiber loops. The shortest loop covers the time separation between two consecutive time bins, whereas the longest one covers the whole duration of the number of modes. The same squeezing source is excited at each pulse, and photon-counting measurements are performed at the end of the evolution. (b) Structure of Borealis.^5,58 The interferometer comprises three consecutive loops of increasing length. (c) Example of a unitary matrix $U_{ij}$ performed by Borealis that shows the limited connectivity among the modes due to the loop structure.

The loop structure of Borealis does not follow the universal layout. Figure 1(b) shows the time-bin interferometer that comprises three consecutive loops with three tunable BSgates $r_{\mathrm{1,2},3}(t)$ and likewise tunable Rgates ${\varphi }_{\mathrm{1,2},3}(t)$. The tunable phase shifters are limited to the range $[-\frac{\pi }{2},\frac{\pi }{2}]$. There are further three static phases ${\tilde{\varphi }}_{\mathrm{1,2},3}$ that cannot be controlled by the user and represent the optical phases of each loop. The time separation between the bins is $\tau =167\mathrm{ns}$. Then the time that covers the evolution of the 216 logical modes plus the 43 ancillary modes of the device is $43.5\mu \mathrm{s}$ and thus coincides with the time needed to obtain one sample from the device. The length of the three loops covers a time delay equal to $\tau$, $6\tau$, and $36\tau$, respectively. The squeezer has four settings for the squeezing parameter, “low,” “medium,” “high,” and “zero.” The squeezing must be set to one of the allowed values, and it cannot be modulated over the temporal modes. The setup ends with the single-photon detection stage. The time bins are translated into 16 different spatial modes by a time-to-space demultiplexer system. Single photons are then detected by 16 photon-number-resolving detectors.⁵

The design of Borealis has been chosen to be a compromise between the need for a programmable device and the requirement to achieve the quantum advantage regime. The architecture does not allow for full connectivity between the modes, with a layout that has been recently shown to reduce the required dimension for simulation via tensor networks,³⁰ while still permitting the cover of high-dimensional spaces with minimized losses and optical resources.¹² For connectivity, we mean the number of modes with which each of the modes interacts directly. In Borealis, for example, the loop structure allows each mode to be connected with at most six-time bins. Thus the three-loop connectivity of Borealis puts some restrictions on the matrices that can be represented, as we can see in Fig. 1(c). The transfer matrix is limited to being zero on the upper half of the diagonal. The moduli of the matrix elements decrease when descending from the diagonal. They increase only every 6 and 36 modes but are still smaller than the values at 6 or 36 modes above. Other restrictions derive from the ranges of tunability of the Rgates that do not allow the encoding of meaningful graphs adjacency matrices inside the device.

4.

Interfacing with Borealis

Borealis was accessed through Amazon Braket,⁵² and the parameters sent were the squeezing values, the parameters of the three-beam splitters, and the parameters of the three-phase shifters for each time mode, as well as the number of shots (i.e., the number of samples). The parameters sent are for 259-time modes, even though the number of outputs that we show is 216. This is because the first 43 modes are used to fill the loops, and thus the output for those first 43 modes is basically background noise. The device thus has 216 “logical” modes, while it has 259 “physical” modes.

When interfacing with Borealis, one also downloads the device certificate, which describes the current calibration of the device. It contains information such as the loop phases, the squeezing values for the various settings, and the various efficiencies. The certificate distinguishes among three types of efficiencies:

The amount of the losses changed on a daily basis. The average number of photons in the output from Borealis was usually around 75% less than that of a lossless simulation with the same parameters. This is due to the slightly different operating conditions of the cloud version of Borealis with respect to the specifications reached in Ref. 5.

The information included in the device certificate is necessary to perform numerical simulations with alternative models, such as thermal samplers, which are used for comparison to assess the operation of Borealis. One aspect to be taken into account regarding the device certificate is that it is measured once a day, right before the 2-h period in which Borealis is available. Thus it could be less representative of the state of Borealis in those runs far away from the calibration. Indeed, the values of these efficiencies seem to oscillate over time, as we will show in Sec. 5. After some discussion with Xanadu’s team, it turned out that these oscillations in the common efficiency are mainly due to variations in temperature. The scale of these variations in efficiency varies daily, with some days being stabler than others.

Another aspect to be taken into account regarding the device certificate is that simulations based on the included parameters, in particular with respect to noise, do not always accurately reproduce some features of the machine output. For instance, in some cases, we found that the difference in the average number of detected photons per shot between the simulations performed according to the certificate parameters and the measured samples was small. This indicated that the loss estimation in the device certificate was accurate. In other cases, the mismatch between simulation and experiment was significantly higher, around 1.5 to 2 photon difference on an average number of photons $\sim 26$, a value not compatible with statistical fluctuations. In these trials, the device certificate did not represent the physical device with sufficient accuracy within the statistical errors of the apparatus. This issue started to be more evident after one maintenance period that was made near the end of January 2023. Such a discrepancy is relevant when trying to compare the experimental samples with simulations of classical samplers in the validation stage of the device, as it is no longer possible to have a truly accurate simulation for those days.

5.

Borealis Data Acquisition and Validation

We will now analyze the results of our runs on Borealis and compare them with some simulations. All of the codes to interface with Borealis, perform the simulations, and analyze the data were written in Python 3⁵⁹ and were based mainly on the libraries Strawberry Fields⁶⁰ and The Walrus.⁶¹

In Fig. 2, we show the distribution of the number of detected photons at the Borealis outputs for the three squeezing levels larger than zero. Our goal is to apply the orbits method to validate our data against classically simulative models, such as thermal, coherent, squashed, and distinguishable SMSV input states. The method of comparing the orbits is effective, even in the quantum advantage regime, as long as the number of modes $m$ is considerably larger than the number of detected photons.⁵⁰ The reason is that in such a condition the probability distribution of the various orbits’ configuration is peaked on a few orbits. This is evident in Fig. 2, which shows the experimental orbits’ distribution for the three squeezing levels. For this reason, all the runs shown in this section were made with the squeezing parameter set to low, which satisfies the condition $m\gg n$. The results of Fig. 2 also motivate the performed choice on the orbits for the application of the validation technique, given that they correspond to highly populated ones.

Fig. 2

Detected photons and orbit probability distributions. The top row: the distributions of the number of detected photons for the three levels of squeezing: (a) “low,” (b) “medium,” and (c) “high.” The bottom row: the distributions of the orbit configurations for the number of detected photons associated with the highest probability, namely, $n=26$ for low squeezing, $n=62$ for medium, and $n=97$ for high squeezing. The $x$ axis represents the various orbits’ configuration ordered as follows, first $[\mathrm{1,1},1,\dots ]$, then $[\mathrm{2,1},1,\dots ]$, $[\mathrm{2,2},1,\dots ]$, etc. The $y$ axis represents the number of samples for that orbit configuration. The size of each sample was 250,000 shots.

To check the stability of the results over time, we ran sampling experiments with the same settings of the circuit for 2 weeks. We executed two runs of 250,000 shots each day from Monday to Friday for these 2 weeks. The parameters we used are the following: all beam splitters are set to a transmissivity of 50%, and all phase shifters are set to 0. We notice that the observed changes in the average number of detected photons in the outputs during the days impacts significantly the probability of orbits $\{O_1,O_2,O_3\}$. This variation is likely mostly due to changes in the common efficiency between days, as mentioned in the previous section. The plot of the three orbits probability for the 2 weeks is shown in Fig. 3, whereas in Table 1, we show the average number of photons per run of the day during the 2 weeks.

Fig. 3

Stability check of Borealis total efficiency over days. Orbit distribution for the 2 weeks of runs with the same circuit settings. The variations in the different days are due to changes in the total efficiency of the apparatus, which is mostly dominated by the variations in common efficiency. The points correspond to the three orbits’ probability for detected photons in the range between 18 and 32 photons. More precisely, each point on the plot is associated with a specific detected photon number. The points at the bottom left corner correspond to 32 photons, while the points on the right are the orbits for 18 postselected photons.

Table 1

Variation in the number of detected photons over daysa

We now discuss the validation of the samples collected from Borealis. Based on the results reported in Ref. 50, we expect that indistinguishable SMSV GBS and indistinguishable thermal samples display orbits on the same hyperplane in the space spanned by $\{O_1,O_2,O_3\}$. Furthermore, in the lossless case, thermal and squeezed light GBS have a similar spread of points in the orbits, while coherent light GBS has a spread of points that is much larger (see also the Supplementary Material). The presence of photon distinguishability results in a change of such a hyperplane for squeezed, thermal, and coherent light. The effect of the photon losses is to change the GBS orbits toward the thermal sampler.⁵⁰ In fact, both balanced (the common efficiency in Borealis) and unbalanced (loop and relative detector efficiencies) losses lead to changes in the orbits but always on the same hyperplane. Indeed, introducing losses to a Gaussian state corresponds to mixing it with the vacuum, thus increasing its thermal component. Hence, this effect does not alter the hyperplane of the orbits. Furthermore, in the presence of unbalanced transmissions between the channels, this introduces also changes to the transfer matrix of the device, which does not move the orbits in different hyperplanes. In Figs. 4(a) and 4(b), we show such behaviors in more detail. The simulation considers thermal light inputs evolving through the transformation performed by the Borealis circuit with different levels of common efficiencies. We set the parameters of the circuit as follows: the transmittance of BSgates was set randomly between 0.4 and 0.6, while the phase shifters were set to random values in their full range $[-\frac{\pi }{2},\frac{\pi }{2}]$. The points that correspond to the same number of detected photons move along lines. More precisely, the orbits somehow expand when the common efficiency decreases and shrink when the efficiency increases. Unbalanced losses, on the other hand, tend to alter the orbits, moving the points away from those lines, but still keeping the orbit on that same hyperplane. Furthermore, from this analysis, other imperfections, such as small errors in the squeezing parameters, beam-splitter ratios, or interferometer phase shifts, are expected to provide only the second-order changes in the orbits. Errors in the squeezing parameter do not alter significantly the class of the Gaussian state, and thus the orbits would move on specific lines of the same hyperplane. Additionally, errors in beam-splitter ratios or phase shifts would correspond to slight changes in the transfer matrix, which do not alter significantly the orbits. We also made an approximate simulation of the orbits of an ideal GBS, as can be seen in Figs. 4(a) and 4(b). From this simulation, we observe that those orbits reside on the same hyperplane as the thermal ones, as expected. However, they do not align with the lines of the thermal sampler.

Fig. 4

Role of losses in the orbit estimation. (a), (b) Orbits for the thermal sampler simulations. All points, except the green and dark gold ones, represent orbits where only the common efficiency has been changed. The lower the efficiency is, the larger the radius of the orbits is. The green points represent the case where one of the detectors was turned off, thus producing an unbalanced loss. The dark gold points are from the lossless indistinguishable SMSV states simulation. The details on the parameters used for the simulations can be found in the Supplementary Material. The lines are obtained by making a linear fit on all points with the same number of detected photons, excluding the points of the orbit with unbalanced losses and the SMSV simulation. Each line represents a different number of photons. We highlighted the number of corresponding photons for some of the lines in this figure. The effect of balanced losses is to move points along the lines. The unbalanced ones move the points away from the lines, but keep them on the same hyperplane. However, such an effect is evident only in the case of a strong amount of imbalance. The lossless SMSV states points are not on those lines but are still on the same hyperplane. We show the projection of the orbits in the plane in (b). (c), (d) Experimental orbits from the runs on Borealis. The lines are the same as the panel above calculated from the simulation with thermal states.

In Figs. 4(c) and 4(d), we compare runs taken in around 3 months, from the end of December 2022 to March 2023, with a break in January for maintenance. The circuit settings were the same that we employed for the numerical simulations with thermal light. The changes in the orbits are compatible with changes in the common efficiency observed during the time of data collection. The other interesting behavior is that the orbits seem to follow closely the lines from the thermal simulations. However, they do to not match exactly those lines. Thus, according to these tests, Borealis appears to behave as a lossy GBS with squeezed states, which suggests that the orbits may be intermediate between a GBS with indistinguishable thermal light and a lossless squeezed light GBS. Furthermore, the orbits still lay on the same hyperplane of the thermal, like the case of a lossy GBS with indistinguishable emitters. This allows us to exclude the notion that the samples derive from a distinguishable particle sampler. Since we have shown that changes in losses do not change the hyperplane of the thermal simulation, from our analysis we can restrict the hypothesis to a scenario where Borealis is either a thermal with losses, a coherent with losses, a squashed with losses, or an SMSV GBS with losses. As a further test, in the following, we will also consider an additional mockup corresponding to the greedy sampler^5,49 that mimics samples with the correct marginals up to a chosen order $k$. More specifically, to further highlight the relevance of losses for Borealis, we considered both samples generated via the greedy algorithm for a lossless device implementing the ideal unitary $U$, with squeezing parameters adapted to lead to the correct average photon number, and a lossy case where the parameters of Borealis device certificate have been used to generate the greedy samples. As explained in Ref. 5, we considered as a relevant case for Borealis the scenario with marginals order $k=2$. We will now proceed to show that it is not a coherent light sampler, and then we will show that it is closer to an SMSV GBS than a thermal or a squashed one and to analyze the greedy sampler scenario.

In Fig. 5(a), we compare data from runs of 250,000 samples from Borealis on different days but with the same average number of detected photons (red points), with simulated samples from coherent (green), squashed (orange), lossless greedy (purple), lossy greedy (cyan), and thermal light (red). In Fig. 5(b), we consider runs performed on the same day, in which the device certificate provided an accurate estimate of the noise in the experiment due to photon losses. The orbits of coherent light are distant from Borealis data in both cases and have a significantly bigger spread than the orbits of Borealis and the thermal simulations. As a side note, we can also see that the orbits of the coherent sampler are distinct from those of a thermal sampler. In fact, even the mean photon number is different between the simulations of coherent and thermal samplers, indicating that losses affect the two classes of states differently, given their different photon statistic. On the other hand, both Borealis and the thermal sampler seem to be affected in a very similar way by the losses, at least in terms of the mean photon number. We also observe that the samples from Borealis can be distinguished from the lossless greedy sampler, while the orbits in the lossy case are closer to the one of Borealis, and thus the discrimination of this mockup needs to be supported by other methods as those used in Ref. 5. This result further highlights the relevance of losses in the analysis of the samples from Borealis. Now, we are left only with the task of showing whether Borealis is closer to a thermal sampler, a squashed sampler, or an SMSV GBS with losses. Squashed states are the closest classical model to the lossy SMSV sampler. Samples drawn from the probability distribution of squashed states are obtained using classical mixtures of coherent states,^43,62 similar to the thermal sampler.⁶³

Fig. 5

Validation of Borealis samples. (a), (b) Comparison of the orbits of Borealis (red), thermal states simulation (blue), squashed states simulation (orange), a second-order greedy sampler with the ideal $U$ transformation and no loss (purple) and with the effective transfer matrix and loss (cyan), and coherent states simulation (green). The coherent, squashed, lossy greedy, and thermal states simulations used the same parameters and device certificate of the Borealis runs. All the points were calculated on a sample size of 250,000 shots with squeezing set to low. (a) Borealis data collected on different days with different $U$ transformations. Such runs display a similar average number of detected photons and device certificates. (b) Borealis data collected on the same day with different unitary transformations. On that day, the fluctuation in the average number of detected photons between runs on Borealis was compatible with statistical fluctuations, and the deviations of the average number of detected photons predicted by the simulations according to the certificate with those measured on the device were small. (c)–(h) Comparison of the covariances of lossy thermal states, lossy squashed states, and Borealis for a given transformation $U$ of the device. The same circuit parameters were used for the simulations and the measurements with Borealis. The same device certificate was used for all the simulations. (c) Covariance of the simulated thermal sampler with losses. (d) Covariance of the thermal sampler with losses calculated from 250,000 samples, to reproduce the additional noise due to a limited number of samples. (e) Simulated covariance of a squashed light sampler with the losses of the Borealis circuit. (f) Covariance of squashed light sampler from a finite sample of 250,000 shots. (g) Simulated covariance of an SMSV GBS with losses. (h) Covariance of Borealis samples.

As we can see in Figs. 5(a) and 5(b), the orbits of Borealis and the thermal and squashed simulations are similar. In particular, the distances between the orbits of the two classical models corresponding to a given number of detected photons are compatible with the ones from Borealis within the typical orbits’ dispersion due to the intrinsic fluctuations of Borealis parameters in time. This means that in this case, the orbits alone are not enough to ascertain whether Borealis is more likely to be a GBS with SMSV input states, a squashed sampler, or a thermal sampler. To further distinguish among these three cases, we look also at their respective sample covariance $\mathrm{Cov}(n_i,n_j)$. These quantities form the matrix of the two-point correlators in the output modes $i,j$, with matrix elements $C_{ij}$ defined as $C_{ij}=\mathrm{Cov}(n_i,n_j)=⟨n_i{n}_j⟩-⟨n_i⟩⟨n_j⟩$. These quantities are informative to help in discerning the nature of the sampler as shown in the theoretical works^44,45 and previous experiments.^2–6,46 Figures 5(c)–5(h) show the simulated covariance matrices and the one calculated from a finite sample of 250,000 shots. From this analysis, the behavior of Borealis is closer to the one of an SMSV GBS with losses, rather than a lossy thermal sampler or a lossy squashed sampler. This conclusion is supported by performing a more quantitative analysis of the two-point correlators. In Fig. 6, we compare the values $C_{ij}$ of samples collected from Borealis with the expected values calculated with the different possible models. We observe in Figs. 6(a) and 6(b), where the data are compared with a thermal and a squashed sampler, that these models do not provide an accurate explanation of the collected samples, as confirmed by the slopes of a linear fit corresponding, respectively, to $\sim 3.76$ (thermal) and $\sim 2.28$ (squashed). Conversely, the comparison of the correlators with the expectations from a lossy SMSV GBS is characterized by a better agreement between data and model, as also quantified by the fitted value of the slope, equal to $\sim 0.94$ and thus close to the ideal value of 1.

Fig. 6

Two-point correlators. Scatter plot of two-point correlators $C_{ij}$ of samples collected from Borealis (experiment) as a function of the corresponding values calculated with different possible models for the device (simulation). Comparison with a (a) thermal sampler, (b) squashed state sampler, and (c) lossy SMSV state GBS. Black dashed line: the expected trend for the ground truth, corresponding to a linear function with a slope equal to 1. Colored dashed lines: trends obtained by a linear fit of the data shown in each figure.

These final observations allow us to conclude that, among the hypotheses tested in the above analysis, Borealis is likely to be performing a genuine SMSV GBS with losses. For completeness, we also analyzed some of the data from Ref. 5, using the method based on graph feature vectors. The results are shown in Fig. 7, where we show in red the experimental orbits for the data corresponding to the lowest squeezing value of Ref. 5. The simulated thermal (blue) and squashed (orange) samples were generated according to the circuit parameters of that time. It is worth noting that both the squeezing values and losses are smaller than the latest run with Borealis shown above. The orbits of the three models seem to be slightly more separated in this set of experimental data, in accordance with the expectation of a GBS device with reduced losses. The device certificates for this case can be found in the Supplementary Material.

Fig. 7

Orbits of the data from the Borealis experiment of Ref. 5. Red, the Borealis data; blue and orange, the simulated data from the thermal and squashed samplers, respectively. The size of the sample was $\sim {10}^6$. The points correspond to the number of detected photons from 18 (from the right) to 32 (to the left). The Borealis data are those of Fig. 3(a) of Ref. 5. The lines are those we obtained from the thermal simulations in Fig. 4.

6.

Conclusions

In this work, we have analyzed the fundamental issue of assessing the operation of noisy intermediate scale quantum devices. More specifically, we have focused on the problem of GBS in a noisy scenario that represents one of the most investigated approaches to demonstrate the achievement of the quantum computational advantage regime. Here, we analyzed the capability of a graph theory-based method to identify and exclude the main noise source in a Gaussian boson sampler. As the test platform, we have assessed the operation of Borealis, a large-scale photonic device that has been made available on the cloud. In particular, we analyzed its working principles and its performance as a sampling machine of indistinguishable SMSV states.

To this end, we first compared the feature vector components of the graph encoded in the device inspired by Ref. 50, namely, the orbits of Borealis, with those of the lossy thermal sampler and a lossless SMSV GBS. This method is effective even in the quantum advantage regime as long as the number of detected photons is much smaller than the number of modes. In such a condition, satisfied in Borealis only for the low level of the squeezing parameters, we investigated the effects of photon losses in the orbits’ estimation. In particular, the main source of noise in Borealis was a result of the amount of balanced loss. Their effect is to move the orbits on the same hyperplane individuated by the data of indistinguishable thermal, squeezed, and coherent light emitters. This allowed us to exclude with high confidence the presence of significant photon distinguishability in the device. The distribution on such a hyperplane of Borealis data shows significant differences also with the points generated by coherent light. The orbits of the simulated thermal and squashed light are still close to the Borealis data instead. The small deviations observed could be compatible with the discrepancy between the device certificate employed to set the parameters of the simulations and the actual experimental conditions. Additional analysis has been also discussed regarding the greedy algorithm mockup, simulating samples having the correct marginals up to a chosen order. As a further benchmarking of the device, we evaluated the two-mode correlation functions^44,45 summarized in the covariance matrix of samples. This last analysis showed a significant deviation from a lossy thermal sampler and a less pronounced, although still present, deviation from the squashed sampler, and thus shows that, among the class of hypotheses tested by the employed approach, Borealis is likely to behave as a lossy SMSV GBS.

On the one hand, our analysis showed the effectiveness of the orbits’ method in noisy conditions and, in particular, its power in highlighting the effect of the various contributions of photon losses in GBS experiments. On the other hand, the observation that the orbits generated by data collected from Borealis and by thermal and squashed samplers underline possible limitations in the use of the device for graph-related problems. More precisely, our results showed that the collected feature vector components can be approximated by a classical simulation with thermal light, and using squashed light the simulation can be even closer to the ones obtained from Borealis. Other recent works open questions regarding the advantage of the use of a noisy GBS device against simulation with classical simulative models for the problem of graph max-clique and densest subgraphs.^64,65 Due to the reduced connectivity of the circuits and the limited range of the phase shifters achievable via the Borealis device, we were not able to find a way to encode a meaningful graph in the device, namely, an adjacency matrix with a submatrix of significant size that is completely nonnegative (or nonpositive), to test the machine in this context. Furthermore, the structure of the device does not permit the encoding of graphs belonging to different classes of isomorphism, thus preventing the use of the graph’s feature vectors, associated with the measured orbits, for graph similarity and isomorphism applications. Future perspectives of our analysis regard a systematic study of noisy GBS devices, such as Borealis for these graph-related problems, and an investigation on whether the method for certification employed here can be extended to the original version of BS, in light of recent studies of its connection with graph theory.⁶⁶