Classical Simulations of Low Magic Quantum Dynamics

TL;DR

This paper introduces classical algorithms for simulating low-magic quantum circuits, especially those with frequent measurements, enabling analysis of larger systems and complex phenomena like measurement-induced phase transitions.

Contribution

The authors develop new classical simulation algorithms tailored for low-magic, measurement-rich quantum circuits, expanding the capacity to study larger and more complex quantum dynamics.

Findings

Simulated all-to-all monitored quantum circuits with low T-gate rates.

Characterized measurement-induced phase transitions in entanglement, purification, and magic.

Demonstrated the utility of algorithms for low-magic, high-entanglement dynamics.

Abstract

We develop classical simulation algorithms for adaptive quantum circuits that produce states with low levels of ``magic'' (i.e., non-stabilizerness). These algorithms are particularly well-suited to circuits with high rates of Pauli measurements, such as those encountered in quantum error correction and monitored quantum circuits. The measurements serve to limit the buildup of magic induced by non-Clifford operations arising from generic noise processes or unitary gates, respectively. Our algorithms also allow a systematic truncation procedure to achieve approximate simulation. To benchmark our approach, we study the dynamics of all-to-all monitored quantum circuits with a sub-extensive rate of T-gates per unit of circuit depth, where we can simulate previously inaccessible system sizes and depths. We characterize measurement-induced phase transitions in the output wavefunction, including in the entanglement, purification, and magic. We outline the utility of our algorithms to simulate dynamics with low magic and high entanglement, complementary to the leading matrix-product state approaches.

Tables

Table 1. Table 1: Summary of efficient simulation regimes for 𝒩 ¯ \overline{\mathcal{N}} from the LRSD algorithm.

Measurement Basis	Gate Type	$T$-gate Rate	Measurement Rate	Scaling	Numerics
$X$-basis	CZ	$p_T=\eta /L^{\beta }$, $\beta \ge 1$	$p_m\gtrsim 0.6$	$\mathrm{poly}\left(L\right)$	Fig. 4
$X$-basis	Random Clifford	$p_T=\eta /L^{\beta }$, $\beta \ge 1$	$p_m\gtrsim 0.6$	$\mathrm{poly}\left(L\right)$	Fig. 18
$Z$-basis	CZ	any $p_T$	any $p_m$	${\left(2L+1\right)}^2$	–

Table 2. Table 2: Critical exponents of Connectivity Transition.

	$p_{Z\vee X}=0$	$p_{Z\vee X}=0.25$	$p_{Z\vee X}=0.5$	$p_{Z\vee X}=1$
$p_m^{cc}$	0.66(1)	0.61(2)	0.61(2)	0.66(2)
$p_c^{ce}$	0.26(2)	0.26(1)	0.25(3)	0.27(2)
${\nu }_c$	2.5(3)	2.3(2)	2.3(2)	2.0(2)
${\tilde{\beta }}_c$	1.80(8)	1.63(5)	1.66(7)	1.48(6)

Equations

$${\rho }_S=2^{r-L}\prod _{i=1}^r\left(\frac{I+g_i}{2}\right)=\frac{1}{2^L}\sum _{g\in \mathcal{S}}g.$$

$$g_i=\underset{j=1}{\overset{L}{⨂}}{P}_{ij},P_{ij}\in \{I,X,Y,Z\}.$$

$$\mathcal{S}(|\psi ⟩)=\{P\in {\mathcal{P}}_L\mid P|\psi ⟩=|\psi ⟩\},$$

$$\text{outerproduct}\psi \psi =\sum _{l\in {\mathcal{L}}_{\mathrm{LRSD}}}{\lambda }_l{\sigma }_l{\rho }_S.$$

$$\text{outerproduct}\psi \psi \to T_i\text{outerproduct}\psi \psi T_i^{†}=\sum _{l\in {\mathcal{L}}_{\mathrm{LRSD}}}{\lambda }_l{T}_i{\sigma }_l{\rho }_S{T}_i^{†}.$$

$$T_i=\alpha I^{\otimes L}+\beta P,P=I^{\otimes (i-1)}\otimes Z\otimes I^{\otimes (L-i)}$$

$$\rho \to \sum _{l\in {\mathcal{L}}_{\mathrm{LRSD}}}{\lambda }_l{\sigma }_l{T}_i{\rho }_S{T}_i^{†}.$$

$${\rho }_S=\frac{1+\bar{P}}{2}{\rho }_S^P,$$

$$T_i\frac{1+\bar{P}}{2}{T}_i^{†},$$

$$T_i{\sigma }_q{T}_i^{†}=({\alpha }^2+{\beta }^2){\sigma }_q+2i\alpha \beta {ϵ}_{qr3}{\sigma }_r.$$

$$\rho =\sum _{l\in {\mathcal{L}}_{\mathrm{LRSD}}}{\lambda }_l{T}_i{\sigma }_l{T}_i^{†}{\rho }_S.$$

$${\sigma }_l\to {\sigma }_l\bar{P}.$$

$$\mathcal{M}(|\psi ⟩)=L-{\log }_2\text{absolutevalue}\mathcal{S}(|\psi ⟩),$$

$${\text{CZ}}_{ij}|s_i{s}_j⟩=(-1{)}^{s_i{s}_j}|s_i{s}_j⟩,$$

$$M_{\alpha ,\pm }^{(j)}|\psi ⟩=\frac{{\Pi }_{\alpha ,\pm }^{(j)}|\psi ⟩}{\Vert {\Pi }_{\alpha ,\pm }^{(j)}|\psi ⟩\Vert },$$

$$T=\left(\begin{array}{cc}1& 0\\ 0& e^{i\pi /4}\end{array}\right).$$

$$|\psi ⟩=\sum _{i=1}^{\chi }{c}_i|{\varphi }_i⟩,$$

$$\rho (t)=\text{outerproduct}\psi (t)\psi (t)$$

$$\mathcal{N}=(2L+1{)}^2+\left|\{{\lambda }_l\}\right|L+\left|\{{\lambda }_l\}\right|$$

$$T^{(j)}\frac{{\Pi }_{z,\pm }^{(j)}|\psi ⟩}{||{\Pi }_{z,\pm }^{(j)}|\psi ⟩||}=e^{i{\varphi }_{\pm }}\frac{{\Pi }_{z,\pm }^{(j)}|\psi ⟩}{||{\Pi }_{z,\pm }^{(j)}|\psi ⟩||},$$

$$S_A=-\mathrm{tr}({\rho }_A\log {\rho }_A).$$

$$\tau \sim L^{z_p},$$

$$\tau L^{-z_p}\sim f_p((p_{m,X}-p_{m,X}^{cp})L^{z_p/{\nu }_p}),$$

$$p_T=\eta /L^{\beta }$$

$$f(L)=a+b{L}^{\gamma },$$

$$\frac{\mathcal{M}(L)}{\log (L)}=\frac{a}{\log (L)}+b.$$

$$\text{outerproduct}\psi \psi =\sum _{l\in {\mathcal{L}}_{\mathrm{LRSD}}}{\lambda }_l{\sigma }_l{\rho }_S,$$

$$\text{outerproduct}\psi \psi \mapsto \sum _{l\in {\mathcal{L}}_{\mathrm{LRSD}}}{\lambda }_l{U}_C{\sigma }_l{U}_C^{†}{U}_C{\rho }_S{U}_C^{†}.$$

$$P=\underset{j=1}{\overset{L}{⨂}}{\sigma }_{ij},{\sigma }_{ij}\in \{I,X,Y,Z\}.$$

$${\rho }_S=\frac{1+\bar{P}}{2}{\rho }_S^P,$$

Full text

I Introduction

Simulating the dynamics of quantum many-body systems is a fundamental challenge, as the Hilbert space grows exponentially with system size, rendering exact state-vector methods intractable for all but the smallest systems. Although tensor network methods such as matrix product states (MPS) can efficiently simulate states with area-law entanglement, they face exponential scaling costs when confronted with the volume-law entanglement characteristic of generic quantum dynamics [36, 44].

Therefore, developing new classical simulation algorithms is crucial for advancing our understanding of key areas like quantum error correction (QEC) and measurement-induced phase transitions (MIPTs). In QEC, robust simulations are necessary to model fault tolerance in noisy, highly entangled circuits [13, 42], while in MIPTs, large-scale simulations are essential for characterizing universal scaling behavior [26, 39, 12, 47]. A powerful tool for studying highly entangled dynamics is the stabilizer formalism, which allows for the efficient simulation of Clifford circuits [16, 1]. However, Clifford circuits are not universal and cannot capture the generic quantum phenomena that arise from non-Clifford operations [10, 9]. Since simulation complexity grows exponentially with the number of non-Clifford gates, there is a strong motivation to develop methods that can efficiently handle circuits that are “close” to the stabilizer manifold [8, 15, 5].

The amount of non-Clifford resources required to prepare a state is known as non-stabilizerness, commonly referred to as “magic”. Magic is zero for stabilizer states and increases monotonically under the application of non-Clifford operations. In recent years, a range of near-stabilizer simulation methods have been developed to simulate quantum circuits with limited non-Clifford resources efficiently. These include Monte Carlo techniques based on stabilizer decompositions [10, 37, 4], mapping to classical statistical mechanics models [41] as well as tensor network methods such as MPS simulations [14] and Clifford-augmented Matrix Product States (CAMPS) [27, 35]. Additionally, hybrid approaches such as stabilizer tensor networks, which combine stabilizer formalism with tensor network contraction and sampling, have been developed  [28]. Collectively, these techniques enable approximate classical simulation of quantum circuits that are not entirely Clifford but remain structured enough to avoid worst-case exponential scaling of simulation complexity with the number of non-Clifford gates.

In the following manuscript, we avoid the sampling overhead of existing near-stabilizer techniques by developing a near-Clifford simulator that uses an exact representation of the state. Instead of sampling over Clifford circuits, our approach updates the stabilizer decomposition dynamically during the circuit evolution. Our approach is designed for dynamic measurement and mid-circuit branching, such as those encountered in quantum error correction and MIPTs [39, 25]. Existing stabilizer rank techniques do not generalize easily to dynamic, measurement-driven evolutions. Our approach can also simulate states with volume law entanglement, for a low number of $T$-gates, unlike MPS approaches. The Clifford+$T$ gate set constitutes the minimal extension of stabilizer circuits to universality, making it the natural and minimal framework for studying the onset of non-stabilizer physics in monitored quantum circuits.

To benchmark our algorithm in a non-trivial setting, we consider random circuit models with two-qubit gates and single-site measurements, but the gates are allowed to act on arbitrary pairs of qubits in the system that are chosen at random. Such “all-to-all” models can exhibit MIPTs, but the usual notion of an area-to-volume law entanglement transition that arises in spatially local models is revised [31]. These models are particularly difficult for MPS-based simulations due to the lack of locality and are therefore useful as a general test case for our purposes.

The remainder of this paper is organized as follows. In Section II, we introduce our algorithm for updating near-stabilizer states under the action of $T$-gates. We also define the measure of magic, stabilizer nullity, which we study in this work. In Section III, we introduce the models of monitored circuits, including entangling gates and protocols for the measurements. In Section IV.2, we test the accuracy of our algorithm, which simulates the full density matrix, by varying a cutoff parameter that controls the accuracy of the algorithm. In Section IV.3, we examine the space complexity of these circuits, namely how the amount of memory to store the state scales with the density of $T$-gates. In Section V, we study the MIPTs in these models, focusing on the magic and purification transitions. In Appendix A, we provide more details on the algorithms used to study these circuits. These include the evolution of the density matrix of the generic state $\ket{\psi}$, $\rho=\outerproduct{\psi}{\psi}$, under measurements and Clifford unitaries, computing the Rényi entanglement entropies, and Bell sampling. In Appendix B, we examine entanglement transitions in more detail. Namely, we locate cluster and all-to-all entanglement transitions. We also provide more details on locating the purification transition. In Appendix D, we include details on the procedure used to locate the magic transition.

