Arresting dynamics in hardcore spin models

We study the dynamics of hardcore spin models on the square and triangular lattice, constructed by analogy to hard spheres, where the translational degrees of freedom of the spheres are replaced by orientational degrees of freedom of spins on a lattice and the packing fraction as a control parameter is replaced by an exclusion angle. In equilibrium, models on both lattices exhibit a Kosterlitz-Thouless transition at an exclusion angle $\Delta_{\rm KT}$. We devise compression protocols for hardcore spins and find that {\it any} protocol that changes the exclusion angle nonadiabatically, if endowed with only local dynamics, fails to compress random initial states beyond an angle $\Delta_{\rm J}>\Delta_{\rm KT}$. This coincides with a doubly algebraic divergence of the relaxation time of compressed states towards equilibrium. We identify a remarkably simple mechanism underpinning this divergent timescale: topological defects involved in the phase ordering kinetics of the system become incompatible with the hardcore spin constraint, leading to a vanishing defect mobility as $\Delta\rightarrow\Delta_{\rm J}$.

Introduction. Within the realm of condensed-matter physics, examples abound of continuum systems whose critical phenomena can be well captured by simplified lattice models. For instance, the liquid-gas transition can, remarkably, be described using the lattice gas model, which then maps onto the Ising model [1]. Similarly, the Edwards-Anderson model and related spin-glass Hamiltonians provide a fruitful avenue toward understanding random impurities in magnetic alloys [2].
In this letter, we take a lattice approach towards what might be called vitrifaction, i.e. the emergence of slow dynamics which has played an important role under the headings of glassiness, freezing, jamming and the like. To do so, we devise a particularly simple yet versatile and tractable family of models of what we have termed hardcore spins. These insulating (non-itinerant) spin models are constructed by analogy to hard spheres/disks, which are the common idealized model systems for the phenomenon of jamming. Much progress has been made to understand the random jamming transition of hard spheres both numerically [3,4] as well as in infinite dimensions [5]. In particular, it is possible to define random jamming without the need to invoke any particular dynamics, the hallmark being a jump of the contact number z as a function of packing fraction at "point J" Φ J .
In contrast, our work takes a different tack: addressing the broad issue of dynamical arrest of an athermal system-one in which thermal fluctuations make a negligible contribution to its dynamics-far from equilibrium, we study whether, and by what mechanism, this phenomenon can arise in a lattice model. Concretely, we address its relation to the physics of phase ordering kinetics, i.e. the extent to which an underlying, but under a given dynamics inaccessible, ordered state and its concomitant topological defects play a fundamental role, an idea that dates back several decades in the glass literature [6][7][8][9][10].
The hardcore spin models are defined by a local constraint, whereby no two neighboring hardcore spins are allowed to enclose an angle smaller than their exclusion angle [see Fig. 1 (a)]. This can be viewed as an orientational, lattice analog to the non-local constraint on the translational degrees of freedom of hard spheres, for which two sphere centers must not be closer than the sum of their radii. The exclusion angle ∆ serves as a tuning parameter for the equilibrium behavior of the system, which we have investigated in Ref. 11. The simplest possible phase diagram of hardcore spins is shown in Fig. 1 (c). At low exclusion angle, the system is weakly constrained and hence in a paramagnetic phase. As the exclusion angle is increased in equilibrium, the system undergoes an entropically driven Kosterlitz-Thouless (KT) transition at ∆ KT .
Crucially, in a constrained system topological defects can not only become confined but their existence can even become strictly incompatible with the hardcore con-straint. This happens at a point ∆ J > ∆ KT , and the defect density vanishes with a power law as ∆ → ∆ J [12].
In this letter, we generalize compression protocols originally developed for hard spheres [13] and study dynamics of hardcore spins far from equilibrium. Importantly, we uncover a remarkably simple mechanism precipitating of what might be called "arresting" or perhaps "jamming dynamics" in our model. In particular, any adiabatic compression protocol, namely one in which the system stays in equilibrium as ∆ is driven through ∆ J , reaches the ordered state with maximal exclusion angle ∆ max , analogous to the close packing. In contrast, a nonadiabatic protocol starting from a random initial state at ∆ = 0 will fall out of equilibrium and reach ∆ = ∆ J with nonzero defect density. As the existence of defects is incompatible with the hardcore constraint for ∆ > ∆ J , the protocol will fail beyond this point. The failure of compression at ∆ J coincides with a doubly algebraic divergence of the relaxation time of compressed states towards equilibrium-as a function of both system size as well as the distance to ∆ J . This can be fully understood in terms of a diffusion-annihilation process of topological defects, with a vanishing defect mobility.
In the remainder of this letter we first present hardcore XY spins on the square lattice as a minimal model exhibiting all the abovementioned features. We then show that on the triangular lattice, the relevant physics is realized in a more complex way: the additional chiral symmetry in the model leads to the presence of two kind of topological defects, that is, domain walls in addition to vortices. This leads to a more structured slow dynamics, but crucially, the behavior of the model can be accounted for by the same fundamental mechanism as before.
Model and phase diagram on the square lattice. We consider a system of XY spins, that is two-component unit vectors {S j }, on a lattice L subject to the constraint ϕ ij := arccos(S i · S j ) > ∆ ∀ ⟨ij⟩ ∈ L, where ⟨ij⟩ denotes an edge between site i and j. For XY spins on the square lattice, the hardcore constraint is illustrated in Eq. (1) (a), where for each spin, we denote in red the orientations forbidden while keeping all neighbors fixed. This model realizes the minimal equilibrium phase diagram sketched in Fig. 1 (a). At ∆ = 0, the system is unconstrained and naturally in a paramagnetic state, while at ∆ = π, the only allowed state (up to global symmetry) is the Néel state. Between these extremes, for any ∆ < ∆ max , long-ranged order is prohibited by an extension of the Mermin-Wagner theorem to constrained systems [14,15]. The system instead undergoes a KT transition at ∆ KT ≈ 0.435π [16], into a phase with quasi-long-ranged (QLR) order, with algebraically decaying correlations. Note that the model has no inherent notion of energy or temperature. Instead, Eq. (1) separates states into two classes -allowed and disallowed -without endowing them either with dynamics or with a notion of (high and low) energy. Thus, nontrivial correlations are of entropic origin, i.e. a form of "order by disorder" (OBD) familiar from the nematic ordering of Onsager's hard rods [17] and its descendants, including a recurrent appearance in frustrated magnetic systems.
It is well known that the KT transition can be understood as an unbinding transition of topological defects, which on the square lattice are vortices, see Fig. 1 (b). Defect unbinding is indeed the mechanism behind the equilibrium transition in our model, but within the QLRO phase, these vortices becomes incompatible with Eq. (1) for ∆ > π/2. This is shown most easily by realizing that on a bipartite graph such as the square lattice, any exclusion model as defined by Eq. (1) maps on an inclusion model defined bỹ where the spins on one of the two sublattices are flipped and ∆ incl = π − ∆. An antiferromagnetic vortex then maps onto a regular ferromagnetic vortex. Such a vortex however must have a core [11], which is a single plaquette with a winding of 2π. Since there are only four sites around a plaquette, this is only possible if ∆ incl > 2π/4 = π/2, which implies ∆ < π/2 in the exclusion model. Jamming hardcore spins. As a first step towards studying our model far from equilibrium, we define a notion of compression, which here will mean an increase of the exclusion angle ∆ while simultaneously evolving the state under some local dynamics. An analogous protocol for hard disks was implemented by Lubachesvky and Stillinger [13]. In this work, the radii of the disks were increased during a molecular dynamics simulation, while keeping the volume fixed. When the compression rate was slow with respect to the time scale set by the molecular dynamics, the final state was approximately close packed. In contrast, when the radii were increased very fast, the system ended up "jammed" in a polycrystalline state (to be contrasted with "maximally random jammed" states [18]), at a density well below close packing.
The Lubachevsky-Stillinger (LS) algorithm is straightforwardly adapted to and implemented for hardcore spins, the most subtle point being the choice of dynamics. Arguably the closest analog to the molecular dynamics studied by LS is given by continuous-time "Hamiltonian" dynamics. For hardcore spins, however, such dynamics are not ergodic even at infinitesimal exclusion angles since they preserve the local vorticity on each plaquette [12] (see also references [19,20] therein). Because of this, in the main text, we focus on a different kind, that is The compression protocol then proceeds as follows. Starting from a random state at ∆ = 0, one alternates N rattle Monte-Carlo sweeps (with one sweep consisting of N Monte-Carlo moves) and an increase of the exclusion angle by δ, where the increment δ is chosen as large as possible without invalidating the current state. Fast compression in this context then means large δ/N rattle . The protocol terminates if the increment δ repeatedly falls below a threshold value.
On a L = 40 square lattice, we run the LS compression protocol for hardcore spins with different values of N rattle to see whether jamming dynamics is observed. Out of 100 runs, with N rattle = 100 all runs terminate at ∆ max = π, that is they reach the Néel state which is our analog to close packing. In contrast, for N rattle = 1 all runs terminate at ∆ J = π/2. For N rattle = 10, runs fall into two classes, with 16 terminating at ∆ J and the rest at ∆ max . The "jammed" states at ∆ J have a finite ordered moment, which increases with N rattle but is always lower than the equilibrium expectation value [12]. To illustrate that the arresting dynamics here is indeed explained by our abovementioned mechanism, we show in Fig. 2 (a) a state at ∆ J for L = 10 and indicate the frozen vortex cores by dashed black boxes. These vortex cores are completely frozen under local dynamics since the constituent spins enclose an angle of exactly π/2. Relaxation dynamics. To corroborate the picture of defect freezing as the mechanism for the failure of nonadiabatic compression beyond ∆ J in our model, we study relaxation dynamics towards equilibrium close to this point in more detail. Generally speaking, one expects such a freezing to appear in conjunction with a diverging relaxation time scale.
While the LS protocol suffices to demonstrate the presence of jamming dynamics in hardcore spins, it is limited in that it allows only moderate compression speeds, resulting in partial equilibration, evidenced by a finite ordered moment of the final states. Because of this, we implement a second kind of compression protocol, originally developed by Xu et al. [22], which utilizes a softened constraint to compress faster and hence avoid partial equilibration. This is done by introducing an energy functional which is zero if Eq. (1) is fulfilled and introduces a quadratic energy penalty if neighboring hardcore spins overlap. The introduction of the soft constraint enables us to use much larger increments δ. This is because one does not have to choose δ such that the current state of the system stays valid, but instead one can choose it such that a conjugate gradient minimization step after the increment recovers a state with zero energy. The softcore compression protocol consistently yields far-from-equilibrium states with an ordered moment close to zero. We prepare such states at a range of exclusion angles ∆ close to ∆ J and show their defect density as a function of Monte-Carlo time in Fig. 2 (b). As expected from scaling arguments [23], it decays diffusively initially, but at long times, this behavior gives way to an exponential decay with a characteristic relaxation time τ defects . This relaxation time as a function of exclusion angle ∆ is shown in Fig. 2 (c), for a range of different system sizes. Evidently, the data is consistent with a doubly algebraic behavior with z = 2 and α = 3.74±0.50. These two power laws are conceptually quite distinct. The divergence of the relaxation times with system size, τ ∼ L z owes its existence to universal long-wavelength physics of phase-ordering kinetics under local dynamics and consequently its value  z = 2 [24] is quite robust. In contrast, the divergence of relaxation time as a function of exclusion angle is related to the defect mobility µ. It is measured by preparing two isolated vortices in an otherwise paramagnetic state. Their distance then follows the time dependence D(t) = √ D 0 − µt, from which µ can be determined by a linear fit. As shown in Fig. 2 (d), it vanishes as a power law as ∆ → ∆ J with, for a given local dynamics, roughly the same exponent α as the relaxation time. However, this exponent can be readily varied by changing the rules of the dynamical evolution even locally. For example, results from studying phase ordering kinetics under Hamiltonian dynamics plus tunneling are consistent with Eq. (5), with z = 2 but α = 1.87 ± 0.02 [12]. Triangular lattice. Finally, we compare the simple picture of the square lattice model to that on the triangular lattice. Its phase diagram, Fig. 3 (a), also has a paramagnetic phase including ∆ = 0, with a single (up to global symmetries) three-sublattice ordered state at ∆ = 2π/3; this comes in two chiralities that are not related by global rotation but instead by exchange of two sublattices. In between, chiral and quasi-long-range order develop either simultaneously or at two separate transitions [25,26] in proximity to ∆ ≈ 0.306π [12] (see also references [27][28][29][30][31] therein).
The properties of topological defects, and in particular their preclusion upon increasing ∆ remain central to the jamming phenomenology. These are richer than in the square lattice case. Vortices, as shown in Fig. 3 (b), become forbidden at a single packing fraction ∆ = π/3. In contrast, domain walls have internal structure (their shape) and become incompatible with the hardcore constraint over a range of exclusion angles. First, kinks in domain walls, shown on the left of Fig. 3 (c) become forbidden at ∆ = 2π/5 while at ∆ = π/2, any domain wall becomes incompatible with Eq. (1) [12].
In Fig. 3 (d), we show a histogram of the angle at which compression of random initial states using the LS protocol at different N rattle terminates. On the square lattice, such a histogram is strictly bimodal, with two narrow peaks at ∆ J = π/2 and ∆ max = π. In contrast, on the triangular lattice, for intermediate compression speed N rattle = 10, we see a continuous distribution between ∆ = 2π/5 with a pronounced peak at ∆ = π/2. This can be understood by considering the same mechanism as on the square lattice: fast compression (that is small N rattle ) fails because it arrives with a finite density of topological defects at a point where these become incompatible with Eq. (1). Now if there are multiple kinds of defects in the system, these can have different relaxation time scales and hence lead to a more complex dependence of the jamming angle on compression speed.
In this picture, it is somewhat surprising that in Fig. 3 (d), there is no peak at ∆ = π/3, where vortices become forbidden. This is just a consequence of the fact that vortices are only well defined in the presence of chiral order, whereas we start from random initial states, that are neither chirally nor QLR ordered and have no welldefined vortex density to begin with. A random initial state will however have a well defined and large domain wall density, which leads to the failure of fast compression (jamming) at ∆ J = 2π/5, which is where kinked domain wall become forbidden. At intermediate speeds, the system then is able to partially equilibrate and reach zero domain-wall-kink density, but still has smooth domain walls, leading to jamming at ∆ J = π/2. Between these two points, in the range 2π/5 < ∆ < π/2, there exist a multitude of local clusters with long but finite relaxation times, leading to jamming of protocols with intermediate compression speeds.
An important, qualitative difference to the square lattice model is that on the triangular lattice the only vanishing defect mobility in the model is that of domain walls at ∆ = π/2. Still, we observe a failure of the compression protocols at other points, where different defects become forbidden without a concomitant freezing. This is not unexpected since a long but finite relaxation time of such defects might still bring compression to a halt. In the language of granular materials, such freezing is fragile in that rattling of a jammed state might unjam it.
Conclusion. In summary, we have provided a detailed phenomenology and comprehensive understanding of arresting dynamics in hardcore spin models, uncovering an intricate interplay of lattice geometry, ordering and defects, and the dynamics at long and short wavelengths. Particularly noteworthy from a conceptual perspective is the role played by the (in)ability to anneal defects under a purely local dynamics. This fundamentally accounts for the phenomenon in this lattice model, with considerable added richness as a result of the variability and multiplicity of defects and their configurations. In addition to lending this model intrinsic interest on its own, the central role of defects to the (slow) dynamics resonates with a possibility previously formulated for the case of the glass transition [9]. The Lubachevsky-Stillinger protocol, inspired by Ref. 13 and described in the main text, was implemented in C ++ . It is given in pseudocode in algorithm 1. In Fig. S1, we summarize the statistical properties of the final states. In particular, we show histograms of the final exclusion angle in the top row and histograms of the structure factor peak value in the bottom row.
The fundamental limitation of this protocol is that what determines the largest possible increment of the inclusion angle ∆ is the smallest angle between any nearest neighbor pair in the system. This value however vanishes in the limit of infinite system size N → ∞, implying that for fixed N rattle also the effective compression speed ⟨δ/N rattle ⟩ → 0 as N → ∞. For L = 40 we find ⟨δ⟩ ≈ 10 −5 .
The softcore protocol, originally devised for hard spheres by Xu et al. [22] was implemented in C ++ as well. The main idea is to use the gradient of the potential to push overlapping spheres away from each other to produce a valid hard sphere state from a state with local violations of the constraint. The authors then use alternating compression and expansion to produce states where spheres are in close contact with each other.
We implement a compression step closely analogous to Ref. 22 for our hardcore spin model as follows, given as pseudocode in algorithm 2. Consider a valid initial state S = {S i } at some initial exclusion angle ∆, that is V (tot) (S, ∆) = ⟨ij⟩ V (ϕ ij , ∆) < ϵ for some fixed tolerance ϵ, which we always take to be ϵ = 2 × 10 −10 . Now, starting with compression, increase the inclusion angle ∆ by some trial value δ > 0, and use conjugate gradient descent to get from the state S, which has now a finite potential energy V (tot) (S, ∆), to a final state S ′ . Depending on the value of V (tot) (S ′ , ∆), we continue with either compression or expansion of the system. If V (tot) (S ′ , ∆) > ϵ, we revert the state of the simulation to its state before the energy minimization and reduce the exclusion angle ∆ by δ (expansion). If V (tot) (S ′ , ∆) < ϵ, we increase the exclusion angle ∆ by δ (compression). In both cases, we then minimize the potential energy again by conjugate gradient descent and repeat the procedure. Every time we switch from compression to expansion or Perform N vice versa, the increment δ is halved. The step terminates if ϵ/2 < V (tot) (S, ∆) < ϵ. Note that although we use the softened constraint in principle, forcing the system to have zero potential energy is equivalent to imposing the hardcore constraint. Allowing the potential energy to be small (ϵ ≪ N −1 ) but nonzero then ensures that the hardcore constraint is almost fulfilled by the state.
Similar to the LS protocol (algorithm 1), we augment the compression step with local Monte-Carlo dynamics as follows. If the difference between successive exclusion angles ∆ is below some threshold δ min , we reduce the exclusion angle ∆ by N exp · δ exp , where N exp is the number of expansion already performed since the last successful update of ∆, before applying N (exp) rattle Monte-Carlo sweeps and applying the compression step (algorithm 2) again. We repeat this procedure with increasing δ exp until either the exclusion angle is finally increased by more than the threshold or we terminate the procedure after a maximum number of expansions N (max) exp . The results of the LS protocol (algorithm 1) can be reproduced by setting the trial increment of the compression step δ 0 to a small value δ 0 = 10 −5 , and by introducing some extra rattling sweeps N rattle between successful increments. These extra rattling sweeps then set the timescale of compression, as N rattle does in the LS protocol. The protocol is given as pseudocode in algorithm 3. Here and in the main text, whenever we refer to states produced by the softcore compression protocol, we use N In Fig. S1, we show histograms of properties of the final states of the two protocols and in the case of the LS protocol also for different N rattle (indicated by color). In the top row, we show the final exclusion angle, and in the bottom row, we show the maximum value of the structure factor in reciprocal space. Since the LS protocol can only realize moderate compression speeds, the system partially equilibrates. This is reminiscent of challenge that arises for monodisperse disks, for which the protocol finds polycrystalline packings that do not pass muster as maximally random jammed states [13,18]. Note that the even though the system in the thermodynamic limit is not long range ordered, for any finite system size the maximum of the structure factor has a finite equilibrium expectation value. In particular, using the reflect algorithm developed in [11], we find max [S(q)] ≈ 0.450 (0.392) for the L = 40 (L = 80) square lattice at ∆ = π/2 and max [S(q)] ≈ 0.384 at ∆ = 2π/5 on the L = 42 triangular lattice. With sufficiently rapid compression, particularly with the softcore protocol, the order parameter stays well below this equilibrium value.

