Modal frustration and periodicity breaking in artificial spin ice

Switchable behaviors and writable defects in artificial spin ice have been extensively studied as nanoscale reconfigurable field sources for applications including neural networks and magnonics. Random walk style chains of correlated excitations have been discussed as potential low energy pathways for such applications, but writing controllable 1D pathways through lattices is generally challenging. We explore the modal behavior of a novel artificial spin ice design that generates magnetic uniaxial chains along a defined direction. In addition, non Ising breakdown of the magnetization within constituent nanoislands is observed, with a rupture of magnetization periodicity is measured. Finally, locally perturbing the vertices in both traditional and non Ising magnetic configurations produces both a selective writing capability, and a 2D network writability dependent on the field strength of the perturbing field point source, and its controlled motion. Such functionalities exhibit a route to high-efficiency magnonic devices via novel lattice designs and controlled magnetic behaviors.

Frustration, the inability of a system to simultaneously achieve a minimum energy within its constitutive parts, yields reconfigurable behaviors in a multitude of natural 1 and synthetic 2 environments. Ice is famously frustrated geometrically as each hydrogen atom is effectively bound to two oxygen atoms, but is not localized symmetrically between them. This creates statistically two possible positional minima for each hydrogen atom depending on the surrounding environment resulting in a residual entropy in the ground state 1,3 . Such frustrated effects can also be seen in complex physical and biological systems, e.g., disordered magnets 4 and protein folding 5,6 . Modern nanolithography techniques offer the means to reproduce atomistic spin arrangements on the nanometer scale using isolated magnetic nanoislands (NIs) with negligible magnetocrystalline anisotropy and preferred moment orientation confined by their shape, acting as "macro-spins". The freedom to position and orientate NIs to high precision has enabled researchers to create artificial models for various atomistic effects, including those with frustration. Examples include realizations of 2D Ising models 7 , engineered spin-glass based Boolean satisfiability 8 , and Artificial Spin Ice (ASI).
ASI are 2D projections of spin ice 9,10 , which are a lithographically patterned into 2D arrays of in-plane ferromagnetic NIs or nanowires. ASI are geometrically frustrated from the intrinsic positional ordering of the constituent components throughout the lattice 2,11 . Similar to spin ice, the magnetization in ASI conforms to a "two-in two out" spin degeneracy in its ground-state. Higher energy, ice rule forbidden states are also easily accessible with applied magnetic fields 12 . ASI have generated attention as potential hardware realizations for more complex systems such as 2D experimental models demonstrating ice-rule dependent colloidal interactions in solution [13][14][15] . In addition, ASI are utilized for novel topical studies, including neural networks 16 , computational logic 17 , and reconfigurable magnonics. In the latter case, ASI designs with nanoscale magnetic channels, arising from topological defects such as Dirac-strings in kagome lattices, have been proposed as functionalized spin-wave (SW) guides that are controlled from their magnetic-field history [18][19][20][21] . As two-dimensional SW transport has a preferred propagation direction perpendicular to the magnetization of in-plane ferromagnetic structures 22 , applied fields can be used to guide the SW vector through a constrained magnetic structure. Locally applied Oersted fields facilitate the "steering" of SWs without impeding their properties, with application to SW-based logic and multiplex devices 23,24 . Where currentcarrying wires cannot be utilized, permanent structures with reconfigurable magnetic modes throughout a lattice can create meta-materials with nanoscale-engineered SW pathways. Hence, alternative ASI-designs [25][26][27] , engineered defects 28,29 and locally "writable" lattices 29,30 are interesting areas of further study.
In this work, we present an ASI design that exhibits an unusual form of long-range magnetic order due to the varying vertex junction shapes and ordering throughout the array ( Fig. 1a-b). Referred to throughout the presented work as the XYrY-lattice, the intrinsic positional ordering demonstrates novel magnetic switchability in the form of distinct, confined uniaxial "channels" of alike energy configurations. In addition, individual NIs can be further excited into unusual non-periodic topological defects that do not propagate through the lattice, acting as artificial null states. We show the switchable behaviors of the lattice in different magnetic arrangements and demonstrate that vertices can additionally be individually written locally utilizing a perturbing point-source field without external applied fields. The experimental conditions are reproduced and understood through micromagnetic modelling, which enables deduction of the equilibrium configurations of the magnetization, the associated energy contributions, and the generated stray field for numerical reconstruction of magnetic force microscopy (MFM) images. This analysis identifies the most energetically favored magnetic states through the applied-field sequence, thus interpreting the obtained experimental results. In Permalloy ASI, the dominating energy term defining the remanent magnetic configuration is the magnetostatic contribution, a reasonable incorporation of the NIs' stray field landscape is therefore important. Thus, in addition to numerically modeling, a calculated vector field for selected remanent states was determined by quantitative analyses of the MFM data (qMFM) as a method to determine their suitability in novel magnonic devices and to quantify probe-sample interactions. Figure 1. MFM images (5 μm ×5 μm) at remanence of mono-modal (a) and multi-modal (b) energy configurations for the XYrY-lattice; topographic outline (white) of the lattice overlaid for reference. (c-d) Energy landscape maps (ELMs) representing vertex energy types seen in (a) and (b), respectively; orange rectangle in (c) denotes the unit cell of the lattice. Accumulative net magnetic moments (yellow arrows) show "magnetic channels" in either modal state for logical sorting or magnonics. Block color in (c-d) represents the magnetic energy configurations depicted in (e (i-iv)): green squares represent type II (T2) energy states, red and blue squares represent type III (T3) states in Y and rY vertices, respectively. The spin states within the vertices are represented as vectors in the ball-diagram, with relative x-and y-components of the magnetic moment (m) defined.

