Entanglement and dynamics of diffusion-annihilation processes with Majorana defects

Coupling a many-body system to a thermal environment typically destroys the quantum coherence of its state, leading to an effective classical dynamics at the longest time scales. We show that systems with anyon-like defects can exhibit universal late-time dynamics that is stochastic, but fundamentally non-classical, because some of the quantum information about the state is topologically protected from the environment. Our coarse-grained model describes one-dimensional systems with domain-wall defects carrying Majorana modes. These defects undergo Brownian motion due to coupling with a bath. Since the fermion parity of a given pair of defects is nonlocal, it cannot be measured by the bath until the two defects happen to come into contact. We examine how such a system anneals to zero temperature via the diffusion and pairwise annihilation of Majorana defects, and we characterize the nontrivial entanglement structure that arises in such stochastic processes. We also investigate simplified"quantum measurement circuits"in one or more dimensions, involving repeated pairwise measurement of fermion parities. We show these can be solved by mappings to classical loop models. The relaxation of the system of Majorana defects to its ground state is analogous to coarsening in a classical 1D Ising model via domain wall annihilation (the classical"$A+A \rightarrow \emptyset$"reaction-diffusion process). Here, however, configurations are labeled not only by the defect positions but by a nonlocal entanglement structure. This"$\gamma+ \gamma \rightarrow \emptyset$'' process is a new universality class for the coarsening of topological domain walls, whose universal properties can be obtained from an exact mapping.


I. INTRODUCTION
Coupling a quantum many-body system to a thermal environment destroys the quantum coherence of its dynamics. Usually, a generic bath coupling leads to dynamics that is essentially classical at the longest timescales, such that the slow degrees of freedom are governed by a classical "hydrodynamics", either stochastic or deterministic. In this paper we exhibit a model that avoids this classical fate, because some of the information in the quantum state is topologically protected [1,2] from the environment. The late-time description is instead a hybrid of quantum and classical. Local degrees of freedom that can be "measured" by the bath become effectively classical, since they cannot remain long in a superposition. But the topologically protected Hilbert spaceassociated with anyon-like defects -allows long-range entanglement to persist. This entangled state in turn affects the slow "classical" degrees of freedom, which here are the positions of the defects.
The model we introduce is for a one-dimensional chain with defect excitations that carry Majorana zero modes. Topologically, these defects are like domain walls between the topological and trivial phases of Kitaev's pwave chain [1]. We model this system during a process of relaxation from a high-temperature state with many defects to a low temperature state where the defects are absent, via the diffusion and pairwise annihilation of defects. During this process, the system loses energy to the surrounding bath (which is taken to be bosonici.e. fermions cannot hop into or out of the system).
The diffusive motion of a single defect is entirely classical thanks to the bath coupling. But when the defects are widely separated -which is mostly the case at late times -the Hilbert space associated with the Majorana modes is protected from decoherence: all bosonic operators on this space are nonlocal, and so are invisible to the bath. Thus describing the state requires one to know both the positions of the defects and their entanglement structure.
One consequence of this hybrid classical/quantum dynamics is that the density of defects decays to zero with a universal amplitude that is different from what it would be if the defects were, say, featureless domain walls in the classical Ising model. The defect positions also exhibit different universal correlation functions. The relaxation of the one-dimensional Ising model to zero temperature under Glauber dynamics [3] occurs by diffusion and annihilation of domain walls -an instance of the "A + A → ∅" reaction-diffusion process [4]. This is a nontrivial stochastic process for which a wide range of universal quantities can be obtained exactly [4][5][6][7][8][9][10]. We find that the density of defects in the Majorana model, decaying via the "γ + γ → ∅ process", is exactly twice that of the purely classical case for a given (late) time.
The crucial feature of the model is a stochastic evolution of the quantum state of the Majorana zero modes γ i . This state is able to remain a pure state, despite the coupling to the environment, but its evolution is entirely nonunitary. In the long time limit of the diffusionannihilation process, the details of the environment coupling become unimportant and a simple universal description emerges. In this limit the environment sim-ply effects projective measurements of the fermion parity (iγ i γ i+1 ) of pairs of adjacent Majorana modes. Specifically, measurement takes place when two Majoranacarrying defects approach each other in space. At such moments, when random, classical diffusion happens to bring two Majoranas into contact, their mutual fermion parity becomes a local operator that can be seen by the bath. The result of this projective measurement determines whether or not the Majorana domain walls are allowed to annihilate back into the ground state.
Importantly, if the measurement outcome is such as to prevent annihilation, the quantum state stores this information until the next encounter. This additional "memory" effect means that the universality class of the dynamics is different from the classical A + A → ∅ universality class. Instead, we show that the γ + γ → ∅ process can be related to two copies of the A + A → ∅ process by an exact mapping.
Separately from this diffusion-annihilation dynamics, we study the process of repeated pairwise measurement of a fixed number of Majoranas, i.e. the evolution of the system in the absence of annihilation. This dynamics, which involves repeated measurement of the parity of adjacent Majorana pairs, is interesting in its own right. In particular, we show that it leads to a nontrivial critically entangled state.
In general, combining local unitary dynamics of a pure quantum state (say, for a lattice of spins) with repeated local measurements yields a phase transition in the dynamics of entanglement growth as a function of the measurement rate [11,12]. Various aspects of this transition are now understood [13][14][15][16][17][18][19][20]. In the generic case this phase transition is continuous, and is characterized by a nontrivial renormalization group (RG) fixed point with connections to unconventional statistical mechanics models [11,16,[21][22][23][24].
The projectively measured Majorana chain that we study here also flows to a nontrivial RG fixed point. In fact, so long as only pairs of Majoranas are measured, the dynamics is critical without the need to tune any parameter. However, the properties of this critical state are very different from (and simpler than) the generic case described above. The dynamics in our Majorana models should also be contrasted with those in models of free fermions subjected to single-site measurements, where a disentangled phase was found [13,17].
Strikingly, the repeated projective measurements organize the system of Majoranas into a random state for which the entanglement of a subregion scales logarithmically with the subregion size. In terms of its entanglement, this state is somewhat like a random singlet state [25,26]: it is described by an "arc diagram" like that in Fig. 1, for which only two-body entanglement is ever generated and the number of arcs crossing the midpoint of a system of length L is of order log L. But unlike the random singlet state, which is a ground state of a random Hamiltonian, the present state is prepared by a stochastic measurement process, without any need to minimize . . , γ12, defined by a choice of pairings (an "arc diagram") and an assignment nij = (iγiγj + 1)/2 = 0, 1 to each arc (connecting i and j, with j > i).
a Hamiltonian, and we show that its universal properties are different from those of the random singlet phase. In fact, these universal properties can be described exactly via an equivalence between the random "updates" to the pairing diagram produced by measurements, and the Temperley-Lieb operations [27] that appear in the context of a 2D lattice model for nonintersecting loops [28][29][30][31][32] and which are equivalent to a stochastic process on pairing diagrams [33][34][35]. We also relate the criticality of the stochastic measurement process to a Lieb-Schultz-Mattis-like theorem for disordered spin chains [36]. The quantum measurement circuit dynamics can be generalized to higher dimensions, and remains exactly solvable by a similar mapping to a classical loop model on a 3D lattice [37]. The diffusion-annihilation dynamics, however, is complicated by the possibility of braiding Majoranas around each other, which can alter the parity of Majorana pairs even when no measurements are performed. The models we study could also be generalized to other types of anyon, but we leave this generalization for the future.

II. MODELS AND SIMULATION
We will discuss two models. The first involves mobile Majorana defects which can diffuse and annihilate: we use this model to discuss a type of relaxational dynamics (Fig. 2). The relaxation we consider is to zero temperature, so that energy always goes down (or stays the same). This process leads to a length scale, the mean interparticle spacing, that diverges like √ t at late times (as in the relaxation of the classical 1D Ising model via the annihilation of domain walls). This diverging length scale ensures that the results are universal.
The second model is simpler: we assume that there is no annihilation between Majoranas. (For example, one can imagine an energy barrier arising from a repulsive short-ranged interaction, which prevents Majoranas from annihilating, but still allows them to approach each other close enough that their parity can be measured.) In this case the evolution of the system is determined only by a random sequence of local projective measurements (see Fig. ??). This is a form of "measurement-only" random quantum circuit. This second model is no longer a model for relaxational dynamics: we study it instead as an example of how a nontrivial entanglement structure can be produced solely by random measurement.

A. Diffusing Majorana defects
We study a model of Majorana defects diffusing in 1D. Let us first describe a coarse-grained model, and then discuss how it can emerge from a microscopic model. One possible physical interpretation of the defects (see Sec. II B) is as domain walls between the two topologically distinct phases of a Kitaev chain [1], with additional interactions to place it at a first order transition [38][39][40] between the trivial phase and the topological phase. When two Majoranas (black points) that are not paired with each other come into contact (horizontal arrows, above), the resulting local pair has fermion number either 1 (left lower configuration) or 0 (lower right), with each outcome having equal probability. In the latter case the pair is removed from the system. In either case the nonlocal pair has a fermion number such that the total parity of the system is conserved modulo 2.
some settings, these two states can be related by symmetry, so that no tuning of parameters is required to access the transition.) The number of defects N t must be even but is otherwise arbitrary, and can decrease with time t. Each defect i = 1, . . . , N t has a position x i (all distinct). For numerical convenience we put these positions on a 1D lattice of size L, with periodic boundary conditions, though one could also imagine a continuum model. To fully describe the state of the system we must specify the number of defects, their positions {x i }, and their quantum state.
Each defect contains a Majorana mode γ i (with anticommutators {γ i , γ j } = 2δ ij ). Any choice of pairings of Majoranas gives a basis for the total Hilbert space: a given pair i, j, with i < j, has a two-dimensional Hilbert space spanned by states with fermion parity iγ i γ j = ±1. We will label the parity by the fermion number operator n ij = (iγ i γ j + 1)/2 = 0, 1.
Thus, a labelled arc diagram such as Fig. 1 completely specifies the quantum state of the system, so long as the instantaneous state happens to be a basis state for some choice of pairing. It turns out, as we explain below, that we only need to consider states of this form (rather than superpositions of such states).
The rules for annihilation and measurement of Majoranas are simple because of universality in the long-time limit; see the remainder of this section for explanations, and Sec. II D for the precise rules of the numerical simulations on the lattice. The dynamics of the coarse-grained model are as follows.
Each isolated Majorana defect performs classical Brownian motion with diffusion constant D. When two Majoranas i and i + 1 come into contact, their fermion parity iγ i γ i+1 is a local operator and is visible to the bath, which makes a projective measurement that assigns a definite value to iγ i γ i+1 . If i and i + 1 are already paired, this measurement produces no update to the pairing. If instead i is paired with some other defect j, and i + 1 is paired with k, then the pairings are updated so that i is now paired with i + 1, and j with k. This re-pairing process is illustrated in Fig. 2. The fermion number n i,i+1 of the newly formed pair i and i + 1 is given a definite value, n i,i+1 = 0 or 1, with probability 1/2 for each option. The fermion parity of the pair (j, k) is fixed by Z 2 conservation of fermion parity: the sum of the values for the initial arcs, n i,j + n i+1,k , must be equal to the sum of the new values, n i,i+1 + n j,k , modulo 2. If we are allowing annihilation in the model, then the adjacent defects i and i + 1 annihilate if and only if their fermion number n i,i+1 is zero (see Fig. 2.) The two outcomes of the measurement of n i,i+1 are equally likely because the reduced density matrix for n i,i+1 is maximally mixed whenever n i,j and n i+1,k have definite values (as they do before the measurement). If the fermion number n i,i+1 is zero the defects can annihilate into the ground state. But if it is unity, conservation of fermion number modulo 2 prevents them from annihilating. In the latter case the defects continue to diffuse. Implicit in these dynamical rules is the hierarchy of energy scales E ∅ < 2E γ < E f , where E ∅ is the energy of the vacuum state, E γ is the energy of a single Majorana defect, and E f is the energy of a single fermion (the local "bound state" of two Majoranas). The first inequality in this chain implies that Majoranas with fermion number 0 annihilate, and the second implies that Majoranas with fermion number 1 remain separate.
The important feature of the diffusion-annihilation dynamics that differentiates it from the purely classical case is that when i and i+1 project into a state with n i,i+1 = 1 and diffuse apart, they remember that they cannot annihilate. Such a pair remains in a state of definite mutual parity until one of its two constituents encounters a third defect. Correspondingly, if i and i + 1 encounter each other again, without meeting anyone else in the meantime, they will again fail to annihilate.
At late times t the typical separation r(t) between defects is large, and the properties of 1D random walks imply that when two defects meet once, they in fact meet a parametrically large number of times (O(r) times) before either of them meets a third party. This large number of meetings justifies having no dependence of the protocol on a microscopic annihilation rate or a microscopic bath coupling. So long as annihilation is allowed by the fermion parity, it will definitely happen in a time short compared to the diffusive timescale (when r is large), irrespective of the microscopic rate. (See Sec. IV for a more detailed discussion.) The irrelevance of the microscopic annihilation rate is well-known in the relaxation of the 1D Ising model, for example: these relaxation processes are "diffusion-limited" rather than "rate-limited" [4].
For this same reason, we can treat the effect of the environment as a simple projective measurement. Even if the microscopic bath coupling is weak, it is amplified by the parametrically large number of repeated encounters.
Thus our model has only a single parameter, which is the diffusion constant D. We will be interested in the Definition of the order parameter Φ(x). Majorana defects (denoted by γ) separate adjacent ground state domains (red/green regions), which have Φ = ±1.
combined evolution of the positions {x i } and the entanglement structure. The latter is represented simply by the geometry of the arc diagram ( Fig. 1). In particular, we will pin down the universal late-time behaviour of the density, and we will compare this with a reference classical model in which the defects are purely classical variables, representing, for example, classical Ising domain walls.

