Phase transitions in random circuit sampling

The computational complexity of quantum systems arises from the exponential growth of the Hilbert space dimension with system size. On near-term quantum processors whose practical complexity is limited by noise, random circuit sampling (RCS) has emerged as the most suitable candidate for a beyond-classical demonstration. The interplay between computational complexity and noise is highlighted by recent RCS experiments with increasing system sizes and fidelity^4,5,6,9,10, although classical algorithms have also advanced substantially^{11,12,13,14,15,16}. RCS is, arguably, the entry point into the realm of classically intractable problems for any experimental quantum processing platform¹⁷. The reason is that RCS circuits can be optimized to maximize the speed of quantum correlations with iSWAP-like gates^2,18,19,20 while preventing potential simplifications in the corresponding classical emulations¹⁷. This intensifying quantum–classical competition motivates two questions. Are there well-defined boundaries for the region where the exponentially large Hilbert space is, in fact, leveraged by a noisy quantum processor? More importantly, can we establish an experimental observable that directly probes these boundaries?

In this work, we provide direct insight to these two questions using RCS on a two-dimensional grid of superconducting qubits. We demonstrate that the interplay between quantum dynamics and noise can drive the system into distinct phases. The boundary between these phases can be resolved using finite-size studies with a fidelity estimation technique called cross-entropy benchmarking (XEB)^2,3,4,21. Reaching the desired phase of maximized complexity requires a noise rate per cycle below a critical threshold whose value is determined by the growth rate of quantum correlations.

The structure of these phases is schematically illustrated in Fig. 1. Driven by the circuit number of cycles or depth, the system first goes through a dynamical phase transition where the output distribution is no longer concentrated in a fraction of bit strings. Anti-concentration is a key ingredient of XEB and of mathematical arguments on the complexity of simulating noiseless RCS^{2,22,23,24,25}. Nevertheless, we will show that this is a necessary but not sufficient condition for global entanglement (Supplementary Information section H), which maximizes the computational cost.

The second transition is driven by noise, specifically the error rate per cycle ϵ × n, where ϵ is the error per gate and n is the number of qubits. As illustrated in Fig. 1, the behaviour of quantum correlations falls into two regimes: when the error rate per cycle is large, the state of the system can be approximately represented by several uncorrelated subsystems (equation (3)). This leaves the quantum system open to spoofing by classical algorithms that represent only part of the system at a time^22,26,27. In the regime where the error rate per cycle is sufficiently low, correlations eventually span the entire system, thus restoring its computational complexity, and the experiment cannot be spoofed (Supplementary Information section F). The boundary between these two phases is determined by the competition between ϵ × n and the convergence of the system to the ergodic state. Reference ²⁸ studied the same phase transition for one-dimensional and all-to-all connectivity.

We found that XEB is a proper observable and can be used to resolve the aforementioned regimes experimentally, as it is sensitive to the nature of the dominant correlations. Specifically, linear XEB is measured as

where n is the number of qubits, \({p}_{{\rm{sim}}}(s)\) is the ideal (simulated) probability of bit string s and the average is over experimentally observed bit strings s. We measured XEB as a function of the number of cycles or depth d for different system sizes to resolve the dynamical phase transition (Fig. 1). The experimental results are shown in Fig. 2a,b for one- and two-dimensional systems, respectively. We found that XEB increased with system size for small d. In this case, the system state was concentrated in a fraction of basis states. However, for large d, XEB decayed exponentially and approximated the circuit fidelity. At an intermediate number of cycles, we observed a single crossing point where all the measured XEB curves intersected and where XEB was approximately independent of the system size. The detailed theory presented in Supplementary Information section E shows that this is indeed a phase transition associated with the role of the boundary.

Having identified the minimum number of cycles at which XEB approximated the system fidelity, we formulated an experimental protocol for locating the transition between the strong- and weak-noise regimes (Fig. 1). A conceptually simple set-up that highlights the underlying physics for this transition is the so-called weak-link model, where two subsystems of size n/2 are coupled through an entangling gate applied every T cycles. In the limit where T = ∞ (no weak link is applied), the subsystems were uncoupled and the overall system converged to a product state ρ_A ⊗ ρ_B, where ρ_A and ρ_B are the pure ergodic states of each subsystem. Adding noise, we assumed the so-called depolarizing channel noise model for the density matrix of each subsystem A and B (ref. ²⁹): F^d/2ρ_A/B + (1 − F^d/2)I_A/B/2^n/2, where I_A (I_B) is the identity matrix and F = e^−ϵn is the fidelity. Direct substitution of this density matrix into equation (1) gives XEB = e^−ϵnd + 2e^{−ϵnd/2}. We use XEB = 1 for ρ_A/B and XEB = 0 for I_A/B/2^n/2.

