Slow dynamics and large deviations in classical stochastic Fredkin chains

The Fredkin spin chain serves as an interesting theoretical example of a quantum Hamiltonian whose ground state exhibits a phase transition between three distinct phases, one of which violates the area law. Here we consider a classical stochastic version of the Fredkin model, which can be thought of as a simple exclusion process subject to additional kinetic constraints, and study its classical stochastic dynamics. The ground state phase transition of the quantum chain implies an equilibrium phase transition in the stochastic problem, whose properties we quantify in terms of numerical matrix product states (MPS). The stochastic model displays slow dynamics, including power law decaying autocorrelation functions and hierarchical relaxation processes due to exponential localization. Like in other kinetically constrained models, the Fredkin chain has a rich structure in its dynamical large deviations - which we compute accurately via numerical MPS - including an active-inactive phase transition, and a hierarchy of trajectory phases connected to particular equilibrium states of the model. We also propose, via its height field representation, a generalization of the Fredkin model to two dimensions in terms of constrained dimer coverings of the honeycomb lattice.


I. INTRODUCTION
The Fredkin spin chain [1,2] is a one-dimensional lattice model with local three-body interactions, whereby hardcore particles can hop to adjacent sites if allowed by constraints involving next to nearest neighbors.This model has been of interest in the quantum many-body community over the last few years for a number of reasons.In its original formulation [1,2], the Fredkin chain can be expressed exactly as an equal superposition of all Dyck paths, i.e. random walk (RW) excursions, with appropriate endpoints, with an entanglement entropy which scales logarithmically in system size, thus violating the area law [1][2][3][4].Furthermore, the model has slow unitary evolution [5][6][7][8] due to dynamical "jamming".With the addition of particular potential energy terms the model features a ground state phase transition between states of bounded and extensive entanglement entropy [9,10].These interesting properties have brought about further studies into generalized Fredkin models [11], including versions which present "quantum scars" [12].
Dynamical constraints, such as those present in the Fredkin model, are responsible more generally for many interesting phenomena in many-body dynamics.A striking example of this are the kinetically constrained models (KCMs) of structural glasses [13,14] -simple lattice models equipped with local dynamical constraints, leading to slow relaxation and dynamical heterogeneity [15,16].Such models can also be considered as systems under closed unitary [17][18][19][20] and open dissipative [21][22][23] quantum dynamics.A recent example of these is the quantum PXP model [24,25] of Rydberg atoms in optical lattices under blockade conditions, which has been shown to exhibit non-thermal eigenstates (often called quantum scars [26]).Another area where dynamical constraints lead to interesting non-equilibrium dynamics is in deter- The top shows the RW representation of the height field, which must always satisfy h > 0. The middle is the corresponding particle representation.The bottom is in terms of Dyck words, where opening "(" must always be matched with a closing ")".
Here we provide such systematic study of both typical dynamics and rare fluctuations.
Classically, the Fredkin model resembles the simple ex-clusion process (or SEP, for reviews see [50,51]).Like the SEP, it describes particles hopping stochastically to neighboring empty lattice sites with at most one particle per site.The key difference is the presence of further local kinetic constraints to motion.These, together with specific boundary conditions, specifically that of an open segment with fixed boundaries, restrict the dimensionality of the state space.For example, for a length N = 2M chain half filled with M particles the dimensionality is the Catalan number rather than the binomial coefficient 2M M .Although the difference in configurational entropy is not extensive, this constrained state space plays an important role in the dynamics, as we explain below.
SEPs and KCMs display interesting dynamical properties which can be studied with large deviation (LD) methods (for reviews see [52][53][54][55]).A central result in the dynamics of these systems is the existence of phase transitions in the space of trajectories, indicated by singularities in the LD functions that quantify the dynamical fluctuations in the long time limit, both in terms of time-integrated currents [56][57][58][59], or dynamical activities [60][61][62][63][64].In the case of the Fredkin model, a preliminary study [49] indicated that it also displays LD transitions.Here we make this finding concrete by studying LD functions using matrix product states.
The paper is organized as follows.In Sec.II we start by defining the model and reviewing its basic properties.We highlight its relationship to Catalan combinatorics and RW excursions [65].In Sec.III we consider the equilibrium states which follow from the properties of the ground state of the quantum problem [1,9].We study the properties of the equilibrium phases in detail by means of numerical DMRG [66].An interesting observation is that there are three distinct equilibrium phases, and a transition between them, despite the fact that this is a one-dimensional system with local dynamical rules.This apparent contradiction with the Landau principle is a consequence of the constrained configuration space of the model.In Sec.IV we study the relaxation dynamics.As in the case of the quantum model [8], the stochastic Fredkin spin chain exhibits slow dynamics.We provide evidence for power law decaying autocorrelations, and for a pattern of hierarchical relaxation when quenched from extremal initial states into the different equilibrium phases.In Sec.V we study the large deviations statistics of dynamical observables by means of numerical MPS.As in other constrained models, the phase transitions at the LD level underpin the slow dynamics and fluctuations seen in typical relaxation trajectories.We reveal the existence of an active-inactive transition, as in other KCMs, but also a hierarchy of trajectory transitions connected to hierarchical relaxation dynamics.In Sec.B we speculate on a possible generalization of the Fredkin model to a two-dimensional setting defined in terms of fully packed dimers on the honeycomb lattice (that is, rhombus coverings of the plane).We give our conclusions in Sec.VI.

