Abstract

The development of large-scale, high-fidelity quantum processors is a fundamental scientific challenge, essential for exploring the boundaries of classical computation and advancing towards fault-tolerant systems. Gaussian boson sampling not only serves as a prominent model for demonstrating quantum computational advantage 1 , 2 , 3 but can also generate bosonic error-correcting codes for fault-tolerant quantum computing 4 , 5 , 6 . However, its scalability has been hindered by significant photon loss in increasingly large and complex encoding circuits. Here we show a programmable photonic quantum processor, Jiuzhang 4.0, which incorporates 1,024 high-efficiency squeezed states into a hybrid spatial–temporal encoded 8,176-mode circuit. By achieving 92% source efficiency and 51% overall system efficiency, the processor produces samples with detection events up to 3,050 photons, representing an order-of-magnitude increase in scale over previous demonstrations 7 , 8 , 9 , 10 . This architecture realizes a cubic scaling of connectivity (16 3 = 4, 096), enabling sampling within a Hilbert space of dimension approximately 10 2,461 . The experimental results are rigorously validated against all current classical simulation methods, especially the matrix product state algorithms recently designed to exploit photon loss 11 . The ability to control thousands of photons in programmable low-loss quantum processors pushes the experimental frontier into a regime far beyond classical tractability and opens a pathway to trillion-qumode three-dimensional cluster states and fault-tolerant photonic quantum hardware.

Main

The transformative potential of quantum computing—the exponential speedup over classical computers for certain problems 1 , 12 —has been increasingly substantiated by recent experimental milestones 7 , 8 , 9 , 10 , 13 , 14 , 15 , 16 , 17 . These efforts have, in turn, sparked substantial improvements in classical simulation algorithms 11 , 18 , 19 , 20 , 21 , 22 , 23 , 24 . This quantum-classical competition drives the development of larger-scale, higher-fidelity quantum processors—particularly by addressing experimental noise sources that classical methods can exploit.

Gaussian boson sampling (GBS) 2 , 3 is an intermediate model for quantum computing with potential applications in various practical domains 25 , 26 , 27 , 28 , 29 , 30 , 31 , 32 , 33 and can be used for generating non-Gaussian states as a building block for universal quantum computing 4 , 5 , 6 . Previous GBS experiments detected 76, 113, 219 and 255 photon events 7 , 8 , 9 , 10 , with outcomes validated against contemporary classical spoofing algorithms 20 , 21 , 22 , 23 , 24 . However, a primary challenge in the optical experiments is photon loss. For example, previous GBS experiments achieved overall efficiencies—incorporating all sources of loss from photon generation to detection—with a value of 0.54 with spatial encoding 8 and 0.32 with time-bin encoding 9 .

Recently, a matrix product state (MPS)-based algorithm 11 quantitatively evaluated the impact of photon loss on intermediate-sized GBS experiments. The algorithm decomposes the output covariance matrix of lossy GBS into a quantum component, defined by an effective squeezed photon number, and a computationally tractable classical component. Consequently, photon loss diminishes the effective squeezed photon number, thereby reducing the computational complexity of GBS.

In this work, we report a new photonic quantum processor, Jiuzhang 4.0, which implements programmable boson sampling with up to 1,024 input squeezed states and 8,176 output qumodes. A schematic of the setup is shown in Fig. 1a. It consists of three main parts: squeezed light sources, a programmable spatial–temporal hybrid encoding circuit and a single-photon detection and coincidence system.

a, The experimental setup. Four OPOs (only two are displayed) generate the SMSSs, which are filtered and sent into a spatial–temporal hybrid encoding circuit. The circuit is composed of three 16-mode interferometers connected by two fibre-delay-loop arrays. b, Evolution of photon population and the resulting transformation matrix structure. Left, trajectory of a photon incident at an arbitrary input port of U1. Right, visualization of the transformation matrix structure. The output population for each input mode is flattened into a 4,096-dimensional array, forming a single row of the matrix. The four independent OPOs within a single input time bin are marked with different colours. The complete matrix structure merges the injected pulses across all time bins.

Single-mode squeezed states (SMSSs) are created from four optical parametric oscillators (OPOs), which are pumped by 1.6 ns-duration laser pulses centred at 775.7 nm (Supplementary Figs. 1–4). To filter out non-degenerate spectral modes, we use three cascaded unbalanced Mach–Zehnder interferometers, which feature a near-unity transmission of 99.8% and >40 dB extinction ratio (Methods). With this design, the system efficiency of our squeezed light source reaches 92%, significantly higher than the 70% reported in ref. 9 .

