key: cord-0702993-xlhqm6lt
authors: Sierant, Piotr; Turkeshi, Xhek
title: Universal behavior beyond multifractality of wave-functions at measurement--induced phase transitions
date: 2021-09-14
journal: Phys Rev Lett
DOI: 10.1103/physrevlett.128.130605
sha: 2ee9d2286af1ae6c29e2654a6f3578dbf471d931
doc_id: 702993
cord_uid: xlhqm6lt

We investigate the structure of many-body wave functions of 1D quantum circuits with local measurements employing the participation entropies. The leading term in system size dependence of participation entropy indicates a model dependent multifractal scaling of the wave-functions at any non-zero measurement rate. The sub-leading term contains universal information about measurement-induced phase transitions and plays the role of an order parameter, being constant non-zero in the error correcting phase and vanishing in the quantum Zeno phase. We provide robust numerical evidence investigating a variety of quantum many-body systems, and provide an analytical interpretation of this behavior expressing the participation entropy in terms of partition functions of classical statistical models in 2D.

We investigate the structure of many-body wave functions of 1D quantum circuits with local measurements employing the participation entropies. The leading term in system size dependence of participation entropy indicates a model dependent multifractal scaling of the wave-functions at any non-zero measurement rate. The sub-leading term contains universal information about measurement-induced phase transitions and plays the role of an order parameter, being constant non-zero in the error correcting phase and vanishing in the quantum Zeno phase. We provide robust numerical evidence investigating a variety of quantum many-body systems, and provide an analytical interpretation of this behavior expressing the participation entropy in terms of partition functions of classical statistical models in 2D.

The evolution of quantum systems driven by the interplay of unitary dynamics and the monitoring action of an environment follows a quantum trajectory through the Hilbert space [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] . While the former generates coherences and spreads information throughout the system, the latter partially resolves the state, which is collapsed according to the Born-Von Neumann postulate. In the many-body framework, this competition gives rise to dynamical phases captured by non-linear functions of the density matrix separated by so-called measurementinduced phase transitions (MIPTs) .

Random quantum circuits, believed to model generic chaotic quantum dynamics [61] [62] [63] [64] [65] [66] [67] , become minimal models to investigate MIPTs when their discrete spacetime is punctured by local projective measurements [68, 69] . Entanglement provides natural measures to explore the phase diagram of these models [70] . Their scaling behavior with system size captures key facets of the dy- We analyze this change using the PEs (1) . The data is for a Haar circuit of L = 12 qubits.

namical phases and justifies their intriguing interpretation in terms of quantum error correction [71] [72] [73] [74] [75] [76] . Specifically, when the measurement rate is tuned, the system undergoes a transition from an error-correcting phase to a quantum Zeno phase, characterized respectively by an extensive ("volume law") and a sub-extensive ("area law") scaling of the entanglement entropy.

The local measurements partially project the state |ψ m of the circuit affecting the structure of many-body wave-function ψ m ( σ) = σ|ψ m , at the quantum trajectory m, in the basis {| σ } (Fig. 1) . The structural changes are probed via the participation entropies (PEs)

where q > 0 and the second equality parametrizes the scaling of S q with system size L by a fractal dimension D q and a sub-leading term c q . The PEs are deeply tied to the concepts of inverse participation ratio (for q = 2) and (multi)fractal dimension [77] . They play a significant role in the Anderson [78] [79] [80] [81] and the many-body localization transitions [82] [83] [84] [85] [86] [87] [88] . Moreover, the PEs distinguish phases of quantum matter [89] [90] [91] [92] [93] [94] [95] [96] [97] [98] without relying on the system specific observables.

In this work, we compute PEs across MIPTs in the stationary state of random stabilizer and Haar circuits. For stabilizer states, we show the absence of multifractal behavior, signaled by S q being q−independent, whereas for Haar circuits, the system state is multifractal for any finite measurement rate, with D q < 1 depending in a nontrivial fashion on q. In analogy to quantum phase transitions [99, 100] , we find that the sub-leading term c q contains universal information about MIPTs. We provide an analytical understanding of those results by mapping the calculation of PEs onto a classical partition function [101] [102] [103] [104] [105] [106] and show that our conclusions apply to a large class of quantum circuits.

We consider random quantum circuits with local projective measurements acting on L qubits in the geometry shown in Fig. 1 (periodic boundary conditions are assumed). The two-qubit gates U j,j+1 are sampled uniformly from the Clifford group for stabilizer circuits and drawn randomly from the Haar distribution on SU (4) for the Haar circuits. Each timestep consists of a layer of projective measurements onto the z component of the spins, performed on each site j = 1, . . . , L with probability p, and of a layer of unitary gates. The system is in an error correcting phase for p < p c and in a quantum Zeno phase for p > p c , with the critical rate identifying the MIPT (p S c = 0.1593(5) for stabilizer circuits [69, 74] , while p H c = 0.17 (1) for Haar circuits [18] ).

We initialize the circuit in the product state |ψ 0 = |0, . . . , 0 and calculate its time evolution for t layers [107] . Stabilizer circuits are simulated in time polynomial in L (up to L ≤ 880) using the package Stim [108] and employing the ideas in [109] [110] [111] . For Haar circuits, we perform exact simulations in the full Hilbert space, up to L ≤ 24 qubits. In both cases, we compute the PEs for each quantum trajectory S q (|ψ m ) specified by the realization m, and consider the mean value. For the stabilizer circuits PEs, the average is over the times 2L ≤ t ≤ 22L and over N ≥ 20000 realizations, while, in the Haar case, we average over 2L ≤ t ≤ 1000, and over N ≥ 4000 (1000) circuit realizations for L < 24 (L = 24).

Participation entropy of stabilizer states. For a stabilizer state |ψ on L qubits there exists L independent Pauli strings g j that stabilize |ψ [112] : g j |ψ = |ψ . The g j generate a group G and can be written as

where X k , Z k are Pauli matrices acting on k-th qubit and n j k , m j k , φ j are equal to 0 or 1. Hence, the stabilizer state |ψ is uniquely determined by the matrices M X = [n j k ], M Z = [m j k ] and the vector of phases Φ = [φ j ]. Its density matrix reads ρ = |ψ ψ| = 2 −L g∈G g and the PEs (1) are

