Dynamically Driven Emergence in a Nanomagnetic System

Emergent behaviors occur when simple interactions between a system's constituent elements produce properties that the individual elements do not exhibit in isolation. This article reports tunable emergent behaviors observed in domain wall (DW) populations of arrays of interconnected magnetic ring‐shaped nanowires under an applied rotating magnetic field. DWs interact stochastically at ring junctions to create mechanisms of DW population loss and gain. These combine to give a dynamic, field‐dependent equilibrium DW population that is a robust and emergent property of the array, despite highly varied local magnetic configurations. The magnetic ring arrays’ properties (e.g., non‐linear behavior, “fading memory” to changes in field, fabrication repeatability, and scalability) suggest they are an interesting candidate system for realizing reservoir computing (RC), a form of neuromorphic computing, in hardware. By way of example, simulations of ring arrays performing RC approaches 100% success in classifying spoken digits for single speakers.


Introduction
Magnetic domain walls (DWs) in ferromagnetic nanowires have been proposed for several von Neumann-type computing devices that store or process binary data, for example, racetrack memory, [1] random access memory, [2] and Boolean logic gates. [3][4][5] These approaches use magnetic domains (regions of continuous magnetization) and the DWs that separate domains to represent data while DW motion in wires and through any wire The ORCID identification number(s) for the author(s) of this article can be found under https://doi.org/10.1002/adfm.202008389.
In nanoscale magnetic systems, emergent behavior has been observed in artificial spin ice (ASI) systems, where simple dipolar/ exchange couplings between neighboring nanomagnets result in a variety of complex and elegant collective phenomena. [15][16][17][18][19] These emergent behaviors have been proposed for use in novel computational paradigms such as reservoir computing [16] or Hopfield networks for pattern recognition. [17] Magnetic fields have been applied to ASIs to simulate the demagnetization from thermal annealing [20] and to study high-speed dynamics, [15,19] but there are few reports of continuous field-driven dynamic complexity in ASIs. [16] To realize the potential of emergent behavior in nanomagnetic systems for device applications, systems must be found in which it is possible to drive these behaviors directly with applied magnetic fields or spin torques.
Here, we demonstrate that arrays of overlapping ferromagnetic ring-shaped nanowires subject to in-plane rotating magnetic field show emergent behavior. There have been many previous studies of individual soft ferromagnetic rings, [20][21][22] including investigations into their suitability for data storage applications. [23] The low magnetocrystalline anisotropy of typical materials, for example, Ni 80 Fe 20 (permalloy), means that the lowest energy magnetic configurations of a nanoscale ring are vortex states, with magnetization oriented around the ring circumference, and a so-called "onion" state, consisting of two magnetic domains separated by DWs at opposite points of the ring ( Figure 1A). Onion states can be formed by DW nucleation under large applied fields while vortex states can be reached by two DWs meeting and annihilating. [20] The onion state domains align to an in-plane magnetic field, and so can be made to move around the magnetic ring by rotating the in-plane field. [21] There have been studies of the magnetic configurations and response to linear field sweeps in two-ring structures, [24] threering structures, [25] and large arrays of connected ring-shaped nanowires. [26] However, these produce relatively simple magnetization dynamics that are well-characterized by standard hysteresis loop measurements. In contrast, the rotating magnetic fields used in this study result in feedback processes, where DWs nucleated within a single element can repeatedly interact with and influence the magnetization configuration of its neighbors. When combined with stochastic DW pinning processes at wire junctions this creates mechanisms of increasing or decreasing DW population, such that interactions of DWs across the array create emergent, non-linear variations of the array magnetization and DW population with the rotating field strength.
This robust, highly non-linear response of the arrays to external stimuli and "fading memory" of previous magnetization states offers the primary properties required to realize a hardware platform for reservoir computing (RC), a form of neuromorphic computing ideal for analyzing complex transient data series. [27][28][29] To illustrate this, we use a phenomenological model of the magnetic ring arrays' behavior to perform RC classification of spoken digits.

Array System and Initialization
We studied square arrays of ferromagnetic (permalloy) rings with 50% wire width overlap between adjacent rings. Figure 1B shows a schematic arrangement for an example 4 × 4 ring array. Arrays were initialized with an in-plane magnetic field pulse of strength H sat to leave all rings in an onion state magnetic arrangement ( Figure 1B). Further details of the arrays are given in the Section 5.   Figure 1C-E show micromagnetic simulations of two overlapping permalloy rings that illustrate the basic DW processes that could occur at wire junctions in larger arrays. The rings had 4 µm diameter, 400 nm track width, 20 nm thickness, and 50% track overlap at the ring junction. Each set of images shows the propagation and evolution of DW structures during rotation of an in-plane magnetic field, H, of a different strength (see Movies S1-S3, Supporting Information for the full simulations). The first sequence ( Figure 1C) shows two rings initially in the onion state, that is , with two DWs in each ring, and H = 80 Oe. The DWs propagated around both rings as the field rotated and were able to overcome the pinning potential presented by the junction; this "propagation" process left the number of DWs unchanged. At lower fields (H = 50 Oe; Figure 1D) and again initialized with two DWs per ring, the first DWs to reach the junction became pinned and were each annihilated by the arrival of another DW, forming vortex states of opposing circulation; hence, this mechanism decreased the system's DW population. At higher fields (H = 100 Oe; Figure 1E), a DW from an onion state ring was able to pass through a wire junction and repopulate an adjacent, initially vortex state ring; thus, this mechanism increased DW population.