II. MEASURING THE DEFECT MOBILITY
In this appendix, we explain how the mobility of vortices, shown in ?? (d) of the main text as a function of exclusion angle ∆, is measured.
In the KT phase, a vortex-antivortex pair at distance D is subject to an entropic attractive interaction potential V (D) ∼ log(D), which makes the pair eventually annihilate. For a single pair, initially at distance D 0 , the distance as a function of time should then satisfy the differential equatioṅ where F (D) is the attractive force between the pair and µ is the mobility of vortices. The above equation has a simple solution We compute the mobility as a function of exclusion angle ∆ numerically as follows. For each value of ∆, we prepare single vortex-antivortex pairs at distance D 0 on a L = 60 lattice, by placing it into a Néel state and evolving the state using 10 6 Monte-Carlo sweeps, but fixing the vortex and antivortex in place. After this, we evolve the system further under the Monte-Carlo dynamics, but without fixing the vortex-antivortex pair. While doing so, we keep track of the distance of the pair as a function of time. We show the result of this procedure in Fig. S2. In (a), we show the distance between the vortex antivortex pair as a function of Monte-Carlo time, averaged over 10 4 initial states, for a range of initial distances at ∆ = 0.48π. As expected, for short times and long distances, it obeys the functional form of Eq. (S2). The mobility can then be extracted from a linear fit to D 2 (t) and is shown as a function of the exclusion angle in (b).
We have also measured the defect mobility for "Hamiltonian" dynamics, results of which are presented in Sec. IV C.

