The effect of temperature and external field on transitions in elements of kagome spin ice

Transitions between magnetic states of one and two ring kagome spin ice elements consisting of 6 and 11 prolate magnetic islands are calculated and the lifetime of the ground states evaluated using harmonic transition state theory and the stationary state approximation. The calculated values are in close agreement with experimental lifetime measurements made by Farhan and co-workers (Farhan et al 2013 Nat. Phys. 9 375) when values of the parameters in the Hamiltonian are chosen to be best estimates for a single island, obtained from measurements and micromagnetic modeling. The effective pre-exponential factor in the Arrhenius rate law for the elementary steps turns out to be quite small, on the order of 109 s−1, three orders of magnitude smaller than has been assumed in previous analysis of the experimental data, while the effective activation energy is correspondingly lower than the previous estimate. The application of an external magnetic field is found to strongly affect the energy landscape of the system. Even a field of 4 mT can eliminate states that correspond to ground states in the absence of a field. The theoretical approach presented here and the close agreement found with experimental data demonstrates that the properties of spin ice systems can be calculated using the tools of rate theory and a Hamiltonian parametrized only from the properties of a single island.


Introduction
Lithographic processing and film growth technologies make it possible to construct complex magnetic systems spanning a wide range in length scale, from nanoscale to microns [1,2]. Of particular interest are artificial spin ice systems which are frustrated and have ground state entropy. Extensive studies, both experimental and theoretical, are currently being carried out on these systems. Interesting physical effects have been observed, such as magnetic monopoles [3], thermally activated changes (even melting) [4][5][6][7][8], and novel thermodynamic phase transitions [9]. The strength of the interaction between the islands, modified by their spatial separation, strongly affects the onset of thermally activated dynamics and transitions in these systems.
Detailed measurements of elements of kagome spin ice, in particular a hexamer ring and a double ring with 11 islands, have been carried out at a temperature of 320 and 420 K with observations of magnetic transitions over a timescale ranging from seconds to days [10]. This is an important model system which is useful for testing and refining theoretical tools for calculating rates of transitions in spin ice systems.
Theoretical work on transitions between magnetic states in overlayers, individual islands and assemblies of islands such as spin-ice has mainly made use of direct dynamical simulations based on the Landau-Lifschitz-Gilbert equation or Monte Carlo simulations [11]. The former has limited applicability because of the short time scale that can be simulated, limited by the fast vibrational motion of the magnetic moments. This problem is analogous to the limitations of direct classical dynamics simulations of atomic rearrangements in chemistry and condensed matter physics [12]. Energy barriers that can readily be overcome on a laboratory time scale are typically prohibitively large for dynamics on such a short time scale. The evolution of systems of interest on a laboratory time scale cannot, therefore, be simulated by direct dynamics. The Monte Carlo simulations on the other hand are purely statistical and give only thermally averaged quantities with no information about time evolution. Furthermore, the Monte Carlo simulations have so far been based on simplified Hamiltonians without direct connection to the microscopic degrees of freedom and the energy landscape of the system. The kinetic Monte Carlo (KMC) method can be used to simulate such systems over an extended time scale, but it requires as input knowledge of the mechanism and rate of the elementary transitions that can take place [10]. While it is possible to guess what such input should be, it is better to have a well defined procedure for deriving such input from the basic properties of the individual islands. This can be done by using tools of rate theory to find the mechanism and rate of the thermally activated transitions. In this way, the KMC method can be made more rigorous and consistent with basic physical properties of the system. Such approaches have been developed for systems undergoing atomic rearrangements [13][14][15], for example dynamics at proton disordered ice surfaces [16,17], but have yet to be developed for magnetic systems.
In the present article, we report results of calculations of magnetic transitions in a 6 island ring and 11 island double ring elements of the kagome lattice. The energy landscape of the system is described by a point dipole representation of each island and is characterized by finding local minima and minimum energy paths (MEPs) for transitions between the local minima. The rate of the magnetic transitions is evaluated using harmonic transition state theory (HTST) for magnetic systems. The calculated lifetime of the ground magnetic state is found to agree well with the experimental observations of Farhan and coworkers [10] when the values of the parameters in the Hamiltonian are taken from experimental estimates as well as micromagnetic modelling of a single island. The effect of external magnetic field on the energy landscape and the transition rates is, furthermore, estimated.

