Magnetic interactions and reversal of artificial square spin ices

Artificial spin ices are nanoscale geometrically engineered systems that mimic the behavior of bulk spin ices at room temperature. We describe the nanoscale magnetic interactions in a square spin ice lattice by an experimentally verified model that accounts for the correct shape of the magnetic islands. Magnetic force microscopy measurements on lithographically fabricated lattices are compared to Monte Carlo simulations of the reversal process of two lattices with different lattice spacings. Lattice node statistics and correlations show significant differences in the reversal mechanism for lattices with different spacings. The effect of structural variations is also compared for the two lattice reversals.


Introduction
The origin of magnetic frustration in bulk spin ices such as Ho 2 Ti 2 O 7 and Dy 2 Ti 2 O 7 arises from the fact that the interaction energy is minimized by two-in, two-out arrangements of spins on a corner-sharing tetrahedral network of magnetic ions [1]. The study of such magnetic spin frustration has attracted much attention, as it results in remarkable properties, such as zeropoint entropy, and unexpected dynamic properties at low temperature [2][3][4]. As a result of magnetic spin frustration, the ground states in these systems are degenerate, and as such it is interesting to study how the system switches. However, since this behavior is only accessible at very low temperatures, and the spin-spin interactions cannot be observed or mapped directly to understand the correlations and formation of such states, there is growing interest in twodimensional (2D) artificial spin ices. An artificial spin ice consisting of single-domain islands patterned in a regular 2D lattice offers significant advantages for the study of frustration behavior, since the nm scale of the islands allows for direct observations of their magnetization state and interactions at room temperature using a variety of observation techniques [5][6][7]. Additionally, lithography patterning provides flexible control of the lattice geometry (e.g. square, Kagome and triangular) and the shape of the individual islands [8][9][10].
The main focus of research on artificial spin ices has been on understanding the groundstate magnetic order. Analyses of the statistics of the different configurations that form at the nodes between islands, and of correlations within the lattice, were carried out by Wang et al [11] and Ke et al [10] for spin ice lattices in which the ground state was achieved through ac demagnetizing protocols. Similar works on connected honeycomb lattices have reported on correlations in the lattice [8] and on the switching fields of individual islands [12]. Similarly, the ground state attained by thermalization has been studied to determine the statistical ordering present in the system [13]. A different approach based on potential theory and the general stray field approach has also been used to study the remanent states in square spin ices with varying lattice spacing [7]. Of particular interest has been the observation of monopole and antimonopole defects, which occur for particular node configurations with a net positive or negative magnetic moment [5,6,12] and which are connected by the so-called 'Dirac strings'. 3 The study of magnetization reversal in artificial spin ices to understand the propagation of defects such as monopole-antimonopole pairs has received attention only recently. In artificial Kagome spin ices the reversal was observed to be mediated by Dirac string avalanches and monopole propagation [5]. However, the analysis was presented for only one specific lattice spacing. In another work, the reversal was studied as a function of interparticle spacing with a focus on determining the energy barriers for magnetization switching of individual islands [14], although the details of the statistics of node formation during the reversal were not presented. The reversal of a square spin ice lattice of fixed lattice spacing was recently studied by Morgan et al [15], with the observation of the formation of chains of flipped magnetization as a result of a magnetic field applied slightly off the lattice diagonal. In addition, all the studies of reversal reported to date have been for pulsed fields rather than for a continuously varying field.
Along with experimental work on artificial spin ices, considerable efforts have been directed at the development of interaction models that provide information about the type of order in the system (short range or long range). One such model based on dipolar interactions was presented in [16,17]. Similarly, a study of the effectiveness of demagnetizing protocols through analysis of node statistics was presented in [18] for square spin ices. Models based on charge interactions have also been investigated using micromagnetic simulations to understand the nature of avalanches seen in the mesoscopic spin ice systems [19]. However, most models assume a Hamiltonian consisting of pure point dipole interactions between the islands without accounting for the individual island shape. It has been shown previously [6] that the shape of the individual islands affects the magnetization within the island, thus altering the magnetic interactions. A recent work combined experimental studies of the coercive fields during the reversal process as a function of applied field angle, with micromagnetic simulations using the true island shape [20]. Results were presented for a fixed lattice spacing, and analysis of the node formation statistics was not included.
In this work, we seek to understand the magnetization reversal process in square artificial spin ices. We develop a model based on a shape formalism in which we derive an analytical expression to compute the interaction energy between island pairs taking into account the actual shape of the island. We then explore the role of defects in the reversal process, their formation and annihilation versus lattice spacing, and their effect on the order/disorder within the lattices, aspects that have not previously been studied extensively. Experimental magnetic force microscopy (MFM) observations of magnetization reversal as a function of lattice spacing are then compared with Monte Carlo (MC) simulations by means of node statistics and magnetization pattern correlations.

