Timestamp boson sampling

Quantum computing, based on the superposition, entanglement, and interference principles of quantum states, has the potential capability to outperform any existed classical computers for certain problems, attracting extensive attention to develop quantum algorithms^1,2 and physically construct quantum computers.³ However, the noise in quantum gates and the decoherence in quantum systems limit the size of the quantum circuits and the number of the controllable qubits, which make the realization of universal quantum computers still challenging and considered as a long-term goal.⁴ As an intermediate but crucial milestone for benchmarking the development of quantum computing is the concept of quantum advantage,^5,6 which means that a quantum machine can perform a specific task demonstrating computational power beyond the state-of-art classical computers. Much effort has been devoted to construct quantum systems^7–10 with a state space large enough to meet the required size of the quantum machine.

Boson sampling, first proposed by Aaronson and Arkhipov,¹¹ is to sample the output probability distribution of n identical bosons scattering in an m-mode (m > n) interferometer, which can be expressed by a Haar random matrix U. The computational complexity of boson sampling is quantified by the calculation of the submatrix permanent of U, which can be further enhanced by the combination factor $(mn)$⁠.^12,13 From the experimental point of view, in the photonic scheme, only single photon sources, linear-optical network and single-photon detection, are needed to realize boson sampling. Due to the theoretically well-defined and experimentally friendly feature, boson sampling is considered to be a short cut to achieve quantum advantage.

Up to now, boson sampling has been realized in many scenarios,^9,14–23 and variant protocols such as scattershot boson sampling^24,25 and Gaussian boson sampling^26–29 have also been conducted. Simulations on supercomputer have been shown that 50 photons or an equivalent size quantum system may beat the current calculation limit, so as to achieve quantum advantage.^30,31 However, most of the experiments are far away from accessing such a system scale, due to the probabilistic nature of single photon sources and inevitable loss in linear-optical network.

Statistical errors caused by Poisson distribution are magnified in the extremely low-flux rate, which is an essential restriction not only for boson sampling but also for all the count-rate-limited experiments. In general, we require hundreds of counts to reduce this statistical error, but in reality, to achieve quantum advantage, output combinations $(mn)$ of boson sampler must be very large, leading to a very low-flux rate for every combination. Thus, most combinations may only be detected a few times, even only once, which prevents a standard boson sampling from realizing in such a situation.

To solve the above problem, we propose a timestamp boson sampling scheme that can effectively reduce the experimental time consumption in any problem size. The experimental procedure has been shown in Fig. 1. Identical photons are injected into the linear-optical network via a polarization-maintaining (PM) V-groove fiber array. The outputs of the network are coupled to a multimode V-groove fiber array and then connected to the avalanche photodetectors (APDs). All the electronic signals are sent to the time-of-flight storage (ToFS) module, which can record all the occurrences of photons at each outputs channel. By setting a reasonable delay between trigger and signal channels, we can extract the simultaneous occurrence of photons at different output channels. Then, we measure the registration time of consecutive occurrences, instead of measured number of occurrences.

We retrieve the registration time τ_i of each combination as the timestamp and calculate the probability based on the photon counting statistics (see  Appendix A for details). We assume that the Fock states generated by SPDC (spontaneous parametric down-conversion) source approximately satisfy the Poisson distribution, then the time interval between two consecutive events satisfy exponential distribution. Based on the relation between the Poisson parameter λ and the mean time τ that $λ=τ−1$⁠, the reconstruction probability p_i of each combination can be calculated by Eq. (1) and the output distribution satisfies the normalization $∑ipi=1$ (see  Appendix A for details).

For recording the complete time information of all the individual modes, we develop a time-of-flight storage (ToFS) technique (see  Appendix B for details), which has been shown in Fig. 1(b). The ToFS is mainly composed of the trigger channel and the signal channels. The delay of 30 signal channels can be set independently and the time information of all 31 channels is recorded simultaneously. The time accuracy of individual channel reaches the picosecond level and such a high time resolution guarantees the precision of the registration time τ_i. We extract the fourfold coincidence events within the 2 ns coincidence time window and use the time information τ_i of the trigger as the timestamp to mark the occurrence of each combination.

