Multi-Step Ordering in Kagome and Square Artificial Spin Ice

We show that in colloidal models of artificial kagome and modified square ice systems, a variety of ordering and disordering regimes occur as a function of biasing field, temperature, and colloid-colloid interaction strength, including ordered monopole crystals, biased ice rule states, thermally induced ice rule ground states, biased triple states, and disordered states. We describe the lattice geometries and biasing field protocols that create the different states and explain the formation of the states in terms of sublattice switching thresholds. For a system prepared in a monopole lattice state, we show that a sequence of different orderings occurs for increasing temperature. Our results also explain several features observed in nanomagnetic artificial ice systems under an applied field.


Introduction
Spin ices have been extensively studied as ideal systems that exhibit geometric frustration effects since not all the pairwise spin interactions can be satisfied simultaneously [1,2,3,4]. Such systems are termed spin ices due to their similarity with a water ice phase in which proton ordering into corner-sharing tetrahedra is frustrated but obeys the "ice rule" of two protons "in" (close to the O atom) and two protons "out" (far from the O atom) in every tetrahedron [3]. In spin ices, the corresponding spin ice rules determine how many spins point toward or away from each vertex in the system, and these rules vary depending on the geometry of the system. For example, in two dimensional (2D) square ice, ice-rule-obeying states have two spins pointing toward each vertex and two spins pointing away, and it is possible for the spins to organize into an ordered ground state configuration [5,6]. For 2D kagome ice, ice-rule-obeying states have two spins in and one spin out out or two spins out and one spin in; in this case, the lowest energy state is not ordered but contains only vertices that obey these ice rules.
There has recently been growing interest in 2D artificial spin ice systems created using nanomagnetic arrays [5,6,7,8,9,10,11,12,13,14,15,16,17,18,19], 2D colloidal particle assemblies [20,21], or vortices in nanostructured superconductors [22]. In the nanomagnetic arrays, the orientation of the magnetic moment of each nanoisland produces the effective spin ordering. The colloidal and vortex systems more closely resemble the water ice system in that the repulsive interactions of the colloids or vortices at vertices resemble the interactions between charged protons, and the ice rules indicate the number of colloids or vortices that sit close to or away from a vertex. Experimental versions of these systems allow for direct imaging of the microscopic spin configurations as well as for careful control of numerous parameters that are not accessible in real spin systems.
In the work of Wang et al. [5] on 2D square ice, as the interactions between neighboring nanomagnets increased, the system was increasingly filled with ice-ruleobeying vertices; however, the predicted ground state was not observed due to quenched disorder effects in the nanomagnet array. Experiments on kagome artificial ice systems soon followed and showed that the ice rules predicted for that geometry were obeyed [14]. In more recent square ice studies using more advanced techniques, large regions of the sample adopted the square ice ground state while non-ice-rule obeying vertices created grain boundaries [9]. Such grain boundary formation was previously predicted in vortex simulations to occur for weak quenched disorder, while for stronger disorder isolated non-ice-rule defects begin to appear [22].
Lack of thermalization is an issue for the experimental nanomagnet systems since thermal activation of the spin configurations is weak or nonexistent. This can be partially mitigated by applying a changing external field [7,8]. It was recently shown that biased states appear under applied external fields for both square and kagome ice systems, while cycled external fields can be used to generate hysteresis curves [12,11,10,13]. For kagome ice, Zabel et al. [11] observed an ordered monopole state where each vertex has either all spins pointing out or all spins pointing in, with the vertices alternating sign in a crystalline arrangement. Other studies have examined the creation and motion of monopole excitations in the presence of a biasing field for both square and kagome ice systems [12,10]. Recent theoretical studies have found that 2D kagome ice can undergo a two-stage ordering transition from a paramagnetic state to an ice-rule-obeying state followed by an ordered monopole state [23]. The thermal disordering or melting of the ordered states can be readily studied in colloidal artificial ices, and recently materials have been created with a spin ordering temperature that is near room temperature, making it now possible to melt magnetic artificial spin ices [24].
Here we show that a variety of different types of orderings are possible in both kagome and modified square ice systems as a function of external field, temperature, and interaction strength for a colloidal artificial spin ice system composed of arrays of elongated optical traps that each capture one colloid. Experimental studies have demonstrated frustration in buckled 2D colloidal layers [21] while numerical studies have shown that charged colloids in elongated trap arrays exhibit the same spin ice rules observed for nanomagnet arrays [20]. Computational and experimental studies of colloidal ordering on 2D square and triangular substrates have revealed numerous ordered states and multi-step transitions between the states [25,26,27,28,29,30,31,32]. The elongated optical traps we consider here permit each colloid to sit in one of two locations at either end of the trap. Arrays of this type of trap have been created experimentally to study various types of ratchet effects, and preferential ordering of the colloids into a state that minimized the colloid-colloid interactions was observed [33]. An advantage of the colloidal system over the nanomagnetic systems is that thermal activation of the effective spin degrees of freedom is possible. The relative importance of the thermal fluctuations can be varied either by adjusting the temperature or by holding the temperature fixed and modifying the strength of the colloid-colloid interactions. It is also possible to bias the colloidal artificial ice system using an electric field.
Simulation-We consider a 2D array of N elongated traps that each contain a single colloid. Each trap has two potential minima where the colloid can sit, as shown in figure  1(a) where we illustrate a single vertex of the square ice system. The effective spin at each site is defined to point toward the end of the trap occupied by the colloid. The colloids interact with each other via a repulsive screened Coulomb or Yukawa potential given by V (r ij ) = r −1 ij exp(−κr ij )r ij where κ = 4/a 0 is the screening length and a 0 is the average spacing between the vertices. Since the colloids are repulsive it is very energetically costly for all the colloids around a vertex to sit close to the vertex in a positive monopole state; on the other hand, the lowest energy state for a vertex is the negative monopole configuration in which all the colloids sit away from the vertex. Due to geometrical constraints it is not possible for all of the vertices to adopt negative monopole configurations. The most cost-effective way to accommodate the colloidcolloid repulsion in the square ice is for each vertex to have two close colloids and two far colloids in a "two in-two out" state that obeys the square ice rules. These ice rule obeying states can be biased, with the close colloids occupying traps oriented at 90 • to each other, or they can be the slightly lower energy ground state configuration shown in figure 1(a) [20]. For kagome ice, shown schematically in figure 1(b), the ice rule obeying vertices have either two in and one out or one out and two in. The colloid dynamics evolve according to the following overdamped equation of motion: Here η is the damping constant and the colloid-colloid interaction force is F cc is the position of particle i(j), Z * is the unit of charge, ǫ is the dielectric constant of the solvent, and q is the magnitude of the charge on a single colloid. The strength of the repulsion between the colloids can be controlled by varying q. The thermal force F T arises from Langevin kicks with F T The substrate force F s i arises from N traps composed of two parabolic ends capping a cylindrical confining area of length l = 1.333a 0 and width d p = 0.4a 0 with a maximum strength of F p and radius r p ; an additional parabolic barrier of height f r is placed at the center of the trap to produce two potential minima at each end of the trap [20]. The external biasing force is given by F ext = F ext [cos(θ ext )x + sin(θ ext )ŷ], with θ ext = 0 for kagome ice and θ ext = 45 • for square ice. For the nanomagnetic system, an in-plane applied external field was used as a biasing force that could align the magnetic moments of the nanoislands. For the colloidal system, the sample can be biased by an in-plane electric field, while for a system of vortices in a type-II superconductor the bias would come from an in-plane applied current. In this work we consider external forces that are strong enough to induce hopping over the central barrier of the traps but not strong enough to cause the colloids to escape from the traps.