With varied pump strength, different squeezing parameters ranging from 0.9 to 1.8 are generated to perform the experiments. The experimentally calibrated photon indistinguishability shows the theoretically expected dependence as a function of the squeezing parameter r (Methods), reaching 97% at r = 1.8. The ground-truth is modelled as a lossy GBS as in the standard practice 7 , 8 , 9 , 10 , a more physical model incorporating partial photon distinguishability is also discussed and yields the same conclusions (Supplementary Information).

The SMSSs are then sent into a spatial–temporal hybrid encoding circuit, which is designed to achieve high connectivity, modularity and scalability simultaneously. The circuit consists of three cascaded fully connected 16-mode interferometers that are interconnected by two delay loop arrays. The first shorter delay loop array ranges from t, t + τ, …, t + 15τ, and the second longer array ranges from t, t + 16τ, …, t + 15 × 16τ. We set the temporal interval τ = 50 ns to match the recovery time of the single photon detection channels.

Figure 1b first traces the evolution of the photon population inside the circuit. Our cascaded architecture harnesses both spatial and temporal interference, so that information encoded in one degree of freedom can be transformed with the other as the photons propagate. In the first step, after U1, the photons are delocalized over all 16 spatial modes. In the second step, the subsequent array of delay loops L1 causes the photons to populate the nearest 16 time bins. A second interferometer, U2, then redistributes the photons across all 16 × 16 spatial–temporal modes. In the third step, the second array of delay loops L2 opens a new temporal dimension by delaying the photons with every 16 time bins, and finally, a third interferometer U3 again mixed all the 16 × 16 × 16 spatial–temporal modes. The experiment implements four independent OPOs operating across 256 time bins for input, which compose a total of 1,024 input squeezed states. An overview of the structure of the entire experimental input–output matrix is given in Fig. 1b (right).

The connectivity of our circuit grows cubically, whereas the physical resources, that is, interferometer size, squeezed light sources, delay lines and detectors, grow only linearly. This favourable scaling law enables unprecedented ultralarge-scale GBS experiments involving thousands of highly connected input and output qumodes.

Our photonic processor offers programmability at both the input-state preparation and the interferometric stages (Supplementary Fig. 8). An electro-optic modulator acting on the seed laser, together with the acoustic-optic modulator that gates the pumping laser, allows us to prepare arbitrary sequences of SMSSs as the experimental input. Downstream along the processor, three 16-mode interferometers incorporate thermal tunability, enabling reconfigurability over the spatial modes. Coupling between temporal modes is governed by fibre delay loops whose phase can be programmed with fast piezo-electric stretchers. Furthermore, the acoustic-optic modulator provides fast phase programmability across time bins within a single experimental run. Taken together, these elements give us control abilities over all components of our photonic circuits shown in Fig. 1a. Although this architecture prioritizes scalability and connectivity over universal gate sets—which would require millions of controls—the connectivity and depth of the implemented network enable sufficient complexity to rule out efficient simulation by state-of-the-art classical methods. Once configured at the start of an experiment, the circuit is fixed by active phase-locking loops, ensuring the machine continuously produces samples from the target ground truth distribution.

The output photons are registered by 16 superconducting nanowire single-photon detectors, with an average detection efficiency of 93% and a recovery time of 43 ns. The overall system efficiency, including detection, is 51%. The phase of the GBS setup, from the generation of the SMSSs to the interferometers, is stabilized to about λ/200 at 1,550 nm (Supplementary Fig. 18).

To characterize the performance of our quantum processor, multiple groups of data are collected under different circuit depths and input scales (Methods): the S64 group with 64 input SMSSs and 4,336 output qumodes; the M256 group with 256 input SMSSs and 5,104 output qumodes; and the L1024 group with 1,024 input SMSSs and 8,176 output qumodes. The three groups differ only in input scale but share the same circuit.

Photon number distribution of the L1024 group is shown in Fig. 2a. We observe a mean photon-click number up to 2,207, which significantly exceeds all previous results by more than one order of magnitude.