At finite yet large T, each subsystem approached the ergodic state in less than T cycles. The application of a two-qubit gate between the two subsystems built global correlations, which, therefore, decreased the XEB term proportional to F^d/2 with some rate λ. This rate λ depends on the gate and was λ = 1/4 for the iSWAP-like gates employed in our experiment (Supplementary Information section D). Therefore, a simplified model for linear XEB is

We justify this equation by formally averaging circuit instances over two-replicas, which results in the so-called population dynamics formalism (Supplementary Information section D). We probed this behaviour by measuring XEB experimentally as a function of d, as shown in Fig. 2c. We employed a noise-injection protocol that effectively changed the gate fidelities in our quantum circuits (Supplementary Information section C2) and show results corresponding to different noise levels. We used the discrete set of single-qubit gates chosen randomly from Z^pX^1/2Z^−p with p ∈ {−1, −3/4, −1/2, …, 3/4} and Z and X the corresponding Pauli matrices. We observed that in the weak-noise regime, XEB converged to the expected fidelity of the entire system, F^d. This was because F was sufficiently high such that F^d dominated the contribution to XEB. On the other hand, XEB was significantly above F^d in the strong-noise regime owing to the dominant contribution of 2λ^d/TF^d/2 to XEB. These results exemplify the competition between the exponential decay of global correlations proportional to ϵn and the entangling rate between subsystems (proportional to 1/T in this example).

The transition between the two different noise-induced phases was more clearly seen by fixing d to a few values and varying the effective noise level (F). We observed that XEB exhibits a non-trivial scaling distinct from F^d (Fig. 2d). In particular, the rate of decay with respect to errors decreased at higher error rates. This is also consistent with 2λ^d/TF^d/2 being dominant at high errors.

To experimentally locate the critical value of the error per cycle (or equivalently, F) where the dynamical exponent of XEB changed, we defined a modified order parameter F^d/XEB, which is asymptotic to a distinct value of 1 (0) in the weak (strong) noise regime. The transition between the two limits became a discontinuity when d → ∞, indicating a phase transition for finite ϵn ≈ κ_c, where the critical value κ_c can be a function of λ and T. In the transition region, we observed the finite-size critical behaviour where F^d/XEB was approximately a function of (ϵn − κ_c)d. This was revealed in the order parameter for different numbers of cycles d crossing at a single point, as can be verified from equation (2) and numerically for the circuits used in the experiment. The experimentally obtained F^d/XEB, shown in Fig. 3a–c,e,f, indeed manifested the expected critical behaviour. For ϵn ≲ κ_c, the order parameter increased as d was increased, whereas for ϵn ≳ κ_c, the order parameter decreased as d was increased. At the critical point ϵn ≈ κ_c, the datasets crossed and the order parameter was approximately independent of d. We attribute the slight drift in the crossing point between different depths to potential systematic errors in the experimental estimation of F.

We extracted the critical noise rate for different link frequencies 1/T experimentally, numerically and analytically (Fig. 3d). When comparing with the functional form \({\epsilon }n\simeq 4/T\log 2\) predicted by the analytical weak-link model, we observed appreciable deviations when the link frequency approached 1/2, which corresponds to the regular one-dimensional chain. The deviation occurred because the weak-link model was no longer applicable in this regime. A more accurate description is provided by a generalization of this model in Supplementary Information section E1.

To explain the nature of the noise-induced phase transition for two and more dimensions, we extended the weak-link model from two subsystems to a number of subsystems equal to the number of qubits (see Supplementary Information section E2 for details). For depths d ≫ 1, the subsystems converge either to an appropriate ergodic state or to a state proportional to the identity matrix, as in the weak-link model. In the noiseless case, there are two stationary configurations, with either all subsystems in the identity matrix or all in the ergodic state. We expanded the noisy state around these two configurations using the so-called dilute flipped spin expansion. This generalizes equation (2) to

where the first (second) term on the right-hand side corresponds to the expansion around the identity matrix (ergodic state). See Supplementary Information section D for the derivation. The index k describes the total number of subsystems not equal to the corresponding global state and k! is its factorial. The first term on the right-hand side corresponds to the low-weight Pauli paths approximation in ref. ²⁷. Finally, λ = 1/4 for the iSWAP gate, as in the weak-link model.

