Dynamics of artificial spin ice: continuous honeycomb network

We model the dynamics of magnetization in an artificial analog of spin ice specializing to the case of a honeycomb network of connected magnetic nanowires. The inherently dissipative dynamics is mediated by the emission, propagation and absorption of domain walls in the links of the lattice. These domain walls carry two natural units of magnetic charge, whereas sites of the lattice contain a unit magnetic charge. Magnetostatic Coulomb forces between these charges play a major role in the physics of the system, as does quenched disorder caused by imperfections of the lattice. We identify and describe different regimes of magnetization reversal in an applied magnetic field determined by the orientation of the applied field with respect to the initial magnetization. One of the regimes is characterized by magnetic avalanches with a 1/n distribution of lengths.


Introduction
Spin ice [1,2] is a frustrated ferromagnet with Ising spins that possesses rather peculiar properties. First, as a consequence of strong frustration, it has a massively degenerate ground state and retains a finite entropy density even at very low temperatures [3]. Second, its low-energy excitations are neither individual flipped spins, nor domain walls, but are point defects acting as sources and sinks of magnetic field H [4,5]. The concept of magnetic charges, while not exactly new [6,7,8], has proven very useful in elucidating the static and dynamic properties of spin ice [9,10,11,12,13]. It is worth noting that these objects are magnetic analogs of excitations with fractional electric charge found in the familiar water ice [14].
Artificial spin ice is an array of nanomagnets with similarly frustrated interactions. The original system made by Schiffer's group had disconnected elongated islands (80 nm by 220 nm laterally and 25-nm thick) made of permalloy and arranged as links of a square lattice [15]. Later versions included a connected honeycomb network of flat magnetic wires [16,17,18,19], in which the centers of the wires form a kagome lattice, hence the sometimes used name "kagome spin ice" [17]. Whereas it had been originally intended as a large-scale replica of natural spin ice, it became clear very soon that artificial spin ice has a number of its own peculiar features. For example, because the magnetic moments in artificial spin ice are extremely large, on the order of 10 8 Bohr magnetons, the energy scale of shape anisotropy due to dipolar interactions, 10 5 K in temperature units [20], effectively freezes out thermal fluctuations of the macrospins meaning that the system is not in thermal equilibrium. Dynamics of magnetization has to be induced by the application of an external magnetic field [15]. Elaborate experimental protocols involving a magnetic field of varying magnitude and direction [21] have been proposed to simulate thermal agitation invoking parallels with fluidized granular matter. It remains to be seen whether the induced dynamics yields a thermal ensemble with an effective temperature. The analogy with granular matter is further reinforced by recent observations of magnetic avalanches in the process of magnetization reversal [18,22].
In this paper we present a model of magnetization dynamics in artificial spin ice subject to an external magnetic field. Two sets of physical variables are used: an Ising variable σ = ±1 encodes the magnetic state of a spin, whereas an integer q quantifies the magnetic charge of a node at the junction of several spins. Magnetization dynamics are mediated by the emission of domain walls carrying two units of magnetic charge from a lattice node, their subsequent propagation through a magnetic element, and absorption at the next node. We specialize to the case of kagome spin ice, in which magnetic elements form a connected honeycomb lattice [18,16,17,19]. The model can be readily extended to other geometries and lattices with disconnected magnetic elements [15,22,23,24]. Some of the results presented here have been outlined previously [25].

Basic features of the model
Our model is specialized toward an experimental realization described previously [17]. That artificial spin ice is a connected honeycomb network of permalloy nanowires with saturation magnetization M = 8.6 × 10 5 A/m and the following typical dimensions: length l = 500 nm, width w = 110 nm, and thickness t = 23 nm. Three nanowires come together at a vertex in the bulk. At the edge of the lattice, a vertex may have one or two links coming in.

