Magnetic reversal of an artificial square ice: dipolar correlation and charge ordering

Magnetic reversal of an artificial square ice pattern subject to a sequence of magnetic fields applied slightly off the diagonal axis is investigated via magnetic force microscopy of the remanent states that result. Sublattice independent reversal is observed via correlated incrementally pinned flip cascades along parallel dipolar chains, as evident from analysis of vertex populations and dipolar correlation functions. Weak dipolar interactions between adjacent chains favour antialignment and give rise to weak charge ordering of ‘monopole’ vertices during the reversal process.


Artificial spin ice
Artificial spin ices are lithographically patterned arrays of single domain nanomagnets [1,2]. The elongated shape of elements on a number of interpenetrating sublattices forms a 2D system of interlinked vertices at which Ising-like dipole moments meet with incompatible interactions. Systems with two sublattices form square ices, whilst those with three form kagome ices. They are analogous to 3D spin ice pyrochlore crystals [3,4], in which rare earth magnetic moments map directly on to the proton ordering of water ice [5], the classic frustrated material. That mapping is in fact so robust that spin ice can be considered a 'magnetolyte' [6], supporting currents of interacting fractionalized Coulombic 'magnetic monopole' charge defects where the energy minimizing 2-in/2-out vertex ice-rules of the macroscopically degenerate ground state (GS) are violated at the ends of sequentially flipped spin chains or 'Dirac-strings' [7][8][9][10][11][12]. Naturally, the 2D artificial analogues have recently become a source of intensive interest [13][14][15][16], as artificial systems allow for vertex interactions to be uniquely controlled. Such systems can also be easily studied using direct imaging via magnetic microscopy.
The artificial systems studied to date are inherently athermal at room temperature, with elemental reversal barriers orders of magnitude greater than the thermal energy k B T . Applied magnetic fields are therefore required to promote reordering and dynamics. Attention has focused largely on the generation of well defined statistical states, and the propagation of magnetic defects. Effective thermodynamics [17] of ac demagnetization can be employed to acquire energy-minimized states [1,14], [18][19][20]. In square ice, an example of which is shown in the scanning electron micrograph (SEM) in figure 1(a), a short range bias towards the twofold degenerate GS (figure 1(b)) with random long range correlation is found in quasi-infinite systems [21]. Finite size effects have recently been reported [22], which may potentially yield much stronger GS order for simplified protocols [23]. True thermodynamics can allow access to long range GS ordering [16], with well defined Boltzmann-weighted monopole-like excitations above this ice-rule obeying background, however, this has only been achieved following a oneshot fabrication-stage anneal process. While indirect evidence exists for charged vertex-vertex The sublattices X and Y are shown as grey and black arrows respectively, defined by their easy axis orientation. Moments on each sublattice form a chessboard pattern of alternating alignment, while energy is minimized by the relative arrangement between sublattices. All (L, P) X,Y -type lines of moments possess an antiferromagnetlike arrangement, and all D X,Y -type moments are aligned with a ferromagnetlike arrangement. (c) The DPS defined by hard polarization of both X and Y. The experimentally applied field direction is shown as a red arrow. (d) The Y-polarized state (YPS) on which Y is hard polarized and X minimizes energy via type-A antiferromagnet-like order. (e) The 16 square ice vertices, V 1−16 , shown as dipole crosses, grouped into increasing energy types T 1−4 . The magnetic north/south poles contributed to each vertex centre are also shown (yellow/red). T 1,2 are magnetically charge neutral overall, while T 3,4 possess normalized charges of ±2 and ±4 respectively. T 2,3 possess a dipole moment, shown as a green arrow.
interactions playing a role in thermal equilibration (such monopole-string excitations having been previously addressed within a theoretical framework [24,25]), it is not clear whether a description in terms of fractionalized objects is truly viable in artificial spin ice systems. Application of dc fields can create a long range polarized state [26], the four-fold degenerate diagonally polarized state (DPS) shown schematically in figure 1(c). This state also obeys the 2-in/2-out ice rules, and its violation also creates monopole-like vertex configurations. A number of reports have recently addressed dc-field magnetic reversal processes in artificial spin ices [13][14][15], [22,27]. A variety of different phenomena have been observed, 4 owing to key differences between the patterns studied, the nature of their inter-elemental interactions, and the overlap between the finite distribution of switching fields, H s (θ), of each sublattice for a given in-plane applied field direction θ . When this overlap is significant, kagome ice reversal mediates via moment flip cascades which nucleate and propagate monopolelike configurations throughout the pattern on a charge-ordered ice-rule (but not GS [28,29]) background, incrementally pinning at sites of higher local H s (θ) [13,15]. A domain wall Coulomb charge interaction model has been shown to qualitatively describe transient states involved in cascading and propagation in nanowire networks [27,30], however, no evidence has been identified for the fractionalized 'vertex object' interactions as envisaged by Castelnovo et al [7]. Simulations have shown that such cascades can convert the square ice DPS into the GS [22] for fields applied along a diagonal symmetry axis. Further, there is current interest in correlating such dipolar reordering processes with observed global properties as measured by magnetometry [31].
In this work, we address experimental magnetic reversal of a square ice pattern, in a dc field applied along an off-symmetry direction. We observe independent sublattice reversal, charged vertex propagation via correlated flip-chain extension, and clearly identify the effects of charged vertex interactions.