For specificity, we choose the eigenbasis of Z k operators (the Z basis) as the basis {| σ }. Then, the matrix element σ|g| σ is non-zero only if the Pauli string g contains no X i operators. Such Pauli strings form a subgroup G ⊂ G, generated by g i = L j=1 g a j i j , with a i being the solutions of M X a = 0 over the field Z 2 .

The independence of g j implies that there are r = L − rk Z2 M X independent generators of G (where rk Z2 M X is the rank of the matrix M X over Z 2 ), hence |G | = 2 r . From (2) we obtain that The operator P G ≡ 2 −r g ∈G g is a projector onto a subspace stabilized by the elements of G . By definition, the elements of G are diagonal in the the basis {| σ }, hence σ| P G | σ can be only 0 or 1, allowing to change the order of sums in (3) . Since the Pauli matrices are traceless, σ σ|g | σ is non-vanishing and equal to 2 L only if g is the identity operator. We conclude that

Crucially, the result (4) holds for an arbitrary stabilizer state and is independent of q, implying that wave-functions of stabilizer states do not exhibit multifractality. This holds true for generic eigenbases of Pauli strings [113] and reminisces about the fact that for stabilizer states all Rényi entanglement entropies are equal [114, 115] . The formula (4) has a simple interpretation. If M X = 0, the stabilizer state |ψ is fully localized in the Z basis and the PE is vanishing. Each linearly independent row of M X corresponds to a generator g j that delocalizes |ψ in the Z basis, incrementing PE by one.

Stabilizer circuits. We now turn to investigation of the MIPT in the Clifford circuits and calculate the PE (4) [116] . The fractal dimension D and the sub-leading term c are obtained from linear fits S(L 0 ) = DL 0 + c. To obtain D(L) and c(L) in a system size resolved manner, we use three chain lengths in the fit: L 0 = L − δL, L, L + δL, where δL = 2 for L < 40 and δL = L/10 for L > 40.

The fractal dimension D is shown in Fig. 2(a) . The wave function is fully delocalized (D → 1) over the Hilbert space only for p → 0. For p > 0, we observe a fractal scaling of the wave function (D < 1), and that the fraction of Hilbert space occupied by the wave function decreases monotonously with the measurement probability p. The numerical results suggest that D collapses on a limiting curve D(∞) with a flex point at p = p c for L → ∞ [113], exhibiting similar behavior to that of the fractal dimension at equilibrium quantum phase transitions [93] .

Crucially, the universal information about the MIPT is completely encoded in the sub-leading term c(L), which acts as an order parameter for the system. Indeed, c(L) approaches a step function, with discontinuity at the MIPT, for p = p S c ( Fig. 2 (b) ). Employing a scaling ansatz c = f [(p − p S c )L 1/ν S ] we find the data collapse on a universal curve with p S c = 0.160(2) and ν S = 1.28(5) in agreement with earlier results based on entanglement measures [69] .

The behavior of D and c in the p → 0 and p → 1 limits, highlighted with dashed lines in Fig. 2 , can be easily understood. In the limit of no measurements, we observe that our local Clifford circuit acts as a global Lqubit Clifford gate [117, 118] for which we analytically obtain S = L n=1 1 − (2 n + 1) −1 = L + c stab , where c stab ≈ −0.7645. In the opposite limit of p → 1, PE is vanishing after the layer of measurement. Subsequently, PE increases due to the layer of L/2 2-qubit gates yielding S = DL with D = 11 /15. (See [113] for details of these derivations, based on counting arguments and on properties of the Clifford group.) Our numerical analysis shows that the sub-leading term c matches these limiting values: c ≈ −0.7645 for p < p S c and c = 0 for p > p S c . Haar circuits. We investigate the MIPT by considering the PEs for q = 1, 2, 3 and focusing on the Z basis. To obtain D q and c q , we consider the linear fits S q (L 0 ) = D q L 0 + c q for three system sizes L 0 = L − 2, L, L + 2. As captured by the non-trivial q−behavior of D q (see Fig. 3 (a)), the Haar circuits exhibit multifractality for any p > 0. Furthermore, similarly to the stabilizer circuits, we observe a mild system size dependence of D q close to the MIPT, suggesting that D q (p) approaches a continuous curve in the L → ∞ limit with an inflection point at p H c . The sub-leading term c q approaches a step function with increasing L as shown in Fig. 3 (b) . The value of c q is non-zero and q dependent in the error correcting phase p < p H c , whereas it vanishes in the quantum Zeno phase p > p H c . The sub-leading term c q collapses onto universal (q dependent) curves upon rescaling p → (p − p H c )L 1/ν . The parameters p H c = 0.166(5) and ν = 1.4(1) are in agreement with the results in [18] obtained from the scaling of the quantum mutual information. As for the stabilizer circuits, the sub-leading term c q acts as an order parameter for the MIPT. The values of D q and c q in the limiting cases p → 0 (p → 1) can be easily understood. For p → 0 and for t > 2L, the PEs are well approximated by the PE obtained replacing the whole circuit by a matrix U drawn with the Haar measure from the SU (2 L ) group. This is expected since a sufficiently deep local Haar random circuit acting on L qubits forms an approximate unitary t-design [119] [120] [121] [122] . Assuming that the initial state |ψ in is the first basis state (which can be done due to translation-invariance of the Haar measure on SU (2 L )), we get that the PEs

, where E U denotes the average with Haar measure. The integral I q N ≡ E U log 2 N j=1 |U j,1 | 2q for U ∈ SU (N ) can be easily evaluated numerically for small N [123] , and for N 1 it is approximated by I q N = (1 − q) log 2 N + log 2 (Γ(1 + q)) [113, 124, 125] . This yields D q = 1 and c q = (1−q) −1 log 2 (Γ(1+q)). In the opposite limit, p → 1, the layer of measurements leads to a collapse of the state of the system onto a product state |ψ m = L i=1 |e i . A subsequent application of a layer of L/2 two-site gates U i,i+1 leads to a state |ψ =

, so that in the p → 1 limit, we get a non-zero D q and a vanishing sub-leading term c q .

Replica approach to participation entropies.

