Emergent Snake Magnetic Domains in Canted Kagome Ice

We study the two-dimensional kagome-ice model derived from a pyrochlore lattice with second- and third-neighbor interactions. The canted moments align along the local $\langle 111 \rangle$ axes of the pyrochlore and respond to both in-plane and out-of-plane external fields. We find that the combination of further-neighbor interactions together with the external fields introduces a rich phase diagram with different spin textures. Close to the phase boundaries, metastable $\textit{"snake"}$ domains emerge with extremely long relaxation time. Our kinetic Monte Carlo analysis of the magnetic-field quench process from saturated state shows unusually slow dynamics. Despite that the interior spins are almost frozen in snake domains, the spins on the edge are free to fluctuate locally, leading to frequent creation and annihilation of monopole-anti-monopole bound states. Once the domains are formed, these excitations are localized and can hardly propagate due to the energy barrier of snakes. The emergence of such snake domains may shed light on the experimental observation of dipolar spin ice under tilted fields, and provide a new strategy to manipulate both spin and charge textures in artificial spin ice.


I. INTRODUCTION
Metastability in physical systems indicates a rough free energy landscape with multiple nearly degenerate local minima. Normally, it requires random disorders for the metastable behavior to occur, with spin glass as the most prominent example [1]. On the other hand, disorder-free geometrically frustrated spin-ice compounds, such as Ho 2 Ti 2 O 7 and Dy 2 Ti 2 O 7 also suffer from spin freezing at low temperature due to the slow dynamics originated from the extensive quasi-degenerate low-energy states, and the topological constraint known as the ice rule [2,3]. For pyrochlore lattices with corner-sharing tetrahedra, the ice rule dictates that the spins in each tetrahedron should obey the two-in-two-out constraint, resulting in a divergence-free flux field on the dual diamond lattice. The ice-rule constraint originated from the effective nearestneighbor coupling of spins can be regarded as a conservation law of emergent gauge fields, giving rise to the anisotropic dipolar correlation [4]. The long-range dipolar interaction in the more sophisticated dipolar spin ice (DSI) model merely provides correction to the projectively equivalent band of lowenergy states [5]; the quasi-degenerate ice manifold remains as a stable low-temperature phase over a wide temperature window [6].
Under a [111] magnetic field, the pyrochlore spin ice exhibits a two-stage magnetization process [23]. The first plateau refers to kagome-ice phase with residual spin degeneracy, and the second plateau is a long-range ordered state which can be interpreted as a magnetic monopole condensate [7]. Neutron scattering experiments reveals the pinch point singularity of kagome-ice phase by applying magnetic field along [111] crystallographic direction on Ho 2 Ti 2 O 7 , yet an unusual critical scattering appears when the field is slightly tilted away from this direction [24]. By using the dipolar spin-ice model with combination of fields, it is shown that the critical scattering is ascribable to the long-range ordered q = X state, and the lack of Bragg peak may be attributed to partial ordering in the real material [25].
In this paper, we study the interplay between furtherneighbor interactions and tilted magnetic field in kagome ice in order to understand possible mechanisms that prevent ordering. Specifically, we study here a two-dimensional classical kagome-ice system with magnetic moments aligning along the local 111 axes originated from the pyrochlore lattice. Due to this canting, the moments can respond to both in-plane and out-of-plane fields, allowing us to explore the rich phase diagrams by tuning both the second-and thirdneighbor interactions, and the fields. In particular, at the boundary between different phases we find metastable winding snake domains emerge with local excitations of monopoleanti-monopole pairs. Unlike in the conventional spin ice, these excitations are immobilized and can not proliferate to destabilize the snake domains. This leads to observable signatures in the neutron experiments and should provide a new perspective on the interpretation of experimental data.
The rest of the paper is organized as follows: In Sec. II we introduce the model Hamiltonian and the geometry of canted moment and external fields. In Sec. III we present the equilib-  I and  III are more favorable than II and IV. rium Monte Carlo results and the phase diagrams of various ground states. To show the relaxation process from saturated state and the emergence of snake domains, we perform the rejection-free kinetic Monte Carlo simulation. The method and results are described in Sec. IV. Then, in Sec. V, we discuss the energetics of snake domain and the localized excitation on its periphery, along with possible experimental signatures in spin-ice crystals or artificial spin ices. We end with concluding remarks in the last section.