Square ice
Square ice can be considered as two identical orthogonal sublattices of moments, defined as X and Y in figure 1(b) in reference to their easy axis alignment with the x and y Cartesian axes respectively. The four-moment cross-shaped vertices formed consequently possess 2 4 = 16 possible configurations, V 1−16 , conventionally grouped into four types, T 1−4 [1], based on increasing energy, as shown in figure 1(e). On a given vertex, those moments lying parallel to each other can take either N-S or N-N/S-S configurations (N = north, S = south). A convenient self-contained way to consider a vertex is therefore four converging magnetic N/S poles. T 1,2 obey the 2-in/2-out ice rules, imposing magnetic charge neutrality. The GS is formed by a chessboard tiling of V 1,2 T 1 sites, while the DPS is a 100% tiling of one single form V 3−6 of the T 2 vertices. T 3,4 are magnetically charged sites, possessing excess N or S, with normalized charge ±2 and ±4 respectively, often referred to as 'monopole' vertices in qualitative analogy with the charge defects of bulk spin ice [7]. T 2,3 have inherent dipole moments.

Methods
A square ice sample was prepared by electron beam lithography, using a Helios Nanolab SEM system, of ZEP520A : acetol (1 : 1) resist spin-coated on a Si substrate, vacuum evaporation, and liftoff. Elements of 100 nm × 250 nm were formed on a lattice of 500 nm lattice constant ( figure 1(a)), over multiple 20 × 20 µm 2 subarrays tiled at a spacing of ∼1 µm on a square grid. A total area of ∼1 mm 2 was formed. A structure of Ta(2 nm)/Ni 80 Fe 20 (25 nm)/Al(2 nm) was deposited. This formed islands with magnetic moments of ∼5 × 10 7 µ B , giving rise to dipolar fields of ∼10 Oe imparted on one island by each of its closest neighbours. The single domain Ising nature of the nanobars was confirmed via magnetic force microscopy (MFM), each element appearing as a dumbbell of light and dark contrast.
In-plane hold-fields of ∼1 s duration were then applied along a direction θ ∼ 10 • offset from the diagonal to the x and y axes, see figure 1  remanent magnetic configuration was followed by MFM, imaging ∼13 × 13 µm 2 areas at the centres of ∼5 closely positioned subarrays after each applied field. It should be noted that exactly the same areas were not reimaged at each field step due to the limitations of positioning the microscope field-of-view, however, varying amounts of overlap do occur between images at successive field steps. While edge effects are potentially important in finite arrays [22,23], strong evidence exists that significant quenched disorder can result in dominant bulk processes [13,15,27]. Initially, a magnetic field of H = +1.5 kOe was applied, sufficiently large to form the DPS at remanence. Increasingly negative hold-fields were then applied from −313 Oe up to −688 Oe, taking the system through a full reversal to the opposite DPS. While applied fields will disturb the Ising-nature of the nanobars, we find, as previous authors report, that we can interpret the remanent states in terms of an ideal Ising dipole picture.