Model for interaction energy
The square lattice exhibits four distinct configurations of islands surrounding a lattice node, as shown in figure 1; the degeneracy of each node type is indicated by a circled number, resulting in 16 configurations. The dominant contributions to the interaction energy arise from nearestneighbor island interactions, and the pairwise magnetostatic interactions between islands are computed using the magnetometric tensor field approach, which properly accounts for the island shape [21]. Traditionally, each island is represented by a single volume dipole and interactions are described by the standard dipolar tensor. At close proximity, however, shape effects become important, and must be taken into account. In our approach, we write the magnetostatic interaction energy between two uniformly magnetized objects of arbitrary shape versus normalized lattice parameter a/2L y for four different node configurations (defined in the inset) and two different island shapes: rectangular prism and stadium. The purple curves were computed using the dipolar approximation, whereas the blue (stadia) and orange (prisms) curves use the shape-dependent interaction energy based on the magnetometric tensor field. as where : denotes tensor contraction,μ i is the unit moment vector of island i, M 0 is the saturation magnetization (identical for all islands), V i is the island volume and N (ρ i j ) is the magnetometric tensor field, defined as the convolution of the shape-shape cross-correlation function and the dipolar tensor field. In this approach, the magnetization in each island is assumed to be uniform and is treated as an Ising variable and the shape anisotropy is included in the magnetometric tensor field via the shape function. Further details of the magnetometric tensor field and a detailed derivation of the interaction energy are given in the appendix. The energy is written as a function of the separation vector ρ i j = r i − r j , where r i locates the center of mass of island i. The complete shape information is incorporated into the magnetometric tensor field. The two shapes considered in this paper, prisms and stadia with dimensions (2L x , 2L y , 2L z ), are shown in figure 1 (inset) for a square lattice with parameter a. The experimental lattice consists of stadia, but we also present results for the prism shape since it is the closest shape to a stadium for which a complete analytical interaction energy can be computed. Figure 1 shows the normalized energy for the independent node configurations as a function of the lattice parameter a/2L y ; the energy plots were computed using the island dimensions of (L x , L y , L z ) = (100, 40, 10) nm for the prisms and (60, 40, 10) nm for the stadia. It should be noted that the rounding of the stadium ends allows for a smaller separation distance between stadia compared to prisms. Figure 1 shows four energy curves for each node type; the purple curves were computed using the dipolar approximation; blue curves represent the exact interaction energy for uniformly magnetized stadia; and orange curves represent the prisms. The interaction energy depends on the square of the volume of the islands; one of the volume factors is used for normalization of the energy on the vertical axis, the other one remains. Since the stadium and prism have different volumes, there are two separate dipolar curves for each shape in the plots. Both blue and orange curves deviate significantly from the dipolar curves with decreasing a/2L y , reflecting the importance of island shape. Only two node types have negative energy; hence, in any stable spin ice lattice made from prisms or stadia, the local node configurations are expected to be of the types I and II. Both these nodes satisfy the spin ice rule of two-in and two-out. However, the difference between the two types is that type I node results in a demagnetized state and type II node results in a saturated magnetization state. Type IV nodes are energetically unfavorable, whereas type III nodes have zero energy due to pairs of mutually perpendicular moments. For lattice spacings larger than about 5a, the dipolar approximation becomes reasonably good, whereas for close separations it underestimates the energy by as much as a factor of two.