a, Photon-click number distribution for the L1024 group across various single-pulse energies, with their mean μ and standard deviations σ listed in the inset. b, A comparison of the photon number distribution of the experimental results (L1024 at 3 nJ), the ground-truth theory and classical mock-ups (squashed state, thermal state, distinguishable photons). c, Bayesian tests against the squashed state, computed for specific photon-click number as a function of subsystem size (30 clicks for 5 nJ and 6 nJ, and 15 clicks for 2 nJ, 3 nJ and 4 nJ). Experimental data from the S64 group are used for the Bayesian test limited by classical computational overhead. The subsystems for validation are randomly selected. Error bars indicate the standard error obtained from 100 resampling groups. d, e, Scatter plot of the second-order correlation function of the IPS sampler (d) and the treewidth sampler (e) in comparison with the experimental results. Coloured dashed lines indicate the linear fit for the results of each sampler. The red dashed line represents the line of unit slope at 45° (ideal sampler) for reference. Both plots use 10 million samples from the L1024 (6 nJ) dataset. f, Comparison of PCC for the second-, third- and fourth-order correlation function between the experimental data and the squashed sampler. We define the normalized PCC excess as (PCC exp − PCC squash )/PCC squash , plotted on the y-axis for visual clarity of the contrast. Data from the L1024 group (10 million samples) are shown for pump settings ranging from 2 nJ to 6 nJ. Error bars are the standard error obtained by randomly and equally partitioning the overall correlators into 100 subgroups.

First, we compare the experimental photon number distribution with the theoretical prediction of the ground truth and the mock-ups, which conveys global information of the ensemble. As shown in Fig. 2b, the experimental distribution overlaps well with the ground truth while strongly deviates from all the classical mock-ups.

Then we continue with the Bayesian test, a robust statistical tool to distinguish the theoretical model of adversary hypotheses from the ground truth based on the likelihood to generate the experimental samples.

We test the classical mock-up states that can maximally approximate the quantum light sources of the experiment among classical Gaussian states: the squashed state 9 , 10 . We start from a subsystem with fewer output bosonic modes and gradually increase the subsystem size 8 . The results in Fig. 2c show that not only are all the Bayesian scores higher than zero, which indicates that the experimental samples are more likely generated from the ground-truth distribution, but there is also a clear trend of increasing as the size of subsystem ramps up, because more complete information of the entire system is included in the test. From this observation, we can infer that greater Bayesian confidence is expected for the full system, although it cannot be directly evaluated.

Another important tool for validating the experimental GBS is the truncated correlation-function benchmarking 8 , 10 . Apart from the correlation coefficient metric, we also use the slope K of the linear fit between the theoretical correlation and the counted sampler correlation to quantify the performance of various samplers. Here we define ΔK = |K sampler − 1| as the metric of deviation from the ground truth. A smaller ΔK indicates better sampler performance. For an ideal GBS sampler, ΔK = 0.

Figure 2d rules out the independent pairs and singles (IPS) sampler designed in ref. 22 to emulate experiments of limited quantum interference. The treewidth sampler proposed in ref. 23 uses the local connectivity of GBS circuits to perform approximate simulation, which has computational complexity exponentially that is dependent on the chosen treewidth parameter. For the largest experiment, we set the treewidth to 161 (a regime beyond practical implementation) to check its approximation performance in terms of the correlation function. The significant deviation from the ground truth in Fig. 2e thus excludes this mock-up. Similarly, we address the greedy sampler 21 , which attempts to spoof GBS samples by reproducing the low-order marginals. We rule out this mock-up by analysing higher-order correlations (Methods).

To provide a comprehensive overview of the data, we use the Pearson test for up to four-order correlation, covering the full range of pumping settings of the L1024 dataset from 2 nJ to 6 nJ. A mode binning method is applied, and the details are described in the Methods. As shown in Fig. 2f, the experimental data consistently yield better correlation values compared with the squashed sampler across all orders and pump settings. We note that the contrast between the experiments and the squashed sampler becomes increasingly pronounced as the correlation order rises. This is because higher-order correlations inherently involve more modes, thereby capturing more complete global information about the system—a trend consistent with observations from the Bayesian test.

Finally, we demonstrate that our experiment vastly exceeds the reach of the MPS algorithm proposed in ref. 11 , which is currently the most powerful classical algorithm to simulate lossy GBS. The MPS algorithm works by first decomposing the output covariance matrix V of a lossy GBS device into a quantum part V p and a classical part W, V = V p + W, where V p represents an ideal but smaller-scale GBS state and W represents random displacement and can be efficiently simulated. The core of the algorithm is to construct a MPS in the Fock basis to approximately simulate the evolution of the quantum part V p . Because photon loss can reduce the effective squeezing, V p typically contains fewer photons than the physical device, allowing the MPS algorithm to simulate lossy GBS more efficiently than exact approaches.

Two crucial parameters determine the accuracy and computational complexity of the algorithm: (1) the bond dimension χ, which controls the truncation error when approximating V p with the MPS. Both runtime and memory cost scale as O(χ 2 ); and (2) the effective squeezed photon number N eff , which accounts for the hafnian calculation cost during tensor construction and evolution. The runtime grows exponentially on N eff .

To exclude an MPS simulation of our experiment, we show that no χ and N eff compatible with present-day resources can reproduce the statistics of our device within experimental uncertainty. Our benchmarking follows two complementary approaches:

1. 1.

   Targeting ideal ground truth distribution. By fixing the target of the MPS algorithm to be the ideal ground truth, we evaluate the truncation error and output statistics error of the MPS simulator with increasing experimental scale. We show that for our largest-scale experiment, the computational resources significantly surpass the abilities of supercomputers, well before the MPS error can be brought down to the level achieved experimentally.

2. 2.

   Targeting modified distribution approximating the ground truth with a tractable N eff . As an alternative investigation, in the second approach, we target the MPS algorithm to approximate the ground truth with a smaller effective photon number, keeping the truncation error sufficiently low and the runtime manageable. We find the output statistics of the sampler deviate sharply from the ideal distribution and fail standard validation tests.

We begin detailing our first approach by investigating the growth of the truncation error at various experimental scales, and consequently how the statistics of the MPS simulator deteriorate. A useful reference is the squashed-state model, obtained by replacing the quantum part V p with vacuum, which represents the classical Gaussian state simulator closest to our ground truth, that is, a baseline from which MPS must improve.

For the S64 dataset, we draw 5.5 million MPS samples across several bond dimensions and benchmark them using the two-order correlation metric ΔK (refs. 10 , 11 ). In Fig. 3a (left), we show the normalized overlineDelta K compared with the truncation error ε, and in Fig. 3a (right), we plot the detailed results for the case of χ = 10 4 . As χ increases from 10 3 to 10 4 , the ε and the simulated overlineDelta K is improved from 0.16 to 0.07 and from 0.44 to 0.34, respectively, converging towards the ground truth (overlineDelta K=0). For this dataset, the MPS outperforms the squashed state baseline (overlineDelta K=0.56). However, it still deviates significantly from the experimental results (overlineDelta K=0.14) even for the largest χ.

a, b, MPS truncation error and ΔK of the two-order correlation function counted from MPS samples, for the S64 experiment (a) and the M256 experiment (b), under progressively increasing bond dimension χ. To facilitate comparisons across different experimental scales, the left panels in a and b display the normalized deviation overlineDelta K=Delta K/Delta K 0. Here, the baseline ΔK 0 is defined as the ΔK obtained with MPS at χ = 0 (where the quantum part V p is ignored, retaining only the thermal displacement in W). The deviation of the squashed state is normalized by this same baseline ΔK 0 to serve as a direct benchmark for the MPS approximation. In a and b, 5.5 million samples from the S64/M256 (2 nJ) dataset are used. c, Truncation error of MPS in approximating the ground truth of the S64, M256 and L1024 experiments, in terms of the chosen bond dimension up to χ = 2 × 10 5 . The effective squeezed photon number N eff of the three groups is also labelled. d, Comparison of PCC for the second-, third- and fourth-order correlation function between the experimental data and the TP-MPS sampler (χ = 10,000). Error bars are the standard error. e, Total variation distance (TVD) of photon number distribution between the ground truth and the experiment/TP-MPS sampler (χ = 10,000). Ten million samples from the L1024 dataset (2–6 nJ) are used in d and e.

The M256 experiment exhibits a different picture, as shown in Fig. 3b. Even at the largest bond dimension of 10 4 , the truncation error reaches 0.86. More importantly, the two-order correlation metric overlineDelta K grows to 0.66, a point that is worse than the squashed-state baseline of 0.51. This is because, at this error level, the MPS method loses so much information that the simulated distribution statistics drift already far from the ground truth.

The experimental size-dependence of the MPS simulator is already evident from the observations above. We further present in Fig. 3c the evaluated truncation error for all three datasets up to χ = 2 × 10 5 . For the L1024 dataset, the truncation error ε reaches 0.9999 under χ = 10 4 , a level at which the MPS-simulated statistics would be essentially meaningless. Yet, even reaching this poor accuracy is already beyond classical computing. For the L1024 dataset, N eff reaches up to 114. Calculating a single MPS element requires exponential time cost hafnian evaluation. For example, running the hardest L1024 sample with χ = 10 4 would require more than 2 × 10 6 years on a 10 5 -GPU cluster such as xAI-Colossus ( ). Consequently, no realistic classical computing resources, to our knowledge, can bring the MPS algorithm anywhere near the accuracy achieved by our experiment.

We next investigate the complementary aspect of whether a simulatable MPS sampler, with smaller N eff and therefore manageable run time as well as negligible truncation error, could approximate the ideal ground truth closer than our experiment does. Concretely, we target the MPS algorithm to simulate a distribution that is interpolated between the ideal ground truth and the squashed state, with the latter being the classical state closest to the ground truth. To construct such a sampler, we artificially reduce the circuit transmission to lower N eff while increasing the input squeezing so that the mean photon number is left unchanged. We refer to such a tailored MPS simulator as a tractable proxy MPS (TP-MPS).