Basic variables: magnetization and magnetic charge
We label nodes of the lattice by a single index i and nanowires connecting adjacent nodes by the indices of its two nodes, ij. In equilibrium, the vector of magnetization M points parallel to the long axis of the wire, so we can encode the two states of a nanowire by using an Ising variable σ ij = ±1. In our convention, σ ij = +1 when the vector of magnetization points from node i into node j. This definition implies antisymmetry under index exchange, σ ij = −σ ji . We define the dimensionless magnetic charge at node i as where the sum is taken over the three neighboring sites j. This definition is quite natural: since magnetic induction B = µ 0 (H + M) is divergence-free, the magnetic charge Q i of node i equals the flux of magnetic field H out of the node, which in turn equals the flux of magnetization M into it: Thus q i is indeed magnetic charge measured in units of M tw. The Bernal-Fowler ice rule [2] enforcing minimization of the absolute value of charge |Q i | is usually justified from the energy perspective: the magnetostatic energy of spin ice can be written as the energy of Coulomb interaction of magnetic charges, The dominant second term-the charging energy of a node-forces minimization of magnetic charges in natural spin ice. The "capacitance" C is determined by the dipolar and exchange couplings energies of adjacent spins [5].
Although we will see below that these energy considerations are not relevant to artificial spin ice within our model, for the moment we will simply adopt the result to it. In honeycomb ice, where the coordination number is 3, dimensionless charge q i can take on values ±1 and ±3. Minimization of node self-energy would select states with Indeed, triple magnetic charges have never been observed in our samples of honeycomb ice. Ladak et al. [18,19] have found nodes with triple charges. The difference is likely due to a higher amount of quenched disorder arising from random imperfections of the lattice [26] in the samples of Ladak et al.
We will find it convenient to use the following notation. A site with a unit charge q i = ±1 has two majority links with σ ji = q i and one minority link with σ ji = −q i . For site i in Figure 1, the minority link is ij.

Basic dynamics: emission of a domain wall
To reverse the magnetization in a nanowire, one must apply a sufficiently strong external magnetic field. The reversal begins when one of the nodes, say i, emits a domain wall (w) into link ij, Figure 2(a). If the link initially has magnetization σ ij = ±1, a domain wall can traverse it from i to j only if it has charge of the right sign, i.e., q w = 2σ ij = ±2. Once the domain wall passes through the link, σ ij changes its sign. Now a domain wall with the same charge q w can only traverse the link in the opposite direction.
The critical field H c , at which a domain wall is emitted from a node, can be estimated as follows [27]. Suppose a node with magnetic charge q i = ±1 emits a domain wall with magnetic charge q w = ±2 [8,25]. Conservation of magnetic charge means that the charge of the site turns to q i = ∓1. The emission process can thus be viewed as pulling a charge q w = ±2 away from a charge of the opposite sign q i = ∓1. The maximum force between the two charges is achieved when the separation between them is of the order of their sizes a, which is roughly equal to the width of the wire w: F max = µ 0 |Q i Q w |/(4πa 2 ). This force must be overcome by the Zeeman force applied to the domain wall by the external magnetic field, F ext = µ 0 |Q w |H ext . Hence the estimate of the critical field, For the system parameters used in our previous work [17] and listed above, this estimate yields µ 0 H c = 18 mT. The critical value observed experimentally [28] is 35 mT. One can envision another possible process, wherein the reversal is triggered when a site with charge q i = ±1 emits a domain wall of charge q w = ∓2 and change its charge to q i = ±3. Considerations along the same lines as above show that the critical field required to pull apart charges q i = ±3 and q w = ∓2 is 3H c . As we will see below, magnetization reversal in samples with low quenched disorder occurs well before the external field has a chance to reach this value. This explains why triple charges are never generated as a result of the emission of a domain wall.
The estimate for the critical field was obtained under the assumption that the external magnetic field H ext is applied along the link into which the domain wall is emitted. When the field makes angle θ with the link, it is reasonable to suppose that only the longitudinal component of the field H ext cos θ pulls the domain wall away from the node. We thus expect the following angular dependence of the critical field: As we will see later in Sec. 3, our educated guess is almost right and that Eq. (6) requires only a minor correction: the angle θ should be measured not from the axis of the link but from a slightly offset direction. This effect is caused by an asymmetric distribution of magnetization around a node, which was missed by the simplified, mesoscopic model of this section.