Preliminary Modeling
The simulation temperature of 0 K in Figure 1C-E meant that no stochasticity was present. DW de-pinning from ring junctions will in practice be probabilistic within a range of fields, as noted above. [8][9][10] In principle, the rate of DW loss will depend upon the number of DWs available, and the rate of DW gain upon both on the number of (empty) vortex states as well as the number of DWs available to repopulate the vortex rings. To explore how these phenomena could result in emergent collective behavior of the arrays, we created a simple analytical model to calculate the equilibrium DW population based on the rates of DW loss and gain in each field cycle. The model assumed a homogeneous DW population in an infinite square array of magnetic rings, with DWs having a uniform probability of passing any junction, and so does not include the detailed interactions at junctions in order to simulate large arrays (see Supporting Information for details of model).
Significantly, the model predicted a minimum in DW population at an intermediate pass probability (Figure 2A), which represents an intermediate rotating field strength, and a significant change in DW population over the range of pass probabilities. This well-defined non-linear response of an array's collective properties resulting from discrete stochastic events at the nanoscale (albeit modeled here as mean-field phenomena) is a clear signature of emergent behavior. Furthermore, these behaviors are both driven and tuned by the rotating applied field strength, an easily accessible parameter both experimentally and in hypothetical applications.

Measurements of Array Behavior
To explore the emergent behaviors predicted by our model experimentally, we patterned several arrays of interconnected permalloy rings on Si substrates (e.g., Figure 2B,C; see Section 5 for further details). Rings had 2 µm diameter and 200 nm wire width, or 4 µm diameter and 400 nm wire width. Their thickness ranged from 10-20 nm and in all cases there was 50% wire overlap between neighboring rings.
We used polarized neutron reflectometry (PNR) to measure the magnetization of a 2 cm x 2 cm sample of interconnected 4 µm diameter magnetic rings (20 nm thick). For each measurement, we applied an initial in-plane field pulse of H sat = 1.9 kOe to saturate the magnetization and then applied an in-plane test field of H during which the sample was rotated through 50 revolutions to drive DW dynamics in the array. The measured reflectivity of oppositely polarized neutrons as a function of wave vector transfer (e.g., Figure S1, Supporting Information) allowed the absolute array magnetization to be measured (see Section 5 for details).
We performed these measurements as a function of test field H in the range 18-150 Oe to manipulate the pinning probability of the DWs at the array's junctions in a manner similar to that in the equilibrium model.
The PNR measurements ( Figure 2D) showed the array magnetization to go through a minimum of close to zero magnetization, and with a similar form of response to the field dependence of DW population predicted by the equilibrium model ( Figure 2A). This agreement in the character of the measurement and model suggests that the basic phenomena at the heart of the model (i.e., an equilibrium between DW loss and gain mechanisms), were likely to be responsible for the functional form of the array magnetization's measured field dependence.
While the PNR measurements offered a quantitative measure of the array's magnetization, the relatively slow acquisition time of each spectrum meant that it was impractical to measure the temporal evolution of the array's magnetization over a continuous rotating field sequence. To achieve this, we used magneto-optical Kerr effect (MOKE) magnetometry, which provided a rapid, qualitative probe of the array's magnetization. The magnetic field protocol for the MOKE experiments involved application of a large bipolar uniaxial pulse to saturate the array, followed by 25 cycles of an in-plane rotating field, H ( Figure S2, Supporting Information). The MOKE signal was recorded throughout the entire field sequence, thus allowing the temporal evolution of array magnetization to be investigated.
MOKE signals were first normalized (to ±1) to signal levels representing states where the array was entirely populated with onion states, and then characterized by the average signal, S av , and maximum difference in signal, S pk-pk , across each field cycle ( Figure S2, Supporting Information). The time-varying signal component in each field cycle, S pk-pk , was due to mobile DWs (in onion states), and so the normalized fraction of these was given by n mobile = S pk − pk /2 . Pinned DWs resulted in static onion states, which contributed a net magnetization and, therefore, Kerr signal; the population of rings with pinned DWs was therefore calculated as n pinned = |S av | . The population of vortex states was assumed to make up the balance of ring states, that is, n vortex = 1 − n mobile − n pinned .
The relative proportions of vortex states and onion states with mobile/pinned DWs after 25 cycles of the rotating magnetic field are shown as a function of rotating field strength in Figure 2E. These show a field region between ≈35 and 50 Oe through which the number of pinned DWs decreased to zero and the number of mobile DWs increased until they filled the array. Within this transitionary region the number of vortex states (or number of rings containing zero DWs) went from approximately zero, through a maximum before returning to zero. The total number of DWs went through a minimum at H = 45 Oe, where the vortex population peaked, and occupied 60% of the array's rings. The field strength-dependence of DW population is remarkably similar to that of the equilibrium model ( Figure 2A) and to the form of the field-dependent magnetization measured by PNR ( Figure 2D).  Figure 3F) and high field (H = 50 Oe; Figure 3H) cases showed little variation from the initial DW population (i.e., rings fully populated) throughout the field sequence. This implies that in these cases all DWs either remained pinned in their initial positions (H = 35 Oe) or were all rotating synchronously with the applied field vector (H = 55 Oe). However, for intermediate fields, such as H = 45 Oe in Figure 3G, a gradual equilibration of ring state populations over many field cycles was observed (as may also  Figure S2, Supporting Information).
Our PNR and MOKE ensemble measurements showed, therefore, that the magnetization and DW population of magnetic nanowire ring arrays subject to rotating magnetic fields underwent minima at particular field amplitudes. The similarity of these observations to our analytical model predictions supports the idea of DW loss and gain processes creating a dynamic equilibrium of DW population and magnetization across a ring array, and their observed non-linear dependence upon rotating magnetic field strength. The gradual approach to intermediate DW population levels seen in MOKE measurements lends further support to this view.