II Near stabilizer simulation algorithm

To set the stage, we begin by reviewing the stabilizer formalism that we will be perturbing with a dilute amount of $T$-gates.

II.1 Stabilizer Formalism

The density matrix of a general mixed stabilizer state on $L$ qubits can be written as a sum over the elements of its stabilizer group $\mathcal{S}$ [1]:

[TABLE]

Here, $\mathcal{S}$ is an abelian subgroup of the $L$-qubit Pauli group $\mathcal{P}_{L}$ that does not contain $-I$. The generators $g_{i}$ of $\mathcal{S}$ are tensor products of Pauli matrices ($X$, $Y$, $Z$),

[TABLE]

An operator $g\in\mathcal{S}$ is said to stabilize $\rho_{S}$ if $g\rho_{S}g^{\dagger}=\rho_{S}.$ For a set of $r$ independent generators, the stabilized subspace has dimension $2^{L-r}$. Pure stabilizer states correspond to a full set of $r=L$ generators, while mixed stabilizer states have $r<L$.

II.2 Non-Clifford Evolution

Given an arbitrary pure quantum state $\ket{\psi}$, the stabilizer group of $\ket{\psi}$ is the set of Pauli operators that have eigenvalue one on that state

[TABLE]

where $\mathcal{P}_{L}$ is the Pauli group on $L$ qubits and $\mathcal{S}$ is a subgroup of $\mathcal{P}_{L}$. We define the low-rank stabilizer decomposition (LRSD) of $\ket{\psi}$ as the representation

[TABLE]

Here, $\rho_{S}$ is a (possibly mixed) stabilizer state that is defined solely in terms of $S$ through Eq. (1). The coefficients $\lambda_{l}=\matrixelement{\psi}{\sigma_{l}}{\psi}$ are real numbers and the “logical” operators $\sigma_{l}$ are elements of the Pauli group that commute with $S$. The set of all logical operators modulo phases $\mathcal{L}_{\mathrm{LRSD}}$ has size $4^{k}$ for a stabilizer group of size $2^{L-k}$ [1]. If $\rho_{S}$ is pure, then $\ket{\psi}$ contains no magic, and Eq. (4) consists of a single term: specifically, $\lambda_{l}=1$ and $\sigma_{l}=I$. In this case, $\ket{\psi}$ is a stabilizer state. More generally, for a broad class of near-stabilizer states with small $k$, the number of logical operators needed in the LRSD remains small, allowing for an efficient classical representation. This approach isolates the non-stabilizer content of the state into a compact, low-dimensional structure, the logical operator set, while leveraging the efficient classical simulability of stabilizer states. This representation builds on the stabilizer rank framework [8, 9], which decomposes near-Clifford states into sums of stabilizer states. The specific structured form used here, expressing a state in terms of the logical operators of a single stabilizer code, was introduced in Ref. [18]. We note that while Eq. (4) is written for the pure-state projector $\outerproduct{\psi}{\psi}$ of a single measurement trajectory, the LRSD naturally generalizes to mixed states $\rho=\sum_{l}\lambda_{l}\sigma_{l}\rho_{S}$ with real coefficients $\lambda_{l}$.

We now present how to update the LRSD representation of the state under the action of $T$-gates. For a $T$-gate on the $i^{th}$ qubit, $T_{i}$, the state evolves as:

[TABLE]

Up to an overall phase, the $T$-gate can be rewritten in the following form:

[TABLE]

Here, $\alpha=\cos\left(\frac{\pi}{8}\right)$ and $\beta=-i\sin\left(\frac{\pi}{8}\right)$. The overall phase is irrelevant because it is canceled when the density matrix is acted upon by conjugation. This form of the decomposition is useful because it separates the action of the $T$-gate into parts that either commute or anticommute with Pauli operators. This simplifies the update rules. In each update with $T_{i}$, the complexity of the LRSD will depend on the commutation relations between $P$ with the stabilizers and logicals of $\rho_{S}.$ Specifically, there are four cases to consider when evolving the state under $T$-gates. In the following, we use $\mathcal{S}$ to denote the stabilizer group of $\rho_{S}$, and $[P,\mathcal{S}]=0$ to mean that $P$ commutes with every element of $\mathcal{S}$ and likewise for $\mathcal{L}_{\mathrm{LRSD}}$.

• •

Case I: If $P\in\mathcal{S}$, then $[P,\mathcal{L}_{\mathrm{LRSD}}]=0$ and $[P,\mathcal{S}]=0$. In this scenario, $\rho$ remains unchanged.

• •

Case II: If $[P,\mathcal{S}]\neq 0$ but $[P,\mathcal{L}_{\mathrm{LRSD}}]=0$, then

[TABLE]

We perform the stabilizer decomposition using Algorithm 3:

[TABLE]

where $[P,\rho_{S}^{P}]=0$ . The action of the $T$-gate on the Pauli string $\bar{P}$ is given by

[TABLE]

which can be deduced from the rule:

[TABLE]

Here, $\epsilon_{qr3}$ is the Levi-Civita symbol and $\sigma_{q}$ is a single-qubit Pauli matrix, where each index runs from $1$ to $3$.

• •

Case III: If $[P,\mathcal{S}]=0$ but $[P,\mathcal{L}_{\mathrm{LRSD}}]\neq 0$, then $P$ commutes with all the stabilizers but is not itself in the stabilizer group. In this case:

[TABLE]

Each term in the tensor product of $\sigma_{l}$ is updated using Eq. (10).

• •

Case IV: If $[P,\mathcal{L}_{\mathrm{LRSD}}]\neq 0$ and $[P,\mathcal{S}]\neq 0$, the anticommuting logicals can be updated so that $[P,\mathcal{L}_{\mathrm{LRSD}}]=0$. If $[\sigma_{l},P]\neq 0$, then

[TABLE]

This transformation leaves the state invariant as $\bar{P}\rho_{S}=\rho_{S},\forall\bar{P}\in\mathcal{S}$. After this update, $[P,\mathcal{L}_{\mathrm{LRSD}}]=0$, which reduces to Case II.

The treatment of Clifford unitaries and Pauli measurements appears in Appendix A. We can now provide some intuition for how the complexity of the classical simulation evolves. As visually summarized in Fig. 1, each $T$-gate application branches the logical operators into new paths, while subsequent measurements can collapse or eliminate them. Each $T$-gate application can increase the number of logicals $|\mathcal{L}_{\mathrm{LRSD}}|$. Specifically, in Cases II and IV, the number of logicals can triple, while in Case III, they can double. Conversely, each measurement can reduce the number of logicals, depending on the commutation relationships. If $[P,\mathcal{L}_{\mathrm{LRSD}}]=0$, or if $[P,\mathcal{L}_{\mathrm{LRSD}}]\neq 0$ and $[P,\mathcal{S}]\neq 0$, the number of logicals remains the same. However, if $[P,\mathcal{L}_{\mathrm{LRSD}}]\neq 0$ and $[P,\mathcal{S}]=0$, the number of logicals can decrease by up to a factor of three. The detailed procedure is provided in Algorithm 3.

Introducing $T$-gates takes us out of the stabilizer manifold. Here, we quantify the non-stabilizerness of a state through its “distance” from the set of stabilizer states [10, 19, 24] and is commonly referred to as “magic”. For a pure stabilizer state, the magic is by definition zero. For a non-stabilizer state, the magic quantifies the degree to which the state deviates from the stabilizer set, and serves as a resource-theoretic measure of the non-Clifford operations required for its preparation. While high magic often correlates with greater circuit complexity [9], the experimental difficulty of preparing a given non-stabilizer state may vary depending on the physical platform and available operations[6, 46].

We now define a measure of magic called the stabilizer nullity [7]. If $\ket{\psi}$ is an $L$ qubit state, the stabilizer group $\mathcal{S}$, of $\ket{\psi}$ is defined in Eq. (3) The stabilizer nullity of $\ket{\psi}$, $\mathcal{M}(\ket{\psi})$, is

[TABLE]

where $\absolutevalue{\mathcal{S}(\ket{\psi})}$ is the size of the stabilizers. For a stabilizer state, $\mathcal{M}=0.$ For non-stabilizer states, $\mathcal{M}>0.$ The stabilizer nullity measures the number of non-stabilizer degrees of freedom in the state. It serves as an efficiently computable upper bound on more general measures of magic, such as stabilizer rank or robustness of magic [9, 20, 38, 40].

III Circuit Model

In all-to-all random circuit models, there is no constraint of spatial locality. Part of the motivation of studying these models is that a graphical description of the MIPT in these circuits leads to a precise connection between the MIPT and mean-field percolation [45]. This is analogous to the connection between two-dimensional classical percolation and the measurement-driven entanglement transition in one spatial dimension [39]. Furthermore, all-to-all models are directly relevant to various experimental platforms. For example, proposals for scaling up quantum computers often involve quantum network architectures, where locally interacting nodes are connected via non-local photonic interconnects, leading to effectively all-to-all connectivity at large scales [22]. At the hardware level, trapped ion quantum computers already realize all-to-all gate connectivity among physical qubits at intermediate scales [33, 29]. These experimental developments further underscore the importance of understanding the fundamental features of MIPTs in nonlocal circuit architectures. We now introduce the circuit protocols considered in this paper. We focus on an all-to-all geometry where the two-site unitary can be applied to any pair of qubits, which is an effective 0D system. Specifically, we consider single-pair all-to-all models, where in a single time step, at most one gate couples a pair of qubits, as shown in Fig. 2. Each time step consists of three probabilistic stages: Clifford unitaries, non-Clifford unitaries, and measurements.

The single-pair all-to-all model, as shown in Fig. 2, is adapted from Ref. [45]. In the first stage, we apply a control-$Z$ (CZ) gate, as defined in Eq. (14) with a probability of $1/2$ to a specific pair of qubits 111This choice ensures that the circuit is in a nontrivial entangling regime: choosing the probability of a CZ gate arbitrarily as $p_{u}$, if $p_{u}$ is too small it drives the system into a trivial area-law entangled phase for all non-zero measurement rates.. This pair is chosen uniformly from all combinations of pairs. After the Clifford unitary gate, we apply a $T$-gate with a probability of $p_{T}$ to a randomly chosen qubit. Note that $p_{T}$ is therefore a per-timestep rate, not a per-qubit rate. The $T$-gate is the simplest non-Clifford gate that, combined with Clifford operations, promotes the gate set to universality. Our model therefore represents the minimal circuit architecture needed to capture the departure from the stabilizer manifold. Finally, we perform a projective local measurement with a finite probability onto a specific qubit chosen randomly.

In particular, we study two distinct models of quantum circuits that are distinguished by the basis in which projective measurements are performed: the $Z$ or $X$ basis. We denote the probability to make a measurement along the $z$-axis as $p_{m,Z}$ and the $x$-axis as $p_{m,X}$. As with $p_{T}$, these are per-timestep rates, not per-qubit rates. In the $Z$ basis model, measurements are diagonal in the computational basis, which allows for an efficient classical simulation for all rates of $T$-gates to remain. However, computing the amount of non-stabilizerness (i.e., the magic) of these states remains a computational challenge, and one of the algorithms we develop allows for its efficient computation (technical details are in Appendix A.3). In contrast, the $X$-basis measurements act in a rotated basis, which can generate additional entanglement and simulation difficulty. Our newly developed algorithm (in Appendix A.1) is able to still efficiently simulate the circuit evolution in the presence of a sub-extensive set of $T$-gates (i.e. $p_{T}\sim 1/L$ or $\sim 1/L^{2}$).

The entangling unitary in our model is the CZ gate in the Clifford group, acting on two qubits $i$ and $j$ as

[TABLE]

where $s_{i},s_{j}\in\{0,1\}$ are the computational basis states. The CZ gate is a Clifford gate, meaning that it maps stabilizer states onto stabilizer states. That is, the set of stabilizing Pauli operators for the output state can be obtained by conjugating the stabilizers of the input state with the CZ gate. This transformation preserves the stabilizer structure and can be tracked efficiently using stabilizer formalism.