In the following, we provide an intuitive interpretation of the extensive scaling of PEs, and of c q as an order parameter for the transition. The entanglement entropy in Haar circuits can be calculated with a replica approach by means of a mapping to a classical statistical model [101] [102] [103] [104] [105] [106] . We generalize the Haar circuit to qudits with on-site Hilbert space dimension d and express the PEs (1) in terms of partition functions of a classical 2D spin model, translating the MIPT to an equilibrium ordering transition [113] . The bulk of the spin model is fixed by the spatiotemporal structure of the circuit (coinciding with the one implemented for the entanglement entropy computation [101] ), and the PEs are fully encoded in the choice of boundary conditions. An immediate consequence is that both PEs and entanglement entropy encode the same critical properties of MIPTs. Furthermore, this mapping provides a direct link between the quantum phase transitions in spin models at equilibrium [99, 100] , and our non-equilibrium setup: in both cases c q is an order parameter for the transition.

The calculation of the PEs simplifies in the limit of large on-site Hilbert space dimension d. At leading order in d and in the limiting cases of small and large measurement rate p, we obtain closed expressions: the fractal dimension is given by D q = log 2 d − p 2 2(1−q) log 2 Γ(1 + q), and the sub-leading term is

. These results are perfectly in line with our numerical findings, indicating the presence of multifractality and showing that D q is a decreasing function of p, whereas c q is constant throughout the error-correcting and quantum Zeno phase. The replica calculation of the PEs provides also an intuition for the order-parameter behavior of c q . For p → 0, we find that the non-zero value of c q originates from the correlation between spins that point in the same direction on all lattice sites. In contrast, for p → 1 the spins on different sites are fully uncorrelated, which results in a vanishing c q . Heuristically, this observation (valid for the limiting cases) can be extended to the whole phase provided the thermodynamic limit is taken. The region p < p c corresponds to a "ferromagnetic" phase of the spin model [102] , in which the correlation between spins throughout the 2D lattice leads to a non-zero value of c q . Instead, p > p c corresponds to a disordered phase of the spin model, in which the spins are uncorrelated beyond a certain length scale, giving rise to the PE strictly proportional to system size L and hence to the vanishing c q .

Universality. To show the generality of our results, we analyze PEs across MIPT in various settings: i) a Clifford circuit with rank-2 measurements onto local Bell pairs, ii) Haar circuit with on-site Hilbert space dimension d = 3, iii) a Floquet circuit with alternating application of U F = e −i j Xj e −i j (Zj Zj+1+Zj ) and measurements M j,± = (1 ± λZ j )/(2 + 2λ 2 ) (j ∈ [1, L]) with probability p = 1 and strength λ. In [113] we detail our numerical results which are summarized as follows. The sub-leading term, c q , plays the role of order parameter for MIPT in each of the cases i)-iii). For the Clifford circuit with rank-2 measurement, we found a collapse of c q with p S2 c = 0.692(2) (in agreement with [68] ) and ν S2 = 1.30(3) consistent with ν S according to the expectation that a local change of the measurement protocol does not alter universality class of MIPT. For ii): Haar circuit with d = 3, c q approaches a step function. The values of c q in the error-correcting phase are given by c q = (1 − q) −1 log 2 (Γ(1 + q)) match the results for d = 2, d → ∞ and, importantly, coincide with random matrix theory prediction [126] for systems with broken time reversal symmetry. The Floquet circuit iii) is composed of U F generated by local Hamiltonians without any randomness and of weak measurements M j,± with strength λ that allows to tune the system across MIPT [69] . Despite the differences with random quantum circuits, the behavior of D q and c q is analogous. The multifractal di-mension 0 < D q < 1 depends in a non-trivial fashion on q for λ > 0 confirming that multifractality of wavefunction can be generally expected in many-body systems with measurements. The sub-leading term c q , shown in Conclusions. We have shown that the PEs of quantum circuits with local measurements allow for a full characterization of MIPTs, alternative to entanglement measures. The expression for the PE of stabilizer states exhibits absence of multifractality. This emphasizes their simplicity as compared to generic many-body wave functions, and provides yet another perspective on the classical tractability of stabilizer circuits and Gottesman-Knill theorem. The multifractal dimensions, the leading terms in system size scaling of the PEs, indicate that the generic wave-functions are multifractal in the presence of projective measurements in the system. The sub-leading term c q encodes universal information about MIPTs, acting as an order parameter with value determined by whether the circuit is generic (e.g. Haar and Floquet circuits) or specifically fine-tuned (stabilizer circuits) and by gross features of quantum dynamics such as the dynamical phase or the presence of time reversal symmetry. This raises general questions about the connection of c q to ergodicity of quantum many-body systems and, in particular, about behavior of c q across many-body localization transition (cf. [83] ). Our results highlight the importance of the structure of many-body wave-functions for MIPTs, and can be used as a reference point in the vividly pursued investigations of MIPTs in many-body systems such as the Bose-Hubbard model [48, 49] , spin systems [50] [51] [52] [53] [54] [55] or free fermions [56] [57] [58] [59] [60] . Recent progress in stochastic sampling of wave functions in atomic as well as in solidstate platforms [127] [128] [129] [130] [131] [132] [133] [134] shows that direct probing of PEs constitutes an experimentally viable alternative to entanglement measures [70] . Random circuit sampling in superconducting quantum processors [135, 136] offers another possibility to study PEs experimentally (see [113] for analysis of an experimental setup). The multifractality at MIPTs was recently investigated also in a different context of: i) correlation functions [137] that exhibit multifractal scalings at MIPT in generic models ii) many-body wave-functions transformed to graph states [138] which can be mapped to Anderson models.

Supplemental Material: Universal behavior beyond multifractality of wave-functions at measurement-induced phase transitions

The Supplemental Material contains:

1. A discussion of the basis dependence of participation entropy for Haar and stabilizer circuits. This includes also a discussion on measuring the PE after a measurement layer, oppositely to the Main Text, where we discuss the PE computed before the measurement layer.

2. The behavior of PEs in stabilizer and Haar circuits in the limits p → 0 and p → 1.

A discussion on self-averaging, and the comparison with the multifractal analysis obtained through inverse participation ratio.

4. The crossover of the fractal dimension in stabilizer circuits toward the limiting thermodynamic form.

5. The replica calculation of the participation entropies, which in turns includes:

• details of the mapping of the circuit to a classical spin model;

• results for the limit of large on-site Hilbert space dimension d;

• cluster expansion results at p → 0 and p → 1.

6. Additional numerical analysis, on:

• Clifford circuits with rank-2 measurements;

• Haar circuits for d = 3 qudits;

• Floquet quantum circuits with generalized measurements.