Basic physics: absorption of a domain wall
Once a domain wall is emitted into link ij, it quickly propagates to the other end of the link, toward node j. Theoretical and experimental studies of domain wall motion in permalloy nanowires [29,30] show that walls move at speeds of the order of 100 m/s in an applied field of just 1 mT. This corresponds to a propagation time of the order of 10 ns, which is too short to be observed in most current experimental setups.
When the domain wall reaches the opposite end of the link ij, its further fate depends on whether the magnetic charge at node j has the same or opposite sign of magnetic charge. We consider the two cases in turn.
If the domain wall and node j at which it arrives have opposite charges, q w = ±2 = −2q j , as in Figure 2(a), the domain wall is attracted to the node. It is easily absorbed by the node, whose charge changes to q j = ±1. A new domain wall with the same charge q w = ±2 may be subsequently emitted into one of the adjacent links jk if two conditions are met: (i) the link has the right direction of its magnetization, q w = 2σ jk and (ii) the external field is sufficiently strong to trigger the emission.
Note that condition (ii) is sensitive to the orientation of the field relative to link jk. It also rests on an implicit assumption that the critical field for a new domain wall is not affected by the just completed absorption of the previous one. This assumption is reasonable if the dynamics of domain walls are strongly dissipative and the energy generated during the absorption process is quickly dissipated as heat. Experiments with domain walls in nanowires indicate that they possess non-negligible inertia [8], and therefore our assumption of strongly overdamped dynamics may not be fully justified. Nonetheless, for the sake of simplicity, we shall assume that the dynamics are strongly dissipative and that the extra energy brought by the arrival of a domain wall does not by itself cause the emission of another domain wall from the same node.
Consider now the other case, where the domain wall and the arrival node have charges of the same sign, q w = ±2 = 2q j , as in Figure 2(b). The two charges now repel and the repulsion grows stronger as the domain wall approaches the node. Under the assumption of overdamped dynamics, the wall stops when the Coulomb repulsion between the charges reaches the level of the Zeeman force from the external field. One might think that this may be an equilibrium situation, but we show as follows that this is not the case. The arriving domain wall generates a strong field at the node, whose magnitude is easy to estimate. Since the domain wall is in equilibrium, the force applied to it by the external field, F = µ 0 |Q w |H c , is balanced by the Coulomb repulsion of the node. By Newton's third law, the domain wall applies an equal force to the node. The field created by the wall at the node is This field is added to the externally applied field H c . The resulting field is sufficiently strong to trigger the emission of another domain wall from the node. (This works for any relevant direction of the applied field.) The charge of node j changes sign, q j = ∓1 = −q w /2, and subsequently absorbs the stopped domain wall.

Basic physics: quenched disorder
Imperfections of magnetic links and junctions create local variations of the critical field H c . If the variations of H c result from a large number of small errors, one expects a Gaussian distribution of critical fields ρ(H c ) with a meanH c and a width δH c given by In the limit of strong disorder, when the distribution width δH c is comparable to the averageH c , nodes with the highest critical fields may fail to follow the scenario shown in Figure 2(b) and remain in a state with a triple charge until the field becomes strong enough. Nodes with triple charges have been observed by Ladak et al. [18,19]. In contrast, other samples have never shown triply charged defects [17], indicating that these samples are in the low-disorder limit, δH c ≈ 0.04H c [28]. The distribution width δH c can be compared to another characteristic field strength, the magnetic field generated by an adjacent node, H 0 = M tw/(4πl 2 ). With the aid of Equation (5), we estimate If H 0 ≪ δH c , the Coulomb fields produced by adjacent and more distant nodes can be ignored to a first approximation. The Coulomb contribution to the net field on a given site is small, but occasionally the redistribution of magnetic charges on nearby sites may trigger the emission of a domain wall if the net field is close to the critical value. See Sec. 4.1 for further details. In the opposite limit, H 0 ≫ δH c , these internal fields must be taken into account. The reversal of magnetization on one link alters the magnetic charges on its ends. The resulting increments of the total magnetic field at nearby nodes, of order H 0 , may be sufficient to trigger the emission of domain walls from them. Samples we studied previously [17,28] appear to be in the regime where H 0 and δH c are comparable.

