{'fraction_non_alphanumeric': 0.04512141496066839, 'fraction_numerical': 0.014535593882835785, 'mean_word_length': 4.202719877986782, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Dynamics and large deviation transitions of the XOR-Fredrickson-Andersen kinetically constrained model

We study a one-dimensional classical stochastic kinetically constrained model (KCM) inspired by Rydberg atoms in their"facilitated"regime, where sites can flip only if a single of their nearest neighbours is excited. We call this model"XOR-FA"to distinguish it from the standard Fredrickson-Andersen (FA) model. We describe the dynamics of the XOR-FA model, including its relation to simple exclusion processes in its domain wall representation. The interesting relaxation dynamics of the XOR-FA is related to the prominence of large dynamical fluctuations that lead to phase transitions between active and inactive dynamical phases as in other KCMs. By means of numerical tensor network methods we study in detail such transitions in the dynamical large deviation regime.


I. INTRODUCTION
Systems with constraints often display interesting cooperative dynamics. This is true both in classical and quantum settings. Broadly speaking there are three classes of constrained systems. One is that of problems where state space is constrained. The canonical example is lattice coverings, for example dimers on a square lattice. In such systems, the constrained nature of configuration space implies constraints in the allowed transitions, making both their classical and quantum dynamics very rich. A second class encompasses systems where constraints in the dynamics are emergent, such as in classical and quantum "fracton" models where the motion of certain effective excitations is severely impeded. A third class comprises systems known as kinetically constrained models (KCMs) with explicit constraints in the allowed dynamical transitions. Here we focus on KCMs.
KCMs were first introduced [1,2] in the 1980s as models of classical glasses. The ones studied most throroughly, such as the Fredrickson-Andersen (FA) [1] and East models [3], are stochastic lattice spin systems with the interesting combination of a trivial thermodynamics and a strongly fluctuating cooperative dynamics (under appropriate conditions -typically low temperatures and/or high densities) due to the constraints. For reviews on classical KCMs see e.g. [4][5][6]. Like their classical counterparts, quantum KCMs also display complex nonequilibrium dynamics, both under closed unitary [7][8][9] or open dissipative [10] evolution.
Modelling dynamics via KCMs can be motivated in many different areas. For example, in classical soft matter, specifically for glasses [11,12], kinetic constraints are meant to encode the local steric interactions of dense fluids. Another application is in the context of ensembles of Rydberg atoms in optical lattices, modelled as a collection of local two-level systems (representing for each atom their ground and some high-lying Rydberg state). When driven on resonance, due to "Rydberg blockade" [13], their dynamics is subject to a kinetic constraint where an atom can change state only if all their nearest neighbours are in their ground state. In a one-dimensional lattice such constraint gives rise to the much studied PXP model [14][15][16][17], the quantum counterpart of the classical "two-spin facilitated" FA model [4].
Here we study a one-dimensional classical KCM which to our knowledge has not been considered in the past. We call it the XOR-FA model to distinguish it from the standard FA model (i.e. the "one-spin facilitated" FA model). The kinetic constraint in the XOR-FA is such that a spin can flip only if one of its nearest neighbours is in the excited state, but not if both are (which is allowed in the FA). Such condition makes the XOR-FA more constrained than the standard FA model. Conversely, the XOR-FA is less constrained than the PXP, whose transitions require the two neighbouring sites to be simultaneously down. The constraint in the XOR-FA model can be motivated by Rydberg atoms in their "facilitated" (or "anti-blockade") regime [18][19][20][21][22][23][24][25][26][27]: when driven out of resonance, specifically when blue-detuned, conditions can be such that an atom may change state only if a single neighbour is in the excited state, but not both.
The paper is organised as follows. In Sec. II we introduce the XOR-FA model and discuss its connection to simple exclusion processes. In Sec. III we consider the relaxation dynamics of the model. In Sec. IV we study the dynamical large deviations by means of numerical tensor networks, and show the existence of a phase transition between active and inactive dynamical phases. In Sec. V we give our conclusions.

II. MODEL
We consider a system of binary variables n j = 0, 1 (we call these states down/up or unexcited/excited) on the sites i = 1, . . . , N of a one-dimensional lattice (with boundary conditions to be specified below). Similarly to other KCMs [4,28] the dynamics will be that of singespin flips subject to a constraint. Specifically, the allowed transitions are 001 → 011 rate = c 011 → 001 rate = 1 − c 100 → 110 rate = c 110 → 100 rate = 1 − c. (1) That is, a site can flip only if both nearest neighbouring sites are in different states. This means that the constraint is a boolean XOR operation on the nearest neighbours of the site that is attempting to flip. We therefore call this model the XOR-FA (short for XOR-Fredrickson-Andersen) to distinguish it from the standard Fredrickson-Andersen (FA) model, where a site can flip if either of its nearest neighbours is up, which in this nomenclature would correspond to the OR-FA (while the PXP would be the AND-FA).
The generator of the continuous-time Markov dynamics is the operator where σ ± j are Pauli operators acting on site j, n j = σ + j σ − j , and the kinetic constraint P j on site j reads, Dynamics with the kinetic constraint Eq. (3) is naturally motivated [25] in quantum many-body systems, specifically in the context of Rydberg atoms in their facilitated/anti-blockade regime [18][19][20][21][22][23][24][25][26][27], whereby an up (down) spin represents an atom in its excited Rydberg (ground) state, and the drive is such that an atom can get excited resonantly only when one of its nearest neighbours is also excited, but not both (as that would make the transition off resonant). The constraint Eq. (3) has also been studied in certain quantum spin chains [29,30] in particular in relation to "quantum scars" (special non-thermal states in constrained quantum systems [17,31,32]) [33][34][35]. Our aim here is consider the classical stochastic dynamics of a system with such a constraint, thus extending the set of known KCMs.
A. Conservation of the number of domain walls and relation to simple exclusion processes The dynamical rules Eq. (1) impose a conservation law in the dynamics, that of the total number of domain walls (DWs) [25,29]. Consider two neighbouring domains of, say, up and down spins · · · 11110000 · · · Due to the constraint Eq. (3) the only allowed changes are to the spins next to the DW, since inside the domains both neighbours to every spin are the same. This means that the possible moves are · · · 11100000 · · · · · · 11110000 · · · · · · 11111000 · · · where we have underlined the sites that changed in each allowed transition.
We can perform a duality transformation to have an explicit DW representation of the problem. We write where X j , Y j , Z j are Pauli operators for the DW between sites j − 1 and j. Notice that this is a canonical (rather than unitary) transformation that preserves the commutation relations between the Pauli operators. The generator in this representation is then where the constraint is and we have used the superscript "DW" to indicate operators in the domain wall representation. Combining the factors we can simplify the generator to The conservation law is now explicit, as the kinetic term simply corresponds to DW hopping. That is, we conserve the quantity N DW = j Z j . For the special case of c = 1/2, the generator simplifies to that of the symmetric simple exclusion process (SEP), [36,37] (10) In the XOR-FA language this is the "infinite temperature limit", where the cost of creating and destroying an excitation is the same. Away from c = 1/2, the dependence of c in the DW-representation generator Eq. (9) encodes the fact that moving left or right a DW depends on whether it is energetically favourable to extend or contract the corresponding spin domain.
For c = 1/2 the generator Eq. (9) corresponds to a SEP with particles with alternating asymmetries in their hopping rates. That is, we have a model where particles (DWs) can hop to neighbouring sites if the sites are not already occupied: the odd particles (DWs) hop left with rate c and right with rate 1 − c, while the even particles hop left with rate 1 − c and right with rate c. Since particles (DWs) cannot cross due to the exclusion, these rates are maintained. This is a special case of the general model introduced in Ref. [38], where each particle is given an individual hopping rate which is maintained under the dynamics.

A. Equilibrium properties
We consider the XOR-FA with N sites, N DW domain walls and periodic boundary conditions (PBC), which formally corresponds to setting n 0 = n N . The dynamics generated by Eq. (2) obeys detailed balance and therefore any initial condition eventually relaxes to an equilibrium state. Since the dynamics conserves the number of DWs, there is one such equilibrium probability for each DW sector. For PBC the number of DWs is even, and the sectors can be classified by the number p of up/down (one/zero) Note that the total length of the domains must be equal to the system size, so in the state above we have and each domain must have at least one site, so that We now define a state which is the translationally invariant superposition of all possible translations of the state above, where the operator T shifts the chain by a single site.
The equilibrium probability vector for the sector with 2p DWs is given by where γ = N − 2p + 1 and N is a normalization constant. One can check that the state Eq. (13) is annihilated by all terms of the generator Eq. (2). This state corresponds to the equilibrium state with non-interacting energy E = j n j (i.e., each up spin costs a unit of energy) at temperature T such that c = e −1/T /(1 + e −1/T ), and subject to the conditions Eqs. (11) and (12). We now study the basic properties of the equilibrium state Eq. (13). In Fig. 1 we show two average observables in equilibrium. The first one is the average excitation density, n = N −1 i −|n i |eq p , where −| = n n| is the flat state and n| = n 1 , . . . , n N |, see Fig. 1(a). We show n for several values of the filling fraction defined as n DW = 2p/N . The symbols are numerical results from standard Monte Carlo simulations. Note that in contrast to the FA or East models [4,28], n is not just equal to c, due to the conservation of the number of DWs. Figure 1(a) shows the agreement of the numerics with an analytical prediction in the N → ∞ limit described in the Appendix.
The second observable coincides with the average dynamical activity (per site) in equilibrium, k . While the activity is an observable at the level of trajectories (see Sect. IV below for further details), its average in equilibrium is given by the average escape rate per site, which is a static observable [39]. The escape rate operator R is (minus) the diagonal part of the generator Eq. (2). Since R is a local operator we can also obtain analytically its equilibrium average in the large size limit, see Appendix. In Fig. 1(b) we show the agreement between k from simulations and the analytic result. Note that n and k are symmetric around c = 1/2 as functions of c as a consequence of the up/down symmetry of the model, cf. Eq. (2) (while there is no corresponding symmetry in terms of DW filling n DW ).

B. Relaxation dynamics
The dynamics of the XOR-FA model is entirely determined by expansion and contraction of the domains (or the movement of domain walls which cannot cross). The system behaves like an "accordion". Depending on the value of c there may be an energetic preference to expand or contract domains of one orientation or the other.  is for c = 0.4, the middle one for c = 1/2, and the bottom one, c = 0.6. The columns correspond to different initial conditions. Column (a) shows equilibrium trajectories, i.e., those that start from an initial condition sampled from Eq. (13). They show pronounced space-time fluctuations in the dynamics associated to the breathing of domains. Column (b) corresponds to the most unfavourable initial state, where DWs are maximally clustered. Relaxation to equilibrium in this case is slow, as DWs within the bulk of the cluster cannot move until the DWs on the outside of cluster diffuse away. Column (c) shows an opposite non-equilibrium initial condition, where DWs are maximally spread out. In this case relaxation to equilibrium is faster, cf. Fig. 2(c). The large space-time fluctuations that are evident in these example trajectories anticipate the large deviation phase transitions that we uncover in the next section.
The different timescales involved in the relaxation of the XOR-FA model can be quantified using timecorrelation functions. In particular we focus on two different correlators in the equilibrium dynamics. The first one is the auto-correlation function, C(t), which mea-sures how many sites that were in the excited state at time 0 are also in an excited state at a later time t. Subtracting the disconnected part, and normalising so that it takes values between 1 and 0, it reads, where the average is over realisations of the dynamics in equilibrium, i.e., starting from a configuration sampled from the equilibrium state Eq. (13) and evolved according to Eq. (2). The second correlator we study is the persistence function, P (t), which quantifies the average probability for a randomly selected site to have not changed state up to time t. We can define it in terms of a local dynamical variable p j (t) at each site j, where p j (t) = 1 if the site has never changed from its initial state at time t, and p j (t) = 0 as soon as it changes for the first time. The resulting aggregate function is then This function is automatically normalised between 1 at the initial time and 0 eventually when all sites flip at least once. In Fig. 3 we show results for time-correlators. We focus mostly on the persistence function as it better captures overall relaxation. Figure 3(a) shows P (t) for various c and two filling fractions of DWs, n DW = 1/4 (top panel) and n DW = 1/2 (bottom panel). For comparison we also show the auto-correlator for c = 0.1 (dashed). We see that decreasing c away from c = 1/2 leads to slower relaxation times. The same can be said for decreasing the density of the DWs. Figure 3(b) shows the same functions as in the top panel of (a) but in a double-logarithmic scale on the ordinate. The change in slope in this representation emphasises the change from exponential decay at short times, to stretched exponential decay at long times.
From the persistence function we can extract a characteristic relaxation time, τ , customarily from the time the function decays to e −1 , that is, P (τ ) = e −1 . These times are shown in Fig. 3(c) for two values of c and as a function of the DW filling. Their behaviour can be understood approximately with simple heuristic arguments.
We first note that for c 1, we can treat the dynamics of the XOR-FA model as small up domains diffusing around a "vacuum" of down domains. To move, the up domain must first expand by exciting a neighbouring spin. This happens slowly at rate c. Following this, the domain then shrinks at rate 1 − c ≈ 1. It can either shrink back to its original position, or shrink such that it shifts by one site across, each happening with equal probability. Thus we say it diffuses around the lattice with diffusion constant D ≈ c/2. The time taken for the system to fully relax can then be estimated as the time it The same functions as in (a) (for nDW = 1/4), but we take the negative log for the y axis, and then also plot on loglog scales. (c) We plot the time it takes for the persistence function to drop to e −1 (markers) for c = 1/2 (blue) and c = 0.1(red) and various nDW. We also show our estimate τ (solid lines) as shown in the main text.
takes for the DWs to diffuse from their original positions around the available space surrounding them, until they hit another DW. On average, the length of each zero domain is given by the average number of down spins split among the number of zero domains. Namely, It then follows that the timescale for the system to relax goes as for c small. As Fig. 3(c) shows, this prediction works well for c small in the whole n DW range, while for c ≈ 1/2 it qualitatively accounts for τ for small DW density [40].
Thinking of the dynamics in this way can also explain the two timescales in Fig. 3(b). At some small time after t 0, the first successful shift of domain(s) will occur. When this happens for c 1/2, the original site is no longer excited, but the site next to it is. In the language of the persistence, this means two sites have flipped from their initial state. For both the persistence and the autocorrelator, this gives a fast initial relaxation, and as these are random uncorrelated events, the initial decay is exponential. Further successive moves of the domain only change at most one more site from its initial state (or in the case of the auto-correlator, will only slightly reduce the probability that the domain may end up in its original position). Thus the rate at which relaxation occurs is reduced, the time is longer, and the decay of the correlators is stretched as the relaxation becomes more collective.

IV. DYNAMICAL LARGE DEVIATIONS AND MATRIX PRODUCT STATES
In this section we study the statistics of trajectories of the XOR-FA model in the long-time regime where we can apply large deviation (LD) methods [6,41,42]. Recent work [43][44][45] has shown the effectiveness of numerical tensor network methods for studying the LDs of KCMs. Here, by means of numerical matrix product states (MPS) we are able to study the LDs of the XOR-FA for large systems to high accuracy. As we show below, the XOR-FA has a trajectory-space phase transition between between dynamical phases with very distinct characteristics, similar to what occurs in several other KCMs.

A. LDs and tilted generators
The dynamical activity [6,[46][47][48] is a trajectory observable which counts the number of configuration changes (in our case the number of spin flips) in some given time. It is the natural trajectory observable to quantify the amount of motion in the dynamics. A question one can ask is what is the probability of observing the activity K for trajectories ω t which run for a total time t. The probability distribution for K is given by where π(ω t ) is the probability of observing trajectory ω t . For long times this obeys a large deviation (LD) principle P t (K) ≈ e −tϕ(K/t) where ϕ(K/t) is the LD rate function [41]. One can also consider the moment generating function which also obeys a LD principle, Z t (s) ≈ e tθ(s) where θ(s) is the scaled cumulant generating function (SCGF) whose derivatives at s = 0 give the cumulants of K, scaled by time [41]. The SCGF plays the role of the thermodynamical free energy and is related to the LD rate function by a Legendre transform θ(s) = − min k [sk + ϕ(k)] [41]. We can deform the generator given in Eq. (2) by multiplying the off-diagonals by a factor of e −s to give the tilted generator, whose largest eigenvalue is the SCGF θ(s) [41]. It has the associated left and right eigenvectors, l s | W s = θ(s) l s | and W s |r s = θ(s) |r s respectively. As the dynamics obeys detailed balance, we can transform the generator into a Hermitian one by using a similarity transformation independent of s [39]. We first define the diagonal matrix Q with matrix elements n|Q|n = ( which has the ground state H s |ψ s = −θ(s) |ψ s . As was done for the generator, one can write the tilted Hamiltonian in the DW representation The eigenvector |ψ s of H s is related to the left and right eigenvectors of W s by where l s (n) = l s |n and r s (n) = n|r s . Thus studying the LDs reduces to diagonilising Eq. (21) to find θ(s) and |ψ s .

B. Matrix Product States
A matrix product state (MPS) is an ansatz for the vector state of a many-body system [49][50][51] . For a chain of N finite dimensional subsystems (of dimension d), it corresponds to states of the form where i k labels the states of the physical basis for the k-th subsystem and each A k is a rank-3 tensor with dimensions d × D × D, with D the so-called bond dimension. Thus the MPS is described by O(N dD 2 ) parameters. Notice that by increasing D, any arbitrary state can be exactly written in the form (24), although this may require up to D = d N/2 . The bond dimension D limits the entanglement within the state. More precisely, in a MPS with bond dimension D, the entanglement entropy for a subchain L (defined as S E = − Tr ρ L log ρ L where ρ L = Tr N \L |Ψ Ψ| is the subchain reduced density matrix) is upper-bounded by S E ≤ 2 log D, independent of the subchain length. This implies that the MPS satisfies an entanglement area law, which is intimately related to their success at approximating relevant physical states [52]. In particular, ground states of local gapped Hamiltonians, which in one spatial dimension are known to satisfy an area law [53], but also of critical models, can be efficiently approximated by MPS [53,54]. Furthermore, MPS constitute the basis of efficient numerical methods, including the celebrated density matrix renormalization group (DMRG) algorithm [55,56] which we use to approximate the ground state of H s .
The DMRG, originally formulated in [55], can be understood as a variational minimization of energy over the set of MPS. By writing the operator Eq. (21) as a matrix product operator (MPO) [57,58], one can perform a local optimization on a single tensor within the MPS to minimize the energy. We iterate through each tensor, applying local updates until convergence. This variational MPS search (vMPS) is well detailed in many reviews (e.g. [51,59]). For completeness, we give a brief description in the appendix.
When applying the vMPS to study the LD statistics of the XOR-FA model, we use open boundary conditions (OBC) which formally corresponds to setting n 0 = n N +1 = 0, as this allows for the most efficient MPS calculations, with computational cost O(D 3 ). In our problem, the number of DWs defines a global conserved quantity, and we need to find the ground state in a particular sector. Although it is possible to encode this symmetry in the tensors [60][61][62], as also local constraints have [63], we opt here for a simpler approach. Namely, we add an energy penalty to the Hamiltonian to favour the desired sector. Specifically, the penalty is λ(N DW − N DW ) 2 where λ > 0 is some Lagrange multiplier and N DW = N i=0 n i (1 − n i+1 ) + (1 − n i )n i+1 is the operator which counts the number of DWs.

C. Results
As we now show, the MPS ansatz combined with the variational search proves to be very efficient for studying the LDs of the XOR-FA model, just like for other KCMs and exclusion processes [43][44][45]. In this way we are able to achieve results for system sizes superior to traditional methods such as exact diagonalisation or importance sampling.

First-order phase transition in the SCGF
A key property of other KCMs such as the FA or the East model is the presence of a first-order dynamical phase transition [47] in the thermodynamical limit N → ∞, manifested as a singularity in the SCGF θ(s) at s = 0. Consequentially there are two dynamical phases -the active phase for s < 0 and the inactive phase for s > 0. We look for evidence for this transition in the XOR-FA model. Figures 4(a,b) show the SCGFs obtained numerically for system sizes N ∈ [20, 100], in linear and log-log scales, respectively. The upper row of Fig. 4 is for c = 1/2 while the lower row corresponds to c = 0.1. For finite size, near enough s = 0 the SCGF should obey linear response, θ(s) ∼ −sk eq , where k eq is the average activity per unit time in the equilibrium state. For the FA and East models, the equilibrium activity is straightforward to calculate exactly (see e.g. [64]). For the XOR-FA, it is more difficult due to the conservation of DWs. In the Appendix we give a way to compute it to a good approximation. As we see from Fig. 4(b), the SCGF we obtain numerically does obey LR close to s = 0.
Still for s 0, beyond the LR regime the SCGF changes behaviour, notably stops scaling with system size, see Figs. 4(a,b). This change in behaviour becomes even more apparent when one considers the activity per unit time as a function of s, k(s) = −θ (s)/N , Fig. 4(c). We see a sudden drop close to s = 0 that becomes more pronounced with system size, a hallmark of a first-order phase transition. From the point of the numerics, this occurs when where there two smallest energy levels of Eq. (21) cross.
The transition point s c can be estimated from the peak of the susceptibility χ(s) = θ (s), see Fig. 4(d). The peaks become higher and sharper with system size. From the numerics we can make a finite size scaling analysis of the critical point. We find that s c (N ) seems to obey s c (N ) ∝ N −α as shown in Fig. 4(e). For n DW = 1/2, we find that the best fit exponent is α ≈ 2.123 for both c = 1/2 and c = 0.1. Furthermore, for n DW = 1/4 we find that α ≈ 2.056 and α ≈ 2.188 for c = 1/2 and c = 0.1 respectively. In each case α is close to the value 2 expected from a diffusive behaviour of the gap in the spectrum of Eq. (21). It could be that these are subleading corrections to N −2 , see Fig. 4(e).

Spatial structure of active and inactive phases
The singularity in the SCGF represents a phase transition in the trajectories of the dynamics at the level of fluctuations: the behaviour at s < 0 corresponding to dynamics with activity that is larger than the typical (equilibrium) one is fundamentally different from that at s > 0 corresponding to dynamics which is less active than typical. This difference is also manifested in the configurations that are predominantly visited by the different trajectories. That is, active and inactive dynamical phases are associated with very different spatial structures.
We mean the following. At s = 0 there is no tilting in the ensemble of trajectories which is the one given by the original dynamics. It corresponds to typical behaviour. Dynamics is ergodic over configuration space, and over long-times the distribution of configurations that are visited -for some fixed value of the number of DWs -is given by the equilibrium probability Eq. (13). At s = 0 the probability of a trajectory is reweighted by the exponential of its activity, cf. Eq. (19). How often configurations are visited in such reweighted ensembles depends on s, and for long-times the distribution over configurations is given by the leading eigenstate of Eq. (20), or equivalently Eq. (21) for the detailed balanceobeying XOR-FA. We have access to these states, |ψ s , from our MPS numerics.
The easiest way to characterise the spatial structure of the characteristic configurations associated with dynamics tilted by s is to study the average local occupations n i s = ψ s |n i |ψ s [65]. In Figs. 5(a,b) we show these local densities as a function of s for two values of c and n DW = 1/4. For s large and negative we see that that domains becomes maximally sized, that is, DWs become maximally spaced apart, maximising the activity, as expected. In contrast, for s large and positive DWs become localised at either edge of the system. When DWs become minimally separated and clustered together, only the sites next to the last DW are allowed to move and activity becomes subextensive and thus minimal.
We can further quantify the average distance between neighbouring DWs by considering the operator which measures the likelihood of observing two neighbouring DWs at distance d apart. The average distance between neighbouring DWs is then given by (as we have N DW − 1 pairs of neighbouring DWs). In Fig. 5(c) we show d s as a function of s for c = 1/2 and n DW = 1/2. It is evident from the plot that the dynamical transition at s c ≈ 0 is also one where there is a singular change in the distance between DWs in the corresponding characteristic configurations, from maximal distance between DWs at s negative, to minimal at s positive. Figure 5(d) shows the same for c = 0.1 and n DW = 1/4. We see that away from the SEP limit of the XOR-FA, there is even richer spatial structure due to the energetic cost associated with domains. At small s < 0 there is an initial plateau in the growth of the average distance between DWs. This is an extension of the equilibrium behaviour, where domains are randomly sized according to the behaviour described in Sec. III. As we move further into the active phase, we observe another plateau. At this point, the excitation density n s = i n i s (as shown in the inset) has not varied much from the equilibrium value. This leads us to believe that this change in structure is due to the excited domains spreading apart and becoming localised as shown in Fig. 5(b). This maximises the activity without having to grow the excited domains as is energetically favourable for c < 1/2. Following this plateau, there is a steady growth to the maximum d s corresponding to the growth of the one domains, such that every DW is on average equally spaced. This is fur- ther supported by the fact that following this plateau, the excitation density rapidly grows.

Maximally and minimally active limits
While we cannot calculate the SCGF analytically in general, there are limits where the calculation becomes tractable (apart from the obvious case of θ(0) = 0). These are the s → ±∞ limits corresponding to the tilting of the dynamics that maximises and minimises the activity. For the former we can easily calculate the ground state of Eq. (21) via vMPS to obtain the rescaled SCGF θ/N = lim s→−∞ e s θ(s)/(N c(1 − c)) using only a small bond dimension of O(10). The numerical data is shown in Fig. 6(a) for various filling fractions. Note that for s → −∞ the dependence on c is irrelevant, and can be scaled out as in our definition ofθ.
We can extrapolate from the numerical results for finite size to obtain an estimate of lim N →∞θ /N as a function of the density of domain walls n DW . This largesize limit of the SCGF can be obtained analytically. For s → −∞, after scaling out the e −s and c(1 − c) factors, the Hamiltonian Eq. (22) becomes that of the XX model. Being non-interacting, the ground state can be obtained exactly by standard means [66], to give: lim N →∞θ /N = π 2 sin(π n DW ). Figure 6(a) shows the agreement between the numerical extrapolation and the exact result. The structure of the state in the maximally active limit for a system with N = 40 sites and N DW = 20 domain walls is shown in Fig. 6(b). In the limit N → ∞ we would expect that lim s→−∞ n s = 1/2, and lim s→−∞ d s = n DW −1 which are both in excellent agreement with the numerical vMPS data.
For the minimally active limit s → ∞, the Hamiltonian given in Eq. (21) becomes a diagonal one, where again the dependence of c falls out as a pre-factor. Each configuration is an eigenstate and one can easily show that the configuration with the least energy is the one where all the DWs are clustered at the edge of the system (which is doubly degenerate). The SCGF at this limit is given by θ(s → ∞) = −c and the exact structure for N = 40 and N DW = 20 is given in Fig. 6(b) for just one of the degenerate states (in practice the vMPS prefers to converge onto just one to keep the bond dimension minimal). Additionally, the excitation density and the distance between DWs are given by lim s→∞ n s , d s = n DW /2, 1 respectively.

V. CONCLUSIONS
We have studied a new classical stochastic KCM, the one-dimensional XOR-FA model. The kinetic constraint in this model is stronger than that of the standard FA model, as spins can flip if only one nearest neighbour is in the excited state, in contrast to the FA where spins can also flip if more than one neighbour is excited. It is also less constrained than the PXP, or 2-spin facilitated FA, model which requires both neighbours to be simultaneously in the same state. As such it shares features with both these models. The constraint imposes a conservation law, that of the total number of domain walls, breaking configuration into a number of disconnected components that scales with system size. This contrasts to the FA model where all but one configuration are dynamically connected, and is closer to the PXP where configuration space is also disconnected. The distinction with the PXP is that in the XOR-FA the conserved quantity is global, while in the PXP there are many local conserved quantities and configuration space is broken into exponentially many disconnected components. This less severe disconnection of state space makes the dynamics within connected components in the XOR-FA still interesting as there is scope for large dynamical fluctuations (in contrast to the PXP).
The XOR-FA is a "thermal" model, in the sense that it obeys detailed balance with respect to the Boltzmann distribution of a non-interacting binary gas, where each excited spin costs a unit of energy, and subject to the conservation of the number of DWs. In analogy to other KCMs like the FA and East models [4], the XOR-FA has a trivial (up to the conservation law) thermodynamics, but complex relaxation dynamics due to the constraint. The conservation law allows the XOR-FA to be represented in terms of the dynamics of its DWs. This establishes a connection to exclusion processes. At infinite temperature the XOR-FA can be mapped via a duality transforma-tion to the SEP [25,29,30], while at finite temperature this mapping leads to a dynamics of hopping DWs with long-range interactions (in contrast to the original spin representation which is always local).
Like in other KCMs [6,39,47,64], we have shown here that the XOR-FA has a trajectory level phase transition between active and inactive dynamical phases. We have proved this to high accuracy via numerical MPS. This adds to recent work [43][44][45] demonstrating the effectiveness of numerical tensor network methods for studying the full counting statistics of stochastic systems. One of the many advantages of this approach is the direct availability of spatial structural information in the various dynamical phases. For the XOR-FA we find spatial structure as expected from its connection to the SEP [67,68]: a change from repulsion of DWs in the active phase, maximising the possibility of motion, through structureless configurations in the typical equilibrium dynamics, to DW clustering in the inactive phase. Away from the strict SEP limit, these structures remain overall, with further richness due to energetics.
Here we have studied the XOR-FA in one-dimension. It is easy to generalise the model to higher dimensions, once again motivated for example by the problem of atoms interacting strongly in Rydberg states. For example, in the blockade regime, when going from one dimension to a, say, two dimensional square lattice, the PXP model becomes the hard square model [69]. Similarly, in the anti-blockade regime of Rydbergs, the XOR-FA would generalise to a model that is less constrained than that of hard squares, but more constrained than the 2-spin facilitated FA model in two dimensions [4]. From the results here it is safe to speculate that such higher dimensional generalisations of the XOR-FA will also display very rich dynamics. Additionally, the structure of the ground state observed at s > 0 is reminiscent of the localized ground states found in quantum KCMs [9], which has dramatic consequences for the quantum dynamics of the model. It may be interesting to investigate whether an analogous localization can be characterized in the ground state of the XOR-FA model in the inactive regime. this criterion for a run of the algorithm at some bond dimension D, we run it again with an MPS with a higher bond dimension, typically using the previous run as our initial guess.

Canonical solutions using MPS
Consider some general probability vector |ν for the configurations of a chain of N spins. The probability of observing some configuration |n = |n 1 , . . . , n N is then given by p(n) = n|ν . Now suppose we have some ob-servableÔ which is also diagonal,Ô |n = O(n) |n and we want to reweigh the probability that any configuration is observed by e −µO(n) as we have done in Sec. IV with the dynamical activity. The new probability vector is then Furthermore, the partition function can be calculated by taking the inner product with flat state, −| = n n|. That is, This allows us to directly find the expectation value Ô µ in the tilted ensemble by taking the derivative and in particular, we can find the expectation value of this observable with respect to |ν (in the non-tilted ensemble) by evaluating Eq. (32) at µ = 0. For the XOR-FA at equilibrium (where one does not worry about the number of DWs), the probability vector |ν eq is defined by wherec = 1 − c. Alternatively, one can just write the product state This can easily be written as an MPS (cf. Eq. (24)) where each tensor is identical, A i = A ∀i. In particular, the MPS has a bond dimension 1 and is given by A 0 =c and A 1 = c. We now want to calculate the observables studied in Sec. III. In particular, we need to calculate the number of DWs, the number of excitations and the escape rate. The corresponding operators are N DW , N and R respectively. Each are a sum of local observables which act on the site i and at most its nearest neighbours. The local operators for each observable are given by r i = cn i−1ni n i+1 +cn i−1 n i n i+1 + (37) cn i−1nini+1 +cn i−1 n ini+1 , wheren i = 1 − n i . The operator which counts each of these quantities exponentially is e −Ô , wherê Notice that each observable has its own re-weighting factor µ i . This is necessary as we want to calculate each quantity, independently of one another. The exponential operator can then be split into the product of many local operators, This is only possible becauseô i = µ 1 (n DW ) i + µ 2ni + µ 3ri are diagonal operators in the configuration basis, meaning they commute with one another. It follows that the MPS |ν O (µ 1 , µ 2 , µ 3 ) = e −Ô |ν eq can be calculated by applying each e −ôi individually. We define e −ôi |n = e −oi(ni−1,ni,ni+1) |n (40) where one can calculate o i (n i−1 , n i , n i+1 ) for each of the possible inputs using Eqs. (35)(36)(37) and Eq. (38). This is a three-body operator, and can be written exactly in MPO form where the tensor for the operator is exp(−ô i ) ni−1,ni,ni+1 mi−1,mi,mi+1 = δ ni−1,mi−1 δ ni,mi δ ni+1,mi+1 e −oi(ni−1,ni,ni+1) .
The delta functions ensure that the operator acts diagonally. We can break the operator down into three separate tensors, one for each site as show on the right of the diagram. The end tensors (black circles) are the threelegged deltas δ ni−1,α,mi−1 = δ ni−1,α δ ni−1,mi−1 , and the middle (square) tensor is δ ni,mi e −oi(α,ni,β) .
The diagram for the new MPS is ... ...
We then group together the indices (j 1 , j 1 ) and (j 2 , j 2 ) to make this a rank-3 tensor as required.
The tensor is then given by the 4 × 4 matrices By considering the inner product of the MPS with the flat state, one can easily show that it also takes the form For large system sizes, the partition function is dominated by the largest eigenvalue Θ(µ 1 , µ 2 , µ 3 ) of B 0 + B 1 . In particular, it obeys a large deviation principle Z N (µ 1 , µ 2 , µ 3 ) ≈ e N ln Θ(µ1,µ2,µ3) .
Before we can calculate the expectation values at equilibrium, we first need to determine the value of µ 1 = µ DW necessary to tilt the probabilities so that on average we have a DW density of n DW . This can be done using Eq. (32) with the derivative of µ 1 (scaled by N −1 ). Thus the average domain wall density is Once we have determined the required µ DW , we can then calculate the average excitation density and the average escape rate per site at equilibrium for the required DW density again using Eq. (32). The expression for each is It is important to note that here we have not enforced that we have exactly the required number of DWs, but we have it on average. In this way we have moved from the micro-canonical ensemble where we conserve the number of DWs to the canonical one. However, in the limit N → ∞, it converges onto the true expected value. In Sec. III we compare these analytical results to numerical ones for large system sizes. Both show excellent agreement.