III. EQUILIBRIUM PROPERTIES
In this appendix, we provide additional details on the computation of equilibrium properties of hardcore spins which go beyond those already presented in Ref. 11. A. The vortex density on the square lattice vanishes as a power law within the KT phase In the main text, we described how the mechanism of jamming of hardcore spins is intimately related to topological defects. In particular, these defects become strictly forbidden at an exclusion angle ∆ J that is the defect density ρ must vanish as ∆ → ∆ J from below. In this section, we investigate how this happens.
Focusing on the square lattice model, vortices become forbidden at ∆ J > ∆ KT , that is inside the Kosterlitz-Thouless (KT) phase. In the range ∆ KT < ∆ < ∆ J , vortices are already confined, but fluctuations lead to a finite density where the sum is over all plaquettes p in the lattice, N p is the number of plaquettes (N p = N on the square lattice with periodic boundaries), and V p is the vorticity of the plaquette p where ∡ S i ,S j ∈ [−π, π] is the signed phase difference between the spinsS i andS j and the tilde denotes the fact that the spins are taken with respect to the Néel state. The vortex density in equilibrium on the square lattice is shown as a function of exlcusion angle ∆ in Fig. S3. The results are compatible with a power law, that is with ∆ J = π/2 and β ≈ 5.36 ± 0.30. The fact that the defect density ρ vanishes continuously (rather than via a first order transition) has nontrivial consequences also for the compression dynamics. As explained in the main text, slow compression is able to reach the exclusion angle ∆ max . This is because under slow compression, the system then stays in equilibrium and, because the defect density vanishes continuously as ∆ → ∆ J from below, has no defects when reaching ∆ J , where defects become strictly forbidden by the hardcore constraint. The protocol can then compress beyond that point, in contrast to a fast compression protocol which falls out of equilibrium and reaches ∆ J with a nonzero density of defects.

