Colloquium: Ice Rule and Emergent Frustration in Particle Ice and Beyond

Antonio Ortiz-Ambriz, 2 Cristiano Nisoli, Charles Reichhardt, Cynthia J. O. Reichhardt, and Pietro Tierno 2, 4 Departament de F́ısica de la Matèria Condensada Universitat de Barcelona, Barcelona, Spain. Institut de Nanociència i Nanotecnologia (INUB) Universitat de Barcelona, Barcelona, Spain. Theoretical Division and Center for Nonlinear Studies Los Alamos National Laboratory Los Alamos, New Mexico 87545 USA Universitat de Barcelona Institute of Complex Systems (UBICS) Universitat de Barcelona, Barcelona, Spain.∗


I. INTRODUCTION
Frustration in life emerges with the impossibility of simultaneously satisfying a set of requirements. Frustration in physics is not very different. A classical example is that of three spins on the vertex of a triangle that want to be antiferromagnetically aligned (Wannier, 1950). This requirement cannot be realized all around the triangle, so at least two spins will display ferromagnetic order and will generate one frustrated bond. More generally, a geometrically frustrated system is one that is subjected to local requirements that cannot be satisfied collectively along certain loops in the system. The notion is therefore intrinsically topological, i.e. invariant under transformations that do not rip those frustrated loops, and indeed it lends itself to abstract, elegant treatments in term of Wilson loops and gauge symmetry, at least in the case of spin systems (Fradkin et al., 1978). In real systems, however, the local requirement is usually the minimization of an energy as a pairwise interaction that is in general geometry-dependent.
Frustration is essential for the understanding of a variety of real materials, such as spin glasses (Mydosh, 2014), water ice (Bernal and Fowler, 1933;Pauling, 1935), spin ices (Bramwell and Gingras, 2001;Harris et al., 1997;Ramirez et al., 1999), and even systems at different length scales such as granular materials (Richard et al., 2005), liquid crystals (Kamien and Selinger, 2001;Lopez-Leon et al., 2011), filament bundles (Hall et al., 2016), coupled lasers (Nixon et al., 2013) and many others (Grason, 2016). There, it is often-but by no means alwaysassociated with slow relaxation, degeneracy, and zero temperature entropy. Indeed, frustration implies compromise, and therefore produces various forms of socalled constrained disorder, which are manifolds whose disorder obeys some non trivial rules, either local or global. The ice rule is an important example of a local rule that constrains disorder and that can have global, i.e. topological, implications.
Typically in an ensemble characterized by constrained disorder, violations of local rules appear as localized excitations of the low-energy states of the system. These local rule violations control the collective dynamics. An example, described further below, is provided by the magnetic monopoles in spin ices (Castelnovo et al., 2008;Ryzhkin, 2005), which are violations of an ice manifold. Another example are the emergent topological charges in certain artificial spin ices (Lao et al., 2018). These concepts are of both practical interest, such as for the measurement of currents of magnetic monopoles (Giblin et al., 2011), and also of theoretical interest. Indeed, we are accustomed to understanding topological defects in terms of alterations of an underlying order, such as misplaced books on a library shelf or dislocations in a crystal, rather than in terms of disorder.
Manifolds of constrained disorder can be rich playgrounds for new physics, and frustration is a fundamental ingredient in the design of artificial systems that can generate novel exotic behaviors which are often not found in natural materials (Heyderman and Stamps, 2013;Nisoli et al., 2013;Wang et al., 2006). Artificial spin ice systems were first created by mimicking natural frustrated geometries and by realizing celebrated models of statistical mechanics (Baxter, 1982;Lieb, 1967b) in settings that allowed characterization at the constituent level, often in real time. However, since artificial materials can be realized in various geometries, a more recent effort (Morrison et al., 2013;Nisoli et al., 2017) has advanced the design of new systems generating a wide variety of new phenomena, including dimensionality reduction, emergent classical topological order, realizations of Pott's models, phase transitions, ice rule fragility, and quasi-crystal spin ices (Barrows et al., 2019;Gilbert et al., 2014Gilbert et al., , 2016aGliga et al., 2017;Lao et al., 2018;Louis et al., 2018;Ma et al., 2016;Östman et al., 2017;Perrin et al., 2016;Shi et al., 2018;Sklenar et al., 2019). Furthermore, many of these ideas proved to be exportable across different platforms, from nanomagnets to trapped colloids, to liquid crystals, and to superconductors (Duzgun and Nisoli, 2019;Latimer et al., 2013;Libál et al., 2009;Ortiz-Ambriz and Tierno, 2016;Wang et al., 2018).

II. THE ICE RULE
We start by briefly recalling the history of the ice rule, and its appearance in magnetic systems either natural or lithographically fabricated. For a more extensive treatment we refer the reader to available reviews and commentaries on the subject: (Bramwell and Gingras, 2001;Gardner et al., 2010;Gilbert et al., 2016b;Heyderman and Stamps, 2013;Nisoli, 2018a;Nisoli et al., 2013;Ramirez, 1994).

From water ice to spin ice
Many fundamental properties of water are still not completely understood (Chaplin, 2006). One of water's early mysteries pertained to its residual entropy and was solved by Linus Pauling in the 1930's. Through a series of carefully conducted calorimetric experiments, Giaque and Ashley (Giauque and Ashley, 1933;Stout and Giauque, 1936) had found that the entropy of ice at low temperature was not zero. Pauling (Pauling, 1935) explained it via the ice rule of Bernal and Fowler (Bernal and Fowler, 1933). In water ice each oxygen atom is bound to four others via hydrogen atoms. Two of these hydrogen atoms are close to the oxygen atom in covalent bonds, and two are further away forming covalent bonds with neighboring oxygens, Fig. 1(a), left. The number of ways in which the hydrogen atoms can be arranged, Pauling showed, grows exponentially with the number of oxygen atoms, leading to a non-zero entropy per molecule. Ice-like systems received renewed interest in the 1990's with the discovery a class of magnetic substances which were named "spin ices" because their magnetic texture at low temperature presents a degenerate low energy state with a residual entropy, consistent with the ice rule model. These rare earth titanates, such as Ho 2 Ti 2 O 7 and Dy 2 Ti 2 O 7 , presented frustrated interactions between the magnetic moments of their constituent at low temperature. Their magnetic cations Ho 3+ and Dy 3+ were known to carry a very large magnetic moment, about ten times the Bohr magneton. At low temperatures such moments could be regarded as binary, classical Ising spins constrained to point along the directions of the lattice bonds forming the pyrochlore lattice, as shown in Fig. 1(a), right. It was then noted (Harris et al., 1997) and confirmed experimentally (Ramirez et al., 1999) that the resulting ferromagnetic interaction favors the ice rule, in which two spins point in each vertex and two out, in analogy with the allocation of protons in the lattice of water ice.
FIG. 1 (a) Ice rule for the water ice I h (left) and magnetic spin ice material (right). The first is characterized by oxygen atoms sharing protons and located at the vertices of a diamond lattice. The spin ice has Ising-like moments pointing along the same lattice, from Castelnovo et al. (2008). (b) A projection of the three-dimensional pyrochlore lattice onto a two-dimensional plane gives rise to a square lattice, from Glaetzle et al. (2014). (c) Left: Magnetic force microscope image of an artificial spin ice composed by a square lattice of permalloy islands with lattice spacing of 400nm. Right: corresponding spin orientations. Bottom: vertex configurations with ice rule highlighted in the violet box. Image edited with permission from (Wang et al., 2006).

Artificial spin ice systems
Between the late 1990's and the early 2000's, a mature effort in the exploration of the magnetic state of nanodots and nanoislands (Bader, 2006) focused on how to obtain exotic magnetic textures and transitions. Meanwhile, a different approach to exotic behavior involved relatively simple elongated, single-domain, magnetic nanoislands, whose magnetization could be described by an Ising pseudo-spin. In this case, complexity is introduced via the mutual interaction of these nanoislands in order to generate possibly interesting collective states of exotic emergent behaviors (Heyderman and Stamps, 2013;Libál et al., 2006;Nisoli et al., 2013;Tanaka et al., 2006;Wang et al., 2006). This approach provided a twofold advantage. First, characterization of the individual degrees of freedom in real space is possible via methods such as Magnetic Force Microscopy (MFM), Photoelectron emission microscopy (PEEM), Transmission Electron Microscopy (TEM), Surface Magneto-Optic Kerr Effect (MOKE), and Lorentz Microscopy. Secondly, the collective behavior of these systems is open to design.