Magnetic force microscopy of reversal
Example MFM images are shown in figure 2, with key configurations boxed and mapped schematically (insets) in terms of elemental Ising dipoles and magnetically charged T 3 vertices. Note that the sequence does not show strictly the same region of square ice elements in each image. In figure 3(a) we track the normalized net digital magnetization  The specific vertices are shown and referenced by colour. Both reversals occur via a peak in complementary positive/negative T 3 pairs, V 10,12 and V 7,13 , which must nucleate to mediate reversal. Note, for clarity, the lower panel of (b) shows the difference between the experimental and random reversal T 2 populations, as the fractional change is small. Further, only the populations of V 10 and V 7 are plotted in (c), within error possessing equivalent populations to their oppositely charged partners, V 12 and V 13 respectively. (This equivalence will be exact at the system level by conservation of charge.) Sequential moment chain reversal acts to transfer vertices within the T 2 group V 3 to V 4 , and V 4 to V 5 . Experimental T 2,3 vertex populations fall significantly far/short of random during reversal, indicating correlated cascade propagation. Initial and final stages of both reversals appears random-like. Weak enhancement (suppression) is observed of T 1 (4) . Strong suppression of all other vertices is also found. The values of field where M X,Y = 0 are marked by vertical orange lines, estimated as where the growing/falling T 2 populations cross. 7 for each step and those of specific interest are shown as solid lines in figures 3(b)-(d). A small number of counting errors occur due to structural defects and tip-sample interactions, however, these occur on ∼0.3% of islands imaged and do not have a significant effect on the system or statistics.
For −375 H < 0 Oe, the initial DPS is observed with m X = m Y = 1, giving a 100 % V 3 T 2 tiling. Full reversal then takes place via two largely independent, but qualitatively similar, stages: initially, reversal of Y occurs between −438 Oe and −563 Oe, followed by reversal of X between −563 Oe and −688 Oe, figure 3(a). For a given applied field angle, θ , a given element on X or Y will have an intrinsic switching field (i.e. the total field required to reverse its Ising state observed at remanence) H X,Y s (θ ). At θ = 0, it would be expected that H X s (0) = H Y s (0) , were X and Y identical under a 90 • rotation of the system. While an intrinsic anisotropy may be possible from patterning, SEM reveals no obvious structural stigmation, however, very small asymmetric artifacts are present on one corner of each element. Slight quantitative differences were observed in behaviour at different locations across the sample, e.g. variation in H X s (θ ) , which appears to vary in the subarray-patterning slow scan direction, presumably the result of a gradual drift in exposure conditions with time. Such variation is small in comparison with the separation in field of the two independent sublattice reversal events. This effect is therefore attributed predominantly to the applied field angular offset, θ ≈ 10 • , meaning that H y > H x and H Y s (θ) < H X s (θ ) . As previously described, we therefore focus here on one self-contained local region, other regions behaving in a qualitatively identical way.
At −438 Oe a low density of isolated flipped-moment events is observed on Y, attributed to elements with lowest H Y s (θ ). These correspond to nucleation events of oppositely charged T 3 pairs, V 10,12 , at the expense of two V 3 T 2 vertices, leading to a reduction in M Y . At −469 Oe a higher density of single flipped moments and associated T 3 pairs is observed, as well as longer sequentially flipped chains further separating oppositely charged T 3 pairs in the y-direction connected by V 4 T 2 vertex chains with polarization rotated 90 • to the initial V 3 DPS. While the DPS is not the GS, we can draw analogy here with charge separation in real and artificial kagome spin ice [7,15], where background and chain are of the same ice-rule-obeying vertex type, with oppositely charged poles propagating in opposite directions.
At −500 Oe, over 50% of Y moments have reversed. M Y = 0 is estimated to be at ∼ −490 Oe, indicated in figures 3(b)-(d) by dashed orange lines at the intersection points of the V 3,4 populations. Qualitatively, the state is very similar to that at −469 Oe, with a substantial population of T 3 vertices and ∼50% V 4 T 2 . The state appears more like a background of V 4 on which chains of unflipped moments now stand out (inset). Due to the coarseness of the field step, the exact reordering processes of specific groups of moments reimaged at successive steps could not tracked, therefore, it is not possible to say whether a chain of 3 moment flips is formed by sequential flipping, or by the oppositely charged ends of two separate chains on the same line of elements meeting and annihilating. At −531 Oe the flip-chains on Y have almost completely propagated out, reversing Y to m Y ≈ −1. A small number of unflipped moments remain and the state appears qualitatively like that at −438 Oe, reflected about the x-axis. A symmetry should be emphasized, resultant of the degeneracy of the DPS. T 3 pair nucleation and annihilation processes appear qualitatively as the inverse of each other. The same is true comparing the extension of a flipped moment chain with the shrinking of an unflipped moment chain. Hence, given no knowledge of field history, the initial DPS background and propagation direction are not revealed from an MFM image. Overall conservation of charge is always maintained, for example, V 10,12 are always observed in equal numbers within error, 8 small discrepancies only occurring over a MFM image due to single poles of such pairs having propagated across the image edge boundary.
Subsequently, for increasing fields, the sublattice X then reverses via the symmetrically equivalent processes-nucleation of V 7,13 T 3 pairs, propagation and annihilation, with V 5 increasing at the expense of V 4 . Full reversal is achieved by −688 Oe. At −625 Oe, a state of M X = 0 has been achieved, which by symmetry should possess a maximum in V 7,13 T 3 . This is also indicated in figures 3(b)-(d) with dashed orange lines. Comparing this value of field with the −490 Oe required to bring about M X = 0 implies that θ ≈ 6.9 • if we assume that X and Y are equivalent in every way, close to our nominal offset angle.
Only a small amount of overlap between the reversal of Y and X occurs, inhibiting the creation of T 1,4 vertices at all fields, as indicated by their small fractional populations. As an example, T 1 vertices can only form if a chain of sequentially reversed moments crosses sublattices, in a manner shown inset in figure 2 at −531 Oe, where the propagation of charged sites from their straight line paths along the y-direction has been diverted in to the x-direction, allowed to occur by the conversion of V 10,12 into V 7,13 . Such infrequent events are negligible in the reversal regime accessed. While the two sublattices do interact via dipolar coupling, by symmetry there is no net effect of the frozen polarized X(Y) sublattice on Y(X) (except negligibly small long range interactions), therefore, the same one-dimensional lines of charge carrying propagating N-N and S-S configurations would be expected if X(Y) were removed. It is likely, however, that the frozen sublattice imparts additional disorder on the reversing sublattice via local variations in dipolar fields.
It is conceivable that the reversal-mediating processes identified could occur in a noninteracting system, the reversal occurring by random flipping events determined by the local values of H X s (θ ) and H Y s (θ ), presumably randomly allocated across the pattern and distributed about H X s (θ ) and H Y s (θ ) respectively. To test this, we generated ideal maps of moments with set fractions in each sublattice flipped at random from the initial DPS, corresponding to the experimental M states. (We average over ten 40 × 40 element vertex maps for each field step, although the results are not sensitive to the absolute size). The resultant random reversal vertex populations are plotted with dashed lines in figures 3(b)-(d). It is clear that during both reversals, the T 2 populations exceed what would be expected for random arrangements of moments, as emphasized by the excess of percentage populations plotted in the lower panel of (b), whilst the T 3 populations similarly fall short of what would be expected. (This is not surprising since the almost total absence of T 1 and T 4 vertices implies that the T 2 surplus should be of the same magnitude as the T 3 deficit.) This indicates that sequential chain flipping is significantly correlated-dipolar interactions bias the system towards reversal via incrementally pinned correlated cascades. It also appears that the initial and final stages of both experimental reversals follow a random-like trend, implying that the nucleation and annihilation events are random-like processes. Y and X reversal appear to be quantitatively different, the latter appearing to have greater suppression of T 3 vertices. Interestingly, the small populations of T 1,4 that occur are slightly enhanced/suppressed with respect to random respectively, reflecting their favourable/unfavourable moment configurations.
Making a visual comparison with a recent report of chain propagation on a kagome ice array [15], while our stronger isolation of sublattice switching strongly suppresses instances of potential 'monopole-trapping' configurations (e.g. figure 2, −531 Oe inset), much weaker cascade correlation is apparent. This can be seen, for example, by comparing our M X = 0 state, figure 2 −625 Oe, with the state following application of 99% the coercive field of the kagome pattern, in which notably longer coherent chain lengths are found. This is attributable to quenched disorder giving rise to a higher density of random-like nucleation and pinning events in our square ice pattern.