Single magnetic island
In order to determine the properties of an individual island, we have carried out micromagnetic simulations. Farhan et al estimated the dimensions of their islands to be´470 170 3.2 nm and the saturation magnetization, M, to be´-2.00 10 A m 5 1 [10]. We have chosen the shape of the islands in our calculations to be as shown in figure 1 with these dimensions, giving a volume of´-2.36 10 16 cm 3 .
The in-plane shape anisotropy constant, K, was evaluated as the difference in magnetostatic self energy when the total magnetic moment is pointing in the direction of the long axis of the island, e l , and when it is pointing in the direction of the short axis of the island, e s . This gives = -K e e l s = 538 J m −3 . Similarly, the out of plane shape anisotropy was evaluated and found to be about 40 times larger than the in-plane anisotropy, 2.3×10 4 J m −3 . As a result of this large difference, the magnetic moments can be assumed to rotate only in plane during the remagnetization transitions. The OOMMF software [18] was used for the micromagnetic calculations. A The magnetic island is taken to be a single-domain, ferromagnetic particle with uniaxial anisotropy. The transition mechanism is assumed to be a uniform, in-plane rotation of the magnetic moments of the atoms in the island, i.e. all atomic magnetic moments in the island are taken to be parallel and lie in the plane of the island [19]. Then, the energy of an island, E, can be expressed as is a unit vector giving its direction, k is a unit vector in the direction of the long axis of the island (see figure 1) and V is the volume of the island. The first term on the right represents the anisotropy energy and the second term represents the Zeeman energy arising from an external field  H . Using the angles illustrated in figure 1, the energy can be rewritten as sin sin cos cos sin sin sin sin cos cos 2 The energy of the island is minimized by pointing the magnetic moment in such a way that f f = 0 so MEPs for remagnetization transitions can be described using only the two polar angles q 0 and θ

