Samples were fabricated by electron beam lithography using a JEOL 6300 system and evaporative metallization, as described previously [17]. Elements were 86 nm by 280 nm in area, possessing a thin film structure of Ti(3 nm)/NiFe(26 nm)/Al(2.5 nm). Each element forms a single domain anisotropic macrospin, with a bistable dipolar moment. Samples each spanned total areas of 2.5 mm × 2.5 mm. Two sample batches were studied: batch 1 consisted of five arrays with lattice constants of a = 500 nm to 900 nm in steps of 100 nm, and batch 2 consisted of two arrays with a = 400 nm and 500 nm.

Magnetic states were prepared by the application of an in-plane magnetic field of magnitude H_a which decreased in a linear fashion from ≈ 1300 Oe to 0 Oe, at a rate r=dHadt, whilst the sample rotated around an axis normal to the sample plane at close to f = 42 Hz. A schematic of the protocol as well as an example of an experimentally measured field profile are shown in Figure 2. As for previous studies [2, 21], this protocol operates within a regime in which the micromagnetic response of the system's constituent nano elements to the time-varying field they experience is expected to occur on time scales significantly faster than the utilized rotational period, hence the dynamic response of an island's magnetization occurs effectively instantaneously. For clarity, note that the following results will be discussed from the sample reference frame. Such a protocol ensures states are initially reset with each application, with bulk array coercive fields of H_c ≈ 600 Oe, and a switching window of ≈ 100 Oe (hysteresis measurements for the a = 500 nm sample of batch 1 have previously been presented [37]; such values are typical for nanomagnets of the dimensions used here [21, 31] and vary little with a [38]). These values are also indicated in Figure 2. Magnetic states were imaged using MFM over ~5 evenly spaced areas each of ≈13 μm ×13 μm. Our images are colored (using the WSxM software [39]) such that each nanomagnet appears as a dumbbell of red and blue contrast, as drawn in Figure 1, indicating the north and south magnetic poles, respectively, allowing moment configurations to be fully inferred. Unless stated otherwise, statistical error bars presented are calculated as the standard error over the collection of images for a given sample and field protocol realization.

The square ice system is well-described by a 16 vertex model, Figure 1B. Four degenerate configuration types can be formed by the four moments converging at a vertex, T_i for i = 1–4 in order of increasing energy [2]. The GS is defined by a chess-board tiling of T₁ vertices, Figure 1, which, crucially, is twofold degenerate. This state presents a background for well-defined magnetic chain defects consisting of magnetically charged “monopole” excitations in the form of T_3,4 vertices, which possess an excess of north or south polarity, connected by energetic strings of T₂ vertices. Such chain defects can exist as isolated excitations flipped out of the GS background [40], or as part of domain wall structure separating two continuous regions of opposite GS phase [3, 17, 32]. Whilst both T_1,2 obey the 2-in/2-out ice rules, they are energetically split by the square ice geometry.

Figure 3 shows a series of MFM images of states prepared for r = 6 Oe/s, with a shown in each case. For clarity, regions continuously tiled with T₁ vertices are boxed and tinted red or blue, indicating whether they belong to one possible GS phase or the other. What is immediately apparent is the variety of states achieved. For largest a, states appear magnetically disordered. As a decreases, small GS-ordered regions emerge, becoming clear to the eye at a = 700 nm, existing within a disordered-looking matrix of higher energy vertices. By a = 500 nm, the order is typified by sizeable continuous GS domains spanning many vertices and separated by narrow domain wall structures the width of a single vertex, much like those observed following thermal ordering [3, 17, 26–29, 32] and in square ice simulations [24]. For a = 400 nm, these features are only enhanced. Here interaction strength is being varied against a close-to-constant intrinsic island property distribution (assuming little variation in patterning as a is changed).

Also quoted in Figure 3 are the normalized digital magnetization components m_{x, y} along the respective x−, y−directions as defined in Figure 1A. Here we assume identical point Ising dipoles. Previous discussion of a similar linear protocol reported a large variation in the net magnetization of states achieved for similar samples and protocol parameters, however, this observation was never explained or understood [31]. Here we observe a tendency for demagnetization to improve as a decreases, attributable to stronger interactions promoting the formation of zero-moment T₁ vertices.