II. MODEL
The Fredkin model is defined in terms of particles hopping on a lattice of N sites with binary occupation, n j = 0 (for empty or down) or 1 (for occupied or up) with j = 1, • • • , N .The system evolves under stochastic continuous-time Markov dynamics with generator where σ ± i are Pauli creation and annihilation operators on site i.The factor in curly brackets in each term is the same as the local generator of the asymmetric SEP (ASEP) [50,51], with rates for hops to the left or right given by c and 1−c, respectively.What distinguishes the Fredkin model from the ASEP is the kinetic constraint which means that hopping between sites i and i+1 is not allowed if n i−1 = 0 and n i+2 = 1, see Fig. 1(a) [67].In Eq. (1) we are considering open boundary conditions on a segment [1, N ] with no injection/ejection of particles at the boundaries.The fixed sites at the edges, which are not acted on by the generator, we fix to be n 1 = 1 and n N = 0. Note that at c = 1/2, Eq. ( 1) is equivalent to the quantum Hamiltonian of the original Fredkin model defined in Ref. [1], up to a minus sign and boundary terms.For c ̸ = 1/2 the generator Eq. ( 1) obeys detailed balance despite the asymmetry in the hopping rates [68], meaning that for all values of c we expect to find an equilibrium stationary state of W. Notice that under a similarity transformation (see below) it becomes equivalent to the "deformed" Fredkin model of Ref. [9].
The model discussed here has various symmetries.The most obvious one is the conservation of the total number or particles (or occupied sites): M = i n i .This property is shared with the SEP.The constraint Eq. ( 2) gives rise to a further subdivision of each subspace of fixed M , which is most easily understood by a representation of the allowed moves in terms of matched brackets [1].In this representation, particles and holes correspond to opening and closing parentheses, and the dynamics respects normal matching rules.Thus the move is forbidden because both sides cannot simultaneously be matched configurations [this forbidden transition is shown in Fig. 1 For concreteness, here we will focus on the case of halffilling by fully matched particles and holes i.e.M = N/2 a = b = 0.In this case the accumulated number of particles starting from the left is never smaller than the accumulated number of holes (that is, the sector that is dynamically connected to having all particles to the left and all holes to the right, see below), cf.[1].We call this sector D.
It is convenient to represent a configuration x = n 1:N also in terms of a height field defined as with boundary condition h 0 = 0, and where Z i = 2n i −1.
For all configurations x ∈ D we have h i (x) ≥ 0 for all i.
If we think of the space direction as "time" and a particle (hole) representing a step up (down), then D is the space of all paths that correspond to random walk excursions [65], that is, random paths that return to the origin while never crossing the horizontal axis.(In contrast, for the SEP in the height representation at half-filling, the space of dynamically connected configurations is that of random walk bridges, which are also constrained to return to the origin but can cross the horizontal axis.).Excursions are also known as Dyck paths.An example configuration is shown in Fig. 1(c) with each of the representations.

III. EQUILIBRIUM STATICS
To determine the equilibrium properties of the model we need to find the state |ss⟩ annihilated by Eq. ( 1).Let us consider as an observable the area under the height profile of a configuration x, It is then easy to see that the the stationary state to the dynamics Eq. ( 1) is given by [9] with N c a c-dependent normalization constant to make ⟨− |ss⟩ = 1, where ⟨−| = x∈D ⟨x| is called the flat state.
The connection to RW excursions means that this probability distribution is related to the Airy function [65,69].The properties of the stationary state at arbitrary c can also be understood from the properties of the ground state of the corresponding quantum model [1,9].That is, if we perform a similarity transformation of Eq. (1) (cf. the ASEP with the same boundary conditions, e.g.[70]) where P 1/2 is the diagonal matrix of the square root of configuration probabilities, we get the Hamiltonian whose ground state is |ψ⟩ = P 1/2 x∈D |x⟩.The transformation to a Hermitian form shows that, despite the asymmetric hopping when c ̸ = 1/2, the Fredkin model obeys detailed balance and consequently the stationary state Eq. ( 7) is an equilibrium one.
The properties of the ground state of Eq. ( 10) are well understood from previous studies [10,11].Here we restate them from the point of view of the equilibrium state of the stochastic model, using matrix product states (MPS, see reviews e.g.Refs.[71][72][73]).
From the connection to RW excursions at c = 1/2 the equilibrium state |ss⟩ can be written exactly as an MPS where n are site dependent tensors, and (i| and |f ) appropriate boundary vectors in the auxiliary (or bond) space of the MPS (we use rounded brackets to distinguish them from vectors in configuration space).
Consider first the slightly simpler problem of the symmetric SEP (SSEP), whose generator is given by an operator like Eq. ( 1) but without a constraint, f i = 1.If we consider the same boundary conditions as for the Fredkin model, but with extra terms in Eq. ( 1) that allow particle hops between sites j = 1, 2 and N − 1, N (no longer prevented in the absence of a constraint), then the SSEP configurations at half-filling are those of RW bridges.If the height field h j describes the position of the RW after step j, the exact transition probabilities at step j for generating bridges of N steps are: for |h| ≤ N +1−j, or zero otherwise.(These are obtained from the naive symmetric RW transition probabilities via a Doob transform, see e.g.[74].)The equilibrium MPS for the SSEP is then given by the (2N + 1) × (2N + 1) matrices with boundaries (i| = (0| and |f ) = |0).It is easy to see that the matrices above satisfy B (j),SSEP 0 = 0 for all j, which means that the MPS Eq. ( 11) with tensors Eqs. ( 13) and ( 14) is annihilated by the SSEP generator.
The construction for the equilibrium state of the stochastic Fredkin chain at c = 1/2 is similar, but the relevant paths are RW excursions.In this case the "Doob" transition probabilities that guarantee an excursion are (cf.e.g.[74]) for 0 ≤ h ≤ N + 1 − j, or zero otherwise.The corresponding matrices have now bond dimension N + 1 and read with the same boundary vectors (i| = (0| and |f ) = |0).
The relevant relations in this case are B = 0 for all j.Given these, one can show that the MPS Eq. ( 11) with tensors Eqs. ( 16) and ( 17) is annihilated by the Fredkin generator Eq. ( 1).In fact, the MPS is annihilated by every local term in the spatial sum that defines Eq. ( 1), so that W can be said to be a parent generator (cf.parent Hamiltonian [73]) of the MPS Eq. (11).
Note that from the definition of the tensors B n above in terms of transition probabilities, the MPS is in "right canonical" form, and Eq. ( 11) therefore satisfies ⟨− |ss⟩ = 1.Away from c = 1/2 we can also write Eq. ( 7) as an MPS if we reweigh the coefficients in Eqs. ( 16) and ( 17) as These reweighed coefficients are not transition probabilities in the height (they do not add up to one), meaning that the resulting MPS is not in canonical form.Finding the normalization N c in this case is non-trivial.

B. Equilibrium phase diagram from numerical MPS
To overcome the difficulty above, in order to study the equilibrium properties for all c we resort to numerical MPS approximations.This we implement with the ITensor library [75], and make use of the density matrix renormalization group (DMRG) [66,76,77] to find the leading eigenvector of Eq. (10).We employ an adaptive bond dimension, which is at most D = 2000 with a truncation cutoff error ϵ = 10 −12 .Furthermore, we exploit the U(1) symmetry which conserves the number of particles and initialize the MPS with a product state which lies in D. We then carefully check the relevant observables to ensure they satisfy the properties associated with D, such as a positive height field.
By looking at various observables at stationarity, it becomes clear that there are three distinct equilibrium phases in the Fredkin model: (i) c < 0.5, (ii) c = 0.5 and (iii) c > 0.5.We denote the expectation value of an observable O with respect to the equilibrium state as ⟨O⟩, with The appropriate order parameter to characterize the equilibrium phases is the average area ⟨A⟩.In Fig. 2(a) we show ⟨A⟩ as a function of c for a range of system sizes N ∈ [20,400].For c < 1/2 the area becomes minimal, while for c > 1/2, the area is maximal.If we consider the area as a function of system size N we find that ⟨A⟩ grows as a power law ⟨N ⟩ ∼ N −β , as shown in Fig. 2(b).This reveals three distinct behaviors: the exponent β takes the values β = 1, 3/2 and 2 for c < 1/2, c = 1/2 and c > 1/2, respectively [9].
For each phase, we show the average of the spatial occupation profile, ⟨n i ⟩, and the average height field, ⟨h i ⟩, in Figs.2(c,d), respectively.For c < 1/2, the particles take an anti-ferromagnetic arrangement, Fig. 2(c) (red circles), thus minimizing the height and therefore the area, Fig. 2(d) (red circles).We sometimes refer to this as the flat phase (in analogy with interacting dimers [78,79]).
At c = 1/2, all configurations occur with equal probability, cf.Eq. (7).In terms of the RW representation of configurations this corresponds to the set of RW excursions.The average occupation, Fig. 2(c) (blue squares) interpolates between 1 and 0, and in the thermodynamic limit, N → ∞, the average occupation density in the bulk is 1/2 [11].In turn, the average height field takes a semi-circular form, Fig. 2(d) (blue squares).Note that this a phase of large fluctuations and this average height field is not representative of typical sample profiles.This is in contrast to the other two phases which are exponentially dominated by extremal area configurations, cf.Eq. (7).
For c > 1/2, the particles (holes) localize to the left (right) edge of the system [10], with a sharp change in average occupation, Fig. 2(c) (green triangles), and with an average height profile in the shape of a tent (with a rounded top, a finite residue of the fluctuations of the c = 1/2 phase), Fig. 2(d) (green triangles).This behavior is similar of that seen in the ASEP in an open segment with fixed boundaries [70].Simple arguments (see Appendix A) give the profile [10] ⟨n j ⟩ = 1 exp ((j − N/2)/λ) + 1 (19) with an inverse localization length λ, We sometimes refer to the c > 1/2 phase as the tilted phase (also in analogy with interacting dimers [78,79]).Note that this shares no connection with the "tilted generator" later introduced in Section V.
An observable which will be of importance later is the dynamical activity ⟨k⟩, which measures the average number of configuration changes per unit time in stochastic trajectories [61,62,80].At equilibrium, it can be measured as the average escape rate, ⟨k⟩ = ⟨−|R|ss⟩, where R is the diagonal part of Eq. ( 1).We show this in Fig. 2(e) as a function of c for various system sizes N ∈ [20,400].It is immediately clear that the dynamical activity scales with system size (up to small finite size effects) for c ≤ 1/2 where occupation is spread out in equilibrium, cf.Fig. 2(c), leading to less constrained and therefore more dynamics throughout the entire system.Conversely, the activity for c > 1/2 is sub-extensive in system size as expected due to the much more inactive conditions given to the localization of the equilibrium state, cf.Fig. 2(c): motion is limited to the center of the lattice (the tip of the "tent") where particle hops are not restricted by exclusion.By fitting the activity with a linear form ⟨k⟩ = a + bN (for c ≤ 1/2), one can extrapolate to infinite size to determine ⟨k⟩ /N in the thermodynamic limit.We show this as the black dashed line, peaking around c ≈ 0.36.Notice that for c > 1/2, the activity goes as O(1) and is surpressed by the scaling in N .The differences in the "active" (c ≤ 1/2) and "inactive" (c > 1/2) dynamics are directly related to the dynamical large deviations in Sec.IV.

C. Localization of the tilted phase
The equilibrium state for c > 1/2 is exponentially dominated by maximal area configurations, that is, configurations in which particles cluster towards the left edge of the system, and holes cluster at the right edge. Figure 3(a) shows the average occupation profile ⟨n i ⟩ for various c > 1/2: for sites beyond the halfway point, i > N/2, we observe an exponential decay of the average occupation, ⟨n i ⟩ ∼ e −i/λ .[Note that the same occurs for the density of holes, 1 − ⟨n N +1−i ⟩, coming from the right, due to fact the generator Eq. ( 1) is invariant under This localization can be further characterized by the density of domain walls (DWs) This is shown in Fig. 3(b): the DW density is close to 1 at the centre of the lattice, and decays exponentially when moving away from it in both directions, Notice that the localization is consistent for increasing system size.As we discuss further in the next sections, exponential localization of DWs at the centre of the lattice has important implications for the dynamics in the tilted phase: particle hops can only occur when there are domain walls, and thus activity is exponentially The line shows the result from the theory, Eq. ( 20), and the blue circles and red crosses the numerically extracted lengths from the occupation and DW profiles, respectively.The numerical data is from DMRG with N = 100.
suppressed away from to midpoint, and is sub-extensive in system size, cf.Fig. 2(e).
The localization length λ decreases with increasing c.We show this in Fig. 3(c) for both particle and DW densities.The agreement with the theoretical prediction Eq. ( 20) is excellent.The numerically extracted lengths here are from DMRG with N = 100.For smaller system sizes, the localization length becomes comparable to system size for c ≈ 1/2, and one might expect to see small deviations from the theoretical prediction.

IV. TYPICAL DYNAMICS A. Dynamics in equilibrium
Figure 4(a) shows representative trajectories in the stationary dynamics of each of the three equilibrium phases (with the initial states sampled from equilibrium).The largest fluctuations occur for c = 0.5.Dynamics in equilibrium can be quantified through the (density) autocorrelation function which provides a measure of the memory of a initial configuration after time t in an equilibrium trajectory.We show C(t) for the three equilibrium phases in Figs.4(b,c).It is apparent from Fig. 4(b) that for c = 1/2 the autocorrelation decays asymptotically as a power law, with a numerically extracted exponent of just under a half.This power law decay can also be observed for intermediate times at c > 1/2 [corresponding to fluctuations of the top of the "tent", cf.Fig. 2(d)], with this intermediate regime becoming longer as c gets closer to 1/2.While at short times decay is exponential, see Fig. 4(c), for long times relaxation is stretched exponential in both the flat and tilted phases.These are signatures of slow dynamics.We can extract a timescale for relaxation of correlations from C(t) from its integral, This is shown in Fig. 4(d) for a range of c.This equilibrium timescale spikes at c = 1/2, as expected from the slow law decay of C(t).Notice that the spike is less sharp for smaller system sizes due to the finite size effects from the boundaries.

B. Relaxation towards equilibrium
Also of significance is the relaxation towards the equilibrium state when starting from non-equilibrium conditions.We explore this behavior by considering dynamics following from an initial state of minimal area, x 0 = 1010 • • • 1010, corresponding to a quench from deep in the flat phase (c 0 ≳ 0) to finite c > 0. When c < 1/2, equilibrium is achieved quickly as the initial state is not far from typical states in the flat phase.Interesting nonequilibrium dynamics occurs when quenching to c = 1/2 or to the tilted phase at c > 1/2.In Fig. 4(e) we show two relaxation trajectories, one for c = 1/2 (left) and another for c = 0.6 (right) [81].The system size is N = 100 and the overall time of trajectories t = 10 5 (note that the time axis is shown on a logarithmic scale).For the case of a quench to c = 1/2, after a slow early regime, equilibrium is reached in reasonable times.
For a quench to c > 1/2, we observe a slow hierarchical relaxation, with a progressive coarsening of clusters of particles and holes.The target state is a tilted one, cf.Figs.2(d), and in the height representation this hierarchical process is the merging of smaller tents in the profile into larger ones.Due to the constraint, Eq. ( 2), local configurations of • • • 0011 • • • , corresponding to troughs in the height field, are locally trapped, and require particles from the right edge of clusters to diffuse away to allow clusters to merge.This process is exponentially scarce in the separation distance, as occupations are exponentially localized, cf.Sec.III C.
The time to complete each stage of relaxation in the tilted phases grows exponentially with the stage.This hierarchy is evident in the growth of the average area normalized by system size, ⟨A(t)⟩ /N , as shown in Fig. 4(f), where we see the area increasing logarithmically in time.This reveals the hierarchical nature of the relaxation process: while smaller systems may have relaxed to equilibrium, larger systems require the merging of larger clusters, and so the growth of the area continues.

V. DYNAMICAL LARGE DEVIATIONS
We now study the statistical properties of the stochastic trajectories ω t = x 0:t of the Fredkin model, in particular the LD statistics of dynamical observables.
If K(ω t ) is a trajectory observable, the probability of it taking a value K is where π(ω t ) is the probability of the trajectory ω t being realized under the stochastic dynamics.For a dynamical observable K that is time-extensive, in the long-time limit, the probability of K obeys a LD principle [52-55] where the function φ(k) is called the LD rate function.
The above asymptotic equality holds as long as the spectral gap is non-vanishing (which it is in the Fredkin model for finite size N [2]).A LD principle also holds for the moment generating function (MGF) where θ(s) is the scaled cumulant generating function (SCGF) whose derivatives at s = 0 give the cumulants of K, scaled by time [52].In analogy with what occurs in equilibrium thermodynamics, the rate function and SCGF are related by a Legendre transform We consider as observable K the dynamical activity.Its SCGF is given by largest eigenvalue of the tilted generator, W s , which for the Fredkin model reads , with s being counting field.As W s is in general non-Hermitian, the leading eigenvalue θ(s) has right and left eigenvectors |r s ⟩ and ⟨l s |.We can write the generator in a Hermitian form with the same similarity transformation as before, Eq. ( 8), with ground state H s |ψ s ⟩ = −θ(s) |ψ s ⟩, related to the leading eigenvectors of W s by where l s (x) = ⟨l s |x⟩ and r s (x) = ⟨x|r s ⟩.
A. Active-inactive trajectory transitions at c ≤ 1/2 From the ground state of Eq. ( 28) we can study statistical properties of the trajectory ensemble of the Fredkin model for long-time trajectories.We do this by means of numerical tensor networks along the lines of similar recent work in KCMs [82][83][84][85][86][87]. Figure 5 shows the LD statistics obtained numerically.The top row gives this for the flat phase at c = 0.4, and the middle row for the c = 1/2 phase.These results are for system sizes in the range N ∈ [20, 400] obtained using DMRG.
Column (a) of Fig. 5 shows the SCGF as a function of s = 0 for a range of sizes.For small s ≳ 0, the SCGF follows linear response (LR), θ(s) ≈ −sk(0), where k(s) = −θ ′ (s) is the average dynamical activity in the ensemble tilted by s, shown in column (b).The LR prediction is shown by the dashed black line for N → ∞, calculated by fitting the dynamical activity for finite sizes at s = 0 with a power law and extrapolating.Notice that at some s c (N ) > 0, which becomes smaller for increasing N , the behavior deviates from LR to one which no longer scales with N (this is most apparent for c < 1/2).The step in the average activity, Fig. 5(b, top and centre), indicates a phase transition between dynamical phases of high and low activity.The change in activity tends to FIG. 6. Structural properties of the LDs.We show observables for each equilibrium phase with (a) c = 0.2, (b) c = 0.5 and (c) c = 0.9.The top row shows the average occupations ⟨ni⟩ s for site i (with differing system sizes and ranges of s).The middle row shows the area scaled by system size ⟨A⟩ s /N for s < 0. Finally, we show the area scaled by system size squared ⟨A⟩ s /N 2 for s > 0 in the bottom row.
a discontinuity with increasing size, indicative of a firstorder transition.
The point s c (N ) at which the crossover occurs at finite size can be estimated from the peak in the corresponding dynamical susceptibility, χ(s) = θ ′′ (s), shown in column (c) of Fig. 5.As the system size is increased, the crossover point shifts towards s = 0 and becomes sharper.The change in dynamics can be seen in the broadening of the LD rate function φ(k) around the equilibrium average, shown in column (d) of Fig. 5.The rate functions show the characteristic Maxwell construct of a first-order transition between two phases, and active one with large k and an inactive one with vanishing k.Note that while the transition in activity looks less sharp for c = 0.5, we expect to recover the usual first order behaviour for increasing system sizes as seen by the broadening of the rate function.
For c ≤ 1/2, the location of the crossover can be fit by a power law s c (N ) ∼ N −α .The upper panel of column (e) in Fig. 5 shows this for c = 0.4 and c = 1/2.The lower panel of column (e) shows the dynamical exponents α as a function of c.When c is far from 1/2, we have approximately α ≈ 1.2.The exponent increases quickly as we approach c = 1/2, to around α ≈ 2.5, a value similar to that found in other exclusion processes [85].It could be that for values close to (but not equal to) c = 1/2, the measured exponent would be lower if we accounted for larger system sizes.

B. Dynamical phases for c > 1/2
Obtaining accurate estimates of θ(s) for c > 1/2 at large system sizes is difficult due to a proliferation of dynamical phases.In particular, it is hard for DMRG to converge to the correct phase due to a large density of states.For this reason, for c > 1/2 we limit our studies to system sizes N = 6, 12, 18 with large c = 0.9 ≫ 1/2 which allows us to effectively study the hierarchy of dynamical phases using exact diagonalization (ED) [88].The bottom row of Fig. 5 shows these results.
Since the typical dynamics (s = 0) of the tilted c > 1/2 phase is itself inactive, cf.Fig. 2(e), we expect transitions to the active phase to occur at s < 0 for finite size systems.In fact, column (b) of Fig. 5 shows several points where the behavior of the SCGF changes.The number of these points seems to increase with system size.In each case, this change in behavior corresponds to transitions in the dynamics.At each of these points we see a sharp drop in the activity, this becoming sharper with increasing N .The values at which the activity plateaus are multiple integers of the equilibrium activity, k(s = 0), and are shown by the black dashed lines.With the limited range of sizes accessible via ED it is not possible to do a finite-size scaling analysis as we did for c < 1/2.From the systems studied we observe that the first away from equilibrium inactive behavior happens at increasing s (that is, getting closer to 0) for increasing N , which shows in the flattening of the rate function, see bottom panel in Fig. 5(d).

C. Structural properties of the dynamical large deviations
The difference in the behavior of the various dynamical phases also manifests in the structural properties of the configurations visited by the trajectories.The eigenvector |ψ s ⟩ obtained from either DMRG and ED contains the probability amplitudes for each configuration, making it easy to calculate averages of configuration observables O(x) in the tilted ensemble [89], In Fig. 6 we show the local occupations ⟨n i ⟩ s (top panels), and the average area ⟨A⟩ s (middle and bottom panels), for (a) c = 0.2, (b) c = 1/2 and (c) c = 0.9.It is clear that the limit of large activity (s < 0 with |s| large) particles spread out in order to maximize the activity.This is evident by the average area ⟨A⟩ s , which scales linearly with system size N , resembling the structures associated with the equilibrium flat phase for c < 1/2.Thus, the active phase for all values of c is also a structurally flat one.Conversely, in the inactive limit for all values of c (large s > 0) particles cluster at the left edge of the system and maximize the area, which scales as N 2 .Thus, irrespective of the equilibrium static phase, the inactive phase of the dynamics is structurally "tilted".
Interestingly, we observe very sharp transitions for c ̸ = 1/2 even at smaller sizes -this is unusual when compared to other constrained models [82,85].This could be due to the sharp transition in activity at equilibrium, cf.Fig. 2. Indeed, for c > 1/2 we notice sharp structural changes at the location of the sharp points of Fig. 5.It is clear that the corresponding structures are related to the assembly of excited at equilibrium (s = 0) obtained by joining multiple ground states of smaller system sizes (compare to what occurs in the excited states of the quantum East model [19]).Of course this makes sense, as despite the scarcity of the configurations associated with these states, they have large lifetimes (as discussed in Sec.III) with impactful consequences on the relaxation behavior.

D. Entanglement entropy
We now consider the bipartite von Neumann entanglement entropy of the MPS approximations to Eq. (28).We partition the system into two subsystems A and B, which denote the spins i ∈ [1, N/2] and i ∈ [N/2 + 1, N ] respectively.The bipartite entanglement entropy between the two partitions is then calculated by where ρ A = Tr B [ρ] denotes the reduced density matrix for A, and ρ = |ψ s ⟩ ⟨ψ s | is the density matrix for the full system.The Hamiltonian Eq. ( 28) exhibits a ground state phase transition in the bipartite entanglement entropy for s = 0.In particular, it scales as S E (0) ∼ log N for c = 1/2 and S E (0) ∼ 1 for c ̸ = 1/2 [1][2][3][4].We now extend this analysis to s ̸ = 0. Figure 7 shows the entanglement entropy for increasing system size N ∈ [20,200] for c = 0.4 (top) and c = 0.5 (bottom) and a range of s.Notice that for c < 1/2, the entropy obeys an area law for all s (although we observe spikes around the transition from active to inactive dynamics).For c = 1/2, we observe for large magnitude s < 0 the states clearly also obey an area law.As s approaches s = 0, the entanglement entropy appears to grow significantly towards S E (0), and looks to scale logarithmically.It is important to note however we only show a small range of system sizes, and it is most likely that for some fixed s < 0, the entanglement entropy will be bounded as N → ∞, and thus obeys an area law.This can be seen by the "branching" behaviour seen in Fig. 7.An important consequence is that for large enough N , one is able to construct a state with arbitrarily high entropy by tuning the value of s towards s = 0.For the inactive phase s > 0, ψ s clearly also obey an area law, again with the entanglement entropy spiking as s approach s = 0.

E. Limits of maximal and minimal activity
The limit of maximum activity is that at s → −∞.In this limit, the diagonal parts of W s (and H s ) are suppressed and only the off-diagonals are left.Notice that for H s , the dependence on c falls out as a prefactor.As the tilting in W s grows exponentially with −s for negative s, we rescale the SCGF as when taking the limit.The (rescaled) eigenvalue θ coincides with the (similarly rescaled) dynamical activity.We show this in Fig. 8 The average area ⟨A⟩ −∞ , see Fig. 8(a) takes an almost constant value, with small fluctuations for small system sizes, lim Notice that the area scales linearly with system size, and is similar to the equilibrium states found for c < 1/2.This is further seen from the average occupations ⟨n i ⟩ −∞ , see Fig. 8(b), showing the antiferromagnetic pattern of the flat equilibrium phase.The opposite limit of s → ∞ gives the most inactive state.In this limit only the diagonal escape rate part of Eq. ( 27) (or Eq. ( 28)) remains and each configuration x ∈ D is an eigenstate.The configurations with the smallest escape rates dominate.Depending on c and N , this is either the maximal area (i.e., fully tilted) configuration, 1111 • • • 0000, which has escape rate R = 2(1 − c), or the minimal area configuration 1010 • • • 1010, which has escape rate R = c(N − 2).The latter dominates if N > 2c −1 , and the former dominates if N < 2c −1 (with degeneracy at N = 2c −1 ).

VI. CONCLUSIONS
Here we have provided a detailed study of the statics and dynamics of the stochastic Fredkin model.Despite being one-dimensional and having local transition rules, this model displays phase transitions between three distinct equilibrium phases.This is a consequence of the constraints in the dynamics which restrict the state space to that of random walk excursions, with these static transitions controlled by the asymmetry in the particle hopping rates.Two of these phases are ordered, one being flat and another one tilted (in terms of the height field representation), with an intermediate disordered and fluctuating phase.This phase behavior is in some ways reminiscent of interacting two-dimensional dimer coverings [78,79,90].
The constraints in the local transitions of the Fredkin model lead to a rich dynamics, both in equilibrium and in the relaxation after a quench.This richness can be seen as a consequence of a non-trivial phase structure of the ensemble of stochastic trajectories.Using numerical matrix product states with DMRG, we compute the large deviations of the dynamical activity and show the existence of active-inactive space-time phase transitions, something that is also observed in other KCMs.The overall picture is one where the static phases extend into dynamical ones, with the flat phase being also a dynamical active phase, and the tilted phase a dynamical inactive one, with first order transitions between them.
There are many possible continuations of the work here.One is to go beyond one-dimension.As an initial step, in Appendix B we propose a two-dimensional generalization of the Fredkin model: by focusing on the fact that Fredkin configurations are random walk excursions, we proposed a two-dimensional model in terms of packed dimers on the honeycomb lattice with constraints in the the dynamics which enforce configurations to be "excursion surfaces".It will be interesting to study this and similar stochastic models in future work.Another interesting area of exploration would be to study the Hamiltonian Eq. ( 28) under unitary dynamics, in analogy with recent work that studied other quantum KCMs.As occurs with the quantum East model [19], we expect the constraints in Fredkin models to provide mechanisms for localization and non-thermal eigenstates.We hope to report on this in the near future.surfaces in two dimensions.In a configuration with an equal amount of the three kind of tiles the height field is pinned at zero at the boundaries (for example in three sites at angles of 2π/3 within a hexagonal region).This is a two-dimensional version of the one-dimensional height field from a lattice of particles and holes at half filling which is bound to return to the origin.
In the one-dimensional case the elementary local move that preserves the filling fraction is to exchange a particle with an adjacent hole.The analogous move for a rhombus tiling is shown in Fig. 9(d) and corresponds to rotating a triplet of tiles forming an elementary hexagon.This move changes the height of the central site by ∆h = ±3.In order to prevent the height field from becoming negative, which is the defining property of the dynamics of the Fredkin model, transitions like those of Fig. 9(d) have to be constrained, cf.Fig. 1(a).Figures 9(e-g) show the corresponding allowed transition in the two-dimensional case: the exchange of tiles is only possible if either of the extra green/blue/red tile as in arrangement (e/f/g), respectively, is present, and not allowed otherwise.This constraint implies that in the transition the height of the site at the centre of the hexagon cannot go below that of the site indicated by a circle.With this dynamical rules it is guaranteed that the height field of the dimer/rhombus arrangement never becomes negative at any point, a two dimensional version of the RW excursions that define the configurations of the Fredkin model.Furthermore, giving different rates to the forward and backward moves in Figs.9(e-g) should lead to flat and tilted phases weighted by the volume under the surface.

FIG. 1 .
FIG. 1. Fredkin spin chains.(a) The local stochastic transition rates for neighboring occupied and unoccupied sites, given by all choices of their neighbors.The fourth transition is not allowed.(b) The disallowed configuration change in the height representation.The "troughs" (• • • 0011 • • • ) are locally immobile.(c) An example configuration in the chosen symmetry sector.The top shows the RW representation of the height field, which must always satisfy h > 0. The middle is the corresponding particle representation.The bottom is in terms of Dyck words, where opening "(" must always be matched with a closing ")".

FIG. 2 .
FIG. 2. Equilibrium properties of the Fredkin model.(a) The average area (scaled by maximum area) ⟨A⟩ /Am as a function of c for various systems sizes N ∈ [20, 400].The dashed line shows the extrapolated value for N → ∞.(b) The average area (symbols) for c = 0.4 (red / dark grey), c = 0.5 (blue / medium grey) and c = 0.6 (green / light grey).The lines show the power laws ⟨N ⟩ ∼ N −β with β = 1, 1.5, 2 respectively.(c, d) The spin occupation ⟨n⟩ i and height profiles ⟨hi⟩ for each equilibrium phase with a system size N = 60.(e) The average dynamical activity (per unit time and system size) as a function of c for various systems sizes N ∈ [20, 400].The dashed line shows the extrapolated value for N → ∞.All results are calculated using numerical DMRG.

FIG. 3 .
FIG. 3. Localization in the Fredkin chain.(a) The occupation profile ⟨ni⟩ of the steady state for c > 0.5 and N = 20.The occupations exhibit an exponential decay for i > N/2.(b) The average domain wall occupations ⟨n DW N/2−i ⟩ for c = 0.75 and N = 20, 40, 60.We see the same exponential decay of domain wall density as we move away from the centre of the lattice.(c) The localization length λ as a function of c.The line shows the result from the theory, Eq. (20), and the blue circles and red crosses the numerically extracted lengths from the occupation and DW profiles, respectively.The numerical data is from DMRG with N = 100.

FIG. 4 .
FIG. 4. Stochastic trajectories and dynamics.(a) Representative trajectories with initial states sampled from equilibrium for c = 0.4 (top), c = 0.5 (middle) and, c = 0.6 (bottom) respectively, for system size N = 100 and time t = 10 3 .(b) The autocorrelation functions Eq. (22) for each of the three distinct equilibrium phases.At large times, the auto-correlator for c = 0.5 decays as the power law t −0.464 (from size N = 100).(c) The same autocorrelation functions plotted on a double-log ordinate scale.For small times they show exponential decay in the three phases.For large times they take a stretched-exponential form for c < 1/2 and c > 1/2.(d) The numerically estimated timescales Eq. (23) as a function of c (for N = 40, 100 and 400).(e) Example trajectories after a quench from the initial state 1010 • • • 1010 for c = 0.5 (left) and c = 0.6 (on a logarithmic time scale).The former relaxes to equilibrium quickly, whilst the latter shows hierarchical relaxation (both panels for N = 100 and t = 10 5 ).(f) The area (scaled by system size) ⟨A⟩ /N after the same quench, for various system sizes N ∈ [20, 100] and c = 0.6.The dashed line shows log t.All data is obtained using continuous-time Monte Carlo.

FIG. 5 .
FIG. 5. Finite size scaling of dynamical LD transitions.The dynamical LD statistics for each equilibrium phase.The top (middle) row of (a-d) shows c = 0.4 (c = 0.5) with N ∈ [20, 400] obtained via DMRG and the bottom row shows c = 0.9 obtained through ED.(a) The SCGF θ(s) as a function of s.The upper and middle panels are scaled by system size, with dotted lines showing the value for s → ∞ and the dashed line showing the linear response prediction (see the main text).(b) The average dynamic activity k(s) as a function of s.The top and bottom panels are shown on log-log scales, whereas the middle one is shown in linear scale.The dashed lines in the bottom panel correspond to integer multiples of k(0).(c) The dynamical susceptibility χ(s) = θ ′′ (s) as a function of s.(d) The LD rate function scaled by system size φ(k)/N as a function of activity k.The black dashed lines show a Poisson distribution with mean k(0)/N in the thermodynamic limit N → ∞, extrapolated from finite-size DMRG data.(e) We estimate the critical point as a function of system size from the peaks of the dynamical susceptibility for c = 0.4, 0.5 in the top panel.The dashed lines shows a fitted power law sc ∼ N −α , with the bottom panel showing the obtained α for various c.

FIG. 8 .
FIG. 8. Extreme active limit.(a) The rescaled SCGF θ/N (top) and the area ⟨A⟩ −∞ /N (bottom) as functions of N measured via DMRG.We fit the SCGF as a + bN −1 allowing us to extrapolate the value in the thermodynamic limit N → ∞ (see main text), shown by the dashed line.The value for the area quickly settles with increasing system size, indicated by the dashed line.(b) The occupation profiles ⟨ni⟩ −∞ for N = 40.
(a) as a function of N ∈ [10, 400] (circles, shown divided by N ), and fit it with the function of aN + b (blue dashed line, shown divided by N ).By extrapolating to infinity, we find that lim N →∞ θ/N ≈ 0.691.

FIG. 9 .
FIG. 9. Two-dimensional generalization of the Fredkin model.(a) Dimer covering of the honeycomb lattice.(b) Equivalent rhombus tiling of the plane.(c) Definition of the height field: given a tiling, moving along the edges of the rhombi the height increases of decreases by one unit as shown.For example, the path in (b) shows that the start and end points have a height difference of +2.(d) The elementary local moves that preserve the "perfect tiling" character (i.e., no tiling defects, or no monomers in the dimer representation) are rotations of a triplet of tiles forming an elementary hexagon.These are the dimer/tiler equivalent of the particle-hole exchange in the SEP.These moves change the height field of the central site by three units.(e-g) Constrained moves: requiring the presence of the extra tile guarantees that the height of the central site (indicated by the filled circle) never goes below the lowest height of the arrangement (indicated by the open circle).These are the two-dimensional equivalents of the allowed moves in the Fredkin chain, see Fig. 1(a).