We consider a single-qubit projective measurement $M_{\alpha,\pm}^{(j)}$ along the $\alpha$ direction in the mid-circuit as

[TABLE]

where $\alpha$ is the axis of the measurement, and $\Pi_{\alpha,\pm}^{(j)}=\outerproduct{\pm}{\pm}_{\alpha}$ project qubit $j$ to the $\pm\alpha$ axis. Here, we consider two bases of measurements, either along $z$ axis ($\alpha=z$), or along $x$ axis ($\alpha=x$), respectively. For the measurement along the $z$-axis, we additionally rotate the measurement outcome for each qubit from $\ket{\pm z}$ to $\ket{+x}$ to keep them entangled under the unitary gates. Our architecture is adapted from the solvable circuit model introduced in [45], which demonstrates a novel phase transition between a fully-entangled phase and a separable phase composed of finite subsystems.

The second type of unitary gate is the single-qubit $T$-gate, which introduces a non-Clifford $\pi/4$ phase to the $\ket{-z}$, defined as

[TABLE]

Applying a $T$-gate to a stabilizer state $\ket{\psi}$ results in a superposition of at most two stabilizer states, $\ket{\psi^{\prime}}=c_{+}\ket{\psi_{+}}+c_{-}\ket{\psi_{-}}$ as depicted in Fig. 2. If $T$ anticommutes with at least one element of the stabilizer group, the stabilizer state is taken out of the stabilizer manifold. This branching increases the complexity of near-stabilizer simulators, which use the stabilizer formalism. A generic pure quantum state $\ket{\psi}$ can be written in terms of stabilizer states,

[TABLE]

where the stabilizer rank is $\chi(\ket{\psi})$, $c_{i}\in\mathbb{C}$, and $\ket{\phi_{i}}$ is a stabilizer state. The minimum number $\chi$ for which such a decomposition exists is the stabilizer rank.

IV Complexity of Simulation

IV.1 Simulation Task

The computational task studied in this work is the trajectory-resolved simulation of monitored quantum circuits. Specifically, for a single measurement trajectory of depth $t\sim L^{2}$, we store and evolve an explicit representation of the full density matrix

[TABLE]

under probabilistic unitaries (CZ), $T$-gates, and projective measurements in either the X or Z basis. The goal is to maintain an exact truncation-controlled representation of $\rho(t)$, thus enabling computation of nonlinear observables such as entanglement entropies, purification times, and stabilizer nullity.

This is the number of entries, or bits and complex numbers, required to specify the state. The stabilizer state is represented by a tableau of size $(2L+1)\times(2L+1)$. Each Pauli operator $\sigma_{l}$ is represented by a $2L$ binary vector, as defined in Section IV. Each phase $\lambda_{l}$ is a complex number. Each logical Pauli operator $\sigma_{l}$ is specified by an $L$-component vector with entries in $\{0,1,2,3\}$, labeling $I,X,Y$, and $Z$, respectively.

These considerations mean that

[TABLE]

entries are required to specify the state. The $(2L+1)^{2}$ factor is the number of entries required to specify $\rho_{S}.$ The $|\{\lambda_{l}\}|L$ factor is the number of entries required to store $\sigma_{l}.$ The $|\{\lambda_{l}\}|$ is the number of entries required to store the complex phases. Although the Pauli indexing is not unique, one could also use a length-$2L$ binary vector to specify $\sigma_{l}$, for example. However, all reasonable choices should give the same asymptotic scaling $\sim L^{2}+\left|\{\lambda_{l}\}\right|L+\left|\{\lambda_{l}\}\right|$. By contrast, a full state-vector description demands $2^{L}$ complex amplitudes, or complex double-precision numbers. Moreover, the state-vector stores exclusively complex numbers, whereas our representation mixes complex numbers, integers, and Booleans, so the count above is a conservative estimate of the true memory advantage of our approach.

We say that the circuit ensemble is efficiently simulable if, for trajectories drawn from the ensemble, the scaling of $\mathcal{N}(L)$ with system size $L$ is polynomial rather than exponential. The parameter regimes in which this condition is satisfied, along with numerical evidence and a detailed analysis of worst-case versus average-case scaling, are presented in Sec. IV.3.

IV.2 Magic Complexity

In this section, we take the limit of $Z$-basis only measurements and examine the single-pair all-to-all model. Since $T$ and CZ are both diagonal in the $Z$-basis, $[T,\mathrm{CZ}]=0$. It can also be shown that

[TABLE]

where $\phi_{+}=0$ and $\phi_{-}=\pi/2$ and we can ignore the global phase factor. Since $T$ and $CZ$ commute, and there are no other gates that mix $Z$-basis states, a $T$-gate on qubit $i$ can be commuted forward in time until it meets a $Z$-measurement on the same qubit, at which point its effect is absorbed into an irrelevant global phase. Suppose a $T$-gate acts on qubit $i$ at time $t$, and the first $Z$-basis measurement after the $T$-gate occurs on qubit $i$ at time $t^{\prime}>t.$ Then, if the last $Z$-basis measurement on qubit $i$ occurs at time $t^{\prime}$, then all the $T$-gates on qubit $i$ at $t\leq t^{\prime}$ are removed from the circuit.

To test the accuracy of the algorithm, with a truncation error threshold $\epsilon$, we evolve the LRSD under the circuit, truncated at each timestep to keep only the coefficients $\lambda_{l}$ such that $|\lambda_{l}|>\epsilon$. If there are no remaining logicals with $|\lambda_{l}|>\epsilon$, the trajectory is invalid and not considered. We pick $\epsilon$ sufficiently small so that this issue does not arise.

In Fig. 3, we plot the dependence of the total number of logicals, $\overline{|\{\lambda_{l}\}|}$ on $\epsilon$ for varying $p_{m,Z}$. The $\overline{\cdots}$ denotes an average over circuit trajectories and $|\{\cdots\}|$ is the set size. We fix the location of the unitary gates and measurement locations to ensure $\overline{\left|\{\lambda_{l}\}\right|}(\epsilon)$ depends only on the truncation error and not the randomness of the circuit-to-circuit fluctuations. As shown in Fig. 3, as the cutoff is varied from $\epsilon=0$ to $\epsilon\approx 10^{-2}$, the number of logicals is unchanged, which demonstrates that for this given circuit, this is the required number of logicals. In contrast, further increasing the cutoff $\epsilon$ to approach 1, $|\{\lambda_{l}\}|$ is reduced by roughly two orders of magnitude, or a factor of 100, in the effective Hilbert space dimension, shown in Fig. 3(a). This greatly reduces the computational cost of Bell sampling at the price of a less accurate wavefunction. We also clearly see the role of $p_{m,Z}$ measurements, which reduces the number of required logicals (a) and we find that (b) the stabilizer nullity grows by no more than five percent across the entire measurement rate window $p_{m,Z}\in\{0.6,0.7,0.8,0.9\}$ as compared to the zero truncation curves (dashed lines). These observations demonstrate that non-trivial observables remain robust to moderate truncation of the LRSD. This allows us to introduce $\epsilon$ as an effective cutoff scheme to reduce simulation time while preserving quantitative accuracy. We note that when operating within the area-law magic phase where the number of logicals scales polynomially, $\epsilon$ can be safely set to machine precision for exact simulation, as was done for our space complexity analysis.

Following the Algorithm in Sec II.2 for applying $T$-gates, each $T$-gate will at most triple $\overline{\left|\{\lambda_{l}\}\right|}.$ This means, in the worst case scenario, $\overline{\left|\{\lambda_{l}\}\right|}=3^{|T|}$, where $|T|$ is the number of $T$-gates in the circuit. In this section, we consider a model where a single layer of $T$-gates is applied at the final timestep of the Clifford circuit. That is, we sweep from $i=1$ to $i=L$ once on a Clifford state. This means $\overline{\left|\{\lambda_{l}\}\right|}$ could, at worst, grow to $3^{L}$. The origin of this $3^{L}$ scaling comes from the fact that each of the $L$ $T$-gates can, in the worst case, triple the number of terms in the LRSD. The rapid growth in $\overline{\left|\{\lambda_{l}\}\right|}$ imposes a challenge in the classical simulation. Therefore, only the state within the area law magic phase can be efficiently simulated using LRSD. In the area law magic phase, we find that the number of logicals, when averaged over circuits, scales with system size as a power law, as shown in Fig. 4.

IV.3 Space Complexity

We now turn to the space complexity in our new simulation algorithm. To quantify the space complexity, we use the metric defined in Eq. (19). To emphasize the classical difficulty of simulating these circuits using the LRSD algorithm versus brute force state-vector methods, we now focus on the circuit-averaged number of entries $\overline{\mathcal{N}}$. The model we examine is the single-pair all-to-all model with $X$-basis measurements, as shown in Fig. 2. To demonstrate robustness to sufficiently generic circuits, we test our LRSD algorithm in Appendix C using completely random Clifford gates in place of CZ gates.

We present the average number of entries $\overline{\mathcal{N}}$ as a function of $L$ in Fig. 4(a), for $p_{m,X}$ varied from $0.6$ to $0.9$ in increments of $0.1$ at fixed $p_{T}=1/L$. This parameter regime lies within the area-law magic phase shown in Fig. 6. Due to the exponential growth of the logical space with the number of $T$-gates, the numerical simulation is constrained to approximately 100 to 200 $T$-gates at low measurement rates: at $p_{T}=1/L$ the total $T$-gate count is $2L^{2}p_{T}=2L$, placing the system sizes shown in Fig. 4 within this window. We provide a comparison of $\overline{\mathcal{N}}=2^{L}$ state-vector methods (dashed black line) with the LRSD (solid blue lines). We find that $\overline{\mathcal{N}}$ exhibits power law growth with some exponent $\alpha$, $\overline{\mathcal{N}}=\beta+L^{\alpha}$, with $\alpha\approx 2$. This is much slower than the exponential growth of state-vector methods. In the CAMPS representation, the computational cost remains polynomial for Clifford circuits and can often remain efficient when a Clifford circuit is “doped” with up to $\mathcal{O}(L)$ $T$-gates, particularly when these are distributed favorably across the system. In this regime, the disentangling procedure can frequently absorb the non-Clifford action without a significant increase in bond dimension [27]. Beyond this point—when the number of $T$-gates exceeds $\sim L$ or their placement generates unavoidable entanglement—the bond dimension typically grows rapidly, leading to exponential cost. For sufficiently sparse $T$-gate densities, such as rates scaling as $p_{T}$ $\sim 1/L^{2}$, the simulations remain tractable, and similar behavior is expected for small prefactors $\eta$ in the scaling form $p_{T}$ $\sim\eta/L$. Within the area-law magic phase, the LRSD algorithm exhibits robust power-law scaling with system size $L$ (see Fig. 4), which plausibly provides an upper bound on the scaling with $T$-gate count at fixed system size.

We place our results in the context of prior work on classical simulation of Clifford circuits with a $T$-gate count $t$ proportional to the system size $L$. Prior works have demonstrated that circuits with $t\sim L$ $T$-gates can admit efficient simulation. It was shown [34] that the cost of exactly computing a chosen output probability scales as $\mathcal{O}(2^{t-r}\cdot t)$, where $r$ is a circuit-dependent quantity, termed the projector rank, that quantifies the number of linearly independent parity constraints induced by the $T$-gates. For random circuits with many Clifford gates between each $T$-gate, $r$ concentrates near its maximum $\min(t,n-w)$ (where $n$ is the number of qubits and $w$ the number of measured qubits), so the effective exponent $t-r$ remains small, enabling efficient simulation. Recent work [32] has introduced a stabilizer tensor network method with magic state injection augmented stabilizer tensor network(MAST) and reported that local expectation values of random $T$-doped Clifford circuits can be computed in $\mathrm{Poly}(L)$ time when $t\lesssim L$, with a sharp phase transition in simulation cost beyond this regime. We note, however, that the $\mathrm{Poly}(L)$ cost in Ref. [32] applies to local expectation values; efficient sampling from the output distribution is explicitly not claimed.