Experimental
Artificial spin ice samples were fabricated by means of electron beam lithography as square lattices of NiFe stadia (2L x = 290, 2L y = 130, 2L z = 20 nm), and lattice spacings of a = 390 (L390) and 700 nm (L700) patterned onto the same substrate. In a previous work [6], we used high-resolution Lorentz transmission electron microscopy to visualize with high spatial resolution the magnetic structure surrounding magnetic monopole defects and their role in the demagnetization process for the L390 lattice. We have used MFM to study the dependence of the magnetization reversal on lattice spacing and defect density. The distance between the sample and the MFM tip was optimized at 50 nm to avoid superposition of topography on the magnetic domain images. A commercial variable field module (Asylum Research) was used to apply and maintain the magnetic field at the desired magnitude during scanning with minimum interference to the Co-alloy-coated tip. This allowed us to study the reversal process in a continuous manner as opposed to applying and removing a field pulse and studying the remanent state. The L390 and L700 lattices were initially saturated by applying a field along a lattice diagonal, forcing all nodes to be of type II. The field was then reversed and MFM images of a selected larger area were recorded at intervals of 5 mT. The images were analyzed digitally to determine the net magnetization within the region of interest, as well as the number statistics of node types, vortices and edge configurations as a function of the applied field. Animations showing the magnetization reversal for the L390 and L700 lattices are available in the supplementary data at stacks.iop.org/NJP/14/075028/mmedia.

Simulations
The total energy of the lattice in an in-plane applied field H can be written as where N is the total number of islands in the lattice. This expression was implemented in an MC simulation to study the reversal behavior of the lattice as a function of applied field, H, and lattice spacing, a. The simulations were carried out for several lattice sizes, and each lattice was terminated along each edge with elements parallel to the edge, to reduce the effect of the lattice edges. Free boundary conditions were used, so that edge elements have only three neighbors instead of four. Simulations were carried out for both the stadium and prism shapes, using either the dipolar approximation or the full shape-dependent interaction energies. We have assumed a Gaussian distribution of reversal field strengths for individual elements around a mean H c , with standard deviation σ : where H p is the component of the applied field along the long particle axis. This distribution was employed in the MC algorithm in the form of its cumulative function [22]: The standard Metropolis algorithm was implemented for individual element reversals. A total of 5000 MC steps (each element is visited once in each MC step) was carried out for each applied field strength, and various average parameters (energy, net magnetization and lattice node types) were computed during the final 2500 MC steps at each field level. Increasing the number of MC steps did not appreciably alter the resulting averages. For the node statistics, all results were averaged over 24 separate MC runs, initialized with different random number generator seeds. All simulations were carried out for the following Permalloy parameters and conditions: saturation magnetization M 0 = 7.958 × 10 −4 A nm −1 and temperature T = 300 K. The starting magnetization configuration was always taken to be saturated in the diagonal direction closest to the applied field direction. Animations showing the results of MC simulations for both the lattices are available in the supplementary data at stacks.iop.org/NJP/14/075028/mmedia. They show the node statistics, applied field, net magnetization of the array and the correlation plots (as described in section 4.3). in L390 the reversal is accompanied by the formation of large numbers of other node types, particularly type I. These results indicate that in the L390 lattice, the net magnetization is reduced by the formation of type I nodes, leading to local regions of zero net magnetization, whereas such regions are not formed in the L700 lattice. Consequently, achieving the ground state of the square spin ice lattice, i.e. a state consisting of two-in, two-out configurations at each node and zero net magnetization, becomes more difficult with increasing lattice spacing. Also, the distribution of node types for the L390 lattice is seen to be narrower as compared to that for the L700 lattice. Figures 2(b) and (c) show the averaged results of 24 MC runs for lattices of 32 × 32 stadia considering only the dipolar energy approximation and using the exact energy of the islands, respectively. The reversal probability of an individual island was assumed to have a Gaussian distribution. The mean value for this reversal field and the standard deviation were varied such that the simulations showed the node counting statistics in reasonable agreement with the experimental results. The values used for the plots in figures 2(b) and (c) are a mean of 44 mT and a standard deviation 2.5 mT for the L390 lattice and 65 and 5.5 mT for the L700 lattice, respectively. For the L390 lattice, there are differences between the node statistics predicted by the dipolar approximation model and the exact energy model. This can be explained because for the L390 lattice with a lattice parameter a/2L y = 3, there are significant differences between the 8 interaction energies of the node types as computed by using the dipolar approximation and exact energy (figure 1). For the L700 lattice however, the dipolar approximation and the exact energy model predict similar node statistics which are in good agreement with the experimental data. This again can be explained because for the L700 lattice with a lattice parameter a/2L y = 5.3, there is not much difference between the interaction energies as computed using the dipolar approximation and exact energy model. The MC simulations also predict a difference in the reversal behavior for the two lattices. According to the exact energy model, the number of type I nodes formed during reversal of the L390 lattice is significantly higher than for the L700 lattice, which agrees with the experimental observations. It should be noted that the MC simulations account for all the interactions within the lattice, not just the nearest neighbor interactions. For the L700 lattice, the percentage of various nodes formed during the reversal are quantitatively close to those observed during the experiments. For the L390 lattice, the predicted percentage of type III and type I nodes is slightly larger than those observed during the experiments. Our results show that the dipolar approximation fails to account for the effect of lattice spacing during the magnetization reversal of spin ices. The exact energy model, however, captures all the essential features of the experimentally observed reversal process.

