Domain dynamics and fluctuations in artificial square ice at finite temperatures

The thermally-driven formation and evolution of vertex domains is studied for square artificial spin ice. A self consistent mean field theory is used to show how domains of ground state ordering form spontaneously, and how these evolve in the presence of disorder. The role of fluctuations is studied, using Monte Carlo simulations and analytical modelling. Domain wall dynamics are shown to be driven by a biasing of random fluctuations towards processes that shrink closed domains, and fluctuations within domains are shown to generate isolated small excitations, which may stabilise as the effective temperature is lowered. Domain dynamics and fluctuations are determined by interaction strengths, which are controlled by inter-element spacing. The role of interaction strength is studied via experiments and Monte Carlo simulations. Our mean field model is applicable to ferroelectric `spin' ice, and we show that features similar to that of magnetic spin ice can be expected, but with different characteristic temperatures and rates.


Introduction
Artificial spin ices [1] are constructed as finite arrays of elongated magnetic dots whose magnetisations are assumed to be well approximated by Ising spins. The stray magnetostatic fields of each island mediate interactions, which are frustrated by geometry. These systems are called 'ices' because minimisation of the magnetostatic energy leads to behaviour resembling that governed by the ice rule for the ground state of solid water, i.e., two-in, two-out spin configurations [2].
The geometry of artificial spin ices can be controlled to a large extent [3][4][5][6][7], and resulting magnetic configurations can be imaged directly using techniques such as MFM [1,[8][9][10][11], PEEM [5,12,13], and electron microscopy [4,14]. Most studies have been made at room temperature with relatively large, thermally stable magnetic elements. This thermal stability is enforced by choosing island volumes large enough that the barrier to magnetisation reversal, primarily associated with shape anisotropy, is much higher than room temperature, so that dynamics can only be induced by applying a magnetic field. Spatial studies are then made with the system in a steady state configuration.
Because the focus to date has been on athermal systems, very little is known about how two dimensional spin ice arrays respond to thermal fluctuations. Recent reports provide evidence that there are magnetic features visible in as-grown arrays that form during early stages of growth [10], namely long range ground state ordering, with well defined domain walls and small clusters representing configurational excitations above the ground state. The origin of the ordering and excitations was argued to be thermal. In another recent work, thermal loss of macro-spin order was reported in a square artificial spin ice patterned δ-doped Pd(Fe) film [15].
From a theoretical perspective, simulation studies [16] of an analogue of square ice constructed from charged colloids in double-well optical traps reveal a freezing transition as the temperature is reduced relative to the height of the barriers between pairs of wells. The same authors also find [17] that strong thermal noise is required to alter the hysteretic behaviour of the system. Recently, Monte Carlo simulations of magnetic square ice have been presented [18], which indicate a sharp peak in the specific heat and a peak in the density of closed loops of spins flipped against the ground state at approximately the same temperature.
In this paper we discuss configurational dynamics of square artificial spin ices in terms of macro-domain growth and boundary movements at finite temperatures. The picture we present is one in which domain boundaries can flow, and also form channels over which magnetically charged vertices move. We find that a domain of ground state configuration should, in the absence of disorder, spontaneously grow until it fills the array in a finite system. We show also that thermal fluctuations can accelerate domain boundary movement, and create clusters similar to those observed experimentally.
We restrict our attention to square array geometries in this work. Spin configurations can be completely specified by vertex arrangements, and for the present work, we use spin or vertex descriptions interchangeably for configurational and energy states. The concept of vertices provides a useful nomenclature for discussing spin configurations [1,8,16,19,20], effective temperatures of field-driven demagnetisation [21,22] and dynamics [11,[23][24][25]. We will use vertex descriptions almost exclusively in what follows. In the square lattice there are sixteen distinct vertex configurations, which can be classified into four groups, Type I-IV on the basis of magnetic charge and total moment. The groups are ordered on the basis of energy. Examples of each type are shown in figure 1(a). The ground state corresponds to a complete chessboard tiling of the two Type I vertices, and is two-fold degenerate. The organisation of the paper is as follows. In the next section we describe general features of domain dynamics in the presence of quenched disorder, using a mean field theory, in reference to conventional magnetic field and zero field cooling experiments. We then discuss effects of thermal fluctuations and strength of coupling using Monte Carlo simulations. These results are compared to experiments. We also discuss how thermal fluctuations can create small cluster excitations above the ground state, which can be understood using statistical arguments. Throughout this work, we see that Type I domain formation and growth are general phenomena, a result which is corroborated by simulation results for artificial 'spin' ice in ferroelectric media, which we discuss in an appendix.

Magnetisation Processes in Mean Field Approximation
Magnetic elements are approximated as block Ising spins, and arranged on a finite two dimensional square lattice, as shown in figure 1(b). The importance of long range dipolar interactions for correlations is a topic of investigation [13,19,26,27], but our simulations suggest that shortcomings of the finite range interaction approximations are significant only in idealised perfect systems. In the models used throughout the main text of this paper, we consider only the three nearest neighbour interactions: J n , J nn , and J nnn in a point dipole model. The square array geometry can be thought of as two sets of parallel lines of elements, with the sets aligned at 90 • to each other. The J nn accounts for interactions along each line, and the J nnn couples elements from adjacent parallel lines. The J n couples elements from the two different sets. These couplings are indicated schematically by dashed lines in figure 1(b). In the appendix, we give results for a ferroelectric ice where all point dipoles in the array are summed over, and find qualitatively similar results to those presented in this section.
The energy of a magnetic element is determined by the local interactions and any applied fields, H. Denoting the magnetic moment at site i as m i , reversal of m i will occur when Here, µ B is the Bohr magneton, m i is the mean field thermal average moment at site i, and ǫ c is an energy barrier to reversal. We suppose that the energy barrier to reversal is proportional to the anisotropy responsible for the Ising like behaviour of the magnetic elements. In the mean field approximation this is where K represents an anisotropy barrier. Each magnetic element is assumed to behave as a single 'soft' macrospin, with a total moment that depends upon temperature. We approximate the temperature behaviour of an element's moment by the Langevin function L(x): β is the inverse temperature 1/(k B T ) where k B is the Boltzmann constant. Previous studies have shown that, in field-driven dynamics, a perfect system allows nucleation of Type I domains on a Type II background through Type III vertex dynamics that start at array edges where moments have fewer neighbours and therefore lower reversal thresholds [24]. Local variations in the coupling or reversal barrier parameters can facilitate nucleation processes to start inside the array, leading to less sensitivity to the array boundaries [20,28,29]. In a simple model for switching barriers, the barrier height is proportional to the element volume. One would expect that during early stages of deposition of material, the element thicknesses are distributed over a range, leading to a spread in switching field values. Furthermore, significant distributions in switching fields have been observed experimentally in fully-grown, athermal artificial spin ice [12,[29][30][31][32] and switching field disorder has been shown to give similar outcomes in numerical simulations as other disorder types [33]. We use this as motivation for describing disorder in terms of the reversal barriers. We assume the barriers vary randomly in an inverval ∆ centred about K, that is, The mean field model allows one to examine the stability of the ground state relative to temperature but cannot capture correlations such as those involved in avalanche processes, and such processes will be discussed later. Initially, all moments are assumed to be of unit magnitude in some specified configuration. In the present model, an element i is picked at random, and the time averaged reduction of its magnetisation is calculated according to (3). If the magnetisation state becomes unstable, that is, if (1) is satisfied, then m i is set to −m i . The process is repeated by choosing a new moment at random, and the iterations continue until a steady state is found such that no m i value changes by more than 10 −4 .
Reduced units used throughout the paper are defined so that |m i | ≤ 1, J n = 1.5J o , J nn = 0.7J o , and J nnn = 0.3J o . The reversal barrier has been chosen to be, in reduced units, K = 10J o . Arrays consist of 640 vertices, and array boundaries are taken to be 'open' boundaries for which the edge elements have only three nearest neighbors, and corner elements have two nearest neighbours, as shown in figure 1 The stability of the Type I ground state can be studied by starting the system at low temperature with a complete Type I tiling, and then increasing the temperature. The critical temperature can be calculated within mean field theory from the average magnetisation of an element in a Type I tiling, which obeys where h loc = 2(2J n − J nn + J nnn ) + 1 2 K. The critical temperature T c , is given by 1/β = h loc /3, to first order approximation. Using the parameters listed above, Because changes in magnetic configurations are driven by instabilities only, in the mean field model, the transition is very sharp with a well defined critical temperature in a perfect system without disorder. Disorder broadens the transition region in temperature, and leads to formation of domains of the two different possible Type I tilings separated by Type II and III chains. Examples are shown in figure 2 for Type I populations n 1 as functions of temperature for different disorders ∆. The population shown is the sum of both Type I tiling possibilities. Other authors also report a broad melting transition in a thermal ice prepared in a magnetised (Type II vertex tiling) state [15].
With disorder, a sharp initial reduction occurs in n 1 near the critical temperature, corresponding to the nucleation and growth of domains. This is followed by a long tail in which the total population of Type I vertices approaches the value n 1 = 1/8, corresponding to the high temperature limit in which each of the 16 possible vertices appear with equal probability. Note that the transition temperature goes to the expected value of 3.4T o as ∆ → 0.
Also shown in figure 2 is the Type I population that evolves from a random vertex configuration starting at low temperature. This is analogous to the 'zero field cooled' state of a ferromagnet that, initially at some temperature above the ordering temperature, has been quenched to low temperature. In the mean field model, a state with large Type I domains immediately evolves at low temperature if K is sufficiently small relative to the J's. The n 1 populations grow with increasing temperature as other vertices become unstable. This growth of domains occurs through the elimination of small domains and clusters. An example is shown in figure 3, where vertex configurations at two steps during the numerical iteration are shown. The temperature is T o , and disorder ∆ = 0.1. Green squares represent Type I vertices. The blue arrows are Type II vertices, and their direction indicate the orientation of the local Type II moment. The light and dark red crosses represent the two charge flavors of the Type III vertices.
The first image (left) is taken after 100 iterations after an initially random configuration of moments. Already a well defined Type I domain structure has formed. The domains are separated by chains composed of Type II and III vertices with net magnetic moment and charge, and the energy of their formation acts as a surface tension on the Type I domains. For this reason, the chains tend towards straight in the limit of no disorder, since corners lengthen the chain and can involve Type III's. The second image (right) shows the evolution after 100 additional iterations. The longest chains have shifted slightly, and positions of the Type III vertices within the chains have changed. Most notably, the smallest domains have vanished, and other small domains are reduced in size. Note that Type IV vertices are highly unstable, and do not survive past the first few iterations.
The reduction in size and eventual annihilation of small domains is the result of an energetic biasing of random domain wall motion towards motion that makes walls shorter by moving closed walls inwards. Dynamics on a Type I/Type II vertex background -such as domain wall motion -can be described in terms of Type III vertex motion [11,23,24]. There are therefore two energy costs associated with dynamics: the Accordingly, the lowest-energy mode of domain wall motion is the propagation of Type III vertices that have formed at random positions during wall creation. These random fluctuations do not favour either growth or shrinking of domains on average. The existence of random fluctuations lead us to speculate that an analogy to creep motion [34] may be possible.
The cost of Type III vertex creation depends on the local spin configuration, and is lowest at corners in the domain walls. For example, in figure 4(c), the easiest spin to flip is the circled spin at the domain wall corner, which can flip to create a Type III vertex pair. The resulting Type III vertices can propagate at approximately zero cost by moving along the domain wall, flipping a diagonal chain of spins on the inside of the wall. This process moves the domain wall inwards, as illustrated in figure 4(d). Such sequences of Type III creation and propagation are the domain wall motion mechanism with second-lowest energy cost. Thus, over time, the wall motion is biased so that despite random fluctuations, all closed domains will disappear.
As a final comment on the mean field results, the local fields along the element magnetisations are weakened within the chains, compared to the fields in regions of lower-energy Type I vertices. In the mean field model this leads to a suppression of the magnetisation in the elements participating in the chains. The magnitude of the a b c d  suppression depends upon the temperature, and can be very large as one nears T c . This effect is shown by the vector magnitude map in figure 5. The magnitude of element magnetisation is shown in grayscale, with black being unity and white being small. The domain pattern corresponds to the configuration shown in the second panel of figure 3. One sees that the moment is considerably reduced within the chains, and especially for small clusters. We note that instabilities occur when the moment of individual elements vanish, and these are most likely to occur in the Type II walls.

Fluctuations and Disorder
The energy differences between configurations are often small, even for configurations that are very different. For this reason, single spin flips can drive avalanches that lead to large configurational changes and strongly modify vertex populations as domains evolve. To capture these dynamics, we turn to Monte Carlo simulations. The model is based as before on the point dipole approximation assumed for the mean field model. The form of the energy used is the same, but evolution is described using a heat bath algorithm, rather than a simple critical field defined by (3). Additionally, we consider the case of a large Curie temperature so that the effective block moments are insensitive to temperature. In this limit the energy used for the heat bath algorithm is where now all magnetisations have fixed magnitude as in an Ising model, with |m i | = 1. The barrier to reversal is determined by the anisotropy constant K. The parameters used are the same as those for the mean field model. The results shown below were made with 10, 000 Monte Carlo steps at each temperature and averaged over 20 realisations of disorder at each disorder strength.
As in the mean field calculation, the role of disorder is largely to broaden the transition. Example heating curves are shown in figure 6, where thermal evolution of ground state and random tilings are shown, for different disorder strengths ∆. The most striking feature, in comparison to the analogous mean field results shown in figure 2, is the smooth behaviour of the population at the transition, and the rounding in the transition region.
The dependence on disorder is also less dramatic than in the mean field case. In both mean field and Monte Carlo simulations, the spatial distribution of switching barriers is weakly connected to the final distribution of domains and chains. An analysis of the distribution of K values shows that there is a propensity for chains to locate on strong pinning sites, where K is large, and for weak pinning sites to lie more frequently inside domains. This result can be expected because domains grow via thermally activated spin flips. Growth is accomplished through distortions of chains, and chains will pin at sites where K values are too large to allow thermally driven element magnetisation flips. Also, as a consequence, the final resulting domain configuration at a given temperature displays rough chains for large ∆ values, and flat segmented chains for low ∆ values. Differences between Monte Carlo and mean field results appear because of how thermal fluctuations drive domain growth. Fluctuations allow the system to explore dynamical pathways different to those dictated by a stability analysis, and are therefore quite significant in determining the final domain configuration. At temperatures below the transition, the dominant paths tend towards the lower energy, single domain state.
An example evolution is shown in figure 7, using the colour scheme defined in figure 3. Four configurations are shown at constant temperature T = T o , starting with a random initial configuration. The first three configurations are taken after 1000, 2000, and 3000 steps. The fourth configuration (in the bottom right hand corner) is taken after 10, 000 Monte Carlo steps. The amount of disorder in this example was small relative to both K = 10J o and J n = 1.5J o , with ∆ = 0.1. One sees the slow growth of one Type I domain at the expense of another, through the motion and deformation of Type II and III chains. The dynamics is equivalent to that seen with the mean field model. Note that here also Type IV vertices do not survive the first few Monte Carlo steps unless interactions are weak (as with large spacings). Blue arrows indicate Type II vertices, and red crosses are Type III vertices. The first three panels show configurations after 1, 000, 2, 000, and 3, 000 steps for a temperature of T o and disorder ∆ = 0.1. The initial state was random. The final panel shows the configuration after 10, 000 steps. Note that small clusters and domains disappear as large domains grow at the expense of smaller domains.
In addition to our mean-field and Monte Carlo simulations, we have imaged vertex configurations in as-grown samples, and found domain structures similar to those in simulations. An example MFM image is shown in figure 8, where domain walls, which carry net moment are clearly visible on the approximately neutral background of Type I domains. The sample imaged here was previously studied in terms of ground state ordering and thermal excitations [10], and consists of an array of nominally 280×85×26 nm 3 Permalloy islands, with a 3 nm Ti buffer and a 2.5 nm Al cap. Full details of sample fabrication are given elsewhere [10], but the key point is that islands are formed by deposition of Permalloy through gaps in a nanopatterned resist mask, so that island thicknesses -and hence barriers to magnetisation reversal -grow over time.
For the simulation, the values shown in the figure were obtained after a finite number of Monte Carlo steps, and the effect of temperature, relative to the local coercive field, is to control the rate at which n 1 changes. In the experiment, the deposition rate relative to temperature probably controls the rate at which n 1 changes. We suggest that domain configurations found in the simulation and the experiment can both bo interpreted as snapshots in time of a thermally driven evolution. We note that a difference between the simulations shown in this example and the experiment does exist, in that the as-grown samples support small clusters of flipped spins. We discuss this further in section 3.2.

Spacing and interaction strength
The preponderance of Type I vertices in domain structures is a consequence of the higher energies needed to form other vertex types, and the corresponding lower probability for their creation. The probability for reversing a moment depends upon the size of local fields relative to thermal energy. This ratio also determines the rate at which domains can grow.
An example from Monte Carlo simulations is shown in figure 9(a). Here the Type I population is shown as function of interaction strength. The interaction strength is scaled by a factor s such that J α → sJ α /J o for each of the three interactions α. Temperature is fixed at T o in all cases. Several different disorder strengths were studied and the simulations were run for 100, 000 steps at each spacing. The n 1 population tends to unity as the interaction strength increases, although the rate at which this occurs depends on disorder strength. For strong disorder, the rate is slowest. We note also that there were no significant Type IV populations observed, although the population did increase with increased disorder. We have also measured statistics of the vertex type populations for as-grown samples with different lattice constants. In these studies, the dimensions of individual elements were held constant at 270 × 115 × 25 nm 3 , but the lattice constant was varied from 400 nm to 600 nm. In this way the interaction strength should decrease, roughly like 1/r 3 where r is the lattice constant. Results are shown in figure 9(b), where the n 1 population is given as a function of 1/r 3 . All samples in the series consist of 0.5×0.5 mm 2 continuous arrays of Permalloy islands grown on a 2 nm Ti buffer, with no capping layer, and were fabricated in a single batch, to ensure that fabrication parameters were consistent across the series and so avoid problems with quantitative reproducibility (such as those mentioned in [10]).
Both figures 9(a) and (b) show that the Type I populations increase with increasing interactions. In principle one should be able to extract estimates of interaction energies as a function of spacing from the experimentally measured n 1 . However, inspection of figure 9(a) shows that one also needs to know the amount of disorder before a prediction for n 1 can be made. The disorder can be expected to be strongly dependent on details of the sample growth and design [32].

Fluctuations and clusters
The role of disorder is subtle. On the one hand, disorder in K leads to randomness in the chains. On the other, disorder also facilitates fluctuations and leads to appearance of small clusters of Type II and III excitations, examples of which are illustrated in figure 10. In Monte Carlo simulations, these fluctuations usually appear and disappear quite rapidly (in Monte Carlo step time), but can also lead to structures that may persist through several Monte Carlo steps. Experimentally, as-grown samples are seen to exhibit frozen-in fluctuations, as seen in figure 8. In previous studies, it was shown that the population of clusters decays approximately exponentially with their energy, which was given as evidence of thermally-driven processes occurring during sample growth [10]. In simulations of field-annealed square ices realised in nanostructured superconductors, small structures can only exist independently of domain walls when sufficient disorder is present [20].
An analytical description of the nucleation, growth and decay of an isolated cluster gives insight into the distribution of experimentally observed clusters. As seen in figure 10, clusters can be thought of as connected lines of spin flips against the ground state. The lines may branch. In our analytical model, the growth and decay of clusters proceeds by single spin flips that extend or shorten these lines. For example, the left hand cluster shown in figure 10 can evolve in to the right hand cluster, and vice-versa. For simplicity, we neglect disjoint clusters: large clusters cannot evolve into two smaller clusters.
Under the assumption that clusters are excitations above the ground state that are nucleated by random thermal spin flips, we start with an initial perfect ground state and study the temperature-dependent evolution of small clusters of up to four spin flips. The clusters may be grouped by their topology (and hence energy). There are twenty five groups of such clusters: apart from the ground state, there is one distinct cluster of one flip, two clusters of two flips, five of three flips and sixteen of four flips. These groups include clusters that contain Type IV vertices, but these high-energy clusters are suppressed relative to those containing only Types II and III vertices.
The processes of cluster growth and decay are transitions between the groups of clusters, and can be described by a master equation for the probability, P (A), of a where m is the island moment and r is the lattice spacing) for cluster excitations of up to 4 spin flips. As indicated in the legend, data are shown from experiment (these were previously published in [10]), the master equations (6) and Boltzmann factors. The labels indicate the mnemonics used in [10]. The labels 4+ and 4t are in parentheses to indicate that those clusters were not observed experimentally in [10]. Error bars on the experimental data correspond to the square root of the population of each cluster type.
cluster having topology A: where the sum runs over all cluster topologies that a cluster of type A can evolve to or from by a single spin flip. G(A → B) is a multiplicity that takes into account the number of 'pathways' from A to B, and depends on the numbers of orientationallydistinct clusters with a given topology and the number of clusters of topology A that can evolve into a particular cluster of topology B. The rate ν(A → B) is given by the Arrhenius-law factor ǫ c − ǫ dip gives the energy cost of a spin flip to grow or shrink a cluster on a Type I background. ǫ c is an arbitrary anisotropy barrier whose exact value does not affect the steady state cluster probabilities. In this model, disorder in anisotropy barriers is neglected, and ǫ c is constant always. ǫ dip is the energy of the spin prior to flipping, due to point-dipole interactions among near neighbours (J n , J nn and J nnn ). The attempt frequency f does not affect the steady state probabilities. The twenty five coupled equations of the form (6) can be solved numerically, with an initial condition that the system is in the ground state: P (∅, t = 0) = 1, P (i, t = 0) = 0, ∀i = ∅, where ∅ represents the ground state. The steady state probabilities of all clusters that do not contain Type IV vertices, normalised to the experimentally observed population of the 1 excitation (that is, a single spin flip), are plotted in figure 11, as a function of their energy difference with the ground state (as given in [10]). The agreement between theory and experiment is remarkable.
The figure shows the solution at k B T = 2.96u, the temperature at which the exponential decay that best fits the theoretical results matches the exponential decay of experimental data. The energy scale u = µ 0 m 2 /(4πr 3 ), where r is the lattice spacing and the moment m = MV , where M = 860 × 10 3 Am −1 is the magnetisation and V is the island volume. The clusters that contain Type IV vertices are not shown, but apart from two exceptions they all have populations less than half of the lowest population shown. This is consistent with experiments [10], where such clusters are not seen and it was concluded such clusters were highly improbable. The exceptions are a Type III -Type IV -Type III line and a Type III -Type IV -Type II -Type III cluster, which both have predicted populations of ∼ 2.
For comparison, the figure also shows populations predicted from Boltzmann factors of the form g exp(−dE/k B T ). The degeneracy g of a cluster topology is determined by the number of distinguishable ways a particular shape can be rotated and reflected, as well as a factor of 2 (equal for all clusters) to account for the possibility of a global spin flip. dE is the cluster's energy above the ground state, taken from [10]. The ratio of interaction strength to temperature has been tuned to match the exponential decay to that of the experimental data, giving a ratio of k B T = 8.18u.
In experiments, configurations are frozen in as the island volume becomes large enough that shape anisotropy barriers suppress island magnetisation reversal. The ratio of temperature to nearest-neighbour coupling can be compared for the master equations and the Boltzmann distribution by estimating the thickness at which freezing-in occurred from the two models. If we estimate the temperature during growth to be 350 K, then the nearest-neighbour interaction is 2.4×10 −21 J for the master equation and 5.9×10 −22 J for the Boltzmann factors. For Morgan et al.'s 280 × 85 nm 2 islands, the estimated interaction strengths correspond to thicknesses of 3.5 nm (master equation) and 1.7 nm (Boltzmann factors). Both these estimates are of the same order of magnitude as the estimate of ∼ 1 nm given in [10].
Interestingly, while the Boltzmann factors and the solutions to (6) both agree well with the experimental data, the Boltzmann factors agree best for the excitations 1, 2L, 3Z, 4Z while the solutions to (6) fit better for the other clusters. For example, the 'pathway' 1 → 2L → 3U → 4O leads to a stable configuration, 4O. In (6), we have truncated the maximum cluster size to 4 flips, making 4O very stable. However, even without this approximation, the energy cost of adding an extra spin flip to 4O is high, and the cluster is stable. Thus, the 4O cluster is in some sense a 'sink' of probability, and has higher probability than would be expected from Boltzmann factors, which is captured by the master equations.
Finally, we note that the master equations (6) explicitly assume an initial ground state configuration, and treats small excitations as having nucleated and grown on this background before being frozen by increasing island switching barriers. The clusters in the as-grown samples of Morgan et al. might also be interpreted as regions which have fallen out of equilibrium during the effective cooling of an initially random configuration.
While we have not explicitly tested such a possibility, the agreement of the master equations with experiment gives further evidence for the picture of nucleating and growing clusters.

Conclusion
Vertices in artificial spin ice are local configurations of magnetisations that can be thought of to some extent as classically emergent objects analogous to quasiparticles, able to move and interact according to well defined rules dictated by the lattice geometry. We have shown that domains of ground state vertices form spontaneously in the square lattice. Domain boundaries are defined by chains of Type II vertices, along which Type III vertices can appear, move and drive growth of one domain at the expense of another.
We have shown that these domain dynamics occur with rates that depend on temperature and the strengths of interactions between elements. Effects of disorder have been studied, and shown to affect mainly the average size and growth rate of domains, while not modifying significantly the fundamental processes involved. A comparison to experimental results was obtained from as-grown samples for which interelement spacing was varied. The experiment and simulations are essentially different in that element thicknesses are changing with time in the experiment whereas temperature is changed in the simulation. Nevertheless, similarities between populations for the Type I ground states were observed between experiment and simulation, suggesting that the asgrown samples are behaving somewhat as frozen snapshots of magnetic ordering within the arrays occurring as described by the simulations presented here.
In addition to domain wall motion, there is also the possibility that an element will flip within a domain. Reversal of a single element in a Type I configuration creates a pair of oppositely charged Type III vertices. Additional reversals lead to clusters of Type II and Type III vertices, that may persist for some time. We have calculated the probabilities of occurrence of small clusters as a function of their energy, and shown that one obtains a distribution in quantitative agreement with that reported recently from experiments [10].
Lastly, we have shown that all features of domain formation and growth can be obtained with a generic mean field model. This leads us to suggest the possibility of creating artificial spin ice using ferroelectric media, and we have provided an example in the appendix using parameters appropriate to bulk PbTiO 3 . In this example we also showed that models using full dipole sums over a lattice of point dipoles produces the same qualitative results as a model using a severely truncated sum.

Appendix: Ferroelectric Ice
Our discussion so far has been within the context of magnetic spin ice, although our results have general applicability to a number of different systems. We illustrate this with a return to the mean field model, and discuss now domain growth within the context of Landau-Ginzburg theory for ferroelectrics.
A second order Landau-Ginzburg model is used to calculate the electric dipole moment of an island i, given the effective electric field E at that point. Each dipole is assumed to lie along only one axis, parallel to the long axis of that island. The ferroelectric free energy density F at island position r i is given by: where P is the electric polarisation, α = A(T − T c ) and β are the Landau parameters for the material. P(i) is assumed to lie parallel with the element long axis. A point dipole approximation is again used with the electric dipole associated with island i denoted as p i = V P(i), where V is the volume of an element and P(i) is assumed to be uniform across the element. Unlike the previous mean field theory, in this model we calculate the electric field at position i by superposing contributions from all other elements (approximated as point dipoles) in the square lattice. The number of elements used is sufficiently small that one can use the simple summation Minimisation of (8) will, for a range of E values, give two solutions for p i with opposite sign. One can show simply from (8) that the energy barrier separating the two solutions disappears when This is the condition for stability of an element's polarisation orientation. When this condition is fulfilled, the polarisation is reversed. As in the previous mean field model, this procedure is iterated through a system of elements at a fixed temperature until a steady state is reached. Parameters used for the simulations are: T c = 1100 K, A = 7.5×10 5 C −2 m 2 N K −2 , and β = 2.4 × 10 9 C 4 m 6 N. These parameters give a polarisation in zero field and at room temperature P = (α/β) = 0.5 C/m 2 , which is typical for PbTiO 3 in bulk [36] and is quite accurate for elements that are over 100 nm long. It should be noted that the Landau expansion (8) should strictly speaking be to sixth-order to most accurately model the ferroelectric phase transition [36] but we neglect the P 6 terms here for simplicity. Results are given for a 20 micron square sample with N = 50 islands along one edge (4900 islands and 2401 vertices). A spacing of 400 nm between island centres is assumed. The starting configuration for an iteration is taken to be islands having polarisation with constant magnitude P i = ( α β ) 1/2 , with random alignments along element axes.
The energy barrier with these parameters for flipping an element polarisation is much larger than the thermal energy, for most temperatures below the Curie point, so reversal via thermal fluctuations occurs only in the vicinity of T c .
The general features observed in the magnetic system are found also for the ferroelectric system. Starting from a macroscopically upolarised state at low temperature, the effect of increasing temperature is to increase the number of Type I vertices. As before, domains of Type I's form, separated by chains of Type II and III vertices. Effects of disorder are also completely analogous.
Finally, we note that domains grow at temperatures approaching T c , as also found for the magnetic system. An example of average Type I domain size is shown in figure  12 for three temperatures. The error bars correspond to the spread in sizes calculated over ten different initial configurations.
It is useful to note that it is possible to define a relaxational dynamics based on the effective field − ∂F ∂p i . The iteration process can thus be related directly to a real time dynamics. In this sense, the numerical iteration method is a time evolution, and these simulations give a feeling for the effect that different cooling rates may have on a system. If the system is cooled faster than the islands' polarisation can respond, then a higher number of unfavourable Type II and Type III vertices will be frozen into the artificial ferroelectric ice. This is interesting because these domain walls form without the presence of disorder and have a density depending on the history of the system.