Dipolar correlations
To further explore these spatial correlations, six correlation functions C(n) = m i · m i±n i can be calculated, with C = (L , P, D) X,Y for nth nearest neighbour dipole pairs, for L, P, or D correlations on sublattice X or Y, as defined in figure 1(a). These are similar to those employed by Ke et al [21] but are here defined on X and Y independently. Perfect alignment/antialignment of all nth pairs yields C(n) = ±1 respectively. Random alignment yields C(n) = 0. Note, C does not distinguish between flipped and unflipped moments, and no correlations between X and Y are considered as their reversals are almost totally unmixed. The observed C(n) are plotted for Y and X in figures 4(a) and (b) (solid lines with symbols), as well as those calculated for the random moment reversal sequence (dashed lines).
To aid discussion, it should first be noted that all random C X(Y) behave identically with state and have no dependence on n, and hence we define and plot simply R X(Y) = L X(Y) (n) n for clarity. For our double reversal process, the collective trend is for R Y to initially decrease from a value of +1 in the initial DPS towards 0 at M Y = 0, then to rise again towards +1 as Y reversal is complete. Subsequently, R X follows the same pattern.
Experimentally, all C X,Y (n) are indistinguishable from random, with the exception of certain short range correlations during the mid-stages of both reversals. During Y reversal, at −469 Oe and −500 Oe, enhanced L Y (1) correlation is observed above random, a greater enhancement for the smaller |M Y | state, indicating the propagation of correlated sequential sublattice Y moment flips along the y-direction. Interestingly, weak suppression of P(1) is observed at −500 Oe. This can be interpreted as P(1) neighbours favouring antialignment, a consequence of solely their direct dipolar interaction (unlike GS correlated configurations). This implies that flip chains on adjacent rows on Y weakly resist nucleating or propagating along side each other, and bias towards a state as shown in figure 1(d) is present, where Y minimizes its energy via order analogous to that for a type-A antiferromagnet [32], under the constraint of the hard polarized uninfluential X sublattice. (As an aside, this state may be expected in a square ice lattice of coupled superparamagnetic moments under an x-(y-)directed applied field.) Again, reversal of X is observed to behave qualitatively like reversal of Y, however, with an apparently stronger enhancement of short range L X compared to L Y during reversal mid-stages, either due to the θ -offset or an intrinsic patterning bias. This is observed distinctly at −625 Oe with enhanced correlation up to n = 3. As this is a demagnetized M X = 0 state, it is expected by symmetry that the effects of dipolar coupling on correlation will be strongest here. Again, suppression of P(1) is observed, now with a true anticorrelation, R X lying at 0. No significant correlation can be seen in D X at any n indicating that chain-chain interactions are weak and local. The agreement of experimental and random correlations during the early and late stages of each reversal confirms that nucleation and annihilation of oppositely charged T 3 pairs are random-like processes.

