Chiral switching and dynamic barrier reductions in artificial square ice

Collective dynamics in lithographically-defined artificial spin ices offer profound insights into emergent correlations and phase transitions of geometrically-frustrated Ising spin systems. Their temporal and spatial evolution are often simulated using kinetic Monte Carlo simulations, which rely on the precise knowledge of the switching barriers to obtain predictive results in agreement with experimental observations. In many cases, however, the barriers are derived from simplified assumptions only, and do not take into account the full physical picture of nanomagnetic switching. Here we describe how the immediate magnetic environment of a nanomagnet reversing via quasi-coherent rotation can induce clockwise and counter-clockwise switching channels with different barrier energies. This barrier splitting for chiral reversal channels can be sizeable and, as string-method micromagnetic simulations show, is relevant for artificial spin ice systems made of both exchange -- as well as magnetostatically --dominated units. Due to the barrier splitting (and further reductions due to non-uniform reversal) transition rates can be exponentially enhanced by several orders of magnitude compared to mean-field predictions, especially in the limit of rare switching events where thermal excitation is less likely. This leads to significantly faster relaxation time scales and modified spatial correlations. Our findings are thus of integral importance to achieve realistic kinetic Monte Carlo simulations of emergent correlations in artificial spin systems, magnonic crystals, or the evolution of nanomagnetic logic circuits.

Collective dynamics in lithographically-defined artificial spin ices offer profound insights into emergent correlations and phase transitions of geometrically-frustrated Ising spin systems. The understanding of experimentally-observed temporal evolution of extended spin ices are often guided and supported by model predictions, for example from kinetic Monte-Carlo simulations. This coarsegrained approach, which disregards microscopic details of the moment reversal, allows to simulate systems with a large number of moments evolving over long time scales, which otherwise would be too computationally-costly to be implemented in full micromagnetic simulations. To obtain correct relaxation time scales and spatial correlations, kinetic Monte Carlo simulations rely instead on the precise knowledge of the rates for individual moment reversal. These rates are determined by the switching barriers which, in many cases, are derived from simplified or approximative assumptions only, which do not take into account the full physical picture of nanomagnetic switching.
In this work, we describe how the immediate magnetic environment of a nanomagnet reversing via quasi-coherent rotation can induce clockwise and counter-clockwise switching channels with different barrier energy. We compare predictions from a perturbative model to switching barriers obtained from micromagnetic string-method simulations for two different -exchange-vs. magnetostaticallydominated -artificial square ice geometries. Taking into account the spatial extension and nonuniform magnetic behaviour, we find further reductions and enhanced barrier splitting, especially in the case of magnetostatically-dominated nanomagnets. These modifications of the switching barriers lead to exponentially enhanced relaxation kinetics, especially in the limit of rare events. From kinetic Monte-Carlo simulations we find that the evolution invoking split barriers yields much faster relaxation time scales and results in different spatial moment configurations compared to the often-employed mean-field transition barriers. Our results highlight how the local magnetic environment can significantly enhance the transition kinetics and affect emergent correlations, even without invoking defects or additional anisotropies. These findings are thus of integral importance to achieve realistic kinetic Monte Carlo simulations of emergent correlations in extended artificial spin systems, magnonic crystals, or the evolution of small-scale nanomagnetic logic circuits.
Artificial spin ice systems are lithographically-created lattices of elongated single-domain nanomagnets, and have been designed to investigate the effect of correlations and the onset of long-range order in frustrated twodimensional magnetic lattices [1][2][3][4].
Of particular interests is the evolution of extended spin ice lattices from a field-saturated state towards an energetically favourable (ground) state, driven by thermallyactivated reversal of individual nanomagnets. Such experiments have been performed mainly using photoelectron emission microscopy, and gave valuable insight on the relaxation process and the formation of spatial correlations [5][6][7][8][9]. These results are often compared to model predictions for the temporal and spatial evolution from kinetic Monte Carlo (kMC) simulations. Their major advantage over full micromagnetic simulations is that they are less computational costly, and thus can be extended FIG. 1. Switching barriers in artificial square ice. (a) Nanomagnets are arranged on a square lattice with periodicity a. The strongest mutual interaction JNN acts between perpendicular nearest-neighbour moments mj. (b) Fully-magnetised double-vertex state #0, with enumeration of moments. The central nanomagnet (black) can rotate either clockwise (red arrow) or counter-clockwise (blue arrow) from left (←) to right (→). (c) Energetically-favourable states, on the left, feature more head-to-tail configurations between the central and neighbouring moments (green arrows). (d) In a pertubative picture, the switching barrier energies can be obtained by adding the interaction energies to the switching barrier ∆E sb of an isolated nanomagnet. If only the energies of the equilibrium configurations (←, →) are taken into account, a mean-field barrier is obtained (gray arrow to cross). In case the environment features a perpendicular magnetisation, i.e.M i,⊥ = 0, the high-energy states (↑, ↓) will split, and thus yields separate transition barriers for clockwise (red) and counter-clockwise (blue) rotation.
sistently overestimates the barriers and underestimates the chiral splitting of the former, and are not applicable even with renormalised parameters in the case where non-coherent reversal modes are possible.
Reductions and splitting of the switching barriers lead to exponentially enhanced transition rates especially in the limit of rare events, as we show with a modified Arrhenius law. Using the rates for the chiral transition channels as input for kMC simulations, we find that the evolution of an extended square ice proceeds much faster, and involves different spatial correlations when compared to a mean-field model. The influence of the immediate environment on the nanomagnetic switching thus is a key ingredient to correctly model the relaxation dynamics of artificial spin ices, as well as of functional magnonic materials and small-scale circuits for computation. We therefore expect our results to be relevant to different communities making use of thermally-driven relaxation of interacting nanomagnets.
This work is structured into three sections: In Sec. I a basic understanding is derived on how the magnetic environment can lead to chiral switching channels in artificial square ice. Sec. II compares point-dipole model predictions to micromagnetic simulations of nanomagnets of different dimensions and material parameters. In Sec. III, ramifications of the modified switching barriers on the switching rates of single nanomagnets and relaxation kinetics of extended artificial square ice are discussed.