The noise-induced phase transition appears in the thermodynamical limit n → ∞ at fixed ϵn and \(d/\log n\) (see Supplementary Information section E for details). We found

The second factor on the right-hand side corresponds to the global fidelity and dominates in the weak-noise regime \({\epsilon }n\ll {\rm{ln}}\,4\). The first factor describes the convergence to anti-concentration (Supplementary Information section E), in accordance with refs. ^22,23, and dominates in the strong-noise regime \({\epsilon }n\gg {\rm{ln}}\,4\). In this case, the first term on the right-hand side of equation (3) prevails, and by taking small k, the state can be approximately represented by several uncorrelated subsystems. On the one hand, as captured by equation (3), this transition corresponds to qualitatively different configurations. On the other hand, both regimes correspond to an ordered phase. The noise-induced phase transition is driven by a control parameter (analogous to a magnetic field in an Ising model) that scales with the system size (number of qubits). This is loosely like Fréedericksz transitions in liquid crystals³⁰ and is qualitatively different from the quantum to classical transition discussed in ref. ³¹. The transition discussed here is a competition between the finite rate of convergence to the overall ergodic state and the fidelity per cycle. The transition in ref. ³¹ is a competition between local interactions and the error rate per qubit and is related to the error threshold of quantum error correction codes.

The two-dimensional experimental results are shown in Fig. 3e,f for a 4 × 4 square grid of qubits and two different circuit structures, whereby the two-qubit gates were applied either in a staggered (Fig. 3e) or an unstaggered (Fig. 3f) fashion. As in one dimension, the 16-qubit system was divided into two halves connected by a single iSWAP-like gate applied every four cycles. For both circuit structures, we observed a similar crossing between F^d/XEB measured at three different cycles, with a higher value of ϵn observed for the unstaggered patterns.

We numerically evaluated critical noise rates for systems of different sizes and circuit structures without weak links, including both the staggered and unstaggered patterns and the ABCD-CDAB pattern used in ref. ⁴. Here, we used Haar random single-qubit gates (Supplementary Information section D).

The vertical line in Fig. 3g gives a lower bound for the noise-induced phase transition. Furthermore, the noise-induced phase transition for the discrete gate set used in the experiment occurred at higher noise rates (Supplementary Information section F2). We compared this lower bound with the 67-qubit and 70-qubit RCS experiments that we will present next. We obtained the error per cycle by fitting the exponential decay of the fidelity (Fig. 4). It is evident that these experiments fell well within the weak-noise regime, satisfying the requirement to fully utilize the computational capacity of the noisy quantum processors.

We demonstrated beyond-classical RCS by performing an experiment on a 67-qubit Sycamore chip (Fig. 4). These random circuits followed the same two-dimensional ABCD-CDAB pattern as ref. ⁴. The single-qubit gates were chosen randomly from the discrete set Z^pX^1/2Z^−p with p ∈ {−1, −3/4, −1/2, …, 3/4}. We show in Supplementary Information section B the fidelity of the elementary operations of the random circuit. On average, we achieved a read-out error of 1.3(0.5) × 10^−2, a dressed two-qubit Pauli error rate of 5.5(1.3) × 10^−3 that can be further decomposed into a single-qubit Pauli error rate of 1.0(0.5) × 10^−3, and an intrinsic two-qubit simultaneous error of 3.5(1.4) × 10^−3. An intrinsic Pauli error rate of 3.5 × 10^−3 corresponds to an average fidelity of 99.72%. We validated the digital error model by looking at patched variations of the random circuit (inset in Fig. 4a), where slices of two-qubit gates have been removed to create patched circuits for which each patch XEB can be verified for a modest computational cost. The total fidelity was then the product of the patch fidelities. Computing XEB over full circuits is, at present, an intractable classical task. We, thus, estimated the fidelity after 32 cycles using the discrete error model, obtaining 0.1%. This increased depth was possible thanks to the substantially reduced errors compared with previous processors. We collected over 70 million sample bit strings for a single circuit at this depth. In Supplementary Information section SC1, we report fidelities for another XEB experiment, SYC-70, done on 70 qubits and with 24 cycles. In Fig. 4b, we verify these extracted fidelities with the Loschmidt echo, where we use the same circuit and its inversion to return to the initial state. We observed good agreement with the XEB experiment and the digital error model.