B. Relation to microscopic models
The Majorana defects in our coarse-grained model live at domain walls separating two types of ground state, which we can label by a local "order parameter" Φ = ±1 (see Fig. 3). The two ground states are assumed to have equal energy density, so that the diffusion of a domain wall is unbiased. Depending on the microscopic model, this degeneracy of the two ground states could either be due to the system being tuned to a first-order transition, or it could be enforced by a symmetry relating the two ground states (so that Φ is the order parameter for the breaking of this symmetry).
For an example of how Majorana defects arise, consider the Kitaev chain, which describes a chain of fermions coupled to a superconducting order parameter. Its Hamiltonian may be written [1]: Here J ± δJ is the Majorana hopping amplitude, and we have labelled the sites x by half-integers for convenience. The "microscopic" Majorana operators η i should be distinguished from the variables γ used above in the coarse-grained model. The former are related to physical fermion operators c x by combining them in pairs, such that c x = 1 2 η x + iη x+1/2 for x ∈ Z. The two phases of the chain can be visualized simply in the limits δJ = J and δJ = −J. In each of these limits the chain is dimerised in one of the two possible ways [1], and the fermion parity of each dimer is zero. A single domain wall between the two phases carries an unpaired Majorana mode.
In the above model the quantum phase transition between topological and non-topological phases, which occurs at δJ = 0, is continuous and described by free fermions. However, adding interactions to the Kitaev chain can make this transition first-order. The details of these interactions are not important for the universal coarse-grained model, but a first-order transition is important in order to have well-defined domains. One method of driving the transition to first-order is by using a four-η coupling (the "Majorana Hubbard model"), as demonstrated in Refs. [38,39]. In this model the firstorder transition occurs at a large interaction strength; a different choice of four-η interaction yields the first order transition for much smaller values of the couplings [40]. The same transition can also be obtained in continuum theories, as proposed in the context of the 1D boundary of a 2D topological superconductor, with spontaneously broken time-reversal symmetry on the boundary [41].
Let us briefly comment on symmetry. In the Kitaev chain, and the interacting versions discussed above, there is a translation symmetry x → x + 1/2 that acts as an Ising symmetry on the "order parameter" Φ distinguishing the topological and trivial phases. In the usual interpretation of the Kitaev chain, where adjacent ηs are grouped to make local electron orbitals, this is not a physical symmetry [42], as is clear from the fact that it relates topologically distinct phases. In this interpretation, a parameter must be tuned to reach the transition point where the topological and trivial phases are degenerate in energy.
However, in other realizations of Majorana defects the two phases can be related by a physical symmetry. In the setup of Ref. [41], for example, where the 1D system comprises the boundary of a 2D bulk, this symmetry is time reversal. If the Majorana modes are realized by a 1D array of defects in a 2D topological system, the translational symmetry of the Kitaev chain can be a physical translation symmetry.
Many realizations of Majorana modes have been proposed (see for example Refs. [43,44] for an introduction), including at the vortex cores of 2D p + ip superconductors [45], as excitations of the Moore Read quantum Hall state [46], in proximity-coupled surfaces of 3D topological insulators [47], and at endpoints of proximatized 1D wires (see for example Ref. [48] for a review). For discussions of physical realizations of stronglyinteracting lattice Hamiltonians for Majoranas, see for example Refs. [49][50][51].
We consider a model with a first-order transition, tuned to the first-order transition point. (In the presence of a symmetry relating the two phases, this simply means we are in the ordered phase for Φ). Now we imagine that the chain is prepared at high temperature, and allowed to relax to low temperature via a coupling to a low-temperature bath. At late times the state will consist of large domains of one phase or the other which in their interior resemble one of the two degenerate ground states. The typical size of these domains grows with time, by the diffusion and annihilation of domain walls (Fig. 2). Since the system is poised right at the phase transition, neither phase is preferred, which means that the diffusion of domain walls is unbiased. (If we tuned slightly away from the transition, then one domain type would be preferred, and on average would grow at the expense of the other.) The model described in the previous section provides a universal description of this class of models in an appropriate limit of temperature and time scale. The temperature of the bath should not be strictly zero: in that limit the defects do not diffuse, and quantum effects other than the ones we consider may become important [52][53][54]. However, the temperature of the bath is assumed to be parametrically small, so that at equilibrium the concentration of thermally excited defects is also parametrically small. We can therefore consider an effective model on length scales below the typical spacing between such thermally-excited defects, at which the density is treated as relaxing to zero: this is the model described in the previous section. An important additional requirement is that the bath should not have low-energy (gapless) fermionic modes that can couple directly to the γs.