7. An experimental setup proposal on state-of-the-art and near term quantum computers.

As we noted in the Main Text, the participation entropies

are dependent on the choice of the many-body basis {| σ }. Here, we describe this dependence for Haar circuits and stabilizer circuits, emphasizing the universal information about MIPT contained in the sub-leading terms c q in system size scaling of participation entropy S q = D q L + c q .

The numerical results presented in the Main Text for the Haar circuits were obtained in the Z basis (the eigenbasis of Z j operators) with participation entropies calculated after the layer of unitary gates. However, within this prescription, the the results remain exactly the same in the generic eigenbasis of local operators U * j Z j U j where U j are arbitrary (fixed at a given site) unitary matrices from SU (2) group. This is an immediate consequence of the translationinvariance of the Haar measure on the SU (4) group from which the 2 qubit gates U i,i+1 are drawn.

If instead the participation entropies are calculated after the layer of measurements, the freedom of local rotations is no longer present since the basis associated with the projectors P ± = (1 ± Z)/2 becomes distinct. The results for that case are shown in Fig. S1 .

The fractal dimensions D q , shown in Fig. S1 (a), (d) coincide only in the limit p → 0, when the wave-function is fully delocalized (D q → 1). Close to the critical measurement probability, we see a crossing point of D q curves for various L. In the Z basis, with the increase of p, the wave function is extended over decreasing fraction of Hilbert space and gets localized (D q = 0) in the p → 1 limit. Localization of the wave-function in the Z basis is equivalent to its full delocalization in X basis, as indicated by the value of D q for p → 1 in X basis. The results for D q show that the behavior of fractal dimensions is strongly dependent on the details of the protocol of their calculation. However, in all of the cases considered, D q demonstrate the multifractal scalings of wave-functions in presence of local projective measurements.

The sub-leading terms c q , shown in Fig. S1 (b-e) are independent of the choice of basis, containing the universal information about MIPT: the finite-size collapses with the scaling form c q = f [(p − p c )L 1/ν ] yield p c and ν that agree, within estimated error bars, regardless of the choice of basis.

The expression

was derived in the Main Text for participation entropy for the Z basis. An analogous formula can be derived for the X basis S2 shows the results of a direct test of the formulas (S2), (S3). Both results show that S q is independent of the parameter q and hence the wave functions of stabilizer states are not multifractal. The equations (S2), (S3) hold in the respective eigenbases of operators X k and Z k that determine the Pauli strings g j specifying the stabilizer state. In general, basis operators related by a Clifford unitary S3 transformation would preserve the absence of multifractality for the wave function of stabilizer states. Instead, generic unitary transformation U * j (e.g. Haar distributed) would induce a non-trivial q dependence and multifractality. 

In this section we provide details of derivations of formulas for PEs in both stabilizer and Haar circuits in the limiting cases of vanishing measurement probability p → 0 and large measurement probability p → 1.

Our derivation in the limit of no measurements (p → 0), is based on the observation that a local stabilizer (Haar) circuits acts as global L-qubit Clifford (Haar) gate. This can be expected for stabilizer circuits since an arbitrary Clifford operation on L qubits can be expressed in terms of a linear in L number of 2-qubit gates in nearest-neighbor architecture [117, 118] . For Haar circuits, such a behavior is also expected since a sufficiently deep local Haar random circuit acting on L qubits forms an approximate unitary t-design [119] . In contrast, for the large measurement probability p → 1 our derivations are based on the fact that for both stabilizer and Haar circuits the state is collapsed onto a Fock state by the layer of L measurements of the Z i operators.

To progress in the case of stabilizer circuits, we note that any operator U C from the Clifford group can be written as U C = F 1 HF 2 [139] where H is a layer of Hadamard gates, and where F 1 , F 2 reshuffle vectors of the Z basis up to a phase. Hence, H is the only term that influences the value of PE: For an initial product state, each Hadamard gate produces an equal superposition of |0 and |1 states, increasing the value of PE by 1. Therefeore, the average PE after the action of U C is equal to the expected number of Hadamard gates in the Clifford operation. If U C is an L-qubit Clifford gate, the probability of having a Hadamard gate acting on qubit n is p n = 2 n /(2 n + 1) [139] . The PE is given by (2))/ ln(2) ≈ −0.7645. (Here ψ q is the q-th polygamma function). In the opposite limit of p → 1, the state is collapsed on a single Fock state and the PE is vanishing after the layer of measurements. Subsequently, PE increases due to the layer of L/2 2-qubit gates. Each of them has the average number of Hadamard gates equal to 22/15 (as implied by the expression p n = 2 n /(2 n + 1)), hence PE is given as S = DL with D = 11/15 and a vanishing sub-leading term c = 0. The numerical analysis from the Main Text shows that the sub-leading term c in stabilizer circuits matches the limiting values for p → 0 and p → 1 throughout the entire quantum error correcting and quantum Zeno phases.

For the Haar circuits, in the p → 0 limit, we calculate the PEs by replacing the whole circuit by a matrix U drawn with the Haar measure from the SU (2 L ) group. Using the translational invariance of Haar measure on SU (2 L ), the PEs are given by S p=0