The linear-optical network is realized by 3D photonic chip fabricated with femtosecond laser direct writing technique^32–36 which has both advantages and limitations (see  Appendix C for fabrication details). The advantage depends on its 3D fabrication capability which can realize large-scale chips with thousands of modes.³⁴ In addition, adding a spatial dimension can more effectively expand the matrix size and experiment scale, which is suitable for multiphoton experiments, such as quantum walks³⁷ and boson sampling. The limitation is that the waveguide bending transformation brings additional losses to the chip and the crosstalk problem induced by phase thermal control. Therefore, the large-scale arbitrary controllable linear networks for the multiphoton experiment is still a challenge. The waveguide structure inside the chip is shown in Fig. 2(a), where six input waveguides are coupled into a 2D waveguide array and then the 3D 5 × 6 structure transforms into a tiled 30 output waveguides array after evolution. The spacing of all the waveguides for input and output is set to 127 μm in order to match the external V-groove structure. To introduce random phase shift,^19,38 we deform the S-bend randomly of each input and output waveguide.

We use single photon scattering and a series of HOM (Hong–Ou–Mandel) interferences to reconstruct the unitary matrix U. The measured results of the amplitude elements and the phase elements are shown in Figs. 2(b) and 2(c) (see  Appendix D for characterization details). In Fig. 2(b), to obtain all the moduli of the matrix, we inject heralded single photons into all six inputs one by one and measure all 30 output modes simultaneously. In Fig. 2(c), the phase relationships between different output modes are achieved by measuring hundreds of HOM interferences using the selected input ports.³⁹ In the reconstruction algorithm, we fix the phase part in the first column and last row as 0, and some other channels have phase value around $π/2$ or $−π/2$⁠. In Fig. 2(c), some values of phases are cluster, but have no obvious periodicity, so it will not affect the results of the experiment.

Aside from the fabrication and full characterization of the 3D photonic chip, we build a six-photon source, which is depicted in detail in Fig. 3(a). The multiphoton states are based on spontaneous parametric downconversion^40–42 (see  Appendix E for details). The generated photons are coupled into PM optical fibers connected to the PM V-groove fiber array. We achieve the temporal overlap of non-classical interference among different down-converted photons by translating linear motorized stages to compensate external delays.

We report the experimental results in two different input scenarios. One is a four-photon experiment with three-photon injection and one external trigger, whose sampling rate is sufficient so that we can detect almost all the output combinations within a moderate accumulation time. The other is a six-photon experiment with three heralded single photons injection, and we can only detect parts of the output combinations due to the low-flux rate. It is worth mentioning that in the six-photon case, we synchronize the coincidence of three trigger signals into one signal with a field programmable gate array (FPGA) module, so that we regard the synthesized triggers as timestamps, but in essence, they are sixfold coincidence events.

We sample the four-photon experiment with both timestamp protocol and counting protocol, respectively. We collect a total of 50 000 s data when the pump power is set as 600 mW. For the counting protocol, we observe 83 455 fourfold coincidence events (4021 combinations) in total with a sampling rate about 6000 per hour. Then, we divide the total time into 6000 intervals and count the event amounts in each time interval. The results shown in Fig. 3(b) are approximately Poisson distribution, which verify the rationality of our theoretical assumptions.

For the timestamp protocol, we measure the average time of each combination occurring five times to reduce the influence caused by the fluctuation of registration time and then 3615 combinations left. We calculate the probability distribution using both counting and timestamp protocol in the subspace (see  Appendix F for processing details). We plot the timestamp reconstruction results (p_i) and counting reconstruction results (c_i) in a five-layer ring,⁴³ as shown in Fig. 3(c). The similarity between the two reconstruction protocols $S=∑ipici$ is up to 98.7%. The total variation distance $D=(1/2)∑i|pi−ci|$ and the value is 0.13. These two reconstruction results are highly similar (see more results with different occurrence number in  Appendix G, Fig. 9), which confirm that the timestamp protocol works properly. Moreover, in our timestamp protocol, each combination only needs to happen once if the registration time accurate enough; therefore, ideally, we can effectively reduce the execution time and speed up the sampling process, indicating a great resource-effective advantage.