With the bond dimension fixed at 10 4 , we scan and choose, for each experimental size, a transmission that yields a favourable truncation error of 0.01. Running this specially tailored MPS simulator on the S64, M256 and L1024 datasets does improve the squashed state baseline, yet its outputs still differ markedly from the ground truth in comparison with our experimental distribution (Methods). Figure 3d presents a comprehensive comparison of up to fourth-order correlations across all pump settings (2–6 nJ) of the L1024 dataset. A consistently positive Pearson correlation coefficient (PCC) excess demonstrates that the experimental data outperform the TP-MPS sampler across all settings and correlation orders. Furthermore, the experiment more accurately reproduces the theoretical photon-number distribution (Fig. 3e). More detailed validation data for all datasets and power settings are provided in Supplementary Information section 8. Furthermore, we have also performed cross-entropy benchmarking (XEB), yielding consistent conclusions (Methods).

These results exclude all possibilities that the MPS algorithm could reproduce the ground truth distribution better than our experiment. We note that the observations here are consistent with the asymptotic criterion in ref. 11 , which states that the MPS algorithm would keep being efficient as long as the transmission η decays as eta =o(1/sqrtK), where K is the number of input squeezed states. In our experimental architecture, the scale-up is achieved at a fixed transmission rate, in which case the MPS is driven outside its efficient regime as the experimental size grows.

Finally, we benchmark the classical computational overhead to simulate our experiments by estimating the computational cost required for the state-of-the-art classical approach—the MPS algorithm 11 —to generate equivalent samples. For simulating an M-mode GBS experiment with a local Fock-space cut-off dimension d, the asymptotic runtime of the MPS algorithm is

where χ(ε) is the bond dimension at a truncation error ε, and N eff is the effective squeezed photon number. We note that the complexity accounts for only the cost of tensor network construction.

The computational cost is dominated by the two exponential factors, 2^N_eff/2 and χ(ε) 2 . We compute N eff exactly from the squeezing parameters and the interferometer matrix (Supplementary Table 2), finding N eff ≈ 114 for our largest experiment. Although this value is lower than the detected photon clicks due to losses, the computational cost remains prohibitive because of the exponential dependence. The reported N eff ≈ 114 also significantly outperforms that of all previous GBS experiments with N eff up to about 10. Computing χ(ε) exactly is infeasible in this regime, because the required tensor calculations exceed any realistic memory budget. Instead, we extrapolate from the tractable regime to estimate the required χ. Following ref. 11 , which shows that χ(ε) = O(polylog(1/ε)), we find our data to be well described by χ(ε) = A(ln(1/ε)) n (Supplementary Fig. 23), with coefficients A and n determined by least-squares fitting. Choosing a conservative truncation error of ε T = 0.05 yields a lower bound of the required χ > 8 × 10 21 for our largest L1024 dataset (6 nJ), as shown in Fig. 4a.

a, Fitted relation between the bond dimension χ and the truncation error ε, and the extrapolation to a target truncation error ε T to estimate the required size of χ. b, Estimated speedup ratio against MPS simulation, which is based on El Capitan, the most powerful supercomputer at present. The results are presented in terms of the two crucial factors, the bond dimension χ and the effective squeezed photon number N eff .

Having determined the least required χ, we then estimate the time cost for the MPS algorithm using equation (1). We choose El Capitan–the most powerful supercomputer at present 34 , as the classical baseline. To obtain the empirical prefactor of equation (1), we first benchmark the runtime of the MPS simulation in a regime in which the direct simulation is feasible 35 . We then scale the measured runtime to the intractable regime according to equation (1), implicitly granting the classical supercomputer unlimited memory—an assumption that favours the classical side. Figure 4b summarizes the projected runtime advantage ratios for all three experimental instances reported in this work, as well as earlier GBS demonstrations. For our largest experiment, El Capitan would require more than 10 42 years to produce a sample, whereas our quantum processor takes only 25.6 μs. This corresponds to a large quantum speedup ratio exceeding 10 54 .

In summary, we have constructed a programmable photonic processor capable of executing one of the largest boson sampling experiments so far. Featuring up to 1,024 input squeezed states and 8,176 output qumodes, the system yields a Hilbert space dimension of about 10 2,461 . This configuration allows us to access regimes that are computationally intractable for the state-of-the-art classical algorithms, including the MPS simulations, even in the presence of photon loss.