Kagome single ring
A single ring element of a kagome lattice consisting of a hexamer of islands was simulated as shown in figure 2.
The distance between parallel islands was chosen to be 1 μm to mimic a system studied by Farhan et al [10]. The interaction between the islands was approximated as dipole-dipole interaction with a single dipole centered on each island. For a hexamer of islands, the energy can then be written as where the vector  r ij is pointing between the centers of islands i and j. This energy expression defines an energy surface as a function of the orientation of the magnetic moments of the islands in the system. A state is defined as a local minimum on this energy surface and is found by setting up a rough approximations of the orientations of the magnetic moments and then minimizing the energy with respect to the orientation of all the magnetic moments. Two equivalent ground states exist corresponding to clockwise and anticlockwise arrangements of the magnetic moments of the islands. The energy along a MEP for Results: Isolated island rotation: Pre-exponential factor, 1 s −19 .9 10 8 Activation energy, eV 0.79 Single ring, first rotation: Pre-exponential factor, 1 s −11 .0 10 91 .0 10 12 Activation energy, eV 0.84 0.925 Single ring, overall transition: Pre-exponential factor, 1 s −12 .5 10 9 Activation energy, eV 0.88 Double ring, overall transition: Pre-exponential factor, 1 s −14 .2 10 8 Activation energy, eV 0.83 the rotation of one of the magnetic moments starting from a ground state is shown in figure 2. The MEP is defined in such a way that at each point on the path the energy is at a minimum with respect to all degrees of freedom except the direction along the path tangent. The calculations of MEPs were carried out using the geodesic nudged elastic band (GNEB) method [20]. The maximum energy along the MEP corresponds to a first order saddle point on the energy surface and gives an estimate of the activation energy of the transition within HTST. The rotation of one of the magnetic vectors creates two points where magnetic moments meet in a way that increases the energy, i.e. two 'defects' are formed (states labeled II in figure 3). All other configurations where additional magnetic moments have been rotated without introducing more defects where generated and the state corresponding to a local minimum in energy found. Then, the GNEB method was used to find MEPs between local minima corresponding to successive rotations. Both clockwise and anticlockwise rotations were studied. The climbing image algorithm was used to converge the highest energy image on the saddle point. The MEPs are the paths of highest statistical weight for the magnetic transitions and show the mechanism of the transitions. The essential part of the energy surface of the system is characterized by the local minima corresponding to configurations with only two defects in the hexamer and the MEPs connecting those minima.
HTST [21,22]  Here, the subscript j refers to the energy minimum corresponding to the initial state and s refers to the saddle point. J stands for Jacobian, ò for eigenvalues of the Hessian, a for coefficients in the expansion of the velocity at the saddle point in terms of the eigenvectors of the Hessian, and E the energy. The sum over eigenvalues at the saddle point excludes the negative eigenvalue corresponding to the unstable mode, i.e. direction for which the energy is maximal at the saddle point.
For a single island, with energy given by equation (3) and parameters listed in table 1, the activation energy for a rotation of the magnetic vector is found to be 0.79 eV. The calculated pre-exponential factor turns out to be small,´-9.9 10 s 8 1 , three orders of magnitude smaller than has been assumed in previous Monte Carlo simulations of this system [10].
The MEP for the rotation of the magnetic vector of one of the islands in a kagome hexamer is shown in figure 2. In this case, clockwise and anticlockwise rotations have equivalent MEPs, with equal activation energy and rate constants. The rate of transition between the two states is therefore twice the rate obtained from equation (5). The activation energy is quite a bit higher than for an isolated island, 0.84 eV, but the preexponential factor is similar,´-1.0 10 s 9 1 . For subsequent transitions, the lowest activation energy is obtained by rotating one of the neighboring islands, (one example shown as state III in figure 3(a)). In those cases, the activation energy for clockwise and anticlockwise rotations is slightly different, by a few percent, and one of the two possible mechanisms dominates. An MEP for a complete transition from one of the ground states of the kagome hexamer to the other A reasonable approximation to the pre-exponential factor for the elementary steps can, therefore, be obtained by using the value for a single island. The highest energy saddle points, corresponding to rotations between states III and IV, and between states IV and V, have energy 0.885 eV above the ground state (see figure 3(a)).
The overall rate from one of the ground states of the kagome hexamer to the other can be calculated by combining the rates of the various elementary steps in a master equation, i.e. a set of differential equations for the time derivative of the probabilities, n i , of finding the system in each of the various energy levels, i (roman numerals are used in figure 3) = -+ = - where = i 2 ... 5. The coefficient 12 multiplying the elementary rate constants W 12 and W 21 takes into account the number of different ways to get between the two levels, and similarly for the coefficient 2 multiplying some of the other rate constants, as illustrated in figure 3. For transitions between states I and II as well as between states VI and VII, clockwise and anticlockwise rotations have the same rate, so there are two equivalent, parallel mechanisms. But, for transitions between the other pairs of states either clockwise or anticlockwise rotation dominates. Transformations back from the final state, i=7, are not included since the goal is to calculate the lifetime of the ground states.
In order to obtain a closed expression for the overall rate, the steady state approximation is invoked.  where = i 2 ... 5. The rate of decrease in the probability of the initial state can then be written as a continued fraction The lifetime of the ground state is then obtained as the inverse of the rate constant, t = k 1 tot . By using the calculated MEPs and HTST to estimate the rate constant for each elementary step given the Hamiltonian in equation (4) and the parameters derived for a single island, a lifetime of τ=14 s is obtained at 420 K. This is in excellent agreement with the experimental measurements of Farhan and co-workers showing a lifetime of ca 11 s at 420 K (see figure 3(b) in [10]). At 320 K, the lifetime is much longer, calculated to be several days.
In order to determine the overall pre-exponential factor and activation energy for transitions between the ground states, the lifetime was calculated for a wide range in temperature, from 80 to 1000 K, using equation (9). The calculated lifetime was found to vary exponentially with 1/T, following closely the Arrhenius law. From this, the effective pre-exponential factor and effective activation energy were found to be´-2.5 10 s 9 1 and 0.88 eV. The effective activation energy turns out to be essentially equal to the difference between the energy of the highest saddle point and the initial state energy. The presence of the intermediate states reduces the preexponential factor for the overall transition as compared with the pre-exponential factor of the elementary steps, but the presence of multiple transition paths increases it. The two effects cancel out to some extent.
Variation in the thickness of the islands has a strong effect on the lifetime of the ground states of the hexamer, as observed by Farhan and co-workers [10]. If the island thickness is reduced to 2 nm, the lifetime becomes -10 s 5 at 420 K, and an increased thickness to 4 nm brings the lifetime to 10 s 6 .