where U ∈ SU (N ) and E U denotes the average with Haar measure over the unitary group SU (N ) can be easily evaluated numerically by noting that |U j,1 | 2 are uniformly distributed on a standard N -simplex ∆ N = {λ ∈ R N : [123] . This can be done with Monte Carlo method by noting that Dirichlet distribution Dir(λ 1 , . . . , λ N ; a 1 = 1, . . . , a N = 1) is a uniform distribution over ∆ N . The calculation simplifies for N 1, when the real and imaginary parts U j,1 behave as uncorrelated Gaussian random variables with a vanishing mean and variance equal to N −1 [124, 125] . Consequently, t j = |U j,1 | 2 are independent random variables distributed with the exponential distribution P (t j ) = N e −N tj . Introducing a variable w = j t q j we note that w, as a sum of N 1 independent random variables, is normally distributed with mean proportional to the Gamma function µ = N 1−q Γ(1 + q) and standard deviation σ = N 1−2q Γ(2q + 1) − Γ(q + 1) 2 . For N 1, we have σ µ, and log 2 (w) ≈ log 2 (µ), which yields the results

This leads to the following expression for the PEs:

which approximates well the results already for moderate size circuits (e.g. L = 8, since 2 L 1). From (S6) we get that the fractal dimension D q = 1 and the sub-leading term c q = (1 − q) −1 log 2 Γ(1 + q). In contrast, as we argue in the Main Text, the PEs in the limit p → 1 are given by S p=1 q = L 2(1−q) I q 4 and the sub-leading term is vanishing. The expressions derived here for the sub-leading term c q in the p → 0 (p → 1) limit match the numerical results in the entire quantum error correcting (quantum Zeno) phase.

Another quantity of employed in the literature to probe multifractal properties of the wave function, particularly relevant for Anderson transitions [79] , is the Inverse Participation Ratio (IPR), defined as

The typical value of IPR, IPR typ q ≡ exp log IPR q (where . denotes the average over ensemble of states |ψ m ) is functionally dependent on the average value of PEs: 1 1−q log(IPR typ q ) = S q . Therefore, the analysis of the average PEs, conducted in the Main Text, is equivalent to an analysis of the typical value of the IPR. In this section we investigate behavior of the average IPR: IPR q . The average IPR can be written as [79] : IPR q = a 0 2 (1−q)DqL whereD q is a fractal dimension, a 0 is a proportionality factor, we have neglected sub-leading terms and taken into account that the dimension of the Hilbert space is 2 L . Taking a logarithm of this expression we observe that

where we have introducedc q = (1 − q) −1 log 2 a 0 . We analyze this expression for stabilizer circuits in system size resolved manner, similarly as for the average PEs in the Main Text. For stabilizer states, using the result that S ≡ S q = rk Z2 (M X ) and (S7), we immediately find that IPR q = 2 (1−q)rk Z 2 (M X ) in the Z basis. We report the numerical results for Clifford circuits in Fig. S3 , where we compare the fractal dimensionD (D) and the sub-leading termc (c) extracted from the IPR (PE). The results are averaged over more than 20000 circuit realizations. Despite mild deviation in the precise value of D and c (compared toD andc) the qualitative features and the critical scaling at the MIPTs are the same. This suggests that the system exhibits a qualitative self-averaging and the information extracted from the average IPR and typical IPR (equivalent to PEs) is the same. Both D q andD q exhibit similar, mild finite size effects at the MIPT as visible in Fig. S3 a) , c). Moreover, the sub-leading termc q approaches a step function at MIPT with increasing L, analogously to c q , as shown in Fig. S3 b) , d). For both cases p c = p S c = 0.160(5) and ν = 1.30(9) (note that the error bars are larger than in the Main Text do tue a smaller interval of system sizes considered). Similar considerations extend also to Haar circuits, but we do not report the numerical results for the sake of readibility. We conclude by noting that the exponential in Eq. (S7) gives rise to larger errors in the data, which is particularly visible in the sub-leading termc q shown in Fig. S3 b) . This justifies our choice of considering the participation entropy instead of the IPR to investigate the multifractal and the MIPTs in the quantum circuits of interest. Figure S3 .

Comparison ofD (a) andc (b) obtained from the system size dependence of IPR: (1 − q) −1 log 2 IPRq = DqL + (1 − q) −1 log 2 a0 with Dq (c) and cq (d) obtained from the PEs: S = DL + c. The information about MIPT encoded both in the IPR and in the PE is the same (although there are minor differences: for instance, the values ofc and c are different in the error-correcting phase.

The numerical results suggest that the fractal dimension D collapses on a limiting curve D(∞) with a flex point at p = p c for L → ∞. This is presented in Fig. S4 , where we analyze the L dependence of D considering D(L) − D(L ref ) and taking the reference system size L ref = 8. Beyond a length scale L * , the value of D(L) is a good estimate for D(∞) (e.g. L * ≈ 40 at p = 0.1). The value of L * increases as p gets closer to p c ; however, outside of the interval of p ∈ (0.15, 0.17), the signs of saturation of D(L) are clearly visible. This suggests that the system size dependence in D(L) are crossover effects which are washed away in the thermodynamic limit and that the D(∞) possesses a flex point at p = p S c , exhibiting similar behavior to that of the fractal dimension at equilibrium quantum phase transitions [93] . However, our numerical results do not allow us to unambiguously exclude the possibility of a non-trivial asymptotic scaling of D(L).

In this section, we detail the computation of the average participation entropies (PEs) for a Haar random circuit acting on d-dimensional qudits. We first introduce the mapping from the random circuit to a classical statistical model and then obtain closed expressions for the PEs in the limit of d 1.

This subsection follows [101] [102] [103] [104] [105] [106] , where the mapping from the Haar hybrid circuit to a classical statistical mechanics model has been already considered. Each quantum trajectory |ψ m = K m |ψ /||K m |ψ || is specified by the action of the non-unitary operator K m describing the circuit realization m on an initial state |ψ . We are interested in the average PEs over the quantum trajectories

which we map below to a free energy cost of classical statistical mechanics model. Notice that here, for convenience, we consider the natural logarithm in Eq. (S9) instead of the log 2 used for d = 2. This will change the result by an overall multiplicative factor ln (2) A key remark is that the bulk of the classical statistical model is fixed by the spatiotemporal structure of the circuit, whereas the observable of interest is fully encoded in the boundary condition of the model. As a consequence, for the problem of interest, the degrees of freedom, the geometry of the lattice and the interaction terms are the same of those presented in, e.g., [101] .

The mapping consists of two steps. First, the PEs of a single trajectory are expressed as a trace over an operator acting on q-copies of the system (Rényi replica). Then the logarithm entering the participation entropy is expanded introducing k copies of the Rényi replicated systems. The final model involves Q = kq + 1 layers, with the additional term coming from the circuit probability density ||K m |ψ || 2 .

We define the boundary operator Λ q such that

(S10)

With this notation, the average participation entropy is given by

(S11)

In the above expression we have defined the partition functions

Pictorially, we can introduce the layer representation of a circuit realization (S14) which recasts Eqs. (S12) and (S13) to .

(S15)

Since Z Λ = Z 0 = 1 in the replica limit k → 0 (Q → 1), Eq. (S11) becomes a free-energy cost associated with change of the boundary condition

where F α = lim k (ln Z α )/(k(1 − q)) are the free energies for the boundary condition α = 0, Λ. Since both measurement and unitary gates are independently and randomly distributed, we factorize the average E m . For the unitary gates we have

where Wg D (s) are the Weingarten symbols, and s, r are permutations over the symmetry group S Q [101] . The operators χ s have matrix elements

δ(σ (r) , τ (s(r)) ).

(S18)

The average over the measurements can be performed in a similar fashion. It is convenient to consider the measurements contracted with the operators in Eq. (S17), leading to (S19) S8 where the square indicates M ⊗Q , and W p (s) = (1 − p)d C(s;S Q ) + pd and C(s; S Q ) is the number of cycles in the permutation s in S Q . The action on the initial (product) state |ψ = ⊗ L a=1 |e a is given by

Instead, the boundary condition α is fixed by the operator acting on the top layer of the circuit

where the contraction signals a sum over the associated degrees of freedom. In particular, we have

where the delta functions in the Λ (pictorially represented in the magenta dots in Eq. (S21)) impose the "book" structure on the k replicae (see Eq. (S15)). After the average is performed, the partition function becomes a classical spin model on an anisotropic honeycomb lattice Ξ

In Eq. (S23), the spins are elements of the symmetric group S Q , and we have introduced the notation Ξ to define the red links, Ξ the blue links, and ∂Ξ the boundary sites. The lattice contains N = L(t − 1)/2 spins for a qudit chain of length L and for a total time t 1.

The expression Eq. (S23) involves no approximation, and is valid for any d. However, the full computation of the partition function Z α (k, q) is involved, due to frustration effects in its weights. To obtain analytic insights we therefore consider the limit of large qudit dimension d 1. At leading order we have

where we have introduced the function

and

S9 and 1 is the identity permutation. The leading order of W bnd Λ (s) is a consequence of the permutation invariance of each book page. The sub-leading corrections in d arise when larger groups are considered. For instance, S 2q ⊗ S ⊗(k−1) q ⊗ 1 would arise when the first two book share the same boundary state, which is equivalent to consider them merged in a single book with 2q pages. Another way to obtain a sub-leading in d correction is to add the Q-th replica to the k-th book (hence forming a new book of q + 1 pages), which gives rise to group S ⊗(k−1) q ⊗ S q+1 . Hence, up to corrections in d −2 , the honeycomb lattice reduces to a square one ,

with 2N links. The expression for the partition function becomes

where the Boltzmann weights can be easily expressed from W p (s) and W bnd α (s), and are given by

(Notice that here we considered Q ≥ 1). This expression shows that | ln p| plays the role of an inverse temperature. Hence, at the leading order the model is a Q!-state Potts model defined on the structure Eq. (S28). Neglecting the constant energy term, which simplifies in the ratio Eq. (S16), we have

Cluster expansion at p → 0 and p → 1

We conclude this section by considering the limit of low and high measurement rate, which correspond, respectively, to the low (β → ∞) and high (β → 0) temperature limits.

At low temperature, the weights are zero, except when the v i,j = 0. Hence, we can approximate

We expand the partition functions in a series in p: q) . A simple combinatoric computation gives