Orderings for Kagome Ice
We first consider a kagome ice array in an external field F ext applied along the x axis, in alignment with the axes of 1/3 of the traps in the sample. In figure 2 we plot the fraction of vertex types N i /N, the effective magnetization M, and the applied external force F ext as a function of time for a sample with noninteracting colloids at q 2 = 0. The vertex types are defined as follows: N 0 has no colloids near the vertex; N 1 has one colloid near the vertex and two far from the vertex; N 2 has two colloids near the vertex and one far; and N 3 has three colloids near the vertex. N 1 and N 2 vertices obey the kagome ice rules, while N 0 and N 3 vertices are negative and positive monopoles, respectively. In figure 2(c), the noninteracting colloids begin in a randomized state with N 0 /N = N 3 /N = 0.12 and N 1 /N = N 2 /N = 0.38. As shown in figure 2(b), F ext is gradually increased from zero. Just above F ext = 1.0, the system jumps into the biased ice-rule obeying state illustrated in figure 3(a). We define the contribution of each trap to the effective magnetization M according to the component of the effective spin direction that is aligned with F ext . Each trap is assigned a value s i = ±1 depending on whether After the sample is fully biased, we decrease F ext as shown in figure 2(b). At F ext = −1.7, the biased state switches into an ordered monopole state where N 0 and N 3 each jump from 0% of the population to 50% of the population while N 1 and N 2 simultaneously drop to zero as shown in figure 2(c). The ordered monopole state, illustrated in figure 3(b), is a lattice of positive and negative monopoles and is the same state observed in experiments on kagome nanomagnetic systems under external drives [11]. The ordered monopole state arises due to the fact that the kagome lattice is effectively composed of two sublattices that have different switching fields under an external x-direction drive. This is illustrated schematically in figure 4 where the darker ellipses indicate the sublattice with a lower switching field, which we term sublattice A. For the case of non-interacting colloids with q 2 = 0, the switching field for sublattice A is simply the maximum force f r of the barrier in the center of the trap. Since the applied drive is not aligned in the direction of the traps of sublattice B, the switching force for sublattice B has a higher value of f r / cos(60 • ). When we apply a sufficiently large F ext in the positive x-direction to cause both the A and B sublattices to switch, the system enters the positively biased state illustrated in figure 4(a). After we reduce F ext back to zero, we begin applying F ext in the negative x direction. Sublattice A switches in the negative x direction first due to its lower threshold, but since the colloids in sublattice B have not yet switched, the ordered monopole state forms as shown in figure 4(b). For q 2 = 0 or finite colloid-colloid interactions, creation of monopoles becomes energetically unfavorable, so for sufficiently high q 2 the switching of sublattice A simultaneously induces a switching of sublattice B, bringing the system directly into the negatively biased state illustrated in figure 4(c) and preventing the formation of the ordered monopole state. For a certain range of nonzero q 2 , however, the ordered monopole state can still be stabilized when F ext passes through zero. Our results suggest that in the experiments of Zabel et al. [11] the magnetic islands are in the weakly interacting or noninteracting limit.
We next examine the thermal disordering of the ordered monopole state for varied colloid-colloid interaction strengths. In figure 5(a) we plot the vertex populations N i /N versus temperature T for a noninteracting sample with q 2 = 0, while figure 5(b) shows N i /N versus T for an interacting sample with q 2 = 0.3. In each case we set F ext = 0 before applying the temperature. In the noninteracting system in figure 5(a), the ordered monopole state persists up to T = 2.5, at which point the N i /N cross over to the values expected for an uncorrelated random arrangement or paramagnetic state. In figure 5(b) for q 2 = 0.3, N 3 /N and N 0 /N drop to zero with increasing T while N 1 /N and N 2 /N increase to 0.5, indicating a crossover from the ordered monopole state into an ice-rule obeying state. This is not the biased ordered state illustrated in figure 3(a); instead, the ice-rule obeying state lacks true long range order, as shown in figure 3(c). As T is further increased, N 3 /N and N 0 /N gradually increase while N 1 /N and N 2 /N decrease to values close to those expected for a random or paramagnetic state of the type illustrated in figure 3(d).
By conducting a series of simulations, we map out the different orderings as a function of temperature and pairwise interactions as shown in figure 6. The line separating the ice rule obeying state from the random state is defined as the point at which N 1 /N and N 2 /N reach values that are within 10% of those expected for a random vertex arrangement. For weak interactions or small q 2 the ordered monopole state melts directly into the random state. As q 2 increases, the ice-rule obeying region increases in extent while the line marking the crossover between the ordered monopole state and the ice-rule obeying state drops to lower T . These results suggest that the ordered monopole states can only be observed for systems with effective spin-spin interactions that are weak compared to the strength of the switching or coercive field. For q 2 > 0.4 the ordered monopole state can no longer be prepared by means of the external field protocol since the repulsion between colloids becomes strong enough that switching of the A sublattice immediately induces switching of the B sublattice. Nevertheless, the ordered monopole state is still stable at lower T for q 2 > 0.4 as can be seen by artificially preparing the lattice in the ordered monopole state and then increasing the temperature.  . Phase diagram for T vs pairwise interaction strength q 2 for the kagome ice system with f r = 1 indicating regions where the ordered monopole, ice-rule obeying, and random or paramagnetic states appear. For weak interactions (small q 2 ) the monopole state disorders directly into the random state, while the width of the icerule obeying region grows with increasing q 2 . Red symbols: beyond the highest value of q 2 where the ordered monopole state can be created with a biased drive protocol, the system can be artificially placed into the ordered monopole state and then melted into the ice-rule obeying state.
Such artificial preparation could be achived in colloidal systems by, for example, tuning the colloid-colloid interactions to a lower value in order to create the ordered monopole state and then increasing the interactions to a higher value; however, it is generally not possible to vary the interaction strength in magnetic systems. The line marking the melting into the ice-rule obeying state obtained by using artificially created ordered monopole states is marked with a dashed line in figure 6.