I. CHIRAL MOMENT REVERSAL
Artificial square ice is a lithographically-designed magnetic metamaterial with identical nanomagnets arranged on a square lattice with periodicity a, see Fig. 1(a).
Due to their shape anisotropy, each stadium-shaped nanomagnet with length l, width w, and thickness t is quasi-uniformly magnetised, and thus behaves like an Ising macrospin. Without applied field or interactions with neighbours, the magnetic moment will align with the long axis, i.e. to the left (←) or right (→). To spontaneous switch between these energetically-degenerate configurations (without loss of generality from ← to →), the moment rotates to overcome a metastable state for which the net moment points along the nanomagnets' short axis, i.e. ↑ or ↓. The difference between the metastable and equilibrium micromagnetic energies E mm gives the single-island switching barrier ∆E sb : The value of ∆E sb depends on the size, shape and material of the individual elements [16][17][18][19]. For the assumption of uniform magnetisation, the barrier is given by the shape anisotropy ∆E shape sb = K shape V , with V being the volume of the nanomagnet. Values of K shape are either tabulated [20] or can be calcuated via magnetostatic simulations [21,22].
The mutual coupling between nanomagnets is given by magnetostatic interactions, which takes the following form between point-dipole moments m i and m j separated by a distance vector r ij = r i − r j : In artificial square ice the strongest coupling, denoted by J NN , occurs between nearest-neighbour nanomagnets at a 90°angle, see Fig. 1(a). Using the lattice periodicity a and the net moment m = |m| = M sat V of a nanomagnet with volume V and material saturation magnetisation M sat , we define a convenient energy scale J dip NN : Due to the pronounced distance dependence, E dip ∝ r −3 , the coupling between nanomagnets meeting at the vertex points highlighted in Fig. 1(a), is dominant over further-range interactions [23,24]. Therefore, we investigate the switching barriers for moment reversal of a central nanomagnet under the influence of its closest neighbours only.

A. Switching environments
For an infinite artificial square ice the environment that influences the reversal of a nanomagnet forms a doublevertex configuration, as depicted in Fig. 1(b). Here, each tip of the central nanomagnet (black) is in close interaction to three other nanomagnet (gray), whose magnetisation remains largely unchanged during the reversal of the central nanomagnet. The extended square lattice then can be considered as an infinite tiling of these motifs.
We denote each magnetic equilibrium configuration with a state number #i determined by the arrangement of the surrounding nanomagnets, and the orientation of the central (switching) nanomagnet, which can point to the left (←) or to the right (→). The state number #i can be obtained from the binary representation of the relative configuration of the six surrounding nanomagnet numbered 1 to 6 according to the scheme shown in Fig. 1(b):