Results
Periodic energy landscapes in ASI lattice. Figure 1 introduces the ASI lattice design and the modal magnetic energy landscape as a function of the magnetic field history. The unit cell is a stretched enclosed hexagon comprised of junctions resembling Xs and Ys in a one-to-two stoichiometric ratio, with one Y rotated 180° (so referred to as rY). The lattice maintains the same spin-ice magnetic energy state degeneracy and nomenclature across all vertex shapes as defined by the simple square ASI 2 . For this design, there is globally one easy axis and one hard axis (i.e. along the y-and x-axes, respectively), confirmed by modeling (Supplementary Information S1). When magnetized along the easy axis with an applied field, By = 100 mT, the moments at all four-point vertices conform to the ice rule allowed configuration, as visualized by MFM in Fig. 1a with a low moment MFM probe (LMP), demonstrating a uniform saturation state at remanence along the unit cell long-axis. Thus, one of the degenerate 'two-in two-out' low energy ground-states is achieved and defined as Type II (T2). The state shown in Fig. 1a (i.e. characterized by a single type of ASI across the lattice) is defined here to be mono-modal.
When the XYrY-lattice is magnetized along the x-axis with a sufficient applied field, e.g. Bx > 100 mT, the NIs orientated at 45° and 135° to the y-axis magnetically switch to align along the field-direction (Fig. 1b), and conserve their magnetic state when the field is released. However, under field application protocol the magnetization in the perpendicularly aligned NIs does not alter, i.e. those aligned along the y-axis. This results in a multi-modal state where the X-shaped vertices are in the T2 magnetic state and the Y/rY-shaped vertices have the higher energy, Type III (T3) state. As a consequence of the geometric ordering of X, Y and rY junctions through the lattice, this creates channels of alike energy states along the x-direction, as shown by the energy landscape maps (ELMs) in Fig. 1(cd). Different vertex energy types are represented by colored blocks explained in Fig 1e, and yellow arrows indicate the net magnetic moment per vertex in respective magnetic configurations, demonstrating its modal characteristics for potential application in magnonics or as logic devices. The structure can be set in either configuration by application of either field vector without memory artefacts. Thus, from the shape design, switchable uniaxial channels are generated depending on its magnetic history. The vector diagrams in (e) represent the summated components of the spin states contributed by constituent NIs (q) at the different vertices, which are further described as x-and -and y-components of the vertex moment. This demonstrates the greater net moment in Y and rY vertices, compared to X-shaped vertices, and explains their bidirectional characteristic in the multi-modal state.