II. THE MODEL
We consider a spin model on the two-dimensional kagome lattice with Ising moments canted along the local easy axes derived from the pyrochlore ice ( Fig. 1), Each magnetic moment can be written as S i = µσ i d κ(i) where σ i = ±1 is the Ising variable and d κ(i) is the unit vector along local easy axis of the ith spin belonging to sublattice κ. Writing the Hamiltoinian in terms of the Ising variables σ i , In the following, we use the effective exchange interactions J 1 , J 2 , J 3a , and J 3b to characterize the phases, and the external magnetic field h is decomposed into out-of-plane and in-plane components h = h ⊥n + h k [26] as depicted in Fig. 1(a). Explicitly, the products of unit vectors and local easy axes read In the following, we set J 1 = 1.

III. GROUND STATES
With only the antiferromagnetic nearest-neighbor interaction J 1 > 0, the system has macroscopic degeneracy in the ground state. Each configuration in the ground-state manifold satisfies the ice rule on the kagome lattice: for both uppointing and down-pointing triangles, the three spins must be either two-in-one-out or one-in-two-out. We can define the topological charge for the kagome ice as [27,28] for up-pointing and down-pointing triangles, respectively. The kagome-ice rule in this charge language becomes Q , Q = ±1 ( Fig. 1(b)). Therefore, the nearest-neighbor Hamiltonian contains only self-energy of topological charges and describes a gas of single charges. The total charge density of the system is calculated through such that a single-(triple-) charge ordered state gives Q t = +1 (+3). L denotes the number of unit cells in each direction and the total number of spins N = 3 × L × L.
The degeneracy of single-charge configurations is lifted by an out-of-plane magnetic field h ⊥ , which acts as staggered chemical potential of triple charges ( Fig. 1(c)). Since each triangular unit contains two majority spins and one minority spin under ice rule, the configurations with majority spins aligned to the field become energetically more favorable in large perpendicular field. For h ⊥ larger than the critical value, the minority spin also tends to be aligned with it, leading to a saturated state of staggered triple charges ( Fig. 2(d)). These triple charges correspond to the magnetic monopoles and anti-monopoles in the pyrochlore spin ice in the [111] external field, and the saturation of charges represents the second plateau of the two-stage magnetization process described in Sec. I. Thus, the canted moment in our model allows us to imitate the magnetic field responses of the three-dimensional spin ice.
The projected component of canted spins on the kagome plane provides a way to perturb the kagome-ice manifold further through a small in-plane field h . The direction of this in-plane field is defined to point toward one vertex of the triangle in the kagome lattice, as shown in Fig. 1(a). Thus, for high-enough h , moments on this sublattice tend to align with h and the other spins are arranged into the local configurations that preserve the ice rule. We denote these fixed spins as kink spins and the others as stripe spins, as shown by different colors in Fig. 1(b). At low temperatures below J 1 , the states satisfying the ice-rule are characterized by polarized kink spins pointing downwards and rows of parallel stripe spins, either pointing to the left or to the right, which can be described by a local Z 2 order parameter. This regime is similar to the sliding phase discussed in other frustrated magnets. As T is further lowered, the second-neighbor coupling J 2 < 0 breaks the sliding symmetry and gives rise to an antiferromagnetic order of the parallel spin stripes. The resultant two-fold symmetric long-range ordered state is called the stripe state, as depicted in Fig. 2 For weaker in-plane fields, the kink spins are no longer pinned and a different long-range ordered state with √ 3 × √ 3 magnetic unit cell emerges. This so-called vortex state can be seen as an aggregation of hexagonal vortices with opposite chirality, as shown in Fig. 2(c). The same ground state can be found in both the equivalent short-range model with in-plane moments [29] and a long-range kagome-ice model with dipolar interactions [27,28]. Figure 2(a) and (b) shows the schematic phase diagrams of the ground states found in equilibrium Monte Carlo simulations. Besides the known ground states mentioned above, we found a new ground state emerging within a small domain of parameter space. The configuration of this state is depicted as Fig. 2(f), and can be seen as a vortex state with columnar defects of reversed kink spins, thus we call it the defected vortex state (d-vortex state). Closer inspection finds one-third of the triangles in this configuration contain triplecharge defects, and the others satisfy kagome-ice rule with staggered single charges. To our knowledge, this is a rare case that the long-range order of single and triple charges coexist.
The temperature and field-dependence of the charge density and various order parameters measured in the Monte Carlo simulations are presented in Fig. 3. In the J 1 -J 2 -h ⊥ model, the critical field can be simply calculated from the energy difference of vortex state and saturated state as h ⊥,c = 12J 1 + 12J 3a + 6J 3b from which we obtain h ⊥,c ∼ 11.96 (in the unit of J 1 ) for our chosen parameters, and it is consistent with the simulation results shown in Fig. 3(a)-(b). For the J 1 -J 2 -h ⊥h model (Fig. 3(c)), the critical field that drives the system into stripe state can be calculated from the energy difference of vortex state and stripe state: mula because both ground states satisfy ice rule. However, in the moderate fields, a mixed phase of vortex and stripe state is observed and numerous metastable states start to emerge and prevent the system from equilibriation. This emergent phase will be discussed extensively in the following sections.