When it is infeasible to detect all the combinations in large-scale standard boson sampling, we should still be able to show a valid validation^44–46 even with a small number of samples. We also sample the six-photon experiment with both counting protocol and timestamp protocol. For the counting protocol, after a collection time of nearly 100 h, we obtain 358 sixfold coincidence events in total, where 321 combinations just happen once, 17 combinations happen twice, and only one combination happens three times. By using the data processing in  Appendix F, only 292 combinations left and each combination occurs once. The calculated probability distribution is equivalent to a uniform distribution and the probability is 1/292. The results cannot pass the row-norm estimator (RNE) test which can distinguish the boson sampler from uniform sampler. The probabilistic of registration events lead to statistical errors in standard boson sampling protocol, especially when the sampling rate is extremely low.

Instead, we use timestamp protocol and validate the reconstruction results^44–48 when each of the 292 combinations occurs once. The process of validation is no different from standard boson sampling. To distinguish our experimental results from the uniform and distinguishable samplers, we perform the row-norm estimator test⁴⁵ and the likelihood ratio test⁴⁶ (see  Appendix H for details). The results are shown in Figs. 3(d) and 3(e). The different growth trend between the experimental data (yellow dots) and the simulated data (red dots) indicates that the reconstruction results are not extracted from a uniform distribution or sampled by distinguishable photons. Timestamp protocol can still work properly even in the extremely low-flux and low-sampling rate, and can be very useful for large-scale validation against classical imposters.

We further mine more information about the nature of the timestamp protocol. In Fig. 4(a), we extract the time interval of both 83 455 events and, as an example, the output configuration (11 022) in the four-photon experiment. The time interval between successive occurrences of events satisfying exponential decay distributions, which also confirm our theoretical assumption.

Based on our initial assumption that the Fock states generated by SPDC source approximately satisfy the Poisson distribution, the registration time measured in the experiment always fluctuates. For getting rid of the influence of time fluctuations, we use the procedure in  Appendix F to process the results, which inevitably discards some data. Therefore, it is important to find a suitable sampling region between the number of occurrences and the number of singular points, see the results in Fig. 4(b). We perform statistics on the data from the four-photon experiment. When each combination only occurs once, a large number of singular points will be discarded. Until the number of occurrences is increased from 5 to 10, the proportion of the discarded singular points drops from 9.96% to 1.77%. We can obtain a good trade-off between the success rate and the execution time in timestamp boson sampling experiment.

Figures 4(c)–4(h) further show the comparison between the experimental (p_i) and the theoretical (t_i) results with the timestamp protocol in the four-photon experiment (the corresponding validation results can be found in  Appendix G, Fig. 10). We use fidelity and total variation distance to quantify the performance of our experimental result, which are defined by $F=∑ipiti$ and $D=(1/2)∑i|pi−ti|$⁠. For a perfect boson sampler, F should be 1 and D should be 0. By adding the number of occurrences one by one, we can observe the fidelities and total variation distance slightly improve, but deviate from the theoretical predictions to some extent. The deviation is mainly due to the differential losses among different channels. The 5 × 6 3D structure is transformed into 1D output layout. Different waveguides have different geometric track and path lengths. Such losses imply the matrix will not be unitary¹⁶ and change the final photon distributions, so deviated from ideal distributions. By using site-by-site addressing technology³⁷ and imaging-based large-scale correlation measurement,⁴⁹ we are already trying to avoid these additional losses and improve this problem.

In summary, we propose and experimentally demonstrate a timestamp boson sampling based on the retrieval of the registration time information of sampling events. The 3D photonic chip provides an integrated and scalable way for multiphoton experiments. The deviation of experiment results from ideal distributions mainly comes from the differential losses caused by the bending transformation of the waveguides. The extra losses also make the large-scale linear networks be challenging. An efficient solution is to avoid unnecessary bending transformation with site-by-site addressing technology or imaging-based large-scale correlation measurement.

