Entanglement dynamics in monitored Kitaev circuits: loop models, symmetry classification, and quantum Lifshitz scaling
Kai Klocke Daniel Simm Guo-Yi Zhu Simon Trebst Michael Buchhold







These two authors contributed equally.
Department of Physics, University of California, Berkeley, California 94720, USA

These two authors contributed equally.
Institut für Theoretische Physik, Universität zu Köln, D-50937 Cologne, Germany

Institut für Theoretische Physik, Universität zu Köln, D-50937 Cologne, Germany
Advanced Materials Thrust, The Hong-Kong University of Science and Technology (Guangzhou), Guangzhou, China


Institut für Theoretische Physik, Universität zu Köln, D-50937 Cologne, Germany

§ ABSTRACT 
    Quantum circuits offer a versatile platform for simulating digital quantum dynamics and uncovering novel states of non-equilibrium quantum matter. 
    One principal example are measurement-induced phase transitions arising from non-unitary dynamics in monitored circuits, 
    which employ mid-circuit measurements as an essential building block next to standard unitary gates.
    Although a comprehensive understanding of the dynamics in generic circuits is still evolving, we contend that
    monitored quantum circuits give rise to robust phases of dynamic matter,
    which – akin to Hamiltonian ground state phases – can be categorized based on circuit symmetries and spatial dimensionality.
    To illustrate this concept, we focus on measurement-only quantum circuits within symmetry classes BDI and D, 
    which are measurement-only circuit adaptations of the paradigmatic Kitaev and Yao-Kivelson models, 
    embodying particle-hole-symmetric Majorana fermions with or without time-reversal.     
    We establish a general framework – Majorana loop models – for both symmetry classes (in arbitrary spatial dimensions) 
    to provide access to the phenomenology of the entanglement dynamics in these circuits, 
    displaying both an area-law phase of localized Majorana loops and a delocalized, highly entangled Majorana liquid phase. 
    The two phases are separated by a continuous transition displaying quantum Lifshitz scaling, 
    albeit with critical exponents of two distinct universality classes.
    The loop model framework provides not only analytical understanding of these universality classes in terms of non-linear sigma models,
    but also allows for highly efficient numerical techniques capable of simulating excessively large circuits with up to 10^8 qubits. 
    We utilize this framework to accurately determine universal probes that distinguish both the entangled phases and the critical points 
    of the two symmetry classes. Our work thereby further solidifies the concept of emergent circuit phases and their phase transitions.



Entanglement Dynamics in Monitored Kitaev Circuits:  Loop Models, Symmetry Classification, and Quantum Lifshitz Scaling
    Michael Buchhold
    September 9, 2024
=======================================================================================================================









§ INTRODUCTION




Many-body entanglement has emerged as a fundamental concept in understanding and classifying quantum phases of matter <cit.>.
In an equilibrium setting, it can serve as a tool to discern ground states of gapped from gapless Hamiltonians, to identify topological order, and to distinguish between generic and non-generic excited states <cit.>. Recently, this notion of many-body entanglement has been extended beyond static Hamiltonian states to characterize dynamically generated quantum states, in particular those arising in monitored quantum circuits. Analogous to Hamiltonian systems, robust entanglement structures, such as for instance area-law or volume-law scaling, arise in quantum circuits and can be used to characterize the generated quantum states <cit.>.


The robustness of these entanglement structures against small variations in microscopic circuit parameters defines entanglement phases of matter, akin to phases of matter in thermal equilibrium, which depend only on global properties, such as underlying symmetries, the range of quantum gates and measurements, and the spatial geometry. This perspective of defining entanglement phases in quantum circuits is intricately connected to the concept of universality in statistical mechanics – stating that macroscopic features of a many-body system are governed by a few global properties, such as symmetries, interaction range, and dimensionality – and it anticipates a symmetry classification of entanglement phases in monitored quantum circuits.


Along these lines it has been recognized that in generic quantum circuits with measurements, the presence or absence of a U(1)-symmetry, associated with the conservation of the total particle number or magnetization, dramatically alters the potential phases <cit.>. However, a comprehensive classification of circuit phases in terms of their symmetries is presently absent. Here we provide a first step in this direction, analyzing two symmetry classes of measurement-only circuits in two spatial dimensions. In addition, we provide a powerful framework for the analytical and numerical analysis of said circuits and symmetry classes in arbitrary spatial dimensions.


In a quantum circuit with mid-circuit measurements, monitored circuit for short, the global symmetries are inherited from the symmetries of  the combined set of generators of the dynamics, i.e., from all measurement projectors and quantum gates that evolve the wave function.  In this work, we focus on systems characterized by particle-hole symmetry, expressed through operators that are even in Majorana fermions. Specifically, our circuits are  measurement-only circuit analogs <cit.>  of the paradigmatic Kitaev model <cit.> (and variations thereof <cit.>), where frustration arises not from non-commuting Hamiltonian terms but from non-commuting two-qubit parity check measurements, see Fig. <ref> for an illustration. 

In terms of symmetry, these “Kitaev circuits" adhere to C^2=1 particle-hole symmetry and, depending on the underlying spatial geometry of the circuit, fall into two categories with regard to time reversal symmetry. For bipartite spatial geometries (such as the honeycomb geometry) the system exhibits T^2=1 time-reversal symmetry and falls within symmetry class BDI of the ten-fold way classification of free-fermion systems <cit.>.
For non-bipartite spatial geometries (such as the next-nearest neighbor Kitaev <cit.> or the decorated honeycomb geometry which Yao and Kivelson have first explored as Hamiltonian model <cit.>), time-reversal symmetry in the ground state is broken and the system belongs to symmetry class D.
For the Hamiltonian systems, this distinction with regard to time-reversal symmetry has crucial impact on their (spin liquid) ground states. The Yao-Kivelson model, for instance, allows for the formation of a chiral spin liquid which is absent in the conventional honeycomb Kitaev model.
Whether a similar distinction arises also for the quantum states stabilized by the corresponding measurement-only circuits in symmetry class BDI and D is a central question that we address in this manuscript.



For both symmetry classes, e.g., the Kitaev honeycomb circuit or its generalizations including next-nearest neighbor couplings or the Yao-Kivelson geometry, we find that when the measurement-induced frustration is large, i.e., when non-commuting operators are measured frequently, each circuit creates an entangled Majorana liquid state. It displays a distinctive subsystem entanglement entropy S(A)∼ Llog(L) for a subsystem A=L× L of linear dimension L. In contrast, when the measurement-induced frustration is low, a weakly entangled area-law state is realized. The two phases are separated by a critical line, which displays a characteristic quantum Lifshitz scaling behavior <cit.>, previously associated with quantum dimer models <cit.>, Dirac fermions and certain (2+1)-dimensional conformal field theories (CFTs) <cit.>.
But while the general entanglement phase diagrams appear almost identical for the two systems, their underlying symmetry does manifest itself, primarily at the phase transition between the two principal entanglement phases which fall into distinct universality classes whose critical exponents we determine.



To unveil this physics, we introduce Majorana loop models in (d+1) dimensions as a joint theoretical framework which covers both the Hamiltonian and measurement cases. These loop models naturally encapsulate the dynamics of Majorana world lines, distinguishing between symmetry classes BDI and D through the orientability or non-orientability of world lines <cit.>. Moreover, they differentiate between Hamiltonian ground states and monitored systems using an internal parameter – the loop fugacity (which, in topological terms, is also referred to as d-isotopy parameter) <cit.>. 
The fugacity n has an intuitive interpretation: it determines whether all world lines appear with equal probability, yielding measurement-only circuits (n=1) or whether world line configurations corresponding to high energy states are suppressed, yielding Hamiltonian ground states (n=√(2)) <cit.>.






The loop framework gives access to both analytical arguments and efficient numerical simulation techniques for monitored Majorana circuits. This allows for a detailed and thorough analysis of the statistical mechanics of circuits in symmetry classes BDI and D in two spatial dimensions and a detailed exploration of Lifshitz criticality in a setting of free, monitored fermions. It is enabled through new, optimized simulation techniques for loop models, which yield a numerical complexity scaling as 𝒪(Nlog(t)) for N∝ L^2 numbers of qubits and circuit depth t. In numerical simulations, system sizes of 10^8 qubits can be realized, reaching the largest system sizes currently accessible for monitored circuits and raising the standard for the quantitative determination of critical exponents and correlation functions in two spatial dimensions. In addition, we simulate the monitored Kitaev and Yao-Kivelson circuits of interest in this manuscript also in the conventional stabilizer representation. We employ state-of-the-art Clifford circuit simulations techniques <cit.> that reach system sizes of 10^4 qubits and critically analyze what kind of qualitative and quantitative insight they afford.

In summary, we discuss measurement-only Kitaev and Yao-Kivelson models as representatives of symmetry classes BDI and D in monitored circuits. We establish a joint framework for both classes – Majorana loop models – and utilize it to accurately determine the long-wavelength properties of a localized area-law and an entangled Majorana liquid phase, and the measurement induced phase transition separating them. This is a rigorous step towards the symmetry classification of  monitored circuits in arbitrary dimensions. In addition, we establish that the measurement-induced phase transition in both symmetry classes is of the quantum Lifshitz-type, uncovering a hitherto unexplored link between loop models in (2+1) dimensions, monitored free fermions, and (2+1)-dimensional CFTs.

The paper is organized as follows: Below we provide an overview of the main results for the universal long-distance behavior of both circuits from the loop model and the Clifford simulations. In Sec. <ref>, similarities and differences between Hamiltonian and measurement dynamics are discussed, introducing the Kitaev and Yao-Kivelson models both in the Hamiltonian as well as in the measurement-only circuit framework. In Sec. <ref>, the framework of Majorana loop models is introduced for class BDI and class D and for Hamiltonians and circuits in arbitrary dimensions. In addition, we provide a dictionary between the loop framework and circuit observables. The analytical and numerical results for the loop model framework are presented in Sec. <ref>. These are complemented by a selection of observables obtained in high-performance Clifford simulations in Sec. <ref>. A discussion and conclusion is presented in Sec. <ref>.







 §.§ Overview of main results


The measurement-induced dynamics of monitored circuits (discussed in this manuscript) share several similarities with the energetics-induced dynamics of their Hamiltonian counterparts.
In both cases, frustration is induced by the presence of non-commuting operators and the statistical mechanics of the stationary state can be mapped to a suitable (2+1)-dimensional loop model. While, however, the physics of Kitaev and Yao-Kivelson Hamiltonians is typically understood in terms of their Majorana band structure, their circuit counterparts feature a different link to many-body quantum physics. As a guiding picture, the Majorana loop model representation of both the measurement-only Kitaev and the Yao-Kivelson model in d spatial dimensions anticipate a link between the measurement-only circuits and the physics of diffusive Majorana metals – highly entangled liquid phases of Majorana fermions undergoing Brownian motion on a lattice – and disorder-induced localization in (d+1) dimensions. The transition between the liquid and the localized phase is induced by the degree of measurement-induced frustration, which we will use as an organizing principle in the following.




  §.§.§ Large measurement-induced frustration



When all allowed bond operators of nearest neighbor qubits (or free Majorana fermions) in the honeycomb or decorated honeycomb lattice are measured with equal probability, i.e., at the isotropic point in the phase diagram, the measurement-induced frustration is maximal. This pushes the wave function into an entangled state with the subsystem entanglement of an A=L× L square growing as 
	
    S(A)∼ Llog(L) ,
 
characteristic of quantum liquids with a nodal (Fermi) surface in two dimensions <cit.>. The liquid state is robust against small variations in the circuit parameters, i.e., away from the point of maximum frustration, defining an extended highly entangled phase at large measurement frustration, as illustrated in the phase diagrams of Fig. <ref>.

For measurements drawn randomly in space-time, this liquid phase is reminiscent of diffusive metals arising in weakly disordered Fermi systems. This picture is strengthened in the loop model framework: large measurement-frustration causes the endpoints of Majorana loops to perform an isotropic random walk in (d+1)-dimensional space-time. This yields an entanglement structure akin to a disordered metal with diffusion constant 𝒟, see Sec. <ref>.




The long-distance properties of the random walkers are governed by a non-linear sigma model in class ℂℙ^n-1 for symmetry class BDI or ℝℙ^n-1 for symmetry class D, see Eq. (<ref>) below. The replica limit is determined by the fugacity n of the underlying loop model, i.e., n=1 for the circuit. 
Both are single-parameter theories and depend only on the diffusion constant 𝒟, which consequently determines the behavior of physical properties at large distances. For instance, the entanglement entropy in the Majorana liquid phase is 

    S(A)∼√(𝒟)· Llog(L) .
 
The effective diffusion constant 𝒟 at large distances emerges from the under renormalization group transformations. Deep in the liquid phase, the renormalization of 𝒟 is negligible and we find reasonable agreement when approximating 𝒟=𝒟_mic, where 𝒟_mic is the microscopic diffusion constant from the measurement probabilities, see Fig. <ref>(a). 

We show that the physical properties in the circuit, such as, e.g., entanglement, mutual information, and correlation functions, can be inferred from the distribution function P(ℓ) of Majorana loops of length ℓ in (d+1)-dimensional space-time. For a two-dimensional system with linear dimensions L_x, L_y in the x- and y-direction, this distribution separates into two parts: a Brownian regime for ℓ<L_xL_y, and a Poisson-Dirichlet tail at space-time volume filling distances ℓ>L_xL_y, see Fig. <ref>. The Brownian regime features a loop distribution P(ℓ)∼ℓ^-5/2, which does not distinguish between symmetry class BDI and D. In contrast, the Poisson-Dirichlet tail differentiates between both symmetry classes: it displays an asymptotic behavior P(ℓ) ∼ℓ^-1(ℒ-ℓ)^θ-1, where ℒ is the maximum distance of a single loop in space-time and θ=1 for class BDI and θ=12 for class D, respectively. Both symmetry classes thus realize space-time filling loops in the Majorana liquid phase. The tails of the distribution are not detectable at the space time boundary and therefore do not discriminate different wave functions at final times. Instead, they reflect different dynamics in the circuit, which can be detected by out-of-time-ordered correlators (OTOCs), observable via circuit-ancilla measurement schemes, see Fig. <ref>. Such measurements allow to universally discriminate the dynamics of Majorana liquids in different symmetry classes. 




  §.§.§ Small measurement-induced frustration


For small measurement-induced frustration, i.e., when a large fraction of measurements can be satisfied simultaneously, the Majorana world lines rarely undergo a measurement-induced motion. Instead, they remain localized in space with an average localization length ξ, akin to disorder-induced localization of a Fermi liquid. The loop distribution P(ℓ) thus acquires an exponential cutoff, P(ℓ)∼ℓ^-5/2exp(-ℓ/ξ), suppressing the buildup of correlations and entanglement over distances larger than ξ. Consequently, the subsystem entanglement of a contractible region L× L obeys an area law 
    S(A)∼ Llogξ - ln 2 + ⋯
 on distances L>ξ. 