Microscopy of Magnetic Configurations
The above experiments offered much insight in considering the collective behavior of magnetic ring arrays. A more complete understanding of the ensemble's magnetization dynamics, however, requires knowledge of magnetic configurations within the arrays and how these evolve over successive field cycles or under different rotating field conditions. To achieve this, we used microscopy techniques to image the magnetic states of nanowire ring arrays that had been subject to different sequences of applied magnetic field.
Magnetic force microscopy (MFM) of the ring arrays allowed the details of their local magnetic configurations to be observed. Figure 3A shows a low magnification MFM image of the same ring array as used in the PNR and MOKE experiments following a saturating field pulse H sat = 250 Oe in the direction shown. Rhombus-shaped regions of contrast at ring junctions where wires are orthogonal to H sat (the y-direction in Figure 3A and highlighted by dashed circles) indicated the presence of DWs, as confirmed by calculations of MFM response from micromagnetically-simulated onion state rings ( Figure S3, Supporting Information). The MFM image in Figure 3A also confirmed that all of the rings in the array were forced in the two-DW onion state following the initialization field pulse, as assumed in the experiments above. The MFM image also showed contrast at other junctions (some highlighted by dotted circles in Figure 3A) where magnetization had been assumed  Figure S3, Supporting Information) demonstrated that this assumption was correct and that the contrast emerged from magnetization deviating slightly from circumferential alignment in these regions.
Higher magnification images of the array following application of 25 cycles of 47.5 Oe amplitude in-plane rotating magnetic field ( Figure 3B,C) showed details of the domain configurations and the DWs themselves. These figures show that the rings were not only in vortex or onion states, but also very often in "¾" domain states where DWs (as revealed by their characteristic contrast signature) are only separated by one quarter of the host ring. Figure 3B shows that the DWs in onion (and ¾) states generally had rhombus-like appearance. However, some other features were also visible at junctions. For example, one ring appeared to contain a 360° DW, to make what is effectively a vortex state with two adjacent but non-annihilated 180° DWs. The domain structure in Figure 3C is annotated in Figure 3D in order to highlight the degree of short-range order, with identical ring configurations appearing across lines of rings. The richness of the magnetic configurations observed in the MFM images suggests that, while the insights provided by the analytical model and MOKE data presented above are very helpful in gaining an intuitive understanding of the arrays' emergent behaviors, the assumptions underlying them are somewhat simplistic. In particular, the assumptions that arrays only contain simple onion and vortex states and that the populations of these are homogeneous do not capture the full picture of complex magnetization configurations that it is possible to obtain in the interacting rings system. Indeed, the ¾ states are a necessary consequence of onion state DWs behaving independently, from which the emergent properties appear.
X-ray photo-emission electron microscopy (X-PEEM) allowed imaging of ring magnetization configurations directly and with a wider field-of-view than MFM. This allowed us to study the magnetic behaviors of whole arrays with greater ease and identify domains directly. Similar to previous experiments, we applied a saturating field pulse, H sat = 160 Oe followed by 30 cycles of in-plane magnetic field with amplitude H. Figure 4A shows the magnetization configuration for a 25 × 25 array of 2 µm diameter rings for a range of H values. At low fields, the array was largely comprised of onion state rings with domain magnetization aligned to the saturating field pulse direction. The magnetic configuration became gradually more complex as H was increased towards H = 26 Oe, with the rings in the array being distributed across onion, vortex and ¾ domain sates. For fields greater than H = 26 Oe the rings increasingly adopted the onion state configuration, until this state was uniformly adopted at H = 42 Oe. The net normalized magnetization of each image in Figure 4A is shown in Figure 4B. This shows the array magnetization was zero at its minimum, in exact agreement with the PNR measurements seen in Figure 2D (the reduced values of H at which the characteristic dip was observed is explained by the reduced, 5 nm thickness of permalloy used for X-PEEM measurements, compared with 20 nm for the other experiments). However, the array's domain configuration was clearly more complex than consisting of the simple vortex or onion ring states, thus explaining how the minimum number of DWs present in the array can be non-zero (as observed in the MOKE measurements and mean-field model) when the net magnetization of the array is zero (as seen in the PNR measurements). The X-PEEM image of an 8 × 8 array of larger (4 µm) diameter rings in Figure 4C (H = 33 Oe) more clearly illustrates the range of possible local magnetic configurations. Onion states (with 180° DW separation) are observable in three orientations (visible as rings either with uniform color or alternating blue/ red quadrants) and vortex states of each circulation are seen. Furthermore, there is a wide range of ¾ ring configurations in different orientations.
In addition to allowing a detailed understanding of how the nanoscale behavior of the array leads to the emergent properties observed at ≈100 µm length scales, the X-PEEM images also offer insight into the manner by which the local interactions give rise to other emergent behaviors at intermediate length scales of a few micrometers. For example, the images in Figure 4A with H ≤ 23 Oe show diagonal regions of continuous magnetization that are reminiscent of Dirac strings in magnetostatically-coupled artificial spin ice (ASI) systems. [18] The low values of H used to obtain these images indicate that DW pinning from junctions on each field cycle had a relatively low probability and that domain growth in the diagonal directions is clearly the energetically preferred route.