Ice rule and topology: conceptual themes
The ice-rule appears an innocuous enough concept, yet it can be understood in most general terms and has profound implications. Consider a lattice or even a directed graph, with binary variables such as Ising spins on each edge, impinging in vertices of various coordination z. The pyrochlore or square geometries introduced before can serve as examples. Then we can define the topological charge of a vertex of coordination z with n spins pointing toward as which corresponds to the difference between the number of spins pointing in and out. This notion is properly topological as it only depends on the topology of the graph. Furthermore, given a spin configuration, any single spin flip will alter the charge in two nearby vertices. Only flipping proper loops of spins will preserve the charge distribution.
In this language, the ice rule can be considered a prescription for a local minimization of |q| at each vertex. In practical systems this is typically (but not necessarily) enforced by the nearest neighboring spin-spin interactions. Any configuration of spins that locally minimizes |q| is said to obey the ice rule, and the subset of the phase space corresponding to such configurations is called an ice manifold. For a lattice of uniform, even coordination, the ice manifold is then characterized by zero charge q = 0 on each vertex. The first violation of the ice rule then corresponds to charges q = ±2, called monopoles. Very often-but not always-the ice manifold represents the lowest energy of the system, and then the ground state is degenerate and disordered, and has a residual entropy. Importantly it cannot be "explored from within." Any single spin flip on an ice-rule configuration creates a pair of monopoles of opposite charge ±2 (Castelnovo et al., 2008;Ryzhkin, 2005), violating the ice rule. Only the coherent flipping of a loop of spins, properly chosen so that they are all arranged head to tail, represents an "update" that does not violate the ice rule.
The simplest, square version of such system has motivated theoretical research in applied mathematics for half a century. In the 1960s, Lieb, Wu, Baxter, Rys, and others began working on simplified, two-dimensional models of mathematical physics, known as vertex models, that captured and also generalized the properties of ice. In these models, different energies are assigned to different vertex configurations on a square lattice, and in many cases the models can be solved exactly, typically via transfer matrix methods (Baxter, 1982;Lieb, 1967a,b;Rys, 1963;Wu, 1969). The six-vertex model in particular (Lieb, 1967b) only admits ice rule obeying vertices on a square lattice, and was meant to represent a solvable, two-dimensional equivalent of water ice. As it forbids monopoles, it is completely embedded into the ice manifold and has strong topological properties. For instance, if the degeneracy of the ice-manifold is lifted by an energetics that selects the antiferromagnetic state, as in the Rys F-model (Rys, 1963), then the corresponding ordering transition is infinitely continuous (Lieb, 1967a), and with an order parameter, also infinitely continuous! (Baxter, 1982).
Beside clarifying the topological nature of the ice rule, these models were also often equivalent to other important statistical mechanics systems, such as dimer cover models, and thus initiated an independent theoretical line of research in mathematical physics that has further evolved in terms of loop representations to describe topological effects at (and of) the boundaries. We will not report here on this half a century long, very interesting developments because it goes beyond the scope of this Review. Indeed the interesting phenomenology in those models arises from a topological structure-given by the ice rule-that in simulations and experiments is always violated. While it is true that such topological structure is present in the ground state, this does not imply that it extrapolates fully to the low energy physics, which is in fact generally a dynamics of topological defects. In a sense, realistic spin ice systems are neither topologically constrained nor unconstrained. They should perhaps be called "Topology Breaking" systems because for any non-zero temperature, violations of the topological constraints (the ice-rule) in form of monopoles substantially changes the physics. For a dramatic example: while the antiferromagnetic ground state of artificial square ice (Morgan et al., 2011;Nisoli et al., 2010;Porro et al., 2013;Wang et al., 2006;Zhang et al., 2013) would seem to be well described by the Rys F-model, the ordering transition of the latter is in the Kosterlitz Thouless class (Kosterlitz and Thouless, 1973;Lieb, 1967a), whereas in "real" square ice the transition is simply second order in the Ising class (Levis et al., 2013;Sendetskyi et al., 2019;Wu, 1969) precisely because monopoles break the topological constraint.
Thus in real spin ice systems it becomes less interesting to eviscerate all the possible topological representations of their idealized ice manifold. It is instead more interesting to see how the topological properties of an ultimately unreachable ground state affect the low-energy physics as the system breaks those topological constraints. This issue has been attacked via the concepts of spin fractionalization into magnetic charges (Castelnovo et al., 2008;Ryzhkin, 2005) and of spin fragmentation (Petit et al., 2016) into a Coulomb (Henley, 2010b) and non-Coulomb magnetization. For instance, in pyrochlore ice the ice manifold can be considered as a classical topological phase of constrained disorder (Castelnovo et al., 2012;Henley, 2010a,b). Such a phase is labeled not by an order parameter, as an ordered phase would, but rather by a fluctuating, solenoidal gauge field M (x), which can be thought of as a coarse grain of the spin magnetization. Then spin fractionalization implies that the low-energy manifold can be described in terms of local violations of its solenoidal nature, such that Thus, the ice manifold corresponds to q = 0 everywhere, and it is easy to show (Henley, 2010b) that it implies algebraic (specifically, dipolar) correlations in the reciprocal space, and thus the observed pinch points in the structure factor of neutron scattering. While this ice phase is critical, this is a purely mathematical abstraction since no transition to the ice manifold exists, and at non-zero temperature monopoles will always be present, however sparsely, and provide a correlation length. The low energy of the system remains however reminiscent of the topological nature of its ground state. Through spin fragmentation (Petit et al., 2016) it is possible to show that the low-energy ensemble can be decomposed into an ice-rule ensemble plus an crystal-charge ensemble. In terms of the coarse grained magnetization this merely corresponds to an Helmholtz decomposition of the magnetization vector in the sum of a divergence free and divergence full part. Because the two components are uncoupled at the lowest order in the effective free energy, dipolar correlations and pinch points survive at non-zero temperature. Artificial realization of topological phases was achieved in the Shakti geometry (Lao et al., 2018) as well as in magnetic square ice (Möller and Moessner, 2006;Perrin et al., 2016) and rectangular (Loreto et al., 2019;Ribeiro et al., 2017) ice. No equivalent phase has yet been realized with colloids, though those magnetic realization provide directions.
For vertices of odd coordination, the cancellation of the charge in Eq. (1) is impossible and |q| is minimized when q = ±1. Thus in a simple system of uniform, odd coordination, such as the kagomé ice (Qi et al., 2008), the ice manifold corresponds to ensembles of disordered spins obeying a 2-in/1-out or 2-out/1-in ice-rule at each vertex. Because there is no charge cancellation, such systems can be regarded as disordered, neutral plasma of opposite topological charges which, depending on realization, can interact, leading to further phases within the ice manifold (Chern et al., 2011(Chern et al., , 2013aDrisko et al., 2015;Möller and Moessner, 2009;Zhang et al., 2013), as we will show below in more detail.
Finally, systems of mixed coordination, say for instance z = 3, 4 host both versions of the ice rule, 2-in/2-out on z = 4 coordinated vertices and 2-in/1-out or 2-out/1 on z = 3 vertices. Here is where the difference between magnetic and particle-based spin ice becomes more dramatic, with the ice-rule breaking down in the latter, as we shall see in Section IV.

Colloids as a model system
Colloidal particles represent a versatile model system to explore collective phenomena in statistical physics and condensed matter, due to their accessible length-scales and to the presence of simple and tunable interactions (Poon, 2004;Yethiraj and van Blaaderen, 2003). Colloidal systems have provided new insight into the glass transition (Weeks, 2017), yielding phenomena (Schall et al., 2007), and the motion of dislocations (Schall, 2004), and have been used to test numerous foundations of both equilibrium (Thorneywork et al., 2017;Zahn et al., 1999) and nonequilibrium statistical mechanics (Martínez et al., 2016(Martínez et al., , 2017. One of the best examples of the use of colloids for both condensed matter and statistical physics studies is in understanding ordering and commensuration effects on surfaces, where the substrate can be one-dimensional (1D) (Bechinger et al., 2001), two-dimensional (2D) (Bohlein et al., 2012;Brunner and Bechinger, 2002;Mangold et al., 2003), quasiperiodic (Mikhael et al., 2008), or random (Deutschländer et al., 2013). For colloids interacting with periodic 1D surfaces, various effects such as transitions among liquid, 2D hexagonal solid, and smectic states appear as a function of increasing substrate strength (Bechinger et al., 2001;Tierno, 2012;Tierno et al., 2008). The next level of complexity is to consider colloids interacting with two dimensional periodic arrays, such as an egg carton (Agra et al., 2004;Bohlein et al., 2012;Mangold et al., 2003;Reichhardt and Olson, 2002;Šarlah et al., 2005), muffin tin (Bechinger et al., 2001), or more complex potentials (Demirrs et al., 2013;Gunnarsson et al., 2005;de las Heras et al., 2016;Loehr et al., 2016a;Massana-Cid et al., 2019;Tierno and Fischer, 2014;Yellen et al., 2005). Such systems can mimic the ordering of atoms on 2D surfaces (Coppersmith et al., 1982), vortices in type-II superconductors with nanostructured pinning (Baert et al., 1995;Harada et al., 1996;Martín et al., 1999), and vortices in Bose-Einstein condensates (Tung et al., 2006) interacting with 2D optical trap arrays. In this case, commensuration effects arise when the number of colloids is an integer multiple of the number of potential minima, giving an integer filling factor f , where at f = 1 each trap captures a single colloidal particle. Experiments (Bechinger et al., 2001;Bohlein et al., 2012) and theoretical studies (Agra et al., 2004;Reichhardt and Olson, 2002;Šarlah et al., 2005) of colloids interacting with 2D substrates reveal a variety of novel orderings, including commensurate colloidal molecular crystals for f = 2, 3...N. When the system is away from commensuration, such as just above f = 1.0, the additional particles act like highly mobile kinks, and both simulations and experiments for colloids on 2D arrays have revealed the motion of these kinks under an applied drive (Bohlein et al., 2012;McDermott et al., 2013;Vanossi et al., 2012). Other studies of colloids have focused on the Aubry transitions that occur as a function of substrate strength (Brazda et al., 2018).
Since various substrates for colloidal particles can be created readily, including systems in which a single trap contains a double well potential (Babič and Bechinger, 2005), a natural question to ask is whether it is possible to create a colloidal system containing effective spin degrees of freedom. For colloids on regular 2D arrays, it was shown that for integer f = 2 and higher, the colloids in each trap act like dimers with a spin degree of freedom that can be mapped to an Ising one, making it possible to use colloids in periodic substrate arrays to model various spin systems (Agra et al., 2004;Reichhardt and Olson, 2002;Šarlah et al., 2005).

Artificial colloidal ice: simulations
Motivated by these ideas, and influenced by the work on magnetic ASI, Libál et al. proposed that 2D array of double well traps for colloids could mimic a square ASI (Libál et al., 2006). In this case the basic unit cell consists of four double well traps arranged as shown in Fig. 2(a), where f = 1 so that each trap captures a single colloidal particle. The colloid-colloid interaction potential is of Yukawa or screened Coulomb form, V (r) ∝ q e exp(−κr)/r, where q e is the electrostatic charge and 1/κ is the screening length.
The simulations for N charged colloids on a 2D array of N traps are performed using 2D Brownian dynamics (BD). The colloid dynamics are overdamped and the where in rescaled units the damping constant is set to η = 1.0 and a 0 is used as the unit of distance in the simulation. The colloid-colloid interaction force is given by, , is the solvent dielectric constant, and κ = 4/a 0 which is the typical screening length for charged colloidal systems. These simulations neglect hydrodynamic interactions between the colloids, which is a reasonable assumption for charged particles in the low volume fraction limit or where the dynamics are dominated by hopping events rather than large scale flows. The effects of thermal noise are captured in the simulation via the force term F T , which represents random Langevin kicks with the properties F T i = 0 and For charged colloidal systems, a bias can be applied using an electric field in order to mimic the effect of an external magnetic field in magnetic spin ice. In the simulation model, this biasing field is represented by the term F ext i . If, for example, the field F ext i = F dc (x +ŷ) is used to bias the colloids along a 45 degree angle from thex axis, the ice state illustrated in Fig. 2 In the absence of a substrate, charged colloids confined to two dimensions will form a triangular lattice (Kusner et al., 1994). For the traps in Fig. 2(a), there are two equally favorable locations. The two resulting configurations can be mapped to a spin degree of freedom in which the spin is defined to point toward the end of the trap occupied by the colloid. The lowest energy state for the single unit cell shown in Fig. 2(a) has all the colloids sitting at the end of the trap that is furthest from the vertex, so that all of the effective spins point away from the vertex. This minimizes the colloid-colloid interaction energy. The highest energy state has all of the colloids sitting close to the vertex, so that all of the effective spins point toward the vertex. This corresponds to a doubly charged monopole state.
As shown in Fig. 2(b), when many unit cells are assembled into a lattice, it is not possible for all vertices to adopt their lowest energy configuration simultaneously since the geometric arrangement competes with the symmetry of the pair interaction. If the electrostatic screening is very strong or the trapping sites are sufficiently far enough apart, the colloids do not interact with each other and the arrangement of the effective spins is random, as illustrated in Fig. 2(b) where the populations of the different possible vertex states matches with what is expected for a purely thermal distribution. For strong colloid-colloid interactions, the system forms the lowest energy collective ground state, in which each vertex has two colloids close to it and two colloids far from it, equivalent to the ice-rule of the square ASI (2-in/2-out), Fig. 2(c). If a biasing field is applied at 45 • along thex axis, the spin ice rule is still obeyed but the system forms a biased state with high energy vertices which still obey the ice rule, Fig. 2(d).
The system can be characterized according to the six possible vertex types, whose energy is listed in Table I. From the lowest to the highest energy, the vertex types are: type I, with all four spins out to create a double monopole; type II, with one spin in and three spins out to create a monopole; type III, the spin ice rule obeying state; type IV, the biased ice rule state; type V, with three spins in and one spin out to create a monopole; and type VI, with all four spins in to create a high energy double monopole. An immediate difference from the magnetic version of square ASI is that in the colloidal artificial ice, the monopole states II and V, and the double monopole states I and VI have different energies. This becomes important in the mobility of the monopoles, and is an indication of the fact that the particle based spin ice system minimizes the global interaction energy rather than the local vertex energy.
Using numerical simulations of this geometry, Libál et al. (2006) examined the evolution of the vertex populations when passing from the weak to the strong interacting limit. Figure 3(a) shows the fraction N i /N of the N vertices that are of type i as a function of colloidal charge q e in a system with fixed κ and distance between the traps. At q e = 0, the fractions of the vertex types match what is expected in a non interacting system, while as q e increases, type I and VI vertices disappear first followed by type II and IV vertices, until for q e ≥ 2.0, only the type III ice rule obeying vertices remain. Figure 3(b) illustrates the fraction of vertex types as a function of the distance d between the traps in a system with fixed q e and κ. When d is small, the colloids are strongly coupled and type III vertices dominate, while as d increases, the vertex distribution gradually shifts back to the random configuration. The results in Fig. 3(b) are very similar to the observations in the initial work on magnetic ASI, where increasing the spacing between magnetic islands produced a more random state (Wang et al., 2006). In a system with fixed q e and d and changing screening strength κ, shown in Fig. 3(c), the interaction between the colloids is weak for strong screening and the vertex distribution becomes random. A sample with fixed d and κ at q e = 0.4 that is subjected to an additional biasing drive F dc at 45 degrees from the x axis is initially in the type IV biased state. As q e increases, Fig. 3(d) shows that a transition occurs into the lower energy type III square ice state. Libál et al. (2006) also found that a transition from an ordered ice state to a disordered state occurs as a function of increasing temperature. These results indicate that a particle based artificial ice system exhibit the same ice rule obeying states as nanomagnetic ASIs, with the advantages of the mesoscopic character of the particles which makes the microscopic degrees of freedom readily accessible.

Experimental realization
The realization of a full optical system, following the original proposition (Libál et al., 2006), remains a challenging task because (a) it would require a large optical power to generate enough double wells and (b) it is difficult to finely tune DLVO based electrostatic interactions between colloidal particles. An alternative realization was demonstrated recently using a combination of different techniques including soft-lithography, magnetic manipulation, and optical tweezers (Ortiz-Ambriz and Tierno, 2016). The experimental system, illustrated in Figs. 4(a,b), features interacting paramagnetic colloids confined by gravity in lithographically generated topographic double wells. Each trap contains two deep wells connected by a small central hill, and these indentations are arranged in a regular lattice, such as the square ice geometry illustrated in Fig. 4(a). Using optical tweezers, each trap is filled with one paramagnetic colloid consisting of a spherical polymer particle that is responsive to magnetic fields. Repulsive and tunable interactions are induced by applying an external magnetic field B = Bẑ, perpendicular to the particle plane. The applied field induces a dipole moment m = πd 3 χB/(6µ 0 ) within the particles, where d is the particle diameter, χ is the magnetic volume susceptibility, and µ 0 is the permeability of the medium (water). All particles interact through magnetic dipolar interactions, Pairs of particles (i, j) with moments m i,j = me i,j and at distance r = |r i − r j | interact through dipolar forces, with an interaction potential given by, . This potential is maximally attractive (repulsive) for particles with magnetic moments parallel (perpendicular) to r. In an unconstrained system, when the applied field is perpendicular to the particle plane, U d reduces to an isotropic repulsion between parallel particles in the same plane, U d (r) = ω/r 3 . The field amplitude is chosen such that a confined particle, when subjected to the dipolar force from a neighbor, can cross the central hill with gravitational potential U g , but never escape from the bistable confinement (U c ), U g < U d < U c , as shown in the small inset in Fig. 4(a). Following the original idea (Libál et al., 2006), one can assign a vector to each particle pointing towards the well occupied by the colloid, Fig. 4(b). This mapping makes it possible to construct a set of vertex rules similar to ASIs. However, in contrast to ASI whose islands have in-plane dipole moments, the colloidal ice features out-of-plane dipoles, and thus the energetic hierarchy is similar to that found for repulsive electrostatic colloids. Within a homogeneous lattice, collective interactions between the particles oppose the local energetics and enforce the ice rule for the colloidal system . Indeed, systematic measurements combined with Brownian dynamics simulations confirm that, for high field amplitudes, the colloidal ice follows the ice rule (Ortiz-Ambriz and Tierno, 2016). The experiments were performed by setting the system in a random configuration using optical tweezers, turning on the magnetic interactions, and waiting for the system to reach thermal equilibrium. The initial disordering process was necessary since, in contrast to the original simulations (Libál et al., 2006), the magnetic colloids experience negligible thermal fluctuations due to the relatively large particle size, d ∼ 10µm. The disordering process can be considered as a way to couple the system with a virtual heat bath, while increasing the strength of the pair interaction represents a cooling procedure. Moreover, in analogy with ASIs, the particlebased ice reaches a unique ground state (GS) filled by type III vertices. The loss of degeneracy results from the different particle distance at the vertex, which makes the type III and type IV vertices energetically different. The most energetically unfavorable vertices with three (type V) or four (type VI) colloids in become topologically connected to the low energy type II and type I vertices, reducing their occurrence at large B, as shown in Fig. 4(c).

Defect dynamics, grain boundaries, and logic gate
The experimental system provides a versatile approach for probing the effect of disorder and topological defects in frustrated lattices. The optical tweezers, which are independent from the collective magnetic coupling, can be used to manually add or remove particles from the double wells. One method for creating an artificial de- The square plaquettes can have a clockwise (orange arrows) or counterclockwise (yellow arrows) chirality (type IV vertices), or be achiral (green crosses). Scale bars are 20µm for both cases. Image edited with permission from Ortiz-Ambriz and Tierno (2016).
fect line is to start with the square ice GS (all type III vertices) and flip particles in a sequence that produces a double line of type IV vertices connecting two high energy type II and V vertices. These defects can be described within the framework of the "dumbbell" model (Castelnovo et al., 2008). To each vertex one can associate an effective topological "charge" q that quantifies the degree of violation of the ice rule. As explained in Section II.1, a spin pointing towards (away from) the vertex center will have a positive (negative) charge q i , and the total charge at each vertex i will be q = i q i . Thus, the GS of the square lattice corresponds to q = 0 (type III), whereas type II (V) vertices have a positive q = +2 (negative q = −2) charge, Fig. 5(a). Using this formalism, it was shown that in a 3D spin ice and at low temperature, topological defects interact only via a magnetic Coulomb law, V (r) ∝ Q 2 /r, where Q is the topological charge and r is the separation distance between the two monopoles. However, numerical simulations demonstrated that for a 2D square ASI, such strings shrink due to an additional line tension term (Mól et al., 2009) that results from the missing degeneracy of this lattice at the single vertex level. In the 2D case, the interaction potential between two defects becomes V (l) = −Q/l + υl + c, where υ is the line tension and c is a constant associated with the creation of defect pairs (Silva et al., 2013). Direct measurements of the average defect line length l showed that the defects shrink following overdamped dynamics described by where γ is a friction coefficient. Eq. 4 lacks the inertial term that is present when describing the motion of topological defects in ASIs (Vedmedenko, 2016). This first order nonlinear equation admits as solution the implicit function, where α = Q/υ is the ratio between the Coulombic and line tension contributions. The solution describes a monotonous shrinkage of the defect line starting from l 0 at t 0 , and was used to fit both the experiments and the simulation data. Such analysis makes it possible to quantify α = 0.0290 ± 0.0014a 2 , where a = 29µm is the lattice constant of the square system. This Coulombic charge is one order of magnitude lower than the value calculated for magnetic ASI (Mól et al., 2009), indicating that the line tension contribution is enhanced in the colloidal ice. Nevertheless, two topological charges that are bound by a defect line in either ASIs or in the colloidal ice interact with a similar Coulombic law.
The modeling of the magnetic colloids is similar to that used for charged colloidal systems since the colloids are confined in double well traps and experience a repulsive pairwise interaction with each other. The colloid dynamics is governed by an overdamped equation of motion given by 1 µ where D is the diffusion constant, µ is the mobility, dt the time step, and N [0, 1] is a Gaussian distributed random number with a mean of zero and a variance of 1. The first term is the thermal force. Magnetization of the colloids in the z direction produces a repulsive particle-particle interaction force F pp (r) = A cr /r 4 with A c = 3 × 10 6 χ 2 V 2 B 2 /(πµ 0 ) for colloids at a distance r apart, where V = πd 3 /6 is the colloid volume. The substrate force F i s traps the colloids in the double wells. Simulations of magnetic colloids have focused on a range of parameters that matches what is used in the experiments.
Since it is possible to tune the interaction forces between the particles by changing the applied field, it is interesting to model the interaction between different types of monopoles in a square ice for increasing field amplitude. In the work of Loehr et al. (2016b) the magnetic colloids in a square ice substrate were initialized in a GS containing two monopole excitations of opposite charge and separated by a distance d. These monopoles are attracted and approach each other at a velocity that increases with the applied magnetic field. At high fields, the monopoles annihilate more rapidly through a nucleation process in which additional monopoles appear and break the defect string which extends between the original two monopoles. This work also showed that different monopole species move at different velocities, in contrast to what would be expected for magnetic spin ice systems. By applying a periodic biasing field, the asymmetry in the monopole mobilities can be exploited in order to create a ratcheting motion of a defect line (Libál et al., 2017).
It is also interesting to investigate grain boundaries (GBs) in the colloidal ice. In general, GBs are ubiquitous in condensed matter and they influence the performance of a variety of systems including high T c superconductors (Graser et al., 2010), organic films (Rivnay et al., 2009), and direct-band gap semiconductors (van der Zande et al., 2013). In the square colloidal ice, as in ASIs, GBs emerge as defect lines that separate unmatched regions of GS, Fig. 5(b). While GBs have been observed in magnetic samples (Morgan et al., 2011;Zhang et al., 2013), monitoring their formation and dynamics remains a challenging task due to the small length scales of the ASIs. In terms of plaquettes and not vertices, the GS of the square ice is composed of a checkboard pattern of loops with alternating chirality, indicated by yellow and orange in the left panel of Fig. 5(b).
For the square ice this GS is twofold degenerate, and a 90 degree rotation around a vertex produces another GS. A domain wall, shown in the right panel of Fig. 5(b), can be constructed by a line of achiral cells that separate two incompatible regions of GS. These are very stable defects because, in contrast to a double line of type IV vertices, it is necessary to modify all the vertices on one side of the domain wall to make them compatible with the GS of the other side. This feature could be used to store binary information. By assigning a value of 0 (1) to a clockwise (counterclockwise) chiral cell, it is possible to use a binary representation to write information in the GS of a lattice of colloidal ice or nanoscale ASI, and later erase the information using a strong field which restores the full GS (Ortiz-Ambriz and Tierno, 2016).
A major interest in nanoscale magnetism and ASIs is the realization of logic circuits based on a dense array of strongly magnetized elements. The particle-based ice could provide a guideline for the realization of simple and resettable logic ports based on magnetic dipolar interactions. Figure 6 shows a proposal for a "NOR" gate, a functionally complete port capable of generating all logical functions (Bartee, 1991). The gate is based on a metastable biased state, filled by high energy type IV vertices, and it can be initialized and reset by an external force F ∼ (B · ∇B) applied along one diagonal, F 1 = F 1 (ŷ −x). In this system, the colloids in two horizontal rows of traps have been replaced with particles that have a higher magnetic susceptibility (pink arrows in Fig. 6). A single particle in the upper left corner is used to trigger the propagation of a defect line. Then a second biasing force F 2 is applied to flip only the two chains containing high susceptibility particles. The input of the gate is applied by switching the line of colloids with higher susceptibility, such that a flipped line means true (1) and an unflipped line means false (0). As shown in Fig. 6(c-d), the defect starts to propagate under a perpendicular field B = Bẑ and is deflected by the flipped input lines. After propagation, the location of the charge in the bottom row gives the output. If it is located in the lower left corner, the output is (1). If either of the two input lines are in a 1 state, the defect line will not be deflected, and the topological charge will end up in a different position, giving an output of (0).

Effect of disorder, doping, and system memory.
A general open question in frustrated systems is how robust or fragile the frustrated states are in the presence of quenched disorder. Numerical simulations of spatially extended samples with significant statistics can shed light on this aspect. In the colloidal artificial ice, one method of introducing quenched disorder is by randomizing the energy of the barrier at the center of each double well. In a square ice geometry subjected to a biasing field, such randomness induces the formation of +1 or −1 monopoles in the background of biased ice rule obeying vertices. Libál et al. (2012) studied the effect of quenched disorder on a colloidal square artificial ice and found that under simulated annealing, instead of a completely ordered state, the system formed a partially disordered configuration containing monopoles and biased ice rule obeying vertices.
Under the application of a cyclic biasing drive, the monopole defects begin to annihilate until the sample reaches a reversible state in which the same effective spin configuration appears during each biasing cycle, so that the system exhibits what is known as return point memory (Pierce et al., 2005). Libál et al. (2012) also showed that the return point memory state is reached much faster in a kagomé colloidal ice than in a square colloidal ice. In the kagomé ice, irreversibility is produced by the motion of individual defects, which can be pinned rapidly by the quenched disorder, whereas in the square ice, the irreversible motion arises from the reconfiguration of grain boundaries, and these extended objects require a longer time to find a pinned configuration.
Return point memory can be quantified by examining the overlap function q o which compares the effective spin configurations from one cycle to the next at the same biasing field, Here S i is an effective spin equal to S i = 1 (S i = −1) if the colloid is sitting in the right or top end of the trap (left or bottom end), and n is the number of cycles that the external force F ext is applied. When q o = 1.0, the spin configuration is perfectly reversible and is exactly the same during each biasing field cycle, while q o = −1.0 would indicate that the spin configuration is exactly the opposite of its previous value. Fig. 7(a,b) shows q o versus F ext for the square and the kagomé colloidal artificial ices containing random disorder in the barrier heights implemented using a Gaussian distribution of width σ. The kagomé system is characterized by a modified version of the ice rule with 2-in/1-out ot 1-in/2-out vertex types. With repeated cycles of the biasing field, q o approaches 1.0, indicating the emergence of return point memory. This process takes about n = 10 cycles in the square ice but only n = 4 cycles in the kagomé ice. The steady state arrangement of artificial spins is not ordered but contains numerous defects; however, the same pattern of defects appears after each biasing field cycle. Return point memory effects have also been observed in ASI by Gilbert et al. (2015). It would be interesting to understand how these effects appear in many of the other proposed ASI geometries (Morrison et al., 2013). Further modeling work on kagomé and square colloidal artificial ices that took into account the effects of temperature showed that when biasing fields are present or when the system is prepared in different metastable states, a series of structural transitions can occur, and that states with monopole ordering appear (Olson . In ASI, each spin is constrained to point along the axis of the magnetic nanoisland, although is possible to remove entire islands. In particle based artificial ices, it is possible to introduce doubly occupied sites, which would be like having spins pointing in both directions on the same nanoisland. Conversely, it is also possible to remove a particle to create a vacancy or dilution of the spin arrangement. Libál et al. (2015) considered the square and kagomé geometries containing doubly occupied wells where the system was initiated in a GS configuration and then subjected to an increasing temperature. The effects of the doping are very different in the two ice geometries. Figure 8 illustrates the effects of doping on square ice (left panels) and kagomé ice (right panels) when the number of dopants is held constant but the temperature is increased. For the square ice at low temperature, N3 (three-in/one-out) vertices form around each doping site in order to screen the doubly occupied wells. As T increases, N1 defects and N2 biased defects begin to appear near the doping sites. Thus in the square ice, the doping sites act like weak spots or nucleation points that generate increased hopping of the colloids in the surrounding sites. In contrast, the kagomé ice GS can readily absorb the doubly occupied sites without creating monopole states, as shown in Fig. 8(d) at low temperature. This is because the kagomé GS consists of an equal number of N1 gs and N2 gs vertices, and the system can shift this balance so that there are slightly more N2 gs vertices, enabling it to incorporate the extra charge of the doubly occupied sites while still remaining in a GS configuration. As the temperature increases, the doped sites in the kagomé lattice reduce rather than increase the amount of hopping occurring in the neighboring vertices.

Kagomé ice and its inner phase
Kagomé systems have been the subject of intense research in nanomagnetism because this simple geometry provides a truly degenerate frustrated system (Qi et al., 2008;Tanaka et al., 2006). The degeneracy arises at the single vertex level, where the three spins have equal distances and interaction energies. The ice rule of this lattice have associated a net topological charge, with q = +1 for 2-in/1-out and q = −1 for 1-in/2-out. Although the kagomé ice does not have a well defined GS, it was shown theoretically in ASIs that the long range tail of magnetic dipolar interactions can further reduce the entropy and give rise to ordered phases. Specifically, multipole expansion (Möller and Moessner, 2009) and numerical simulations (Chern et al., 2011) predicted that lowering the temperature should favor the transition from a short range ordered phase ("ice I") to a long range ordered state ("ice II"). In the latter case the spins on the hexagons define chiral and achiral loops along the lattice. Finally, the system can transition into a completely ordered "spin solid", or chiral phase, with chiral and achiral loops alternating in a chessboard-like pattern. Direct visualization of these exotic states with their associated relaxation dynamics remains elusive. An indirect signature of the ice I phase has been reported based on magnetotransport measurements of the Hall signal in a cooled ASI sample (Branford et al., 2012). Further, the ice II phase was visualised via MFM on a lattice of permal- loy nanoislands heated above the Curie temperature of the constituent material (Zhang et al., 2013). More recently, low energy muon spectroscopy was used to probe the existence of peaks in the muon relaxation rate that can be identified with a critical temperature associated with a phase transition that bridges such phases (Anghinolfi et al., 2015). Numerical simulations of the magnetic colloidal system in the kagomé ice geometry show the emergence of different states as the colloid-colloid interaction strength changes . For weak interactions, the system is in a paramagnetic state, as shown in Fig. 9(a). As the interaction strength increases, the system first enters a state where the ice rule is obeyed in an ensemble of disordered spins, as illustrated in Fig. 9(b). This is followed by the appearance of the topologically ordered charged state shown in Fig. 9(c), and finally, in the limit of very strong interactions, a three fold degenerate ferromagnetic state emerges as shown in Fig. 9(d). Unlike the phases predicted and experimen-tally observed with ASI, the colloidal kagomé ice forms a ferromagnetic state at high interaction strength. This ferromagnetic state arises due to interactions between non-nearest-neighbors, rather than simple vertex energetics, and disappears when the interactions are short ranged, explaining why it was not observed previously.

Other geometries
In addition to the square and kagomé artificial ices, other types of frustrated systems have been realized with microscale colloids. Chern et al. (2013b) proposed a system of colloids interacting with a honeycomb array of optical traps that contain three wells instead of two, Fig. 10(a). This system is a realization of a fully packed loop (FPL) model and Baxter's three-coloring problem (Baxter, 1970). As a function of temperature and interaction strength, a series of phases appear in the triple well system, including a stripe state, stripes with sliding symmetries, random packed loop states, and disordered states containing broken loops. Figure 10 FPL models have been used to explain a broad class of phenomena in magnetism, optics, and polymer physics (Blöte and Nienhuis, 1994;Duplantier, 1998;Jaubert et al., 2011;O'Holleran et al., 2008), but physical realizations are rather scarce. An experimental realization of the colloidal version has been achieved using a honeycomb lattice of triangular shaped magnetic minima, as shown in Fig. 11(a). The complex magnetic potential arises when a uniaxial ferrite garnet film (FGF) is subjected to a constant magnetic field B z = Bẑ. The FGF is composed of a triangular lattice of magnetic "bubbles," i.e., cylindrical ferromagnetic domains that are uniformly magnetized (Tierno et al., 2007(Tierno et al., , 2009. Paramagnetic colloids dispersed above this honeycomb magnetic lattice self-assemble into interacting microscopic dimers. To anneal the lattice into its minimum energy state, the system is subjected to a precessing magnetic field consisting of a combination of the perpendicular field and an in-plane rotating field B xy ≡ B xy [cos (ωt)x − sin (ωt)ŷ]. The resulting field B = B z + B xy performs a conical precession around theẑ axis with angular frequency ω, and sets the dimers into rotational motion. Depending on the field parameters, the resulting dimer arrangement can be mapped to a long range striped phase or to a random FPL state, as illustrated in Fig. 11(b). The mapping of the dimers to the FPL model is achieved by using the arrow representation originally introduced by Elser and Zeng for the spin-1 2 kagomé antiferromagnet (Elser and Zeng, 1993). Since the dimer sits on one of the three sides of a triangular minimum, an arrow can be defined that points from the dimer center to the free corner of the triangle. This is shown by a green dot with a blue line in  Another recent realization of a particle ice composed by topographic double wells but in a different geometry is the highly coordinated triangular lattice (z = 6) presented by Lee and Tierno (2018). Although not directly inspired by the water or spin ice materials, the triangular geometry reflects the spin disposition in several magnetic materials with moments lying on weakly coupled parallel planes (Dublenych, 2017;Nakatsuji et al., 2005). The triangular order reflects the natural arrangement of repulsive particles in the absence of the substrate. However, such ordering can be frustrated by the presence of the central hill in the double wells. In this geometry, collective interactions between the particles lead to a unique GS characterized by vertices with three colloids pointing inward and three outwards, similar to what was predicted for ASIs by Mól et al. (2012) and Rodrigues et al. (2013). It was also found that the use of a bias force that magnetizes the system allows the GS to be accessed easily via a structural, martensitic-like transition characterized by the coherent sliding of one particle at each vertex. FIG. 11 (a) Magnetostatic energy landscape of a ferrite garnet film calculated under a static, perpendicular magnetic field B = 3.6mT. The inset shows a 3D view of one triangular minimum. (b) Random FPL phase observed with a monolayer of paramagnetic colloids above the FGF after an annealing performed with a precessing magnetic field, from .

Nature of the ice rule in colloidal ice
While the ASI and the particle-based ice have often been considered as equivalent, their frustration and energetics are fundamentally different. At the nearest neighbor level, a vertex in the particle-based ice has a lowest energy unfrustrated configuration in which all particles are far from the vertex. Indeed the nearest neighbor energy of a magnetic spin ice vertex with n spins pointing toward the vertex is proportional to the square of its topological charge, thus favoring the ice rule, which minimizes the topological charge 1 . In contrast, the energy of a vertex in particle-based ice scales as where n(n − 1) is the number of repulsive interactions among n particles. This energy favors n = 0 and n = 1 states, which correspond to large negative charges according to Eq. (1). In particular, the energies of the particle-based ice vertices do not posses a Z 2 symmetry. Thus the local energetics actually work against the ice rule. The ice rule in particle-based ice is only recovered in the thermodynamic limit and its origin is collective: not all vertices can simultaneously have all particles located away from the vertex. In a finite size particle-based ice system of uniform coordination, the energy can be reduced by pushing particles onto the boundaries. The total charge in the bulk is = 1 2 + 1 2 FIG. 12 Schematics showing that a colloidal ice, here in a random configuration, can be decomposed into a spin ice and a background exerting a geometric field, from Nisoli (2018b) then bounded by the flow of the pseudo-spin field through the boundary. Since the pseudo-spins have an absolute value of one, the density of charge in the bulk scales at least as the reciprocal of the length of the boundaries, leading to an ice manifold in the thermodynamic limit.
In a system of multiple coordination, however, such as z = 4, 3, we can consider the z = 3 vertices as an "internal boundary" onto which the z = 4 sub-lattice can push topological charges. This would lead to the emergence of negative charges on z = 4 vertices, in violation of the ice rule. Note that this behavior, which has been predicted (Nisoli, , 2018b and verified numerically and experimentally  points to an essential difference in the origin of the ice rule in magnetic and particle-based spin ices. The ice-rule in magnetic ice is known to be robust against decimation (Morrison et al., 2013), mixed coordination (Gilbert et al., 2014(Gilbert et al., , 2016a, and dislocations (Drisko et al., 2017), and is present even in finite sized clusters . In particle-based ice, however, the local energy in Eq.(9) opposes the ice rule, which is regained only as a collective compromise. Other differences appear in the kinetics. For example, when defect lines in square colloidal ice are driven with a field, the two monopoles at the ends of each line have different mobilities since, unlike the monopoles in magnetic spin ice, they have different energies (Libál et al., 2017). The similarities and differences of the magnetic and particle-based ices can be quantified exactly. A mean field approach provides useful quantitative predictions that are fit well by experimental and numerical results (Libál et al., 2017;. Starting with the approximation of a gas of decorrelated vertices, a constraint must be introduced so that the total topological charge is conserved. This takes the form of a Lagrange multiplier φ that modifies the original vertex energies from those in Eq. (9) tõ For a lattice of coordination z, the choiceφ ∼ (z − 1) ensures conservation of topological charge, giving a spinice-like effective energetics,Ẽ n ∼ q 2 n , and therefore a spin-ice behavior. Eq. (10) can be interpreted as follows: the collective effect of the particle-sharing vertices can = FIG. 13 Schematics showing that a particle-based ice obtained by removing traps from a regular lattice is equivalent to a spin ice stuffed with negative charges at the locations of the missing traps. This leads to violations of the ice rule through formation of positive topological charges (blue circles), from Nisoli (2018b) be subsumed into a emergent field φ, which modifies the energetics of the individual vertex. Indeed, in a better approximation, one can allow the constraint-enforcing field φ(x) =φ + η(x) to fluctuate in space, since it mediates an entropic interaction among the topological charges to which it is coupled. This results in a familiar Debye picture for purely entropic screening with a correlation length ξ 2 ∼ T /Q 2 , where Q 2 is the charge fluctuation of the manifold. For instance, in a fully degenerate square ice, Q 2 = 0 in the ice manifold and thus the correlation length is infinite, as mentioned in the section on the topological properties of the ice rule. In contrast, Q 2 ≥ 1 in kagomé ice since each vertex has a charge of at least ±1 and the ice phase is never critical.
Consider now a lattice of multiple coordination. It is impossible to find a value of φ that can return a Z 2 invariant, ice-like effective energetics in Eq. (10) for both sublattices. One of the lattices must break the ice rule. How this happens can be understood by a geometrically intuitive, exact treatment (Nisoli, 2018b). The energy of the system is given by where ψ(r) is an isotropic repulsive interaction and y labels the position of the particles in the traps. These positions represent a binary variable that we can represent as -or -. Indeed Eq. 11 does not look like a spin ice Hamiltonian. Then we can ascribe a positive charge to the real particles and we can consider the empty locations y − of the traps as virtual negative charges , which repel (attract) other negative (positive) charges. We can fractionalize a trap on an edge x as i.e. a positive dumbbell -(a trap doubly occupied by positive charges) plus a dipole of negative and positive charges represented by a spin σ = -located in x, the center of the trap. Then the energy in Eq. (11) can be rewritten as The first term is the interaction between dumbells and is clearly a spin-ice Hamiltonian. J ii (x) is a tensor field and the background field B mediates the interaction between dipoles and the positive dumbbells, both of which can be reconstructed from the particle-particle repulsive interaction ψ. Fig. 12 illustrates this decomposition. It follows that a particle-based ice becomes equivalent to a magnetic spin ice when the second term in Eq. (13) is zero. That is certainly true if a lattice has point reflection symmetry in the middle points {x} of each edge, and explains why the kagomé and square particle-based ices follow the ice rule, as found previously numerically and experimentally (Libál et al., 2006;Loehr et al., 2016b;Ortiz-Ambriz and Tierno, 2016).

Decimated systems
The analogy between the colloidal ice and nanoscale ASI does not hold in general for more complex geometries. Imagine decimating a simple lattice such as the kagomé lattice by removing some traps. This breaks the reflection symmetry and the background field will perturb the spin ice energetics since the saturated dumbbells are replaced by negatively saturated ones at the locations of the missing traps. As shown in Fig. (13), the decimated system is equivalent to a spin ice stuffed with negative charges. The charges polarize the dumbbells close to the decimated vertices, breaking the ice rule. In these geometries with mixed coordination z, the conservation of the topological charge holds only at the global level, not at the sub-lattice level.
In the case shown in Fig. 14(a), the square system (z = 4) is decimated with optical tweezers to a lattice of mixed coordination, z = 3, 4. The decimation process is equivalent to a partial, random dimer covering of the edges, as illustrated in Fig. 14(b). Randomly chosen edges of the square lattice are covered by dimers under the constraint that each vertex contains at most one dimer. Removing an edge between two "dimerized" vertices with z = 4 gives rise to only vertices of z = 3. The process avoids formation of z = 2 vertices, and thus simplifies the vertex hierarchy of the mixed coordination system (Gilbert et al., 2014), although similar effects can also be observed in the presence of z = 2 vertices.
Combined experiments and numerical simulations show that in such geometry, the ice rules are spontaneously yet selectively violated via the formation of negative topological monopoles of charge q = −2 on the vertices of high coordination z = 4. The low coordinated vertices z = 3 still obey the ice rule (2-in/1-out or 2-out/1-in); however, the relative ratio of q = 1 to q = −1 charges changes in order to compensate the negative charges that accumulate around the z = 4 vertices. Even in this situation, the total topological charge of a system of "dipoles" remains zero.
A quantitative analysis of the experimental and theoretical results provides additional insight into the decimation process. Figure 14(c) shows the relative fraction of z = 4 vertices, n z4,q , versus the ratio between the two vertex coordinations η = N z3 /N z4 , grouped by topological charge. The number of negative q = −2 charges generated around the z = 4 vertices increases with increasing η, quantifying the strength of the ice rule violation. Above a critical decimation threshold, the entire system disorders due to the spontaneous appearance of entropy-driven negative monopoles which induce topological charge transfer between the sublattices. As a result, the colloidal ice has a "fragile" low-energy manifold that is produced by an energetic compromise between locally excited vertices. This is in contrast to magnetic ASI, which are structurally "robust" ices. These observations prompt different exciting ideas. Since ice rule fragility is associated with topological charge transfer among sub-lattices, these new phenomena can be exploited for domain wall engineering, in which membranes that are semipermeable to the topological charge of de- fects are structurally designed. These results also apply to lattices of different coordination or to ASI systems of nonzero residual entropy.

Discrete versus continuum models
The interactions between the pseudospins in the colloidal ice are given by a generalized dipolar interaction constructed from the effective interaction V (r) between pairs of colloids (which in the case of magnetic colloids is V ∼ r −3 ). This implies that colloids with different repulsive interactions exhibit different phases within their ice manifold, and in particular, that electrically charged colloids are expected to more faithfully reproduce the phases of magnetic spin ices. Inspired by these results, in a theoretical study, Le Cunuder et al. investigated the energetics and phase transitions that occur in the kagomé geometry by varying the temperature rather than increasing the interaction strength. Numerical calculations of the total magnetostatic energy for a system containing N particles showed that the chiral state ("spin solid") always has a lower energy than the ferromagnetic one, a result which is robust regardless of the system size, Fig. 14(c). This calculation was performed by assuming that the particles are fixed at the bottom of each trap, without considering any relative displacement.
The apparent discrepancy between the analytic re-sults and those obtained in the simulations and experiments described above was resolved by performing two types of Monte Carlo simulations (Le Cunuder et al., 2019). In the first, discrete model, particles are only allowed to jump from one side of the bistable traps to the other, while in the second, continuous model, particles are allowed to make these same jumps and also to move continuously around the lowest point of the bistable traps. While the discrete simulations reproduced the well-known phase structure of the dipolar spin ice, the continuous Monte Carlo simulations showed a single phase transition from charge disorder directly to the ferromagnetic state. Furthermore, using the continuous model, it was shown that by modifying the strength of the traps, the GS of the particle-based ice can be changed from the chiral ordering to the ferromagnetic one.

Future directions in particle ice
The experimental realization of particle ice represents a starting point for exploring the novel physics that emerges from the mesoscopic character of the particles. A variety of avenues for future research are now available.
Relaxation and dynamics: Apart from the annihilation of simple defects, the kinetic mechanisms behind the relaxation dynamics of the particle ice remain largely unexplored. Since direct visualization of the dynamics is possible in colloidal ice, correlation functions for relatively long times can be extracted from the particle positions, and rearrangement events (spin flips) can be identified as the system ages from a metastable state. An interesting question is whether the colloidal ice could exhibit glassy behavior and kinetic arrest. The presence of disorder in the system, such as an inhomogeneous distribution of the hills within the double wells, or induced decimation, may make the energy landscape more rugged and induce glassy behavior. There have been reports of evidence of glassiness in frustrated systems from both theoretical models and experiments. For example, a spin glass phase was predicted to occur for a decimated square ice in an ASI (Sen and Moessner, 2015), and experiments with buckled colloidal monolayers  showed evidence of dynamical arrest.

Thermalization effects:
Using smaller particles with large thermal fluctuations would allow spontaneous particle switching within the double wells, a feature that is currently absent in the experimental system. This platform could be used to explore thermalisation effects in the colloidal ice. Such an investigation could be performed by using optical tweezers to prepare a square or kagomé ice lattice in the lowest energy state while keeping the particle-particle interaction strength high. Subsequently lowering the magnetic field would then favor thermal disordering starting from the GS and make it possible to determine the lifetime, fraction, and dynamics of the emerging, thermally induced "Dirac" strings.
Complex geometries: The lithographic approach imposes no limit on the type of two-dimensional structures that can be engineered. While most of the works on particle based artificial ices have focused on square or kagomé lattices, numerous other geometries that have been proposed for ASI (Morrison et al., 2013) could also be realized in particle based systems. Since the particle ice system minimizes the global energy rather than the local vertex energy, it often exhibits fragility , so it is likely that these systems could have very different behaviors than their ASI counterparts. In addition, aperiodic structures or even disordered lattices with hyperuniform properties can easily be designed and implemented with the colloidal ice.

Annealing procedures:
Frustrated systems can easily be trapped in a metastable state when the material or sample is cooled to reach the GS. Different annealing protocols have been developed to date, such as thermalisation during sample growth (Morgan et al., 2011) or above the Curie temperature of the constituent material (Zhang et al., 2013), and rotational demagnetization (Ke et al., 2008;Nisoli et al., 2007;Wang et al., 2018). However, most of the theoretically predicted ordered phases in highly frustrated systems are still far from being realised since these techniques do not allow for direct system visualization during annealing, making it impossible to monitor in situ how close or far from the GS the system is as a function of time. The colloidal ice could provide a way to directly follow the annealing process in-situ. The current realization relies on static fields; however, it would be possible to introduce a time dependent field that spins the particles in order to induce tunable anisotropic or even time-averaged attractive in-plane interactions.
Three dimensions: A major effort of the frustrated magnetism community is devoted to the realization of a 3D artificial version of the pyrochlore lattice, such as by using magnetic wires or nanobars oriented at determined angles in order to replicate the tetragonal order (Fernández-Pacheco et al., 2017). Recently, an alternative approach to this goal has been proposed (Mistonov et al., 2013) based on the realization of an inverted colloidal opal filled with cobalt via electrochemical crystallization. This opal was obtained using the colloidal crystal template technique, which is a well-established method in material science (Zakhidov et al., 1998), based on replicating the long-range order of an assembled colloidal crystal in a solid matrix.
Although this type of approach may restrict the lattice geometry , or lead to artifacts due to domain wall pinning within the magnetic network, it represents a rather simple, fast and versatile method of obtaining a highly ordered 3D porous structure. Future attempts at the realization of a 3D colloidal ice system may exploit similar self-assembly techniques. On the other hand, progress in optical manipulation has made it possible to create 3D optical traps for colloids (Grier and Roichman, 2006), and thus it should be possible to create new types of fully 3D spin ice geometries using optically confined colloids. In principle, such geometries could be used to model water ice more accurately, or possibly as a method for realizing deconfined phases.

V. OTHER PARTICLE BASED FRUSTRATED SYSTEMS
Frustrated configurations have also been realized for colloidal systems in the absence of a substrate (Han et al., 2008). For example, when colloids confined to a 2D plane are allowed to buckle into the third dimension, as shown in Fig. 16, frustration emerges since upward and downward buckling are equal energy states. In the triangular lattice naturally formed by the colloids, the buckling process produces a structure similar to an antiferromagnetic Ising model on a triangular lattice (Shokef et al., 2013), and interesting disordered and stripelike patterns appear. While the original Ising antiferromagnet at low temperature features extensive entropy (Wannier, 1950(Wannier, , 1973, the buckled system displays subextensive one, which points towards the presence of glassy dynamics and kinetic arrest. Indeed glassiness in such system has been recently reported , with other features as the emergence of an order by disorder transition (Leoni and Shokef, 2017;Shokef et al., 2011). Other colloidal systems that can exhibit frustration effects include either particles with complex shapes (Brown et al., 2000), or isotropic ones deposited above deformed surfaces (Bausch et al., 2003;Soni et al., 2018). In the first case, the frustration can be tuned by designing different shapes which compete with the confinement and impede crystallization, often producing a disordered state Mason, 2009, 2015). Tuning the shape of the confining surface may also lead to frustration effects between isotropically repulsive colloids. The surface topology will induce the formation of lattice defects as disclinations, releasing energetic stresses arising form the packing on the curved confinement Irvine and Vitelli, 2012;Nelson, 2002).
On a macroscopic scale, geometric frustration has been addressed by arranging classical bar magnets. One significant case is shown in Fig. 17 where an ensemble of interacting millimeter-size magnets is arranged into a kagomé lattice. The unit base of the system are ferromagnetic  Han et al. (2008) rods attached to planar rational units which allows only out of plane angular motion (polar angle α) but not in plane one (the angle θ is fixed). The initial state of the system was prepared by polarizing these "macroscopic spins" along the perpendicular plane (ẑ) via a strong static field. Switching off the field induces a relatively fast (∼ 2s) reorganization process, and the magnets stabilized into a equilibrium pattern filled by vertices with 2-in/1-out and 1-in/2-out. The rotors relaxation process occurred in three steps; namely a relatively fast inertial reorientation where the rotors break the axial symmetry, than a sequence of dissipative librations followed by a final damped oscillations which leads to a nearest neighbor spin correlation s i · s j = 1/3. The authors also showed the exciting possibility to extend the system towards 3D, by stacking different plates composed by the rotors in a tetrahedral-like configuration. Such demonstration man- ifests the ubiquitous character of geometric frustration, which transcends length scale.
Numerous other condensed matter systems can be described as an assembly of particles with repulsive pairwise interactions. Examples include vortices in type-II superconductors (Libál et al., 2009), ions (Pyka et al., 2013), dusty plasmas (Morfill and Ivlev, 2009), skyrmions (Ma et al., 2016), and Wigner crystals (Reichhardt et al., 2001). Any of these systems, when coupled to the correct substrate geometry, could exhibit effective spin degrees of freedom and ice rule.
One of the first proposals for such particle-based artificial ices involved vortices in type-II superconductors (Libál et al., 2009), where considerable progress has already been made in creating various types of pinning arrays to control the vortex ordering. In the vortex system, a series of double well traps can be fabricated by placing two pinning sites very close together, as illustrated in the top panels of Fig. 18. When the sample is nanostructured in this fashion, the thinner parts of the superconductor have lower vortex condensation energy, and the vortex will preferentially sit at the highest points of the double well traps. For N p double well traps, the number of vortices N v is directly proportional to the magnetic field, so that when N v /N p = 1/2, the system is equivalent to the colloidal square ice at half filling. Particle based modeling of vortices in pinning arrays is performed using similar techniques as those described for the colloidal systems; however, the pairwise interaction between the vortices is a modified Bessel function for vortices in a bulk superconductor, and has a ln(r) form in thin film superconductors. Many superconducting samples contain a significant amount of intrinsic random disorder, so Libál et al. (2009) explore the effects of additional quenched disorder on the GS in a square ice pinning arrangement at N v /N p = 1/2, Fig. 18. As the disorder strength δ increases, an increasing number of non-ice rule obeying vertices appear in the GS in the form of grain boundaries, but individual monopoles that are not associated with a grain boundary do not begin to appear until δ ≥ 0.5. The appearance of defects arranged in grain boundaries was subsequently observed experimentally in a square ASI (Morgan et al., 2011).
The superconducting vortex configurations can be observed directly using imaging techniques, but it is also possible to probe the stability of the vortex configurations by applying a current in order to depin the vortices. When the vortex arrangement is highly ordered, the depinning threshold or critical current I c is higher than when the vortices are disordered. Latimer et al. (Latimer et al., 2013) studied the critical currents of superconductors with square ice geometries in experiments. They find a series of peaks in the critical current I c as a function of vortex density. When the temperature is lowered, a prominent critical current peak appears at N v /N p = 1/2, as shown in Fig. 19 for I c versus H/H 1 , FIG. 19 The critical current of the vortex artificial square ice sample, normalized by the zero magnetic field critical current Ic(0), as a function of the magnetic field H at T = 4.8 K (red circles) and T = 4.9 K (black squares). Top inset shows the results obtained at T = 4.72 K. The zero magnetic field critical currents are: Ic(0) = 802 µA at T = 4.72 K, Ic(0) = 613 µA at T = 4.8 K, and Ic(0) = 335 µA at T = 4.9 K. The lower inset shows contour plots of the simulated local magnetic field h for the half matching field, which shows the ground state vortex configuration at T = 4.8K. White circles indicate the positions of vortex-free holes and the dashed red square shows the unit cell of the simulation area. Reprinted with permission from Latimer et al. (2013).
where H 1 is the first matching field. The inset shows a simulation of the flux configuration at the f = 1/2 filling where the square ice rule obeying state is highlighted. For H/H 1 > 1.0, the vortices start to become doubly quantized, and another peak in the critical current appears at H/H 1 = 3/2, which again corresponds to an ice rule obeying GS.
In a series of experiments, Trastoy et al. (2014Trastoy et al. ( , 2015 examined spin ice pinning arrays for the high temperature superconductor YBCO, where thermal effects are important. They found that changes in resistance correspond to vortex configurations that are dominated by a geometrically frustrated energy landscape that favors ice like ordering and frustration, and they find a series of peaks in the transport response similar to what is found in ordered square pinning lattice arrangements. For vortices in a kagomé pinning arrangement, various ordered and disordered arrangements have been observed in imaging experiments, and it has been argued that the long-range vortex-vortex interactions are insufficient to lift the degeneracy of the different vortex configurations in the strong pinning sites (Xue et al., 2017). As the magnetic field is increased, additional vortices become trapped in the interstitial regions surrounding the pinned vortices, which could be a step toward the realization of a stuffed artificial spin ice (Xue et al., 2018). Imaging experiments for vortex configurations on the kagomé ice (Wang et al., 2018) showed the predicted kagomé ice rule obeying states (Libál et al., 2009) both at H/H 1 = 1/2, as expected, and also at the higher field of H/H 1 = 3/2, where the additional vortices occupy the interstitial regions between pinning sites. Interestingly, the ice rule state at H/H 1 = 3/2 is even more ordered than the state at H/H 1 = 1/2, suggesting that the interstitial vortices may play a role in annealing defects in the ice GS. Wang et al. (Wan, 2016) examined vortex ordering and dynamics on a switchable artificial spin ice array that was motivated by the realization of a magnetic rewritable artificial spin ice . In this system, frustration effects can be switched on or off by changing the magnetic charge ordering of the underlying reconfigurable magnetic ice substrate. Different magnetic charge orderings can produce disordered or frustrated vortex configurations, as well as non-frustrated configurations that produce ordered vortex crystals, as illustrated in Fig. 20 for different possible ordered states at B/B φ = 1.0 and the f = 2.0 filling. Here types I, II and III correspond to different magnetic charge ordered states of the magnetic substrate (Wang et al., 2018). Other work on vortex states with rewritable magnetic ice substrates showed that the different patterns can be switched with an applied current .
Square and kagomé spin ice geometries have also been proposed for skyrmions in magnetic dot arrays, where the dots are arranged such that the skyrmion can sit on either side of the dot. A recent work demonstrated that, when the skyrmions are strongly localized or particle like, the behavior of the system is very similar to what is observed in both superconducting vortex and colloidal artificial ices (Ma et al., 2016). When the magnetic field is varied, however, the skyrmions can become elongated and fill the entire dot, destroying the frustration, so the system can exhibit a transition from an ice rule obeying state to square or triangular ordered states (Ma et al., 2016).

VI. OUTLOOK
In this section we provide further general directions that cover different disciplines across condensed matter systems not limited to particle based systems.
Stuffed Spin Ice: In particle based artificial ices, additional particles can be placed in the spaces between the traps of the substrate as well as within the traps themselves. These interstitial particles could potentially modify the ice rule for the particles that are sitting in the trap sites. One study has already provided evidence that the interstitial particles can actually enhance the ordering of the ice state (Xue et al., 2018). It would also be interesting to explore the dynamics of interstitial particles for different ice substrate configurations. For example, an ordered square ice background might produce increased or decreased diffusion of the interstitial particles compared to a degenerate GS ice background.
Driven Dynamics: In ASI the spin degrees of freedom are permanently localized on the magnetic nanoisland. In contrast, in particle based artificial ices it is possible for the particles to hop from one ice substrate trap to another, making it possible to explore the dynamics of a frustrated system compared to that of a non-frustrated one. Protocols that could be considered include sliding dynamics under a dc drive, shaken dynamics under external forces, or simply applying a drive to only one side of the sample but not to the other. This would be equivalent to subjecting a portion of an ASI to a magnetic field while the other portion of the sample experiences no field.
Dynamic Substrates: If optical traps are used to create the ice substrate, than a new type of dynamics involving the traps themselves could be implemented. For example, the trap could be flashed on and off, or the orientation of the double well potential could be rotated. Such protocols could provide new methods to reduce frustration effects. Alternatively, it may be possible to induce new types of frustration effects using dynamical protocols that bring repulsive particles together. Such protocols would create local high energy states, but due to the frustration, there may not be an easy direction in which one of the particles could move away in order to reduce the interaction energy. In particle based systems, dynamic protocols could be introduced readily by using dynamic traps (Bhebhe et al., 2018;Curtis et al., 2002).

Frustrated Active Matter Systems:
Active matter composed of self-propelled particles is a rapidly growing research field, with an explosion in the experimental realization of artificially self-propelled colloids (Bechinger et al., 2016). An interesting area for future exploration is frustrated systems on lattices where the fluctuations are active rather than thermal. Some studies have already shown the organization of active matter into vortex states (Beppu et al., 2017) and vortex lattices (Nishiguchi et al., 2018;Wioland et al., 2016), and it would be possible to arrange these vortices into a frustrated configuration. Additionally, if active and passive particles are mixed together, a frustrated ice geometry might reduce the diffusion of the active particles, while an ordered ice state could produce flow paths for the active particles. Other systems such as coupled gyroscopes or spinners could be placed in frustrated geometries which could lead to interesting dynamics (Nash et al., 2015).
Deformable Substrates: In most of the ASI studied to date, the substrate or confinement is fixed and cannot react to the forces exerted on it by the particles. In soft matter systems there are many ways to create a substrate that could react and be deformed elastically. If this reaction is introduced into an ASI system, it could lead to new methods for reducing frustration effects. In addition, excitations such as monopoles could create additional long range strain fields in the soft substrate that could be attractive or repulsive for other monopoles, leading to emergent effects such as excitons and polarons.
Nonlocal Frustration Effects: Numerous experiments have shown that large numbers of colloids can be addressed individually in optical feedback experiments (Lavergne et al., 2019). Thus, another route for studying frustrated systems would be to introduce non-local effects in which the forces experienced by particles in one part of the system are correlated with the positions of other particles that are far away. This would make it possible to test whether phase transitions can be induced by increasing the strength of the nonlocal interactions, to study the effects of small world versus random interactions, and to determine whether nonlocal effects change the nonequilibrium properties of the frustrated state. It would also be possible to engineer nonlocal frustration effects in which a reduction of frustration at one location in the system increases the frustration at a distant point, and to study how the system would relief such effects.

Continuous Spin Directions:
Another avenue in 2D soft mater systems is to allow the effective spin degree of freedom to point along multiple directions or have some freedom to move in the third dimension in order to create a frustrated Potts model, Heisenberg like spin models, and/or XY models.
Nondissipative Dynamics: In the particle based models considered so far, the motion of the particles is overdamped. It is, however, possible to realize particle ices in which nondissipative effects are important, such as in underdamped systems including dusty plasmas, trapped ion crystals, or even milimetric magnets on 3D printed substrates. The inertial effects could produce phononic excitations, and it would be possible to study the difference in the phonon modes for crystalline versus degenerate or frustrated ice states, as well as the propagation of solitons or shock waves.
Functionalized Colloids for Frustrated Geometries: There are many examples of functionalized colloidal systems in which the interactions between the colloids can be tailored to be anisotropic (Glotzer and Solomon, 2007), reconfigurable (Ortiz et al., 2014), or patchy (Bianchi et al., 2011), or even to mimic chemical bonds (Wang et al., 2012). Using such systems, it would be possible to create 2D or 3D colloidal assemblies that naturally exhibit frustration, or even a colloidal system that mimics the structure of water ice (Pauling, 1935).
Quantum Effects: A final potential avenue would be to place cold atoms or cold ions on trap arrays arranged in an ice structure. This system could be realized at a size scale and temperatures for which quantum effects could become important (Glaetzle et al., 2015;Mazurenko et al., 2017). support from from MINECO 344 (FIS2016-78507-C2), DURSI (2017SGR1061) and Generalitat de Catalunya under Program "ICREA Académia".