B. Switching barriers from a point-dipole model
In the following, we discuss a perturbative approach to switching barriers in artificial square ice. Here, the single-nanomagnet barrier ∆E sb , Fig. 1(d), is modified due to energy contributions arising from the interactions with the immediate magnetic environment.
To discuss a specific example, we focus on the fullymagnetised double-vertex environment #0 depicted in Fig. 1(b). The energy of the configuration where the central moment points (exactly) to the left, ←, is lower than when it points to the right, →, where three magnetic charges meet at each vertex point. Using simplified assumptions and symmetry arguments, one can derive a mean-field switching barrier ∆E i dip ← to → (as the average of the barriers for clockwise and counter-clockwise rotation, see Appendix A in Ref. [15]). Its value depends solely on the energies of the equilibrium configurations, as indicated by the gray arrow in Fig. 1 This mean-field barrier, however, is missing a crucial point, as independent relaxation pathways via clockwise and counter-clockwise rotation of the central nanomagnet need to be considered, i.e.
These barriers are not necessarily equivalent, as shown in Fig. 1(c): Due to their staggered spatial arrangement, the central moment in a fully magnetised environment #0 will preferably align ferromagnetically with its neighbours, thus forming a head-to-tail flux-closure configuration which reduces the dipolar energy term in Eq. (2). Therefore, transitions via counter-clockwise rotations (blue) of the central nanomagnet will be largely favoured over those via clockwise rotations (red). From the dipolar energies for all environmental states and orientation of the central moment, Fig. 3(a), one can obtain the respective clockwise and counter-clockwise switching barriers, Fig. 3(b). Under the assumption that ↑ and ↓ align perfectly along the short axis of the nanomagnets, the energy splitting between the configurations is symmetric around the mean-field barrier ∆E i dip ← to → (marked by crosses), and equals the three distinct values shown in Fig. 3(c).
Barrier splitting occurs for all environments which feature a perpendicular effective fieldM i,⊥ generated by the neighbouring nanomagnets, that acts on the central nanomagnet [the indices b j are as defined in Eq. (4)]: The normalised perpendicular magnetisationM i,⊥ can take the values of 0 (black and purple in Figs. 2 and 3), ±2 (orange), and ±4 (red and blue) only. Thus, we can modify Eq. (5) to include an additional energy term: For a central moment initially pointing to the left (←), the second term (derived in Appendix A) is subtracted for clockwise, and added for counter-clockwise reversal.
In conclusion, with a simplified point-dipole model the switching barriers are modified by the choice of clockwise vs. counter-clockwise rotation, if the moment interacts with an effective perpendicular stray field generated by its environment. As shown in Fig. 2, 40 out of the 64 double-vertex configurations feature a finite perpendicular magnetisationM i,⊥ = 0 acting on the central nanomagnet. In particular, we expect a maximum chiral barrier splitting for the fully-magnetised environments (marked in red), which are commonly-used initial states for thermal relaxation studies of artificial spin ice [7].
These observations are by no means a curiosity, since nanomagnetic switching will occur predominantly via the more favourable pathway. We thus expect profound consequences on the switching rates and transition kinetics when taking into account the chiral splitting.

II. MICROMAGNETIC SWITCHING
The point-dipole switching barriers represent a perturbative approach parametrised by two parameters only: First, ∆E sb describes the switching barrier of an isolated nanomagnet, and implicitly depends on its shape and size [20]. Second, J dip NN quantifies the energies of equilibrium configurations due to the interactions between nanomagnets placed on the square lattice, and modify the switching barrier of individual nanomagnets. The mean-field barrier energy in Eq. (9), however, does not take into account possible non-coherent moment reversal. In particular, it does not describe the influence of material parameters such as the saturation magnetisation M sat and the exchange strength A ex , thermal fluctuations, and the magnetic environment. Due to these effects the net moment can be dynamically reduced during reversal, e.g. via non-uniform buckling modes, vortex creation, or domain formation [25][26][27].
To have a nuanced look on how the barrier energy depends on (1) the material parameters, (2) the nanomagnet shape and size, and (3) the interactions with neighbouring nanomagnets, we now turn to a full micromagnetic simulation of the reversal barriers.
A. Implementation of string-method simulations Contrary to simulations employing the Landau-Lifschitz-Gilbert equations to explicitly solve the timedependent evolution of the nanomagnetic reversal, we determined the associated barrier energy using a timeindependent string method. In this approach, starting from a coherent moment reversal, the moment configurations are iteratively optimised to a minimum-energy path through configuration space, and thus yield the lowest barrier energy associated to that reversal process. As in our previous work [15], which also discusses further simulation details, we implement here the simplified and improved string method (SISM) [28] using the finite-element micromagnetic code magnum.fe [29].
We consider two artificial square ice geometries with distinct choices for M sat and A ex , representing different regimes. Meshes discretizing the considered geometries, i.e. an individual nanomagnet and the double-vertex configurations, were created with the software gmsh [30].
First, we consider a geometry largely dominated by exchange interactions, which favour a coherent reversal, and thus may resemble the macrospin model derived in Sec. I: Nanomagnets with dimensions l × w × t = 150 nm × 100 nm × 3 nm are placed on a square lattice with periodicity a = 240 nm. The material parameters M sat = 790 kA/m and K = 0 correspond to bulk permalloy (Fe 0.2 Ni 0.8 ) values at 300 K. The exchange stiffness A ex = 13 pJ/m was obtained from a temperaturedependent scaling [31][32][33][34] where A ex (0) = 950 kA/m and M sat (0) = 18 pJ/m denote the respective permalloy bulk parameters at 0 K. Second, we consider a system for which we expect sizable magnetostatic effects leading to non-uniform magnetic configurations during reversal: Nanomagnets with dimensions l×w×t = 470 nm × 170 nm × 3 nm are placed on a square lattice with periodicity a = 600 nm. This geometry, or choices close to it, have been used in several experimental studies [7,13,35,36]. The saturation magnetisation M sat = 350 kA/m corresponds to the value given in [7], and from the scaling in Eq. (10) we obtain A ex = 3.25 pJ/m. Although the saturation magnetisation is lowered significantly, the exchange stiffness is even more reduced, thus we expect the energetics dominated by magnetostatic interactions. We again assume a vanishing magnetocrystalline anisotropy, K = 0.