In light of these results, we emphasize that our contribution is not the observation that certain $t=\mathcal{O}(L)$ regimes can be tractable, but rather the identification and characterization of an average-case polynomial scaling regime for trajectory-resolved full density matrix evolution in measurement-driven ensembles. Unlike the approaches above, which target output probabilities [34] or local expectation values [32], the LRSD maintains the complete density matrix $\rho(t)$, enabling the computation of nonlinear observables such as entanglement entropies and stabilizer nullity. This represents a fundamentally different and more demanding simulation objective.

To rigorously justify the average-case polynomial efficiency of the LRSD algorithm, it is critical to address the potential for a “hard trajectory” problem. This scenario is common in near-Clifford simulators where a small fraction of quantum trajectories incur exponential simulation costs, thereby dominating the ensemble average. As defined in Eq. 19, the space complexity is determined by the number of logical operators $|\{\lambda_{\ell}\}|$. While the worst case cost for a single trajectory is theoretically bounded by $3^{|T|}$, our numerics confirm that the ensemble is not bottlenecked by rare, exponentially costly instances. As detailed in Appendix C, the trajectory-resolved probability distribution of $\mathcal{N}$ is tightly concentrated around its mean, with a rapidly decaying right-hand tail. This rapid supression of high-cost trajectories ensures that the polynomial scaling of the ensemble-averaged complexity $\mathcal{N}$ in Fig. 4 is accurately representative of a typical circuit realization. Consequently, the simulation remains efficient in the area-law magic regime without being derailed by statistical outliers.

We note that the $X$-basis model—initialization in $\ket{+}^{\otimes L}$, evolution under diagonal unitaries (CZ and $T$), and measurement in the $X$-basis—is structurally related to Instantaneous Quantum Polynomial-time (IQP) circuits [11], for which classical sampling from the output distribution is believed to be intractable under standard complexity-theoretic assumptions. Since maintaining the full density matrix is a strictly more demanding computational task than sampling, worst-case hardness of the LRSD simulation in this circuit class is expected. This complexity-theoretic expectation is consistent with the $3^{|T|}$ worst-case bound discussed above. Crucially, the mid-circuit projective measurements—absent in standard IQP circuits—provide the mechanism that collapses the LRSD and concentrates the trajectory-resolved cost around its polynomial mean, enabling average-case efficient simulation to coexist with worst-case hardness.

We summarize the efficient simulation regimes in Table 1. For $X$-basis measurements with CZ gates, $\overline{\mathcal{N}}$ scales as $\mathrm{poly}(L)$, as demonstrated numerically in Fig. 4. To verify that this scaling is robust to circuit structure rather than an artifact of the CZ gate, we also test $X$-basis measurements with random Clifford gates replacing CZ (Fig. 18 in Appendix C), finding the same polynomial scaling. For the $Z$-basis model with CZ gates, $[T,\mathrm{CZ}]=0$ and $T$-gates commute through $Z$-measurements (Sec. IV.2), so the LRSD evolution is purely Clifford and $\overline{\mathcal{N}}$ reduces to the stabilizer tableau size $(2L+1)^{2}$; no numerical verification is needed.

V Measurement-induced phase transitions

There are different patterns of multipartite entanglement in the all-to-all steady state, which can be distinguished by various metrics. In monitored circuits, three distinct entanglement phase transitions that can occur as the measurement rate $p_{m}$ is varied: (i) an entanglement transition at $p_{ce}$ which separates the volume- and area-law scaling, (ii) a connectivity (percolation) transition at $p_{cc}$ that fully disentangles the state, and (iii) a purification transition at $p_{cp}$ where a maximally mixed ancilla purifies in a time that becomes system size independent. Typically $p_{cp}\leq p_{ce}\leq p_{cc}$.

V.1 Purification Transition

Here, we focus on the purification transition in the single-pair all-to-all model with $X$-basis measurements and provide numerical evidence for its critical point and scaling. We defer detailed definitions and diagnostics for the connectivity and entanglement transitions to Appendix B. Typically $p_{cp}\leq p_{ce}\leq p_{cc}$. Numerically, we find that $p_{ce}=p_{cp}=0.27(1)$ and $p_{ce}<p_{cc}=0.66$.

The dynamical purification phase transition can be efficiently probed by studying how the system preserves entanglement over time with a single reference qubit. After preparing $\ket{\psi_{o}}=\ket{+x}^{\otimes L}$, the reference qubit is entangled to a randomly selected system qubit to form a Bell pair. The entangling operation is followed by a scrambling stage, which consists of $\sqrt{10}L$ random two-qubit Clifford gates applied to random system qubit pairs. We choose $\sqrt{10}L$ gates because the system’s all-to-all entanglement, Eq. (79) in Appendix B.2, is saturated at this point. This ensures that the ancilla is maximally entangled with all the qubits in the system, and not just a single qubit. This reduces the finite-size effects. We evolve the circuit up to a final time $t_{f}=2L^{2}$ in the single-pair all-to-all model with $X$-basis measurements. This is shown in Fig. 2(a). This timescale is chosen to ensure that the ancilla entanglement entropy has entered its long-time regime, where it decays exponentially with time for all system sizes and measurement rates. The decay of $\overline{S}_{Q}$ with $t$, for $0<t<2L^{2}$ is shown in Fig. 14 at $p_{m,X}=p_{m,X}^{cp}$.

When we use $X$-basis measurements, Eq. (20) no longer applies. Recall that we used Eq. (20) in Sec.IV.2 to argue that the $T$-gates constitute a single layer of single-qubit unitaries at the end of the Clifford circuit. The $T$ gates are commuted to the end of the circuit using the fact that they are diagonal in the $Z$ basis and commute with both the $CZ$ gates and $Z$-basis measurements. A final layer of single-qubit unitaries cannot alter the entanglement of a pure state $\ket{\psi}$. In this case, $\ket{\psi}$ is a stabilizer state, which has a degenerate entanglement spectrum. However, since Eq.  (20) no longer applies, the argument at the beginning of Section IV.2 is invalid. The $T$-gates are therefore interspersed throughout the circuit, which means the entanglement is no longer the same as the stabilizer state. Instead, the final state is a generic non-stabilizer state without a degenerate entanglement spectrum. We must therefore use the LRSD rather than the stabilizer formalism to simulate the state and calculate the entanglement.

We track the entanglement entropy between the system and the ancilla qubit, $\overline{S_{Q}(t)}$. $S_{Q}(t)$ is the von Neumann entanglement entropy, $S_{A}$, for a bipartition $A$ and $B$ of the system:

[TABLE]

The purification time $\tau$ can be estimated from the decay of the ancilla entanglement, namely, $\overline{S_{Q}(t)}\sim\exp(-t/\tau)$ (see Appendix B.3.1 for details on the fitting procedure used to obtain $\tau$). Therefore, by fitting the late-time decay and plotting $\tau$ as a function of system size, we can identify the purification transition, as shown in Fig. 5(b) and (c), leading to a phase diagram in Fig. 5(a) in the limit of $X$-basis only measurements, as a function of $p_{m,X}$ and $p_{T}$. Recall that $p_{T}=\eta/L^{\beta}$ in Eq. (24) is the probability of injecting a single $T$-gate per timestep at a randomly chosen qubit (Sec. V); we choose $\eta=1$ and $\beta=2$ in Eq. (24). Each timestep contains $\eta/L^{\beta}$ $T$-gates on average, and accounting for $2L^{2}$ timesteps, this means that there are a total of $2L^{2}p_{T}=2\eta L^{2-\beta}$ $T$-gates in the circuit. For $L=64$, the number of $T$-gates is sufficiently low so that it fits in the constraints of our LRSD simulation method. In Section A.1, we explain how to perform Clifford operations on the LRSD.

In the phase diagram Fig. 5, we find that the critical measurement rate $p_{m}^{cp}$ is unaffected by the rate of $T$-gates within error bars. This is consistent with prior work [14], which has shown that the entanglement transition critical point $p_{c}^{ce}$ is insensitive to the rate of $T$-gates for $\beta\leq 1$, and we expect $p_{c}^{ce}=p_{c}^{cp}$ in this model.

In Fig. 5(b), we show the decay time as a function of system size. At small measurement rates $p_{m}<0.25$, the initial state takes an exponentially long time to purify, indicating a mixed phase; whereas at large measurement rates $p_{m}>0.25$, the decay time saturates to an $L$-independent value, indicating a pure state. The critical measurement rate $p_{m,X}^{cp}=0.26(1)$ can thus be identified by locating the inflection point of the decay time, where the decay time grows algebraically with system size, i.e.,

[TABLE]

where $z_{p}$ is the dynamical critical exponent for the purification transition. We find $z_{p}=0.22(2)$, which is consistent with the value found in Ref. [39] for the Hartley entropy of Haar random circuits with measurements. With the estimate of the critical measurement rate $p_{m,X}^{cp}$, and dynamical exponent $z_{p}$, we can also perform the finite size scaling of $\tau L^{-z_{p}}$ as a function of $p_{m}$ as shown in Fig. 5(c) using the following ansatz,

[TABLE]

where $f_{p}(x)$ is a universal scaling function and $\nu_{p}$ is the critical exponent for the correlation length in the purification transition. The estimated critical exponents are $\nu_{p}=0.45(5)$ and $z_{p}=0.22(2)$ with the data collapse as shown in the inset of Fig. 5(c). Our estimate of the critical exponents is consistent with Ref. [33] with long-range $XX$ gates, where $XX=\text{exp}(i\frac{\pi}{4}X\otimes X)$, implying the same universality class as mean-field percolation.

For the $n=0$ Rényi entropy of Haar random circuits with measurements, also known as the Hartley entropy, there is an exact mapping to a percolation problem in the circuit geometry. In an all-to-all geometry, this mapping predicts $z_{p}=1/5$ [31].

Our primary focus is on the purification transition, which is about temporal correlations. For completeness, however, we also locate the entanglement transition, which is about spatial correlations. In the limit of no $T$-gates, we find in Appendix B.1 that $p_{ce}=0.27(1)$, coinciding with the entanglement transition. We use an all-to-all definition of the entanglement suited for 0D systems.

The connectivity transition is analytically solvable for Haar random circuits in 1+1 dimensions, $p_{cc}=1/2$, due to the mapping onto 2D percolation [21]. However, this mapping no longer accurately describes the connectivity properties of the state. Instead, the geometric structure of the entanglement correlations in the system is completely contained in the graph state $\ket{G}.$ The graph $G$ describes clusters of entangled qubits, which we use as the basis of our percolation analysis in Appendix B.1.

V.2 Magic Transition

Here, we analyze a magic transition in the single-pair all-to-all model with only $Z$-basis measurements. Following the argument around Eq. (20), $Z$-basis measurements constrain the placement of $T$-gates in a structured way. This leads to a density matrix representation that is amenable to Bell sampling, as discussed in Appendix A.3. We then apply the algorithms developed in Appendix A.3 to efficiently compute the magic. The phase transition is between a (sub)-extensive and area law scaling of magic controlled by the rate of measurements $p_{m,Z}$ and probability of $T$-gates $p_{T}$. To find a phase transition in the magic, we study states that interpolate between states produced by Clifford and non-Clifford circuits in a controlled way. We set the probability to sample $T$ gates as

[TABLE]

Previous work in a bricklayer geometry has shown that if the density of non‑Clifford $T$ gates is kept constant as the system size $L$ grows, $\beta=0$, the total magic follows a volume law for any measurement rate $p<1$ [14]. Magic is additive and can be produced locally, whereas entanglement can only be generated by gates whose support straddles the boundary between two subsystems. Because each site receives a fixed injection of magic from a finite gate density, the total magic therefore scales extensively with $L$ for $\beta=0$. For $\beta=1$, $\eta$ represents average number of $T$ gates per time step, independent of $L$. The resulting density is $p_{T}=\eta/L^{\beta}$; choosing $\beta=1$ drives $p_{T}\to 0$ in the thermodynamic limit while still providing a finite reservoir of magic.

