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

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\times 170\times 3.2\,\mathrm{nm}$ and the saturation magnetization, M, to be $2.00\times {10}^{5}\,{\rm{A}}\,{{\rm{m}}}^{-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\times {10}^{-16}$ cm³.

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}_{{\rm{l}}}-{e}_{{\rm{s}}}$ = 538 J m⁻³. 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⁴ J m⁻³. 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 summary of the parameter values used in our simulations as well as the values used in the KMC simulations of Farhan et al is given in table 1.

Table 1.  Parameter values used in the present calculations as well as values used in the calculations of [10], and results obtained for an isolated island, a single ring and a double ring element of the kagome lattice.

Parameter:	Value used	[10]
Magnetization, ${\rm{A}}\,{{\rm{m}}}^{-1}\cdot {10}^{5}$	2.00	2.00
Volume, ${\mathrm{cm}}^{3}\cdot {10}^{-16}$	2.36
Anisotropy, ${\rm{J}}\,{{\rm{m}}}^{-3}$	538	618
Results:
Isolated island rotation:
Pre-exponential factor, 1 s−1	$9.9\times {10}^{8}$
Activation energy, eV	0.79
Single ring, first rotation:
Pre-exponential factor, 1 s−1	$1.0\times {10}^{9}$	$1.0\times {10}^{12}$
Activation energy, eV	0.84	0.925
Single ring, overall transition:
Pre-exponential factor, 1 s−1	$2.5\times {10}^{9}$
Activation energy, eV	0.88
Double ring, overall transition:
Pre-exponential factor, 1 s−1	$4.2\times {10}^{8}$
Activation energy, eV	0.83

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

$$\begin{eqnarray}&&E/V=-K{(\hat{m}\cdot \hat{\kappa })}^{2}-(\vec{M}\cdot \vec{H}),\end{eqnarray} \tag{ 1 }$$

where $\vec{M}$ is the magnetic moment of the island, $\hat{{m}}=\overrightarrow{{M}}/\parallel \overrightarrow{{M}}\parallel$, is a unit vector giving its direction, $\hat{\kappa }$ 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 $\vec{H}$.

Using the angles illustrated in figure 1, the energy can be rewritten as

$$\begin{eqnarray}E/V & = & -K{\cos }^{2}\theta -\parallel \vec{M}\parallel \parallel \vec{H}\parallel (\sin {\theta }_{0}\sin \theta \cos {\phi }_{0}\cos \phi \\ & & +\sin {\theta }_{0}\sin \theta \sin {\phi }_{0}\sin \phi +\cos {\theta }_{0}\cos \theta )\end{eqnarray} \tag{ 2 }$$

The energy of the island is minimized by pointing the magnetic moment in such a way that $\phi ={\phi }_{0}$ so MEPs for remagnetization transitions can be described using only the two polar angles ${\theta }_{0}$ and θ

$$\begin{eqnarray}&&E/V=-K{\cos }^{2}\theta -\parallel \vec{M}\parallel \parallel \vec{H}\parallel \cos ({\theta }_{0}-\theta ).\end{eqnarray} \tag{ 3 }$$

Zoom In Zoom Out Reset image size

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

$$\begin{eqnarray}E/V & = & -\displaystyle \sum _{i}^{6}K{({\hat{m}}_{i}\cdot \hat{\kappa })}^{2}-\displaystyle \sum _{i}^{6}({\vec{M}}_{i}\cdot \vec{H})\\ & & -\displaystyle \frac{V}{2}\displaystyle \sum _{i\ne j}^{6}\left(\displaystyle \frac{({\vec{M}}_{i}\cdot {\vec{M}}_{j})}{{\parallel {\vec{r}}_{{ij}}\parallel }^{3}}-\displaystyle \frac{3({\vec{M}}_{i}\cdot {\vec{r}}_{{ji}})({\vec{M}}_{j}\cdot {\vec{r}}_{{ij}})}{{\parallel {\vec{r}}_{{ij}}\parallel }^{5}}\right),\end{eqnarray} \tag{ 4 }$$

where the vector ${\vec{r}}_{{ij}}$ is pointing between the centers of islands i and j.

Zoom In Zoom Out Reset image size

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 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.

Zoom In Zoom Out Reset image size

HTST [21, 22] is used to calculate the rate constant for each elementary transition $j\to k$