Orderings For Square Ice
In order to create an easily accessible ordered monopole state for a square ice lattice, we propose a simple extension of the standard square ice geometry where we now have two sublattices as illustrated in figure 7. Each sublattice has a different strength of the central barrier, f 1 r and f 2 r , and we take f 2 r /f 1 r = 2. In figure 7(a), the darker traps have the weaker barriers. Such a geometry could be produced in the nanomagnetic system by using two different sizes or two different materials for the nanoislands, producing two different coercive field values. We apply a biasing force along a line oriented at 45 • to the x axis. In the square ice array, the vertex types N 0 , N 1 , N 2 , N 3 , and N 4 have 0, 1, 2, 3, and 4 colloids close to the vertex, respectively. The N 2 vertices are further    figure 7(c). In analogy to the kagome ice system, in the two-sublattice square ice sample the monopole state forms when sublattice A, with its weaker barrier, switches before sublattice B. If the colloid-colloid interaction strength is small enough, the particles in sublattice B remain in their original positions, producing the monopole state shown in figure 7(c). As F ext is further increased in the negative direction, sublattice B switches at F ext = −1.3 to form a negatively biased ice rule obeying state aligned with the drive. We note that in our previous work on square ice with only a single trap sublattice, monopole formation was rare due to its high energetic cost. We next consider thermal effects on the ordered monopole state. To prepare the sample, we sweep F ext to a value at which the monopole state appears and then switch off F ext . If the pairwise interactions between colloids are not too strong, the monopole state remains stable even without a biasing drive. Figure 9(a) shows the vertex populations N i /N versus temperature T for a system of noninteracting colloids with q 2 = 0. Here the monopole state illustrated in figure 10(a) disorders into the paramagnetic state illustrated in figure 10(e), with N 0 /N and N 4 /N monotonically dropping from 0.5 to around 0.06 but never falling below this value. In samples with interacting colloids, such as the system with q 2 = 0.2 shown in figure 9(b), N 0 /N and N 4 /N both drop from 0.5 to 0 with increasing T indicating that all monopoles disappear from the system. As T is increased further, N 0 /N and N 4 /N gradually increase back to the paramagnetic value of 0.0625. For 1.9 ≤ T ≤ 2.3, where N 0 /N and N 4 /N are dropping to zero, the values of N 1 /N and N 3 /N pass through a peak. This is followed by a window 2.3 ≤ T ≤ 3.3 where N 1 /N and N 3 /N drop again while the value of N b 2 /N peaks. A weak peak in N gs 2 /N appears in the vicinity of T ∼ 3.5, while at higher temperatures the vertex populations approach the values expected in a completely random configuration.
The same trends are more clearly evident in figure 9(c) for a sample with larger pairwise interactions, q 2 = 0.4. Here the highest value reached by N 1 /N +N 3 /N is nearly 0.9, indicating that large portions of the sample are filled with N 1 or N 3 vertices. The corresponding colloid configuration is illustrated in figure 10(b), where local ordering of the N 1 and N 3 vertices into a lattice can be seen. We term this state a biased triple state; it forms due to the preferential switching of colloids in the weaker sublattice A traps. As T increases, N 1 /N and N 3 /N drop while N b 2 /N peaks at a value of 0.8 corresponding to the thermally induced biased ice rule obeying state, illustrated in figure  10(c). At higher T , N b 2 /N decreases again while N gs 2 /N reaches its peak value of 0.7 corresponding to the thermally induced ground state shown in figure 10(d). For still higher T , the N i /N approach the values expected for a random configuration and the sample enters a paramagnetic state, illustrated in figure 10(e).
The sequence of ordered states we observe for finite pairwise interactions under finite temperatures can be understood by considering the interaction of the metastable monopole state with the two trap sublattices. Our protocol of sweeping F ext into the monopole regime and then switching off F ext at zero temperature causes the sample to be trapped in this metastable state for low temperatures. As T increases from zero, the system first lowers its energy by moving one colloid out of each N 4 vertex into an N 0 vertex, creating N 1 and N 3 vertex pairs and generating the biased triple state. This produces the peak in N 1 /N and N 3 /N near T = 1.6 in figure 9(c), while figure 10(b) indicates that due to the finite temperature, some non-N 1 or N 3 vertices still exist in the sample. In the biased triple state, the lower energetic cost of an N 3 vertex compared to an N 4 vertex means that a higher thermal activation energy is required to destroy N 3 vertices compared to N 4 vertices. The difference in thermal activation energies determines the size of the temperature window where the biased triple state remains metastable. Once the temperature is large enough, colloids jump out of N 3 vertices into neighboring N 1 vertices, resulting in the formation of N 2 vertices. This hopping is more likely to occur in the weaker traps of sublattice A, resulting in the preferential formation of biased N b 2 vertices. As a result, N b 2 /N peaks near T = 2.3 in figure 9(c) and the sample enters the biased ice rule obeying configuration (with thermal defects) shown in figure 10(c). Since the N b 2 vertices have a slightly higher energy than the ground state N gs 2 vertices [20], at even higher temperatures additional hopping produces an increased fraction of N gs 2 vertices, placing the system in the thermally induced and thermally defected square ice ground state configuration illustrated in figure 10(d). For sufficiently high T , the correlations between neighboring traps are lost and the sample enters the paramagnetic state shown in figure 10(e).
In figure 11 we plot a phase diagram of the regions in T and q 2 space where the different ordered states appear. The lines indicate crossovers and are determined based on when the different N i /N cross threshold values. As q 2 increases, the destruction of the ordered monopole state drops to lower values of T . The biased triple state and biased ice rule obeying state maintain roughly constant widths in T as q 2 increases, although both phases shift to lower ranges of T . In contrast, the square ice ground state, which appears only for larger values of q 2 , increases in width as q 2 increases. Above q 2 = 0.45  Figure 11. Phase diagram for T vs q 2 for the modified square ice sample with f 1 r = 1.0 and f 2 r = 2.0 indicating regions where the ordered monopole, biased triple, biased ice rule obeying, square ice ground state, and random states appear. As q 2 increases, the width of the ordered monopole state decreases while the transition to the random state shifts to higher T .
we can no longer stabilize the ordered monopole state using our biasing protocol. Our results indicate that a simple modification of the square ice geometry can produce a wealth of orderings in the artificial spin ice system.
We note that the phase diagrams for the kagome and modified square ice presented here do not exhibit true phases since we have driven the system into a metastable state. In many of the nanomagnetic systems that exhibit hysteresis, most of the observed states are also metastable unless a successful annealing protocol has been applied. All of our samples contained no quenched disorder. We expect that many of the states we observe would be gradually washed out if increasing amounts of quenched disorder are added; however, there should still be observable correlations in the different states that can be detected by analyzing the vertex population densities.

Summary
We have described the creation of different types of orderings in artificial kagome and modified square ice systems using external fields, temperature, and particle-particle interaction strength. We show that multiple-step ordering-disordering transitions occur that can be identified by counting the vertex populations. In the kagome ice, an ordered monopole state can be induced by applying a biasing field to samples with pairwise interactions that are not too strong. With increasing temperature, the system crosses into an ice rule obeying state and then into a high temperature paramagnetic state. For stronger pairwise interactions, the monopole state is unstable to the formation of the ice rule obeying state. Our results suggest that recent experimental observations of monopole ordering in kagome nanomagnetic systems were performed in a weakly interacting or noninteracting regime. For square ice, we make a simple modification to the geometry by introducing two sublattices of traps that have different switching fields and show that under the application of an external biasing field, this system can form a checkerboard ordered monopole state. The modified square ice exhibits multi-state ordering as a function of temperature starting from the monopole state, passing through a biased triple state, a thermally induced biased ice rule obeying state, a thermally induced square ice ground state, and a disordered or paramagnetic state. Although our results are obtained specifically for colloids on periodic trap arrays, the behavior should be generic to other artificial ice systems under external drives, including nanomagnetic array systems.