Microscopic basis for the model
To test the basic model of magnetization dynamics presented in Sec. 2, we performed numerical simulations of magnetization dynamics in a small portion of the honeycomb network by using the micromagnetic simulator OOMMF [31]. The typical numerical experiment involved a junction of three permalloy magnetic wires of length l = 500 nm, width w = 110 nm, and thickness t = 23 nm [17]. We used the two-dimensional version of the oommf code with cells 2 nm × 2 nm × 23 nm. (The lateral size of the unit cell should not exceed the minimal length in the micromagnetic problem, the magnetic exchange length obtained from exchange and dipolar couplings. In permalloy, it is about 5 nm [29].) The magnetization field M(r) was allowed to relax to an equilibrium state with magnetic charge q = ±1 at the junction, Figure 3. An external magnetic field was then applied in a fixed direction and its magnitude was slowly increased keeping the system in a state of local equilibrium. Eventually, the magnet reached a point of instability when a domain wall was emitted from either the central junction or one of the peripheral ends of the wires, depending on the direction of the applied field. The wall then propagated to the opposite end of the link reversing the link's magnetization. Using those orientations of the field for which a domain wall is emitted from the junction, we determined the dependence of the critical field H on the angle θ between the field and the link in which the reversal occurs, Figure 4.
Two features of the angular dependence in Figure 4 stand out. First, H(θ) is not an even function of the angle θ, and contrary to our expectations, the critical field is not at its lowest when the field is parallel to the link. Second, the critical fields for three different links in the experiment have the same shape but differ in the overall scale H c .
We have traced the physical origin for the asymmetric dependence of the critical field H(θ) to an asymmetric distribution of magnetization at the junction, Figure 3. The energetics of the emission process shown in the figure can be described in the language of collective coordinates [32]. The soft mode associated with the emission of   a domain wall into the vertical link is the domain wall displacement X along the link.
To the first order in the applied field H and to the second order in X, the energy is where Q xx , Q xy and k are phenomenological constants. Generally speaking, the offdiagonal component Q xy does not vanish unless the magnetization distribution is symmetric under the reflection y → −y. The equilibrium position of the wall depends on the direction of the applied field H = (H cos θ, H sin θ, 0) as follows: where the offset angle α and effective chargeQ are defined through According to Equation (10), the relevant component of the magnetic field H is found by projecting the field onto the easy axis of a (majority) link, which is rotated through angle α toward the minority link. These considerations suggest the following modification for the postulated field dependence of the critical field (6): As

Numerical simulations
The heuristic considerations of Section 2 and the micromagnetic simulations of Section 3 suggested a coarse-grained model of magnetization dynamics in which the basic degrees of freedom are Ising variables of magnetization σ ij on links and magnetic charges q i on sites of the honeycomb lattice. Each link has its own fixed critical field H c . The critical fields form a Gaussian distribution (7) of width δH c around the meanH c . The average,H c = 50 mT, was chosen on the basis of our micromagnetic simulations, whereas the relative width was set to δH c /H c = 0.05, a value inspired by our experimental observations [28]. Simulations were performed in a rectangular sample with 937 links. The edge consisted of "dangling" links with no other links attached to their external ends. We choose the initial state with a maximum total magnetization that can be obtained by placing the system in a strong magnetic field along one set of links, Figure 5(a). Simulation details are described in Appendix A. Following initialization, the external field is switched off and reapplied along a different direction, at an angle θ to its initial orientation, Figure 5(b-d). To stimulate magnetization dynamics, the rotation angle must be large enough so that H would have a negative projection onto at least some of the majority links. When |θ| is between 30 • and 90 • , only one of the three sublattices of links will reverse. Two sublattices reverse when |θ| is between 90 • and 150 • . The entire lattice undergoes a reversal when |θ| > 150 • .
Aside from the number of active sublattices, there are marked differences in the dynamics of the reversal. For small angles of rotation, |θ| < 131 • , the reversals occur in a gradual and uncorrelated manner, with individual links switching when the applied field reaches the link's critical field. For larger angles, |θ| > 131 • , we observed avalanches in which long chains of links reverse magnetization simultaneously. This kind of switching happens when the sublattice whose magnetization is most antiparallel to the applied field cannot switch first because it consists entirely of minority links and must wait for one of the other sublattices to begin its reversal. If that happens in a higher field, the former sublattice acts like a loaded spring, making the reversal nearly instantaneous. A diagram depicting different regimes as a function of the field rotation angle θ is shown in Figure 6.
In the simplest case, the reversal of magnetization in a link occurs when the magnetic field reaches the critical value for that link. The links thus reverse on an individual basis, largely independently from the others (but see below). To be more precise, the two ends of a link have different critical fields and the reversal begins from the end with the lower critical field and stops at the other end. The effective probability density of the critical fields thus changes from a Gaussian distribution to It can be seen in Figure 6 (right panel) that the resulting distribution is very close to a Gaussian with renormalized mean and width, In our simulation, the renormalized values areH ′ c = 48.7 mT and δH ′ c /H ′ c = 0.042.