Charge density functions
Qualitative evidence for the effects of the P(1) neighbour interactions can be seen in the formation of configurations such as that shown inset in figure 2 for −625 Oe. A number of oppositely/like charged pairs of T 3 V 7,13 vertices, propagating via flip chains on adjacent lines on sublattice X, appear to be pinned/antipinned at adjacent sites, attributable to the weak resistance of the chains to propagate along side each other. This is an exciting idea as it works towards validating understanding of the system in terms of coupled vertex objects, rather than its underlying dipoles. To this end, we focus on the M X = 0 state achieved at −625 Oe, and calculate the average fractional vertex type density ρ −± (n V ) of positively or negatively charged T 3 vertices relative to all other negative T 3 vertices at a separation of n V vertex sites in the y-direction (perpendicular to propagation), and compare these to ideal random values ρ R −± (n V ) (inferred by symmetry), shown in figure 5. As a hard-polarized Y sublattice is demanded in this state, only 4/16 of the vertices are compatible with this condition and, therefore, allowed to occur: V 4,5,7,13 . For a random allocation, producing a state with M X = 0, it would be expected that at any given distance relative to a negatively charged T 3 vertex (or indeed any other reference vertex) both positively and negatively charged T 3 vertices will have an average density of ρ R −± (n V ) = 0.25, except at n V = 0, where ρ(n V ) must always be 0 or 1 respectively. Our randomly reversed maps agree closely with these ideal values. Experimentally, the general suppression of both positive and negative T 3 vertices relative to the random reversal state by virtue of correlated chain propagation is reflected in the average density of ρ −± ∼ 0.17, at all |n V | > 1, agreeing closely with that observed in figure 3(c). For n V = ±1, however, it is observed that the suppression of ρ −+ is notably weaker, maintaining a value close to ∼0.23, while ρ --is much more strongly suppressed with a value ∼0.09. (Calculations for ρ +± yield a similar result for like/opposite charges.) From this we can directly infer that during the reversal process,  Figure 5. Fractional density ρ −± (n V ) of positive or negative T 3 vertices relative to negatively charged T 3 vertices in the M X = 0 state compared to ideal random values for nearest vertex neighbour distance n V . The schematic shows the propagation direction of a negatively charged T 3 vertex, and its n th V y-direction neighbours. The general suppression of average T 3 density from 0.25 to ∼0.17 for n V > 1 is understood by the enhanced cascade propagation processes. At n V = 1, however, a higher/lower density of unlike/like poles is observed, indicating weak charge ordering by attractive/repulsive pinning/antipinning. oppositely/like charged monopole-like T 3 vertices do indeed couple, at least over short ranges. While such vertices also possess a N-S dipole moment, figure 2(e), this component is identical for both the oppositely charged T 3 V 7,13 vertices allowed during the reversal of sublattice X and is part of the uninfluential sublattice Y. Therefore, this weak T 3 vertex coupling, which produces weak charge ordering, can be attributed predominantly to the ±2 magnetic charge carried by their N-N/S-S components.

