Unidirectional subsystem symmetry in a hole-doped honeycomb-lattice Ising magnet

We study a model of a hole-doped collinear Ising antiferromagnet on the honeycomb lattice as a route toward the realization of subsystem symmetry. We find nearly exact conservation of dipole symmetry verified both numerically with exact diagonalization (ED) on finite clusters and analytically with perturbation theory. The emergent symmetry forbids the motion of single holes -- or fractons -- but allows hole pairs -- or dipoles -- to move freely along a one-dimensional line, the antiferromagnetic direction, of the system; in the transverse direction, both fractons and dipoles are completely localized. This presents a realization of a `unidirectional' subsystem symmetry. By studying interactions between dipoles, we argue that the subsystem symmetry is likely to continue to persist up to finite (but probably small) hole concentrations.

Introduction.Fractons are the newest addition in the line of exotic quasiparticles in condensed matter physics.Recent years have witnessed tremendous progress in understanding fractons , which unmasked connections with other areas of physics including topological order [2,4,8,9], gauge theory [12,13], quantum computing [4], glasses and soft matter [21,22], These exotic properties result from unusual mobility constraints whereby a fracton is a charge-like excitation which has restricted mobility when in isolation, but which, nonetheless, can easily move in a subdimension of space when bound to an oppositely charged fracton in a dipolar bound state [41].
The realization of these aberrant mobility constraints and subsystem symmetry in physical systems is, however, not naturally available, stimulating many interesting theoretical proposals [42][43][44].One particularly appealing approach is based on the idea that defect motion-induced frustration in an ordered background leads to emergent immobility constraints on the motion of defects themselves which in turn gives rise to fractonic quasiparticles, with hole-doped Ising antiferromagnets [45][46][47][48][49][50][51][52][53][54] being a prime example of this mechanism [27,28].Unlike other constructions proposed to realize fracton topological order (e.g.[55]), this approach does not require an extensive number of locally conserved quantities (which results in an extensively degenerate ground state), making the physical realization of fracton conservation laws (without topological order) in quantum magnets significantly more accessible.However, the emergence of fracton-like quasiparticles in fully two-dimensional (2D) hole-doped Ising antiferromagnets has been demonstrated only in the asymptotic limit of t J (where t is the hole hopping and J is the Ising coupling) [27], where already at the leading order in t 2 /J, the direction of the dipole moment is not conserved even though its magnitude is because hole pairs can rotate as they move in the 2D plane.
In this letter we overcome this limitation and pro- pose an essentially exact physical realization of fractonic quasiparticles with subsystem dipolar symmetry in a 2D non-degenerate, ordered spin system with local two-spin interactions.The key ingredient in our proposal is the collinear antiferromagnetic order in which defect motioninduced frustration completely prevents hole pairs or dipoles from moving in the perpendicular direction, resulting in exact conservation of both the magnitude and direction of dipole moment in the asymptotic limit t J.This further endows the system with a subsystem symmetry since dipole motion is restricted to a one-dimensional (1D) submanifold of the system.This symmetry manifests only along the antiferromagnetic direction, thus we denote it as 'unidirectional' subsystem symmetry.Importantly, we find via exact diagonalization (ED) that this symmetry continues to hold quantitatively away from the perturbative limit when t J.We further suggest that interaction between hole pairs is weak implying continuity of these immobility constraints to small, finite concentrations.
Model.Using a coupled spin chain construction to ferromagnetically couple antiferromagnetic Ising chains we construct a 2D ordered Ising magnet, in which we study doped holes.Specifically, we consider a model of holes doped into an Ising collinear antiferromagnet on the honeycomb lattice given by where t denotes the hopping amplitude, c † (c) are fermionic creation (annihilation) operators, and σ ∈ {↑, ↓ } is the fermion spin.The coordinates r i are defined by the sites of a Bravais lattice given by r j = m j a 1 + n j a 2 where a 1 , a 2 denote the primitive lattice vectors of the honeycomb lattice.The basis has two sites referred as A and B. The vectors δ x , δ y , δ z connect spins on nearestneighbor sites in the three directions of the honeycomb lattice, see Fig. 1.A no double occupancy constraint, is enforced on the Hilbert space.The Ising spin Hamiltonian H Ising is constructed by ferromagnetically coupling alternating sites on antiferromagnetic chains (this can be viewed as a model of a striped antiferromagnet in a brick-wall lattice).The ferromagnetic spin couplings are taken to be along the δ z direction and antiferromagnetic spin couplings along the δ x and δ y directions: where S z ri denotes the spin-z operator.One ground state of H Ising is given by |ψ . A hole can be created (annihilated) on a given site by annihilating (creating) an electron on that site.For any given hole density, the total magnetization is conserved.Thus, we associate the removal (addition) of a fermion with spin σ with the creation (annihilation) of a hole with spin −σ, as either amounts to a total net change of the magnetization of the entire system by −σ.Therefore the hole creation operator is given by f † ri,σ = c ri,−σ .A hole can move to a neighboring site along the antiferromagnetic direction if the electron with antialigned spin on that site moves to the hole's original site.One can view this as a spin flip operation at the original hole site accompanied by hopping of the hole to the concerned neighbor site.The original hole site with a flipped spin is now in a "wrong" orientation with respect to its two remaining neighbors and thus we view this as defect creation.A hole dressed by such bosonic (spin wave) defects forms a magnetic polaron.We can represent a misaligned spin as a bosonic magnon defect for the sites r i .(aA unique characteristic of this model on the honeycomb lattice is that each site has two antiferromagnetic neighbors in the x and y directions and one ferromagnetic neighbor in the z direction.This means that in order for a hole to move from a site to another in the x or y direction it must create a defect on the site of its departure or annihilate a defect on the site of its arrival.In contrast, hole hopping in the z direction will not involve any defect creation.Thus, in the language of the hole and magnon defect operators our model in Eq. ( 1) can be recast as In Eq. ( 3) a no-double occupancy constraint at each lattice position is imposed, as there can be either a defect or a spin in a given lattice site.We ascribe an effective charge degree of freedom to the dressed hole in terms of its spin density: , where τ z is the Pauli-z matrix in the hole spin flavor space.The total charge Single hole.We first consider the single-hole sector We find an effective Hamiltonian for the single hole in the limit t/J 1 by treating the hole-boson coupling term in Eq. ( 3) perturbatively.The effective Hamiltonian, up to second order, is given by The first term in this expression reflects the localization of single holes where 2t 2 /J is the formation energy associated with creating an immobile polaron.This can be understood as follows.Since hole motion in the x and y directions frustrates the antiferromagnetic bonds, the hole must retrace its path and return to its original site in order to heal the background [27], thus a single hole cannot move along the antiferromagnetic direction.This is represented pictorially in Fig. 2 where we demonstrate that a single hole can hop only between two sites in the z direction via the second term in Eq. ( 4), but cannot hop in the x and y directions without creating misaligned spins even following a move in x direction, and thus the hole is localized on the z bond connecting the two sites.A single hole can only move away from its original site via virtual processes involving closed loops known as Trugman loops [45].In contrast to the case of a square-lattice antiferromagnet where Trugman loops first appear at sixth order in perturbation theory, Trugman loops in our model first appear only at fifteenth order.Thus, we expect fracton physics to persist in a parama-trically larger regime in t/J.Furthermore, away from the The insets depict the momenta resolved by the simulation cluster.The symbols of the energy levels denote different point group representations.We observe a flat spectrum in both momentum directions for n h = 1 and a flat spectrum in the ky direction only for n h = 2.This implies full localization of the single hole and localization of the two holes along the δ z direction only.
perturbative limit, we expect these closed loops to be energetically much more costly and thus their contribution to hole motion to be negligible [54].To confirm this behavior beyond the limits of applicability of perturbation theory, we use ED [56] for the model with J/t = 0.4 on N = 24 and N = 32 site clusters.Figure 3(a), (c) show the energy spectrum as a function of momentum in the Brillouin zone for n h = 1.We observe nearly exact degeneracy of the spectrum at all momenta.This fully flat degeneracy in momentum space indicates localization in real space, as expected from our arguments for the spin polaron [57].
Two holes.Next, we consider the two-hole sector n h = ri,σ f † ri,σ f ri,σ = 2.We derive an effective Hamiltonian for two holes perturbatively, which, up to second order in t/J is given by: Here, the first and second terms correspond to single-hole processes derived in Eq. ( 4), while the third and fourth terms correspond to near-neighbor repulsive densitydensity interactions between magnetic polarons, and the fifth and sixth terms correspond to pair hopping [58,59] along the x − y direction.This effective Hamiltonian shows that two neighboring holes oriented along the x − y line can move in a bound state via a second-order process in which one hole hops by creating a bosonic defect which is then absorbed by the second hole only if it follows its partner along the x or y directions.In contrast, hole pairs cannot move along the z -direction because any such motion frustrates the antiferromagnetic bonds, in which case the holes have no option but to retrace their paths to their original non-frustrating configuration.One can visualize these processes pictorially in Fig. 2 which shows that two holes can move only together and only in the antiferromagnetic x -y direction, but not through the ferromagnetic z direction [60].This picture is confirmed via ED for the model with J/t = 0.4 on N = 24 and N = 32 site clusters.Figure shows the energy spectrum as a function of momentum in the Brillouin zone for n h = 2.We observe that the lowest-energy states in the different momentum sectors are nearly exactly degenerate only along the k y -direction of the Brillouin zone, indicating localization of states solely along the δ z direction in real space [61].
The phenomenology of one-and two-hole states in our model implies an emergent unidirectional subsystem symmetry along the antiferromagnetic direction with conservation of both charge and dipole moment (defined as D = i q ri r i ).Here a single hole forms a spin polaron which is almost perfectly localized up to very high order in perturbation theory mimicking a fracton, while two near-neighbor holes form a spin bipolaron which moves with ease along an antiferromagnetic 1D submanifold of the system while conserving dipole moment exactly like a dipole.The effective Hamiltonian in Eq. ( 5) manifestly conserves the total charge.The conservation of dipole moment, [D, h 2nd 2h ] = 0 [62], can be seen from Eq. ( 5) which shows that hole pairs can only move together whilst conserving their relative separation and cannot rotate [63].We note that unlike the case of the square-lattice antiferromagnet [27], there is no need to impose external energetic constraints to realize fracton physics in our model, which appears to hold robustly even away from the perturbative limit t/J 1 as seen in ED.Furthermore, our results will continue to hold even if the magnitude of J along the antiferromagnetic direction is different from that along the ferromagnetic direction, making physical realization more readily accessible.
Finite density of holes.Having established unidirectional dipole symmetry, we address the question of dipole-dipole interactions at finite hole concentrations.We study numerically via ED n h = 4 holes in a N = 24site cluster (large clusters are beyond reach).Figure 4 shows a nearly degenerate spectrum along the k y momentum direction corresponding to the ferromagnetic direction δ z in real space suggesting that dipole conservation may persist at finite hole densities.To investigate interactions between dipolar pairs we compute the pair-pair binding energy E pair defined as For the N = 24 site cluster, at J/t = 0.4 we find E pair /t = 0.13, suggesting that dipole-dipole interaction is repulsive.However, it is not clear to what extent these results are sensitive to finite-size effects.We argue that current indications suggest that dipole symmetry may play an important role at finite hole concentrations, and hope to address this in future work.Conclusion.We considered the physics of fractons and dipoles emergent in a hole-doped collinear antiferromagnets on a honeycomb lattice.By means of analytical arguments and ED, we showed that individual holes are completely localized in the 2D system, while near-neighbor hole pairs form dipolar lineons which can move freely only along the antiferromagnetic direction.These observations reflect an emergent quasi-exact unidirectional subsystem symmetry along the antiferromagnetic direction.These results were obtained for an Ising magnet, but, based on perturbative arguments [27], we expect a sufficiently small spin exchange J ⊥ to not affect our results significantly [64,65].Our results indicate that dipole symmetry is a robust feature in the limit of a single and a pair of doped holes and may persist to finite hole concentrations where dipole-dipole interactions are nontrivial and have implications that will be the subject of future work.Another promising future direction involves using a coupled plane construction analogous to our approach to engineer a three-dimensional model with exact subsystem symmetry along lines or planes.We hope that such approaches will enable the study of the exotic properties of fractons such as their unusual dynamical behavior in simple, potentially accessible spin systems.
We acknowledge useful discussions with M. Bukov, A. L. Chernyshev, J. Romhanyi, C. Xu, and especially with M. Berciu, V. Calvera, K.-S.Kim, A. Prem, K. Wohlfeld and H. Yan. S. S. acknowledges support from Science and Engineering Research Board (Department of Science and Technology) Govt. of India, under grant no.SRG/2020/001525 and an internal start up grant from Indian Institute of Science Education and Research, Tirupati.J. S. acknowledges support from the Gordon and Betty Moore Foundation's EPiQS Initiative through Grant GBMF8686 at Stanford University and the National Science Foundation (NSF) Materials Research Science and Engineering Centers (MRSEC) program through Columbia University (where this work was initiated) in the Center for Precision Assembly of Superstratic and Superatomic Solids under Grant No. DMR-1420634.J. S. also acknowledges the hospitality of the Center for Computational Quantum Physics (CCQ) at the Flatiron Institute.The Flatiron Institute is a division of the Simons Foundation.Exact diagonalization calculations were performed using the Hydra software package [66].

FIG. 1 .
FIG. 1. Ground state of the Hamiltonian HIsing (Eq.(2)) with one electron per site.The three bonds of the honeycomb lattice are denoted as x , y and z respectively (x and y are the Cartesian directions).The Hamiltonian describes an Ising magnet on a honeycomb lattice with antiferromagnetic exchange along the x and y bonds and ferromagnetic exchange along the z bonds.The vectors δ x , δ y and δ z connect near-neighbor sites in the x , y and z directions, respectively.The unit cell of the honeycomb lattice with A and B sublattices is shown in the dotted region, and a1 and a2 are the primitive lattice vectors.

FIG. 2 .
FIG. 2. (a)An isolated hole can move by one site only in the ferromagnetic direction but cannot move in the antiferromagnetic direction without frustrating the antiferromagnetic bonds.(b) A pair of holes on neighboring sites can move only along the antiferromagnetic x − y direction.Note that a pair of holes on neighboring sites connected by a bond in the z direction cannot move without frustrating the background and as such will be localized.
and b † ri+δ z = σ − ri+δ z , and for the sites r i .(a 1 + a 2 ) ∈ 2Z + 1 as b ri = σ + ri , b ri+δ x = σ − ri+δ x , b ri+δ y = σ − ri+δ y , and b ri+δ z = σ + ri+δ z , where σ ± ri are the Pauli Ladder operators.In contrast to motion along the antiferromagnetic direction, a single hole can move -only by only one sitein the ferromagnetic z direction since there no wrongly aligned spin.

FIG. 3 .
FIG. 3. Energy spectrum as a function of momentum from ED on an N = 32 ((a) and (b)) and N = 24 ((c) and (d)) cluster for t = 1.0 and J = 0.4.Results for a single hole are shown in (a) and (c), and results for two holes are shown in (b) and (d).The insets depict the momenta resolved by the simulation cluster.The symbols of the energy levels denote different point group representations.We observe a flat spectrum in both momentum directions for n h = 1 and a flat spectrum in the ky direction only for n h = 2.This implies full localization of the single hole and localization of the two holes along the δ z direction only.

4 .
FIG. 4. Energy spectrum as a function of momentum for n h = 4 holes from ED on an N = 24 cluster for t = 1.0 and J = 0.4.The momenta Γ and M1 are degenerate up to a difference of ∆/t = 1.4 • 10 −3 , indicating localization along the ky direction.