30 • < θ < 131 • : gradual reversal
With the field rotated through θ = 120 • , two sets of links have a negative projection of magnetization onto the field. In Figure 5(b), they are the horizontal minority links and the majority links parallel to the field. Because emission of a domain wall into a minority link requires a very high field, it is the majority links that undergo magnetization reversal first. The field makes an angle α ≈ 19 • with their easy axes, so the reversal is expected to occur around the field H 1 =H ′ c / cos (−19 • ) = 51.5 mT. Magnetization reversal in the links parallel to the field alters the magnetic charges on all sites, Figure 5 magnetization measured along the applied field is expected to be a superposition of error functions: The three terms reflect the contributions of the three groups of links with different orientations. The simulated dependence M (H) is shown in Figure 7 along with the theoretical curve (15) that takes into account the renormalization of the Gaussian distribution parameters (14).
A close inspection of the simulated curve M (H) shows that on occasion several adjacent links reverse simultaneously due to a positive feedback during the reversal. When magnetization of a link is reversed, magnetic charges at its ends are switched. The net magnetic field on an adjacent site, projected onto its easy axis, increases by The extra field is not negligible on the scale of the critical-field distribution width δH ′ c = 2.0 mT. It can help to stimulate the emission of a domain wall at an adjacent site if that site's critical field is not too high. This kind of positive feedback causes avalanches, in which magnetization reversals occur nearly simultaneously in links residing along a one-dimensional path determined by the orientations of easy axes. For example, an avalanche occurring in the background of a fully magnetized state of Figure 5(b) would travel along the vertical direction. In the limit of small feedback, ∆H ≪ δH ′ c , the distribution of avalanche lengths is exponential. Indeed, if the link starting an avalanche of length n has a critical field H, n − 1 of its neighbors must have critical fields in the range between H and H + ∆H. The probability to find such a collection of links is for a Gaussian distribution of critical fields (7). The distribution of avalanches seen in the simulation is shown in the inset of Figure 7 along with the theoretical distribution (17). These results can be directly compared to the experimental reversal curve measured in the same geometry [28], Figure 7 (right panel). Although the overall scale of the magnetic field is substantially lower, the data are well fit by Eq. (15) with H 1 = 35.9 mT and H 2 = 45.9. The ratio of the reversal fields, H 2 /H 1 = 1.28, agrees well with the theoretical value H 2 /H 1 = cos (−19 • )/ cos 41 • = 1.25. The relative widths are δH 1 /H 1 = 0.037 and δH 2 /H 2 = 0.046.
Overall, it appears that our model provides a reasonably good description of magnetization reversal when the field is reapplied at θ = 120 • to the direction of initial magnetization. In this regime, the reversal proceeds in two well-defined stages, each involving one subset of links. During each stage, links reverse largely independently, although sometimes the reversal in one link changes the field on a nearby site and triggers magnetization reversal there. The reversal fields are given approximately by the equations The reversal follows the two-stage scenario as long as For larger field rotation angle θ, the reversal proceeds in a very different manner.