Finally, we estimated the equivalent classical computational cost of RCS with the tensor network contraction method^{7,11,12,13,14,15,16,18,32,33,34}. With this method, ref. ¹⁵ classically sampled the RCS experiment performed in 2019⁴ in 15 h using 512 GPUs. Reference ¹⁶ also performed this task with a similar cost. Furthermore, another team computed⁷ the corresponding XEB, which confirmed the predictions of ref. ⁴, which is a harder computational task than noisy sampling. Hence, these notable improvements in classical algorithms significantly raised the beyond-classical threshold. For completeness, in Supplementary Information section H, we also study matrix product states, a popular tensor network variational representation of one-dimensional quantum states with limited entanglement^35,36,37. We found that given current supercomputer memory constraints, matrix product states failed to reach the experimental fidelity and offered worse performance than tensor network contraction.

We report improvements in tensor network contraction techniques that resulted in lower estimated computational costs for simulated RCS (Supplementary Information section G). Figure 4c shows the time complexity or number of floating-point operations (FLOPs) (the number of real multiplications and additions) as a function of the number of qubits and cycles required to compute a single amplitude at the output of a random circuit without memory constraints, as we are optimizing only the contraction ordering of the underlying tensor network. The time complexity in this context is rigorously defined as the minimum contraction cost achievable in terms of FLOPs per amplitude¹², which we approximated by optimizing the contraction ordering. This served as a proxy lower bound for the hardness of both sampling and verification. For a fixed number of qubits and increasing number of cycles, there was a notional crossover in the scaling of the time complexity from exponential to linear. Given a noisy experimental set-up, this implies an optimal number of cycles for the trade-off between computational time complexity and fidelity. Beyond the crossover, fidelity decreases faster than the time complexity increases. The crossover number of cycles is consistent with a scaling \(\sqrt{n}\). A line following the functional form \(d=A\sqrt{n}\), where n is the number of qubits, has been added as a guide to the eye. Note that this is related to the number of cycles at which ‘typical’ entanglement is achieved (Supplementary Information section H) and is a stronger requirement than the anti-concentration of the output distribution. For both 67 and 70 qubits, 24 cycles was deep enough to saturate the exponential growth of this time complexity. Figure 4d shows the growth of the time complexity without memory constraints of the largest RCS experiments run over the last few years.

A practical estimate of the computational resources to simulate RCS needs to take into account the finite FLOPS (FLOPs per second) delivery of a supercomputer as well as its memory constraints and other limitations such as finite bandwidth. Table 1 shows estimates of the runtime for the approximate simulation of the largest instances of RCS from refs. ^4,5,6 and the current work when using the state-of-the-art methods discussed in Supplementary Information section G. For these estimates, we considered sampling 1 million uncorrelated bit strings at a fidelity like that of the experiment using the current top-performing supercomputer, Frontier, which performs 1.7 × 10¹⁸ single-precision FLOPS of theoretical peak performance spread across GPUs with 128 GB of RAM each. This required the computation of 10 million approximate probability amplitudes of uncorrelated bit strings to be used in rejection sampling¹⁸. Despite notable progress in achieving tensor contraction algorithms that are embarrassingly parallelizable over each GPU^14,32, we found that this technique broke down for the much deeper SYC-67 circuits with 32 cycles given the tight memory constraints. As a result, an estimated lower bound for the sampling cost became substantially more prohibitive. Assuming a distributed use of all RAM and under the unrealistic assumption of negligible bandwidth constraints, the computational cost was around 1 × 10⁴ years (Table 1). In the untested case in which we expanded working memory to all secondary storage and still ignored the bandwidth, we obtained an estimate of 12 years.

In conclusion, our experiment provides direct insights on how quantum dynamics interacts with noise. The observed phase boundaries lay out quantitative guidance to the regimes where noisy quantum devices can properly leverage their computational power. In addition, we present new RCS results with an estimated fidelity of 1.5 × 10^−3 at 67 qubits and 32 cycles or 880 entanglement gates, corresponding to more than doubling the circuit volume compared to ref. ⁴ for the same fidelity. Global correlations dominate XEB in the weak-noise phase, which protects RCS against ‘spoofing’ attacks, in contrast to boson sampling³⁸, where all known metrics for recent experiments^39,40,41 are dominated by local correlations⁴². Looking forward, despite the success of RCS in quantifying the available coherent resources, finding practical applications for near-term noisy quantum processors still remains an outstanding challenge. Certified randomness generation^43,44 could be a promising candidate for such an application (Supplementary Information section I).