When we add $T$-gates to the circuit in this model, the entanglement entropy does not change for each circuit, because the $T$-gates constitute a single layer of single-qubit unitaries at the end of the circuit. Any sequence of single-site unitaries acting on an initial state cannot alter the entanglement across any bipartition; the entanglement in the final state is identical to that of the initial state. This means that the entanglement entropy of the circuit with $T$-gates can be calculated by simulating the stabilizer circuit with the $T$-gates removed. It also means that the magic phase transition cannot be detected by the entanglement, since it is degenerate in the number of $T$-gates. We emphasize that this follows from the fact that the $T$-gates, CZ gates, and $Z$-basis measurements are all diagonal in the computational basis and therefore mutually commute. As a result, the effective action of the $T$-gates can be represented as a single layer applied at the final stage of the circuit.

The observed magic transition depends solely on the locations of measurements and $T$-gates, which are chosen stochastically. For each circuit realization, the amount of magic is determined by the measurement history, which dictates how many $T$-gates survive until the end of the circuit. The final state along a given trajectory is a stabilizer state dressed by a measurement-conditioned static layer of $T$-gates.

The entanglement dynamics in this model are governed entirely by the underlying monitored Clifford circuit. The $T$-gates do not modify the scrambling-measurement competition between the Clifford entangling operations and measurements. Conversely, the influence of $T$-gates on magic is completely independent of the Clifford unitaries. Magic depends exclusively on which $T$-gates survive.

Across many trajectories, the survival probability of $T$-gates depends on the measurement rate, giving rise to a phase transition. There is no dynamical interplay between entanglement growth and the propagation of non-Clifford resources. Instead, the distribution of surviving non-Clifford resources, and thus the scaling of the stabilizer nullity, remains a dynamical property of the monitored circuit ensemble. Although $T$-gates commute forward through the circuit, their survival probability depends on the stochastic and dynamical measurement history.

We locate the magic transition using two methods. First, we fit the data to the function

[TABLE]

with the constraint that $\gamma>0$ corresponds to (sub)-extensive scaling, and $\gamma\leq 0$ to area-law behavior. If the best fit yields $\gamma\leq 0$, we classify the data as area-law; otherwise, we interpret it as (sub)-extensive. Details of the fitting procedure are provided in Appendix D. The values of magic for the trajectories plotted in Fig. 6 empirically indicate a (sub)-extensive/log/area-law scaling of the magic with $L$. This is consistent with the scaling behavior of the magic, which was previously studied in bricklayer random Clifford circuits with a vanishing density of $T$-gates [14]. At the critical point $\beta=\beta_{c}$, we therefore invoke logarithmic scaling, i.e., $\mathcal{M}(L)=a+b\log L$ or

[TABLE]

For $L\rightarrow\infty$, $\mathcal{M}/\mathrm{log(L)}\rightarrow b$, so we expect a crossing at $b$ in the thermodynamic limit.

In Fig. 6(a), we identify four distinct phases in the magic and entanglement. The phases are area-law magic, volume-law entanglement (Phase I), area-law magic, area-law entanglement (Phase II), power-law magic, volume-law entanglement (Phase III), and area-law entanglement, power-law magic (Phase IV). In this model, the entanglement and purification transitions coincide. The exact exponent $\beta$ of the scaling plays a crucial role. Consider parameters $p_{m}$ and $\eta$ in regions III or IV. For some finite $\beta$, these points lie in the volume-law magic phase, while for $\beta=\infty$, no magic is present. We show numerically in Fig. 6(a), that for sufficiently low $p_{m}$, $\beta_{c}>1$ and decreases to $\beta_{c}=1$ for increasing $p_{m}$. This suggests that the $\beta=\infty$ “zero-magic” fixed point is stable. Conversely, for $p_{m}$ and $\eta$ in regions I or II, the system is in the area-law magic phase for some finite $\beta$, and is in the volume-law magic phase at $\beta=0$. For sufficiently high $p_{m}$ and low $\eta$, it is likely that a critical point $\beta_{c}<1$ exists.

We leave to future work as to whether, as $p_{m}\to 0$ and $\eta\to\infty$, the critical value $\beta_{c}$ also diverges. In other words, whether the $\beta=\infty$ “zero-magic” fixed point is unstable, implying that the system resides in the magic phase for any positive $\beta$. Similarly, it remains an open question whether, for sufficiently large $p_{m}$ and small $\eta$, $\beta$ tends to zero, placing the system in the area-law magic phase for any finite $\eta$.

VI Discussion and conclusion

In this work, we introduced a simulation method based on the LRSD that stores and evolves the full density matrix of a near-stabilizer state. This capability allows access to quantities nonlinear in the state, such as entanglement. Using this approach, we identified, in the limit of $Z$-basis-only measurements, a MIPT in magic that is distinct from the MIPT in entanglement.

To compute stabilizer nullity for large systems ($L=256$), we developed a Bell sampling method. From a resource-theoretic perspective, we analyzed when LRSD offers a computational advantage, with the number of entries $\mathcal{N}$ providing a practical bound on simulation hardness as the $T$-gate rate varies. Introducing $T$-gates alters entanglement depending on the initial state and circuit structure. In $Z$-basis measurements within our CZ-based circuits, they can be commuted to the end and leave entanglement unchanged. With $X$-basis measurements, they can actively modify entanglement.

For the $Z$-basis measurement model, we find four distinct phases in entanglement and magic. Phase III begins generating sub-extensive entanglement and magic—resources essential for universal quantum computation—while Phases II and IV remain low-entanglement and can be efficiently simulated with MPS. In the area-law magic regime (Phases I and II), the LRSD can further explore entanglement MIPTs in various circuit geometries and probe dynamical features such as operator growth and information scrambling via out-of-time-order correlators.

A limitation of our approach is that efficient simulation of the near-stabilizer $\rho$ for large ($L>32$) system sizes is possible either near the critical measurement rate, which separates the area and (sub)-extensive magic phases, or within the area law magic phase. This constrains the values of $\beta$ and $\eta$ in the $T$-gate rate $p_{T}=\eta/L^{\beta}.$

Previous work has used MPS to study the magic transition in 1+1 dimensions using Clifford unitaries doped with $T$-gates. However, MPS fails for the nonlocal gates in the all-to-all model. Our simulator allows us to characterize the magic transition in this geometry. Previous work has computed the entanglement transition in Haar-random all-to-all circuits [31] and the purification transition in Clifford all-to-all circuits [33]. Both of these works find similar scaling behavior, suggesting that the two universality classes are the same. Our results find critical exponents, $\nu=0.5$ and $z=1/5$, as these work for various rates of $T$-gates, and confirm that the entanglement and purification transitions coincide in this all-to-all model.

In summary, our work presents new methodologies for simulating quantum circuits near the stabilizer limit and characterizing distinct MIPTs. To study the magic transition, we implemented a scalable Bell sampling algorithm capable of computing stabilizer nullity for large system sizes. The LRSD allowed us to study the purification in these circuits, which is independent of the magic transition. In other geometries, the purification transition may be nontrivially affected by the rate of $T$-gates and be related to the magic transition. Future research directions include extending this framework to non-Clifford resources beyond $T$-gates, and exploring MIPTs in other circuit geometries, where the universality class may undergo nontrivial flows.

Data Availability

The data used in this work is available from the authors upon request.

Acknowledgements.

We thank David Huse, Grace Sommers, and Bryan Clark for useful discussions. This work was partially supported by the Army Research Office Grant No. W911NF-23-1-0144 (K.A.  and J.H.P.), the US-ONR grant No. N00014-23-1-2357 (H.P. and J.H.P.). K.A. and M.J.G. acknowledge support from the Defense Advanced Research Projects Agency (DARPA) under Agreement HR00112490357 and the NSF QLCI award OMA2120757. We are grateful to the OSG Consortium [Sfiligoi2009Pilot, PordesOpen2007, OSG, OSPool] for providing us computational resources to run simulations. This work was performed in part at the Aspen Center for Physics, which is supported by National Science Foundation grant PHY-2210452 (J.H.P.) and at the Kavli Institute for Theoretical Physics (KITP), which is supported by grant NSF PHY-2309135 (M.J.G., J.H.P.)

Appendix A Algorithms

In this section, we provide an algorithm to update the LRSD in Eq. (4) under Clifford operations. These are the Clifford unitaries and Pauli measurements. We also explain how to compute the Born probabilities of measurement outcomes, and update the state after the measurement. We also explain how to compute the Rényi entanglement entropy from the LRSD. In the final part, we develop algorithms to perform Bell Sampling on a specific class of non-Clifford states. We also present the method of computing the stabilizer nullity.

A.1 Clifford Evolution

In this section, we describe how to update the LRSD of a quantum state as it evolves under Clifford unitaries and Pauli measurements, see Eqs. (14) and  (15). Throughout, we write the density matrix in the form of Eq. (4) of the main text,

[TABLE]

where the coefficients $\lambda_{l}$ are real, the operators $\sigma_{l}$ are logical Pauli strings drawn from the set $\mathcal{L}_{\mathrm{LRSD}}$, and $\rho_{S}$ is a stabilizer density matrix. Clifford unitaries and Pauli measurements are resource‑free in the sense that they do not generate additional non‑Clifford resources, such as magic. Nevertheless, they alter both the logical operators and the stabilizer structure, so the LRSD must be updated after every such operation.

Let $\mathcal{S}=\expectationvalue{g_{1},\dots,g_{r}}$ denote the abelian stabilizer group of $\rho_{S}$, generated by the commuting Paulis $g_{i}$. For a mixed stabilizer state, $r<L.$ We then define a corresponding set of destabilizers $\mathcal{D}=\{d_{1},\dots,d_{r}\}$ with the properties $[d_{i},d_{j}]=0$ and $[g_{i},d_{j}]=0$ for $i\neq j$, while each pair $(g_{i},d_{i})$ anticommutes, $\{g_{i},d_{i}\}=0$. The set $\mathcal{D}$ complements the stabilizer generators but, unlike $\mathcal{S}$, does not form a group.

A.1.1 Clifford Unitaries

Equation (28) gives the update rule for a Clifford unitary $U_{C}$. Because any Clifford maps Pauli operators onto Pauli operators by conjugation, each logical operator transforms as $\sigma_{l}\mapsto U_{C}\sigma_{l}U^{\dagger}$, and the stabilizer state transforms as $\rho_{S}\mapsto U_{C}\rho_{S}U_{C}^{\dagger}$. Consequently the LRSD retains its form after the update: the set of stabilizer generators $\{g_{i}\}$, the destabilizers $\{d_{i}\}$, and the logical Paulis $\{\sigma_{l}\}$ are merely mapped onto other Pauli operators by the conjugation, while the coefficients $\lambda_{l}$ remain unchanged. The LRSD is then updated as:

[TABLE]

A.1.2 Pauli Measurements

We now discuss how to perform a projective measurement of the Pauli operator $P$, where

[TABLE]

The effect of the measurement depends on how $P$ commutes with elements of $\mathcal{L}_{\mathrm{LRSD}}$ and $\mathcal{D}.$ We use the notation $[P,\mathcal{S}]=0$ to mean that $P$ commutes with all elements of $\mathcal{S}$, and $[P,\mathcal{S}]\neq 0$ to mean that $P$ anticommutes with at least one element of $\mathcal{S}.$ We recall that $P$ either commutes or anticommutes with every other Pauli operator. The algorithm then branches into three cases. When $[P,\mathcal{L}_{\mathrm{LRSD}}]=0$, the stabilizer state is updated directly. When $[P,\mathcal{L}_{\mathrm{LRSD}}]\neq 0$ but $[P,\mathcal{S}]=0$, certain logical operators are dropped. When $[P,\mathcal{L}_{\mathrm{LRSD}}]\neq 0$ but $[P,\mathcal{S}]\neq 0$, a stabilizer decomposition is performed before applying the measurement.