B. Forbidden defects on the triangular lattice
The phase diagram of hardcore spins on the triangular lattice is shown in Fig. S4 (a). As on the square lattice, there is a low-exclusion-angle paramagnetic and a large-exclusion-angle quasi-long-range ordered phase. However, the situation on the triangular lattice is richer because in addition to the O(2) symmetry of the XY spins, there is a chiral (Z 2 ) symmetry which the noncollinear ground state of the model breaks. Because of this, there are two kind of topological defects: vortices as well as domain walls. The former are defined with respect to a state of fixed chirality, whose noncollinearity leads to vortices and antivortices having very different free energy costs associated with their creation. The latter even become incompatible with the hardcore constraint at different exclusion angles depending on the their local shape [ Fig. S4 (b)].
To understand all of this in detail, we must get a sense of what kind of spin textures the topological defects correspond to. We begin by considering domain walls of chirality. For this, note that at ∆ max = 2π/3, not all allowed states can be transformed into each other via global rotations. Instead, there are two incarnations of the "120-degree" or " √ 3 × √ 3" state, which are related by the exchange of two sublattices. These two states are shown in the top right and bottom left of the right of Fig. S4 (b). There, we also show how they can be separated by a domain wall, indicated in pink. To distinguish the two domains quantitatively, we define a vector chirality κ t on each elementary triangle t of the lattice In the 120 degree state, the chirality is extremal on all triangles, that is |κ t | = 3 √ 3/2, and it has opposite sign on upwards and downwards triangles. This is also shown in Fig. S4 (b), where we indicate the sign of κ for each triangle as + or −. We call a domain wall as shown on the right of Fig. S4 (b) "smooth," because the domain wall takes no sharp turns or "kinks" as shown on the left of the same figure. Inspecting the smooth configuration, we recognize the 90-degree structure familiar We also indicate the succession of forbidden defects. The system exhibits extremely slow relaxation dynamics from 2π/5, with the dynamical transition thought to occur as ∆ → π/2. from the vortex cores on the square lattice. We thus conclude that these kind of domain walls become forbidden at ∆ = π/2. Turning to the domain wall kink, note that such a configuration necessitates three adjacent triangles to have the same sign of the vector chirality κ. This however already implies that the winding number along the loop encircling the three triangles is 2π, which becomes incompatible with the hardcore constraint for ∆ > 2π/5.
Turning to vortices, these are defined exactly as on the square lattice, but with respect to a 120-degree state of fixed chirality. In particular we define a vorticity V t where ∡ S i ,S j ∈ [−π, π] is the signed phase difference between the spinsS i andS j and the tilde denotes the fact that the spins are taken with respect to one of the 120-degree states. Note that this implies that vortices are only well defined in the presence of chiral order. For a uniform choice of chirality, we show a vortex (right) and antivortex (left) in Fig. S4 (c). As in (b), we highlight the bonds that minimize the enclosed angle between spins in red. In contrast to the square lattice case, the spin textures of vortices and antivortices on the triangular lattice are qualitatively different. Still, as we will show in the following, both become forbidden by the hardcore constraint at the same exclusion angle, that is ∆ = π/3. First, note that a perfect (meaning unperturbed) antivortex core is incompatible with the hardcore constraint for any nonzero exclusion angle. What is still possible is a state in which the configuration on the center triangle of the antivortex is perturbed. For this, two of its spins are rotated in opposite directions, away from the third. For the center triangle on the right of Fig. S4 (c), one possibility of this would be to rotate the bottom right spin of the triangle clockwise and the bottom left spin anticlockwise. Via such a rotation by an angle ϵ, we can produce a valid configuration on the center triangle for any ϵ < ∆. However the vorticity on the central triangle for such state is only equal to −1 as long as |ϵ| < π/3, which means that for |ϵ| > π/3 it is not an antivortex core. Since an antivortex needs a core, its existence is hence only compatible with the hardcore constraint for ∆ < π/3. Second, the vortex core as constructed on the right side of Fig. S4 (c) is just the 120-degree configuration of opposite chirality, which might seem to suggest compatibility with the hardcore constraint for all exclusion angles up to ∆ max . However, a vortex is not just a core but an extended object. The next largest closed loop to consider is of length six and indicated in red on the left of Fig. S4 (c). The spin configuration of this loop in a perfect vortex is clearly only compatible with the hardcore constraint for ∆ < 2π/6 = π/3. At this point, only three spins with any available phase space are those at the three corners, which however can only be deformed in a way that already flips the vorticity of the whole loop.
The exclusion angles at which the different topological defects become forbidden by the hardcore constraint are summarized in Fig. S4 (a).