Beyond foundational interest, this architecture also enables immediate applications in the NISQ era, such as image recognition 33 and cryptographic one-way function 36 . Moreover, the demonstrated technological advances, including the high-efficiency squeezed source, the hybrid spatial–temporal encoding with full connectivity and distributed active phase stabilization for modularity, open up immediate applications in generating large-scale three-dimensional cluster states and implementing bosonic error-correcting codes 37 , 38 , 39 , 40 . We anticipate that the inherently low-loss nature of this platform will enable further investigations into non-Gaussian resource generation and error-corrected quantum information processing 4 , 5 , 6 , 41 , 42 , 43 .

Methods

Squeezed light sources and spectral filtering

The optical parametric oscillators (OPOs) are pumped by rectangular optical pulses (1.6 ns duration, 775.7 nm centre wavelength) to produce down-converted photons at 1,551.4 nm. The pump pulse duration is selected to satisfy the single temporal mode condition imposed by the OPO design, whereas a 50 ns repetition interval aligns with the delay loops. The OPOs are continuously locked by the Pound–Drever–Hall (PDH) technique, with a phase-modulated rectangular probe light pulse with 1 μs duration. These probe pulses share the polarization and centre frequency of the pump pulses but are temporally interleaved to prevent cross-talk (Supplementary Fig. 4).

To achieve narrowband spectral filtering with minimal photon loss, we implement cascaded unbalanced Mach–Zehnder interferometers (MZIs). By introducing an optical delay difference L between the MZI arms, the transmission at one output port exhibits a sinusoidal frequency dependence,

By setting L equal to the cavity length L c and aligning the transmission maximum ν′ to the zeroth-order spectral mode ν 0 of the cavity, frequencies at nu 0 +/- fracc2L c ,nu 0 +/- frac3c2L c ,nu 0 +/- frac5c2L c ,... acquire a π phase. This results in destructive interference that suppresses the odd-order spectral modes (±1, ±3, ±5, …) of the cavity. More generally, an unbalanced MZI with L=frac2n+12 n L c filters out the spectral modes of order ±2 n , ±3 × 2 n , ±5 × 2 n ….

We implement three cascaded unbalanced MZIs for each OPO with optical path differences of L=L c ,(3)/(2)L c ,(5)/(4)L c (Extended Data Fig. 1). This configuration suppresses the most non-degenerate spectral modes, except for orders ±8, ±24, ±40, …. Each MZI consists of two balanced beam splitters and a mirror, forming a specially designed triangular interference structure. The incidence angle of the mirror is set at 20.6°, ensuring the beams strike the beam splitters at the Brewster angle (55.3°). Consequently, the p-polarized beam, as used in the experiment, experiences negligible reflection loss without anti-reflection coating. Residual spectral components are removed by three subsequent layers of 0.4 nm narrow-band-pass filters.

With careful mode matching tuned to >99.95% and an average transmission of the narrow band pass filters >98.5%, we achieve a >40 dB rejection ratio of all high-order spectral modes and an total transmission efficiency above 95%.

Cross-OPO interference calibration

As our setup uses multiple OPOs for generation of squeezed light, ensuring the indistinguishability of photons from distinct OPOs is important. We quantitatively calibrate the interference indistinguishability of the sources by unheralded second-order correlation measurements. These measurements are performed using a Hanbury–Brown–Twiss (HBT) experiment on thermal states generated by interfering two SMSSs (see Supplementary Information section 1.3 for details).

We conduct these measurements in two scenarios: interference from a single OPO using a delayed self-interference scheme and the direct interference between two independent OPOs. The results for both cases are nearly identical and agree well with theoretical simulations (Extended Data Fig. 2).

Details of the photonic circuit

The circuit architecture, shown in Extended Data Fig. 3, comprises three cascaded ultralow-loss, fully connected 16-mode interferometers interconnected by two delay loop arrays. Each interferometer uses a multi-layer Mach–Zehnder-type interference structure, which features both compactness and compatibility with external thermal tunability. Specifically, each unit consists of two perpendicularly oriented four-layer triangular interferometers (Extended Data Fig. 3, bottom left), providing full connectivity across the 4 × 4 = 16 spatial modes. The transmission efficiency of each 16-mode interferometer is 99% per mode, with an average wave-packet overlap exceeding 99.7%. The delay loop arrays form arithmetic progressions with common differences of τ for the shorter array and 16τ for the longer array, where τ = 50 ns is chosen to match the recovery time of single-photon detectors. Photons traversing the first interferometer and short delay loop array are split into 16 sub-beams, arriving at time bins [t, t + τ, …, t + 15τ]. These sub-beams then enter the second interferometer and the long delay loop array, where the ith loop directs beams to time bins [t + i × 16τ, t + (1 + i × 16)τ, …, t + (15 + i × 16)τ]. This configuration enables each input mode to couple densely to the subsequent 16 2 temporal modes and redistribute across all 16 spatial modes using the third interferometer.