Phenomenological Simulation
The length scale of the ring arrays makes modeling with existing approaches challenging: the tens-of-micrometers lateral dimensions of the arrays are too large for micromagnetic modeling to be used, and the thermally activated pinning events occur over much longer timescales than can typically be simulated by micromagnetic approaches. However, our microscopy results have shown that local magnetization configurations are critical to the arrays' behavior, and thus very simple macroscale models (such as presented above) are unlikely to capture their behaviors accurately.
To overcome these difficulties, we have developed a phenomenological model that maps the progress of each DW in a square-lattice ring array of arbitrary numbers of rings (see Section 5 for further details). DWs were allowed to propagate through continuous wire regions by alignment to a rotating field vector. When DWs arrived at wire junctions, further propagation was tested against a field-strength-dependent probability, P, that was characterized by the (zero-Kelvin) deterministic switching field H 0 , zero-field energy barrier E 0 , and a factor α to scale the effect of applied magnetic field on the energy barrier (see Section 5). Alternatively, the model used different probabilities for one (P 1 ) and two DWs (P 2 ) at a junction to reflect the different de-pinning fields seen in micromagnetic modeling ( Figure 1C,E). P 1 was calculated from H 0 , E 0 and α, and P 2 as a simple proportion of P 1 . H 0 and E 0 values were assigned to each junction individually, so that arrays could have uniform or distributed property values. The model was performed to simulate room temperature; included DW annihilation and repopulation events; treated each DW individually so that vortex, 180° onion, and ¾ states were all available; and included realistic edge and corner rings, which have fewer junctions than the "bulk" rings within the arrays. The model was conducted to replicate the field sequence used in the above experiments. The initial magnetic configuration consisted of all rings being in onion states with domains aligned in all rings to represent relaxation from a saturating field pulse (e.g., Figure 5A for a 6 × 6 ring array). This configuration has the largest possible net array magnetization. The field direction was then rotated anti-clockwise and the DW dynamics allowed to evolve over multiple field cycles. We fitted the modeled values of M after 50 field rotations for a 26 × 26 ring array to the magnetization measured by PNR ( Figure 2D) by iterating the model parameters to minimize the mean-square difference in magnetization values. The fits shown in Figure 2D were obtained from averages of 100 simulations operated at a frequency of f = 0.2 Hz, which is similar to that used in the PNR experiment. For the simple model, H 0 = 162.5 ± 15 Oe (one standard deviation), E 0 = 1.09 eV and α = 1.45. The advanced model that distinguished between one and two DWs at a junction gave a best fit with H 0 = 142.5 ± 12.5 Oe, E 0 = 1.05 eV, α = 1.1, and P 2 = 0.75P 1 .
The two-probability model gave a normalized mean square error (NMSE; normalized to variance) to the PNR data of 6.2% compared with an NMSE of 13.0% for the simple single-probability model ( Figure 2D). This indicates that it is important to consider the independent behavior of the different mechanisms shown in Figure 1C  We note that for values of H close to that for the minimum value of M in Figure 2D the model produced complex domain configurations (e.g., Figure 5B and Movie S8, Supporting Information) containing local arrangements strikingly similar to those seen by X-PEEM and MFM. Notably, X-PEEM ( Figure 4C), and phenomenological modeling ( Figure 5B) images of ring arrays following application of intermediate strength rotating field showed array edges having a more uniform magnetic configuration. This appears to have arisen due to the lack of junctions to pin DWs along the array edges and highlights that the response of arrays, or at least the field range over which magnetization varies, is likely to be sensitive to array geometry.
We used the optimized phenomenological model to calculate the time-dependent magnetization, M, in the initialization direction during field rotation. The modeled time-dependent magnetization shown in Figure 5C was from a 26 × 26 ring array, as used in our PNR studies, and for H = 60 Oe. M oscillated with the change in applied field direction as DWs rotated around the array's rings. Over and above this, there was a changing dynamic in M, with a sharp decrease over the first few cycles due to DW pinning and annihilation, its recovery as DW repopulation increases over approximately the next ten cycles, and its equilibration after ≈15 cycles in total. The noise present in the amplitude of M arose due to the dynamic nature of the equilibrium and, of course, the probabilistic nature of the local DW pinning events that create the equilibrium.
The relatively good quantitative fit between the model and PNR data ( Figure 2D) offers encouragement that the model is a good description of the magnetization dynamics of the ring arrays. The model may, in fact, be a powerful tool for investigating both detailed and ensemble ring array behavior, while being relatively fast to run (for example, simulating the dynamics within a 26×26 ring array during 50 rotations of magnetic field takes 56 s on a 3.1GHz i5 processor).