B. Influence of environment on reversal
We select representative environment states for each of the five classes introduced Fig. 2(b), for which we expect no (black and purple), intermediate (yellow), and high barrier splitting (red and blue), respectively. For these states, shown at the top of Fig. 4, we compare the chiral switching barriers obtained from a macrospin approximation (left) and micromagnetic string-method simulations (right) for the two different square ice geometries (schematics on the left are shown to scale). Clockwise and counter-clockwise barriers are marked in red and blue, respectively, and the difference (barrier splitting) is plotted as black bars.
The macrospin model considers variation of the single-moment barrier ∆E sb (calculated from the shape anisotropy of a uniformly-magnetised nanomagnet ∆E shape sb = K shape V ) due to point-dipole-like interactions quantified by J dip NN ∝ (M sat V ) 2 /a 3 , as derived in Sec. I B. We also plot the mean-field barrier from Eq. (5) (black horizontal line).
The chiral barriers from the string-method simulations are calculated from the energy difference between the micromagnetic net energies of the metastable barrier configuration (↑, ↓) and the initial configuration (←): For the exchange-dominated square-ice geometry, Fig. 4(a-e), the mean-field barrier always overestimates the lower of the two micromagnetic barriers (for the chosen cases corresponding to counter-clockwise reversal marked in blue), as already discussed in our previous work [15]. Compared to the point-dipole model, micromagnetic simulations consistently give lower chiral switching barriers. The difference of the barrier energies, i.e. the chiral splitting, is enhanced, however, it remains proportional to the perpendicular moment |M i,⊥ | generated by the environment. The reversal process is still governed by an almost-coherent rotation of the central moment (Appendix C 1). Therefore, one can obtain reasonable switching barriers from the perturbative decomposition of Eq. (9) by using a lower single-moment barrier ∆E string sb and stronger interactions J mm NN (Appendix C 2). For the magnetostatic-dominated square-ice geometry, Fig. 4(f-j), the differences are even more pronounced. This is because the reversal behaviour is no longer uniform, as discussed in Appendices C 1 and C 2. Interestingly, non-coherent reversal is particularly pronounced in environments with vanishing perpendicular magneti-sationM i,⊥ = 0 (black and purple). This leads to large reductions of the switching barrier, e.g. more than −50 % in the case of state #22 in Fig. 4(f), when compared to the point-dipole model. In general, the point-dipole barriers overestimate all micromagnetic barriers by a considerate margin, and underestimate the chiral barrier splitting, as is evident when comparing the black bars shown in Fig. 4(h-j). For environments with a finite perpendicular magnetisation the barrier splitting is of similar magnitude, Fig. 4(h-j), and thus does not follow the proportionality |E i,↑ − E i,↓ | ∝ |M i,⊥ |. Therefore, the sim-plified decomposition of switching barriers of Eq. (9) is no longer valid.
In general, by taking into account the micromagnetic nature of the moment reversal, we observe both a reduction of the switching barriers as well an enhanced separation of the chiral barriers for environments witĥ M i,⊥ = 0. As lower barriers are easier to overcome for thermally-induced reversal, we therefore expect significant enhancement of the kinetics of artificial square ice, which may yield different relaxation time scales and emergent correlations.