Single ring in an external field
An important question is how an external magnetic field affects the properites of a spin ice system. Calculations were carried out for the hexamer ring in a field applied as shown in figure 3 with magnitude of 2.0 and 4.0 mT. The effect of the field on the energy landscape is illustrated in figures 3(b) and (c). Even such a small field clearly has a dramatic effect on the energy landscape. At a field of 2 mT, states I and VII no longer correspond to the ground state, but rather states III and V. When a larger field of 4 mT is applied, several of the of the local minima on the energy surface disappear and the corresponding states no longer exist. This is an important consideration in experimental measurements were external magnetic field may be present, on purpose or inadvertently.

Kagome double ring
Similar calculations were carried out for a double ring element of a kagome lattice consisting of 11 islands, shown in figure 4. Using the same values of parameters as for the single ring hexamer and a fit to the calculated rate over an extended temperature range, the effective activation energy for transitions between equivalent ground states of the double ring is found to be lower than for the single ring, 0.83 versus 0.88 eV. Again, this is essentially the difference in energy between the highest first order saddle point along the MEP and the ground state energy. Futhermore, there are more intermediate states for the double ring, each one offering the possibility of a return to the initial ground state. The calculated lifetime of a ground state of the double ring within the steady state approximation turns out to be longer than for the single ring, 22 s at 420 K, as the lowering of the prefactor has a larger effect than the lowering of the activation energy in this case.
Again, excellent agreement is obtained between these calculations and the experimental measurements of Farhan and coworkers which indeed observed a longer lifetime for the double ring than for the singe ring, 29 s versus 11 s [10]. The results on the double ring further illustrate that close agreement can be obtained with experiment using the present theoretical approach and parameters determined from the properties of a single island.

Discussion
The theoretical approach presented here and the good agreement found with the measurements of Farhan et al [10] on the lifetime of the ground state of the single and double ring kagome elements, demonstrates that the properties of spin ice systems can be calculated using the tools of rate theory and a Hamiltonian parametrized from the basic properties of individual islands. The parameter values used here were not obtained by fitting to the measurements of the kagome ring or double ring, but were obtained from estimates obtained from micromagnetic calculations and experimental measurements on a single island. This is a powerful approach that opens the door for rigorous simlations of melting and other interesting phenomena in extended spin ice systems. In particular, the pre-exponential factor is evaluated from HTST and found to be three orders of magnitude smaller than what had previously been assumed [10]. The calculated effective activation energy is correspondingly smaller than had previously been extracted by fitting to the experimental data while using the postulated value of the pre-exponential factor. This difference between the calculated and fitted values of the activation energy together with the difference between the calculated and previously assumed values of the preexponential factor will lead to different estimates of the rate of transitions at temperature outside the narrow range used in the fitting.
We note that the methodology used here for calculating pre-exponential factors in the Arrhenius rate expressions for magnetic transitions, HTST, has been shown to give results in good agreement with experimentally measured prefactors for remagnetization of Fe islands on a W(110) surface [23] and simulations of hysteresis loops of spring magnets [24]. Also, the calculated activation energy for annihilation of a magnetic skyrmion obtained with this methodology is in close agreement with direct Landau-Lifshitz-Gilbert simulations carried out at a relatively high temperature (where the life time is short enough for such direct simulations) [25]. Full transition state theory, without the harmonic approximation, where both the location and orientation of the transition state is variationally optimized (as has been done for atomic rearrangement transitions [26]) followed by dynamical corrections involving direct simulation of trajectories starting at the transition state (the two step WKE procedure, see review in [12]), will, however, give more accurate estimates of the transition rates and may be required for systems where the energy landscape is more complex and has low second order saddle points in comparison with the first order saddle points.