Simulation of Reservoir Computing
We have further used the phenomenological model to investigate using the magnetic rings arrays for performing "reservoir computing" (RC). [27][28][29] RC is a form of recurrent neural network that is particularly well suited to analysis of complex, timedependent signals. Conventional recurrent neural networks consist of networks of neurons that receive an external input stimulus as well a weighted stimulus from other neurons, and therefore their activity evolves in time. Weighted connections to specific output nodes are used to extract information to classify data or predict outputs. Since the activity depends on the previous states of the network, training such a network requires the expensive storing and unrolling of the system states, a technique known as back propagation through time. To avoid this complex process, in RC, the network is replaced by a fixed dynamic "reservoir" of neurons connected via random fixed weights. The reservoir performs a non-linear transformation of the input signal that effectively simplifies the problem to be solved and limits the required training to the weights from the reservoir to the output layer. Since training does not involve the input weights or the weighting of internal dynamics or connectivity, RC training is extremely efficient. The recurrent nature of RC creates a "fading memory" to new input, and from which RC's power in time-varying signal analysis is derived. RC can be achieved in materio by replacing the neural network reservoir with any physical or algorithmic system that shows a reliable non-linear response to a driving stimulus, [29] with signal output taken from different regions, as desired and as convenient. In materio RC has been demonstrated with a wide range of systems (see [29] for a recent review), including electronic, spintronic, mechanical, memristive, carbon nanotube, and even liquid systems as the hardware reservoir. The robust emergent behavior of the magnetic ring array system here makes it an interesting candidate. Given the phenomenological model appears to be a good analogue for the physical system, we used simulations to perform a voice recognition classification task of spoken digits zero to nine by eight different speakers in the TI-46 dataset. [30] Initial utterance data (e.g., Figure 6A) was split into 13 spectral bands that were concatenated to form input stream, u(t), of 1200-1500 data points, depending on the utterance data. These values were transformed into N magnetic field amplitudes where k is an arbitrary scaling value. Before each measurement, a 26 × 26 ring array with values of P 1 and P 2 equal to those above (H 0 = 142.5 ± 12.5 Oe, E 0 = 1.05 eV, α = 1.1 and P 2 = 0.75 P 1 ) was initialized into onion ring states throughout. Application of the N field cycles with successive amplitudes  Figure 6B). x(t) was then averaged to produce a reduced data stream of 52 values, x i (i: 1-52; see Section 5 and Figure S4, Supporting Information). Training to a subset of utterance data X generated a trained weight matrix W out , allowing the ring Figure 6. A) Raw audio waveform for an utterance of the digits zero (left), and three (right). [30] B) Reservoir's response to inputted signal with field center of 55 Oe and unity input scaling, time averaged over 52 windows for digits zero (left), and three (right). C) Activation on each of the 10 output nodes after linear combination of time-averaged reservoir output x t ( ) and trained weight matrix W. The bars represent activation values for the "Zero" (orange) and "Three" waveforms shown in (B). D) Network accuracy for classifying test samples across all spoken digits from eight different speakers and average for each speaker trained separately using time-averaged reservoir outputs x t ( ) (red), and time averaged raw inputs u t ( ) (blue). E) Heatmaps of accuracy across all digits for a single speaker for different central magnetic fields and input scaling factors, using a single (left) and ten (right) training samples per digit.  reservoir's DW population characteristics from an unseen "test" subset of data to generate activation values against each digit being tested for ( Figure 6C), and identify the most likely digit being spoken.
RC with the ring array reservoir achieved 99.4% accuracy for a single speaker ( Figure 6D) and up to 89% accuracy for all eight speakers (averaged over all speakers). This outperformed a control dataset without the reservoir transformation, and very significantly for all speakers and few training utterances. The convergence with greater levels of training of the control classification with that using the reservoir seen in Figure 6D is similar to previous examples of reservoir computing being used in relatively simple classification problems. [28] The degree of convergence was also seen to be far less strong previously [28] when the initial data was more complex with less-effective preprocessing filtering applied. This gives confidence that Figure 6D demonstrates the simulated ring array was performing a meaningful RC transformation.
The heat maps in Figure 6E show the magnetic ring array reservoir performance with a single speaker for a range of values of H centre and k, and different levels of training. This demonstrates reservoir viability over a range of operating conditions of rotating magnetic field amplitudes drawn from the non-linear region of Figure 2D.