Mono-/Multi-modal duality statistics.
The transition between the mono-and multi-modal states was studied systematically using statistical analysis of MFM images (10 μm × 10 μm) of ASI in the center of the patterned area. The sample was initially demagnetized along the easy-axis (y) by an AC-field, and subsequently the magnitude of the external x-orientated magnetic field ( x) was increased step-wise. Between field-steps (from 0 to 213 mT in nine steps, see S2 in the Supplementary Information) the ASI was imaged by MFM under zero-field to study variations in energy states within the presented structure. To compare between the MFM images, the different frustrated vertices were binned into energy configurations, and compared using the χ 2 statistical test: here o and e are the observed and expected population frequencies for each junction's energy state, respectively.
The statistical purpose of the χ 2 test is to evaluate the degree of independence between categorical datasets by comparing the acquired dataset to an expected one; the latter of which is defined by a null hypothesis 31 . In the present case, the null hypothesis is defined as all junctions in the lattice imaged are energetically homogeneous, i.e. the same vertex energy type. Therefore, e is essentially the total population per vertex junction within the scan area, and o is the absolute difference between the populations of the different energy states. This means that χ 2 measures the divergence from homogeneity for magnetic states per vertex shape across a population-based dataset as a function of Bx. Here, the function log10(χ 2 + 1) is used to accentuate the change in frequency of individual magnetic states, and remove mathematical errors associated with the logarithm of zero (Fig. 2a). Thus, any deviation from log10(χ 2 + 1) = 0 demonstrates a global transitional state where two or more energy states are observed per unit vertex, i.e. the lattice is neither mono-or multi-modal. In these named cases all vertices are uniform when subcategorized by junction shape. In Fig. 2a(i), at Bx < 30 mT, all three junction shapes are in a uniform T2 energy state, i.e. mono-modal [X(T2) Y(T2) rY(T2)]N structure, across the MFM scan area. This is the magnetization state as set by the y-axis orientated pseudodemagnetization procedure (the 'pseudo' process notifies that a complete demagnetization was not achieved). Minor variance in the initial data-point (Bx = 0 mT) is due to topographical defects, but was statistically insignificant to the dataset. In this regime, the minimum switching field is not reached for any individual NI in the lattice, thus no switching was observed. For 44 < ≤ 119 mT ( Fig. 2a(iii)), multi-modal energy states are observed. Specifically: all X-and Y/rY-vertices are solely in T2 and T3 states, respectively (i.e. [X(T2) Y(T3) rY(T3)]N).
By considering the same external field sequence applied in the measurements (details are reported in the Supplementary Information (S3)), the magnetization of the XYrY-lattice as a function of Bx is studied by micromagnetic modeling. After the demagnetizing process along the y-axis, the periodic monomodal magnetization is seen (Fig. S2ii (b)), confirming the experimental results. According to the modeling, the remanent multi-modal state appears after applying 75 mT field along the x-axis direction, this is also matched closely to the experimental values (Fig. S2ii (c)).
In between the mono-and multi-modal regions exists a transition period in the experimental data that causes a "spike" in Fig. 2a(ii) where a heterogeneity in energy states across the individual vertices is detected, most prominently in the X-shaped junction. Fig. 2b demonstrates the gradual transition in 1 mT iteration steps from the mono-to multi-modal state by plotting the population frequency (P) of T3 states. X-shaped vertices undergo a stochastic two-step transition from T2 (vertically aligned) to T3, and back to T2 (horizontally aligned) as seen diagrammatically and experimentally in Fig 2c and d, respectively, Full transition between mono-and multi-modal occurs within a narrow field window of ~20 mT. Interestingly, there is a slight asymmetry/hysteresis in the switching profiles of Y and rY junctions, as the latter demonstrates a supposedly weaker switching field. This is likely from a small angular offset between the sample and the field (~0.7°).
The multi-modal state remains the dominant energy configuration of the lattice until Bx > 119 mT, as an exponentiallike increase in log10(χ 2 + 1) in Fig. 2a for the Y and rY states is identified because the periodicity across the lattice vertices deteriorates. The energy states in the MFM images at remanence show the Ising states breakdown within the constituent NIs, resulting in an NI with weaker MFM contrast, and thus stray-field; this suggests a more complex magnetization structure within the parallel longitudinal NIs. From replicating the experimental conditions in modeling, the magnetic breakdown of the longitudinal NIs appears after the application of a field of 275 mT, far greater field than the equivalent experimental dataset. This discrepancy between simulated and experimental values may originate from lattice defects and/or lattice-field misalignments, which can induce a break in uniformity and decrease the switching energy barrier.
Non-periodic Landau states. The model revealed non-Ising NIs had a magnetization pattern akin to a vortex/Landau structure, in agreement with the MFM results (Fig. 3) expected for the patterns. From the experimental study, in 40% of cases one NI in the longitudinal pair breaks down into a Landau state (Fig. 3b). This defect is preferable to maintaining the ASI degeneracy and switching into an antiparallel Ising magnetization (Fig.  3c) that occurs in only 4.2% of cases, which produce a multimodal state consisting of T1 and T2 energy types at the Y-rY vertices. Fig. 3d demonstrates a state (2.8%) where one NI in the pair breaks down into the Landau state and the other flips Ising-like, whereas Fig. 3e demonstrates an occurrence of a double Landau state (1.4%). After this field history, 51.6% of vertices remain in the original configuration (Fig. 3a).
In this population frequency analysis, Landau chirality and location in the parallel pair were not considered. However, a more detailed analysis of these variables is demonstrated in the Supplementary Information (S4). In that case, a greater field (Bx ≈ 600 mT) was used. The higher field excites a greater amount of longitudinal NIs into a Landau state, increasing the frequency of energy states described by Figure 3b and e. The analysis demonstrates an approximately equal split between clockwise (CW) and counterclockwise (CCW) Landau states. However, there is an apparent preference for the CW and CCW vortices to appear on the left and right of the parallel pair NIs, respectively. From the micromagnetic modeling the magnetic behavior of the ASI nanostructure was determined to be significantly influenced by the previous magnetic history, exhibiting strong hysteretic properties. As an example, when starting from saturation along x-axis, the remanence state obtained by gradually reducing the field to zero is characterized by a chaotic magnetization configuration, with the simultaneous presence of different magnetic states, as illustrated in Fig. S3ii (d) in the Supplementary Information.
To link simulation and experimental results and to quantify the stray magnetic field emanating from the ASI, MFM measurements were performed with calibrated probes and analyzed as described in Methods and Refs. [32][33][34] . The tip transfer function (TTF) quantifies the stray field derivative profile of the MFM probe below the physical tip apex at a distance ztot, which corresponds to the tip-sample distance during the calibration measurement on a Co/Pt multilayer reference film. The TTF of the used LMP is largely identical when evaluated before and after imaging the ASI. It confirms that the probe's imaging characteristic did not changed throughout the experiment. Details of the calibration measurements are given in Supplementary Information S5. An estimate of the expected MFM contrast is obtained for the three remanent states from calculating the surface charges through micromagnetic modeling and periodically extending them laterally to the dimensions of the scan area. Subsequently, the surface charge map was convolved with the averaged TTF, resulting in a calculated MFM contrast (Fig. 4a), which shows the typical blurring of the surface charges due to the volume character of the probe at tip-sample separation. While the surface charges are concentrated at the island ends and reproduce the NIs curvature, the calculated MFM contrast resembles isotropic poles with about 140 nm full-width half-maximum.
Selected line profiles through experimental and calculated data show very good qualitative agreement (Fig. 4b), and for the mono-modal state also good quantitative agreement. Here, simulated and experimental values match perfectly when a small correction factor of 1.3 is applied to the simulated data. To corroborate this finding, qMFM measurements have been performed with two additional probes (see S6 in the Supplementary Information). Despite their largely different quantitative characteristics and resulting MFM signal, calculated MFM profiles matches the experimental data well, when the corresponding TTFs are employed. In these experiments, agreement is within 10% deviation. Quantitative evaluation of the experimental and simulated MFM images for the multi-modal state (Fig. 4b (ii)) is comparatively inferior. A larger correction factor of 1.45 is required to match the experimental MFM image with the simulated values. In addition, an underestimation of the attractive regime is measured because the measured MFM contrast at the two negative poles is both stronger and laterally more confined. Contrary to this, along a profile through two neighboring positive poles, the corresponding double peak is neither increased substantially nor narrowed (see Supplementary Information S5f). This finding is attributed to an unfavorable magnetic interaction between the probe and the sample because the former experiences a strong net stray field when positioned above the T3 energy vertex. When the probe experiences a repulsive force, the MFM contrast is less affected. However, in the attractive regime contrast enhancement in MFM measurements is often observed due to a focusing of the probe's magnetization configuration in the sample's stray field or vice versa. This focusing effect is even an advantage in MFM measurements with superparamagnetic probes 35 .
The qualitative appearance of the non-Ising Landau state is convincingly estimated by the calculated MFM data, however a quantitative comparison of the profile is only possible in the close (~2 μm) proximity of the experimentally measured Landau state. Experimental and modelled profiles coincide both in peak height and shape when a correction factor of 1.5 is applied (Fig. 4b (iii)). In total, the micromagnetic simulation, which is based solely on geometry data and materials properties, is a valid description of the experimental ASI realization. The quantitative discrepancy most likely originates from a combination of two factors: a possible reduced saturation magnetization in the patterned permalloy compared to the bulk literature value used, as discussed in by Farhan et.al in ref. 36 (and references 31-32 therein); and influential probe-sample interactions described above which are not accounted for in the qMFM calculations.
Such probe-sample perturbations have been utilized across literature to locally write pathways into ASI lattices, because they have been shown to have relevant applications in creating local low-energy logical pathways 29,30 . The correction factors described above provide an approximate figure-of-merit for the perturbing probe-sample interaction. In Figure 5, we assess the ability to locally write into the lattice when set in its non-periodic state using an MFM probe with a >10-fold greater magnetic moment 33 . The XYrY-lattice was magnetized by Bx ~600 mT and scanned with a high-moment probe (HMP) in single-pass AFM mode. By changing the scanning angle, we can alter the dynamics of the perturbation field with respect to the individual islands and therefore the switching field of the NI changes according to the Stoner-Wohlfarth model 37,38 . Thus, local pathways can be established and selectively written into the lattice.
Figure 5a (i) is an MFM image of the initial magnetic configuration of the lattice, taken with the LMP before writing with the HMP in semi-contact (tapping) mode. To aid with interpretation, the MFM image was converted into an ELM in Fig. 5a (iii)), which displays the different energy states at each vertex as a colored square (Fig. 5a (ii)). The HMP scan area is displayed as the black outline at an approximate slow axis scan angle of 4.5° from the y-axis of the lattice, this angle was chosen because it is slightly off the easy-axis of the lattice/NIs thus reducing the NIs' switching field. Subsequently, the lattice was imaged again with the LMP (Fig. 5a (iv)), where a change in block color from 5a (iii) indicates a change in vertex energy. In the switching occurrences, only parallel-magnetized longitudinal NIs switched (red and blue blocks in Fig. 5a (iii)), predominantly into an antiparallel configuration, i.e. lower energy T1 (yellow) and T2 (green) states at the Y and rY vertices, respectively. Specifically, the right-hand NI in the longitudinal pair switched into the antiparallel configuration in 57% of instances, whereas the left NI switched in 14% of instances. The remaining switching events (29%) were right-hand Ising state NIs evolving into non-Ising Landau states (white blocks). The HMP did not change the magnetization of any non-Ising NIs and the two remaining parallel longitudinal pairs are at the edge of the scan area, where the probe has not passed over the full area of the NI. This contrasts Fig. 5a (v-vi), where the fast axis scan direction was along the lattice hard axis. In this instance, no transitions were measured because the required switching field along this axis was not reached. type key (ii) and energy landscape maps (ELM) before (iii) and after (iv) high-moment probe (HMP) writing step. Black outline indicates the area scanned with the HMP in single-pass tapping mode. (v-vi) Equivalent ELMs scanned before and after writing step in tapping mode, but the scan angle was along the lattice hard axis, resulting in no switching events. (b) ELMs from MFM images with LMP before (i) and after (ii) HMP writing step in constant contact mode, black outline indicates scan region of HMP. Red highlighted area indicates a T3 chain in X vertices Figure 5b represents the magnetization in the lattice before and after (i and ii, respectively) HMP writing in contact single-pass AFM mode as ELMs. Within the scan area marked by the black line, there is no specificity in writing magnetic configurations at the vertices, the states are uniformly written into the monomodal magnetization, following the direction of the scan's slow axis. This provides a non-selective "rewrite" function, which brings selected areas back to a homogeneous energy state. This contrasts areas at the edges of the HMP-scan area, where repeatable T3 energy states in the X-shaped vertices, trace the bottom of the scan area, creating a "logical pathway". When the scan is shifted up one unit cell on the right-hand side of the image, we see that the pathway follows the red-pathway. This demonstrates an ability to produce a singular 2D pathway through a complex medium as is desirable in the field of reconfigurable magnonics, complementing the 1D pathways produced by globally applied field.