131 • < θ < 180 • : reversal with avalanches
When the field is rotated through θ = 170 • relative to the direction of magnetization, the theory described in Section 4.1 no longer applies. Because H 2 , the reversal field of horizontal links, is lower than H 1 , these links should reverse first. However, that is impossible because in the initial, fully magnetized state, Figure 8(a), these are minority links whose critical field is roughly 3H c (Section 2.2), i.e., much higher than H 2 = H c / cos 11 • ≈ H c . For this reason, a horizontal link does not reverse until one of its neighbors, a majority link, reverses and in the process alters the charge at one of the horizontal link's ends. This converts the horizontal link into a majority link enabling it to reverse magnetization. It turns out that this mode of reversal is accompanied by long magnetic avalanches.
In the simplest scenario, the dynamics begins with the reversal of the weakest link with the critical field near H 1 , Figure 8(b). The reversal turns the horizontal link next to it into a majority link, which is now ready to reverse since the applied field exceeds its critical field: H ≈ H 1 > H 2 . A q = −2 domain wall emitted from its left end travels to the right end where it encounters a site with charge −1, Figure 8(b). As discussed in Section 2.3, the arriving domain wall induces the emission of another domain wall into an adjacent link, Figure 2(b). The magnetization of that link gets reversed, bringing us to the state shown in Figure 8(d). The cycle repeats creating an avalanche. In effect, we have a q = +2 charge moving along a zigzag path parallel to the applied field and reversing magnetization of the links along the way. The process continues until the moving charge reaches the edge of the system so that an avalanche extends from edge to edge. A different scenario may take place if the system has "weak" links that trigger the reversal when the applied field is at or below H 2 . These can be links at the edge of the system or some sort of defects. Their reversal converts one of the horizontal links (critical field H c1 ) to the majority status, as shown in Figure 8(b). When the applied field reaches a value sufficient to induce the reversal of that link, an adjacent link is also reversed as described above, Figure 8(c-d). The next horizontal link down the line (critical field H c2 ) will switch immediately if H c1 > H c2 . The switching will continue until the avalanche comes to a stubborn link whose critical field exceeds H c1 . Its reversal will happen in a higher applied field, possibly triggering another avalanche.
If the first reversal occurs in a link whose critical field H c1 is at the lower end of the critical field distribution, the first avalanche will be short because it is unlikely that a large number of subsequent links will have even lower critical fields. As further avalanches get terminated at links with higher critical fields, their lengths will tend to increase. Toward the end of the reversal, avalanches will begin with links whose critical fields are near the higher end of the distribution. These avalanches will be particularly long. The last avalanche in a given string of links will terminate at the edge or will meet an avalanche traveling in the opposite direction. These qualitative considerations anticipate a wide distribution of avalanche lengths. Indeed, we show in Appendix B that the avalanches have a power-law distribution of lengths, Remarkably, this result applies to any distribution of critical fields, not just a Gaussian one, and numerical simulations confirm this picture. As can be seen in Figure 9, magnetization reversal begins in an applied field H ≈ H 2 − δH 2 = 47 mT, where H 2 is given by Eq. (18). At that point, the reversals include single pairs of links from two sublattices. Long avalanches, involving as many as n = 10 and more links, are observed by the time the applied field reaches H ≈ H 2 + δH 2 = 51 mT. The length distribution is well fit by the power law (19) as can be seen in the inset of Figure 9.
The third sublattice reverses in much higher fields, H ≈ H 3 = 77 mT, where This stage of the reversal proceeds in the gradual manner described previously.