Discussion
DWs in patterned nanowires have shown hugely promising functionalities, for example, for shift register memory, [1] random access memory, [2,27] and logic; [3][4][5] but the stochastic nature of DW interactions with potential barriers has been a major obstacle to their development into commercial applications. Our experiments transform DW stochasticity from being a barrier to technical operation into the underpinning basis of emergent behavior, with reliable, field-dependent array magnetization, and DW population resulting from the probabilistic nature of DW de-pinning at ring junctions. Field rotation frequency will have an effect on the field strength of observed DW de-pinning [31] but we found similar emergent behavior across the limited range of field rotation frequencies used here (0.2-27 Hz, as appropriate for each experimental technique).
We found that emergent behavior in the nanoring arrays was observed with similar character in arrays of different width wires (2 and 4 µm) of different thickness (5 and 30 nm), although these differences affected the magnetic fields at which the minimum in DW population or magnetization appeared. This demonstrates that the array behavior is remarkably robust to changes in sample geometry, although sample design could be used in future to tune the emergent response to appear at convenient field strengths. For example, micromagnetic simulations of a single bi-ring junction (200 nm wire width, 5 nm wire thickness) for different junction geometries ( Figure S5, Supporting Information) show DW de-pinning fields from 30-80 Oe ( Figure S5C, Supporting Information), depending on the size of wire overlap at the junction. The simulations also showed pinned DWs to be drawn into the wire junction for small wire overlap (e.g., Figure S5D, Supporting Information for 10 nm overlap) but remain outside the wire junction until the de-pinning field is reached for larger wire overlaps (e.g., Figure S5F, Supporting Information). It is likely that these differences in magnetic configuration will also lead to differences in the energy barrier for DW de-pinning and thus on the emergent behavior of the array as a whole. Permalloy was chosen in this study due to its near-zero magnetocrystalline anisotropy and the relative ease with which it allows DW propagation. Use of other materials and an overall in-plane magnetocrystalline anisotropy (e.g., from using an in-plane magnetic field during deposition) is likely to create differences in the energy barriers associated with wire junctions along orthogonal directions. DWs in very close proximity, typically closer than a wall's width, have been observed to experience dipolar interactions. [32] Hence, although dipolar interactions will clearly be significant in the formation of the magnetic configurations at junctions, they are unlikely to be seen between DWs across a ring.
Our simulations of harnessing the arrays' emergent behavior to perform RC classification of spoken digits suggests that magnetic ring arrays are an interesting candidate system for in materio computation. The nanowire ring arrays show the characteristics required to perform in materio RC: representation of data (by magnetization or DW population); a nonlinear response characteristic (field-response of magnetization [ Figure 2D, Figure 4B] and DW population [ Figure 2E]); and fading memory to new input (seen here via MOKE measurements as the gradual approach to equilibrium from an initial state [ Figure 2F-H]). The ring system also offers repeatable fabrication and scalability (easily achieved with lithographic patterning), which most other in materio RC systems lack, for example, carbon nanotubes or waves on water. The square array arrangement used here was simply due to fabricational simplicity for investigation of emergent properties and was not optimized for RC uses, or even designed originally with RC in mind. Future work exploring different, perhaps non-uniform, ring arrangements is likely to result in improved RC performance.
The further requirements of an ideal in materio RC system of electronic integration of input/output also appear to be achievable: the continuous nature of nanowires means that magnetoresistance effects could provide CMOS-compatible output; and the nanowire junctions could be addressed (e.g., with current-induced field) to tune the local DW pinning probability (the effect of differences in pinning is seen in the magnetic configurations within rings at the edge or center of arrays in Figures 4C, 5B). Few other candidates for in materio RC hardware offer this breadth of ideal characteristics. The demonstration here of spoken digit recognition offers encouragement for future experimental RC demonstrations with magnetic ring arrays and detailed characterization of magnetic ring array RC performance. These future experiments may highlight unforeseen limitations in the description offered by the phenomenological model but the model's ability to predict experimental results seen elsewhere in this article give sufficient confidence that magnetic rings array RC holds promise.

Conclusion
We have demonstrated that the stochastic nature of magnetic domain wall (DW) propagation in arrays of interconnected soft magnetic nanorings when driven by rotating magnetic fields leads to emergent behavior being observed in the arrays, and that this behavior can be used to perform reservoir computing (RC). The stochasticity arises from the probability of moving DWs overcoming the potential barriers presented by ring junctions. DW interactions mean DW pinning and de-pinning create DW loss and gain mechanisms. The interconnected nature of arrays leads to dynamic equilibria being established in array magnetization and DW population that are robust nonlinear functions of field strength. By way of example, a phenomenological model demonstrated use of the nonlinear emergent behavior of magnetic nanoring arrays to perform RC of voice recognition tasks highly effectively. There are many other suggestions for in materio RC systems based on nonlinear response but the magnetic ring arrays appear to be particularly appealing as they offer the other features required of an ideal RC platform: clear representation of data, fading memory, scalable, and repeatable fabrication, and proven electronic integration. Future work on the sensitivity of the ring arrays' emergent properties to wider experimental parameters (e.g., ring dimensions) will be helpful to understand how flexible they may be in applications such as RC.