The results in the four-photon experiment show that the timestamp protocol works properly and effectively reduce the execution time, which can improve the scale and state space of the experiment (see  Appendix I for details). In the six-photon experiment, the validation results verify that the experimental events sampling by timestamp protocol are still from a genuine boson sampler even with limited datasets. In addition, we give a reasonable region between the number of occurrences and the number of discarded data. By using higher-efficiency and higher-purity quantum source such as quantum dots, the time fluctuations and discarded data can be reduced. Moreover, we show that it is possible to obtain faithful results from very limited events by retrieving their time information. Such a mechanism can be generalized to a wide range of quantum optics, such as multiphoton quantum walk,⁵⁰ noisy intermediate-scale quantum (NISQ) technology,⁴ and multiphoton entanglement,⁸ driving our capacity of controlling quantum systems to an entirely new accessible realm.

This research is supported by the National Key R&D Program of China (2019YFA0308700, 2019YFA0706302, and 2017YFA0303700); National Natural Science Foundation of China (NSFC) (11904229, 61734005, 11761141014, 11690033); Science and Technology Commission of Shanghai Municipality (STCSM) (20JC1416300, 2019SHZDZX01); Shanghai Municipal Education Commission (SMEC) (2017-01-07-00-02-E00049); China Postdoctoral Science Foundation (2020M671091). X.-M.J. acknowledges additional support from a Shanghai talent program and support from Zhiyuan Innovative Research Center of Shanghai Jiao Tong University.

The authors have no conflicts to disclose.

Wen-Hao Zhou: Data curation (lead); Formal analysis (lead); Investigation (lead); Methodology (lead); Project administration (lead); Software (lead); Validation (lead); Visualization (lead); Writing – original draft (lead); Writing – review & editing (lead). Xian-Min Jin: Conceptualization (lead); Formal analysis (lead); Funding acquisition (lead); Investigation (lead); Methodology (lead); Project administration (lead); Supervision (lead); Writing – original draft (lead); Writing – review & editing (lead). Jun Gao: Data curation (equal); Methodology (equal); Software (equal); Validation (equal); Visualization (equal); Writing – original draft (equal). Zhi-Qiang Jiao: Formal analysis (equal); Methodology (equal). Xiao-Wei Wang: Data curation (equal); Formal analysis (equal); Software (equal). Ruo-Jing Ren: Data curation (equal); Formal analysis (equal); Methodology (equal). Xiao-Ling Pang: Software (supporting). Lu-Feng Qiao: Software (supporting); Visualization (equal). Chao-Ni Zhang: Software (supporting). Tian-Huai Yang: Software (supporting).

The data that support the findings of this study are available from the corresponding author upon reasonable request.

The time-of-flight measurement and the single photon detection have been researched in single-photon imaging.^51–53 The record of the detection times can be described by the Poisson process. We expand the timestamp method to multiphoton experiment, but the experimental aim and result are different. In the single-photon imaging, we are concerned about the number of photons detected on each pixel. In the boson sampling experiment and other multiphoton experiments, we use timestamp protocol to detect and record the coincidence events.

We use the standard statistical techniques to explain the timestamp protocol. The photon statistics of the coherent laser is Poissonian statistics, and the photon statistics of the squeezed light is sub-Poissonian statistics.⁵⁴ 

In our experiment, we assume that the Fock states generated by SPDC source approximately satisfy the Poisson process. The Poisson process can be described by Poisson distribution which is used to show the number of times in an interval of time or space as shown in Fig. 5(a). Each blue ball means that a certain combination is detected once, and we assume that N counts are detected after time T. We can divide time T into n intervals and the time interval is t₁. For the convenience of description, we define t₁ as 1 s. Based on the fundamental assumptions above, the number of events detected per second satisfies the Poisson distribution $P(x=k)=λkk!e−λ,k=0,1…$⁠. The mathematical expectation $E=∑k=0∞kP(k)=λ$⁠, where λ represents the average number of events detected per second. We can also calculate the variance $D=∑k=0∞k2P(k)=λ$⁠, so the error bar of each measurement is $λ$⁠.