Conclusion
A novel ASI lattice design was presented demonstrating two switchable modal states with large windows of field stability, separated by a meta-stable state appearing in a narrow field range. The ability to globally set uniaxial 1D magnetic channels through a lattice with separable energy configurations may prove useful in reconfigurable magnonics and magnetic logic. In addition, non-Ising frustrated states, understood through experiment and modeling, may inspire future works from their potential as magnetic "null" states from their inert behavior compared to the rest of the lattice. Future studies on controllable local generation of non-Ising states, opposed to the nonperiodic generation demonstrated here, could provide additional dimensions to ASI-based read/write applications. The revelation of non-negligible probe-sample perturbation in MFM imaging presented a unique selective-writing property of the lattice dependent on the scanning parameters, the energy state of the vertices and constituent NIs. This was contrary to contact mode, which universally wrote the scan region with a distinctive edge boundary, resembling a 2D string which could be used as a highly localized reconfigurable pathway through a complex matrix for magnonic applications.

Methods
Fabrication. ASI nanostructures were patterned using e-beam lithography with a combination of lift-off and Ar ion etching from a sputterdeposited film with a Si/SiOx(300 nm)/Py(25 nm)/Pt(2 nm) film architecture. The individual NIs are 380 nm ×100 nm laterally in a 100 μm ×100 μm array with a stadium geometry.