III. TRANSITION KINETICS
The effect of barrier splitting on the net transition rate can be generalised by using an average barrier and a factor f , which describes the symmetric splitting of the clockwise and counter-clockwise barriers around the average barrier, as depicted in Fig. 5(a), For the barriers derived from the point-dipole picture in Sec. I B, the average barrier corresponds to the meanfield barrier of Eq. (5), i.e. ∆E dip avg = ∆E i dip ← to → . The splitting factor vanishes, f = 0, for environmental states without a perpendicular magnetisationM i,⊥ = 0 acting on the switching moment. From the results of micromagnetic simulations discussed above we obtain non-zero values of f between a few percent up to about 20%.

A. Modified Arrhenius law for barrier splitting
The temperature-dependent transition rate for spontaneous switching over an average energy barrier ∆E avg can be obtained via the Arrhenius law [37,38], with the attempt frequency ν 0 and the Boltzmann constant k B = 8.62 × 10 −5 eV/K: The attempt frequency ν 0 depends on the shape, size, and material of the nanomagnets. Typical values of ν 0 are in the order of 10 9...12 Hz [38,39], and even faster time scales have been discussed [40,41]. We need to consider the clockwise and counterclockwise reversal as parallel and independent channels of relaxation, and the rates associated to each of the two barrier energies, i.e. ∆E avg (1−f ) and ∆E avg (1+f ), need to be added to obtain an effective transition rate. Therefore, for the definition of the rate in Eq. (15), we explicitly included a pre-factor of two, to account for degenerate The sum rate νsum = νmax + νmin of the two parallel relaxation channels can be significantly enhanced compared to the rate expected from overcoming an average barrier νavg, shown by the solid lines indicating νsum/νavg. The colours denote different kinetic regimes given by the ratio of the average barrier ∆Eavg compared to the thermal energy kBT , as indicated by the numbers on the right. Dashed lines denote the respective rate enhancement νmax/νavg associated with transitions via the lower barrier only, which underestimates the rate by a factor of two in the limit of f → 0. Within the shaded area νmax exceeds 90% of the net rate νsum. The rate enhancement is particularly large for rare events where the reduced energy, as indicated by the numbers on the right, is large, i.e. ∆Eavg/kBT 1. However, it hardly matters in the limit of superparamagnetic fluctuations where ∆Eavg/kBT → 1.
clockwise and counter-clockwise relaxation channels over the average barrier.
Due to the pronounced non-linearity of the Arrhenius law, the summation of rates leads to an effective increase of the net transition rate ν sum of thermally-activated switching when compared to the rate ν avg associated with the average barrier (derivation in Appendix B): The ratio of the joint rate compared to the averagebarrier rate, i.e. ν sum /ν avg , are compared in Fig. 5(b) for different splitting ratios f and reduced energies ∆E avg /k B T . The exponential rate enhancement is particularly pronounced in the limit of rare events with ∆E avg /k B T 1 where a splitting of barriers can increase the (albeit low) transition rates by several orders of magnitude (purple lines). In contrast, in the limit of superparamagnetic fluctuations, i.e. ∆E avg /k B T → 1, bar-rier splitting increases the net rate only moderately (blue lines).
For the assumption that transitions occur predominantly via the lower barrier only, we have to consider the transition associated to the smaller barrier, i.e. (1 − f )∆E avg , see dashed lines in Fig. 5 giving the ratio ν max /ν avg . For high splitting ratios f and ∆E a /(k B T ) 1 the rate ν max will approach ν sum , as transitions by the higher-lying barrier become irrelevant. The shaded area of Fig. 5 marks were ν max exceeds 90% of the value ν sum . Within this regime, transitions via the lower-lying barrier might be a good approximation of the net reversal rate. In the limit of f → 0, however, where both barriers are equal, using the rate ν max will underestimate the net rate by a factor of two, i.e. (ν max /ν sum ) f →0 = 1/2.
In the case of artificial square ice, and depending on the kinetic regime given by the relation between the energies ∆E sb , J NN , and k B T , approximating the transition rates via the mean-field barrier [7,8,13,36] or the minimum barrier [42,43] thus may significantly underestimate the speed of evolution.