C. Details on the numerical estimation of ∆KT on the triangular lattice
Here, we give the details on our numerical computation of the equilibrium phase diagram of hardcore spins on the triangular lattice. We use a variant of the Wolff algorithm, in which a single Monte-Carlo step consists of a global reflection of a cluster of spins. The axis of reflection is chosen at random and a "pocket" is initialized to contain a random seed spin. Each spin in the pocket is reflected about the chosen axis, and if doing so leads to a violation of the hard constraint with its neighbors, these neighbors are added to the pocket. The move terminates when the pocket is empty, i.e., when the lattice has returned to an allowed spin configuration. By definition, the "cluster" consists of all spins that are added to the pocket at some point in the move. We call this the reflect algorithm.
While the reflect algorithm is provably ergodic on bipartite lattices [11], we consider it likely that the triangular lattice hosts exceptions. One indication of this is that the average cluster size near the critical point scales as L 2 , whereas the square lattice enjoys a superior scaling of L 2−η ′ , with exponent η ′ in close agreement with the critical exponent η of the staggered spin susceptibility [11]. Nevertheless, several other metrics-including the decay of the autocorrelation function of the spin susceptibility, the ability of cluster moves to unwind vortices and extended defects, and the convergence of different initial states to the same time-averaged expectation values-support the use of the algorithm as an approximate probe of the equilibrium phase diagram, especially at low ∆.
Based on studies of the thermal model [25], the spin and chiral transitions are expected to at very close exclusion angles, at or below π/3, where vortices become forbidden. Finite-size scaling on lattices up to L = 54 at angles ranging from 0.3π to 0.4π support this hypothesis, suggesting a phase transition around 0.307π.
Looking first at the spin transition, since the associated U (1) symmetry is continuous, only QLRO is permitted for ∆ < ∆ max , with a continuously varying exponent η in the critical phase. This exponent is determined from the scaling of the spin susceptibility, or structure factor, defined as the maximum in the Fourier spectrum of the spin correlation function: Triangular lattice exclusion models in the QLRO phase exhibit a peak at q s = (4π/3, 0) corresponding to the "120-degree" or " √ 3 × √ 3" order. The susceptibility χ s = S(4π/3, 0) then scales as a power law with the finite system size χ s ∼ L 2−η , as shown in Fig. S5 (a).
If the spin transition belongs to the Kosterlitz-Thouless universality class, the critical exponent will assume the value η = 1/4 at the critical angle ∆ KT . Precision numerical studies are split on whether the magnetic transition in this model and the closely related fully frustrated XY model is in the KT universality class [27? , 28] or of a non-KT character [25,26]. At the system sizes considered here we cannot draw a firm conclusion regarding the universality class for the constraint-only model, but the data in Fig. S6 (a) are broadly consistent with a KT transition. First, η s = 1/4 (the value at the KT transition) around ∆ ≈ 0.307π, which is also the angle below which the power-law scaling of χ s becomes a poor fit, indicating a transition to the paramagnetic phase. In this neighborhood of ∆, we also examine the Binder ratio (middle panel), which should also take a universal value at the critical point, and become asymptotically independent of L above this angle. A variety of normalizations exist in the literature, but we use the definition [29]: where M s is the K-point magnetization [25]. Finally, we measure the second-moment correlation length (bottom panel), defined as [25]: ξ 2nd (L) = 1 2 sin(π/L) S(q s )/S(q s + 2π Lx ) − 1. (S10) To leading order, ξ 2nd /L becomes independent of system size in the QLRO phase and takes a universal value in the thermodynamic limit at the KT transition [30]. Turning to the chiral phase transition, we measure analogous quantities to Eqs. (S8), (S9) and (S10), with the chirality (Eq. (S6)) serving as the order parameter in place of the K-point magnetization. In this case, since the associated symmetry is discrete, the transition is to a long-range-ordered phase. Consequently, in the thermodynamic limit, the connected correlation function is algebraically decaying only precisely at the critical point. In terms of the finite-size scaling of the chiral susceptibility, this means that the power-law scaling with L is only exhibited in the close vicinity of the transition, as shown in Fig. S5 (b). The system-size independence of the Binder ratio and second-moment correlation length, as an additional signature of critical behavior, also appears only in the vicinity of the transition. These signatures are shown in Fig. S6 (b), again indicating a critical point in the range of ∆ ≈ (0.306 ± 0.001)π. Consistent with the expectation on grounds of symmetry that the chiral transition belongs to the Ising universality class, the value of η deduced from a power-law fit to χ is close to η Ising = 1/4. Moreover, the value of U c attained at the intersection of different system sizes is consistent with the Binder ratio at the Ising critical point on a triangular lattice geometry [31], a universal value that has also been observed at the chiral transition of the antiferromagnetic XY model [? ].
To summarize, within uncertainties, the data are suggestive of coincident chiral and spin transitions, with the latter falling in the KT universality class. The possibility of two separate transitions cannot be ruled out; in the thermal model, the difference between the two transitions only manifests above a threshold system size of L ≈ 300, so the same effect may be present here [25]. Regardless, the existence of a transition below ∆ = π/3 indicates that, as on the square lattice, there exists a range of exclusion angles where defects are geometrically allowed but entropically disfavored, leading to quasi-long-ranged order by disorder. In terms of the dynamics, this means that the jamming point and its concomitant doubly algebraic relaxation dynamics occur within the QLRO phase.