The negative subleading term is a universal topological correction contributed by the gauge flux, as the complementary part of the fractionalized degrees of freedom from the qubits. Such a topological state is in the same universality class as the toric code wavefunction <cit.>, and the circuit can potentially serve as a (dynamical) quantum error correction code – when the frustration is tuned towards zero, the state resembles an instantaneous state of the Hastings-Haah Floquet code <cit.>. 





  §.§.§ Critical measurement-induced frustration


The Majorana liquid and the area-law phase are separated by a phase transition at a line of critical measurement-frustration. In the Majorana loop framework, this is a localization-delocalization transition, which takes place when the diffusion constant in the non-linear sigma model renormalizes to zero. Precisely at the transition, the Majorana fermion entanglement obeys an area law but with a peculiar subsystem-dependence that reveals quantum Lifshitz criticality. Despite obeying an area law, correlation functions and the loop distribution function are not bound by a finite correlation length but display algebraic, scale invariant decay, which is faster than the conventional ∼ 1/(distance)^2. 

The fact that the is a single-parameter theory, with the only microscopic parameter entering being the microscopic diffusion constant 𝒟_mic, allows us to identify the location of the critical line of the phase transition. It is located at 𝒟_mic=𝒟_mic,c for some critical value 𝒟_mic,c.
In symmetry class BDI, we can extract the value of 𝒟_mic,c=3/16 exactly from the high symmetry points of the phase diagram (the yellow circles in Fig. <ref>(a)). This enables us to analytically determine, without any free parameter, the phase boundary in the Kitaev honeycomb geometry, see Fig. <ref>(b), and, tracing the same value of 𝒟_mic,c, also allows us to determine the phase boundary of the  Kekulé Kitaev model, see Fig. <ref>(c).
This connection to the universal value of 𝒟_mic,c in these two BDI models not only illustrates the universality of the loop framework, but it also explains the previously obtained 
circular structure of the phase diagram for the honeycomb geometry <cit.>.


The universal critical behavior at long wavelengths is determined by the ℂℙ^n-1 (ℝℙ^n-1) non-linear sigma model for symmetry class BDI (D) with the replica limit n→1 in d+1 dimensions. The corresponding universality classes have been discussed in the context of dense (fully-packed) polymers <cit.>, localization in gapless spin-singlet superconductors with disorder <cit.>, quantum magnets and spin-ice systems <cit.>.  We numerically extract four different universal exponents using both the loop model framework with simulations of up to 10^8 qubits and, using the stabilizer framework, for Clifford circuits with up to 10^4 qubits. These results are summarized in Tab. <ref>, where we also compare our estimates with previous loop model simulations of both symmetry classes (which reached system sizes equivalent to 10^6 qubits <cit.>). 