By a stabilizer decomposition, we mean that

[TABLE]

where $\rho_{S}^{P}$ is a new stabilizer state such that $[\rho_{S}^{P},P]=0$ and $\bar{P}$ is a Pauli operator defined analogously to Eq. (29).

The proof of the evolution of $\rho$ in Case III is as follows. We give the proof for a projection onto $\frac{I+P}{2}$; the case for $\frac{I-P}{2}$ follows the same argument. First, using Algorithm 3, we can decompose $\rho_{S}=\frac{1+\bar{P}}{2}\rho_{S}^{P}$, where $[\rho_{S}^{P},P]=0$ and $\{\bar{P},P\}=0$. Then,

[TABLE]

There are two cases for each term labeled by $l$ in the sum. If $[P,\sigma_{l}]=0$, then we replace that term by

[TABLE]

Note that $\frac{I+P}{2}\rho_{S}^{P}$ defines a new stabilizer state $\rho_{S}.$ If $\{\sigma_{l},P\}=0$, then the term is

[TABLE]

Substituting Eqs. (32) and  (33) into each term of Eq. (31), we see that we arrive at the mapping in Case III. To evaluate whether $\{P,\mathcal{L}_{\mathrm{LRSD}}\}=0$ in Algorithm 5, and all other commutation relations, we use the symplectic inner product. A generic Pauli operator $P_{i}$ can be written as

[TABLE]

where $x_{ij},z_{ij}\in\{0,1\}$, where $i$ labels which Pauli operator it is. The symplectic inner product for two Pauli operators $P_{1}$ and $P_{2}$ is defined as:

[TABLE]

where $\oplus$ denotes addition modulo 2. This product equals 0 if $[P_{1},P_{2}]=0$ and 1 if $\{P_{1},P_{2}\}=0$. We now present an approach that leverages the stabilizer structure to efficiently compute the Born probability of projective measurements. The algorithm computes the probability of obtaining the outcome $+1$ when performing a projective measurement of operator $P$ on the general pure state $\ket{\psi}.$ The routine projects the LRSD‐represented state $\outerproduct{\psi}{\psi}$ onto the $+1$ eigenspace of a Pauli $P$. First, each logical Pauli $\sigma_{l}$ that commutes with $P$ is rewritten as a difference of projectors, allowing the LRSD sum to be expressed entirely in terms of projector products. The trace of this projected operator yields the probability $P_{\uparrow,P}$ of obtaining outcome $+1$. $P_{\uparrow,P}$ can be interpreted as a sum of terms, where each term is a conditional probability of obtaining two consecutive measurements with fixed outcomes. Summing over the terms reconstructs the total measurement probability.

We now present the Algorithm for performing a stabilizer decomposition as in Eq. (30). Firstly, $\rho_{S}^{P}$ is obtained from $\rho_{S}$ by removing a stabilizer generator $\overline{P}$ satisfying $\{\bar{P},P\}=0.$ To see how this decomposition is obtained, we have

[TABLE]

By a change of basis, we can choose at most one $g_{i}$ to anticommute with $P.$ The details of the procedure are given in Algorithm 3. Without loss of generality, suppose that $\{g_{1},P\}=0$ and $[g_{i},P]=0$ for all other $i.$ Then,

[TABLE]

Identifying $g_{1}=\overline{P}$ and

[TABLE]

we have Eq. (30).

We now explain how the state $\rho_{S}^{P}$ is stored in stabilizer simulations and is obtained from the representation of the original state $\rho_{S}.$ For a mixed stabilizer state $\rho_{S}$, the stabilizer formalism includes not only stabilizers and destabilizers, but also a set of logical operators denoted by $\mathcal{L}_{S}=\{l^{x}_{1},\dots,l^{x}_{L-r},l_{1}^{z},\cdots,l^{z}_{L-r}\}$. We group the logical operators into two blocks: the first $L-r$ operators denoted with $x$ and their anticommuting partners denoted with $z$. These logical operators satisfy the commutation relations $[l^{x}_{i},l^{x}_{j}]=[l^{z}_{i},l^{z}_{j}]=[l^{x}_{i},l^{z}_{j}]=0$ for $i\neq j$ and $\{l^{x}_{i},l^{z}_{i}\}=0$. Additionally, the logical operators commute with all stabilizers $s_{j}\in\mathcal{S}$ and destabilizers $d_{j}\in\mathcal{D}$, i.e., $[l^{x}_{i},s_{j}]=[l^{z}_{i},s_{j}]=[l^{x}_{i},d_{j}]=[l^{z}_{i},d_{j}]=0$ for all $i,j$.

A.2 Partial Trace and Rényi Entropy of LRSD

We now present algorithms to compute the Rényi entanglement entropies using the representation of the state in Eq. (4). We present an algorithm for computing the reduced density matrix for an arbitrary state, from which we can compute the entanglement. We also present an algorithm which is valid only for $[\mathcal{L}_{S},\mathcal{S}]=0.$ These methods allow us to compute the full entanglement spectrum for an arbitrary $n$, where:

[TABLE]

and $\rho_{A}$ is the reduced density matrix of some region $A.$

Our method to compute the entanglement relies on converting the LRSD representation of $\rho$ in the compact binary form into a $2^{|A|}\times 2^{|A|}$ matrix in the full Hilbert space. This restricts the size of $|A|$ for which we can compute the entanglement. However, non-trivial quantities such as correlation functions between two qubits and the ancilla entanglement entropies can be computed.

As a concrete example, we demonstrate how to compute correlation functions between an arbitrary pair of qubits using the LRSD representation. Given any two qubits $A$ and $B$, we trace out all remaining qubits $k\notin\{A,B\}$ (using Algorithm 4) to obtain the reduced density matrix $\rho_{A\cup B}=\mathrm{Tr}_{k}(\rho)$, which is a $4\times 4$ matrix. From this reduced density matrix, one can evaluate any two-qubit observable, including standard two-point correlators $\langle O_{A}O_{B}\rangle$. A particularly useful measure is the mutual information,

[TABLE]

which captures the total correlation—both quantum and classical—between any two arbitrary sites in the system. Here, $S_{n}$ denotes the Rényi entropy defined in Eq. 39. The mutual information provides an upper bound on all connected two-point correlation functions via $\langle O_{A}O_{B}\rangle-\langle O_{A}\rangle\langle O_{B}\rangle\leq\sqrt{2I_{1}(A\cup B)}$, and thus serves as a comprehensive diagnostic of two-site correlations accessible within our framework.

A.2.1 Partial Trace and Rényi Entanglement Algorithms

We now present a general-purpose algorithm to compute the Rényi entanglement entropies for arbitrary $n$. Given region $A$ with qubits $\{i_{1},\cdots,i_{k}\}$ and region $B$ with qubits $\{i_{1},\cdots,i_{L-k}\},$ the first step is to perform a partial trace of the state in Eq. (4) to obtain the reduced density matrix

[TABLE]

where $\rho_{S}^{A}$ is the reduced density matrix of $\rho_{S}$ and $\sigma_{\ell^{\prime}}^{A}$ are the logical operators restricted to region $A$. The size of the new logical set $\{\lambda_{\ell^{\prime}}\}$ satisfies $|\{\lambda_{\ell^{\prime}}\}|\leq|\{\lambda_{\ell}\}|.$ Each of the logical operators $\sigma_{\ell}$ and $\rho_{S}^{A}$ can then be converted to the full Hilbert space, represented as $2^{2k}\otimes 2^{2k}$ matrices, respectively. Using matrix operations, the matrix for the full state $\rho_{A}$ can be computed, from which Eq.  (39) can be easily evaluated.

The first step brings the tableau for stabilizer state $\rho_{S}$ into row-echelon form using the procedure of Ref. [3], which applies permutations and scalar multiplications to the generators. This form is particularly well-suited for partial trace operations and simplifies handling both pure and mixed stabilizer states.

We now present a method to perform a partial trace of $\rho$, which is represented using the LRSD in Eq. (4). This method extends the technique in Ref. [3], which shows how to compute the entanglement in stabilizer states. First, we represent the stabilizers of $\rho_{S}$ of Eq. (8) in an $r\times L$ matrix, where columns $1$ to $L$ represent qubits $1$ to $L$ and each row represents a different stabilizer. Each entry of the matrix is either the Identity or Pauli operator: $I$,$X$,$Y$,$Z$. Each logical operator $\sigma_{l}$ is represented analogously. The first step is to bring the tableau of $\rho$ into a reduced row echelon form (RREF), following the procedure in [3]. The general structure of the RREF can be described in a recursive fashion. There are three cases:

[TABLE]

The symbols “$\vdots$” and “$\cdots$” denote a number of repeated rows and columns. This number can be zero. The symbol RREF denotes a possibly empty sub-array that is also in RREF form. The symbol $*$ denotes either a Pauli operator or $I$. Furthermore, $\sigma$, $\sigma_{1}$, and $\sigma_{2}$ are Pauli operators and $\sigma_{1}$ and $\sigma_{2}$ anticommute. We denote the sequences of $*$ operators by $g$, $g_{1}$ and $g_{2}$, respectively. The first column in each of the three cases represents a qubit to be traced over. The stabilizer in RREF’ defines a new stabilizer state $\rho^{\prime}$ which is the reduced density matrix obtained after tracing over the first qubit. In each of the three cases, we can decompose $\rho$ into its action on the qubit to be traced over and the rest of the system. In Case I, recalling that we can substitute the generators $g_{i}$ from the tableau into Eq. (1), we see that:

[TABLE]

In Case II, we can also use Eq. (1) to observe that

[TABLE]

Finally, in Case III, we have that

[TABLE]

We must have that $g_{1}\neq g_{2}$, since $g_{1}=g_{2}$ requires $\sigma_{1}=\sigma_{2}$ to preserve the property that all stabilizers commute. However, $g_{1}=g_{2}$ and $\sigma_{1}=\sigma_{2}$ leads to duplicate stabilizer generators. For each logical $\sigma_{l}$, we can simply decompose it into

[TABLE]

where $\sigma_{r}$ is the action of $\sigma_{l}$ on the rest of the system. Now, we can present our algorithm to calculate the partial trace of the LRSD.

Step 1. Using column permutations, bring the columns of the qubits to be traced out to the first position. Do this for each of the logicals $\sigma_{l}$ and the tableau of $\rho_{S}$.

Step 2. Bring the stabilizer tableau into RREF*′*.

Step 3. For each qubit to be traced over, apply Algorithm 4.

In Algorithm 5, we present an algorithm for updating the logical set $\mathcal{L}_{S}$ and the destabilizer set $\mathcal{D}$ when performing a measurement of the Pauli operator $P$ in Eq. (29). We now describe how to perform a measurement of the Pauli operator $P$, Eq. (29), with the projectors $(I\pm P)/2$, on the stabilizer state $\rho_{S}.$ After measurement, $\pm P$ becomes a stabilizer of $\rho_{S}$, and must therefore be to the stabilizer group. The list of stabilizer generators must then be updated so that they all commute with $P$. Furthermore, the list of destabilizers must be updated to preserve the required commutation relations.

If a single stabilizer $\bar{g}_{i_{1}}$ anticommutes with $P$, we replace $\bar{g}_{i_{1}}$ with $\pm P$ and $\bar{d}_{i_{1}}$ with $\bar{g}_{i_{1}}$. When multiple stabilizers anticommute with $P$, we pair each anticommuting stabilizer and destabilizer into products and perform the updates $\bar{g}_{i_{1}}\mapsto\bar{g}_{i_{1}}\cdot\bar{g}_{i_{k}}$, so the stabilizers now commute with $g_{i_{1}}$. The formation of a list of these products is equivalent to constructing a tree graph, shown in Fig. 7, where the nodes represent the anticommuting stabilizers and the edges represent the choices of pairs to multiply.