On the other hand, a Poisson process is a renewal process in which the interarrival intervals have an exponential distribution function.⁵⁵ Therefore, we can revisit this measurement from the timestamp perspective as shown in Fig. 5(b). We detect the execution time t₁, t₂, $t3…$ and obtain the time intervals t₁, $t2−t1,t3−t2⋯$⁠. Statistically, these time intervals satisfy the exponential distribution. The probability density $f(x)=λe−λx,x≥0$ and the distribution $F(k)=1−e−λk,k≥0$⁠. Therefore, the mathematical expectation $E=∑x=0∞xf(x)=1/λ$ and the variance $D=∑x=0∞x2f(x)=1/λ2$⁠, where $1/λ$ represents how long the certain combination takes to happen once.

Obviously, counting protocol is inversely proportional to timestamp protocol according to their mathematical expectation. We know the counting probability distribution p_i is defined as $pi=λi/∑iλi$⁠. Corresponding to the timestamp protocol, we change the λ_i as $τi−1$⁠, which is Eq. (1) in the paper.

We use Fig. 6 to define the τ_i. As shown in Fig. 6, the blue parts represent the trigger channel and red parts represent the signal channels. We define $τi(j)$ as the time passed between the beginning of the experiment and the coincidence detection of a certain combination j. Since there are delays among different channels, we uniformly define the trigger registration time as $τi(j)$⁠. When a combination occurs twice, the registration time are $τ1(j)$ and $τ2(j)$⁠. The time interval is $τ1(j)$ and $τ2(j)$ − $τ1(j)$⁠. If we extract the time information of each combination occurs twice to reconstruct the probability distribution, there is no difference between using absolute time and time interval because we take the average time information $τ2(j)/2=(τ1(j)+(τ2(j)−τ1(j)))/2$⁠.

The main part of the ToFS has been shown in Fig. 1. The output part of the photonic chip is connected to 30 avalanche photodetectors (APDs) with a multimode V-groove fiber array. The detector signals are linked to the 30 channels time-of-flight storage module. We set the delay of each signal channel relative to the trigger channel, respectively. We scan the delay between different signal channels and the trigger, which identify the differential delays by maximizing their coincidence counts. The results have been shown in Fig. 7. The coincidence time window is set as 2 ns, and then, we get fourfold coincidence events through post-processing. In order to mark each fourfold coincidence, we use the time information τ_i of each trigger signal as the timestamp.

The 3D photonic chip is fabricated by a femtosecond laser with the repetition frequency of 1 MHz, the pulse duration of 290 fs, and the wavelength of 513 nm. After being reshaped by a cylindrical lens, the writing laser is focused into a borosilicate glass substrate with a ×50 objective lens (NA = 0.55). The sample is placed on a high-precision three-axis translation stage moving at a velocity of 15 mm/s. We obtain a photonic chip with 1D six-input fan-in, 30-output fan-out, and 3D 5 × 6 random coupling structure. All the inputs and outputs are set at 170 μm underneath the surface of the chip substrate. In order to introduce the randomness to the scattering matrix, the six input waveguides have different S-bend structures (geometric tracks) during the transformation process. Therefore, the coupling length and coupling strength between these waveguides are not uniform. This further brings randomness into the evolution. The input waveguides are coupled to a polarization-maintaining (PM) V-groove fiber array and the output waveguides are coupled to a multimode V-groove fiber array.

To estimate the efficiency of the chip, we couple the coherent laser (780 nm semiconductor laser) into the six inputs one by one and measure the output intensity distribution. The chip efficiencies of the three chosen input 2, 4, 5 are 26.8%, 30.5%, and 32.5%, respectively. Due to the different bending transformation of different paths, the efficiency of each input is also different. The average chip efficiency of all six inputs is about 25.3%. The length of the 3D chip is 5 cm, so the losses of the chip is about 1.2 dB/cm. Considering the losses of end-face coupling, the actual loss of the chip itself will be smaller.

The scalability of the experiment depends on the efficiency of the system, including the quantum source, 3D photonic chip, and detection. For the 3D photonic chip, femtosecond laser direct writing technique provides a promising platform. The 3D fabrication capability can easily expand the matrix size of the linear networks. However, additional losses induced by waveguide bending transformation limit the sampling rate in the multiphoton experiments. To improve this problem, we can use 2D site-by-site addressing technology to reduce the bending losses. Furthermore, we can also use single photon sensitive camera to directly measure the coincidence in the future. In addition, we can also improve the brightness of the quantum source to expand the scalability by using quantum dots or integrated quantum source.⁵⁶ 