The critical exponents and their interpretation in the loop model and in the quantum circuit can be summarized as:

    
  * The correlation length exponent ν determines the average loop length ⟨ℓ⟩ 
    and the spanning number between the final and the initial state in the loop framework. 
    In the circuit, the former sets the average entanglement entropy in the localized phase 
    and the latter controls the late-time residual entropy during dynamical purification, see Fig. <ref>.
    
  * The anomalous dimension η describes the power-law scaling of loop correlations in the bulk of space-time. 
    	In particular, for the two-point watermelon correlator we have
	
    G_2(r) ∼{[       exp(-r/ξ) localized phase;        r^-(1+η)  critical point;          const.    liquid phase ]. .

          In the circuit this corresponds to the scaling of a two-point OTOC, see Fig. <ref>. 
          Alternatively, the anomalous dimension can be extracted through the hyperscaling relation τ = (11-η)/(5-η) 
          through the Fisher exponent τ. 
          The latter determines the loop length in (d+1)-dimensional space time  P(ℓ) ∼ℓ^-τ, see Figs. <ref> and 
          <ref>.
    
  * The fractal dimension d_f, which is related to η,τ by hyperscaling relations d_f = (5-η)/2=3/(τ-1), 
    describes the structure of the loops in space-time. As with η, these correlations are reflected in circuit OTOCs in the bulk of space-time. 
    See Fig. <ref> for the scaling at the transition and Figs. <ref>, <ref> 
    for scaling in the critical phase.
    
  * The order parameter exponent β determines the probability that a given space-time point falls on a macroscopic loop. 
    Analogous to quantum magnets, extensive loops correspond to macroscopic ordering of the spins, e.g., in a ferromagnet. 
    This is determined by the two-leg watermelon correlator in the limit of large separations, G_2(r ≫ 1, Δ) ∝Δ^2β, 
    where Δ is the distance to the critical point in the liquid phase, see Fig. <ref>. 
    This determines circuit OTOCs at long space-time distances.


Together with the universal OTOC ratios in the Majorana liquid phases, this provides a complete characterization of the two symmetry classes of measurement-only circuits.




§ KITAEV MODELS, CIRCUITS, AND SYMMETRIES




To set the stage we briefly review the fundamental traits of Kitaev physics – the fractionalization of quantum mechanical degrees of freedom and the subsequent formation of long-range entanglement <cit.>. We do this first in the context of the well-known honeycomb Kitaev model, i.e. in the Hamiltonian model in the parlance of this manuscript, where the ground states are quantum spin liquids with a characteristic and defining entanglement structure. 
In a second step we turn to the quantum circuit variant of interest in this manuscript. We then turn to an elementary symmetry classification of these systems and the conclusions one can draw from this with regard to the formation of long-range entanglement.







 §.§ Kitaev physics






  §.§.§ Hamiltonian model


The Kitaev Hamiltonian can be defined on any lattice geometry that allows for a tricolorization of its bonds. This is naturally the case for trivalent lattice geometries such as the honeycomb lattice, but can also be expanded to more complicated lattice geometries with higher coordination number 
<cit.>. 
On each bond, an Ising-like interaction couples two neighboring spin-1/2 degrees of freedom along the X, Y, Z directions (in spin space) 
associated with the bond's color

    H = ∑_⟨ jk⟩_μ K_μσ_j^μσ_k^μ ,

where μ =x,y,z indicates the easy-axis of the bond-directional exchange.

For any tricoordinated lattice geometry, this model can be exactly solved via rewriting the model in a Majorana basis <cit.>.
To this end, each spin-1/2 is mapped onto four Majorana fermions {b^x, b^y, b^z, c} satisfying a gauge constraint b^x b^y b^z c = 1, see also the illustration in Fig. <ref> above and Fig. <ref> below. This new basis allows us to reexpress the Pauli operator as Majorana bilinears σ_j^μ=ib_j^μ c_j^μ. The bond-dependent Ising interaction is then equivalent to a gauged Majorana hopping bilinear σ_j^μσ_k^μ =  i u_kjc_jc_k, where u_kj≡ i b_k^μ b_j^μ=± 1 denotes the gauge connection on the bond.
Note that the gauge field u is not a gauge-invariant physical observable, but instead one needs to consider the Wilson loop plaquette operator, as a product of the Ising interactions surrounding a plaquette, which evaluates the gauge flux penetrating the plaquette. Importantly, these gauge fields remain entirely static and assume a fixed low-energy configuration (favored by a Majorana fermion mediated interaction between the gauge fluxes)
that can typically be predicted via Lieb's theorem <cit.>.
We will discuss this gauge physics in more detail below for the two principal lattice geometries of our manuscript.

The virtue of this reformulation, which at first sight might be perceived as simple operator algebra, is that it faithfully captures the low-energy physics
of the Hamiltonian model: The spin-1/2 degrees of freedom fractionalize into an itinerant Majorana fermion coupled to a static lattice gauge theory in the background –  key characteristics of quantum spin liquid ground states <cit.>. The precise nature of these spin liquids depend on the underlying lattice geometry and the relative strength of couplings, but generally include gapless, gapped, and chiral spin liquids which are accompanied by long-range entanglement structures <cit.>. In fact, it is this formation of long-range entanglement in a quantum magnet that is the defining characteristic of a quantum spin liquid state and allows to unambiguously distinguish them from ordered ground states or finite-temperature paramagnetic states <cit.>.





  §.§.§ Quantum circuits


In a new twist on Kitaev physics, several groups have recently started to look into monitored Kitaev models <cit.> or “Kitaev circuits". 
In these circuits one replaces the bond-directional exchange terms of the Hamiltonian by bond-directional joint two-qubit measurements, also known as parity checks. These are illustrated in Fig. <ref>, along with their circuit implementation using an ancilla qubit.
On a theoretical level, these bond-directional parity checks implement a non-unitary Kraus operator

    M_jk≡exp(-τ/2s^μ_jkσ_j^μσ_k^μ) / √(2 cosh(τ)) ,

where μ = x, y, z indicates the type of bond-directional parity check, s=± 1 is the measurement outcome of the ancilla qubit, 
and τ∈[0,+∞) characterizes the strength of measurement that is controlled by a unitary entangling gate parameter t∈[0,π/4]: 
tanh(τ/2)=tan(t). 
The entangling CNOT gates of Fig. <ref> correspond to t = π/4 rotations, resulting in strong, projective measurements.
Here we focus on this Clifford regime of projective measurements and randomly sample the two-qubit parity checks in every microstep without evoking any spatio-temporal ordering (see Fig. <ref> above).
On a global scale one thereby transitions from (continuous) Hamiltonian dynamics where each exchange term is a projector onto a single, energetically favored two-spin state to a (discrete) measurement-only dynamics <cit.> where each projective measurement results in two possible, random outcomes.









While this might sound like a somewhat obscure idea, one should be reminded that such a connection between Hamiltonian and measurement-only circuit is also deeply embedded in the concept of the toric code <cit.> – it points to the implementation of a topological quantum memory via 
stabilizer measurements, i.e. a measurement-only quantum circuit that implements rounds of four-qubit measurements 
(whose outcome is interpreted as a syndrome and fed into a decoder to perform quantum error correction). 
Of course, the toric code is also very well known in its Hamiltonian version which serves as an exactly solvable spin model that gives rise to a topological spin liquid ground state. 
A key distinction between the toric code and the monitored Kitaev models of interest here is the commutativity of operators. 
While all Hamiltonian terms/measurement operators commute with one another for the former, they do not for the latter. 
This elevated level of measurement-induced frustration has important consequences on the entanglement structures 
that can be stabilized in these circuits as we have discussed in the introduction and which, in contrast to the toric code, 
are also distinct from the entanglement in the ground states of the corresponding Hamiltonians. 

Returning to the monitored Kitaev models at hand, it is important to note that running such a Kitaev circuit on a given initial state 
will give rise to a multitude of possible trajectories, depending on the sequence of measurement outcomes. 
It is the ensuing ensemble of pure states,
created by running a Kitaev circuit multiple times, that is the object of interest in the following. It allows to
define (in a statistical sense) typical states with a characteristic entanglement structure which we wish to classify
and compare to the ground-state entanglement stabilized by Hamiltonian dynamics 
(see Appendix <ref> for some Hamiltonian phase diagrams).

To understand the formation of long-range entanglement in the circuit setting, one can again rely on the language of fractionalization
which we have introduced above in the context of the exact solution of the Hamiltonian model in a Majorana basis. To do so, let us
perform a completely analogous rewriting of the 2-level qubit states in terms of four Majorana fermions and a subsequent regrouping into one
c-Majorana fermion per site and gauge connections (arising from Majorana bilinears on every bond) as illustrated in Fig. <ref> below.
In this language, every two-qubit parity check has the effect of measuring the local Majorana fermion parity, 
while effectively teleporting single Majorana fermions <cit.> – for a schematic illustration see panels (c) and (d) in Fig. <ref>. Repeated local measurements can then quickly lead to teleportation over large distances and the formation
of long-range entanglement, see Fig. <ref>(e).




 §.§ Time-reversal symmetry and gauge physics



A crucial ingredient in the classification of entanglement emerging either in the ground states of Kitaev Hamiltonians or in typical states stabilized by Kitaev circuits is the role of time-reversal symmetry and gauge physics. 

To set the stage, let us recall that, for a bipartite lattice/circuit geometry, the time reversal transformation can be defined as 𝒯=iσ^y 𝒦 where 𝒦 is the complex conjugate, which leads to σ^x(y)(z)→ -σ^x(y)(z). The combination of time-reversal symmetry 𝒯^2=+1 with the intrinsic particle-hole symmetry of Majorana fermions puts both the Hamiltonian model and the quantum circuit in the BDI symmetry class.

The situation is fundamentally different when the symmetry class is changed from BDI to class D, which can be realized by moving to a non-bipartite geometry as first noted by Kitaev <cit.> and later explored in the context of a lattice model by Yao and Kivelson <cit.>. 
Even though on the microscopic level of the Hamiltonian or quantum circuit description the system remains time-reversal invariant, 
the low-temperature ground state or dynamically stabilized typical state is not.
Instead, time-reversal symmetry is spontaneously broken in this non-bipartite setting as we will discuss in the following.




  §.§.§ Hamiltonian model



Let us first inspect time-reversal symmetry in the context of Hamiltonian models, where its role in stabilizing ground-state flux patterns 
and its spontaneous breaking for non-bipartite lattice geometry is well established.
As a first example, consider the seminal honeycomb model with its bipartite lattice geometry. 
Here the Wilson loop, indicating the flux through one elementary hexagonal plaquette, can be expressed as 
W_p = σ_1^xσ_2^yσ_3^zσ_4^xσ_5^yσ_6^z=± 1, 
with a clock-wise ordering of the six spin-1/2s around the plaquette. 
This flux is invariant under time reversal symmetry, i.e., 𝒯^-1 W_p 𝒯 = W_p. 
In addition, the Wilson loop operator commutes with the Hamiltonian, resulting in a static flux pattern, which by virtue of Lieb's theorem <cit.>
is the completely flux-free state in the ground state. As a consequence, the ensuing spin liquid is primarily characterized by the band
structure of the free Majorana fermions. For isotropic coupling, K_x = K_y = K_z in Eq. (<ref>), this is the celebrated
quasi-relativistic band structure with two gapless Dirac cones. 
In terms of entanglement, this is an area-law entangled state albeit one with a subleading O(ln L) correction <cit.>

    S(A) ∼ L + ln L + … .


Now, if we move to the Yao-Kivelson model on its non-bipartite lattice, note that while the flux through the large, dodecagonal plaquette (of length 12)
is invariant under time reversal, the triangle flux is not – W_ = σ_1^xσ_2^yσ_3^z changes sign under time reversal. In the low-temperature regime of the Hamiltonian (i.e. below the thermal crossover at which the spin degrees fractionalize <cit.>), the Majorana fermions mediate an effective Ising interaction between the triangle fluxes <cit.>. This leads to a flux-ordering transition at a finite transition temperature <cit.> at which the system spontaneously breaks time reversal symmetry. Simultaneously, the Majorana fermions respond to the chiral triangle flux and form a p+ip superconductor which falls into two topologically distinct phases <cit.>: a topological phase at weak pairing and a trivial phase at strong pairing. The transition between the two states can be driven by the relative strength of the triangle versus hexagonal coupling (J and K in the overview of Fig. <ref>).

In the language of quantum magnetism, the time-reversal symmetry breaking, topological state is called a chiral spin liquid <cit.>,
with an entanglement structure 

    S(A) ∼ L  - ln 2  ,

where the topological nature of this phase is reflected by a non-vanishing, subleading topological correction γ_ topo = ln(2) to the area-law entanglement <cit.>. 





  §.§.§ Quantum circuits



Let us now turn to a discussion of the role of time-reversal symmetry and gauge physics in the measurement-only circuit analogs of the Kitaev model. 

If starting from a featureless, maximally mixed state, the purification dynamics due to measurements can be separated into two contributions <cit.> – a gauge flux purification and a Majorana purification. 
Since the gauge flux commutes with the measurement operators, it is purified with a constant rate leading to an exponential decay of its entropy in time.
The Majorana purification, in contrast, depends on the degree of frustration – it is exponentially fast for the small frustration phase, but algebraically slow for the large frustration phase. Its behavior is tightly related to the pure state entanglement entropy that will be discussed in the following. 
To summarize the main results, the purification time for the gauge flux always scales logarithmically with the system size (ln N), while the purification time of the Majorana fermion takes (ln N) in the less frustrated and less entangled phase, but (N^p) power-law scaling with the system size in the highly frustrated and entangled liquid-like phase. 
Such two-stage purification dynamics was found to be a rather generic feature for measurement-only circuits built from subsystem codes <cit.>.






  
*Gauge purification.—

Let us first discuss the gauge flux purification dynamics in more detail. Note that the flux is not directly measured in our protocol but instead indirectly via the Majorana fermion parity checks. 
Every single Majorana parity measurement enforces a Majorana hopping (see Fig. <ref>), and when a set of measurements encircles a plaquette the moving Majorana detects and collapses the gauge flux to one of its eigenstates. 
As a consequence, a particular block in the block-diagonal matrix representing the Gaussian Majorana fermion state is picked up

    (
    ρ(W)     0     0    ⋯
    
    0    ρ(W')     0    ⋯
    
    0     0    ρ(W”)    ⋯
    ⋯   ⋯   ⋯   ⋯
    )  ,

where each block corresponds to a fixed gauge flux configuration W, W', W”. 
Over time, this block matrix is being purified towards a pure state, akin to lowering the temperature of a free fermion system. 
Due to the Clifford nature of the Kitaev circuits discussed here, the probability distribution for any possible flux configuration is equal, such that the state is purified into a random flux configuration, distinct from the flux-ordered Hamiltonian ground state <cit.>. Nevertheless, the states associated with any gauge flux configuration share exactly the same entanglement entropy. This is due to the generic property of Clifford stabilizer states: the entanglement of a Clifford stabilizer state ρ =  ∏_j=1^N1± O_j/2 where {O_j} are the stabilizer generators, does not depend on the eigenvalues being positive or negative. A minimal example is that all four possible Bell pairs |↑↑⟩+|↓↓⟩, |↑↓⟩+|↓↑⟩, |↑↑⟩-|↓↓⟩, |↑↓⟩-|↓↑⟩ share one bit of entanglement entropy. 

Since the Majorana parity checks are performed stochastically, there is a constant rate of purifying the gauge flux (denoted as r in the following), 
which leads to an exponential decay of the flux entropy in time  <cit.>

    S_flux(t) ∼ N_p e^-rt ,

where N_p is the number of the plaquettes of the same type. 
This decay is illustrated in Fig. <ref> 
for various points in the phase diagram of the monitored honeycomb Kitaev model – the initial decay of the state entropy towards a plateau at S(t)/N=12 signals the flux purification, which is independent of the circuit phase.
In general, the purification rate r depends on microscopic details, such as the number of edges of the plaquette 
– the larger the plaquettes, the slower the gauge purification. 
Apart from these microscopic details, the time scale to purify the flux to a pure configuration can be estimated as t ∝ln N_p / r
– the purification evolution here is akin to lowering the temperature to approach the ground state flux configuration in the Hamiltonian <cit.>. 

For plaquettes with an odd number of edges <cit.>  there is an obstruction to purification (that renders r=0) due to the time-reversal symmetry of the protocol, as proved in the following: Any stabilizer of the time-evolved state O is composed of a subset of the local quantum measurement operators K_j

    O = ∏_j K_j  .

Thus, when every local measurement respects a global symmetry such as time reversal, so does any stabilizer operator of the evolved state

    𝒯^-1 O 𝒯 =∏_j ( 𝒯^-1 K_j𝒯) = O  .

For example, a hexagon or dodecagon plaquette can be decomposed into its surrounding six bond operators, while a triangle plaquette operator is time reversal odd and cannot be generated by any two-body time-reversal-even operator measurements. Nevertheless, an even product of triangle plaquettes are time reversal even and can be purified, for example considering the sequence of stabilizer generators depicted in Fig. <ref> below.


This implies that a product of an even number of triangle plaquettes can still be purified exponentially fast, while a product of an odd number of triangle plaquettes remains always undetermined, respecting the global time-reversal symmetry. Thus, one global bit, corresponding to the action of a chiral Ising time-reversal transformation, remains conserved, leaving one bit of entropy, the same as the time reversal symmetric mixed state of the Yao-Kivelson Hamiltonian at low temperature. 

Crucially, for the quantum circuits, the entanglement dynamics does not depend on the individual flux assignments of the gauge configurations. 
Thus we have two strategies to break time-reversal symmetry: 
(i) we directly measure the 3-body next-nearest neighbor (n.n.n.) Kitaev interaction on the honeycomb lattice, 
which explicitly breaks the time reversal symmetry by fixing the ±π/2 fluxes on this triangle spanned by the n.n.n. spins; 
(ii) we move to the decorated honeycomb lattice with triangle plaquettes (a la Yao and Kivelson) and fix a single triangle flux from the initial state. 
Both strategies yield, as we will show in the following, the same entanglement phases. 




  
*Majorana purification.—

The purification dynamics of the Majorana fermions (the second stage of purification in Fig. <ref>), in contrast, primarily depends on the amount of frustration, given by the relative weight
of non-commuting operators, and not the specifics of the underlying lattice geometries.
For small frustration, i.e. a large bias towards certain commuting measurement patterns (such as the toric code stabilizers), 
the Majorana fermion purifies into the preferred dimer states exponentially fast, thereby stabilizing only a few, exponentially long-lived, 
global topological states (potentially serving as logical qubits). 
For large frustration, i.e. closer to the point of isotropic measurement strengths, in the liquid-like phases of the phase diagram, 
the Majorana fermion generally purifies in a power-law, governed by the physics of the statistical loop model to be introduced in Section <ref>. 




 §.§ Relation to Hastings-Haah code



This section has introduced and discussed Kitaev circuits as measurement-only quantum circuit analogs of the Kitaev model with an eye on the underlying lattice/circuit geometry. Before closing we also want to draw a few connections to the Hastings-Haah code <cit.>. 
The latter has attracted considerable interest from the quantum information community for its ability to dynamically stabilize a topological qubit,
while also being intimately connected to the Kitaev physics discussed above. 

The Hastings-Haah code imposes two-qubit parity checks with a strong spatio-temporal ordering. In space, it assigns a Kekulé pattern 
of couplings to the honeycomb geometry (see Fig. <ref>) and then performs time-periodic measurements of 
all XX, YY, ZZ parity checks. This Floquet dynamics periodically swaps the instantaneously stabilized state among the three corners of the
Kekulé Kitaev circuit phase diagram (see Fig. <ref>), 
which like their conventional honeycomb counterpart, are toric code phases. This is what allows it to create a topological qubit state. 

Here we take the broad interest in the Hastings-Haah code as motivation to study the measurement-only circuit variant of the 
Kekulé Kitaev model <cit.> along with its conventional honeycomb siblings.







§ MAJORANA LOOP MODELS IN  DIMENSIONS




Loop models arise quite generically in systems of quadratic Majorana fermions <cit.>. 
Here we describe how the loop picture emerges, both in Hamiltonian dynamics and in quantum circuits.
We briefly review the statistical mechanics of loops in 3D and then discuss the universality classes of loop models, 
based on symmetry and loop fugacity, and name representative models therein.
Finally, we define important observables and entanglement quantities for loops and Majorana circuits, providing a dictionary between the two.





In a quadratic Majorana model all information about the state is encoded in the two-point Majorana correlation functions ⟨ i γ_l γ_m ⟩. 
If one asks, in addition, that the quadratic Majorana state is also a stabilizer state, then all Majoranas have definite pairings ⟨ i γ_l γ_m ⟩ = 0, ± 1. 
Such a stabilizer state can be graphically represented by a pairing diagram wherein nodes l and m are connected by an open arc if and only if 
i⟨γ_l γ_m ⟩ = 1. 
Any generic Gaussian state can then be written as a linear combination of stabilizer states or, alternatively, of pairing diagrams. 
Evolving the state over time leads to rearrangement of the Majorana pairings, with fermion world lines tracing out loops in space-time as schematically illustrated in Fig. <ref>.
The entanglement and mutual information between two distinct subregions in space A,B are given by the number of arcs connecting A and B <cit.> (e.g., the blue loop in Fig. <ref>),
independent of the sign of the Majorana parities [Diagrammatic representations including the parity sign are possible but follow a more complicated loop algebra than required here <cit.>.].
When the number of Majoranas is conserved at all times (e.g., as for a spin-chain after a Jordan-Wigner transformation), the space-time is fully packed with loops.

Generic models of fully-packed loops may be described by a braid-monoid algebra <cit.>.
Loop configurations and transformations between them are then generated by a sequence of local operations (i.e., generators of the algebra) which rearrange the world lines.
The allowed local operations acting on Majoranas γ_l and γ_m are depicted in Fig. <ref>(a).
The identity 1_lm leaves the pairings unmodified, amounting to uninterrupted propagation forward in time.
By contrast, the Temperley-Lieb (TL) generators e_lm result in spatial propagation of loops.
On the Majoranas, this generator acts as a projector onto the local fermion parity, 

    e_lm∝_lm = 12(1 + iγ_lγ_m)  .

In particular, for initial pairings (k,l) and (m,n), e_lm implements a loop surgery to yield new pairings (l,m) and (k,n).
Finally, the braid operator b_lm acts as the Majorana swap operator 

    _lm = 1√(2)(1 + γ_l γ_m)  ,

exchanging the two Majoranas and causing the world lines to cross over one another.

Loop models corresponding to a braid-monoid algebra have three key properties – ambient isotopy, scalar fugacity (d-isotopy), and twist relations – which are depicted in Fig. <ref>(b-d).
Ambient isotopy reflects the ability to smoothly deform loops.
The scalar fugacity endows closed loops with a weight n.
Lastly, the twist relations allow braids b_lm (b_lm^-1) to be undone by incurring a factor ω (ω^-1) related to spin statistics.
All three leave the connectivity of the loop endpoints unchanged and may thus be viewed as equivalence relations (up to a scalar) between Majorana trajectories with identical entanglement in the initial and final states.





Throughout this work, we are interested exclusively in the dynamics of entanglement and mutual information in Majorana models, such that only the connectivity ⟨ iγ_lγ_m⟩ matters.
Braiding phases can thus be neglected and we fix ω = 1.
This yields the equivalence relation b_lm∼ b_lm^-1 on the space of loop configurations, reducing the braid-monoid algebra to the Brauer algebra (n) <cit.>.

By decomposing the Majorana dynamics into a sequence of operations {1_lm, e_lm, b_lm}, one obtains a transfer matrix representation of the evolution, i.e. a brickwall quantum circuit as illustrated on the left-hand side of Fig. <ref>. 
We may now distinguish two principal scenarios where Majorana loop models arise: Hamiltonian and circuit dynamics.








  §.§.§ Hamiltonian dynamics

Suppose that the local transfer matrix takes the form 

    T ∼1 + τ_lm .

Such a transfer matrix arises naturally from Trotterizing the imaginary time evolution e^-β H of a quadratic Majorana Hamiltonian H = i∑_⟨ l m ⟩ A_lmγ_l γ_m.
In order to satisfy the loop isotopy condition e_lm e_mn e_lm = e_lm, we take normalized projectors e_lm = √(2)_lm akin to defining the generators of the TL algebra <cit.>.
This implies the relation e_lm^2 = √(2) e_lm and identifies the loop fugacity 

    n = √(2)

for the imaginary time evolution of a Gaussian Majorana Hamiltonian.
The details of the Majorana adjacency matrix A_lm≠ 0 fix the geometry and the symmetries for the model.



  §.§.§ Circuit dynamics

Alternatively, the transfer matrix may arise in a measurement-only circuit consisting of projective measurements of parities iγ_lγ_m.
Here the sign of the measurement outcome poses a challenge to specifying a loop representation.
In particular, the two possible outcomes correspond to _lm and _ml.
Distinguishing such phase information in the graphical language of loops would require non-trivial twisting relations (ω≠ 1), thereby complicating the relative weights in the trajectory-averaged loop ensemble.
If, however, we are interested only in entanglement and other quantities which depend exclusively on the modulus ⟨ i γ_l γ_m ⟩, then the sign of the measurement outcome may be neglected.
This amounts to imposing a twist relation ω = 1 and thus an equivalence relation _lm∼_ml, with a loop configuration now representing an equivalence class of circuits.
Unlike the Trotterized imaginary-time evolution from before, here the projector _lm is always accompanied by a normalization of the state (i.e. |ψ⟩→_lm√(⟨_lm⟩)|ψ⟩) yielding a transfer matrix of the form

    T ∼1 + p1-p_lm .

Due to the normalization, projective measurements automatically satisfy the loop isotopy condition.
We may thus identify the generator e_lm as the operator which acts on any state |ψ⟩ as the normalized projective measurement _lm.
In addition, the idempotence of projectors _lm^2 = _lm fixes the loop fugacity 

    n=1

for each individual measurement in the circuit.




  §.§.§ Universality classes of loop models



The behavior of correlation functions and observables at large distances in loop models is determined by two ingredients: the loop fugacity n and the underlying symmetry of the loop model. Despite realizing a different loop fugacity, the transfer matrix of both the Hamiltonian and measurement-only circuit are generated by the operators 𝒫_lm, which are quadratic in Majorana fermions. Thus each Hamiltonian symmetry class  has a corresponding symmetry class counterpart in the measurement-only circuit. The loop fugacity then yields a fine structure of universality classes for each symmetry class.

For quadratic Majorana theories, particle-hole (PH) symmetry is always present and we thus distinguish two different symmetry classes based on whether time-reversal symmetry is present or not. In class BDI time-reversal is present, which enables a bipartition of the Majorana lattice into sublattices A and B, such that each operator 𝒫_lm acts exactly on one Majorana fermion on sublattice A and on one Majorana fermion on sublattice B. The paradigmatic example of such a model is Kitaev's honeycomb model in two spatial dimensions <cit.>. In the loop model framework, this symmetry is known as “orientability”. It allows to assign a unique orientation to each lattice site, either forward or backward in time, such that each loop carries this orientation through space-time.

In the absence of time-reversal, no bipartition of the lattice can be found, translating to the absence of orientable loop configurations. In the loop framework, the latter is known as non-orientable loops models or completely packed loop models with crossings (CPLC). This is symmetry class D and paradigmatic Hamiltonian models include the Yao-Kivelson model <cit.>, the Kitaev honeycomb model with next-nearest neighbor interactions as well as in the presence of a magnetic field. In the following, we will synonymously apply the language of Hamiltonian symmetry classes and loop model symmetry classes, i.e., we will refer to symmetry class BDI for orientable loop models and symmetry class D for non-orientable ones. 




 §.§ Statistical mechanics of loops in 3D



A d-dimensional Majorana circuit (d=2 in this work) corresponds to a d+1-dimensional loop model in space-time. 
Each lattice site l with Majorana fermion γ_l at a given time t when measurements are performed represents a vertex in a regular space-time lattice.
When all Majorana fermions are included in at least one measurement 𝒫_lm during the circuit evolution, the loop model is fully packed, meaning that every vertex in the lattice has coordination number z=4.
Then the space of possible loop configurations is generated by routing loops through each vertex in one of the three ways depicted in Fig. <ref>(a).

Such a model is described by a partition function 

    Z = ∑_𝒞 W(𝒞) n^N(𝒞) ,
 
where the sum is over all possible loop configurations 𝒞.
Each term in the sum involves two components: (i) a “local” part W(𝒞) and (ii) a “non-local” part n^N(𝒞).
The local term W(𝒞) is a product of Boltzmann weights associated to the choice of loop connections at each vertex, giving the probability of performing the unique series of measurements (or Hamiltonian evolution steps) that yields the loop configuration 𝒞.
This is the product of local measurement probabilities (or Hamiltonian matrix elements A_lm). 
The non-local term n^N(𝒞) accounts for the loop fugacity n associated to each of the N(𝒞) closed loops in configuration 𝒞.

As a side remark, we note that the weight n^N(𝒞) poses a challenge to Monte Carlo simulations of loop models for any fugacity n≠1 since it cannot be represented by a local update rule for evolving a single time slice.
This singles out measurement-only circuits with n=1 as a particularly attractive example for simulating novel entanglement phases.
For larger integer fugacity n ∈ℤ_+, simulation is possible by more complicated Monte-Carlo methods involving the whole space-time lattice <cit.> or additional ancilla degrees of freedom <cit.>, but this is beyond the scope of this work.




  §.§.§ Non-linear sigma models


A field theory formulation of the loop model partition function Z is provided by Z=∫ D[Q]exp(-S[Q]) with the non-linear sigma model () action

    S[Q] = 1/2g∫^d x  [(∇ Q)^2 ] + (topological terms)

Here, Q takes values in for symmetry class BDI and for symmetry class D. 
The field Q is parametrized by a vector z, which is complex for class BDI and real for class D, such that Q^αβ = z^αz̅^β - δ_αβ is a traceless n × n Hermitian (orthogonal) matrix with normalization condition z^† z = n. 
At short distances, the coupling constant g scales with the microscopic diffusion constant 𝒟_mic∼ g^-2, which is derived in Sec. <ref>.
The measurement-only limit n=1 is obtained by taking n→ 1^+ as a replica limit or alternatively by considering a supersymmetric formulation <cit.>.

Both the and the non-linear sigma models appear in a variety of different applications. 
Prominent examples are disordered fermion systems and quantum Hall layers in two dimensions. 
Although the loop models considered here are in (2+1) dimensions, the relation to disordered fermions remains prominent: independent of the loop fugacity and the symmetry class, but depending on the parameter g and the topological terms, the support: (i) a long-loop phase with 𝒟≠0 at the largest distances, yielding long-range entangled Majorana pairs akin to a disordered metal for which the topological terms are irrelevant,
(ii) a short-loop phase with a diffusion constant 𝒟 flowing to zero under renormalization group transformations, yielding area-law entanglement akin to localized fermions states but with distinct topological properties [We note that this is different from the (1+1)-dimensional case, for which class BDI is always localized, except for at the critical points separating different topological phases, while class D may support a metallic state akin to weak anti-localization, known as Goldstone phase.].

While for n=2 the metallic phases of the action S[Q] in (2+1) dimensions are known from the Hamiltonian case, the n=1 limit describes loops undergoing Brownian motion with diffusion constant 𝒟∼ g^-2. For n=1, the intermediate distance behavior, i.e., on distances ℓ≤ L_x L_y of the spatial volume, of Brownian random walkers is the same for both symmetry classes. This, however, manifestly changes on larger distances ℓ≥ L_x L_y: here, the Brownian random walkers have non-zero probability of forming space-time volume filling loops which occupy a non-zero fraction of the space-time lattice. These macroscopic loops are sensitive to the symmetry class and yield a means to distinguish symmetry classes BDI and D in a measurement-only quantum circuit, see Sec. <ref>.

At the transition from the metallic to the localized regime, a phase transition occurs with universal scaling behavior depending on the symmetry class and the fugacity n <cit.>, which we explore below.




 §.§ Loop-Circuit Dictionary


Here, we connect the key observables of the loop model framework to the Majorana circuits, providing a mapping between the two pictures.
In passing from loops in (2+1)-dimensional space-time to the dynamics of Majorana fermions in two dimensions, one dimension in the loop picture is designated as time.
The state of the Majorana circuit at a given time corresponds to the surface of the loop model along one direction.
There are two such temporal surfaces: (i) that at t=0 where boundary conditions on the loops correspond to a choice of initial state for the Majoranas, and (ii) that at time t=T corresponding to the final state produced by the circuit.
Correlations in this output state amount to boundary correlations in the loop model.
By contrast, loop correlations in the space-time bulk require the ability to peek at the state midway through the circuit evolution.

To track information regarding both the surface and the bulk of the loop model, we keep track of not only the loop connectivity, but also their integrated path lengths.
Any open arc is labeled by a tuple (l, m, ℓ), where l and m represent the endpoints and ℓ is the path length. 
In a numerical simulation this information is readily incorporated into the update scheme:
Consider a state with initial pairings (k, l, ℓ_1) and (m, n, ℓ_2) where l and m are nearest neighbors.
Measuring the parity iγ_lγ_m yields a state with new pairings and loop lengths (k, n, ℓ_1 + ℓ_2 + 1) and (l, m, 1).
When k = m, such a measurement closes the loop with total length ℓ_1 + 1 which is then recorded in a histogram for the trajectory.




  §.§.§ Loop length distributions



The central observable for the loop models we consider here is the probability distribution for the length of loops, both along the (temporal) surface and through the space-time bulk. 
The surface loop length between two points l,m at fixed time T describes the distances between two entangled Majorana fermions γ_l,γ_m in the circuit. 
The distribution of surface loops thus yields complete information on the quantum state at time T with regard to its entanglement and mutual information. 
The bulk loop length yields information on out-of-time ordered correlators (OTOCs), explored below.

For an open arc connecting Majoranas γ_l and γ_m on the t=T boundary, we define the surface length ℓ_surf^α along the α-direction as follows.
Let ℓ_x x̂ + ℓ_y ŷ be the vector connecting sites l and m in the two-dimensional plane. 
Using periodic boundary conditions (PBC), the surface lengths of the loop are defined as ℓ_surf^α≡min(ℓ_α, L_α - ℓ_α).
We define the corresponding length distributions P_surf^α(ℓ) as that obtained in the stationary state of the circuit when taking a pure initial state with only local pairings.

Characteristics of the entanglement dynamics during the evolution are encoded in the bulk loop length distribution P_bulk(ℓ). 
It is the probability distribution for the total path length of a randomly selected loop.
Here we need to make a distinction between open arcs and closed loops, as these will generically exhibit distinct length distributions.
During the circuit evolution, many closed loops are formed, but the number of open arcs is fixed by boundary conditions and remains constant.
If we impose PBC in time, all loops are closed and there is no need to make this distinction.
In practice, temporal PBC can be implemented with the assistance of ancilla measurements, which is discussed in Sec. <ref> together with the subtle but important role of boundary conditions for the bulk loop statistics.





  §.§.§ Surface observables: Entanglement entropy and mutual information



Observables described by the loop distribution of a temporal surface at time t=T correspond to the entanglement content of the Majorana fermion state ρ(T). 
We consider the von Neumann entanglement entropy S_A ≡ -[ρ_A log_2 ρ_A ] for a subsystem A, where the reduced density matrix ρ_A ≡_A̅[ρ] is obtained by tracing out the complementary system A̅,
and the mutual information I_2(A,B) ≡ S_A + S_B - S_AB between subsystems A and B.
For Majorana stabilizer states, these measures of entanglement are determined wholly by the non-vanishing two-point Majorana correlations ⟨ iγ_l γ_m ⟩ = 1.
Thus the connectivity encoded by the loop representation contains all necessary information to specify S_A and I_2(A,B).
The mutual information I_2(A,B) therefore corresponds precisely to the number of loops connecting subsystems A and B, while the von Neumann entanglement entropy for a pure state is S_A = 12 I_2(A, A̅), reflecting that each Majorana arc carries a half qubit of entanglement.
Multipartite entanglement measures such as the tripartite mutual information I_3(A,B,C) vanish due to the additivity of mutual information in the loop framework.

Throughout this work we will frequently consider the entanglement entropy of a cylindrical subsystem of length ℓ, with boundaries encircling the system in the x̂ or ŷ directions.
For such a partition, the entanglement entropy can be computed from the surface loop length distribution <cit.>.
In particular, the average entanglement entropy for subsystem A in a cylindrical region of length ℓ  along the x̂ direction is given by

    S^x_ℓ = L_y ∑_ℓ'min(ℓ, ℓ') P_surf^x(ℓ')  .

The mutual information between such cylindrical subsystems can be computed in a similar manner.




  §.§.§ Spanning loops and purification



A single loop connecting the two temporal surfaces at t=0 and t=T correspond to Majorana fermions, which remain unpaired during the entire circuit evolution. Such loops are called spanning and reflect the dynamical purification under measurement: Consider a maximally mixed initial state ρ(t=0) ∝1, corresponding to each Majorana being unpaired at t=0. The purity of the Majorana state ρ at time t=T is given by the total entropy S_L(T) = -[ρ(T)log_2ρ(T)]. Each unpurified qubit at time T corresponds to two unpaired Majorana fermions and thus to half the number of spanning loops n_s(T), i.e., S_L(T) = 12 n_s(T).

In the loop model framework, the spanning number, i.e., the dynamical purification, serves as an accurate quantifier for the different loop model phases and the critical point separating them <cit.>.
To do so, we assume fixed aspect ratios T / L_y and L_x / L_y.
Then for small distance δ from the critical point, we employ a scaling ansatz,

    n_s(δ, L)     = f(x)[1 + βL^ ]  ,           

    x                  = L^1/ν δ[1 + αδ]  ,


where is the dimension of the leading irrelevant operator and α≠ 0 gives a polynomial correction.





  §.§.§ Bulk loops, watermelons, and OTOCs




Beyond the distribution of bulk loop lengths P_link(ℓ), bulk properties of the loop model are captured by the watermelon correlators G_k(r), the probability that two points in space-time separated by distance r are connected by exactly k distinct arcs.
For example, G_2(r) gives the probability that two points with separation r lie along the same loop.
When the points do not lie along the same fixed-time slice, then G_k(r) takes the form of an out-of-time order correlator (OTOC) for the Majoranas.
As was previously pointed out in Ref. Klocke_2023, watermelon correlators can be accessed in the circuit by an ancilla-based measurement scheme.
For example, the two-leg watermelon correlator G_2(r) can be measured by considering space-time points (t_1, l) and (t_2, m) separated by distance r.
Consider ancilla A and B, each consisting of a pair of definite parity, e.g. ⟨ i γ_A,1γ_A,2⟩ = 1.
After evolving the circuit to time t_1 (t_2), a “marked” loop can be effectively inserted into the circuit by coupling the ancilla and bulk via measurement of iγ_A,2γ_l (iγ_B,2γ_m).
The history of the world lines which had been passing through these points in space-time is then effectively stored in the remaining ancilla γ_A,1 and γ_B,1.
We then evolve the circuit to late time t_f ≫ t_2 and impose boundary conditions on the bulk qubits (e.g., by measurement) such that the only remaining open arcs are those through the ancilla.
Had we taken only individual Majoranas for each ancilla, then the watermelon correlator G_2(r) would explicitly take the form of an OTOC, ⟨γ_A(t_f) γ_A(0) γ_B(t_f) γ_B(0) ⟩ - ⟨γ_A(t_f) γ_A(0) ⟩⟨γ_B(t_f) γ_B(0) ⟩, which for stabilizer states coincides with the mutual information between the ancilla I_2(A,B).
By construction, the mutual information I_2(A, B) counts the number of loops connecting points (t_1, l) and (t_2, m) in space-time such that G_2(r) = 12 I_2(A, B).
For higher order correlators G_k(r), additional ancilla are required, but the general procedure remains the same.




§ DYNAMICS OF MAJORANA LOOPS



In this section we present results from large-scale simulation of the measurement-only circuits in the loop framework. As representatives for the symmetry classes, we consider the honeycomb Kitaev model for class BDI and the next-nearest neighbor honeycomb Kitaev model for class D. Numerical results for the Yao-Kivelson circuit are presented in App. <ref>.
By varying the relative probability for measuring different bonds, one may tune the circuit between distinct measurement-induced phases, including several topologically distinct area-law phases and an extended metallic phase with a logarithmic area-law violation and long space-time loops.
In symmetry class BDI, i.e., in the measurement-only Kitaev circuit, the area-law phases arise in the toric-code limit where one flavor of bond measurement dominates.
Area-law phases also appear in symmetry class D, i.e., in the next-nearest neighbor Kitaev geometry and the Yao-Kivelson circuit by appropriate choice of parameters, see Fig. <ref>.
Besides the various area-law phases, circuits in both symmetry classes host an entangled liquid phase, reminiscent of a weakly disordered Majorana metal phase, which we characterize in Sec. <ref>.
Then in Sec. <ref>, we shed light on the universality of the entanglement transition separating metallic and area-law phases. Finally, we discuss the shape of the phase diagram for class BDI and the properties of the metallic phase based on the picture of random Brownian walkers.




 §.§ Long loops and the Majorana liquid phase



When the non-commuting measurements introduce a sufficient degree of frustration, the endpoints of the Majorana world lines no longer remain confined to a finite area but start to undergo a random Brownian motion.
This generates long loops in the circuit and stabilizes a Majorana liquid or metallic phase with characteristic Llog L entanglement. 
The entanglement structure in the stationary state is determined by the surface loop distribution in the loop model. 
In the liquid phase, this distribution does not distinguish between the two symmetry classes BDI and D, which we discuss below. 
A universal distinction between BDI and D, however, is observable in the bulk loop distribution.








  §.§.§ Surface loops and entanglement



The entanglement structure of a Majorana state ρ(T) at time t=T is captured by the probability distribution P^x,y_surf(ℓ) of open loop arcs along the temporal boundary of the loop model.
In the liquid phase, the loop length distribution exhibits power-law scaling 

    P^α_surf(ℓ) ∼ℓ^-2 .

For a cylindrical subsystem of length L and circumference L_y, Eq. (<ref>) yields an entanglement entropy S(L, L_y) ∼ L_y log(L) (see Appendix <ref> for a numerical confirmation). 
Such a logarithmic violation of the area-law of fermions in 2D is often associated with the presence of a Fermi surface <cit.>, but it is also expected for metallic states with non-vanishing conductivity in the presence of weak disorder <cit.>. 

The scaling of P^α_surf(ℓ) and of the entanglement entropy is a general feature of the liquid phase for both class BDI and D.
Let us denote c̃_α for α = x,y as the coefficient for the entanglement S=c̃_α L_αlog(L) for a cut along the α-direction. Away from the critical point, this coefficient is generally not universal but rather varies continuously with microscopic details of the model (i.e. the lattice geometry and the bond weights). This may be seen analogously to a smooth deformation of the Fermi surface yielding a continuous modification of the logarithmic entanglement scaling in the Hamiltonian setting, an analogy, which 
we develop further below. 

The effect of orientability here is a rather trivial one.
Suppose that we fix the length ℓ_x.
In the orientable case, there are fewer allowed Majorana pairings for any ℓ_x, as correlations between sites on the same sublattice must necessarily vanish.
As a result, the typical total loop length ℓ = √(ℓ_x^2 + ℓ_y^2) is longer, and the probability density P_surf(ℓ_x) is reduced by a constant factor.





  §.§.§ Bulk loops and universal circuit correlations



The two symmetry classes BDI and D are distinguishable by their bulk loop statistics. On short distances ℓ≤ L_x L_y and away from the critical point, Brownian motion in three space-time dimensions yields a mean-field decay 

    P_bulk(ℓ) ∝ℓ^-τ

with τ=5/2 and a fractal dimension d_f=3/(τ-1)=2 for both symmetry classes,  depicted in Fig. <ref>(b). However, Brownian walkers in 3D have a non-zero probability of never returning to their starting point, and to form “macroscopic” loops which occupy a non-zero fraction of the space-time volume.
The statistics of such extensive loops provide an unambiguous means for distinguishing between the two symmetry classes.

For finite system sizes, macroscopic loops are sensitive to the spatiotemporal boundary conditions. 
Here, one distinguishes absorbing boundary conditions, yielding open Majorana world lines, such as, e.g. a mixed initial state, from reflecting boundary conditions in space-time [Strictly speaking, macroscopic loops would emerge also for periodic boundary conditions in time, which are, however, challenging to implement in a circuit.]. 
For the latter, all world lines are closed including at the final and initial time. 
This is realized by starting with a pure initial state at t=0 and by terminating the circuit at t=T with a fixed set of measurements that closes all world lines.
With reflecting boundaries all loops are closed, and macroscopic loops manifest in the distribution P_bulk(ℓ) giving rise to a distinct scaling regime at distances ℓ≥ L_x L_y, shown in Fig. <ref>(b).
Here, macroscopic loops follow a Poisson-Dirichlet (PD) distribution <cit.>
ℓ·P_bulk(ℓ) = θ/ℒ(1 - ℓ/fℒ)^θ- 1,   (ℓ≫L_x L_y)

where ℒ is the total number of links in space-time (i.e. the volume). Whereas f is a non-universal quantity, 
the parameter θ is a universal parameter which depends only on the loop fugacity n and the symmetry class.
In the most general case, one finds θ = n for symmetry class BDI, while for symmetry class D, θ = n / 2.

For the Majorana circuits at hand with fugacity n=1, this leads the normalized bulk loop length distribution to approach a constant 

    ℓ· P_bulk(ℓ)  ∼ℒ^-1   (orientable)

for ℓ→ℒ in class BDI, 
while class D displays a square-root divergence 

    ℓ· P_bulk(ℓ)  ∼ℒ^-1/2(ℒ-ℓ/f)^-1/2   (non-orientable) .

Such a divergence reflects the tendency for a concentration of probability density in larger macroscopic clusters, due to the greater mobility of loops without the BDI constraint of time-reversal.
The results of numerical simulations, shown in Fig. <ref>(b), confirm the PD scaling regime for macroscopic loops, with the parameter θ consistent with loop fugacity n=1.





In order to access the universal part of the PD distribution, i.e., the parameter θ, one faces two challenges: 
(i) implementing measurements of loop quantities in space-time and on distances ℓ≫ L_x L_y much larger than the spatial extent and
(ii) eliminating the dependence on the non-universal parameter f in Eq. (<ref>). 
Let us start by solving the latter.

Consider m points r⃗_i with i=1,…,m in space-time which are separated by macroscopic distances r⃗_i - r⃗_j≳ L_x,L_y.
In the thermodynamic limit, if any two points r⃗_i and r⃗_j are on the same loop, then that loop is macroscopic.
The probability that all m points lie along the same loop is <cit.>
P_m ≡P(r⃗_1, …, r⃗_m) = f^m Γ(1 + θ)Γ(m)/Γ(m + θ).
    
Here, Γ is the Γ-function and the non-universal quantity f appears as a prefactor.
Thus, f can be eliminated by taking appropriate ratios P_m_1^k_1 / P_m_2^k_2 for m_1 k_1 = m_2 k_2, leading to universal quantities which distinguish the two symmetry classes <cit.>.

We propose a way to measure the probabilities P_m directly by generalizing existing ancilla schemes <cit.>. 
Since P_2 is simply the large-separation limit r→∞ of the two-point watermelon correlator G_2(r), it may be measured via the protocol described earlier in Sec. <ref>.
For higher-order quantities P_m, we use m ancillae to insert marked loops at all m points in space-time, as depicted in Fig. <ref>.
Once again we evolve the circuit to late times and then measure all bulk Majoranas, leaving open loops only terminating on the ancillae.
The mutual information between groups of ancillae counts precisely the probability that the space-time points where we measured the bulk and ancilla had been along the same loop, so a quantity like P_3 amounts to requiring I_2(A,B) = I_2(B, C) = I_2(C,A) = 1.
In practice this can be obtained by comparing the magnitude of I_2(A,B) before and after tracing out ancilla C, where the latter case simply reduces to P_2.
Here we consider such a scheme with four points r⃗_i at separations on the order of L_y and measure P_2, P_3, and P_4 to obtain the following ratios,

    P_2^2/P_4    = Γ(1+θ) Γ(4+θ)/6Γ(2+θ)^2    = 1,    θ = 1 
    35/36,    θ=1/2,       
    P_2^3/P_3^2    = Γ(1+θ)Γ(3+θ)^2/4 Γ(2+θ)^3    = 9/8,    θ = 1 
    25/24,    θ=1/2.

Importantly, both give values which differ between the two symmetry classes.
In Fig. <ref>, we show that the ratios indeed converge toward the universal value as system size is increased.
This then offers a practical scheme for distinguishing between the long-loop phases which would otherwise not be possible from the steady-state entanglement.

We note that when implementing (open) absorbing boundaries in space-time, any macroscopic loop will be cut into an extensive number of individual open loops.  
Then the loop length distribution matches the first passage time distribution for a random walk in one spatial dimension with absorbing boundaries 

    P_bulk(ℓ) ∼ P_FP(ℓ) ∼ℓ^-5/2 e^-αℓ
 
with α^-1∼ L_x L_y.
In Fig. <ref> (see Appendix), we verify this expected scaling for open loops in the metallic phase of both symmetry classes.




  §.§.§ Analytical entanglement scaling and phase boundaries  from the diffusion picture



Viewing the dynamics in the metallic phase as Majorana world lines undergoing Brownian motion in space-time allows us to establish a direct link between (i) the non-linear sigma model in Eq. (<ref>), (ii) the entanglement entropy and, (iii) for class BDI, the microscopic measurement probabilities. 
Note that the latter enables accurate analytical predictions, for symmetry class BDI, for the entanglement scaling and the location of the transition line between the metallic and the localized phases for arbitrary geometries. 

We start by considering loop endpoints undergoing Brownian motion in space-time with an effective diffusion constant 𝒟 at large wavelengths. 
The entanglement entropy at a fixed time t=T is then determined from open arcs terminating at the temporal boundary. 
Designating one of the arc endpoints as the “start”, we may view the arc as the path of a random walker which starts near and eventually terminates upon an absorbing boundary at t=T.
Then the total bulk length of the arc is equivalent to the first passage time τ of the random walker
[More precisely, the movement of the loop endpoint along the temporal direction should be viewed as a persistent random walk, where the direction of propagation is reflected whenever the Majorana is involved in a measurement. This yields a mild correction to the timescales, but does not alter the universality.].
The distribution of first passage times for one absorbing boundary in three dimensions is known to be P_FP(τ) ≈√(2/π)τ^-3/2. 
During this time, both endpoints undergo Brownian motion in the two-dimensional spatial plane, yielding a distribution for the distances between endpoints P_2D(r,τ) ≈2r𝒟τe^-r^2/𝒟τ.
Averaging with respect to the time τ, yields the expected distribution of spatial displacements

    P_surf(r) = ∫τ  P_FP(τ) P_2D(r,τ) = √(2𝒟)/r^2

with r^2=ℓ_x^2+ℓ_y^2. 
The entanglement entropy along a cut depends on the distribution of loop lengths projected along a fixed axis.
For isotropic Brownian motion in space, this gives

    P_surf^x(ℓ_x) = 4∫_0^∞ℓ_y  P_surf(r)/2π r = 2√(2𝒟)/πℓ_x^2.

From this distribution we infer that the logarithmic entanglement scaling has coefficient c̃_x = 6√(2𝒟)π.
Varying the relative measurement probabilities of different bonds will generally introduce some anisotropy to the Brownian motion such that 𝒟^x ≠𝒟^y and c̃_α∝√(D^α). This scaling of the logarithmic entanglement proportional to the square root of the diffusion constant appears also in monitored fermion systems with unitary dynamics <cit.>, and is seemingly generic for measurement-induced fermion liquid-type phases. 

The effective diffusion constant 𝒟 emerges at long wavelengths from the renormalization group (RG) flow of the in Eq. (<ref>). The is a single-parameter theory, initialized with the microscopic short-distance diffusion constant 𝒟_mic, which depends on the lattice geometry. The RG flow only depends on the fugacity n, the symmetry class and the dimensionality d. For n=1 and d=3, each symmetry class thus has a one-to-one correspondence 𝒟_mic↔𝒟. For bipartite lattices, i.e., for class BDI, the random walker is always on the same sublattice after two steps and the microscopic diffusion constant is readily inferred from the measurement probabilities.
Consider the measurement-only Kitaev circuit, where bonds are measured with probabilities K_x + K_y + K_z = 1. The diffusion constants may be estimated by considering the mean-squared displacement of a random walker with jump rates set by K_α and bonds of unit length. Setting 𝒟^z_mic as the diffusion constant for the direction parallel to the ZZ-bonds and 𝒟^⊥_mic as the diffusion constant perpendicular to it yields

    𝒟_mic^z     = 9/8 K_z(K_x + K_y)  ,  𝒟_mic^⊥  = 1/6 (2𝒟_mic^z + 9K_x K_y)  ,                                                       
    𝒟_mic   = 1/2(𝒟_mic^z + 𝒟_mic^⊥) =3/4 (K_x K_y + K_x K_z + K_yK_z)   .

The total diffusion constant 𝒟_mic is invariant under continuous rotation of K_x/y/z with respect to the isotropic point 
– this explains the circular symmetry of both the metallic phase and its phase boundary, as numerically found in Refs. <cit.>.

We emphasize two particular regimes: (i) In the vicinity of the isotropic point K_α=1/3, i.e., deep in the metallic phase, we expect the RG flow corrections to 𝒟 to be weak. Thus setting 𝒟=𝒟_mic proves to be a good approximation in the vicinity of the isotropic point, see Fig. <ref>(a). Since under RG, the diffusion constant will generally flow towards a smaller value, 𝒟_mic at the isotropic point serves as an upper bound for the average entanglement growth in each direction 𝒟≤ 1/4. (ii) In the vicinity of the critical line marking the transition between the metallic and the localized phase, the renormalization of 𝒟_mic is strong and 𝒟_mic≫𝒟. Nevertheless, the one-to-one correspondence 𝒟_mic↔𝒟 [Neglecting the effect of anisotropy in the microscopic values.] allows us to extract the position of the critical line for any bipartite geometry. At K_x = K_y = 12 and K_z = 0, the location of the critical point is exactly known: the system decouples into disconnected one-dimensional strings each of which displays percolation critical behavior <cit.>. Using this point as an estimate for the critical value 𝒟_mic,c provides

    critical diffusion constant: 𝒟_mic,c=3/16.

For the Kitaev honeycomb lattice, this yields a circular contour in parameter space (see Fig. <ref>(b)) which coincides exceptionally well with numerical estimates for the entanglement transition. 
To further underpin the strength of the diffusion picture, we provide the numerically obtained phase diagram for the measurement-only Kekulé-Kitaev geometry illustrated in Fig. <ref>(b). This system also belongs to symmetry class BDI but its phase diagram exhibits a different shape of the critical line as displayed in Fig. <ref>(c). Using Eq. (<ref>) again provides a rather accurate estimate for its phase boundary based on a calculation of the diffusion constant as
 𝒟_mic  = 3/8 [K_g K_r (2-3K_gK_r) + K_b(K_g + K_r)(2 + K_g K_r) - 3 K_b^2 (K_g^2 - K_g K_r + K_r^2) ], 
 where again we have considered the mean-squared displacement of a random walker after two steps along lattice edges, averaging over starting points in the unit cell.





Lastly, let us discuss a more rigorous approximation scheme for the entanglement entropy. We note that even deep in the metallic phase, the microscopic diffusion constant 𝒟_mic is larger than 𝒟 determined from the entanglement scaling by a factor of approximately 3/2.
This difference reflects the fact that the microscopic calculation treated the loop endpoints as non-interacting Brownian walkers.
In actuality, due to the Pauli principle, the loops are constrained to not overlap with themselves or one another, corresponding to an effective excluded volume.
To account for this, consider the role of a nonzero density ρ of closed loops in the bulk of space-time.
These act as obstructions to the Brownian motion of open arc endpoints, suppressing diffusion.
At leading order, this may be treated as a uniform probability ρ that at any step the random walker remains in place, sending 𝒟→𝒟(1-ρ). The density ρ is hard to extract analytically but easily accessible numerically. 
We find that the fraction ρ varies continuously from ρ≈1/3 at the isotropic point to ρ=1 in the area-law phase.
In Fig. <ref>(a), we show that the rescaled microscopic diffusion coefficient 𝒟_mic(1-ρ) better matches the fitted 𝒟, reproducing the magnitude deep in the metallic phase and going to zero continuously at the transition.
Thus, accounting for the excluded volume effect accurately captures the interactions neglected in 𝒟_mic.
Alternatively we might arrive at the same correction by considering a self-avoiding walk on a lattice with coordination number z and connectivity constant μ.
Here the quantity 1-ρ corresponds to the ratio μ/z, which falls in the range μ/z ∈ [1.2, 1.5] for various 3D lattices <cit.>.


The direct relation between the entanglement scaling, i.e., the prefactor c̃ of the metallic Llog(L) term, and the effective diffusion constant 𝒟 at distance L allows us to numerically track the RG flow of 𝒟 from its microscopic starting point 𝒟_mic towards its asymptotic value as illustrated in Fig. <ref>(a) for linear system sizes L=16,64,1024. This is remarkable, since in d=3 dimensions, the is no longer perturbatively controlled and the RG flow of 𝒟 is not accessible analytically. 
We emphasize the connection between c̃, the diffusion constant for Brownian motion of loops, and their relation to the logarithmic entanglement scaling found in weakly disordered Fermi liquids.
In the ground state of a disordered metal with disorder strength γ, the diffusion constant is related to the Fermi velocity v_F and the disorder by 𝒟∼ v_F^2 / γ, implying c̃∼ v_F / √(γ).
Indeed, such proportionality is expected from calculations of the entanglement entropy from a Fermi surface <cit.>. 




 §.§ Quantum Lifshitz criticality and universal scaling behavior at the entanglement transition




In both symmetry classes, the Majorana liquid and localized phases are separated by a critical line exhibiting quantum Lifshitz scaling of the entanglement entropy. At this critical line, the effective diffusion constant vanishes, resulting in an area law for the entanglement entropy. The distribution of Majorana loop lengths, however, is not yet localized but displays universal algebraic scaling behavior that clearly differentiates the two symmetry classes. This unique combination of universal algebraic scaling and the area law leads to the distinctive subleading entanglement scaling characteristic for quantum Lifshitz criticality.

In the Kitaev honeycomb model for symmetry class BDI, the critical line is reached by tuning the probabilities K_α. Quantum Lifshitz scaling is observed, except at the three points where one of the K_α = 0, reducing the model to an array of one-dimensional chains. For class D, we set the next-nearest neighbor bond measurement probabilities to J_x = J_y = J_z = J and tune J, K_α through the critical regions. For J > 0, quantum Lifshitz scaling is generally present.





  §.§.§ Quantum Lifshitz entanglement and surface loops




The entanglement entropy in the loop framework is computed from Eq. (<ref>). 
At a genuine two-dimensional transition we observe an algebraic decay of surface loop lengths, P_surf(ℓ) ∝ℓ^-γ_surf. 
Away from the percolation points with one K_α=0, we find that γ_surf≈3, shown in Fig. <ref>(a). 
This is faster than in the Majorana liquid phase, where γ_surf = 2 leads to an Llog(L)-growth.
As a result, the transition exhibits a subsystem dependence of the entanglement entropy S(ℓ, L) = aL + b J(ℓ / L). 
Here J(u) is the “Lifshitz” scaling function [For unit aspect ratio of the spatio-temporal volume.]

        J(u)    = log(θ_3(iλu)θ_3(iλ(1-u))/η(2iu)η(2i(1-u)))  ,
    
    
where θ_3 is the Jacobi theta function, η is the Dedekind eta function, and λ is a fitting parameter <cit.>.
In Fig. <ref>(b) we show the collapse of the entanglement data onto a curve accurately described by this Lifshitz scaling function.
For both symmetry classes, we find λ = 3.4 ± 0.1 for cuts along either direction of the system.

Although derived initially for a quantum Lifshitz model <cit.>, this torus entanglement scaling has been observed in a variety of (2+1)-dimensional CFTs <cit.>.
Notably, the Rényi entropy takes precisely the form of Eq. (<ref>) also for quantum dimer models and resonating valence bond (RVB) states <cit.>.
Given the close link between models of dimers, loops, and height fields <cit.>, one may expect that this entanglement scaling ought to extend also to the case of loops with fugacity n=1 (relevant to the quantum circuits at hand).
Indeed, higher fugacity loop models have been employed in numerical studies to verify the Lifshitz scaling of entanglement in quantum magnets <cit.>. 
Similarly, the bond-length distribution in the antiferromagnetic Heisenberg model on the square lattice exhibits the same scaling exponent as P_surf(ℓ) in our circuit model <cit.>.
We note here that quantum Lifshitz scaling is compatible with the simultaneous presence of conformal symmetry and an area law of the entanglement in (2+1) dimensions. 
It is remarkable that this entanglement scaling, though apparently quite universal, has been reported previously only in the context of interacting circuits <cit.>.
For fermion-parity symmetric circuits, the mapping to an interacting dimer model suggests that Lifshitz entanglement scaling is quite generic.
More broadly, requiring conformal invariance at the entanglement transition imposes strict constraints such that one might generically expect scaling of the form Eq. (<ref>) except at fine-tuned points.




  §.§.§ Correlation length exponent ν









To accurately determine the critical point K_c and the correlation length exponent ν, we consider the spanning number n_s for fixed aspect ratios L_x = 2L_y (i.e. L_x × L_x plaquettes) and circuit depth t=L_y.
We then perform a finite-size scaling (FSS) analysis using an ansatz generically of the form in Eq. (<ref>) and a cubic B-spline to fit the scaling function.
Moreover, by taking sufficiently large system sizes and a narrow window around the critical point, we find that any irrelevant scaling variables may be dropped from the ansatz.

In Fig. <ref> we show the resulting data collapse for the transition in (i) the (orientable) honeycomb Kitaev circuit along the line K_x = K_y = (1 - K_z)/2
and (ii) the (non-orientable) honeycomb Kitaev circuit with next-nearest neighbor parity checks along K_x = K_y = J = (1 - K_z)/3.

The results are summarized in Tab. <ref>. We find consistent values of ν in both symmetry classes for several other lattice geometries (see Appendix <ref>).
These results thereby highlight that the entanglement transitions in (2+1)-dimensional measurement-only Majorana circuits have distinct universality depending on the presence or absence of orientability, i.e., depending on the symmetry class. We fix the critical points K_c at the values reported here and use them in all further analysis.










  §.§.§ Fisher exponent τ



At the transition and throughout the extended critical phase, the bulk loop length distribution is expected to follow a power-law P_bulk(ℓ) ∝ℓ^-τ, where the Fisher exponent τ is related to the fractal dimension d_f via the hyperscaling relation τ = (d / d_f) + 1.
In the extended critical phase, the Brownian nature of the loops yields d_f = 2 and thus τ = 5/2.
By contrast, a mean-field analysis <cit.> predicts η = 0, corresponding to d_f = 5/2 and τ = 11/5 at the transition.
Tuning to the critical point determined from FSS of the spanning number, we obtain the exponent τ by fitting a power-law to the bulk loop length distribution P_bulk(ℓ). 
The result is provided in Tab. <ref>.

These results deviate very slightly from previous numerical studies on 3D classical loop models, where hyperscaling gave τ = 2.184 ± 0.003 for the orientable case <cit.> and τ=2.178 ± 0.002 for the non-orientable case.
In both symmetry classes, the decay of P_bulk(ℓ) is slower than expected for τ = 11/5, indicative of small but finite renormalization due to fluctuations beyond the mean-field prediction.
However, the difference in τ between the two symmetry classes is rather subtle, making it challenging to distinguish them even at the large system sizes accessible in the loop framework.






Finally, we note that a direct observation of the Fisher exponent τ might be practically infeasible as it would require access to the bulk lengths of loops in space-time.
As an alternative, it might be far more practical to employ ancilla probes <cit.> to instead extract the anomalous dimension η, from which we may determine τ via hyperscaling relations, see below.






  §.§.§ Anomalous dimension η





The anomalous dimension η can be determined by examining the finite-size scaling behavior of the total length of spanning loops ℳ,

    ℳ L^-(5-η)/2 = f(L^1/ν (K - K_c)).

Moreover, this provides a second, independent estimate of the correlation length exponent ν.
In Fig. <ref>, we show the scaling collapse of ℳ in both symmetry classes. The results are displayed in Tab. <ref>. We note that the estimates for ν obtained in this finite-size scaling analysis are independent of the values extracted from the spanning number collapse, Tab. <ref>.
Both approaches yield almost perfect agreement, which is a strong internal consistency check of our numerical approach.



An alternative way to obtain an estimate for the anomalous dimension η is to probe the two-leg watermelon correlator G_2(r).
At the transition, this correlator decays algebraically with the distance r, G_2(r) ∝ r^-(1 + η).
Analogously to the correlators P_n discussed in Fig. <ref> which we considered earlier, G_2(r) can be measured in an experiment by using an ancilla scheme: G_2(r) is the average mutual information between two ancillae that were entangled with the circuit at space-time points separated by a distance r.






Here, one would entangle ancillae at identical spatial positions but at times t_1, t_2 such that r = t_1 - t_2.
Then, for a fixed system size, the scaling of G_2(r) can be studied by simply varying the temporal separation (i.e. the circuit depth).
As shown in Fig. <ref>(a), we observe a power-law decay of the watermelon correlator with respect to the space-time separation that allows one to extract the anomalous dimension for both symmetry classes, displayed in Tab. <ref>.
We note that this approach is highly sensitive to subleading corrections and small deviations from the exact local of the critical point. We consider it thus generally less accurate to extract η in this manner.





  §.§.§ Order parameter exponent ß



In the loop framework, the Majorana liquid phase acquires a non-zero ferromagnetic-type order parameter, represented by the probability that two infinitely far separated space-time coordinates are connected by a single closed loop, i.e., P_2∼ |K-K_c|^2β. 
The critical exponent β associated to this order parameter can be detected in an ancilla scheme similar to Fig. <ref>. 
In practice, for systems of linear dimension L, one measures G_2(r) for large space-time separations r∼ L. 
This introduces a subleading correction G_2(r) = A r^-x_2 + P_2, where x_2 = 2 for Brownian loops <cit.>.
For even moderate system sizes L, this subleading correction becomes negligible.
We may then treat finite-size effects for G_2 similarly to those for the order parameter in percolation,

	G_2(r∼L, K)    = L^-2β/νf(K̃),

	K̃    = (K - K_c)L^1/ν[1 + A(K-K_c)].
	
    
In the liquid phase, this scaling function reduces to the simpler form G_2(r∼ L, K) ∝K_c - K^2β on distances L much larger than the correlation length ξ.


We determine β by fitting the ansatz in Eq. (<ref>) to the data in a narrow window enclosing the critical point.
In order to avoid systematic drift in the fitted exponents, we fix ν and K_c at the values found previously from examining the spanning number.
The resulting scaling collapse is shown in Fig. <ref>(b). 
This yields accurate estimates of the critical exponent β from the ancilla scheme which clearly distinguish both symmetry classes and which are summarized in Tab. <ref>.






  §.§.§ Hyperscaling and Discussion



Having obtained several independent estimates of various critical exponents, we now examine their consistency with one another via hyperscaling relations.
We further comment on the reliability of the different approaches for determining critical exponents and justify the choice of values reported in Table <ref>.

Given the Fisher exponent τ, the anomalous dimension is given by the scaling relation η = 5 - 6τ-1.
The resulting values are given in Tab. <ref>.
For class BDI, scaling gives a value of η which lies between the two estimates determined earlier in Tab. <ref>.
In symmetry class D, scaling suggests a slightly smaller but not inconsistent magnitude of η than previously found.
The anomalous dimension η is notoriously difficult to determine directly, even for the large system sizes accessible with our loop model simulations.
Computing η via the hyperscaling relation is highly susceptible to small deviations in the Fisher exponent τ and thus even smallest uncertainties in the critical point K_c.
By contrast, determining η directly via FSS of ℳ leaves K_c as a fitting parameter and is generally found to be more robust.
This is supported not only by the near perfect scaling collapse shown in Fig. <ref>, but also by the close agreement with the value of ν found by FSS of the spanning number.
As such we opted to report the corresponding value of η in Tab. <ref>.

Hyperscaling allows for two alternative ways of estimating the order parameter exponent β: 
(i) given exponents η and ν from the FSS of ℳ, one finds β = ν(η + 1)/2, and
(ii) using the correlation length exponent ν and Fisher exponent τ from the FSS of ℳ and P_bulk(ℓ), one finds  β = 3ν(τ - 2)/(τ-1).
We find that both approaches give values consistent with one another and clearly distinct in the two symmetry classes.
Moreover, hyperscaling gives values of β which are in very good agreement with those found from the FSS analysis of the watermelon correlator G_2(r,K) (see Tab. <ref>).
For consistency we also report β as determined by the hyperscaling relation β = ν(η+1)/2 in Tab. <ref>.

Lastly, let us comment on how the exponents we determine here compare with previous numerical studies.
The orientable loop model was studied previously in Ref. Ortuno_2009.
Our results for the measurement-only Kitaev model yields critical exponents perfectly consistent with these earlier results 
(see Tab. <ref>), confirming the universality of the entanglement transition.
With our loop model simulations alloing us to reach appreciably larger system sizes (by two orders of magnitude) 
we are able to refine the estimated critical exponents to greater precision. 



By contrast, the non-orientable loop model was studied only relatively recently in Ref. Serna_2021.
Unlike the orientable case, the critical exponents we report for this symmetry class seem to differ appreciably from 
these previous estimates.
In particular, whereas we find ν=0.9403± 0.0006, Ref. Serna_2021 suggested a smaller value ν=0.918 ± 0.005.
This in turn leads to a systematic difference in all other critical exponents found by hyperscaling relations or FSS which involves ν [
For instance, consider the anomalous dimension η as found by FSS of ℳ.
From the RG flow <cit.>, one expects that ν^-1 - 2 + η is a constant.
Thus, a small error in the fitted value of ν is compensated by a corresponding shift in η.
For Δ(ν^-1) = (0.9403^-1 - 0.918^-1), one finds Δη = 0.026, almost precisely matching the difference in reported values of η.
Since the fractal dimension was determined in Ref. Serna_2021 via the scaling relation d_f = (5-η)/2, we find that the discrepancy in d_f can also be accounted for using a similar reasoning.
Finally, the assumption that ν^-1 + η is constant in tandem with the scaling relation β=ν(1 + η) allows one to show that Δβ=Δν = 0.0223, consistent with the values reported in Tab. <ref>.
].
To reconcile this apparent discrepancy, we have simulated two alternative circuit/lattice geometries, including the L-lattice studied in Ref. Serna_2021, which has allowed us to track down its source: In our approach we have expanded the FSS ansatz by non-linear terms, which we find to give superior data collapses over wider parameter ranges and system sizes. As discussed in Appendix <ref>, this also allows us to smoothly interpolate between our estimates and the ones found in Ref. <cit.> as a function of (i) the range of included system sizes (ii) the choice of scaling variable, and (iii) the density of spline points. This leads to conclude that an improved FSS ansatz for the numerical finite-size data of Ref. Serna_2021 would give a somewhat higher estimate of the critical exponent ν, perfectly consistent with what we report here.

Despite the discrepancy with earlier works, there is strong reason to believe that the critical exponent values we report here represent the most accurate and precise determination for this symmetry class to date.
Our numerical simulations reach system sizes which are much larger than previously studied, enabling greater control over subleading corrections in the FSS analysis.
This is further bolstered by the internal consistency between the two independent approaches for estimating ν.
In Appendix. <ref> we further demonstrate the universality of ν≈ 0.94 by examining two additional lattice geometries, namely the measurement-only Yao-Kivelson model and Cardy's 3D L-lattice.





 §.§ Simulation setup




Let us close this Section by highlighting the benefits of the loop model representation in numerical simulations. 
In short, the loop model representation allows for a substantially more efficient numerical simulation of Majorana circuits <cit.> reaching system sizes of up to 10^8 qubits far exceeding those accessible via standard Clifford circuit methods (employed and discussed in the next Section). 
To understand this efficiency gain of about two orders of magnitude in number of simulated qubits remember that, in the tableau representation of conventional Clifford simulations, 
an N qubit stabilizer state requires 𝒪(N^2) memory, and each Clifford operation (e.g. a measurement) requires 𝒪(N^2) time.
By contrast, in the loop model representation,  the loop connectivity compactly encodes the state in only 𝒪(N) space and can be updated after a measurement in 𝒪(1) time.

This substantial computational advantage of the loop model representation is further enhanced by noting that a circuit of depth t acts as a transfer matrix on the loop connectivity.
Two such depth t circuits can be efficiently composed in time 𝒪(N), yielding a new circuit of depth 2t.

The loop model simulations are then carried out as follows:
We begin by preparing a “pool” of N_p shallow circuits of depth t=1.
In particular, we measure N bonds from the Majorana honeycomb lattice chosen at random with relative probabilities K_α for nearest-neighbor bonds and J for next-nearest-neighbor bonds.
Then all trajectories are evolved to depth T by log_2(T) rounds of concatenation.
When concatenating, the circuits are randomly translated with respect to one another in order to enlarge the accessible configuration space as well as to reduce spurious correlations arising from finite N_p.
From this pool, N_s independent sample trajectories are drawn by choosing two random elements, translating them with respect to one another, and concatenating.
Due to self-averaging in large systems, random translations before concatenating allows one to take N_s ≫ N_p.

For each such sample, quantities such as the entanglement entropy and the spanning number may be directly computed from the loop connectivity after imposing appropriate boundary conditions at the temporal boundaries.
Moreover, for PBC, the entanglement entropy averaged over all cuts can be efficiently computed in one 𝒪(N) shot.
By contrast, computing entanglement entropy in the tableau (Clifford) representation requires 𝒪(N^3) time for each subsystem, making it far more costly.
To probe bulk loop statistics, we track not only the loop connectivity, but also the total length of each loop, storing the lengths of closed loops in a histogram during the evolution.

Our simulations are completely based on self-implemented loop model codes (as opposed to codes using highly optimized open-source community packages as in the case of the Clifford simulations, see next Section). These codes are run on national high-performance computing resources, specifically the AMD EPYC (v3 Milan)-based Noctua2 cluster at the Paderborn Center for Parallel Computing (PC^2), the AMD EPYC (v4 Genoa)-based RAMSES cluster at RRZK/University of Cologne, and the Intel XEON Platinum 8168-based JUWELS cluster at FZ Julich.
Limited by the amount of system memory available, we initialize pool sizes between 120 (L=128) and 20 (L=4096) and generate 10^5 samples from a given configuration pool. Depending on the system size, we generate between 50 and 100 independent pools.

Let us provide an example of the resulting computational cost: at each set of parameters in Fig. <ref>, generating 10^6 samples consumes 3 core-hours at L=128 and 1000 core-hours at L=4096. For both systems, across all system sizes and 100 different parameters, a budget of less than 200,000 core-hours was sufficient to generate the data shown in Fig. <ref>. Using a budget of about 2000 core-hours per data point, this has allowed us to simulate system sizes up to N ∼ 10^8 qubits in the loop representation (versus N ∼ 10^4 in the Clifford tableau representation). An important technical detail is that our loop model code is highly parallelizable and limited in performance primarily by memory bandwidth, inviting GPU acceleration in future simulations.

 



§ CLIFFORD SIMULATIONS


In this final section, we provide numerical results from conventional Clifford simulations relying on stabilizer tableau algebra. 
Such simulations often serve as a starting point for a quick exploration of entanglement phase diagrams before indulging into
more complex numerical simulations (such as the loop model numerics discussed in the previous Section) or analytical 
considerations. It is the purpose of this Section to demonstrate what kind of qualitative and quantitative insight can already be
afforded when performing such Clifford simulations using state-of-the-art codes <cit.>, 
with which one can readily simulate systems of the order of 10^4 qubits. 






 §.§ Phase diagrams from tripartite mutual information





One particularly powerful observable to quickly map out the entanglement phase diagram of a given Clifford circuit is the tripartite mutual information

    I_3(A:B:C) = I_2(A:B) + I_2(A:C)-I_2(A:BC)                        
    
                   =S_A + S_B + S_C - S_AB - S_BC - S_AC + S_ABC .
 
where the subregions A,…, D are equal partitions of the torus on which the lattice is defined. 

Already for a fixed, finite linear system size, this observable allows to distinguish the entanglement structure of the localized and liquid phases,
for which is assumes values <cit.> of I_3(A:B:C) = +1 (localized) and I_3(A:B:C) = -1 (liquid), respectively.
As we show in Fig. <ref> for the original honeycomb Kitaev circuit and its Kekulé variant, a heat map-style visualization of I_3 readily maps out the entanglement phase diagram already with reasonable accuracy (even though the critical point does not necessarily coincide with I_3 = 0) – compare, e.g., with the phase diagrams in the overview of Fig. <ref>. 





A more precise quantitative estimate of the phase boundaries can be obtained via a finite-size scaling analysis of the tripartite mutual information, which exhibits an easy-to-locate crossing point for different system sizes at the phase transition K=K_c.
Performing a finite-size data collapse as shown in Fig. <ref> one can further extract the correlation length exponent, though with far less accuracy as what is achievable with the loop model framework discussed in the previous section.


The values of the correlation length exponent ν are clearly distinct between the two symmetry classes, even at these small system sizes. 
However, the specific values do not agree with our previous results from the loop simulations, even within margin of error. In the Clifford simulations, 
there are several limitations at play which affect the accuracy of the extracted correlation length exponents and which are likely to lead to this disagreement. This includes 
(i) strong finite-size effects arising from the limited system sizes available in the Clifford simulations (which are hard to overcome with the current algorithms/numerical platforms), and
(ii) a necessity to extend the first-order scaling ansatz x = (K - K_c) L^1/ν to higher orders (which is not straight-forward).
One manifestation of these limitations is a system-size dependent drift of the tripartite mutual information within the liquid phase, 
making it impossible to properly collapse the data. The origin of this drift can be rationalized as follows: even deep in the liquid phase, there is always a non-zero probability that the system will not be long-range entangled at a given time. In essence, a series of measurements on the K bonds can localize the Majorana pairs and temporarily destroy long-range entanglement. Alternatively, this can be phrased as a random walk in configuration space with a small fraction of area-law states, which vanishes in the thermodynamic limit. In taking the average over a large number of samples, this small fraction of localized (area-law) states then manifests itself as a small, finite offset in the tripartite information and is most pronounced for small system sizes and more so for the orientable circuit (due to the reduced number of paths for long-range strings). This is indeed visible in Fig. <ref>, where I_3 never actually saturates to -1 (in the liquid phase). 
In our finite-size scaling analysis we have worked around this by performing the optimization of the fitting cost function only for data points at K > K_c.










 §.§ Entanglement entropy and Lifshitz criticality


Another way to explore the entanglement phase diagram of a Kitaev circuit and to zoom in on the phase transitions is, 
of course, to directly explore the entanglement structure of the ensemble of stationary states. 
To this end, one might want to look at “entanglement arcs",
i.e. the bipartite entanglement entropy for subsystem cuts of varying size 0 ≤ l/L ≤ 1. This is shown in Fig. <ref>
for parameter scans through the phase diagram of the honeycomb Kitaev circuit as an example of an orientable circuit as well as the non-orientable version with next-nearest neighbor parity checks. One can clearly distinguish the liquid region where the entanglement entropy S/L spans an arc akin to the well-known Calabrese-Cardy result <cit.>  

    S(l/L)/L = a  logsin(π l/L) + b  ,

for conformally invariant systems (red traces), versus the localized phase whose area-law scaling reveals itself as a wide plateau of constant entanglement entropy (blue traces). 



When moving towards the phase boundary within the entanglement phase diagram, a close inspection of the scaling of the bipartite entanglement entropy reveals the characteristic shape of Lifshitz critical scaling, marked by the pink trace in Fig. <ref>, which shows a clearly distinct behavior from the “entanglement arcs" in the liquid regime and the almost constant behavior in the localized phase. This deviation from the liquid-like scaling is also reflected in an enlarged χ^2 value of the respective fit, as shown in the right panels of Fig. <ref>. 
Zooming in on the critical behavior and plotting the entanglement entropy for different system sizes in Fig. <ref>, we find already for the system sizes achievable in our Clifford simulations 
an indisputable signature of the Lifshitz scaling, consistent with our discussion in the context of the loop model simulations in Section <ref> above.



These findings support the usefulness of Clifford simulations to not only map out the entanglement phase diagram 
of a given monitored circuit with relatively modest numerical resources (for instance via scans of the tripartite mutual information), but also to extract qualitative features of the phase transitions (such as a distinction of ∼ Llog L versus 
Lifshitz scaling) and even get a rough estimate of critical exponents. At the same time, high-quality fits of critical 
exponents require substantially larger system sizes than attainable in Clifford simulations (∼20,000 qubits), leaving 
an opening for more advanced numerical schemes (such as the loop model simulations) that allow to go beyond 
these limitations.




 §.§ Simulation setup




Let us close also this Section on Clifford simulations with a detailed review of our numerical approach. 
To obtain the stationary state as efficiently as possible, we stabilize all conserved quantities at initial time. Depending on system size, we take between 5 × 10^2 and 4× 10^4 samples for each set of parameters K_x,y,z, J, which determine the probability to projectively measure a bond (see Fig. <ref>). We use a highly efficient implementation <cit.> of stabilizer tableau algebra. While it is not possible to improve upon either the 𝒪(N^2) scaling for performing a single measurement (and thus 𝒪(N^3) per circuit layer) or the 𝒪(N^3) scaling for canonicalization, we perform several HPC optimizations to reduce memory access, since tableau-based simulations are always bound by memory bandwidth. This includes a hybrid parallelization strategy based on the physical layout of the hardware used, which in this case is a 50 × 2 × AMD EPYC 9654 cluster.
Furthermore, we disregard phase information and measurement outcomes – ensuring that measurements always scale ∼𝒪(N^2) – and make use of the fact that there are 3L/4 possible partitions of the honeycomb lattice, allowing to perform self-averaging. Even though computing the entanglement entropy of a stabilizer state always requires to canonicalize the tableau (𝒪(N^3)), it is possible to get a slight advantage when the only quantity of interest is the scaling of S(l/L) by converting to the clipped gauge. While this also requires an 𝒪(N^3) operation, it allows one to extract the entanglement entropy for all subsystem sizes in one shot. This trick, however, is of little use for the tripartite mutual information, since it requires tracing out non-contiguous regions – which is not possible in the clipped gauge representation. In the end, we are able to generate high quality data for up to 2× 96^2=18,432 qubits, albeit at exploding computational cost. 
As an example, for a sweep of 100 different parameters, generating 192 trajectories with 150 samples each consumes close to 2,000 core hours at L=32 while generating only 96 trajectories with 5 samples each takes slightly more than half a million core hours at L=96 – the majority of which is however spent on performing canonicalization.







§ DISCUSSION



Summarizing our main results in brief, we have investigated measurement-only circuit variants of the Kitaev and Yao-Kivelson model to elucidate that 
(i) monitored quantum circuits give rise to robust entanglement phases, i.e. dynamic phases of matter that are stable to local variations of the circuit parameters (measurement probabilities in our case), 
(ii) the notions of symmetry and dimensionality can be used, akin to their Hamiltonian counterparts, to define universal behavior of the entanglement phases and the quantum criticality of phase transitions between them, and
(iii) that the particle-hole symmetric Majorana circuits, relevant to these monitored Kitaev models, can be recast in terms of Majorana loop models.
These Majorana loop models are a powerful resource as they not only allow to recast the entanglement phases in terms of localization physics of loops and map their phase transitions to non-linear sigma models, but they also provide a numerical framework allowing to simulate up to 10^8 qubits. This has allowed us to clearly discriminate the universality classes of orientable versus non-orientable loop models in (2+1) dimensions, relevant to Majorana circuits in symmetry classes BDI and D, respectively and determine their critical exponents to unprecedented precision.

Furthermore, our work suggests a broader perspective on the tenfold way classification of free fermionic wave functions. The connection between orientable and non-orientable loop models offers a unified framework for embedding both free Majorana circuits (fugacity n=1) and ground states of free Majorana Hamiltonians (fugacity n=√(2)) within symmetry classes BDI and D. This raises intriguing possibilities for future research, such as exploring whether other symmetry classes allow similar embeddings or investigating higher loop fugacities, potentially realized in quantum circuits coupled to classical dynamical agents <cit.>. This approach could refine the tenfold way classification by distinguishing universality not only by symmetry classes but also by Hamiltonian dynamics versus monitored circuits. Moreover, this offers a tractable handle on broader emerging questions regarding symmetry as an organizing principle in open quantum systems <cit.>.

One unexpected discovery in our study is the presence of quantum Lifshitz scaling at the phase transition of monitored Majorana circuits—a feature previously associated only with interacting circuits. This finding suggests a new and unexplored connection between loop models in (2+1) dimensions, monitored free fermions, and (2+1)-dimensional conformal field theories, warranting further investigation.

There are many future directions one might pursue in studying universal behavior in monitored quantum circuits and measurement-only Kitaev models. 
This includes variations of dimensionality, such as going to higher spatial dimensions and studying three-dimensional monitored Kitaev models <cit.>, which is also closely related to the 3D Floquet code <cit.>, or exploring quantum systems with arbitrary qudit dimension <cit.>. 
One might also think of generalizations of monitored Majorana systems to more general monitored non-Abelian anyon systems, such as a monitored Fibonacci anyons in a measurement-only variant of the golden chain <cit.>.

This work also raises a number of questions regarding the relation between the universality of entanglement transitions in Majorana circuits and those in more generic (2+1)-dimensional Clifford circuits.
Moving away from the Gaussian limit, we obtain random quantum circuits which can be related to a broad class of interacting loop or dimer models and for which new strong-coupling fixed-points may emerge.
Clifford operations which are even (odd) under fermionic-parity will preserve (break) the orientability of the loop model <cit.>.
Exploration of the parity-symmetric case has been initiated already in Ref. Zhu23structuredVolumeLaw, where the one finds volume-law entanglement and a correlation length exponent consistent with the loop model in class BDI.
Far from the Gaussian limit, large-scale numerical simulation of generic (2+1)-dimensional Clifford circuits give compelling evidence that the measurement-induced transition is in the same universality class as 3D percolation <cit.>.
However, the exponent ν=0.876 ± 0.004 <cit.> for 3D percolation is markedly smaller than what we find here for non-orientable Majorana circuits via loop simulations.
As is apparent from our own Clifford simulations in class D, finite-size effects lead to an underestimation of ν even for the largest system sizes.
It would be of interest to better resolve whether this apparent difference is due to finite-size effects or if perturbing away from the Gaussian Majorana limit leads to a genuine renormalization of the critical exponents.





Let us close this outlook by pointing to some of the physics awaiting us when going beyond the (random) Clifford regime 
of strong, project measurements.
To do so, one can consider non-random but weak, non-projective measurements 
–  as realized by moving the rotation t∈[0,π/4] in the unitary entangling gates, drawn as CNOT gates in Fig. <ref>, 
away from a perfect π/4-rotation. 
One crucial implication here is that moving away from a perfect t=π/4 rotation is akin to reducing τ, 
the effective inverse temperature in the Kraus operator (<ref>) above, to a finite value.
This connection between weak measurements and finite temperature physics has been demonstrated <cit.> in rigorous terms 
for commuting two-qubit checks, resulting in Nishimori physics that subsequently has been probed in experimental devices <cit.>. 
In the context of the Hastings-Haah code it allows one to explore the physics of qubit fractionalization
 (and intermediate phases) as one goes out of the Clifford limit of projective measurements to weak measurements <cit.>. 
 It will be interesting to investigate such weakly monitored Kitaev circuits in more depth in the future. 
 For instance, note that analogous to cooling a Hamiltonian in equilibrium, tuning the measurement strength from weak to strong can purify all the time reversal invariant plaquettes or the product of plaquette operators, leaving a Kramers doublet as a time reversal symmetric mixed state. 
 Applied to the non-orientable circuits of interest in this manuscript this poses an interesting question for future research – is the spontaneous breaking of time-reversal symmetry in such non-orientable circuits captured by a hitherto unexplored chiral Nishimori critical point?






    KK was supported by an NSF Graduate Fellowship under Grant No. DGE 2146752 and by the U.S. Department of Energy, Office of Science, National Quantum Information Science Research Centers, Quantum Science Center.
    
    The Cologne group acknowledges support from the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy Cluster of Excellence Matter and Light for Quantum Computing (ML4Q) EXC 2004/1 390534769, and by the DFG Collaborative Research Center (CRC) 183 Project No. 277101999 - projects B01 and B02.
    
    GYZ additionally acknowledges support from a start-up fund from HKUST(GZ). 
    
    The numerical simulations have been performed on the JUWELS cluster at Forschungszentrum Jülich,
    the Noctua2 cluster at the Paderborn Center for Parallel Computing (PC^2), and the RAMSES cluster at RRZK, University of Cologne.








Data availability.– 
All code for generating the figures shown and the underlying numerical data is available on Zenodo <cit.>.













§ HAMILTONIAN PHASE DIAGRAMS



Here for completeness we show, in Fig. <ref>, a comparison of the Hamiltonian phase diagrams of the original honeycomb Kitaev model and the Kekulé-Kitaev model. 

For the isotropic honeycomb Kitaev model, a unit-cell consists of two spin sites, and in the momentum space it hosts two linearly dispersive Majorana Dirac cones at the Brillouin zone corners <cit.>. The Dirac cones are stable in the presence of translation and C_2 𝒯 symmetry (a combination of lattice inversion and time reversal transformation). The perturbation away from the isotropic point gradually shifts the location of the two Dirac cones, which merge and annihilate and gap out at the phase boundary between the gapped and gapless phase. 

In contrast, the Kekulé variant of the model tricolors the bonds and triples the size of the unit-cell into 6 spin sites, which induces a periodic potential with large wave vector that can scatter one Dirac cone to the other, which immediately gaps out the Dirac cones when perturbed away from the isotropic point <cit.>. As a result, the gapless phase shrinks to a singular point in the phase diagram. 




The distinct topology of the phase diagram of the Kekulé-Kitaev model allows for an adiabatic path connecting the three gapped phases. This feature is exploited in the Floquet code <cit.>, which measures all the plaquette operators by four rounds of bond measurements without collapsing the global Wilson loops as logical qubits, and inducing an effective Hadamard transformation for the logical qubit due to the winding around the gapless point. 






§ FRUSTRATION GRAPHS



For generic protocols involving only Majorana parity checks, we can generate the frustration graph by a duality transformation. The original matter Majorana fermion spans a graph where each node hosts a matter Majorana fermion, and each edge corresponds to a parity check. To dualize the graph, we turn every edge to a node, and every Majorana site into a all-to-all connected cluster, see Fig. <ref>. In the dual graph, the node refers to the measured operators, and the edge connects two operators that anti-commutate each other, due to the fact that they overlap with a Majorana fermion. 

If the original graph has uniform coordinate number, say, every matter Majorana fermion has z neighbors, then in the dual graph, each node as an operator has 2(z-1) neighbors. For example, in the (decorated) honeycomb lattice z=3 and the frustration graph has 4 neighbors. For the special case of next-nearest neighbor measurement only, the Majorana fermions form two disconnected triangle lattice with z=6, and thus the coordinate number of the dual graph is 10. 








§ UNIVERSALITY AND LATTICE GEOMETRY IN NON-ORIENTABLE SYMMETRY CLASS D



Here we provide additional evidence for the universality of the critical exponents reported in Tab. <ref> for the non-orientable symmetry class.
We consider two additional lattice geometries: 
(i) the measurement-only Yao-Kivelson model, and
(ii) Cardy's L-lattice, as studied in Ref. Serna_2021.
Results are summarized in Tab. <ref>, showing a striking agreement between the three geometries as expected for universal
critical behavior.











 §.§ Measurement-only Yao-Kivelson Model



Let us now consider the decorated honeycomb lattice of the measurement-only Yao-Kivelson model depicted in Fig. <ref>(d).
We take a lattice of L × L large (dodecagonal) plaquettes such that the corresponding squashed brickwall geometry has dimensions L_x = 6L and L_y = L.
Moreover, we fix parameters K_α = K and J_α = J = (1 - K) for simplicity such that for 0 < K < 1 the loop model is non-orientable.
Then varying K allows one to tune the system through the localization transition for symmetry class D.





  §.§.§ Correlation length exponent ν






In Fig. <ref> we show the data collapse for the measurement-only Yao-Kivelson circuit in the vicinity of the critical point.
From this FSS analysis we find correlation length exponent ν = 0.9401 ± 0.0015, which is in close agreement with the value found both on the honeycomb lattice with NNN couplings and on Cardy's L-lattice.
This result further bolsters the reliability of our reported value of ν and demonstrates universality of the transition independent of the underlying lattice geometry.




  §.§.§ Fisher exponent τ




In order to find the remaining critical exponents, we determine the Fisher exponent τ from the bulk loop length distribution P_bulk(ℓ) at the transition.
Fitting the power-law scaling P_bulk(ℓ) ∼ℓ^-τ shown in Fig. <ref> yields an estimate of τ which is consistent with the corresponding result on the other two lattice geometries.
As for Cardy's L-lattice, we report the critical exponents η, d_f and β as found by hyperscaling relations in Tab. <ref>.






  §.§.§ Bulk loop statistics







Finally, we highlight universal aspects of bulk loop statistics in the metallic phase.
In Fig. <ref> we show the distribution of bulk loop lengths P_bulk(ℓ) at the isotropic point K = J deep in the critical (metallic) phase.
These results mirror those presented in Fig. <ref> for the honeycomb lattice, revealing Brownian and Poisson-Dirichlet scaling regimes.
Moreover, the asymptotic scaling in the Poisson-Dirichlet regime is consistent with Eq. (<ref>) for both symmetry classes.
We further validate the universality of the metallic phase by examining ratios of the long-distance loop connectivities.
In Fig. <ref>, we show that these ratios rapidly approach the values predicted from the universal quantity θ, consistent with the result in Fig. <ref> for the honeycomb lattice.









 §.§ Cardy's L-lattice




As noted in the main text, the value of the critical exponent ν reported here differs notably from that reported previously in Ref. Serna_2021, which employed loop simulations on Cardy's L-lattice <cit.>.
Here we re-examine the non-orientable loop model on the L-lattice in order to validate our results from the honeycomb lattice with NNN couplings.

The lattice geometry and loop model parameterization are defined in Fig. <ref>.
Taking a time slice through the unit cell, we see that only four of the loops may contribute to the spanning number.
We thus define the cross-sectional area of the unit cell to be 2 × 2 rather than 4 × 4 as was used by Ref. Serna_2021.
For linear dimension L, we take the 3D lattice to have unit aspect ratio, corresponding to size L × L × L.
Here we focus on the line p=q, along which the loop model is always non-orientable.









  §.§.§ Correlation length exponent ν








We compute the spanning number n_s for a set of 100 points in close proximity to the transition for linear dimensions L = {64,128,256,512,1024,2048,4096}.
For each parameter p=q, we average over 100 independent “pools”, drawing 10^4 circuit trajectories from each such pool.
In order to fit the scaling function f(x) we again employ a cubic B-spline.
Here we find it useful to take the spline points more densely clustered near x=0 where the curvature of f(x) is greatest.
In particular, we follow Ref. Ortuno_2009, taking spline points at spacing Δ x = 0.1 for x ∈ [-0.5, 1.5] and Δ x = 1 for x ∈ [2, 5].
While this choice of spline points gives a more robust fit, taking uniformly distributed spline points yields consistent results.
In Fig. <ref>, we show the scaling collapse resulting from two different ansatze.
For scaling variable x = L^1/ν(p-p_c), we obtain a good collapse with ν=0.9318±0.004 when restricting to a narrow window about x=0 and system sizes L ≥ L_min = 256.
When allowing for a non-linear correction to the scaling variable x = L^1/ν(p-p_c)(1 + A(p-p_c)), we find a nearly perfect collapse with ν=0.939 ± 0.001 for all system sizes and a much wider interval around the transition point.
Note that this value of ν also nicely matches the one found for the non-orientable circuit on the honeycomb lattice with NNN couplings.
Moreover, this FSS analysis is consistent with the earlier observation that a non-linear term in the scaling variable is necessary for robust scaling collapse.

To further validate the reported value of ν, we carry out a careful examination of the fitting procedure.
In Fig. <ref>, we show the quality of fit and optimal value of ν as a function of (i) the minimum system size L_min, (ii) the choice of scaling variable, and (iii) the density of spline points.
We find that the simpler (linear) ansatz consistently underestimates the exponent ν relative to the non-linear ansatz except in the limit of very large system sizes.
Moreover, we see that the non-linear ansatz consistently results in a robust scaling collapse which is largely independent of the density of spline points.
This analysis provides strong support for the universality of critical exponent ν = 0.9403 ± 0.0006 in non-orientable loop models, independent of the underlying lattice geometry.









  §.§.§ Fisher exponent τ



As an additional check, we determine the Fisher exponent τ at the critical point from the power-law scaling of the bulk loop length distribution P_bulk(ℓ) ∼ℓ^-τ.
Fixing the critical point from the FSS of the spanning number, we plot the distribution P_bulk(ℓ) in Fig. <ref>.
The results here may be compared with those for the honeycomb lattice shown in Fig. <ref>(b), showing comparable values of the exponent τ.
Taken together, the values of ν and τ determined here allow for the remaining critical exponents to be determined by hyperscaling relations.
In particular, we report the values of critical exponents η = 5-6τ-1, d_f = 3/(τ-1), and β=3ν(τ-2)/(τ-1) in Tab. <ref>.
All such exponents are in close agreement with the values found for the honeycomb lattice in the main text.









§ SUPPLEMENTAL NUMERICAL DATA





 §.§ Entanglement scaling in the metallic phase


Throughout the metallic phase of both the orientable and non-orientable models, we observe logarithmic entanglement scaling S(L,L_y)   L_y log(L) for a cylindrical subsystem of length L along the x̂ direction and circumference L_y.
This is accompanied by a power-law scaling of the surface loop length distribution P_surf^α(ℓ) ∼ℓ^-2.
In Fig. <ref>, we verify the expected scaling of P_surf^α(ℓ) and S(L, L_y) at the isotropic point.
Notably, the prefactor of the logarithmic entanglement scaling is larger in the non-orientable symmetry class as alleviating the orientability constraint yields greater mobility for the random walk undergone by loop endpoints.
Moreover, the relative magnitudes of the logarithmic term in the symmetry classes is comparable to that predicted by the ratio of the mean-squared displacement ⟨r⃗^2⟩ of the random walk.







 §.§ Bulk length of open arcs


As noted in Sec. <ref>, an open boundary is absorbing for the random walk undergone by loop endpoints.
Thus in the metallic phase, the bulk length of open arcs is dictated by the first passage time distribution P_FP(ℓ) ∝ℓ^-3/2 of the random walk along the temporal direction.
This is shown clearly in Fig. <ref> for both symmetry classes.
For finite system size L, the first passage time distribution acquires an exponential cutoff which is reflected in Fig. <ref> for loop length ℓ comparable to the cross-sectional area A.