C. Quantum circuit model (no annihilation)
We may also consider the dynamics in Sec. II with diffusion but no annihilation. We will find that such dynamics lead to a critically entangled state. This is interesting to study even though it is not directly related to the relaxation process above. In this limit the defects diffuse around, being projectively measured whenever they bounce into each other. We present numerical simulations of this process in Sec. III A.
Such a model invites an even further simplification, in which the Majoranas γ i are taken to be at fixed posi-tions i ∈ Z and are subjected to repeated measurements of randomly chosen adjacent pairs. This procedure gives a quantum-circuit-like model for random projective measurement. The simplest way to set this "circuit" up is as in Fig. 4, so that at even time steps even-numbered links have the opportunity to be measured, with a probability p for a measurement to occur, and at odd time steps the odd links are measured with the same probability.
Note that this model is for the dynamics of a pure state: i.e. the measurement outcome is chosen randomly, but we do not average over it (which would give instead a mixed state density matrix). Conceptually, we can imagine that all the measurement outcomes are recorded by an observer, who therefore has access to the pure state. See Refs. [11][12][13][14]18] for discussions of these issues.
By representing the worldlines of the Majoranas in spacetime, the dynamics in this circuit model maps to a well-studied classical 2D loop model [29,30]. Equivalently, if we neglect the fermion parity labels, the dynamics on pairing diagrams in this model map to a known stochastic process derived from the loop model [33][34][35]. Exact results are available for these problems that determine the universal properties of the entanglement structure generated by the circuit (Sec. III C).

D. Simulation method
We implement a computer simulation of randomwalking Majorana defects on a discrete 1D lattice with periodic boundary conditions. The system is intialized with a large number N 0 of defects uniformly distributed across a system of size L, and the initial pairings are such that each Majorana is paired with a nearest neighbor (a dimerized configuration) and all pairs have initial parity n = 1.
The system evolves by random unit hops (left or right) of randomly-selected Majoranas, such that during each time step there are N t hops, where N t is the number of Majoranas in the system at the beginning of time step t. If a hop brings two Majoranas i and i + 1 onto the same site, then the hopping Majorana is reflected back to its initial site and the parity of the two Majoranas is updated according to the rules in Sec. II A: If i and i + 1 are already paired, then no updates are made to the pairings in the system; If i and i+1 are not already paired, they become paired upon contact. Their fermion number is set randomly to either n i,i+1 = 0 or n i,i+1 = 1 with equal probability. The previous partners of i and i + 1 are joined into a (nonlocal) pair, whose parity is fixed by fermion parity conservation (conservation of n modulo 2). This is shown in Fig. 2.
Below we study two cases for the evolution of the number of Majoranas. In the first case (Sec. III), the number of Majoranas in the system is fixed. In this case the arc labels n are irrelevant to the dynamics of the pairing diagram, since the rules for pairing Majoranas do not de- pend on the parities of the arcs involved. In the second case (Sec. IV), two Majoranas are annihilated immediately if they come into contact and are measured in the n i,i+1 = 0 state. Details about the simulation parameters, including the system size and run time, are listed in the corresponding figures.