Here the combinatoric constant in Z Λ (k, q) comes from the number of permutations satisfying the constraint Eq. (S25) (i.e. s ∈ S ⊗k q ⊗ 1) and from the condition of having the δ-function on each link. We note that this condition enforces a ferromagnetic-like phase, with the bulk degrees of freedom fixed by the boundary ones. The first sub-leading correction is second-order in p and is given by

Therefore, in the limit p → 0 we have

where we reabsorbed the ln 2 from Eq. (S9) for comparison with the Main Text. Importantly, we note that the value of the sub-leading term c q = log 2 (Γ(1 + q))/(1 − q) is not affected by the increase of the measurement rate p (within the range of validity of the second order approximation in p). This is consistent with the numerical evidence presented in the Main Text, supporting c q = 0 within the error-correcting phase.

In the opposite limit of infinite temperature (p → 1), it is convenient to express the Boltzmann weights as

an expand the partition functions in a series in (1 − p) as Z α (k, q) = 2N n=0

(1−p) n n! Z (n) α (k, q). In this limit, the spin system is in a paramagnetic phase: the permutations on different lattice sites are independent of each other. Hence, the partition function is approximated by an equal weight sum over all the possible degrees of freedom. A combinatorical computation gives

Similarly, it is possible to compute the first sub-leading corrections, which are order

Hence, at high temperature, the PEs reads

where, again, we reabsorbed the ln 2 from Eq. (S9) for comparison with the Main Text. Our computation is consistent with the numerical analysis for d = 2, and in particular support the stability of c q = 0 within the quantum Zeno phase p > p c .

In order to test the universality of the behavior of the sub-leading term c q in scaling of PEs, in this section we numerically investigate three additional setups: (i) rank-2 measurement stabilizer circuits, (ii) d = 3 qudit Haar circuits, and (iii) Floquet quatum circuits with weak measurements. Throughout this section, we consider the PEs in the Z-basis computed after the layer of unitary gates. Figure S5 . Scheme for the stabilizer circuits with the rank-2 measurements P

(2) ± . The unitary gates Ui,i+1 are drawn from the Clifford group as in Fig. 1 of the Main Text.

We consider the stabilizer circuit in Fig. S5 . The difference with respect to the setup in Fig. 1 of the Main Text is the substitution of the rank-1 measurement P (1) i,± ≡ P i,± = (1 + Z i )/2 by rank-2 measurements P (2) i,± = (1 ± Z i Z i+1 )/2. These operators project onto the following Bell states

We note that the rank-2 measurements have less disentangling power compared to the rank-1 ones, as the residual entanglement of the Bell states in Eq. (S42) is ln 2. Therefore the critical point of this system is expected to occur to higher measurement rates. This setup has been previously considered in Ref. [68] . We note that the exponent ν estimated in that paper (ν 1.75) does not fit later and more accurate studies for the stabilizer circuits with rank-1 measurements (ν = 1.30(5) [70] ). (The reason for this discrepancy is due to a smaller accuracy of the entanglement entropy collapses considered in Ref. [68] in comparison to the entanglement measures considered in [69, 70] .) Thus, our analysis has an additional benefit. It extracts a precise location of the critical point and an accurate value the correlation length exponent which have been not yet reported for stabilizer circuits with rank-2 measurements.

The numerical results, averaged over more than 20000 circuit realizations, are given in Fig. S6 . We see that, similarly to the rank-1 case, the fractal dimension D converges at large system sizes, to a universal curve, whereas the sub-leading term c approaches a step-function for increasing system size L, with c = c stab −0.76 for p < p c and c = 0 for p > p c . The estimated critical point is p c = 0.693(2) and the correlation length critical exponent is ν = 1.30 (3) .