IV. HAMILTONIAN DYNAMICS
Here we outline the main results of the study of Hamiltonian dynamics on the square lattice. In these dynamics (algorithm 4), each spin is initialized with an angular velocity drawn from a Gaussian distribution, and we work in the frame in which the net angular momentum is zero. The system is evolved until a pair of neighboring spins "collides," i.e. saturates the hard constraint, at which point the two colliding spins swap velocities.

A. Conservation Laws
The Hamiltonian dynamics are more constrained than the stochastic Monte Carlo dynamics in a few crucial ways.
First, the vorticity on each plaquette with respect to the antiferromagnetic state, is conserved. To see that this is true, denote the four spins belonging to a plaquette as ϕ 0 , ϕ 1 , ϕ 2 , ϕ 3 ∈ [0, 2π), where without loss of generality we can define the x axis such that ϕ 0 = 0 at a given instant in time. Denoting spin orientations with respect to the Neél state byS j = (−1) j S j , the enclosed angle betweenS andS ′ is defined in the interval [−π, π]: (S12) Plugging this into the expression for the vorticity yields: where H(x) is the Heaviside function. Since the relative ordering of ϕ 1 , ϕ 2 and ϕ 2 , ϕ 3 is preserved by the dynamics, V p (t) depends only on the initial conditions. Thus, there is one local conserved quantity per plaquette. A square lattice of N spins has N plaquettes, so this is enough to make the system integrable [? ]. In addition, periodic boundary conditions impose that the net topological charge on the lattice is zero.
An additional global constraint comes from the fact that the colliding spins behave like equal-mass objects in 1D, which simply exchange velocities upon collision. This implies that the set of velocities is fixed by initial conditions. In the problem of hard rods confined to a 1D box, this conservation law enables an explicit solution [19,20]. First, positions are redefined such that the rods have zero length. Then, the positions of these point particles at later times are given by free trajectories labeled {x i (t)}, but with each intersection of trajectories inducing a swap of the labels. However, the difference between translational degrees of freedom in real space and orientational degrees of freedom in spin space complicates attempts to apply a similar solution strategy to the problem of hardcore spins. Rather than pursuing this strategy here, in keeping with the main focus of this letter, we will lift some of the constraints to recover a form of "jamming."