Experimental Section
Samples were fabricated using electron beam lithography with patterns defined using a RAITH Voyager system. A positive resist was spin-coated onto a Si (001) wafer substrate with a native oxide present. Metallization using permalloy (Ni 80 Fe 20 ) powder was performed in a custom-built (Wordentec Ltd) thermal evaporator (base pressure < 10 -7 mbar) before lift-off completed the lithography process. Samples for PNR and MOKE measurements covered 2 cm x 2 cm in blocks of 26×26 rings in squarelattice arrays, with rings 4 µm in diameter, 400 nm in wire width and 20 nm thick, and with 50% (200 nm) overlap between neighboring rings. Samples for X-PEEM were in various sized square arrays of up to 25 × 25 rings. The rings were made from 5 nm permalloy covered with 2 nm Al, had 50% overlap between neighboring ring wires and had either 200 nm wire widths with 2 µm diameter rings, or 400 nm wire widths with 4 µm diameter rings.
Polarized neutron reflectometry (PNR) data were acquired using the Offspec beamline at the ISIS neutron and muon source (Rutherford Appleton Laboratory, U.K.). The incident beam was polarized so that the neutron spin was aligned either parallel or antiparallel to the magnetic field, which was applied in plane, perpendicular to the direction of travel of the neutrons and also defines the neutron quantization axis. For each sample state, measurements for both incident spin-states, "up" and "down", were acquired but no polarization analysis was performed. As a consequence, the measured reflectivities R+ and R-consisted of sums of the non-spin flip and spin-flip components for a given incident spin-state and were sensitive to the magnetization along the neutron quantization axis. In analyzing the data both the scattering from the magnetic rings and the exposed silicon had to be taken account. The second contribution further splits into contributions from the exposed areas within and between the rings and the contribution from the surrounding Si wafer, away from the patterned area. The in-plane coherence length for neutron reflectometry is of the order of 100 µm [33] which is significantly longer than the relevant length scales for the contributions of the rings and the Si in and between them but not longer than the length scale, of the order of mm, of the surrounding Si on the sample wafer. Our model therefore coherently summed the ring and inter-ring signals and then incoherently summed those with the signal of the surrounding Si wafer. All model fitting was carried out using the GenX package. [34] As the rings were individually too small to resolve their magnetization configuration with specular neutron reflectometry, each ring contributed its net magnetization projected onto the applied field axis to the overall signal. A vortex state would therefore always contribute zero magnetization, and an onion state aligned with the magnetic field will contribute the maximum possible magnetization achievable at moderate applied fields. As the neutrons flooded the entire sample, of the order of cm 2 , it was possible to obtain information of the collective state of the array. All measurements were performed with an in-plane magnetic field of 18 Oe as a compromise to maintain neutron polarization but without affecting sample magnetization significantly.
Magneto-optical Kerr effect (MOKE) magnetometry was performed on a system similar to one described elsewhere. [35] Briefly, light from a 532-nm-wavelength continuous-wave laser (Vortran Stradus 532-40) was expanded to ≈8-mm diameter, passed through a Glan-Laser prism (giving p-polarized light at the sample) and directed onto a sample at a 45° angle of incidence. A half-wave retardation plate was used to remove ellipticity in the reflected beam before an analyzing Glan-Laser polarizer was used to convert Kerr rotation in the optical polarization into intensity change, which was detected using a Si photodiode. A quadrupolar electromagnet around the sample allowed computerdefined magnetic fields to be applied to the sample. The time-dependent optical signal and magnetic field information were recorded on a digital oscilloscope (Agilent infiniium 54832D MSO) before being sent to a personal computer for analysis.
X-ray photo-emission electron microscopy (X-PEEM) was performed at beamline I06 at the Diamond Light Source. X-PEEM images were obtained with an Elmitec SPELEEM-III microscope on beamline I06 at Diamond Light Source. I06 is a soft X-ray beamline with two APPLE-II undulators that allow full control of the light polarization. The magnetic domains images were obtained by averaging a series of X-ray absorption (XAS) images on and off the Fe-L 3 resonance and with left and right X-ray circular polarization in order to generate contrast by X-ray magnetic circular dichroism (XMCD). Samples were mounted on cartridges with a quadrupole magnet, to provide an in-plane magnetic field with arbitrary direction, designed at the CIRCE XPEEM beamline at the Alba Synchrotron. [36] Orthogonal field directions were each driven by a current power supply and used to generate in-plane rotating magnetic field at a frequency 8 Hz. An initial saturating field pulse was orthogonal to the measurement sensitivity direction, hence, the images for 10 and 42 Oe field amplitude in Figure 4A both show saturation but in orthogonal directions. Analysis of images from intermediate field strengths required removal of regions that showed the initial saturation before interpreting the changes in magnetization.
Magnetic force microscopy (MFM) was performed using a Digital Instruments Veeco Multimode instrument with Bruker MESP-LM-V2 CoCr-coated tips. Magnetic fields were applied before samples were placed in the MFM system using the electromagnet from the MOKE magnetometer.
Micromagnetic modeling was performed using the open-source MuMax3 [37] finite difference package to solve the Landau-Lifshitz-Gilbert equation of motion for dynamic magnetic systems. The two-ring system in Figure 1 (and images for Figure S3, Supporting Information) was described using 4 nm (width) × 4 nm (length) × 20 nm (thickness) cells. The parameters used here were those for permalloy: saturation magnetization, M s = 715 × 10 3 A m −1 (based on ferromagnetic resonance measurements of our films), exchange stiffness constant, A = 13 × 10 −12 J m −1 . The Gilbert damping constant, α, was set to one to allow rapid equilibration of magnetization states and allow the dynamics of relatively slow (Hz-scale) field rotations to be simulated. Magnetic field was rotated in angular steps of 15°. Each step was held for 12-15 ns to allow equilibration to occur. Simulation of MFM images from micromagnetic images was performed using the built-in MFM MuMax3 function. [37] The half bi-ring system in Figure S5, Supporting Information was modeled using cell sizes of 4 nm (width) × 4 nm (length) × 5 nm (thickness), α = 0.5 and other constants to represent permalloy. The geometry of the sample design is shown in Figure S5A, Supporting Information. Magnetization was allowed to relax from an initially saturated state in the direction shown in Figure S5B DW in two of the four wire arms. A magnetic field was applied in the saturation magnetization direction during the relaxation process and was then rotated anti-clockwise by 5° every 12 ns of simulation time. The magnetic field strength was varied across successive simulations in steps of 10 Oe until the DWs were able to pass through the wire junction.
Phenomenological modeling of square-lattice ring arrays with arbitrary numbers of rings recorded the position and motion of DWs through rings, with each divided into 16 sections of equal length. A field vector was rotated in 22.5° steps (i.e., 16 steps per revolution) and DWs allowed to move in order to align neighboring domains to the new field direction, unless a junction was encountered. In this case, a probability, P, for DWs to remain pinned was calculated using the approach by Wernsdorfer et al. [38] by P (t) = e −t/τ , where t is time and τ is the average time of de-pinning. The latter is given by τ (T, H) = (1/f 0 ) exp(E(H)/k B T), where f 0 is the attempt frequency ( =10 GHz here), E is the magnetic-field (H)-dependent energy barrier, T is temperature and k B is Boltzmann's constant. The energy barrier is calculated as E (H) = E 0 (1 − H/H 0 ) α , where E 0 is the zero-field energy barrier, H 0 is a characteristic switching field, and α ( =1.5) is an analytically-derived exponent. [39] E 0 and H 0 were fitted to the PNR data, as described in the Results section.
The simulations of performing RC with magnetic ring arrays used the phenomenological model with different de-pinning probabilities for one or two DWs at a junction. The array was treated as the reservoir; data were input by modulating the amplitude of a cyclic magnetic field and the reservoir response was measured as the normalized DW population. We used the method in ref. [40] to perform the RC simulation, and terms below are taken from this source. 26 vocal utterances of each digit from zero to nine by different speakers [30] were used as input datasets. These were separated into "training", "validation", and "test" sets: the training set consisted of N train utterances of each digit used to calculate a weight matrix W out ; the validation set used ten utterances of each digit used to train a regularization factor, the hyperparameter λ (see below); The remaining utterances were used as the test dataset to measure the model's classification performance. Pre-processing was used to transform the raw audio waveforms ( Figure 6A, Figure S4A, Supporting Information) into Mel-frequency cepstrum coefficients ( Figure S4B, Supporting Information), as is typical for speech recognition tasks. [41,42] The frequency response for the 13 Mel-frequency bands in each time window τ w were concatenated to produce an input vector u(t) for each utterance ( Figure S4C, Supporting Information). Each datum u i from u(t) was transformed into a single rotation of a magnetic field with amplitude H i input given by where k is an input scaling factor and H center is an offset field strength.
The reservoir response x i (DW population) for each input was measured, forming reservoir state vector x(t) ( Figure S4D, Supporting Information). This vector was then averaged to produce 52 features in output vector x t ( ) for each utterance ( Figure 6B, Figure S4E). Each utterance had a target output vector y j , where j is an integer from 0-9; y j was zero for all values of j except for the value corresponding to the digit, when it was equal to unity, that is , for an utterance of the digit "7", y 7 = 1 and y j ≠ 7 = 0. The matrix connecting the averaged reservoir state to the target vector was calculated using ridge regression, and W out calculated by where X and Y denote matrices made up of all the vectors x and y in the training sets, respectively, X T is the transpose of matrix X, and I the identity matrix. Multiplication of vector x of a given utterance in the "test" set with the trained weight matrix W out gave a vector of output activation nodes. The node with the highest activation was taken as the winner and provided the predicted class of the utterance ( Figure 6C). If the prediction matched the label, the reservoir had correctly classified the spoken digit. The input data was also trained against without the reservoir transformation ( Figure S4F, Supporting Information) to provide a control measurement.

Supporting Information
Supporting Information is available from the Wiley Online Library or from the author.