Phase stabilization of the system

Phase stabilization is important for GBS experiments. To ensure the stability of the entire interference network, our system incorporates three active servo loops for phase locking.

The first phase-lock module stabilizes the phase from the initial splitting of pump pulses to the first interferometer, as shown in Extended Data Fig. 4a. A mirror mounted on a piezoelectric transducer (PZT) is positioned in the path of the pump beam to each optical parametric oscillator (OPO). Because the pump beam co-propagates with the squeezed light through the cavities, it serves as an ideal probe for feedback. A dichroic mirror, transparent at 1,550 nm and reflective at 775 nm, separates the pump light, which is then interfered with a reference beam. The resulting error signal, obtained by balanced homodyne detection, is used for phase locking by adjusting the PZT displacement.

The second and third phase-lock modules stabilize the phases of the two fibre delay loop arrays, as shown in Extended Data Fig. 4b. Taking the first array as an example (the same method applies to the second), a PZT-mounted mirror is positioned in the path of each spatial mode for active phase compensation. A counter-propagating probe beam, frequency-matched but orthogonally polarized to the squeezed light, is split into 16 beams and coupled to each spatial mode at the input of the second interferometer by a polarization beam splitter. To minimize stray counts, this probe beam is pulsed using an acousto-optic modulator and temporally offset from the squeezed light pulse train. At the output of the first interferometer, the probe is separated by another polarization beam splitter and interfered with a reference beam to generate an error signal that drives the PZT. To address significant phase drift in long fibres due to temperature fluctuations, which often exceeds the compensation range of a single PZT, a thermoelectric cooler feedback system is implemented to mitigate slow thermal drift.

The phase-locking precision is estimated to be 1 nm for the OPO source array, 3 nm for the first delay loop array and 10 nm for the second delay loop array, and is continuously maintained over tens of hours. The integration of high-efficiency squeezed light sources working at 20 MHz repetition rate and spatial–temporal hybrid encoding architecture—underpinned by robust, long-term phase stabilization—establishes a viable pathway towards the control of trillion-qumode 3D cluster states.

Multiple groups of experiments with different circuit depth, input scale and unitary transformation

We collected multiple groups of data under different circuit depth and input scales to fully characterize and validate the performance of our quantum processor, as shown in Extended Data Fig. 5. We first finished a group of experiments under the miniature setup with a shallower circuit consisting of the first two interferometers, which contains 3 × 16 = 48 input SMSSs and 16 × (16 + 15) = 496 output qumodes. For the whole setup comprised of all the three interferometers, experiments under three groups of progressively expanding input scale were conducted: the S64 group with 4 × 16 = 64 input SMSSs and 16 × (16 + 255) = 4,336 output qumodes; the M256 group with 4 × 64 = 256 input SMSSs and 16 × (64 + 255) = 5,104 output qumodes; the L1024 group with 4 × 256 = 1,024 input SMSSs and 16 × (256 + 255) = 8,176 output qumodes.

To experimentally demonstrate the programmability of our system, two datasets were collected with the respective distinct unitary transformations. The first dataset (dataset A), is thoroughly examined through the paper and used for QCA benchmarking; the second dataset (dataset B), which produces up to 3,050 photon-detection events, is presented in Supplementary Information section 9, together with supportive validation results.

Bayesian test

The Bayesian test score ΔH between the ground-truth model H 0 and the classical adversary model H 1 on a set of n-photon-click samples is defined as

where N is the number of samples, vecs i represents the ith sample, P (0/1) (vecs i ) represents the event probability of samples vecs i under the H 0/1 hypothesis, and P (0/1) (n) is the coarse-grained probability for the n-photon click under hypothesis H 0/1 . ΔH > 0 indicates that the experimental samples are more likely from the ground truth GBS rather than the mock-up. A higher score indicates higher confidence. As shown in Extended Data Fig. 6, we observe that ΔH is independent of the detected photon-click number k across the calculated range for various datasets (S64, M256 and L1024). Consequently, when investigating the scaling behaviour of ΔH with respect to subsystem size, we may specify a fixed detected click number k without loss of generality.

In Fig. 2c and Extended Data Fig. 6, we observe a decrease in the Bayesian score as the pumping power increases. This trend arises from two compounding factors. First, owing to computational limits, we can validate only a fixed-size subsystem. As the global experiment size grows (up to about 3,000 photons), the validated subsystem represents a greatly shrinking fraction of the whole (down to about 1%), effectively introducing severe truncation loss that reduces the distinguishability of the signal. Second, regardless of photon loss, the squeezed state and its classical squashed-state counterpart become intrinsically harder to distinguish at higher pumping powers. Numerical simulations confirm these mechanisms, demonstrating strong agreement between our experimental observations and theoretical expectations.