IV. FIELD-QUENCH PROCESS
Since the transitions between the various phases in Fig. 2(a) and (b) are mostly first-order, mixed-phase states with complex structures often emerge at the phase boundaries especially at finite temperatures. Coupled with the underlying geometrical frustration, metastable states with unusually slow dynamics were observed particularly at the phase boundary between the stripe and the vortex states. Since the relaxation time scale of these metastable states is far beyond the capability of equilibrium Monte Carlo with conventional single-spinflip or loop algorithms [6], we applied the rejection-free algorithm called the waiting time method (WTM) to investigate the slow dynamics from given initial states [30]. Instead of us- ing acceptance ratio in the Metropolis algorithm, the waiting time method associates to each spin a waiting time as a local indicator of relaxation time, specifying how long a given spin has to wait before flipping. The flipping time for each spin is given by where the X i is a random number from the uniform distribution, X i ∈ (0, 1] and τ i is the waiting time from the Boltzmann weight, where ∆E i is the energy cost of flipping the i th spin. Therefore, the local flipping time t i is a stochastic variable from the exponential probability distribution associated with average τ i . The global time t g is assigned to the system and updated as t g = t i in each step after the i th spin with lowest t i is flipped. Then, the new waiting time of the flipped spin and its affected neighbors are calculated and added to the local flipping times where the prime in superscript denotes that the value is generated after flipping the target spin. To present the simulation results, the global time t g is used as the time scale of evolution. Each simulation run is executed at constant temperature and may suffer from metastability at low temperatures. Therefore, averaging measured observables over numerous runs is necessary for obtaining general dynamical behaviors. By using WTM, we study the J 1 -J 2 -J 3a -h ⊥ -h model on a L = 72 kagome lattice with 15552 spins. The fieldquench process is built up by assigning the initial configuration as saturated state that imitates the high-field condition (h ⊥ 11.96), and then quenching to a lower perpendicular field that favors the stripe state. The magnitudes of magnetic fields at t g ≥ 0 are chosen to be (h ⊥ , h ) = (11.76, −1.68), along with third-neighbor interaction J 3a = 0.05 owing to match the energetic phase diagram in Fig. 2(b). We found that it is enough to reach final stage of each simulation by 10 5 Monte Carlo steps, and the results are the average value over 100 samples.
In Fig. 4 we show the time evolution of the system in the field-quench scenario. In addition to charge density, we measured another two ad hoc order parameters to track the evolution of spin textures. The stripe-and kink-order parameters are defined as where η = ±1 for different sublattices of stripe spins, and ν = ±1 for odd and even rows of stripes. The symbols and denotes the stripe spins while ↓ represents kink spins.
As shown in the figure, the relaxation process from the saturated state can be divided into five stages. In stage I, the charge density drops drastically by rapid pair-annihilation of triple charges. The stripe ordering remains imperceptible since the emergent domains of single charges are randomly distributed. In the next stage, larger striped domains are constructed by further annihilation and diffusion of triple charges. However, in most simulation runs the system is stuck in stage III with partial stripe ordering and shows no further evolution. Only in very rare cases, the system in our simulation can escape from the long-lived metastable states and reach the stripe ground state in stage V.