We remark that both the value of the sub-leading term c in the error correcting phase, and of the exponent ν at the critical point match those for the rank-1 measurements discussed in the Main Text. This shows that the stabilizer circuits with rank-1 and rank-2 measurements lie in the same universality class, and that c, as a property of the phase, only depends on gross features of the model (such as the state being a stabilizer or not). (S43) S12

To the best of our knowledge this setup has not been previously discussed in the literature, therefore it constitutes a novel test for the universality of Haar circuits. Importantly, it is expected that the qudit dimension d is a relevant perturbation for the MIPT in random Haar circuits, since the limit d = ∞ differs from the finite d [13] .

In our numerical simulations we consider systems of up to L = 16 sites. This corresponds to Hilbert space dimension N = 3 16 = 43046721, larger that the largest Hilbert space dimension for d = 2 case considered by us or in [18] . Our findings are gathered in Fig. S7 . The results are averaged over more than 2000 disorder realizations. While, we believe that the system size is insufficient to extract a precise position of the critical point and of the exponent ν, we obtain a good collapse of the data for the sub-leading term c q with p c = 0.31 (5) and ν = 1.4 (3) . We note that the exponent ν is compatible with the numerical result for d = 2 given in the Main Text.

Furthermore, we note that the value of c q in the error correcting phase matches those estimated for d = 2 and for d → ∞ in our paper: c q = (1 − q) −1 log 2 Γ(1 + q). Hence, we conclude that also for Haar circuits the characterization of the phase through c q depends only on the gross features of the circuits. This value of c q precisely coincides with the result for Gaussian Unitary Ensemble (GUE) of random matrices [126] , appropriate for systems with broken time reversal symmetry. It is, however, important to note that the GUE result D q = 1 and c q = (1 − q) −1 log 2 Γ(1 + q) describes our systems only at p → 0. For p > 0, the wave-functions in random circuits with measurements become, in general, multifractal (D q < 1) but the sub-leading term persist to be equal to c q = (1 − q) −1 log 2 Γ(1 + q) in the whole error-correcting phase p < p c . Since there are no works in the literature discussing the Haar circuits with d = 3, it is interesting to compare the results obtained from the analysis of PEs with the results obtained through entanglement measures. Specifically, we consider the bipartite mutual information I 2 and the tripartite mutual information I 3 , which are considered as precise observables to locate the critical point in random circuits [18] .

These quantities are both defined in terms of the entanglement entropy. Given a state |ψ and a bipartition M ∪ N , the entanglement entropy between the subsystems M and N is given by

We consider a partition of our system (with periodic boundary conditions) into four adjacent (and not overlapping) subsystems A, B, C, D, each of length L/4, so that the subsystems A and C are antipodal. We calculate the bipartite quantum mutual information between subsystems A and C:

where S(A) is the entanglement entropy (S44) with M = A and N = B ∪ C ∪ D. The tripartite mutual information is given by

We compute Eq. (S45) and Eq. (S46) after the unitary gate layer, and consider the average value over the trajectories. Our numerical results are given in Fig. S8 , results are averaged over more than 400 disorder realizations. For the considered partitions, the critical point is identified by a crossing point [18] . (Differently from the choice given in Ref. [69] , where the transition is identified by a peak. The different behavior stems from the scaling of the subsystems size with the total system size.) From the data, we find that the crossing point for the bipartite mutual information is shifting toward smaller values of p with increasing system size L, suggesting large finite size effects. On the other hand, the tripartite mutual information exhibits a crossing point in the vicinity of p 0.24, which is roughly compatible with our estimate through the PEs. However, the limited sizes do not allow for a proper finite size scaling. Hence these results serve only as a consistency check. 

Lastly, we consider Floquet quantum circuits with weak measurements. This framework provides two tests. On one hand, we see if the Haar circuits are representative of generic quantum evolution. On the other hand, we also test if weak measurements change the physics of the system. This framework has been previously considered in Ref. [69] ; however our numerical analysis gives better estimates of the position of the critical point and of the exponent ν, as we extend the numerical simulations over larger system sizes and consider finer observables for the transition.

The circuit has a similar architecture to that given in Fig. 1 of the Main Text, but here the layer of two-qubit unitary gates is replaced by the Floquet operator

while the projections P ± are substituted with M ± = (1 ± λZ j )/(2 + 2λ 2 ) which act on each qubit (measurement rate p = 1) with strength determined by the parameter λ.

We consider both the PEs as well as the tripartite mutual information (S46) with the same partition of the system as for the Haar circuits. Our results, averaged over more than 1000 circuit realizations, are given in Fig. S9 and in Fig. 3 of the Main Text.

The estimates for the critical measurement strength λ c and for the critical exponent ν both for the sub-leading term c q and for the tripartite quantum mutual information are compatible and are given by λ c = 0.320 (8) and ν = 1.3 (1) . Therefore, we conclude the Floquet circuits with weak measurements also belong to the same universality class of Haar circuits with d = 2 qubits and projective measurements. Furthermore, we observe that in the error correcting phase, for λ < λ c , the sub-leading terms in scaling of PEs are given by c q = (1 − q) −1 log 2 Γ(1 + q), and are the same as for the Haar random circuits. The value of c q mathches the GUE prediction (note that the time reversal symmetry is broken in the Floquet circuits by the presence of measurements). This suggests that c q plays the same role in this phase as for the Haar error correcting phase.

In this section we simulate an experiment that observes the change of the sub-leading term c q in system size scaling of PEs across MIPT. This allows us to establish the experimental relevance of our findings by arguing that systems of size L = 10 seem to be well in reach of our method on superconducting quantum processors (this system size should be compared with recent experiment performed in systems of up to L = 8 sites [44] ). The following analysis demonstrates also that the experimental resources needed to observe MIPT through the PEs scale exponentially with the size of the system.

In numerical calculations, the many-body wave function ψ(σ) = σ|ψ of state |ψ in many-body basis {|σ } is directly stored in the memory of computer and the calculation of PEs