We now present an algorithm for computing the Rényi entanglement entropy, which is valid only for the case $[\mathcal{L}_{\mathrm{LRSD}},\mathcal{S}]=0.$ In this case, we can write the density matrix $\rho$ as a sum over pairs of stabilizer states, where the states in each pair are orthogonal to each other. A partial trace over the environment is carried out for each stabilizer term. After tracing out the environment, we obtain a mixture of reduced density matrices. The purity, $\mathrm{tr}(\rho_{A}^{2})$, can be computed from overlaps of these density matrices using standard stabilizer techniques, from which $S_{n=2}$ can be easily computed. To compute the entanglement for arbitrary $n$, the reduced density matrix can be obtained from a partial trace over each term.

A.3 Bell Sampling and Magic Computation Algorithms

A.3.1 Bell sampling

To compute Eq. (13), we need to Bell sample the Pauli strings from the state $\ket{\psi}$, as illustrated in Fig. 8. To define the Bell samples, we first define the Pauli matrices $\sigma_{00}=I$, $\sigma_{01}=X$, $\sigma_{10}=Z$, and $\sigma_{11}=Y$. The $4^{L}$ Pauli strings are $L$-qubit tensor products of Pauli matrices, which can be written as $\sigma_{\bm{r}}=\bigotimes_{j=1}^{L}\sigma_{r_{j},r_{j+L}}$, where $\bm{r}\in\{0,1\}^{2L}$. We first prepare two copies of the state, $\ket{\Psi}=\ket{\psi}\otimes\ket{\psi}$, and then measure pairs of qubits $(j,j+L)$ in the Bell basis, which is defined as,

[TABLE]

Following Ref. [23], the probability of obtaining the sequence of outcomes $\vec{r}$ is:

[TABLE]

where $\ket{\psi^{*}}$ is the complex conjugate of $\ket{\psi}$  [30]. For a stabilizer state $\ket{\psi}$ with stabilizer group $\mathcal{S}\subset\mathcal{P}_{L}$, there exists a unique Pauli operator $g$ such that $g\ket{\psi}=\ket{\psi^{*}}.$ $p(\vec{r})$ is non-zero if and only if $\sigma_{\vec{r}}$ belongs to a single coset $g\mathcal{S}$ of the stabilizer group $\mathcal{S}$ of $\ket{\psi}.$ Bell sampling produces a uniform distribution over all $2^{L}$ elements of the coset. All operators inside that coset are sampled with equal probability $1/2^{L}$, while those outside have zero probability. The Bell samples lie in a vector space of dimension $G^{\prime}=L$, so $\mathcal{M}=0$ for stabilizer states. For a non-stabilizer pure state, the Bell samples are not in a single coset. The Bell samples form an unstructured subset of $\mathcal{P}_{L}$, which has dimension $G^{\prime}>L$, so $\mathcal{M^{\prime}}>0$ for non-stabilizer states.

We develop a Monte Carlo method for Bell sampling from the distribution $p(\vec{r}),$ in Appendix A.3, building on methods introduced in Ref. [23]. Bell sampling can be intractable for large $L$ as the number of $\sigma_{\vec{r}}$ for which $p(\vec{r})=0$ is exponentially growing in $L$. Our method always generates samples $\sigma_{\vec{r}}$ with $p({\vec{r}})\neq 0$ at each sampling step. A sampling step means drawing one independent sample from $p(\vec{r})$. We compute $\mathcal{M}$ after each sampling step, and find that $\mathcal{M}$ converges in fewer than $\tau<2L$ sampling steps for all measurement rates, provided $L\leq 256$ and $0<\beta<1$, at $T$-gate rates given by $p_{T}=12/L^{\beta}$, as shown in Appendix A.3. We also find that our Monte Carlo sampling method converges to the exact value of the magic obtained by exhausting all possible $4^{L}$ samples in polynomial time. A reason we use the stabilizer nullity to quantify magic is that it admits an exactly correct result in a very low number of samples, shown in A.3.

A.3.2 Bell Sampling Algorithm

We now introduce a Monte Carlo sampling algorithm that efficiently produces Bell samples for a specific subclass of states that include, but are not limited to, stabilizer states. The sampler only visits Bell samples, $\sigma_{\vec{r}}$, which have nonzero probabilities $p(\vec{r})>0$, where $p(\vec{r})$ is defined in Eq.  (58). This algorithm is adapted from Ref. [23]. The non-stabilizer state must be able to be written as $\ket{\psi}=\mathcal{T}\mathcal{C}\ket{\psi_{o}}$, where $\ket{\psi_{o}}$ is a stabilizer state, $\mathcal{C}$ is a Clifford circuit, and $\mathcal{T}$ is a layer of $T$-gates. The non-Clifford layer factorizes as

[TABLE]

where $N_{i}$ is the number of $T$-gates applied to qubit $i.$

Using $T^{2}=S$, $T^{4}=Z$, and $T^{8}=I$, each power $T^{N_{i}}$ can be decomposed into a product of Clifford and non-Clifford gates:

[TABLE]

Here, $\left\lfloor\cdots\right\rfloor$ is the floor function and $S=\mathrm{diag}(1,i).$ The explicit Clifford expansion allows at most one residual non-Clifford gate to act on each qubit. The rest of the $T$-gates are absorbed into the Clifford part of the circuit, so the number of $T$-gates remaining is $k=\sum_{i=1}^{n}a_{i}.$

We now consider how the density matrix $\rho_{S}$, for the initial state $\mathcal{C}\ket{\psi_{o}}$, evolves under the $T$-gate layer in Eq. (59). The final state $\rho$ can be represented as a product of a mixed stabilizer state and a product over $k$ terms, each of which is written as a sum over three projectors:

[TABLE]

where $\alpha_{1}=\frac{1}{2}$, $\alpha_{2}=\frac{1-\sqrt{2}}{2}$, $\alpha_{3}=\frac{\sqrt{2}}{2}$. Expanding the product over the $k$ stabilizer triplets yields a weighted sum of $3^{k}$ stabilizer states. The $L-k$ operators $g_{1},\cdots,g_{L-k}$ stabilize all of the stabilizer states. We call these isotropic stabilizers using the terminology in Ref. [23]. Along with a set of isotropic destabilizers $d_{1},\cdots,d_{k}$, they define a mixed stabilizer state. There are additional $k$ triplets of stabilizer generators $(g_{L-k+1},-g_{L-k+1},h_{L-k+1}),\cdots,(g_{L},-g_{L},h_{L}).$ Each of these generators is called a primary symplectic stabilizer (PSS) following Ref. [23]. Each PSS stabilizes $1/3$ of the stabilizer states. We require the commutation relations $[g_{i},g_{j}]=0,[h_{i},h_{j}]=0,\{g_{i},h_{j}\}=0$ for $i=j$ and $[g_{i},h_{j}]=0$ for $i\neq j$. The state is written in this form, as Bell sampling can be performed efficiently. The Bell samples can be quickly produced from products of the isotropic stabilizers and PSS following Algorithm 10.

A.3.3 Magic Computation

In this section, we explain how to compute the magic of a circuit using the representation in Eq. (61). The key idea is to use the stabilizer tableau, a compact binary encoding of the stabilizers, to bring the stabilizer generators into a compact form as in Eq. (62). The procedure to compute the magic is summarized in Algorithm 7. We recall the explicit form of the stabilizer tableau for $\rho_{S}.$ We first recall again that a Pauli operator $P$ can be expressed as a $2L$ binary vector $(\mathbf{x}|\mathbf{z})$ as discussed at the beginning of Section IV. The state $\rho_{S}$ has generators $\{g_{1},\cdots,g_{L}\},$ where $g_{i}=c_{i}P_{i}$ and $c_{i}=\pm 1.$ Then, the tableau is

[TABLE]

where $\mathrm{sgn}(c_{i})=\frac{1-c_{i}}{2}.$ This can be more compactly represented as: $[X|Z|\vec{c}]$, where $X,Z\in\{0,1\}^{L\times L}$ and $\vec{c}\in\{0,1\}^{L}$.

After generating samples $\sigma_{1},\cdots,\sigma_{\mathcal{N}_{s}}$ from the Bell measurement procedure in Step 3 of Algorithm 7, we can proceed to estimate the stabilizer nullity $\mathcal{M}$ of the state. This estimation is performed by computing the rank of pairwise differences between the strings.

We now present an algorithm for Bell sampling from the state in Eq. (61). Recall from Section A.3.1 that a Bell sample $\sigma_{\vec{r}}$ is obtained from $\ket{\psi}\otimes\ket{\psi}$ with probability $p(\vec{r})$. This is equivalent to measuring each qubit pair in the Bell basis and obtaining the outcome sequence $\vec{r}$ with Born probability $p(\vec{r})$. We now define $\sigma_{\vec{r}_{o}}$ such that

[TABLE]

and

[TABLE]

We recall that a Bell sample $\sigma_{\vec{r}}$ from $\ket{\psi}\otimes\ket{\psi^{*}}$ occurs with probability [30]

[TABLE]

From Equations (66) and  (67) we see that obtaining a Bell sample $\sigma_{\vec{r}}$ from $\ket{\psi}\otimes\ket{\psi}$ is equivalent to obtaining a Bell sample $\sigma_{\vec{r^{\prime}}}$ from $\ket{\psi}\otimes\ket{\psi^{*}}$, and then shifting it by $\vec{r}_{o}$: $\sigma_{\vec{r}}=\sigma_{\vec{r^{\prime}}\oplus\vec{r_{o}}}.$ We then present in Algorithm 9 how to compute $\sigma_{\vec{r}_{o}}$ and then in Algorithm 10 how to obtain Bell samples from $\ket{\psi}\otimes\ket{\psi^{*}}.$

We recall from Ref. [23] that, if $\pm\sigma_{r}$ is an isotropic stabilizer, then

[TABLE]

where $L$ is the number of qubits.

If $\pm\sigma_{r}$ is a product of $m$ Pauli-string stabilizers (an $m$-PSS), the probability of sampling any $m$-PSS is

[TABLE]

where $n_{m\text{-PSS}}$ denotes the total number of such objects.

There are $k$ disjoint PSS pairs. To form an $m$-PSS, we (i) choose the $m$ pairs that actually appear (giving the factor $\binom{k}{m}$), (ii) choose one stabilizer from each selected pair ($g_{n-k+i}$ or $h_{n-k+i}$, giving $2^{m}$ choices), and (iii) multiply by the $2^{n-k}$ isotropic stabilizers that commute with all pairs. Thus

[TABLE]

Substituting (70) into (69) and using (68) gives

[TABLE]

To test the accuracy of Algorithm 10, we fix the location of the unitary gates and measurement locations for many circuits. We sample all $4^{L}$ possible Pauli operators on $L$ qubits. We then compute $\mathcal{M}_{\mathrm{exact}}$ using Eq. (64) with these samples. We also obtain a numerical estimate $\mathcal{M}_{\mathrm{numerical}}$ for the same set of circuits, using Algorithm 10 with $\tau$ samples. Then, define $|\delta\mathcal{M}|=\mathcal{M}_{\mathrm{numerical}}-\mathcal{M}_{\mathrm{exact}}.$ The absolute value $|\delta\mathcal{M}|$ will then be positive and decrease to 0 with increasing $\tau$, as shown in Fig. 9.

We now demonstrate our algorithm on the single-pair all-to-all model with $Z$-basis measurements. From the discussion at the beginning of Section IV.2, we observe that the state can be factorized with the $T$-gate layer in Eq.  (59) acting at the final timestep.

Fig. 10 illustrates the convergence behavior of the stabilizer nullity $\mathcal{M}$ and the number of distinct Bell samples $\mathcal{N}_{s}$ as a function of the number of samples $\tau$, for different measurement rates $p_{m}.$ In the top row, we observe that $\mathcal{M}$ exhibits rapid growth followed by saturation for all values of $p_{m}.$ The saturation value decreases as $p_{m}$ increases, reflecting expected suppression of non-stabilizerness by projective measurements. Additionally, saturation occurs earlier for larger $p_{m}$, which correlates with the reduced complexity of the state. The bottom row, $\mathcal{N}_{s}$, the number of distinct Bell samples obtained during sampling, is independent of the measurement rate $p_{m}.$ It grows at a rate of approximately one per sampling step, indicating that it is rare for the same Pauli to be sampled twice in the same Monte Carlo run.