V. THE SNAKE DOMAINS
The vortex state of the kagome ice is characterized by a broken Z 3 symmetry on top of the staggered ±1 charges. On the other hand, an Ising-like order parameter describes the anisotropic stripe phase. The incompatibility of these two order parameters indicates a first-order phase transition between these two ordered states. As discussed above, complex mixed-phase states often emerge at the first-order phase boundary. Here we show that the coexistence regime of the canted kagome ice exhibits an intriguing snake-like structure as well as localized "bound" states of monopole/antimonopole pairs at the edges of the snake domains. The appearance of this novel snake structure is related to the emergent gauge-like sliding symmetry of the stripe phase at finite temperatures. A typical structure of a snake domain is illustrated in Fig. 5(b). Due to the presence of further-neighbor interactions, a snake excitation is lower in energy than the string excitation ( Fig. 5(a)), and can be seen as a composite of single string and a series of hexagon-loop excitations. Since these hexagonal loops are non-winding loops and do not change the magnetization, snake and string excitation share the same topological sector and are indistinguishable by their winding numbers [31,32]. In contrast to the string excitation, the width of one snake domain can be any odd number greater than or equal to three. Moreover, the snakes with different widths are degenerate in energy. The energy change of a single snake from the background stripe state can be written as where L is the lattice dimension along one edge and ∆E s denotes the energy difference per unit length. Note that for the simulation result presented in Fig. 4, the parameters are chosen to be (J 2 , J 3a ) = (−1/3, 0.05) and h = −1.68. Therefore, the Zeeman energy of in-plane field is the only term that raises the energy from the stripe state. Since the energy increase of forming a snake from the stripe state is proportional to the vertical dimension L, it is expected that a perfect snake domain is unlikely to be seen in the thermodynamic limit. However, according to our simulation, a broken snake can only propagate through the migration of single monopoles at the end points. Thus, the kinetic bottleneck of moving one monopole along a highly anisotropic path results in the slow relaxation of these domains even in larger system sizes. At low temperatures, once a uniform snake domain is formed and wind through the periodic boundaries, spins become very robust against thermal fluctuation except for the kink spins on the edge. When a snake emerges from the stripe state, the ice-rule constraint forces the kink spins on the edge to align opposite to the in-plane field. Flipping those spins via thermal fluctuation can lower the Zeeman energy but simultaneously break the ice rule and create triple-charge pairs. In a large perpendicular field as in the WTM simulations, cancellation of spin-spin interaction and the Zeeman energy allows the edge spins to be more easily thermally flipped, resulting in isolated triple-charge pairs. Interestingly, although the edge spins are easily excited, these charges can hardly diffuse to any adjacent positions. In a sense, these triple charges are not only confined but also localized, partly attributed to the stable structure of the snake domains. The excitation energy of each edge spin can be shown explicitly as h . (11) With sufficiently large h or J 3b , all the edge spins can be excited and aligned with the in-plane field, as shown in Fig. 5(d).
For a metastable state with snake domains, there are two types of disorders. Firstly, both the number of the snake domains and their widths can fluctuate. This contributes to the destruction of the stripe order in our simulation. Secondly, although thermal fluctuation creates localized edge excitations, the total number of these edge spins is governed by the number of snakes. In general, for a given configuration of n snakes and m edge excitations as shown in Fig. 5(e), the total energy of this configuration is where n and m are non-negative integers with the upper bounds n ≤ L/3 and m ≤ nL. Two distinctive conditions are worth considering. For n = L/3 and m = 0, all the snake domains have minimum width and all edge spins are not excited. This fully compact configuration of snakes is exactly the vortex state. In a similar way, for n = L/3 and m = nL it is equivalent to a defected vortex state. Therefore, these metastable states can be regarded as mixture of the three ground states. Using the order parameters we measured in the simulations, we can calculate the exact values of n and m via As a remark, combining the two magnetizations M ⊥ and M one can also reconstruct the values of n and m, hence the ad hoc order parameters are dispensable for this task. The temperature and J 3a dependence of n and m are presented in Fig. 5(f). In almost all the cases, the number of excited edge spins increases monotonically as temperature rises, indicating the isolated and localized nature of these spins. In contrast, the number density of the snakes does not show the same trend and is changed by third-neighbor perturbation. On account of fixed lattice size in our simulations, the number density of snake is inversely proportional to the average size of snake-like domains. For J 3a = 0.1, the zero-temperature ground state is the vortex state ( Fig. 2(b)), and the domain structure is quite robust against thermal fluctuation, showing the saturated value in the number density of snakes (0.33) at all temperatures considered in Fig. 5(f). However, thermal fluctuation does excite the edge spins locally, leading to the monotonic increase in the number of excited edge spins.
On the other hand, for J 3a = −0.05, larger snake domains are easier to form since stripe state is the ground state. In this case, we notice that at intermediate temperatures thermal en-ergy can reduce the number density of the snakes and expand the domain width, moving the system closer to the stripe state. Finally, when J 3a is closer to the vortex-stripe phase boundary, the number density of the snakes shows a gradual increase with temperature, indicating that the nucleation process of the snakes is thermally enhanced. Therefore, although the stripe state is the zero-temperature ground state for J 3a = 0.05, at intermediate temperatures the domain structure is actually closer to the vortex state. At higher temperatures, the domain structure is eventually destroyed by thermal fluctuations.
Finally, we remark on the possible experimental signatures of the snake phase in kagome-ice systems. From the equilibrium Monte Carlo simulations, we calculate the magnetic structure factor for different phases at finite temperatures. On one hand, the Bragg peaks of vortex states appear at the corners of the Brillouin zone, indicating the six-fold symmetric configurations ( Fig. 6(a)). On the other hand, the presence of h breaks the rotational symmetry and the Bragg peaks only appear at the center of two edges with oval-shaped diffuse scattering ( Fig. 6(b)). Note that this similar pattern is observed in the neutron-scattering experiment of Ho 2 Ti 2 O 7 single crystal under tilted [111] magnetic field [24], and later explained as the presence of long-range ordered q = X state, which is the three-dimensional analogy of our stripe state [25]. However, the Bragg peaks of q = X state is not observed in the neutron experiment. For the snake phase in our model, the structure factor shows similar horizontally smeared scattering along the top and bottom edges, but the highest intensity no longer appears at the stripe position due to the presence of snake domains that breaks the periodicity of the stripe spins ( Fig. 6(c)). In one simulation, the points with highest intensity are located at somewhere between the Bragg peaks of two long-range ordered states, indicating the frozen topological sector of kagome-ice manifold.
In the artificial spin ice systems, thermally active kagomeice nanoarray of permalloy is successfully synthesized and studied recently. The equilibrium thermal behaviors such as phase transition from kagome-ice phase to long-range ordered vortex state is detected by low-energy muon spin relaxation (µSR) [33], while the non-equilibrium field-quench relaxation is also investigated through time-dependent spatial correlation of moments and charges [34]. Most importantly, through the X-ray photo-emission electron microscopy and magnetic force microscopy, both the magnetic domains and local triplecharge pairs can be imaged in the artificial kagome-ice system [34,35].