is straightforward (but still requires to sum the number of terms that is exponentially large in system size L). In turn, in atomic and solid state platforms [127] [128] [129] [130] [131] [132] [133] [134] as well as in superconducting quantum processors setups [135, 136] one can measure local observables (e.g. spin-z components Z i ) which leads to a collapse of the state |ψ of the system on one of the basis states |σ with probability | σ|ψ | 2 . Multiple repetitions of such a set of measurements on the state |ψ allow to estimate the probabilities | σ|ψ | 2 for all basis states |σ and to calculate the PEs with the formula (S48). Notably, this procedure of sampling of the final state of quantum circuit |ψ was employed in the so called cross-entropy benchmarking in experiments that demonstrated quantum advantage [135, 136] . In those experiments samples of over 10 7 output states |σ were collected, giving us a rough estimate of the number of samples that could be collected in a hypothetical experiment that aims to observe the change of the c q term across MIPT. In the following, we describe how the PEs (S48) can be accurately estimated in such an experiment with superconducting quantum processor.

Our simulation of the hypothetical experiment consists of the following steps: In our simulation the 2-qubit unitary gates are fixed for each set of local measurements encoded in the index m (an experimentally more convenient route would be to consider a fixed set of 2-qubit gates). Each measurement layer consists, on average, of pL local spin measurements, giving rise to approximately 2 pLt different possible states |ψ m at the output of the circuit. Hence, after the first run of the experiment (steps 1.-3.), the probability to end up with the same state |ψ m is 2 −pLt . To gather information about the state |ψ m it is necessary to perform multiple destructive measurements on it. This exponential in system size barrier in preparation a given state |ψ m is the so-called problem of post-selection. The problem is common for all experimental setups that aim to observe MIPT (one possible solution is to restrict the class of circuits to circuits with space-time duality [24] ). As we show below, the number of measurement outcomes |σ for each state |ψ m required to accurately estimate PEs scales exponentially with system size L. This scaling, combined with the fact that the number of possible final states |ψ m also scales exponentially with L results in the overall exponential scaling with system size of the number of repetitions n s of the experiment.

We now assume that by the measurements in step 3. we gathered, for each state |ψ m , in total N MC basis states |σ ("wavefunction snapshots"), distributed in such a way that state |σ occurs i σ times in the gathered output (note that this procedure is actually a Monte Carlo sampling of the wave-function, hence the subscript "MC"). Then, the probability | σ|ψ | 2 ≈ i σ /N M C , and the PEs (S48) can be approximated as

The approximations of S q from our numerical experiment, averaged over 10000(100) final states |ψ m for L = 6, 12 (L = 20) are shown as a function of the number of snapshots N MC in Fig. S10 a) , b), c). We observe that S q (N MC ) is an increasing function of N MC and that the increase persists even for N MC significantly larger than the Hilbert space dimension N = 2 L . However, the functional dependence S q on N MC is rather simple. For small system sizes (L = 4, 6, see Fig. S10 a) ) we observe that S q (N MC ) ∼ 1/N MC . Fitting a first order polynomial in 1/N MC in interval N MC ∈ [250, 1500] and extrapolating it to N MC → ∞ allows us to reproduce the exact value of S q (calculated with the full knowledge of many-body wave-function with formula (S48)) with accuracy better than 0.1%. For larger system sizes (L ≥ 8) we observe that S q (N MC ) ∼ 1/ √ N MC in the interval from N MC ≈ N /8 to N MC ≈ N and that for larger N MC we again have S q (N MC ) ∼ 1/N MC . This observation motivates the following extrapolation scheme.

• The curve S q (N MC ) is fitted by a function f (N MC ) = a/ √ N MC + b (where a, b are fit parameters) in an interval N MC ∈ [N min MC , N max MC ]. • For N MC > α q N we define a function g(N MC ) = a 0 (1/N MC − 1/(α q N )) + f (α q N ), where a 0 is chosen in such a way that dg(N ) dN | N =αqN equal to df (N ) dN | N =αqN ; α q is fixed for all system sizes for given q (we take α 1 = 2 and α 2 = 0.7).

• The desired value of PE is obtained by the extrapolation: S q = lim NMC→∞ g(N MC ).

The extrapolation scheme is illustrated in Fig. S10 b) , c) where the fits and the functions f (N MC ) and g(N MC ) are denoted by dashed lines. The scheme allows to estimate the value of PEs with accuracy of about 0.1% for all system sizes considered. The experimental resources needed for the determination of S q are proportional to N max MC which scales exponentially with L as visible in Fig. S10 but remains bounded from above by the Hilbert space dimension N = 2 L . A rough estimate of the total number of repetitions n s of steps 1.-3. can be obtained as follows: there is approximately ≈ 2 pLt possible final states |ψ m (assuming there is no randomness in the 2-qubit gates), for each of them one needs to perform at least ≈ N max MC ≈ 2 L measurements. This gives n s ∼ 2 pLt+L . For instance, assuming that for system size L = 10 and for measurement probability p = 0.1 the circuit of depth t = 10 is sufficient to reach the saturation value of S q , we get n s = 2 20 ≈ 10 6 . For larger values of p, the steady state value of S q is reached for a circuit of lower depth (for growing p the value of t needed to reach the steady state decreases) and the number of samples n s ≈ 10 6 seems to be a reasonable estimate for any value of p. Such a number of samples can be gathered easily with superconducting quantum processors (more than 10 7 samples or output were gathered in [136] for much larger system size). This should be compared with systems of up to L = 8 sites considered in the recent experiment [44] . However, n s increases very quickly with L: an analogous estimate for p = 0.1 and L = 14 gives n s ≈ 10 10 . If the problem of post-selection was eliminated (note that this problem does not occur at p ≈ 0), system sizes beyond L = 20 could be reached with present day superconducting quantum processors.

An open systems approach to quantum optics

The theory of open quantum systems

Quantum noise: a handbook of Markovian and non-Markovian quantum stochastic methods with applications to quantum optics

Quantum measurement and control

Since we are interested in the stationary state, our results are independent of the choice of initial states

Quantum computation and quantum information

In the main text we consider the former case and in [113] we show that differences involve only non-universal quantities

quantum circuit consisting of t = L unitary and measurement layers acts on the state (we employ the Haar circuit for d = 2 studied in the Main Text, however, other circuits, consisting of gates more straightforwardly accessible in superconducting quantum processor setups

spin operators Z i are measured in the final state |ψ m (where m specifies the circuit realization) collapsing the state to basis state |σ

the steps 1.-3. are repeated multiple times, and the basis states |σ are collected for each final state |ψ m