A. Overview and numerical results
We begin by considering the simpler model in which annihilation of Majoranas is not allowed. When the number of Majoranas in the system is fixed, there is still a nontrivial evolution of the arc diagram, i.e. of the wavefunction's entanglement structure. Starting from the initial state with only nearest-neighbor pairs, longer-ranged pairs form over time via the process illustrated in Fig. 2, as well as by the Brownian motion of endpoints. One can then ask: what distribution P ( ) of arc lengths arises from this dynamics, and over what characteristic time scale does it develop?
The arc length distribution is directly related to the amount of entanglement in the state. The von Neumman entanglement entropy S A of a subregion A is where ρ A is the reduced density matrix for the fermion orbitals c x in the subregion (see Sec. II B). Each arc that connects A to its exterior contributes half a bit of entanglement entropy. (In general there is also an order 1 contribution from the boundary of A, which we neglect.) If we choose half a bit as our unit of entropy, i.e. if we take the logarithm in Eq. (2) base √ 2, then S A is simply the number of arcs that leave region A.
Let A be a section of R contiguous lattice sites in an infinite chain, and let ρ be the concentration of Majoranas in the chain (in units of the inverse lattice spacing), so that Rρ is the mean number of Majoranas in A. The mean entanglement of A with its surroundings is equal to the mean number of arcs exiting the region, which is Here, is the integer length of the arc in units of the lattice spacing. Notice that, if P ( ) is a power law P ( ) ∝ 1/ α for arc lengths ρ −1 , then the value α = 2 is a critical case. When α > 2 the above formula gives only order 1 entanglement for Rρ 1, i.e. a "boundary law". On the other hand, a normalizable distribution with 1 < α < 2 gives S ∝ R 2−α . Precisely at α = 2, the entanglement entropy scales logarithmically with the region size.
We find numerically that the pairing dynamics is critical: α = 2. (Simulations in this subsection correspond to random-walking Majoranas; we will show below that the universal constants characterizing the entanglement are the same in the circuit model.) Figure 6 suggests that Therefore the entanglement scales logarithmically with a coefficient K. Figure 7 shows the results of a direct measurement of this entanglement, yielding the consistent estimate This logarithmic scaling of entanglement is due to each "scale" contributing order 1 bits of entanglement (as at a conformally-invariant 1D quantum critical point with The coefficient K in Eq. (5) is dimensionless and universal. This universality can be seen by mapping the quantum circuit model of Sec. II C, in a worldline representation, to a well-known classical loop model, or equivalently to a stochastic process based on the Temperley Lieb algebra [33][34][35]. We explain this mapping in detail in Sec. III C. One can then argue that the universal properties survive when the effects of diffusion are included (Appendix B). An exact result for the loop model in Ref. [32] (see also [56,57]) equates to with √ 3/π = 0.551..., in close agreement with our numerical result. Recall that we are measuring entanglement in units of half-bits, so that S A is precisely the mean number of arcs leaving S A . Equation (6) is for a region with two endpoints; for a finite chain of size L, split into two equal halves that meet at a single point, the entanglement is 1/2 times the value in Eq. (6), with R = L.
The correspondences mentioned above also determine the dynamical exponent z for the evolution of the entanglement: (This dynamical exponent is also explained via a simple heuristic picture for the dynamics in the beginning of the following subsection.) In Fig. 6 we see the gradual approach of P ( ) to the stationary ensemble with time, with longer arcs taking longer to be generated. At a given time t there is a length scale ξ(t) beyond which P ( ) decays exponentially. This length scale grows roughly linearly with time, consistent with the z = 1 scaling ξ(t) ∼ t 1/z = t. The structure of one arc per length scale, along with the dynamical exponent z = 1, implies that the temporal fluctuations of the entanglement, S(t), produce "1/f noise", such that |S(f )| 2 ∝ 1/f , where f is the frequency. Equivalently, in the time domain (and in the limit of infinite region size) This 1/f -noise is demonstrated using simulation data in Fig. 8. It can be understood by noting that the entanglement of some region with large size R has statistically equal contributions from each (logarithmically spaced) length scale up to R. Since the number of arcs of size changes on a timescale of order , this implies an equal contribution to the fluctuations from each time scale, which is equivalent to 1/f noise. 1 In its logarithmic entanglement and two-body pairing structure, the critically entangled state found above resembles random-singlet wavefunctions. These arise 1 The equal contribution from each temporal scale, i.e. from each df |S(f )| 2 is independent of k. This gives the 1/f scaling. as ground states of 1D chains with quenched disorder [25,26,[58][59][60][61][62][63]. In the context of disordered spin-1/2 antiferromagnets, arc diagrams like Fig. 5 are used to represent frozen long-range spin singlets in the ground state. Similar states arise as ground states [62] (and even highly excited states [64,65]) of the Kitaev chain [Eq. (1)] with quenched disorder, with arcs representing pairs of Majoranas of definite fermion parity.
Interestingly, however, the universal properties of these random singlet states are distinct from those of the state obtained dynamically here. This difference can be seen from the universal coefficient K in the entanglement entropy. In the random singlet states (both for spin-1/2 fermions and for Majorana fermions) the number of arcs leaving a region of size R in an infinite chain is [61] 1 3 ln R, in contrast to Eq. (6). Crudely speaking, this difference in coefficients means that the random singlet ground state of a Majorana chain is less entangled than the state obtained from the repeated pairwise measurement process. The scaling properties of this dynamical state are related to a two-dimensional nonunitary conformal field theory: this structure is presumably quite different from that governing the random singlet phase.

B. Note on translation symmetry
Let us comment on the role of translation symmetry in these ensembles of wavefunctions. It is interesting to note that the value of the exponent α = 2 exhibited by both of them -both the random states obtained by measurement, and those obtained as ground states of random Majorana Hamiltonians -is the largest that is possible under the assumption that the ensemble does not spontaneously break "translation symmetry" in the index i. 2 This symmetry is a property of our dynamics with diffusion and measurement, 3 and it is also a statistical symmetry for the random Majorana Hamiltonians (i.e., a symmetry of the ensemble of Hamiltonians) and similar random spin-1/2 models.
The basic point is that in order to avoid breaking translational symmetry, the pairing diagram must have many large pairs and an average entanglement that diverges with the system size. This statement follows from a much more general mathematical theorem in Ref. [66], and for random singlet states it follows from a more general result for random spin-1/2 chains [36]. But the result for pairing diagrams is a simple special case and follows directly from a basic property of such diagrams [67,68], as we now discuss.
If all pairs are short-ranged, the parity P i of the number of arcs "passing over" a given link (i, i + 1) between Majoranas i and i+1 alternates as a function of i for any state in the ensemble [67]. (This alternation can easily be seen by drawing a diagram.) Therefore -again assuming all pairs are short-ranged -this parity acts as a local "dimerization" order parameter, implying a spontaneous breaking of translation symmetry.
The key question is then about the locality of this order parameter P i . We define spontaneous dimerization to be present in the pairing diagram if there exists a quantity O i (the translated version of O 0 ) which depends only on pairing information within a finite widow around i, and whose two-point function O i O j has nondecaying period-2 oscillations at arbitrarily large |j − i|. (Note that we do not demand that O i O j be expressible as a correlation function of quantum operators of the form ψ| O qm i O qm j |ψ in the Majorana chain, and in general it will not be.) P i does not automatically satisfy the above definition, since it may include contributions from arbitrarily large arcs. But it is straightforward to argue that, if the mean number of arcs crossing a given link in the infinite chain is finite -i.e. if the mean entanglement is finite -then contributions from large arcs are rare, and P i may be approximated arbitrarily well by an operator O i that is local in the above sense.
In other words, statistical translational symmetry guarantees long-range entanglement in both the random singlet phase and in the stochastic measurement process studied here. Translation symmetry can also be used to constrain the entanglement for various other examples of measurement dynamics (Sec. VI).
For completeness, let us note that the criticallyentangled ensembles mentioned above are qualitatively different from a simpler ensemble in which we choose a pairing diagram at random, with uniform probability, from all possible non-crossing pairing diagrams. 4 This ensemble has the arc length distribution P ( ) ∝ 1/ 3/2 (reviewed in Appendix A). A typical arc diagram from this uniform distribution is presented in Fig. 9. Note that this state looks qualitatively different from the result of our stochastic dynamics (Fig. 5), and has many tightly-nested large arcs.  (Fig. 5), we show a configuration from the equally-weighted ensemble of all noncrossing configurations. This ensemble has arc length distribution P ( ) ∝ 1/ 3/2 , and so is statistically very different from the configurations generated in the Majorana model.

C. Analytical treatment
The entanglement structure, encoded in the pairing diagram, undergoes a stochastic evolution that is nontrivial because of its nonlocality: when defects i and i + 1 come into contact and become paired, their previous partners j and k also simultaneously become paired, irrespective of their spatial separation.
This nonlocality means that large arcs can grow much faster than they would if their endpoints simply diffused. For example, if an L-sized system of O(L) defects is initialized with only short-range pairs, then arcs of order L size are generated in a time of order L z = L. This value z = 1 for the dynamical exponent can be understood using a known mapping from pairing diagrams [27,28] and their stochastic updates [33][34][35] to isotropic 2D statistical mechanics models [29,30]. We will review this mapping below. First, however, we give an intuitive argument for the exponents governing the entanglement dynamics.
Assume to begin with that the stationary arc length distribution P ( ) ∝ 1/ α at ρ 1, where α is some exponent to be determined. Now imagine labeling a large arc in the system, and following the evolution of its length. We adopt the rule that after an update to the pairing, the label attaches to the larger of the two  Fig. 4 to a loop model (Fig. 11). The two possibilities, absence/presence of a measurement, map to two ways of connecting up the worldlines. arcs produced. The endpoints of the labelled arc move both by diffusion and by coalescence of the labelled arc with other arcs, whose lengths are distributed according to P ( ). Through this latter process the length of the labeled arc can take "steps" much longer than the mean interparticle spacing. So long as α < 3, the dynamics of the arc growth is dominated by rare long steps, and is therefore more like a Levy flight than a random walk. Considering the longest expected step in a time interval t gives a typical displacement x ∼ t 1/(α−1) . This sets the relationship between length and time scales: One can now fix α = 2 by arguing that either α < 2 or α > 2 leads to a contradiction. Briefly, either inequality leads to a dynamic instability that can be seen by considering the number of arcs with length ∈ (L/2, L) in a system of large size L. When α < 2, there are many such "L-sized arcs", which results in a close nesting of L-sized arcs (like that in Fig. 9, for example). Such arcs quickly annihilate each other by encounters of their endpoints, on a time scale much faster than they can be generated. 5 On the other hand, when α > 2, the expected number of L-sized arcs is much smaller than unity. In this case there is nothing preventing the longest arc in the system from growing to become comparable to the system size, so that again the system is unstable dynamically.
In short, the distribution P ( ) ∝ 1/ 2 is the critical case with an order-unity number of arcs of length of order in any subsystem of length ρ −1 , as mentioned above. Such a distribution provides enough long arcs to keep the upward growth of short arcs in check (via the nested annihilation process), but not so many that long arcs are forced to nest tightly inside each other, which would lead to them rapidly annihilating. This structure with "one arc for every scale of length" implies the logarithmic entanglement scaling in Eq. (5), and the resulting Levy flight of arc endpoints guarantees z = 1.
A more precise picture can be obtained by applying results from 2D "loop models" and related stochastic processes. To make this mapping, we simplify the diffusionplus-measurement process to the quantum measurement circuit model shown in Fig. 4, where Majorana modes γ i live at fixed positions i ∈ Z, and adjacent pairs (γ i , γ i+1 ) are measured with probability p in either even or odd time steps, depending on whether i is even or odd. At the end of this section we argue that the universal results also apply in the model with diffusion.
If we neglect the fermion parity numbers, the "updates" to the pairing diagram are precisely those in the Temperley-Lieb transfer matrix representation of a twodimensional lattice model for nonintersecting loops related to percolation [27,28]. In the loop model, pairing diagrams, representing loop connectivities, can be used as a basis for the transfer matrix. It was pointed out in Refs. [33][34][35] that this transfer matrix can be interpreted as the transition matrix of the stochastic process on pairing diagrams. 6 (The transfer matrix of a lattice model cannot, in general, be re-interpreted as a stochastic process, but in geometrical models like percolation, with uncorrelated local degrees of freedom, it can be. 7 ) The relation to 2D loop configurations can be seen directly as follows. In the circuit diagram of Fig. 4, each block is either a measurement or a non-event, with probability p for the former. We represent each of these outcomes by replacing each horizontal bar with one of the two diagrams in Fig. 10: measurements become the reconnections shown in yellow in the loop diagram, and nonevents correspond to continuation of the vertical lines (blue). Doing this for every horizontal bar, we produce a configuration of loops, plus open strings, like that in Fig. 11. At the bottom boundary of this diagram, representing the initial time, the strands are connected in a pattern determined by the initial pairing of the Majoranas, which in the figure we have taken to be a dimerized configuration. The open strings in the configuration terminate on the top boundary, where they connect the sites i pairwise. Therefore, the loop configuration specifies a pairing configuration. It is easy to check that this pairing configuration is precisely the one resulting from the measurement dynamics. As noted above, the relationship between the 2D loop configurations and pairing diagrams is well-known [27,28], including the formal relationship with a stochastic process [33][34][35]: here we provide a physical realization of this stochastic process in terms of measurement in a Majorana circuit.
When the measurement probability is p = 1/2, the 2D loop configurations generated above are drawn from a simple ensemble where every allowed configuration of loops is equally likely. This ensemble is scale-invariant. The loops have the same fractal structure as cluster boundaries in critical percolation [31], as can be seen from an explicit lattice construction (see Refs. [29,30] for reviews). The ensemble is isotropic in 2D, and this isotropy is preserved in the scaling limit even if the measurement probability is different from 1/2, so long as we define the unit of time appropriately. Therefore the dynamical exponent for the pairing dynamics is z = 1 [33]. The entanglement entropy maps to the number of strands that connect section A of the top boundary to other parts of the top boundary. This number has been computed exactly in Ref. [32], giving the universal coefficient quoted in Eq. (6). (This quantity has also been considered in other settings in Ref. [56] and Ref. [57], where a heuristic similarity with entanglement was noted.) The length distribution P ( ) is (twice) the probability that two boundary sites at a distance are connected, which is a boundary two-point function in the conformal field theory language that decays as 1/ 2 , consistent with the discussion around Eqs. (3) and (4).
Above, the probability of measurement was the same on the even and odd bonds of the Majorana chain. As an aside, we note that if the probability is "dimerized" so that measurements are more frequent on (say) odd bonds, this drives the loop model into a non-critical phase with short loops. The corresponding correlation length exponent is that of percolation in two dimensions, ν = 4/3. The fact that the model with equal measurement rates is critical is related to translation symmetry in the Majorana index, as discussed in the previous section. By contrast, the 2+1D models we introduce in Sec. III C (and a variant 1+1D model discussed at the end of that section) can be long-range entangled even in the absence of such a symmetry.
Finally, to complete the analysis, we must confirm that reintroducing the diffusive motion of the Majoranas (but not the annihilation process) does not change the univer-sal properties obtained for the circuit. This follows from a simple renormalization group argument and is confirmed by our numerics.
Specifically, in the case with diffusion we have a loop ensemble similar to the one shown above, but rather than being defined on a regular lattice, its structure is determined by the encounters of the diffusing particles. In Appendix B we argue that this difference introduces only renormalization-group-irrelevant perturbations of the critical theory described above for the circuit. Therefore we do not expect the coupling to the diffusive density fluctuations to affect the leading scaling, though it will certainly contribute to subleading corrections. We now address the dynamics of the γ + γ → ∅ diffusion-annihilation process. In this model, Majorana defects that come into contact and have zero parity are removed from the system (Sec. II A). In the language of the microscopic quantum Hamiltonian discussed in Sec. II, this annihilation locally heals the ground state by removing two opposite domain walls. Thus, the model with annihilation describes a type of dissipative dynamics, in which the system is initiated in a high-energy state with many domains, and gradually relaxes to a ground state with no domain walls.
As mentioned above, this dynamics differs from a standard classical model for relaxation of the Ising model [3][4][5][6][7][8][9][10] by the fact that Majoranas with fermion number n = 1 are prevented from annihilating: in such cases the domain walls cannot heal upon contact (see Fig. 2).
Let us first review the standard classical Ising problem, which is equivalent to the "A+A → ∅" reaction-diffusion process. In this problem, A single species of particles undergoes diffusion on the real line, and when two particles come into contact they annihilate with some probability θ. For this process one can show that, regardless of the value of θ, and of the initial value of the particle density, the density ρ = N t /L at sufficiently long times t follows where D is the diffusion constant. This exact result, in which 1/ √ 8π is a universal amplitude, has been obtained in various ways [5][6][7]: one approach is to map the Markov transition matrix to an integrable non-Hermitian quantum spin chain [7]. However while the value of the universal amplitude requires a more sophisticated analysis, the scaling ρ ∝ 1/ √ Dt can be understood simply as follows. When the density of particles ρ 1, the typical spacing between nearest neighbors is 1/ρ 1. Thus, the typical time required for one particle to diffuse into its nearest neighbor is ∆t ∼ 1/(ρ 2 D). In the limit where particles are very sparse, two particles that come into contact have many opportunities to collide and annihilate before one of them can diffuse away and encounter a third party. (In the absence of annihilation, two diffusing particles which collide go on to collide ∼ 1/ρ times before either of them can diffuse a distance ∼ 1/ρ.) Thus, at small ρ, two particles that collide will almost certainly annihilate before either can escape the other, regardless of the value of θ. This means that the value of θ is irrelevant at late times (one can say that it is renormalized to unity). Therefore an order 1 fraction of the particles are annihilated in the characteristic diffusive time ∆t ∼ 1/(ρ 2 D) 1. Equivalently, the change ∆ρ in this time is of the order of ρ itself. Relating (∆ρ)/(∆t) to the derivative dρ/dt gives dρ/dt ∼ −Dρ 3 , whose solution is ρ ∼ 1/ √ Dt. (Note that this relation for dρ/dt should be distinguished from the mean field rate equation, which would give a different, incorrect, scaling; see e.g. Ref. 4.) In our problem of diffusing Majoranas, the power law ρ ∝ 1/ √ Dt still holds, since the same time scale ∆t is associated with diffusion over the nearest-neighbor distance, and again an order one fraction of particles annihilate in this time window. However, Eq. (10) no longer holds: the fraction of particles that annihilates in the time window is smaller, and the universal coefficient in Eq. (10) increases.
Unlike in the classical case, where two colliding particles will almost certainly annihilate if t is large, two colliding Majoranas that happen to be projected into a state of parity n = 1 remember this parity when the two particles collide again: no subsequent collision between them can produce annihilation, without first involving a third party. One should therefore expect the particle density at long times to be larger than Eq. (10). Figure 12 shows our numerical result for ρ(t) in the Majorana model with annihilation. For comparison we also plot the density from a numerical simulation of the classical model with θ = 1/2. The latter model would apply if, instead of remembering the fermion parities, we were to re-choose them at random for a pair each time they collided. The classical model has a density ρ cl (t) given at long times by Eq. (10).
The inset of Fig. 12 shows that the Majorana model produces a density which at long times is precisely twice that of the classical reaction-diffusion problem: We explain this surprising relation in the next subsection. It is worth emphasizing that this factor of 2 difference from the classical model cannot be understood simply as a probability for annihilation of two defects that is reduced by (say) a factor of 2 on average by the parity constraint. Indeed, as noted above, the microscopic annihilation probability θ has no effect on the density at late times in the classical model.

B. Mapping to two copies of A + A → ∅
Remarkably, the universality class of the Majorana diffusion-annihilation process can be characterized exactly by a mapping to an auxiliary classical diffusionannihilation process that does not involve a nonlocal pairing structure. Instead, this classical model contains additional "fictitious" labels on the particles. These labels have no meaning in the original quantum system, but by averaging over them, we reproduce the exact dynamics of the defect positions in the original quantum system. This classical model reduces to two independent copies of the A + A → ∅ process at late times.
The classical model that we define involves two particle types, which we label A and B, and the allowed annihilation processes A + A → ∅ and B + B → ∅. The classical model undergoes a stochastic process that is fully local (i.e., unlike the Majorana process, it does not require any pairing information). The mapping between the two models is simple: loosely speaking, each Majorana is replaced randomly with either an A or B particle, but in a way that is constrained by the fermion parities in the pairing diagram. This is illustrated in Fig. 13. If a pair of Majoranas has n = 0, the two Majoranas are replaced by two classical particles with the same label, while a pair with n = 1 is replaced by two classical particles with opposite labels. We define the mapping precisely below.
For concreteness, let us set up the Majorana dynamics on the 1D lattice as follows, in continuous time. In each infinitesimal time window dt, a given bond of the lattice has a probability Γdt of receiving an "update", where Γ is a rate. 8 Updating a bond means the following takes place: (i) If neither of the sites on the bond is occupied, nothing happens.
(ii) If one of the two sites on the bond is occupied, the occupying Majorana defect hops across the bond.
(iii) If the two sites are both occupied, and the corresponding Majoranas are in a "1" pair, nothing happens.
(iv) If the two sites are both occupied, and the Majoranas are in a "0" pair, they are annihilated.
(v) If the two sites are both occupied, but the Majoranas are not paired with each other, they are either paired into a "1" pair, or annihilated, with probability 1/2 for each option; their previous partners are now paired with the appropriate fermion parity.
We now introduce a different stochastic model. We refer to it in this section as the "classical" model, to contrast it with the "quantum" dynamics, i.e. the dynamics in the Majorana model that involves the pairing information. Again we consider particles that occupy the sites of a 1D lattice, with a given site occupied by at most one particle. In this classical model there is no pairing of particles, or any nonlocal structure. However, each particle now carries a new label, which is either A or B. Again, bonds of the lattice are updated at rate Γ. When a bond is updated, the following occurs: (i) If neither of the sites on the bond is occupied, nothing happens.
(ii) If the one of the two sites on the bond is occupied, the occupying particle hops across the bond.
(iii) If the two sites are occupied by two particles with opposite labels, then the outcomes are AB or BA with probability 1/2 each: i.e. the particles exchange places with probability 1/2.
(iv) If the two sites are occupied by particles with the same label (AA or BB) these are annihilated.
We now describe the mapping between the Majorana model and classical model. As mentioned above, this involves replacing each Majorana with either an A or B particle (Fig. 13). The assignments are made separately for each pair in the pairing diagram (note that the two constituents of the pair may be spatially distant). If the pair has n = 0, its two Majoranas are replaced with classical particles with the same label: both A, or both B, with probability 1/2 for each option. If the pair has n = 1, its two Majoranas are replaced with classical particles having opposite labels, with a probability 1/2 for the left member of the pair to be the A particle and a probability 1/2 for the right member to be the A particle. Note that this mapping does not change the positions of the particles: it only replaces the pairing structure with a random assignment of As and Bs. To describe things more formally, consider the evolving probability distribution in each of the processes. Let P qm t be the probability distribution for the state of the Majorana process at time t. This distribution tells us the probability that there are defects at given positions, with a given pairing structure, and given fermion parities. Let P cl t be the probability distribution for the state of the classical process at time t: this distribution tells us the probability that there are particles at given positions, with given labels.
The prescription in the paragraph before last allows us to map a probability distribution in the Majorana model to a probability distribution in the classical model. We call this map M. Abusing notation slightly, we write: Above, we defined a mapping for a single state of the Majorana model: this corresponds to the special case where P qm t is "delta function" supported on a single state. Since a general distribution is a linear superposition of such delta functions, the extension to a general P qm t is straightforward, with M being a linear map.
The key point is that this map is compatible with the time evolution we have defined for each of the models. Formally, this means that the following diagram commutes: That is to say, it does not matter whether we make the "quantum-classical mapping" at the very beginning, for the initial conditions, and then run the classical dynamics until the final time; or whether we instead run the Majorana dynamics until the final time, and make the classical mapping at the end. In either case we get the same final probability distribution in the classical problem.
Furthermore, the mapping preserves the positions of the defects. Therefore, as long as we choose initial conditions correctly in the classical model, the probability distribution of the total density is identical in the two models for all t.
Proving that the diagram (13) commutes is a matter of checking the various cases, corresponding to the various possible updates on a bond.
First note that, by the linearity of all the maps involved in (13), it is sufficient to check the case where the Majorana model is in a definite initial state (i.e. where P qm t is supported on a single state). Also, it is sufficient to check commutativity for the update operation on a single bond. We just have to go through the possible initial states for defects on the bond.
For brevity, we discuss a couple of illustrative cases below. It is straightforward to check the remaining cases, which completes the proof.
If the updated bond is unoccupied, or occupied by a single 'particle' which then hops, the result is immediate, so we consider the cases where the updated bond is occupied by a pair of particles.
First consider a case where the two particles occupying the bond are already paired with each other, say into the '1' state. After the mapping M, we get two classical states with equal probability, which we denote If we update the left-hand side, in the Majorana problem, nothing happens, since the two Majoranas are in an n = 1 state and cannot annihilate. Therefore, for consistency with Eq. (13), the right-hand side of Eq. (14) should also be invariant under the fictitious classical time evolution we have introduced. This is easily seen to be the case using the update rule (iii) for the classical model given above.
Next, consider a case where the particles on the updated bond are not initially paired with each other. For example, let each be in a '0' pair. We also draw their partners: After the update in the Majorana model, which acts on the two central particles, the configuration on the left becomes (again, the coefficients are classical probabilities, not wavefunction amplitudes!): Mapping this final state [the right-hand side of Eq. (16)] to the classical model using M gives: For consistency, we must obtain the same state, Eq. (17), by applying the update in the classical model to the right hand side of Eq. (15). Using rules (i, iii, iv) for the classical process, we can easily check that this is the case. The other initial configurations may be checked similarly.
Having established this mapping, the problem is now reduced to understanding the dynamics of the classical model. This classical problem almost reduces to two separate diffusion-annihilation processes occurring in parallel, one for the As and one for the Bs. This is not quite the case, however: the update rule for a bond where the particle content is AB or BA shows that when As and Bs come into contact their hopping is correlated. (Roughly speaking, As and Bs interact through the restriction that they cannot simultaneously occupy the same lattice site.) However, this effect is irrelevant at late times when the particles are dilute. This diluteness means that an A particle is adjacent to a B particle only a parametrically small fraction of the time, so that the coarse-grained dynamics of an A particle (in the absence of other As) remains simple Brownian motion, with a diffusion constant D = Γa 2 set by the hopping rate Γ and the lattice spacing a.
Thus, at late times we effectively have decoupled A + A → ∅ and B + B → ∅ diffusion-annihilation processes. Each of these contributes a density given by Eq. (10) at late times (recall that this result is independent of the initial density). Summing these contributions gives Eq. (11).
So far we have assumed that the initial state is either a state with definite fermion parities for some choice of pairing, or a classical mixture of such states. One may ask how the density evolves if, instead, the initial state is a quantum superposition that cannot be written as a single pairing state, or an otherwise arbitrary density matrix. In Appendix C, we argue that the universal late-time results for the density and correlation functions (presented in the following subsection) are valid for any initial state. For example, one could take the initial state to be the ground state of the critical noninteracting Kitaev chain, representing a quantum quench for an open system in a certain limit.

C. Exact results for correlation functions
A considerable amount is known about the classical A + A → ∅ process, or equivalently about the relaxation of the 1D Ising model via domain wall annihilation, and these results may be translated to our present problem using the mapping above. We now discuss some examples.
In Ref. 8, Bray calculated the correlation function C cl (r; t) = S(x)S(x + r) in the 1D Ising model at late time t during the domain wall annihilation dynamics, where S(x) = ±1 is the Ising spin at position x. The result is a universal scaling form where D is the diffusion constant and erf is the error function. In terms of the diffusing particles (domain walls), S(x)S(x + r) is equal to (−1) m , where m is the number of domain walls in between x and x + r.
Returning to the microscopic models of Sec. II B, let us consider the Ising-like order parameter that distinguishes the two local ground states (for example of the interacting Kitaev chain). We have denoted this order parameter by Φ(x), with Φ(x) = ±1 for the two ground states (see Fig. 3). The correlation function is again the expectation value of (−1) m , where m is the number of Majorana domain walls in the interval (x, x + r). Making the classical mapping, this is equal to the expectation value of (−1) m A × (−1) m B , where m A is the number of A particles in between (x, x + r) and equivalently for m B . Since the A and B particles are independent at late times, we obtain the correlation function in the Kitaev chain as the product of two classical correlation functions: One can also consider non-equal-time correlation functions. For these, it turns out that the critical exponents differ between quantum and classical models. Let Previously we defined a map between the instantaneous probability distributions in the quantum and classical problems. But in fact one can go further and show that the probability measure for a complete history of the particle positions, from time 0 up to some final time t, is the same in the classical and the quantum problems. 9 Then a simple extension of the above argument relating the quantum and classical correlation functions 10 gives again The result for the classical Ising model [8], at zero spatial separation and with t ≤ t , is Therefore, when t /t 1, so that the power-law for the the long-time decay (t → ∞) is different in the two cases. Equations (20), for the equal-time correlator, and (22) for the non-equal-time correlator, are verified directly using simulation data in Fig. 14.
Finally, another interesting critical index is the "persistence exponent", Θ [9,10]. The probability that, during the dynamics up to time t, the order parameter Φ(x) at a given location has never flipped (i.e. the probability that no domain wall has ever traversed x) scales as t −Θ . The exact result in the classical Ising model is Θ cl = 3/8 [10]. Similar considerations to those above show that in the quantum model this exponent is doubled: Θ qm = 3/4.

D. Entanglement in the γ + γ → ∅ process
In the diffusion-annihilation model, the length scale associated with entanglement is inevitably large at late times simply because the interparticle separation grows as √ t. Even a "dimerized" pairing state would involve pairs of O( √ t) length. One can then ask: does the system contain pairs much longer than √ t, spanning many interparticle separations, at late times? It turns out that the answer to this question depends on the initial state.
In the following discussion we continue to assume that the initial state is either a state of definite pairing, or a statistical ensemble (mixture) of such states. This assumption excludes initial states that can only be written as quantum superpositions of different pairing states: we comment on these in Appendix C.
Let s ∈ {1, 3, 5, 7, . . .} be an arc "size" measured in terms of the Majorana index, so that a arc of size s pairs Majoranas i and i + s for some i. This definition removes In both plots the inset shows the same data in logarithmic or semi-logarithmic scale. All data in this figure is taken from a simulation with 10 5 lattice sites and an initial density ρ(0) = 1/5, and is averaged over 10 4 random realizations.
the trivial √ t scaling of the physical length which comes from the large interparticle separation. We will see in a moment that the leading effect of annihilation on the entanglement structure is to "advect" the arc size distribution, P (s), by carrying weight from large s to smaller s. This advection implies that if the initial state has a finite typical arc length ξ, no matter how large, the typical value of s at late times is only of order 1.
However if the initial state has ξ = ∞ -for example if the initial state is an equilibrated state of the model without annihilation, which has P (s) ∼ s −2 -then ξ always remains infinite. This latter result is consistent with the general fact that is impossible to have a short-range P (s) unless the pairing ensemble breaks "translational symmetry" i → i + 1 in the Majorana index (see Sec.III B). If we draw the initial state from the critical P (s) ∼ s −2 ensemble, which does not break this symmetry, then by causality the symmetry must remain unbroken at any finite time, because symmetry breaking requires correlations over arbitrarily large distances. Therefore the tail of the size distribution must at all times decay slower than s −(2+ ) for any > 0.
To understand the "advection" of the arc length distribution, consider the evolution of a particular large arc with s 1. (As in the heuristic argument at the beginning of Sec. III C, in order to give the arcs a welldefined identity after a collision, we follow the larger of the two arcs produced by the collision.) The characteristic timescale for an encounter between this arc and another is the diffusive timescale ∆t ∼ 1/(Dρ 2 ), where ρ is the instantaneous density.
In an encounter, the chosen arc instantaneously grows or shrinks by an amount ∆s. If s 1, then typically ∆s s. On the other hand, on the same timescale ∆t, an order one fraction 1 − f of the interior Majoranas (lying underneath our chosen arc) collide and annihilate with their neighbours. This annihilation changes s by s → f s, since we measure arc sizes by counting Majoranas. This simple rescaling dynamics can be argued to dominate over the effect of collisions involving the long arc so long as α > 2. We can conclude that an arc of size s 1 will eventually shrink to a size s of order 1, and that the exponent for the power-law tail of the distribution is preserved (App. D).
If the initial distribution has a cutoff ξ on the largest value of s, then at late times all arcs will have s of order 1, and the typical entanglement length scale in the system will be set by the nearest-neighbor spacing. In our simulation with a dimerized initial condition we can directly measure the evolution of the arc-size distribution P (s), where s − 1 is the number of Majoranas underneath a given arc. We find (see Appendix D) that on a timescale of order 1/(Dρ 2 0 ), where ρ 0 is the initial density, the distribution becomes stationary, with a tail P (s) ∝ e −s/s0 at s 1 (s 0 ≈ 2.1). In order to examine the contrasting situation where the initial state has ξ = ∞, we have also run two-step simulations in which we first "equilibrate" the system using the dynamics without annihilation. This yields a power-law distribution P (s) ∼ s −2 that is cut off at large s only by the total number N 0 of Majoranas in the initial state. We then switch on annihilation and observe the evolution of P (s). What we see is roughly consistent with the crude picture above, in which the part of the distribution with s 1 is advected to smaller s at a rate set by the decay of the particle density (see Appendix D for simulation results). The distribution retains the initial power law in the range 1 s N t , where N t is the total num-ber of Majoranas at time t. This means that the longest arc in a system of physical length L system always has a physical length of order L, so long as a nonzero number of Majoranas remain in the system, although this arc's s value (which is of order N t ) is decreasing with time.
It is interesting to note that while the late-time density decay is insensitive to the initial state (see Sec. IV B and Appendix C), the late-time entanglement structure is not. This dependence of the entanglement on the initial state also implies that, in the case with annihilation, the quantum state of the Majorana defects that survive at time t in general never becomes a pure state if the initial state is a mixed one [18] (except in the trivial regime when t L 2 /D and no more defects remain). By contrast the measurement-only dynamics of Sec. III does eventually purify the initial state. In the loop picture (cf. Appendix C), the full system is completely purified when no loops connect to the initial time boundary. In the model with annihilation, this disconnection never happens while there are an extensive number of defects remaining. In the model without annihilation, it happens after a time of order L.

V. MAJORANAS IN TWO DIMENSIONS
There is intense interest in two-dimensional systems with Majorana-like defects, including superconducting systems with pointlike vortices that carry Majoranas [45,47], and topologically ordered systems (e.g. fractional quantum Hall states [46] and spin liquids [72]) with nonabelian "Ising" anyons, which are closely related to Majorana modes. Here we do not attempt to make direct connections with these systems, since it would require many other features to be taken into account. Instead we discuss the simplest possible two-dimensional generalizations of our models.
Perhaps surprisingly, the entanglement structure of a 2+1D quantum measurement circuit, in which Majorana modes are subjected to a random sequence of repeated projective measurements, remains exactly solvable in the scaling limit. We describe this exact solution below.
However, 2D versions of the diffusion-annihilation process are complicated by the nontrivial braiding of Majoranas. A convenient way to handle braiding is to imagine attaching branch cuts to each Majorana defect [73] (the geometry of these branch cuts is arbitrary). One then flips the sign of the operator γ i whenever the ith Majorana defect crosses the branch cut associated with another defect.
This procedure implies that the parity of a given pair depends on the path it has taken. As a result, the classical mapping of Sec. IV B does not carry over directly to the 2D case, and we do not attempt to treat the 2D γ + γ → ∅ problem here.
We note, however, that the simplicity of the braiding rule for Majoranas means that the 2D diffusionannihilation process could be studied numerically with low computational cost. The precise formulation of the problem will depend on the system being modelled. For example, defects may come in two species, representing vortices and antivortices, leading to an A + B → ∅ process.
Two dimensions is the upper critical dimension for such reaction diffusion problems [4], which exhibit logarithmic decay of the density with a universal constant. For example, in the case of the A + A → ∅ process [5,74] ρ(t) = 1 16π It would be interesting to obtain the corresponding universal constants for the quantum models.
Let us now consider the measurement-circuit dynamics for a fixed 2D lattice of Majoranas γ r . Here r runs over the sites of a 2D lattice that we are free to choose. We will find that there is a slight distinction between bipartite and non-bipartite lattices: we consider first the square lattice, as an example of the bipartite case. 11 To avoid having to choose an arbitrary lattice structure in 2+1D spacetime, we may imagine applying projective measurements to arcs at random in a Poissonian fashion (so that iγ r γ r+δ for each pair of adjacent sites is subjected to measurements at a rate Γ; of course only the sequence of measurements matters, not the precise times). Using a 2+1D spacetime lattice, in analogy with Fig. 4, would make no difference to the long-distance properties.
First, numerical results for the square lattice are shown in Fig. 15. The top panel shows a typical configuration of the arcs at late time, and the bottom panel shows the length distribution for a randomly chosen arc. We find that it fits well to at large , where is the Cartesian distance between endpoints of the arc. This power law is fortuitously the same as the 1D case, but the underlying universal structure is different. The power law can again be explained using a mapping to a classical statistical model of nonintersecting loops [37], now in 3D spacetime. This mapping goes through in complete analogy to the 1+1D case. We obtain loop configurations like that shown in Fig. 16, with a loop "turnaround" for every measurement event. 11 If we wish to think of the circuit model as a simplification of dynamics of mobile Brownian defects undergoing random spatial encounters, then the non-bipartite case is appropriate.
A field theory for 3D loop models of this kind was derived in [37,75]. The field theory is a "replica limit" of a nonlinear sigma model with continuous symmetry, in the phase 12 where the continous symmetry is spontaneously broken. (The symmetry is SU(N ) and the target space is complex projective space, CP N −1 . The replica limit is N → 1.) However, the details of this field theory need not concern us here, because the leading scaling turns out to be simple. The fact that the ordered phase can be described with weakly interacting Goldstone modes implies that a long loop resembles a Brownian path at large scales [37], somewhat like a polymer in a melt [76], with simple Brownian exponents.
This fact can be used to compute the length distribution P ( ) for the arcs in the stationary state. P ( ) is the probability that, if we follow a Brownian curve that is initiated at a point on the final time surface of the 2+1D region, its first return to this surface is at a Euclidean distance from the starting point. This simplifies because the various components of the Brownian path, corresponding to the three spacetime axes, are statistically independent. A simple calculation gives Eq. 26 at large . This is the same power law as in the 1D case, but at least in this approach that seems to be a coincidence (since in 1+1D a loop is not Brownian). In 2+1D it also leads to a different scaling of the entanglement. For a compact 2D region, say a disc of size R in an infinite 2D system, an appropriate integral of Eq. (26) shows that the entanglement scales as (a is the lattice spacing) with an O(1) prefactor that we do not determine here. This R ln R entanglement scaling coincides with that for a free fermion ground state with a Fermi surface [77][78][79], though the wavefunction (a translationally-invariant Slater determinant) is of a very different nature there.
The case with a non-bipartite spatial lattice is similar. Strictly speaking, the loop model is in a different universality class as compared to the bipartite case 13 [37,80]. This difference leads to a different field theory (the RP N −1 model with N → 1). However, there is again a phase where the loops are Brownian at long scales, leading to the same power law for P ( ).
A difference between the two cases appears if we dimerize the measurement rates, in analogy to the 1+1D discussion in Sec. III C. This dimerization can be used to drive a phase transition into an area law state with only short-range pairs. Unlike in the 1+1D model in Sec. III, the critical value of the dimerization is nonzero, so that even in the absence of translation symmetry there is a stable phase with the scaling in Eq. (27). The phase transition out of this phase is in a different universality class in bipartite and non-bipartite models.
These results in fact generalize immediately to spatial dimensions higher than two, where there is again a phase where the spacetime loops are Brownian, the arc length distribution is P ( ) ∼ −2 , and the entanglement of a d-dimensional ball is scales as (R/a) d−1 ln(R/a).
To complete the discussion of universality classes, let us describe a model that yields a different 1+1D universality class to that discussed in Sec. III. Such a model can be obtained by making the one-dimensional system effectively non-bipartite. The simplest way to do so is to randomly replace some of the measurements (of adjacent Majoranas) with "swap" operations. These are unitary transformations that exchange the Majoranas γ i and γ i+1 . 14 The loop model mapping of Sec. III C goes through as before, but now with an additional node configuration to be added to the two shown in Fig. 10. This configuration is one where the worldines of the two γs cross, which leads to a 2D loop model with crossings [71,81,82]. The crossings are a renormalization-grouprelevant perturbation, leading to a new universality class of loop gas. In this model we expect the half-system entanglement to scale not as ln L but as (ln L) 2 . 15 If we further perturb the model by "dimerizing" the measurement probabilities, so that even links are measured more frequently than odd links, then there is a continuous phase transition into a short-range entangled state at a nonzero critical value of the dimerization [71].

VI. OUTLOOK
In a system of nonabelian anyons or Majorana modes, some quantum information is protected from a decohering environment. We have shown that this protection leads to new universality classes for relaxation that are not entirely classical, since the active degrees of freedom at late time include a nonlocal entanglement structure. The γ + γ → ∅ coarsening process, which is solvable by the mapping in Sec. IV, is perhaps the simplest case, but we expect there are other interesting examples.
Our analysis of the dynamics of entanglement in this process has also led us to study related models involving measurement only. These models do not directly describe relaxation in many-body systems, but they illustrate that random measurements can lead to a nontrivial, criticallyentangled pure state. The late-time statistics of the state 14 Recall that the phase induced by the unitary operation is not important for the dynamics we discuss. 15 See the discussion of the "spanning number" in Ref. [71].
can be understood exactly, including a logarithmic violation of the area law for entanglement in any number of dimensions.
Let us suggest some directions for future research.
An interesting feature of the 1+1D measurement-only models studied here is the role of statistical translation symmetry by one Majorana. Of course, the individual states generated by the dynamics are not translationinvariant, because a particular realization of the measurement process yields random, non-translation-symmetric outcomes. However, as noted in the text, the translation symmetry of the ensemble of states still imposes a constraint on the entanglement structure, guaranteeing greater-than-boundary-law entanglement. Similar observations can be used to constrain the phase diagram of various dynamical protocols involving both measurements and unitary operations.
For example, consider a spin-1/2 chain subjected to unitary/measurement dynamics that is invariant under spin SU(2) symmetry (i.e. involving measurements only of spin-singlet operators) and under statistical translation symmetry. The stationary ensemble of states in the spin-zero sector must then have super-area-law entanglement, by a Lieb-Schultz-Mattis-like constraint on wavefunctions for spin-1/2s [36].
For the Majorana chain, heuristic arguments also suggest that statistical translation symmetry is enough to guarantee super-area-law entanglement for more general unitary/measurement dynamics, whether free or interacting (we will discuss this elsewhere). It would be interesting to generalize such statements to other settings.
The models considered here can be extended to other types of anyon, for example Fibonacci or SU(2) k anyons [65,83,84]. This kind of extension could give interesting stochastic models for entanglement structures represented by fusion trees. We expect that longrange-entangled states can be obtained by measurement here too, and that there are again interesting diffusionannhilation processes. Two-dimensional versions of these problems are also interesting, and have the additional ingredient of braiding, with braiding rules that become more complex when one goes beyond Majoranas. An important distinction will be between models where braiding [2] and measurement [85][86][87] together yield a computationally universal set of operations and models where they do not, since in the former case the full Hilbert space for a set of anyons can be explored, while in the latter case only a discrete subset can be accessed.
Even for Majoranas, which have the simplest braiding rules (such that the measurement-only entanglement dynamics remains tractable in 2D) braiding has a nontrivial effect on the 2D diffusion-annihilation process, which would be interesting to solve. operators that are irrelevant even when they are added to the action with uniform couplings [30,89,90].
Therefore we do not expect the coupling to the diffusive density fluctuations to affect the leading scaling of the entanglement, though it will will contribute to subleading corrections.
for all t. The crude advection picture in Sec. IV D can be argued to be self-consistent for α > 2, and for large s it gives: where f t = ρ(t)/ρ(0) is the fraction of particles that have not been annihilated. This suggests that if P 0 (s) is a power law s −α at large s then the exponent is preserved under the dynamics (and in the special case α = 2 the amplitude is also preserved).
In this Appendix we provide simulation results to support these claims by examining the cases where P 0 (s) is either short-ranged or given by the critical distribution P 0 (s) ∼ s −2 . Figure 17(a) shows the distribution P (s) as a function of time, with time t = 0 corresponding to a dimerized state, such that s = 1 for all pairs. As time increases, P (s) settles into a limiting distribution P (s) ∝ exp(−s/s 0 ), with the characteristic range s 0 ≈ 2.1. Figure 17(b), on the other hand, shows the time evolution of P (s) starting from the fully-equilibrated state of the model without annihilation. At time t = 0, the annihilation is turned on, and the number of Majoranas begins to decline. However, the distribution P (s) retains its functional form ∼ 1/s 2 , and is cut off only by the total number of Majoranas remaining in the system. Each curve is labeled by the corresponding simulation time. Data corresponding to t > 100 produce a curve P (s) that is indistinguishable the one labeled t = 100. (b) The total number of arcs having a given size s is plotted in double-logarithmic scale as a function of time, for a system that starts in the fully-equilibrated state of the model without annihilation. At time t = 0 the annihilation is turned on, and the total number of arcs starts to decline. However, the distribution of arc sizes retains the dependence 1/s 2 (illustrated by the dashed line). All data in this plot correspond to a simulation of 5000 Majoranas at time t = 0, and is averaged over 1000 independent realizations.