To address the impact of subsystem size, we simulated ideal samples for a moderate GBS system (10 SMSS input, 100 modes output) solvable with available computational power. We selected 6-click samples and performed a Bayesian test while progressively decreasing the subsystem size. As shown in Extended Data Fig. 7a, ΔH drops significantly as the subsystem size decreases.

Finally, we examine the relation between the Bayesian test score and the pump strength (or the input squeezing). We generate ideal samples using the exact GBS simulation method for a small GBS system (5 SMSS input, 10 mode output) with varying input squeezing. The simulation is conducted under two scenarios: one with the same photon loss rate as the experiment and the other without any photon loss. As shown in Extended Data Fig. 7b,c, ΔH decreases with increased squeezing, reproducing the experimental trend in both the lossy and lossless cases. This confirms that the observed decrease stems from the inherent convergence of the target and squashed models at higher squeezing levels, rather than from experimental noise or degradation in source quality.

Correlation function benchmarking

The k-order correlation function is recursively defined as

where X k represents the experimental measurement operator at the kth output mode and P k represents all the partitions of 1, 2, …, k except the universal set.

We implement a mode-binning technique that aggregates multiple physical output modes into a single binned port—a detection event in any of the constituent physical modes registers as a click on the binned port (details provided in Supplementary Information section 6). This method is particularly advantageous for analysing correlations beyond the second order, as it mitigates the prohibitive costs on samples to suppress statistical fluctuations. In the processing of Figs. 2f and 3d, each consecutive 8 modes are binned together. For the second- and third-order correlations, all correlators are included (except for a small part of the data points whose theoretical values are strictly zero). For the fourth-order correlation, the total 8,176 modes are partitioned into four disjoint parts [1, 2,048], (2,048, 4,096], (4,096, 6,044], (6,044, 8,176] to count the interior correlation statistics separately and merge the results as a whole. The standard error in Figs. 2f and 3d is obtained by randomly and equally partitioning the overall correlators into 100 subgroups.

Extended Data Fig. 8 shows a comparison of the third-order correlation between our experimental samples and the mock-up samples from the so-called greedy sampler proposed in ref. 21 , which aims to spoof GBS samples by reproducing the low-order marginals. The mock-up is shown not able to capture higher-order correlations, and thus is conclusively ruled out.

Tractable proxy MPS

As the standard MPS is computationally intractable for large-scale experiments, we investigate whether a proxy MPS sampler—one with a reduced effective squeezed photon number (N eff ) and negligible truncation error—can better approximate the ideal ground truth than our experimental hardware. Concretely, we define a trajectory of states interpolating between the intractable ground truth and the squashed state. This is achieved by scaling down the circuit transmission to reduce complexity, while simultaneously increasing the input squeezing to maintain a fixed mean photon number.

Extended Data Fig. 9a quantifies the computational cost in terms of truncation error along this trajectory. With the bond dimension fixed at χ = 10 4 , we observe that the truncation error (y-axis) remains prohibitively high near the experimental operating point (left dashed line). To define a valid comparison, we select a transmission rate for each system size that yields a manageable truncation error of ϵ ≈ 0.01. As the system size increases, this requirement forces the algorithm to revert significantly towards the squashed state. Consequently, although this tractable proxy MPS (TP-MPS) simulator improves on the squashed state baseline (Extended Data Fig. 9b–d), its output distribution still deviates significantly from the ground truth compared with the experimental data.

Cross-entropy benchmarking

Apart from the correlation function presented in the main text, the cross-entropy benchmark (XEB) is another tool to compare which set of samples is closer to the ground-truth model. The relative XEB between two samplers are defined as

where N is the number of samples, vecs i represents the ith sample, P(vecS i ^exp(c)) represents the event probability of samples from the experiment (classical mock-up sampler), calculated with the GBS ground-truth model. A positive ΔH XEB indicates that the experiment outperforms the corresponding classical mock-up sampler in achieving higher cross-entropy with respect to the ground truth. The relative XEB results between experimental samples and the classical samplers (squashed state and MPS) considered in Fig. 3 are shown in Extended Data Fig. 10. Our experiment surpasses both the squashed mock-ups and the MPS algorithm under this test. For the large experiments containing thousands of photons, the XEB benchmark is infeasible, even taking the subsystem for validation, because the resultant subsystem is too small to output informative results due to the excess artificial dropping out of the global information.