B. Tunneling
Among the several quantities that exhibit power-law scaling near jamming in the Monte Carlo dynamics is the vortex mobility, which vanishes with power α as ∆ → ∆ J [??(d)]. Since algorithm 4 imposes zero mobility at all angles, this constraint must be lifted to observe a similar scaling in the Hamiltonian dynamics. A minimal modification that achieves this end is a change to how the spins are updated in a collision. In place of line 12, in which the colliding spins exchange velocities, the update in algorithm 5 is applied: one of the two spins is chosen at random to "tunnel" through the other, i.e., is reflected about the axis of the other spin (lines 1-2). If the proposed reflection causes the tunneled spin to violate the exclusion constraint with one of its other neighbors, the move is rejected and velocities are exchanged instead (lines 5-9), as in the original algorithm.

Algorithm 5: Hamiltonian Dynamics Tunneling
Step The reflection flips the sign of the phase difference between S j and S k , so by definition this tunneling move creates a vortex-antivortex pair on the two plaquettes that have (i, j) as common vertices. For exclusion angles above ∆ > π/2, defects are forbidden, so the tunneling moves are rejected with probability one. Likewise, in stochastic Monte Carlo dynamics, if a proposed spin flip S i → S ′ i induces a change in the sign of the phase difference with a neighboring spin, then it is guaranteed to violate the hard constraint for ∆ > π/2. In both cases, the dynamics just move through the sector with zero vorticity on each plaquette.
For ∆ < π/2, on the other hand, tunneling moves are accepted with nonzero probability. As the exclusion angle approaches the ∆ J = π/2 from below, the acceptance probability vanishes, which we expect to be accompanied by a vanishing of the vortex mobility.