Magnetization reversal of interior islands
The experimental and exact energy model results indicate that magnetization reversal in the L390 and L700 lattices occurs via different mechanisms. This can be understood in the following manner: reversal in both lattices is initiated by the formation of type III nodes because they are easily created from type II nodes by magnetization reversal of only one island. In the L390 lattice, the reversal then proceeds by propagation of type III nodes or monopole defects. As these defects move through the lattice, they lead to the creation of type I nodes with local zero magnetization, leading to an increase in N I . The smaller lattice spacing of L390 leads to stronger interactions between neighboring islands. Thus, when one island reverses its magnetization, it forces the neighboring islands to reverse as well, causing a cascading propagation of monopole defects. In the L700 lattice, due to the larger spacing, the interactions between the neighboring islands are weaker. The islands reverse independently of each other; although this results in an increased N III , there is no propagation of monopole defects and a cascading effect is not seen, resulting in low N I values.
Further, a Gaussian function was fitted to the experimental node statistics using a curve fitting program and the mean and sigma of the function were calculated as: the mean and sigma were determined to be 49 and 3 mT for the L390 lattice and 58 and 8 mT for the L700 lattice. This agrees with the parameters used in simulations where the sigma for the reversal field in L390 is smaller than that of L700. The differences can be understood by taking into consideration that all the islands in the experimental lattice are not exactly the same. There are small structural variations which occur as a result of the patterning process, which affects the distribution of the reversal fields of the islands. In the L390 lattice reversal, as a result of the stronger interactions, the structural variations are dominated by the cascading effect, leading to a narrow distribution of interior nodes. In the L700 lattice, however, since there are no strong interactions, the variations in the island structure dominate the reversal, leading to a wider distribution.