Summary
In the propagation regime accessed by our protocol, X and Y square ice sublattices reverse independently via L-type neighbour flip chains. The ends of these chains are the N-N/S-S soliton-type objects [33] possessed by T 3 vertices, which carry the vertex charge component along adjacent 1D channels. While we do not address their interactions quantitatively, we show that the underlying coupling of the dipole lattice manifests as vertex charge ordering, adjacent like/unlike charges being weakly repelled/attracted. While the observed charge coupling is only evident at n V = ±1, long range interactions also can play an important role in such systems [20,23], which may manifest as long range vertex interactions, and accessing a strong coupling/low disorder regime is desirable to aid their observation. Pairs of moments possess an interaction that weakens as they get further separated, however, the distribution widths of properties imparted by quenched disorder (e.g. the H X,Y s (θ) distributions) are constants of the system. Interactions between dipolar neighbours of increasing separation therefore become quickly swamped by quenched disorder. Longer range dipolar correlations, and therefore vertex correlations, resultant of direct interactions consequently do not sum coherently. A full quantitative understanding could require a model incorporating charge and dipolar vertex components.
It should be noted that modification of ideal Ising dipole behaviour, e.g. due to internal or external fields, is not explicitly revealed in this study, however, non-uniform island magnetization is possible in elements of similar dimensions [22] and symmetry-breaking field history dependent remanent vertex configurations have been reported in kagome networks of significantly longer NiFe nanowires [27], formed by pinned transverse domain wall charge distributions with width comparable to the wire-width. Such quasiparticles have also been shown to couple on continuous adjacent wires [34] when sufficiently close. It is possible that our islands, by virtue of their size, reverse via a more coherent rotation of magnetization, as opposed to domain wall propagation, which could impart different transient magnetic dynamics. This would certainly be the case for elements of lesser volume, of dimensions comparable with a domain wall width, which is well within the capabilities of modern nanofabrication. Nevertheless, the details of the reversal mode of individual elements are not expected to affect the remanent states of the whole system studied here.
The low [13][14][15] and high energy [14] charge ordered states previously reported in kagome ice patterns are both results of field-condensation. While the latter is resultant of independent sublattice reversal, neither occur by virtue of charge or dipolar interactions and would be observed in an uncoupled system. Furthermore, alternative reversal regimes may be accessed by altering θ . Specifically, decreasing θ to more closely align the field with the lattice diagonal would allow for increased overlap between the distributions of H X s (θ) and H Y s (θ), potentially enhancing T 1 GS tile formation [22]. Understanding the myriad of methods by which to mediate order could allow for realization of frustrated patterns as magnetic information processing devices [35,36].