C. Vortex mobility and diverging timescale
To measure the vortex mobility, as with the Monte Carlo dynamics discussed in the main text, we can measure the diffusion-annihilation process by tracking the separation of a vortex-antivortex pair as a function of time. The initial states used are the same as those in the main text, prepared by placing a single defect pair into a Neél state at a distance D 0 and evolving the state with 10 6 Monte Carlo sweeps. These states, whose structure factor is suppressed below the equilibrium value due to the presence of the defects, are then fed into the Hamiltonian dynamics + tunneling algorithm, where the pair is no longer fixed. The results are shown in Fig. S7. As anticipated from the theoretical considerations of phase-ordering kinetics, the vortex separation decays as D 2 ∝ D 2 0 − µt, where the mobility µ depends only on the exclusion angle ∆. As ∆ approaches π/2, µ scales as a power law, but with a different exponent from the Monte-Carlo dynamics.
The vanishing mobility of defects coincides with a diverging timescale for relaxation towards equilibrium. In the main text, this is probed by tracking the defect density as a function of time starting from far-fromequilibrium initial states [??(b-c)]. For the Hamiltonian dynamics, we found the structure factor to be a more reliable probe of the diverging timescale. The structure factor approaches its equilibrium value exponentially, with a dynamical exponent z = 2, just as for Monte Carlo dynamics. For a given system size, the relaxation time diverges as a power law ∼ (∆ J − ∆) −α with a similar exponent α to that of the vortex mobility. More robust results for α are obtained by returning to the initial states containing a single vortex-antivortex pair. After an initial transient, S(π, π) eq − S(π, π)(t) exponentially decays with a relaxation time that increases with the initial separation [ Fig. S8(a)]. As a function of ∆ J −∆ [ Fig. S8(b)], τ diverges with the same exponent, within uncertainties, as the vanishing mobility; the latter exponent is estimated to be α = 1.87 ± 0.02.  S8. (a) Exponential approach of the structure factor S(π, π) to its equilibrium value, for L = 60 using the same initial states as in Fig. S7. Hamiltonian evolution with tunneling is performed at exclusion angle ∆ = 0.48π, fit to an exponential C exp(−t/τ ) (black dashed lines). (b) The relaxation time τ behaves as a power law in π/2−∆, with the same exponent as the mobility. Consistent results are obtained for different initial separations; the reported value of α is an average over the best fit values for each D0.

V. SOFTCORE COMPRESSION WITH FORBIDDEN DOMAIN WALL KINKS
In this appendix, we give details on and present the results of softcore compression on the triangular lattice in the absence of domain wall kinks. A domain wall kink is a configuration as on the left of Fig. S4 (b), that is two domain walls meet at an angle of π/3 or, equivalently, three adjacent triangles have the same vector chirality. We begin by describing the adjusted preparation of initial states and subsequently discuss the results of the softcore compression protocol (algorithm 3) on the kink free model. In contrast to softcore compression in the presence of kinks, we here find (Fig. S10) a broad distribution between 2π/5 and π/2 of the exclusion angle ∆ J at which the protocol terminates for different runs, superimposed with pronounced peaks at rational fractions of 2π. The weight of the distribution shifts towards larger ∆ as N (exp) rattle is increased, without ever exceeding ∆ = π/2. We interpret this as evidence for the fact that we expect the first infinite relaxation timescale of the system to occur at ∆ = π/2, where domain walls become strictly incompatible with the hardcore constraint. We attribute the presence of a broad distribution of ∆ J for small N (exp) rattle to the physics of small finite clusters and a general slowing down of dynamics due to the large ∆ in conjunction with the large connectivity of the triangular lattice.

A. Generation of kink free initial states
We generate kink-free initial states by starting with a random state and subsequently removing kinks by a combination of a special "kink-removing" Monte-Carlo update and regular single-spin-randomization dynamics. The "kink-removing" Monte-Carlo update is shown in Fig. S9. It consists of a swap of the two spins on the shorter of the two parallel edges of the parallelogram defined by the three adjacent triangles of equal chirality. Performing multiple sweeps of such updates on a random state at ∆ = 0 removes most, but not all domain wall kinks. We therefore alternate such sweeps with regular single-spin-randomization Monte-Carlo sweeps, but with the additional constraint that no single update should increase the total number of kinks in the system. This procedure after a few iterations successfully removes all domain wall kinks from the state, without inducing any obvious form of chiral or orientational order.

B. Distribution of ∆J
Using the kink-free initial states at ∆ = 0, we now run the softcore compression protocol (algorithm 3) with the additional constraint on the dynamics that no update is allowed to introduce a domain wall kink into the system.
In Fig. S10, we show a histogram of the exclusion angles ∆ J at which the protocol terminates for 10 3 initial states. First, in contrast to the model with kinks (see Fig. S1) here we find a broad distribution of ∆ J . We attribute this somewhat surprising result to the fact that for ∆ > 2π/5, the system finds a multitude of ways to exhibit slow dynamics with contributions from both, the nonuniversal physics of small clusters as well as the general slowing down at these exclusion angles. Note that both of these phenomena are more important on the triangular lattice than on the square lattice because of the larger connectivity of the former. The pronounced peaks on top of the broad distributions are located at rational fractions of 2π and are a consequence of what we call Ouroboros configurations. These are configurations on a closed loop C of length L on the lattice, such that ϕ j = 2πnj L , j ∈ C (S14) for some integer n. They are local minima of the energy functional in Eq. (4) of the main text and become forbidden at exclusion angles ∆ = 2πn/L. At these values of ∆, they introduce long but finite relaxation timescales into the system (see also discussion in the main text). These attributions are evidenced by the fact that as N (exp) rattle is increased, that is compression is performed more slowly, the broad, continuous distribution gives more and more away to a collection of peaks. Second, the weight of the distribution shifts to larger ∆ J as N (exp) rattle is increased and for N (exp) rattle = 10 5 , 189 of the 1000 runs even reach the 120-degree state at ∆ = 2π/3. However, no run ever fails between π/2 and 2π/3. This implies that smooth domain walls of chirality are indeed the last defects to become forbidden by the hardcore constraint.