Magnetization reversal of edge islands
We also investigated the statistics of the edge nodes during magnetization reversal. Figure 3  reversal process, and the inset shows the different types of edge nodes. Similar to the results obtained for the internal nodes, significant differences are observed in the edge nodes formed during the reversal of the two lattices. In L390, e1 and e3 nodes (shown in the inset of figure 3) form before the reversal onset of the interior nodes, and their number decreases slowly well after the interior reversal is complete; for the L700 lattice these edge nodes are formed before the onset of reversal of the interior nodes, and their reversal is complete before the completion of the reversal of the interior nodes. In addition, the width of the distribution of edge nodes for the L390 lattice is significantly broader than that of the L700 lattice, which is opposite to what is observed for the interior nodes. Figures 3(b) and (c) show the node statistics for the edge nodes as predicted by the MC simulations using the dipolar approximation model and using the exact energy model, respectively. The mean and sigma of the Gaussian distribution of the reversal probability were taken to be the same as those for the interior nodes. The dipolar approximation model predicts the same reversal behavior for the edge nodes of the L390 and L700 lattices, which is not what is seen experimentally. However, the exact energy model predicts a difference in the width of the distributions for the L390 and L700 lattices, which is qualitatively similar to the experimental observations. It can be noted that the field values at which onset of reversal of the edge nodes occurs are lower in the experimental data than those predicted by the exact energy model. This is because the experimental lattice has different types of edges; two of the edges of the lattice terminate with islands having only three nearest neighbors, while the remaining two edges terminate with islands having four nearest neighbors. During the experiments, only one edge terminating with islands having three nearest neighbors was used for accounting for the reversal statistics.
The differences between the behaviors for edge and interior nodes can be understood from the fact that the edge nodes have only three interacting islands rather than four. Thus the field required to reverse the magnetization of one edge island is expected to be smaller than in the lattice interior, causing the reversal process to initiate at the lattice edges for both the L390 and L700 lattices, in agreement with MFM observations. However, complete reversal of the edges in L390 does not occur until after all interior nodes have reversed. The widening of the distribution of the edge nodes is a result of the absence of the cascading effect at the edges combined with more structural variations in the islands along the edges. In the L700 lattice, since there is no strong correlation between the islands, the edge islands reverse in the same manner as the interior islands, i.e. in a random, uncorrelated way. Thus the edge reversal begins and ends at slightly lower field values as compared to the interior nodes. The mean and sigma for the experimental distribution of edge nodes was determined by fitting a Gaussian function as 46 and 12 mT for the L390 lattice and 31 and 6 mT for the L700 lattice. Comparing these values with those for the interior nodes, it can be seen that for the L390 lattice, the mean reversal field is almost the same for interior and edge islands but the sigma is larger for the edge islands. For the L700 lattice, the mean reversal field for the edge islands is approximately half of the interior islands, but sigma is approximately the same.

Visualization of lattice correlations
It is also important to understand the type of order (long range or short range) that develops during the magnetization reversal of the spin ice lattices, and this can be determined by analysis of the lattice correlations. We have used a Fourier transform-based method to visualize these correlations during the magnetization reversal, similar to methods used in surface magnetization characterization [23]. In this method, binary magnetization maps were generated from the simulations as well as from experimental MFM images in which the islands were assigned white color if the magnetization of that island was pointing along the positive xor y-direction and black color if it was pointing along the negative xor y-direction. A 2D correlation plot was then computed from these binary images using Fourier transforms. Figure 4 shows the simulated and experimental plots at different field values during reversal of the L390 lattice. In the saturated state (field applied along the diagonal) the plots show peaks (bright regions) at regular intervals extending throughout Fourier space, indicating perfect correlation in the saturated lattice. As the field strength decreases towards zero, these sharp peaks begin to disappear because the correlation length decreases as other node types are formed. Of particular importance is the correlation plot at close to zero net magnetization: two clear correlation peaks are present in the direction of the applied field (along the diagonal) which indicate positive correlation and are indicative of stronger interactions in the lattice and a finite correlation length. The experimental plot shows a few more peaks because it does not correspond to the exact point of the net magnetization going to zero in the lattice. Similar plots computed for the L700 lattice show the absence of any peaks, indicating that the magnetization reversal process in this lattice occurs in a more random manner with an almost zero correlation length.

Summary
We have developed a model for the magnetostatic interactions in square artificial spin ices based on the correct shapes of the individual islands. This model was then used in MC simulations to understand the experimentally observed magnetization reversal. Both experiments and exact-energy simulations indicate that the reversal mechanism depends on lattice spacing. In closely spaced lattices, the reversal occurs primarily by the formation and propagation of monopole defects via a cascading effect that leads to the formation of local regions of zero net magnetization. For lattices with larger lattice spacing, monopole defects are formed during reversal, but there is no propagation and reversal occurs randomly throughout the lattice. We have shown that the model based on a dipolar approximation for interaction energy fails to account for these differences in reversal mechanism. We have also studied the reversal of the edge islands and compared it with the reversal of the interior islands. In closely spaced lattices, the structural variations in the experimental lattices do not strongly affect the reversal behavior of the interior islands but do affect the reversal of edge islands. In a lattice with larger spacing, the reversal of both the interior and edge islands is equally affected by such structural imperfections.