$$\begin{eqnarray}&&{W}_{{jk}}^{{\rm{HTST}}}=\displaystyle \frac{1}{2\pi }\displaystyle \frac{{J}_{{\rm{s}}}}{{J}_{j}}\sqrt{\displaystyle \sum _{i=2}^{D}\displaystyle \frac{{a}_{i}^{2}}{{\epsilon }_{{\rm{s}},i}}}\displaystyle \frac{{\prod }_{i=1}^{D}\sqrt{{\epsilon }_{j,i}}}{{\prod }_{i=2}^{D}\sqrt{{\epsilon }_{{\rm{s}},i}}}{{\rm{e}}}^{-({E}_{{\rm{s}}}-{E}_{j})/{k}_{B}T}.\end{eqnarray} \tag{ 5 }$$

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\times {10}^{8}\,{{\rm{s}}}^{-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 pre-exponential factor is similar, $1.0\times {10}^{9}\,{{\rm{s}}}^{-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 is shown in figure 3(a) and the possible paths of elementary transitions indicated in the inset. HTST is used to calculated the rate of each of the elementary steps using equation (5). The pre-exponential factor for the elementary steps in the MEP shown in figure 3(a) is found to range from $9.7\times {10}^{8}$ to $1.0\times {10}^{9}\,{{\rm{s}}}^{-1}$. 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\,\mathrm{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)

$$\begin{eqnarray}\displaystyle \frac{{\rm{d}}{{n}}_{1}}{{\rm{d}}{t}} & = & -12{W}_{12}{n}_{1}+12{W}_{21}{n}_{2},\\ \displaystyle \frac{{\rm{d}}{{n}}_{i}}{{\rm{d}}{t}} & = & -(2{W}_{ii+1}+2{W}_{ii-1}){n}_{i}+2{W}_{i-1i}{n}_{i-1}+2{W}_{i+1i}{n}_{i+1},\\ \displaystyle \frac{{\rm{d}}{{n}}_{6}}{{\rm{d}}{t}} & = & -(2{W}_{67}+2{W}_{65}){n}_{6}+2{W}_{56}{n}_{5},\end{eqnarray} \tag{ 6 }$$

where $i=2...5$. The coefficient 12 multiplying the elementary rate constants W₁₂ and W₂₁ 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. There, it is assumed that the probability of intermediate states does not change with time

$$\begin{eqnarray}\displaystyle \frac{{{\rm{d}}{n}}_{i}}{{\rm{d}}{t}} & = & -2({W}_{ii+1}+{W}_{ii-1}){n}_{i}+2({W}_{i-1i}{n}_{i-1}+{W}_{i+1i}{n}_{i+1})=0,\\ \displaystyle \frac{{{\rm{d}}{n}}_{6}}{{\rm{d}}{t}} & = & -(2{W}_{67}+2{W}_{65}){n}_{6}+2{W}_{56}{n}_{5}=0,\end{eqnarray} \tag{ 7 }$$

where $i=2...5$. The rate of decrease in the probability of the initial state can then be written as a continued fraction

$$\begin{eqnarray}&&\displaystyle \frac{{{\rm{d}}{n}}_{1}}{{\rm{d}}{t}}=-12{n}_{1}\left({W}_{12}-\displaystyle \frac{{W}_{12}{W}_{21}}{{W}_{21}+{W}_{23}-\tfrac{{W}_{23}{W}_{32}}{{W}_{32}+{W}_{34}-\tfrac{{W}_{34}{W}_{43}}{{W}_{43}+{W}_{45}-\tfrac{{W}_{45}{W}_{54}}{{W}_{54}+{W}_{56}-\tfrac{{W}_{56}{W}_{65}}{{W}_{65}+{W}_{67}}}}}}\right)\end{eqnarray} \tag{ 8 }$$

and the overall rate constant for transitions from a ground initial state expressed as ${k}_{{\rm{tot}}}=\tfrac{1}{{n}_{1}}\tfrac{{{\rm{d}}{n}}_{1}}{{\rm{d}}{t}}=\tfrac{A}{C}$, where A and C are a combination of the elementary rates, W_ij, given by

$$\begin{eqnarray}A & = & 12{W}_{12}{W}_{23}{W}_{34}{W}_{45}{W}_{56}{W}_{67},\\ C & = & {W}_{21}{W}_{32}{W}_{43}{W}_{54}{W}_{65}+{W}_{21}{W}_{32}{W}_{43}{W}_{54}{W}_{67}+{W}_{21}{W}_{32}{W}_{43}{W}_{56}{W}_{67}\\ & & +{W}_{21}{W}_{32}{W}_{45}{W}_{56}{W}_{67}+{W}_{21}{W}_{34}{W}_{45}{W}_{56}{W}_{67}+{W}_{23}{W}_{34}{W}_{45}{W}_{56}{W}_{67}.\end{eqnarray} \tag{ 9 }$$

The lifetime of the ground state is then obtained as the inverse of the rate constant, $\tau =1/{k}_{{\rm{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\,{\rm{K}}$. This is in excellent agreement with the experimental measurements of Farhan and co-workers showing a lifetime of ca 11 s at $420\,{\rm{K}}$ (see figure 3(b) in [10]). At $320\,{\rm{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\times {10}^{9}\,{{\rm{s}}}^{-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 pre-exponential 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}^{-5}\,{\rm{s}}$ at 420 K, and an increased thickness to 4 nm brings the lifetime to ${10}^{6}\,{\rm{s}}$.

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\,{\rm{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\,{\rm{mT}}$, states I and VII no longer correspond to the ground state, but rather states III and V. When a larger field of $4\,{\rm{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.

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. As can be seen by comparing the energy landscapes in figures 3(a) and 4, the intermediate states of the double ring are not as high above the ground state as of the single ring. For example, state II in the double ring is up in energy by 0.052 eV while for the single ring the enegy increase is 0.091 eV. The presence of the junctions where three islands meet in the double ring make it easier to accomodate one of the defects formed when an island rotates from the ground state configuration. This smaller increase in the energy of the intermediate states of the double ring is also reflected in the energy of the first order saddle points.

Zoom In Zoom Out Reset image size

The pre-exponential factor for the elementary transitions in the double ring is similar as for the elementary transitions in the single ring, between 9.8 × 10⁸ and 1.0 × 10⁹ s⁻¹. From these considerations, one would expect that the lifetime of the ground states of the double ring to be shorter than that of the single ring. However, the effective prefactor is significantly smaller for the double ring, 4.2 × 10⁸ s⁻¹. The reason for this is essentially the reduced symmetry of the double ring. While there are 12 equivalent ways for transitions out of the ground state of the single ring, there are only 4 equivalent ways for transitions out of the ground state of the double ring. 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.