VI. CONCLUSION
In conclusion, we study a kagome-ice model with canted spins that respond to both in-plane and out-of-plane mag-netic fields. The inclusion of second-and third-neighbor spinspin interactions give rise to complex texture in both spin and charge degrees of freedom. The various competing phases at low temperatures lead to novel snake-like structure in the coexistence regime between the vortex and stripe phases. The combination of geometrical frustration and the anisotropic nature of these snake domains give rise to robust metastable states with extremely slow relaxation dynamics after the field quench. We provide detailed structural and energetic characterizations of these snake domains in the mixed-phase states. In spite of its long-lived nature, spins on the edge of each domain can be easily excited by thermal fluctuation, leading to frequent creation and annihilation of triple charges. These topological charges are localized since they can hardly overcome the energy barrier to migrate into the interior of the snake domain.
The structure of the mixed-phase states also plays an important role in the kinetics of the first-order phase transitions. For example, in the nucleation and growth scenario, the surface tension of domain walls enclosing the stable phase determines the size distribution of the nuclei as well as the nucleation speed. The elongated structure of the snake domains implies a rather anisotropic nucleation process as well as relaxation dynamics, which should be rather different from the isotropic case and is an interesting subject for future studies. From the experimental perspective, similar metastability and domain structure may be relevant to the absence of q = X Bragg peak in neutron scattering experiments [24]. Also, the recently discovered kagome-layered compound Dy 3 Mg 2 Sb 3 O 14 [14] confirms the charge-ordering transition of kagome ice and can potentially demonstrate the metastability and charge texture when external field is applied.