Appendix. Derivation of interaction energy
The magnetostatic interaction energy between two uniformly magnetized objects of arbitrary shapes can be written as whereμ i is the unit magnetic moment vector of island i, M 0 is the saturation magnetization (taken to be identical in both islands), V i is the island volume and N (ρ i j ) is the magnetometric tensor field. The energy is written as a function of the separation vector ρ i j = r i − r j between the two islands, where r i is the position vector of the center of mass of island i. The magnetometric tensor field is defined by the convolution of the overlap function C(ρ i j ) and the standard dipolar interaction tensor D(ρ i j ). The overlap function depends on the shape of individual islands as described by their shape function, D(r), which equals 1 inside the shape and 0 elsewhere. The overlap function is computed as a cross-correlation C(ρ i j ) = D i (r) ⊗ D j (−r)/V i V j , where ⊗ denotes the convolution product. When the shape function and its Fourier transform, the shape amplitude D(k), are known analytically, the overlap function can sometimes be expressed analytically. For more complex shapes, the magnetometric tensor field must be computed numerically, to obtain the interaction energy (A.1). The geometries considered in this paper are shown in figure 1(inset). For each case, there are two interaction types to be considered: between pairs of islands (1)-(1) and (2)-(2) and between pairs of type (1)- (2). Due to the square symmetry of the lattice, interactions of the type (1)-(1) and (2)-(2) are energetically equivalent. The shape amplitudes D(k) are needed for the computation of the overlap function; for a rectangular prism, ( p), with edge lengths (2L x , 2L y , 2L z ) along the Cartesian axes, the shape amplitude is given by [24] where V = 8L x L y L z is the volume and sinc(x) = sin(x)/x. For the stadium shape (s) in the same orientation, with lengths (2(L x + L y ), 2L y , 2L z ), the shape amplitude is given by The y-integral does not have a simple analytical solution, so that the magnetometric tensor field for the stadium shapes requires a numerical approach. For the rectangular prisms, however, the computations can be carried out analytically. The interaction energy for (1)-(1) prism island pairs (with uniform magnetization along the x-direction) is given bȳ Normalizing the coordinates by the length of the island, 2L x , and introducing two dimensionless parameters, τ y ≡ L y /L x and τ z ≡ L z /L x , the overlap function and the dipolar tensor can be written as C(ρ x ,ρ y , 0) = 8L 3 x τ z C (1)−(1) (ρ x ,ρ y ) with C (1)−(1) (ρ x ,ρ y ) ≡ (1 − |ρ x |)(τ y − |ρ y |) (A.5) and D x x (ρ x −x,ρ y −ȳ, 0 −z) = 1 32π L 3 x −2(ρ x −x) 2 + (ρ y −ȳ) 2 +z 2 ((ρ x −x) 2 + (ρ y −ȳ) 2 +z 2 ) 5/2 , (A.6) whereρ = ρ/L x . Substituting the terms back into equation (A.4), integrating over z and denoting p 2 = (ρ x −x) 2 + (ρ y −ȳ) 2 and q = (ρ y −ȳ) 2 − 2(ρ x −x) 2 , the final expression for the prism interaction energy of type (1)-(1) can be written as Similarly, for the type (1)-(2) prism interaction, with the type (2) elements uniformly magnetized along the y-direction, the overlap function can be determined as (in normalized coordinates) C(ρ x ,ρ y , 0) = 2τ z L 3 x C (1)−(2) (ρ x ,ρ y ), (A.8) with and the relevant dipolar tensor element as After integrating over z, the interaction energy for type (1)-(2) prisms is given bȳ