Discussion
The dynamics of magnetization in artificial spin ice is a complex problem. In this paper, we have presented a simple model for this system in terms of coarse-grained physical variables (Figure 1), Ising spins σ ij living on the links of the spin-ice lattice and magnetic charges q i residing on its sites. Inspired by our earlier studies of magnetic nanowires [33,32], where magnetization reversal is mediated by the propagation of domain walls, we have expressed the magnetization dynamics in spin ice in similar terms. Magnetization reversal in individual links of the lattice proceeds through the emission, propagation, and absorption of domain walls with magnetic charge q w = ±2. Coulomb-like interactions between the magnetic charges of the walls and lattice sites play a major role in the dynamics. For example, the magnitude of the critical field, required for the emission of a domain wall, is set by the strength of magnetostatic attraction between a domain wall and the magnetic charge of the lattice site. These heuristic considerations have been confirmed and refined through micromagnetic simulations of a small portion of the spin-ice lattice containing a few links. Quenched disorder is another major element affecting the magnetization dynamics. Small imperfections of the artificial lattice are expected to produce a Gaussian distribution of critical fields. The experimentally measured curve [28] is consistent with a Gaussian shape and width δH c /H c ≈ 0.05.
The dynamics of magnetization reversal strongly depends on the direction of the external magnetic field. If the field is applied at a small angle relative to the magnetization of a (fully magnetized) sample, θ < 131 • for the parameters we used, the reversal proceeds in a gradual way, with links reversing more or less independently of each other, when the strength of the applied field exceeds the threshold of a given link. For larger angles of rotation, the reversal proceeds in one-dimensional avalanches that can easily span the entire length of the system. The reversal in one link with a critical field H triggers the reversal in several others along the chain. The avalanche stops when it encounters a link whose critical field exceeds H. In this regime, avalanche lengths are distributed as a power law, P n = C/n.
It should be pointed out that we model the magnetization dynamics in artificial spin ice as a purely dissipative process, in which the system moves strictly downhill in the energy landscape. Such a picture is very different from an earlier approach extending the notion of an effective temperature to these far-from-equilibrium systems [20,34]. Whereas energy of a microstate plays a major role in the effective thermal approach, our method puts the focus on energy gradients, or forces between magnetic charges.
This study has a limited scope. We focus on a continuously-connected honeycomb network realized in several experimental studies [18,16,17] and cover only the basic regimes of its magnetization dynamics, Figure 6. Interesting phenomena arise at the boundaries between different regimes, particularly when the field is completely reversed, θ = 180 • . In this case, avalanches lose their unidirectional character and become random walks. As the magnetization reversal proceeds, avalanches can begin to intersect and block one another.
Our method can be easily extended to connected networks with other geometries such as square spin ice [15]. Budrikis, Politi and Stamps used a similar heuristic approach to study the dynamics of disconnected magnetic islands [35].
program checks whether the net field has a negative projection H e = H cos (θ − α) onto the link's easy axis, Eq. (12). If H e < 0, the program calculates the weakness of the site and link, W = |H e | − H c . The site and link with the largest W in the sample are considered to be the weakest. As the applied field increases, the largest W becomes positive, triggering the emission of a domain wall from the weakest site into the weakest link. The domain wall propagates to the other end of the link where it is absorbed, either immediately or after the emission of another domain wall as described in Section 2. Once the reversal process that started with the weakest site is complete, the program looks for the next weakest site. The process is repeated until there are no positive W in the system. Spin ice rules are satisfied at each site at all times. No thermal effects are considered.

Appendix B. Statistics of avalanches in the presence of a weak link
Here we derive the statistics of avalanches discussed in Section 4.2. In this case, the reversal begins in Link 1 (critical field H) and spreads to consecutive Links 2, 3, . . . , n (of the same sublattice) as long as their critical fields are lower than H. The avalanche stops when it encounters link n + 1 whose critical field exceeds H. The probability density of the critical-field distribution is ρ(H) and the cumulative probability distribution is Consider an avalanche beginning on link k with a critical field between H and H + dH. The k − 1 preceding links must have critical fields less than H. If the avalanche has length n then links k + 1, k + 2, . . . , k + n − 1 must have critical fields less than H, whereas link k + n must have a higher critical field. The probability of such a distribution is The probability to find an avalanche of length n is found by summing this distribution over the initial position of the avalanche k and integrating over the critical field H.
Performing the sum first, we find Note that F n is an expectation number of avalanches, not a probability distribution normalized to 1. Observing an avalanche of length n does not exclude the possibility of observing an avalanche of a different length n ′ in the same chain during the same reversal process. The probability that an avalanche will have length n is P n = F n / L n=1 F n ∼ 1/(n ln L) (B.6) for large L.