Quantum advantage, benchmarking the computational power of quantum machines outperforming all classical computers in a specific task, represents a crucial milestone in developing quantum computers and has been driving different physical implementations since the concept was proposed. A boson sampling machine, an analog quantum computer that only requires multiphoton interference and single-photon detection, is considered to be a promising candidate to reach this goal. However, the probabilistic nature of photon sources and the inevitable loss in evolution network make the execution time exponentially increasing with the problem size. Here, we propose and experimentally demonstrate a timestamp boson sampling scheme that can effectively reduce the execution time for any problem size. By developing a time-of-flight storage technique with a precision up to picosecond level, we are able to detect and record the complete time information of 30 individual modes out of a large-scale 3D photonic chip. We perform the three-photon injection and one external trigger experiment to demonstrate that the timestamp protocol works properly and effectively reduce the execution time. We further verify that timestamp boson sampler is distinguished from other samplers in the case of limited datasets through the three heralded single photons injection experiment. The timestamp protocol can speed up the sampling process, which can be widely applied in multiphoton experiments at low-sampling rate. The approach associated with newly exploited resource from time information can boost all the count-rate-limited experiments, suggesting an emerging field of timestamp quantum optics.

## I. INTRODUCTION

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.^{3} 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.^{4} 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,^{11} 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.

## II. EXPERIMENT

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 $\lambda =\tau \u22121$, the reconstruction probability

*p*of each combination can be calculated by Eq. (1) and the output distribution satisfies the normalization $\u2211ipi=1$ (see Appendix A for details).

_{i}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

*τ*of the trigger as the timestamp to mark the occurrence of each combination.

_{i}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.^{34} 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^{37} 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.^{39} 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 $\pi /2$ or $\u2212\pi /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.

## III. RESULTS

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*) in a five-layer ring,

_{i}^{43}as shown in Fig. 3(c). The similarity between the two reconstruction protocols $S=\u2211ipici$ is up to 98.7%. The total variation distance $D=(1/2)\u2211i|pi\u2212ci|$ 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^{45} and the likelihood ratio test^{46} (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*) 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=\u2211ipiti$ and $D=(1/2)\u2211i|pi\u2212ti|$. For a perfect boson sampler,

_{i}*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

^{16}and change the final photon distributions, so deviated from ideal distributions. By using site-by-site addressing technology

^{37}and imaging-based large-scale correlation measurement,

^{49}we are already trying to avoid these additional losses and improve this problem.

## IV. SUMMARY

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,^{50} noisy intermediate-scale quantum (NISQ) technology,^{4} and multiphoton entanglement,^{8} driving our capacity of controlling quantum systems to an entirely new accessible realm.

## ACKNOWLEDGMENTS

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.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**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).

## DATA AVAILABILITY

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

### APPENDIX A: TIMESTAMP PROTOCOL

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.^{54}

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*_{1}. For the convenience of description, we define *t*_{1} as 1 s. Based on the fundamental assumptions above, the number of events detected per second satisfies the Poisson distribution $P(x=k)=\lambda kk!e\u2212\lambda ,k=0,1\u2026$. The mathematical expectation $E=\u2211k=0\u221ekP(k)=\lambda $, where *λ* represents the average number of events detected per second. We can also calculate the variance $D=\u2211k=0\u221ek2P(k)=\lambda $, so the error bar of each measurement is $\lambda $.

On the other hand, a Poisson process is a renewal process in which the interarrival intervals have an exponential distribution function.^{55} Therefore, we can revisit this measurement from the timestamp perspective as shown in Fig. 5(b). We detect the execution time *t*_{1}, *t*_{2}, $t3\u2026$ and obtain the time intervals *t*_{1}, $t2\u2212t1,t3\u2212t2\cdots $. Statistically, these time intervals satisfy the exponential distribution. The probability density $f(x)=\lambda e\u2212\lambda x,x\u22650$ and the distribution $F(k)=1\u2212e\u2212\lambda k,k\u22650$. Therefore, the mathematical expectation $E=\u2211x=0\u221exf(x)=1/\lambda $ and the variance $D=\u2211x=0\u221ex2f(x)=1/\lambda 2$, where $1/\lambda $ 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=\lambda i/\u2211i\lambda i$. Corresponding to the timestamp protocol, we change the

*λ*as $\tau i\u22121$, which is Eq. (1) in the paper.

_{i}### APPENDIX B: TIME-OF-FLIGHT STORAGE TECHNIQUE

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 $\tau 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 $\tau i(j)$. When a combination occurs twice, the registration time are $\tau 1(j)$ and $\tau 2(j)$. The time interval is $\tau 1(j)$ and $\tau 2(j)$ − $\tau 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 $\tau 2(j)/2=(\tau 1(j)+(\tau 2(j)\u2212\tau 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.

### APPENDIX C: FABRICATION OF 3D PHOTONIC CHIP AND THE SCALABILITY OF THE EXPERIMENT

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.^{56}

### APPENDIX D: CHARACTERIZATION OF 3D PHOTONIC CHIP

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.

### APPENDIX E: PREPARATION OF QUANTUM LIGHT SOURCE

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.

### APPENDIX F: TIMESTAMP RECONSTRUCTION PROTOCOL

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\tau $, which is defined as $N\tau =T/\tau $. *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*fluctuates but most often stays in the range of $[Na\u2212Na,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\tau <2Na$. We normalize the data in the subspace and compare it with the distribution obtained from theory and standard boson sampling.

_{a}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.

### APPENDIX G: MORE EXPERIMENTAL RESULTS

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.

### APPENDIX H: VALIDATION OF BOSON SAMPLING DATA

To distinguish from a uniform sampler, we use the row-norm estimator (RNE) test. For each event *k*, we introduce an estimator $Pk=\u220fi=1n\u2211j=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.

### APPENDIX I: THE ADDITIONAL ENHANCEMENT OF THE TIMESTAMP PROTOCOL

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)$,^{57} 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,^{12} *n* between 20 and 30,^{11} to *n* = 50.^{58} 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*^{2} 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^{17}

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\u22121n)$ 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^{59} and the moderate number of photons is about 15.^{9} We suppose all other efficiencies equal to 1. We assume sampling rate remains and set the benchmark at $\eta =0.6$ and *n* = 15. Then, we get the following equation:

We set $n\u2032$ from 7 to 19 and calculate the corresponding efficiency $\eta \u2032$. 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\u22121$ and the required time for timestamp protocol is $\gamma SR\u22121$ (*γ* represents the reduced proportion). Therefore, the expected reduction time is $(1\u2212\gamma )SR\u22121$. The execution time $SR\u22121$ 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\u2212\gamma )SR\u22121$ 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\u2032$ and $m\u2032$. For example, we set *n* = 15, *m* = 30, $\eta =0.6$, and calculate $n\u2032\u224817,\u2009m\u2032\u224834$ (taking the nearest integer). In Fig. 11, the “two-photons four-modes” and “three-photons six-modes” represent the enhancement ($n\u2032\u2212n,m\u2032\u2212m$) 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.

## References

^{14}-dimensional Hilbert space