B. Temporal evolution of extended square ice
To illustrate the consequences of barrier splitting, we now turn to the evolution of extended square-ice arrays. In many thermal relaxation studies of artificial spin ice, the main interest lies in the onset of phase transitions and formation of emergent correlations. In many experiments, the system evolves from a field-set fully magnetised state, for which we predict a particularly strong barrier splitting. We therefore expect that the initial demagnetisation of a magnetic-field-saturated artificial square ice array will be particularly affected by the modified transition kinetics.
To model the relaxation, kinetic Monte Carlo (kMC) simulations are often employed [44][45][46]. The kMC algorithm provides a numerical solution to the master equation, which is a system of linear differential equations describing the evolution of the probabilities for Markov processes in systems that jump from one state to another in continuous time [47]. Using this method, both the equilibrium expectation values of populations and their dynamical evolution during a thermalization process can be retrieved.
In this work, kMC simulations are performed using a custom-written code [48], with a system of 50 × 50 moments and periodic boundary conditions. The initial configuration is uniformly-magnetised, with the net magnetisation being parallel to a diagonal direction of the array. The demagnetisation due to spontaneous moment reversals is tracked over time for 125 × 10 3 kMC steps, and averaged over 20 individual simulation runs. We use the point-dipole energy barriers as shown in Fig. 3(b) with parameters J dip NN = 0.178 eV and ∆E sb = 1.327 eV (i.e., using values for the exchange-dominated square-ice ge- ometry), and calculate the environment-dependent transition rates ν avg and ν sum at a temperature of T = 300 K as input parameters for the kMC simulations. Fig. 6(a) compares the time evolution of the net magnetisation of square ice for rates from the mean-field barriers (dashed black line) to the model taking into account chiral barrier splitting (solid red line). The time is measured in multiples of the inverse attempt frequency, τ = ν −1 0 . We find that the onset of demagnetisation for the split-barrier model (red) happens two orders of magnitude earlier than for the average-barrier model (black). In the case of the average-barrier model, the demagnetisation involves bouts of rapid evolution interrupted by phases with little change, indicating avalanche-like dynamics [49][50][51][52]. In contrast, the split-barrier model shows a smooth demagnetisation, with a rate (solid thick lines indicate evolution from from 90% to 50%) which is about three orders of magnitude faster compared to the meanbarrier model.
When assessing the emergent spatial correlations, we find that the evolution for the mean-field model is governed by the propagation of a string of ground-state vertices (in blue) wrapping the system (due to the periodic boundary conditions), as shown in Fig. 6(b). The final state has a magnetisation of about 16% of its initial value. The snapshots of the spatial configuration of the splitbarrier model at M = 0.9 and M = 0.5 (with M normalised to the initial field-set magnetisation) in Fig. 6(c) appear somewhat similar to that of the mean-field case. There are more possible transitions for the system to explore, however, and the final state of the evolution corresponds to a multi-domain state with almost vanishing magnetisation, M ≈ 0.
Thus, our kMC results show that the modified hierar-chy of transition barriers due to the chiral barrier splitting may have subtle, but relevant, consequences: In certain cases, the kinetic relaxation pathways are not simply dictated by equilibrium-energy arguments. This will modify the emergence of spatial correlations, which needs to be explored in a systematic study and compared to experimental results [1,9,10,49,[53][54][55][56][57][58][59][60].

IV. CONCLUSIONS
To realistically model the temporal evolution of artificial spin ices or small-scale nanomagnetic circuits it is necessary to know the switching barriers for the singlemoment reversal. In this work, we quantified how magnetostatic interactions with neighbouring nanomagnets modify the switching barriers in artificial square ice in absence of extrinsic effects such as defects or spurious fields. We found that for environments which feature a finite perpendicular magnetic field acting on the switching nanomagnet clockwise and counter-clockwise moment reversals need to be considered independently. The resulting barrier splitting can be sizeable. In the case of exchange-dominated nanomagnets supporting coherent rotation modes the splitting can be predicted from a modified point-dipole model. Taking into account the finite size of the nanomagnets and the influence of material parameters, further barrier reductions were obtained from micromagnetic simulations. These reductions are particularly strong for magnetostaticallydominated nanomagnets embedded in environments that do not promote reversal via a distinct chiral switching channel.
The splitting and reduction of transition barriers ex-ponentially increase the transition rates when compared to a mean-field average barrier. Depending on the dynamical regime, which depends on the relationship between the average barrier energy, barrier splitting and temperature, we found that transition rates are especially enhanced in the limit of rare events. We modelled the evolution of extended artificial square ice with kinetic Monte Carlo simulations, and compared a meanfield model with the model that takes into account the barrier splitting. We found that the onset and speed of evolution is largely enhanced in the latter case. Furthermore, while mean-field barriers are solely dictated by equilibrium-energy arguments, the chiral switching barriers depend on the kinetics of reversal. Thus, more and different relaxation pathways are accessible, which modifies the emergent spatial correlations and routes towards the ground state.
Our results are a step towards a deeper understanding of the single-moment switching of nanomagnetic systems, highlighting how faster time scales of relaxation can be caused via intrinsic interactions with the magnetic environment. These findings are relevant to the field of artificial spin systems, and can be extended from artificial square ice to other moment configurations, such as kagome ice [6,49,50,61], and square-ice-like tetris, shakti, and brickwork lattices featuring asymmetric moment coordinations [62][63][64]. We also expect that these concepts are relevant for the utilisation of magnetic metamaterials for magnonics [65][66][67][68][69][70][71] and nanomagnetic computation [13,36,43,72,73]. We assume that all moments are strictly parallel (←, →) to the nanomagnet long axis, and remain static during the reversal of the central moment. We furthermore assume that the configuration of highest energy corre-

FIG. 7.
Definition of interactions JNN and J * NN . Nearest-neighbour dipolar interactions for (a) the equilibrium and (b) the high-energy configuration. In the latter, the staggered arrangement of moments gives rise to a preferred ferromagnetic head-to-tail arrangement in the excited configuration, with a modified dipolar nearest-neighbour interaction strength J * NN = JNN/3. Interactions of the central moment with horizontal moments via J * = 0 vanishes. spond to those with perpendicular central moment (↑, ↓). These approximations are approximately valid for weak interactions only. In general, however, they are a gross oversimplification, as due to the pairwise couplings the macrospins may rotate away from the local symmetry axis. This would result e.g. non-symmetric splitting for clockwise and counter-clockwise transitions (i.e. ∆φ = π). Nevertheless, the strict limitation of moment direction allows to employ the anti-symmetry of the dipolar interaction energy under moment rotations of π, The mean-field barrier ∆E i dip ← to → is the average of clockwise and counter-clockwise energy barriers. Under the above assumptions, and as derived in Appendix B of Ref. [15], it is determined by ∆E sb and the energy difference between the equilibrium states before and after switching. It does not depend, however, on the energies of the intermediate high-energy configuration.
The barrier splitting ∆E dip split between clockwise and counter-clockwise rotation, see Eq. (9), can be calculated from the energy difference of the high-energy states. Using the anti-symmetry argument of Eq. (A1), i.e.
Here, E dip, * i,↑ takes into account the dipolar interaction terms with the central moment only (couplings between other moments remain unchanged by the reversal, and thus fall out of the energy difference).
As shown in Fig. 7(b), in the high-energy state interactions with the horizontal nanomagnets, i.e. moments 5 and 6 in Fig. 1(b), will vanish as J * = 0. Due to the staggered arrangement of the moments, a ferromagnetic head-to-tail alignment of the central nanomagnet to its perpendicular neighbours is favourable, whereas the opposite orientation is penalised, see Fig. 1(c). Thus, the net perpendicular magnetisationM i,⊥ = Σ 4 j=1 (−1) bj is relevant to the splitting only. With a modified pair-wise interaction J dip, * Appendix B: Arrhenius law for barrier splitting If the transition barriers split symmetrically by a fraction f around the average barrier ∆E avg to values ∆E cw/ccw = (1 ± f )∆E avg , the joint effective rate ν sum = ν cw + ν ccw can be expressed as follows: Here, C = ∆E avg /(k B T ) denotes the reduced average switching barrier energy. We assume that the attempt frequencies ν 0 are independent of the energy of the saddle point, i.e. ν cw 0 = ν ccw 0 = ν 0 . The transition rate ν avg (∆E avg , T ) is defined in Eq. (15).
The maximum of the clockwise and counter-clockwise switching rates is associated to the lower-lying barrier energy (1 − f )∆E avg . In the limit of f → 0, ν max is a factor of two smaller than the rate ν sum (f = 0), and approaches the value of ν sum for large splitting f or large reduced energy C 1: Appendix C: Additional simulation results

Magnetisation during reversal
To quantify the uniformity of the magnetic reversal, Fig. 8 shows the averaged moments for a single (noninteracting) nanomagnet and each of the representative double-vertex configurations presented in Fig. 4. Here, the average magnetisation of each nanomagnet is plotted for every step of the string-method minimum-energy path. The horizontal coordinate roughly corresponds to the rotation angle φ of the central moment, with end points denoting initial and final equilibrium states.
In general, we find that in equilibrium the net moment |M| (magenta lines) is very close to one. Therefore, the static nanomagnets assume an almost saturated configuration, with the magnetisation largely aligned with the long axis of the nanomagnet and limited edge bending (M x , red lines). Interactions with neighbouring moments, however, can induce sizeable perpendicular moment contributions (M x ⊥ M y , blue lines) in environ- ments that feature a finite perpendicular magnetisation M i,⊥ = 0, i.e. configurations #0, #2, and #16.
In the case of the geometry with small islands dominated by exchange interactions the magnitude |M | remains largely constant, Fig. 8(a). The reversal thus represents a quasi-uniform rotation of the central moment. During reversal, the magnetisation components of the neighbouring nanomagnet can vary, allowing the system to evolve via the most efficient pathway.
For the geometry with large islands dominated by magnetostatic energy, shown in Fig. 8(b), the switching of the non-interacting nanomagnet (left) involves a reduction of the net moment to 91%, and thus does not conform to a uniform moment rotation. For environment states #0, #16, and #2 withM i,⊥ = 0 the reduction of net magnetisation is similar to that of an individual nanomagnet. In contrast, for environment states #6 and #22, witĥ M i,⊥ = 0, we observe a pronounced reduction of magnetisation in the high-energy configuration to less than 70% of the net moment, indicating non-coherent reversal. This leads to a reduction of the switching barrier as well, as discussed in Sec. II. The magnetic configuration of neighbouring nanomagnets varies less during reversal when compared to the small-island geometry. This is because the relative volume fraction of the large magnets meeting at the vertex point, where the spin structure varies the most, is smaller.

Micromagnetic barriers
The switching barriers obtained from micromagnetic string-method simulations for the environment states #0-#31 are summarised in Fig. 9 with (a-c) showing the results for the exchange-dominated, and (d-f) the magnetostatic-dominated geometry.
We compare the clockwise and counter-clockwise barriers, large and small circles in Fig. 9(a,d), to a modified mean-field model. The predictions are based on Eqs. (5) and (9), but instead of energies derived from point-dipole calculations we use those obtained from micromagnetic simulations, as follows: First, the switching barrier ∆E sb of an isolated nanomagnet simulated with the string-method is used, as opposed to the shape anisotropy calculated for a uniformlymagnetised nanomagnet. For the exchange-dominated nanomagnet geometry we obtain ∆E shape sb = 1.540 eV compared to ∆E string sb = 1.327 eV, which is a reduction of −14%. For the magnetostatic-dominated nanomagnet geometry, the barrier reduction is even bigger, with ∆E shape sb = 2.153 eV and ∆E string sb = 1.691 eV (−22%), as the reversal is no longer coherent [see Fig. 8(b)].
Second, E mm i,← and E mm i,→ correspond to the micromagnetic equilibrium energies of the static configurations. Together with ∆E string sb , a modified mean-field barrier can be calculated from Eq. (5) [crosses in Figs. 9(a,d)].
Third, to estimate the barrier splitting, the nearestneighbour interaction J mm NN is rescaled. The re-scaling is motivated by the point-dipole model, which predicts an energy difference of (16 − 24 √ 5 125 )J dip NN between the lowestlying ground state and highest monopole state: (C1) For the exchange-dominated small-island geometry we find that the modified mean-field barrier gives a passable estimate for the average switching barrier, as the small differences in Fig. 9(c) show. The chiral barrier splitting is well-described by the re-scaled energy J mm NN , albeit with a small reduction for fully-magnetised environments marked red and blue in Fig. 9(c).
For the magnetostatic-dominated reversal in the largeisland geometry, Fig. 9(d-f), the mean-field approach fails: The barrier splitting is both overestimated for environments with |M i,⊥ | = 4 (red and blue) as well as underestimated in the case of |M i,⊥ | = 2 (yellow). In particular, the mean-field barrier predictions fails for environments with |M i,⊥ | = 0 (black and purple), where reversal via non-uniform modes are favoured, as discussed before. The reductions compared to the mean-field barrier seem to be particularly strong for environments which feature "X" and "C" configurations of the perpendicular moments, see Figs. 9(f,g). This highlights the importance of considering the magnetostatic interactions with neighbouring nanomagnets during reversal to obtain the correct barrier energies.
With the exception of ∆E string sb , which is a result of the string-method simulation, the energies E mm i,↔ and J mm NN can be obtained from static equilibrium micromagnetic simulations, e.g. using OOMMF [21] or MuMax3 [22]. This makes this approach attractive to estimate more realistic switching barriers based on a pertubative decomposition of a single-nanomagnet behaviour plus a correction term due to interactions with the neighbouring moments. This approach seems valid for relatively small nanomagnets favouring reversal via uniform modes, but fails if more complex reversal mechanisms are accessible.