MFM.
Measurements for the qMFM and the statistical study were performed using the Bruker Dimension Icon ® scanning probe microscope and were conducted at zero-field. The scan location was consciously not fixed between scans to assure effects were globally observed. The detailed gradual field study was conducted with the NT-MDT NTEGRA Aura with in situ in-plane electromagnet. MFM scans were conducted in the same lattice location at zero-field after application of the field states. The MFM probes used for imaging the ASI structures was an NDT-MDT MFM_LM probe whereas the writing step was conducted with Nanosensor PPP-MFMR high-moment MFM probe. The quantitative evaluation of the sample's stray field pattern by MFM followed the procedure outlined by Vock et al. 32 and as applied in Puttock et al. 33 . Reference measurements on a flat (Co/Pt[100]) multilayer sample were conducted before and after measurements of the ASI structure. All experimental parameters were identical to those used during the ASI measurement. MFM measurements were quantitatively predicted by a forward convolution of the stray field patterns calculated via micromagnetic modelling and TTF. The reliability of the analysis has been tested by using three different magnetic probes, which result in largely different MFM measurements (see Supplementary Information S6), but converge into a quantitatively comparable description of the ASI once the tip calibration procedure is applied.
Micromagnetic modeling. Micromagnetic modelling was performed by means of a GPU-parallelized solver for the solution of the Landau-Lifshitz-Gilbert equation in large patterned magnetic films 39 . The code implements a geometric integration scheme based on Cayley transform for the magnetization 40 and a Fast-Multipole based approximation for the evaluation of the magnetostatic field. This was separated into a short-range term, which describes the interactions between close NIs via Green integration, and a long-range term, which took into account the contributions from far NIs via a multipole expansion approximation 41 . The exchange field was calculated with a finite difference method able to handle non-structured meshes and thus suitable for the discretization of NIs with curved boundaries, without introducing fictitious shape anisotropy effects 42 . Micromagnetic modelling was performed to individuate the easy and hard axes of the ASI lattice, according to the analysis reported in the Supplementary Information (S1), and to calculate the equilibrium magnetization configurations and the associated energy contributions, following the experimental applied-field sequence. From the obtained magnetization spatial distributions, we also calculated the stray field patterns for the derivation of the MFM maps. The micromagnetic simulations were performed by considering a saturation magnetization of 860 kA/m, an exchange constant of 13 pJ/m and negligible contributions from magnetocrystalline anisotropy and thermal noise. To accelerate the reaching of equilibrium states, the damping coefficient is set at 0.1, according to the analysis reported in Ref. 40 .