Figure 4 shows the corresponding vertex type populations, n_i for i = 1–4, as a function of a. Populations are non-random over all a, as n_i are not given by their multiplicities q_i, plotted as dashed lines in Figure 4A. Throughout this series, vertex populations mainly shift from n₃ to n₁, being respectively, equal to 0.36 ± 0.01 and 0.24 ± 0.02 for a = 900 nm, and 0.073 ± 0.003 and 0.628 ± 0.009 for a = 400 nm. n_2,4 stay approximately constant at n₂ ≈ 0.33 and n₄ ≈ 0. It is clear that the observed behavior with a is qualitatively similar to that of thermally annealed square ice systems: the GS domains are comparable in form [3, 41]. Furthermore, the data for a = 400 nm shows the strongest GS order yet reported following field-driven dynamics.

The maximum n₁ population yet reported following ac demagnetization of quasi-infinite square ice is 0.35 [2]. Whilst the ac demagnetization protocol was further optimized in later studies [14], this work addressed exclusively pairwise correlation between elemental magnetizations. If a correlator C is defined as ±1 depending on whether two first-nearest-neighbor elements (as illustrated in Figure 1A) are favorably or unfavorably aligned, respectively, at a = 400 nm we find a value of 〈C〉 = 0.618 ± 0.009, averaged over all such pairs imaged. This compares with roughly 0.5 ± 0.02 for the previously reported ac demagnetization study [14]. Furthermore, whilst rotating and non-rotating field demagnetization experiments on finite area square ice patterns have also been shown to generate significant n₁ populations [21, 35, 42], finite size edge effects are expected to have played a significant role in the onset of coherent ordering [22], which will be absent from our “bulk” measurements.

For completeness, note that our two a = 500 nm samples possess similar statistics, however, they are statistically distinct, indicating the subtle difference between sample batches as mentioned in previous reports [17]. The protocol was also repeated two or three times for each sample, yielding consistent results (not shown in Figure 4 for clarity).

To directly characterize the observed domain structure, we define the normalized mean areal domain size, N_d, as the average number of T₁ vertices in continuous contact within a domain (calculated cumulatively over all MFM images for a given parameter set, along with a standard error). A vertex contributes a unit area a² to a domain, as shown by a box in Figure 1A. Shown in Figure 4B, N_d steadily increases from 1.9 ± 0.1 at a = 900 nm to 11.3 ± 0.9 at a = 400 nm. The value of N_d for a random state equals 1.32 to 3 significant figure, as indicated in Figure 4B by a dotted line, which would be expected for a demagnetized non-interacting system in the limit of a → ∞ (This was calculated by averaging the statistics of pseudo-randomly generated square ice states, 200 × 200 vertices in size, over 30,000 realizations). Note, we follow recent studies of dipole domains in nanomagnet arrays [5], neglecting the finite field of view in the calculation, which will act to reduce the measured size of edge-straddling domains. While our domains could be interpreted in terms of dipoles which belong to either phase of the GS, the vertex picture is more appealing as the T₁ vertex is the natural object of interest in this study.

It is interesting to consider the means by which the observed order should form during the field protocol. The evolution of order will be governed by whether or not the net local field experienced by a given macrospin overcomes its intrinsic switching barrier as H_a evolves in time. This is indeed a non-trivial many-body problem requiring the attention of simulations or “real-time” microscopy studies, beyond the scope of our present report. Despite this, we can make general inferences based on the state-of-the-art. As we will discuss, this is a problem involving the balance of three magnetic field scales defined by interaction strength (controlled by a), QD, and r.

Two trivial field regimes must exist: one where H_a is so large that all macrospins track it, the other where H_a is too low to allow any dynamics at all. In the crossover between these two regimes the non-trivial regime of interest exists, both interactions and QD influencing its width. This region is illustrated in the plot of Figure 2. The most influential interaction is that between first nearest neighbors, as indicated in Figure 1A. The field strength H₁ can be calculated assuming point Ising dipoles with moments μ = VM_s, a prismatic volume V and the saturation magnetization of NiFe M_s = 860 × 10³ A/m. H₁ varies from ≈ 40 to 4 Oe as a increases from 400 to 900 nm, respectively. For QD, recently a width σ = 1.25H₁ was estimated for a Gaussian distribution of intrinsic switching fields in similar patterns with a = 400 nm [21], hence σ ≈ 50 Oe. This estimation is expected to be representative of all a studied, as similarly fabricated patterns studied by MOKE magnetometry were reported to show little variation in the width of the field range over which their net magentizations reversed as a function of a [38], implying that extrinsic perturbations to magnetic reversal of the constituent nanoelements, such as pairwise interaction disorder, are a much lesser contribution to this width.