The unitary scattering process of the 3D photonic chip can be described by the reconstruction of a unitary matrix U, and each element in the matrix is a complex number. In order to determine the moduli of the elements, we inject heralded single photons and measure the intensity of each output mode, and then we can easily obtain these values. Considering the chip efficiency and randomness, we choose port 2, 4, and 5 as the input ports. The determination of the argument can be divided by two steps: first, we determine the absolute value by measuring the visibility of two-photon HOM interference, second, we determine the sign of argument by scanning and comparing a large number of HOM interferences between different channels. All the HOM interferences curves have been shown in Fig. 8. Finally, a 3 × 30 matrix can be reconstructed by using the characterization data obtained above.

The setup of six-photon source has been shown in Fig. 3(a), a mode-locked Ti:sapphire laser emits a femtosecond pulse (140 fs) with the center wavelength of 780 nm and then is frequency-doubled by a lithium triborate (LBO) crystal. The generated 390 nm laser is split into two spatial modes by a balanced UV beam splitter. The reflected laser passes through two beta barium borate (BBO) crystals and generates two pairs of correlated photons via spontaneous parametric downconversion. The BBO1 crystal is type-II phase matching in a beam-like scheme, while the BBO2 crystal is type-II collinear phase-matched configuration. The transmission laser passes through the BBO3 crystal which is the same as the BBO1. The 3 nm bandpass filters (BPFs) are employed to guarantee the spectrum of each path approximately the same. The half waveplate (HWP) and quarter waveplate (QWP) are used to initialize the polarization. At last, photons are coupled into PM optical fibers connected to the PM V-groove fiber array. We achieve the temporal overlap of non-classical interference among different down-converted photons by translating linear motorized stages to compensate external delays.

Taking the four-photon experiment as an example, we illustrate our reconstruction steps when the number of occurrences equals to 1. Statistics show that 4021 out of combinatorial $(30,3)=4060$ happen at least once.

Step I. We reconstruct the probability distribution of all 4021 combinations by using their time information τ denoted by the detection time of the trigger.

Step II. We calculate the theoretical occurrence number $Nτ$⁠, which is defined as $Nτ=T/τ$⁠. T represents the total time 50 000 s and τ represents the timestamp of each combination. We also count the actual occurrence number N_a of each combination in total time.

Step III. Based on our assumption in  Appendix A, the mathematical expectation of the Poisson distribution is λ, which represents the average number of events detected per second N_a. The error bar of each measurement is $Na$⁠. Therefore, we can assume N_a fluctuates but most often stays in the range of $[Na−Na,Na+Na]$⁠. For low-flux rate scenario, we use a more tight but appropriate bound to reshape our data and only keep the events in the range $Na/2<Nτ<2Na$⁠. We normalize the data in the subspace and compare it with the distribution obtained from theory and standard boson sampling.

We process the four-photon and six-photon experimental data with the same procedure. The only difference is that when the number of occurrences increases, we replace the single timestamp with the average timestamp of multiple events, which makes the data processing more precise.

We show more experimental results in this section, including the experimental comparison of the timestamp reconstruction and the counting reconstruction with different n_o (the number of occurrences) in Fig. 9, and the validation which is the supplement to Figs. 4(c)–4(h) in Fig. 10.

To distinguish from a uniform sampler, we use the row-norm estimator (RNE) test. For each event k, we introduce an estimator $Pk=∏i=1n∑j=1n|Ai,j|2$⁠, where A is the submatrix related to the event. Then, we can update the counter from C = 0 with the following rules:

For a boson sampler, the counter will be above zero after a few events, which is distinguished from the uniform sampler.

To distinguish from a distinguishable sampler, we use the likelihood ratio test. We also introduce an estimator $Lk=pkind/qkdis$⁠, where $pkind$ and $qkdis$ represent the event probabilities of different samplers. We update the counter from C = 0 with the following rules:

In our experiment, we set $a1=0.9$ and $a2=1.5$⁠. After collecting a few events to update the counter, for a boson sampler, the value of C tends to increase, but for a distinguishable sampler, the value of C tends to decrease.

The timestamp protocol can help to reduce the experimental execution time and increase the experiment scale. The computational complexity of the boson sampling problem is related to the calculation of the submatrix permanent of U. As a #P-hard problem, even with the best known classic algorithm, the computational steps still need $O(n2n)$⁠,⁵⁷ where n represents the number of photons. In fact, there has been no conclusion about the minimum n to achieve quantum advantage, from n = 7,¹² n between 20 and 30,¹¹ to n = 50.⁵⁸ If we consider the contribution of all combinations $(mn)$ to the computational complexity, it will take $O(n2n)(mn)$ equivalent computational steps, where m represents the number of modes. In general, it is commonly believed that m needs to be much larger than n² in order to meet the permanent-of-Gaussians conjecture (at least $m=n2$⁠). However, another view is that this condition can be extended to $m=O(n)$⁠, such as $m=2n$⁠, so we choose a quite stringent condition $m=2n$ for our estimation.

We first introduce our estimation steps. The sampling rate of n-photon boson sampling can be expressed by¹⁷ 

where R_pump represents the pumping repetition rate of the single-photon source, η represents the total efficiency of each channel in the device, including the single-photon excitation rate, fiber coupling efficiency, detection efficiency, etc., n represents the number of photons and m represents the number of modes.

Step I, as P. Clifford and R. Clifford pointed out, the fastest classical algorithm requires the calculation of $(m+n−1n)$ permanents of n × n matrices and each one takes $o(n2n)$ time. For our experiment, we measure the results in the collision-free dominant regime, so the combinations reduce to $(mn)$⁠. Then, we get the computational steps $o(n2n)(mn)$ and we assume $m=2n$⁠.

Step II, considering the state-of-art parameters, the excitation rate of the best known quantum dot source is about 0.6⁵⁹ and the moderate number of photons is about 15.⁹ We suppose all other efficiencies equal to 1. We assume sampling rate remains and set the benchmark at $η=0.6$ and n = 15. Then, we get the following equation:

We set $n′$ from 7 to 19 and calculate the corresponding efficiency $η′$⁠. Then, we obtain the different experimental scales (n and m) corresponding to different system efficiency and further calculate the computational steps, as blue histograms in Fig. 11 show.

Step III, in the main text, we have explained that the timestamp protocol can be applied in the “low-flux regime,” which means each combination only needs to occur a few times even only once. Therefore, the sampling rate and execution time required for the experiment can be effectively reduced.

With Eq. (H1), the experimental execution time is related to $SR−1$ and the required time for timestamp protocol is $γSR−1$ (γ represents the reduced proportion). Therefore, the expected reduction time is $(1−γ)SR−1$⁠. The execution time $SR−1$ depends on the system efficiency η and the system size n. So, for a certain system efficiency η, we can always reduce the experimental execution time for any size system and the absolute value $(1−γ)SR−1$ of the expected reduction time increases exponentially as the scale of the system increases.We use 100 times to estimate the actual situation and establish the following equation:

For each system scale from 7-photons 14-modes to 19-photons 38-modes, we assume sampling rate and system efficiency remain and calculate $n′$ and $m′$⁠. For example, we set n = 15, m = 30, $η=0.6$⁠, and calculate $n′≈17, m′≈34$ (taking the nearest integer). In Fig. 11, the “two-photons four-modes” and “three-photons six-modes” represent the enhancement (⁠$n′−n,m′−m$⁠) induced by timestamp protocol. The results indicate that the sampling rate and experiment execution time required for a 15-photon standard experiment and a 17-photon timestamp experiment are close. The purpose is to exchange time resources for larger system scale.

In Fig. 11, we present the required total efficiency η under different photon number and mode number. We show the enhancement of photon number and mode number further enabled by our timestamp protocol. Meanwhile, we give all the computational steps with the formula $O(n2n)(mn)$⁠. We can find that the timestamp protocol can enhance the computational complexity and reduce the system efficiency required by the experiment to a certain extent.