Appendix B Purification, cluster, and all-to-all entanglement transitions

B.1 Cluster

We begin with the connectivity transition, which is a percolation transition from a connected state (when all qubits are entangled) to a disconnected state (when the state effectively breaks into multiple disentangled clusters). The quantum state factorizes as

[TABLE]

where each $C_{i}$ denotes a disjoint subset of qubits, or cluster, and $\ket{\psi_{C_{i}}}$ is an entangled pure state supported on that cluster. The full state $\ket{\psi}$ is a product across these clusters with no entanglement between disjoint $C_{i}$.

Therefore, we can probe the critical connectivity (cc) measurement rate $p_{m}^{cc}$ by tracking the trajectory-averaged number of qubits in the largest cluster $n_{\text{max}}\equiv\absolutevalue{C_{\text{max}}}$. Inside the connected phase, the size should be $\mathcal{O}(L)$, as all qubits are entangled, whereas in the disconnected phase, the number should be $\mathcal{O}(1)$. Numerically, we convert each stabilizer code into an equivalent graph state to identify entangled clusters efficiently [43, 2]. This graph-state conversion is only valid for stabilizer states; when non-Clifford gates such as $T$-gates are present, it is not yet known how an analogous mapping can be performed. Because every stabilizer state is local‑Clifford equivalent to a graph state, we first apply single‑qubit Clifford rotations that bring the stabilizer generators $g_{i}$ into the canonical graph form:

[TABLE]

where $\mathcal{N}(i)$ is the neighbor set of qubit $i$. We then define a graph $G=(V,E)$ with one vertex per qubit, $V=\{1,\dots,L\}$, and an edge $(i,j)\in E$ whenever $j\in\mathcal{N}(i)$; equivalently, each edge corresponds to a $\mathrm{CZ}$ gate in the standard preparation circuit:

[TABLE]

The connected components of $G$ give precisely the entangled clusters $C_{i}$.

We study the connectivity transition in the single-pair all-to-all model with both $X$-basis and $Z$-basis measurements. The probability of performing a measurement in the $X$-basis given a measurement is performed is denoted by $p_{x\vee z}$. We also turn off the $T$-gates in the circuit, $p_{T}=0$. The numerical phase diagram is shown in Fig. 12, where we find that the estimated phase boundary between disconnected phase and connected phase is roughly always around $p_{m}^{cc}\approx 0.66(1)$, independent of $x$-basis measurement $p_{m,x}$.

This phase boundary of the connectivity transition is detected by directly estimating the number of qubits in the largest cluster defined as [17]

[TABLE]

where $C_{i}$ represents the $i$-th cluster and $|C_{i}|$ is the number of qubits in that cluster. As shown in Fig. 11(b), we find that the trajectory-averaged number of qubits in the largest cluster $\overline{n_{\text{max}}}$ grows extensively in the connected phase while it remains $\mathcal{O}(1)$ in the disconnected phase.

This motivates us to define the following finite size scaling ansatz for the $\overline{n_{\text{max}}}$

[TABLE]

where the estimated critical exponents are shown in Table 2.

A close variant of this model is introduced in [45], which exhibits an exactly solvable phase transition between a fully entangled volume-law phase and a separable phase composed of finite subsystems. Because these specific unitary dynamics and projective measurements can be recast as a measurement-outcome-dependent unitary transformation, the model allows for an exact analytical treatment that maps the separability transition to mean-field percolation. Furthermore, this model serves as a compelling starting point since this entanglement transition coincides with a sharp phase transition in the computational hardness of classical simulation.

B.2 Entanglement Transition

To study the entanglement transition, we need to generalize the concept of entanglement entropy and the corresponding volume/area-law phase to the all-to-all model, as there is no clear notion of surface and bulk in this effective 0D system. Therefore, we generalize the entanglement entropy to the following max-min form [17]

[TABLE]

where $A$ is a half of a bipartition of cluster $C_{i}$, and $S_{n}$ is the $n$-th Rényi entropy. This max-min procedure is designed to find the most entangled cluster and then the least entangled bipartition within it, providing a robust measure that avoids being skewed by small, disconnected clusters. Definitions of entanglement which are based on averaging over bipartitions tend to yield $p_{ce}=p_{cc}$ on all-to-all circuits. To distinguish the entanglement transition from this connectivity transition, we introduce a minimization procedure over bipartitions. The generalized entanglement entropy should be extensive in the volume-law phase, and sub-extensive in the area-law phase.

In Fig. 13(b), we show the max-min entanglement entropy as a function of system size. At small measurement rates, $p_{m}<p_{m}^{ce}$, $S_{mm,1}$ scales extensively in system size. At large measurement rates, $p_{m}>p_{m}^{ce}$, $S_{mm,1}$ scales sub-extensively in system size. The critical measurement rate $p_{m}^{ce}$ can thus be identified by locating the inflection point of $S_{mm,1}$ on a linear-log plot, where $S_{mm,1}$ grows logarithmically in system size, i.e.,

[TABLE]

We find $p_{m}^{ce}=0.27(1)$.

B.3 Purification Transition

This section provides additional data supporting the purification transition results reported in Sec. V. We give more details on how we locate the purification transition in Section V. We present the data used to numerically extract the phase boundary of the purification transition in Fig. 5(a). We also explain the fitting procedure used to extract the decay times $\tau$ shown in Figs. 5(b) and (c). We obtain the decay times $\tau$ by fitting the late-time exponential regime of $\overline{S_{Q}}(t)$.

B.3.1 Fitting Procedure

We do not fit all of the data, $\overline{S}_{Q}(t)$ from $t_{i}=0$ to $t_{f}=2L^{2}$ to the exponential decay ansatz $\overline{S}_{Q}(t)=Ae^{-t/\tau}.$ Instead, we fit over all possibilities of $t_{f}$ and $t_{i}$, where $t_{f}-t_{i}\geq L^{2}/8$, and $\frac{\text{err}[S_{Q}(t_{f})]}{\overline{S_{Q}(t_{f})}}<0.35.$ The first condition ensures that an artificially small region is not being fit. The second condition ensures that the error of each datapoint is sufficiently small relative to the datapoint value. We compute $\tau$ using the region that gives the largest coefficient of determination, $R^{2}$, for the fit.

Figures 15 and  16 show the purification transition in the single-pair all-to-all model with measurements in the $X$-basis only. In this model, where the entanglement is no longer degenerate in the rate of $T$-gates, the purification transition remains unaffected.

B.3.2 Additional Result on Purification Transition

We now present data supporting the phase boundary of the purification transition (Fig. 5) in the single-pair all-to-all model with $X$-basis only measurements. Figures 15 and 16 show the scaling behavior of the decay time $\tau$ as a function of both system size $L$ and measurement probability $p_{m}$ for different $T$-gate densities. These results provide further confirmation that the purification transition persists even in the presence of non-Clifford gates, and allow for an independent extraction of the dynamical exponent $z$ and the critical purification rate $p_{m}^{cp}$.

Appendix C Statistical Properties of the LRSD Space Complexity

In this appendix, we provide an extended statistical analysis of the space complexity $\mathcal{N}$ to supplement the ensemble-averaged scaling results presented in Section  IV.3. Here, we present the full cost distributions to confirm the absence of exponentially hard trajectories, analyze the second moment of these distributions to rigorously bound the variance, and provide additional scaling data utilizing random Clifford gates to demonstrate the robustness of this efficient scaling regime across different circuit structures. Furthermore, to characterize the scaling complexity and demonstrate robustness to generic circuit structures, we completely replace the CZ gates with random Clifford gates in Fig. 18.

C.1 Distribution of Simulation Costs

Since the number of logical operators $|\{\lambda_{\ell}\}|$ directly dictates the memory and computational overhead of the LRSD algorithm [Eq. (19)], evaluating its trajectory-resolved probability distribution, $p(|\{\lambda_{\ell}\}|)$, serves as a direct proxy for the simulation cost.

In Fig. 17, we plot $p(|\{{\lambda_{\ell}}\}|)$ for varying system sizes, $L\in\{16,48,96\}$, across a range of sub-extensive $T$-gate rates from $p_{T}=0.25/L$ to $p_{T}=1/L$. The empirical distributions demonstrate a distinct unimodal structure. As expected, as both the $T$-gate injection rate and system size increase, the mean of the distribution shifts to higher values and the variance broadens.

The probability weight on the right-hand side of the mean exhibits a power-law decay $p(|{\{\lambda_{\ell}\}|)}\sim|{\{\lambda_{\ell}\}}|^{-\alpha}.$ This is observed in for $p_{T}=0.25/L$ to $p_{T}=1/L$ at $p_{m,X}=0.6$. In contrast to regimes with exponential suppression, this behavior suggests that the likelihood of encountering highly demanding trajectories decreases more gradually. As a result, even at low injection rates, such extreme events can retain non-negligible probability weight. This observation is consistent with the presence of “hard trajectories” in this regime, indicating that these rare but costly realizations may play an important role in shaping the average-case sampling dynamics. Importantly, we emphasize that this feature does not restrict our claims regarding efficient simulability within the parameter regimes studied here as we reach $L=96$ despite this heavy-tailed behavior.

Appendix D Magic Analysis

We fit the data, $\mathcal{M}(L)$ to the function $f(L)=a+bL^{\gamma}.$ For the power-law fit, we impose the constraint $\gamma>0.$ For the area-law fit, we impose the constraint $\gamma\leq 0.$ If the best fit is $\gamma\leq 0$, then the scaling is classified as area-law. Otherwise, it is classified as power-law scaling.

To assess the quality of each fit, we compute the chi-squared fit parameter using:

[TABLE]

Here, $\overline{\mathcal{M}}_{\mathrm{ss}}(L)$ denotes the steady-state stabilizer nullity, and $\sigma(L)$ is the standard error of each datapoint. The degrees of freedom (dof) are the number of data points minus the number of fitted parameters.

To identify a proxy for the critical $T$-gate rate $\beta_{c}$, we plot $\chi^{2}/\mathrm{dof}$ as a function of $\beta$ for three different scaling ansatze: the (sub)-extensive and area-law power-law forms, and the logarithmic scaling function $f_{\mathrm{log}}(L)=a+b\ln L$. In Fig. 19, we plot $\chi^{2}/\mathrm{dof}$ against the scaling parameter $\beta$ for the (sub)-extensive and area law scaling fit functions. In addition, we also plot $\chi^{2}/\mathrm{dof}$ for fits of the logarithmic scaling function.

We define $E_{\mathrm{ext}}(\beta)=\chi_{\mathrm{ext}}^{2}/\mathrm{dof}$ and $E_{\mathrm{area}}(\beta)=\chi_{\mathrm{area}}^{2}/\mathrm{dof}.$ $\chi_{\mathrm{ext}}^{2}$ and $\chi_{\mathrm{area}}^{2}$ are the extensive and area law fits, respectively. We then define $\beta_{c}$ as the value at which $\mathrm{ln}[F(\beta_{c})]=0.$ To identify a transition in the measurement rate $p_{m}$ at a fixed $\beta$, we follow the same procedure.

We estimate the uncertainty $\sigma_{\beta_{c}}$ as

[TABLE]

with two contributions: To compute $\sigma_{A}$, we find the values $\beta^{\pm}$ where $\mathrm{ln}[F(\beta^{\pm})]=\pm 1$, and then take $\sigma_{A}=(\beta^{+}-\beta^{-})/2.$ The second contribution is $\sigma_{B}=\left|\beta_{c}(L_{\max}=256)-\beta_{c}(L_{\max}=192)\right|$. Here, $\beta_{c}(L_{\max}=192)$ denotes the estimate of $\beta_{c}$ computed by excluding data for the largest system size $L=256$.