In dc field reversal, a picture of propagation of defected vertices across a background state via the sequential flipping of underlying moments is often invoked [12, 35, 36, 43]. In our protocol, we can consider there to be a background of “loose” spins which track the field. As H_a drops through the critical field window, increasing numbers of moments will lock out due to local stabilization. Initial nucleation of order will occur at random sites defined by QD. Correlated order can then form around such nucleation sites due to dipolar interactions. Simulations have showed that at optimal values of H_a spin flip dynamics increment every 1/4 turn in constant-magnitude rotating field protocols [21, 22], proceeding via the creation, propagation and annihilation of oppositely charged T₃ vertex pairs, and the same is anticipated for our linearly decreasing field profile. The field step every 1/4 turn, ΔH_a, is hence a third relevant scale, influencing which moments might or might not freeze. For r = 6 Oe/s and f = 42 Hz, ΔH_a = 0.036 Oe. Tuning r is hence expected to be a means of controlling the order in a given sample, crudely defining an “annealing time.”

To further study these ideas we explore the effects of parameter r on the states formed. We test this for two samples from batch 1, the densest and sparsest patterns, with a = 500 nm and 900 nm, respectively. The crucial difference between these two samples is the interaction strength. Both should, however, possess a similar intrinsic distribution of moment properties, hence, the effects of QD should be much more prominent for the latter sample, in which interelemental coupling is relatively weaker. Rates of r = 0.3 Oe/s, 3 Oe/s, 6 Oe/s, 60 Oe/s, 300 Oe/s, and 600 Oe/s were used.

Figure 5 shows MFM images for a selection of four values of r for the a = 500 nm sample, and two images for the a = 900 nm sample with smallest and largest r, respectively. All states were found to be adequately demagnetized, as quoted, with a consistently higher moment for a = 900 nm. Again, a shorter a has resulted in better demagnetization. Figure 6 shows (A) n_i and (B) N_d for both samples.

For r = 0.3 Oe/s and a = 500 nm, the strongest GS order and largest average domain size are observed, again with narrow domain walls. Here n₁ = 0.62 ± 0.09 and N_d = 11.3 ± 1.4 vertices (statistically similar to the state of a = 400 nm and r = 6 Oe/s discussed previously). As r is increased, the domains shrink in size. Whilst continuous T₁ domains are never fully suppressed to a random-state size, a messy matrix of T_2,3 structure develops by r ≈ 600 Oe/s. Here, n₁ = 0.41 ± 0.01 and N_d = 3.4 ± 0.2.

This observation is consistent with the picture that more active rotations within the critical field window allows for increased incremental propagation of T₃ vertices, which has allowed them to travel further and find opposite charges with which to annihilate, increasingly forming coherent GS domains in their wake. Once domains are established, they may further evolve or coarsen by motion and annihilation of T₃ vertices along their walls [23, 32]. As r increases, there is less opportunity for T₃ vertices to travel before becoming trapped by the decreasing of H_a. Interactions are still clearly important at r = 600 Oe/s, as few T₄ vertices are found, however, the system is simply not active for long enough to allow for strong GS order formation. This state is not dissimilar to the a = 600 nm and 700 nm states for r = 6 Oe/s discussed previously. Tuning ΔH_a appears to influence the extent to which interactions can generate correlated GS order.

For the a = 900 nm sample only very weak effects are observed. For all r, GS order is weak. Whilst there is some evidence of GS order enhancement on reducing r, this is far from conclusive. At r = 600 Oe/s, ΔH_a = 3.6 Oe which is ≈ H₁(a = 900 nm) hence local variations in field due to interactions may not be allowed to play a significant role in ordering. Interactions are, however, not irrelevant as n₄ is suppressed for all values of r. Below r = 6 Oe/s, even though ΔH_a « H₁, QD must dominate the behavior via the local variation in macrospin properties as σ » H₁. The effect of QD can be emphasized by comparing states for parameters sets (a = 500 nm and r = 3 Oe/s) and (a = 900 nm and r = 0.3 Oe/s). ΔH_a/H₁ are very similar in these two cases, however, the final states achieved are clearly very different. Neither state is random, however, the relative scale of QD is much greater in the latter.

It is further interesting to compare these results with previous reports of ac demagnetized square ice states. It was reported that ac demagnetization was notably superior to linearly ramped field demagnetization [31], despite few details of state statistics being given. Here, using a linear field, we report a closer approach to the GS than those reported for ac demagnetization [2, 3, 14]. Whilst this at first appears to be contradictory, this can also be understood by considering the effects of QD. The samples used in our current study which have the strongest GS order possess less QD relative to the interaction strength, meaning that stronger coherent ordering has been allowed. Evidence for this can be found by comparing various electron microscopy images of experimental square ice patterns previously reported [41, 44].

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.