Perpendicular-anisotropy artificial spin ice with spontaneous ordering: a platform for reservoir computing with flexible timescales

TL;DR

Researchers created nanomagnet arrays with perpendicular anisotropy that self-organize at room temperature, offering a new platform for computing with adjustable dynamics.

Contribution

First demonstration of strongly coupled 2D out-of-plane Ising spin lattices with spontaneous ordering at room temperature.

Findings

Nanomagnet arrays with perpendicular anisotropy spontaneously order at room temperature.

System enables tunable temporal dynamics suitable for reservoir computing.

Strong magnetostatic coupling opens new possibilities for artificial spin ice geometries.

Abstract

Arrays of coupled nanomagnets have wide-ranging fundamental and practical applications in artificial spin ices, reservoir computing and spintronics. However, lacking in these fields are nanomagnets with perpendicular magnetic anisotropy with sufficient magnetostatic interaction. This would not only open up unexplored possibilities for artificial spin ice geometries but also enable novel coupling methods for applications. Here, we demonstrate a method to engineer the energy landscape of artificial spin lattices with perpendicular magnetic anisotropy. With this, we are able to realize for the first time strongly magnetostatically-coupled 2D lattices of out-of-plane Ising spins that spontaneously order at room temperature on timescales that can be precisely engineered. We show how this property, together with straightforward electrical interfacing, make this system a promising platform for reservoir computing. Our results open the way to investigate the thermodynamics of out-of-plane magnetostatically coupled nanomagnet arrays with novel spin ice geometries, as well as to exploit such nanomagnet arrays in unconventional computing, taking advantage of the adjustable temporal dynamics and strong coupling between nanomagnets.

Aleksandr Kurenkov and colleagues demonstrate nanomagnet arrays with out-of-plane anisotropy that spontaneously order at room temperature. This system offers a new platform for unconventional computing through tunable dynamics.

Figures

Dependence of the ordering timescale on the energy balance in a 2D Ising square lattice.a Schematic of the system with 11 × 11 spins. It is initialized to a uniform magnetic state with all spins pointing up and then released at *t* = 0. The lattice then relaxes to lower energies for *t* > 0. This system size was chosen because it matched the size of our fabricated systems discussed later. b Evolution of the local antiferromagnetic parameter *q*_NN_ and the average magnetization *m*_avg_ for different effective anisotropy *E*_EA_ and magnetostatic coupling *E*_MC_. 20 simulations (fine lines) a

Phase diagrams of the average magnetization *m*_avg_ and the local antiferromagnetic parameter *q*_NN_ as a function of effective anisotropy *E*_EA_ and magnetostatic coupling *E*_MC_ at *t* = 1000 s.Five regions can be distinguished: I – frozen state, II – transient states, III – state with domains and domain boundaries, IV – checkerboard ordering, V – superparamagnetic state. The white (black) dotted line corresponds to *m*_avg_ = 0.05 (*q*_NN_ = 0.95). Labels 3c-3f correspond to the energy landscapes shown in Fig. [3](#Fig3)c–f. The three symbols correspond to the experimental lattices with

Dependence of the energy landscape of multilayered nanomagnets with perpendicular magnetic anisotropy on the material stack and lateral dimensions.a Magnitude of the magnetic flux density for a permalloy nanomagnet (dimensions: 100 × 300 × 20 nm^3^) with in-plane anisotropy (left) and 200 nm-diameter nanomagnets with perpendicular magnetic anisotropy consisting of 3 Co layers of 0.5 nm thickness (centre) and 8 Co layers of 1.45 nm thickness (right). The Co layers are separated by 0.8 nm of vacuum. b Schematic showing the geometrical parameters used in energy landscape calculations. c–g Energy

Observing the dependence of the ordering timescale on nanomagnet separation and nanomagnet shape with MFM.All systems were initialized to a uniform magnetic state with an out-of-plane magnetic field before the measurements. a MFM images of lattices with *t*_Co_ = 1.45 nm, *D*_NM_ = 170 nm and increasing nanomagnet separation *S*_ASI_ observed at *t* ~ 1000 s after initialization. Black nanomagnets are in the initial magnetic state while white nanomagnets have switched. Image frame colours match those of the corresponding curves in d and frames in e. b Array of 170 nm-wide square nanomagnets se

Dependence of ordering on the Co layer thickness and the nanomagnet size observed with MFM.All systems were initialized to a uniform magnetic state (black contrast) with an out-of-plane magnetic field before the measurements. a MFM images of lattices with *D*_NM_ = 200 nm, *S*_ASI_ = 20 nm and increasing Co thickness *t*_Co_ observed at *t* ~ 1000 s. Image frame colours correspond to those of the curves in b. White contrast indicates that the nanomagnets have switched. b Energy landscapes calculated for the lattices shown in a. c Average magnetization *m*_avg_ and local antiferromagnetic param

Monte Carlo simulations of the performance of an 11 × 11 lattice reservoir for a sine to sawtooth signal transformation.The parameters used for the simulations are *E*_MC_ = 2.5 *k*_B_*T*, *E*_EA_ = 20 *k*_B_*T* ± 5% and *B*_0_ = −0.35 mT. a Schematic of the lattice used as a reservoir. Blue readout lines indicate the row-by-row averaged magnetic state readout. b Inverse mean squared error 1/MSE as a function of the input magnetic field frequency and amplitude *B*_1_. Higher (lower) values in yellow (purple) indicate better (worse) performance. The black contours correspond to 1/MSE ~ 30, high

Full text

Introduction

Coupled nanomagnets organized on the sites of various lattices are widely-explored in the field of artificial spin ice^1–3^ since they exhibit a large variety of fascinating phenomena including collective dynamics^4–6^, frustration^7–9^, dynamic chirality^10^ and phase transitions^11–13^. Furthermore, these properties can be exploited for novel forms of conventional^14–19^ and unconventional^20–24^ computing. The single-domain nanomagnets typically have an in-plane magnetic anisotropy, which results in strong magnetostatic coupling between the magnets due to extended demagnetizing fields in the lattice plane. In contrast, dipolar-coupled nanomagnets with perpendicular magnetic anisotropy have rarely been employed for artificial spin ice because they have not shown spontaneous ordering–so they are not “thermally-active”–at experimentally measurable timescales^25–30^. This lack of spontaneous evolution makes such magnets a common choice for high-density bit-patterned media^31,32^.

Creating lattices comprised of magnetostatically coupled nanomagnets with perpendicular anisotropy that show spontaneous ordering on an experimentally-accessible timescale would open the way to study the thermodynamics of out-of-plane artificial spin ices based on a large variety of lattices. Among them are the lattices associated with the crystal planes of many bulk magnetic systems^33–35^, the triangular lattice^36,37^, and the canonical two-dimensional (2D) Ising model on the square lattice^38–40^. While any of these lattices can be created with the method described in this work, we will focus on the latter due to its broad importance and universality. Being universally complete, so that any other statistical model^41,42^ or Boolean circuit^43,44^ can be derived from it, the 2D Ising model is of fundamental importance for statistical physics and has been used to model numerous physical, mathematical and biological processes^45–47^. The modelling is typically carried out by selecting a particular variant of the model, setting an initial state, and then observing how it evolves into a lower-energy state over time. The nature of this evolution provides a means to solve non-deterministic polynomial-time (NP-) hard problems with only polynomial overhead (P-hard) by mapping them onto a corresponding Ising system^48,49^. Since any NP-hard problem can be formulated as an Ising problem, the 2D Ising system is of great interest for a variety of combinatorial problems^46^. In addition, energy minimization and state dynamics in 2D Ising lattices bear similarities to equivalent processes in the human brain^50^, providing a basis for computational models of the brain^50,51^ and several types of neural networks such as Boltzmann machines^52^ or Hopfield networks^53^.

Here, we create for the first time a lattice of magnets with out-of-plane anisotropy that spontaneously orders at room temperature. We focus on a 2D square lattice of Ising spins (Fig. 1a) because of its importance for fundamental science and for applications. By fabricating such lattices with nanomagnets with various diameters and separations, and made from precisely tuned multilayer films with different number and thicknesses of the layers, we have been able to determine both theoretically and experimentally how the energy landscape of the system defines its ordering dynamics. Utilizing this information, we have engineered the timescale of the system response, which provides a means to tailor the arrays for specific computing applications.Fig. 1. Dependence of the ordering timescale on the energy balance in a 2D Ising square lattice.a Schematic of the system with 11 × 11 spins. It is initialized to a uniform magnetic state with all spins pointing up and then released at t = 0. The lattice then relaxes to lower energies for t > 0. This system size was chosen because it matched the size of our fabricated systems discussed later. b Evolution of the local antiferromagnetic parameter qNN and the average magnetization mavg for different effective anisotropy EEA and magnetostatic coupling EMC. 20 simulations (fine lines) and their average (bold line) are shown for each case. The green-shaded region highlights the second stage of ordering where the slope of qNN decreases and involves the reversal of spins within the domain boundaries. The simulations are stopped after 1000 s and the resulting magnetic states are shown in c. c Magnetic states at t = 1000 s. Here white (black) contrast corresponds to the up (down) spins. Panels 1–9 correspond to those of b.

Results

Relating the energy landscape to the magnetic relaxation timescale with Monte Carlo simulations

We begin with a theoretical exploration of the degree and timescale of spontaneous ordering in a square 2D Ising system and how it is influenced by the energy landscape of the spins. For this, we performed kinetic Monte Carlo simulations using the “Hotspice” package^54^. The system simulated was a lattice of 11 × 11 magnetostatically (dipolarly, in this case) coupled Ising spins (Fig. 1a) in contact with a heat bath of temperature T. The lattice was first initialized in a state with all spins pointing up (upper panel of Fig. 1a) and allowed to relax for 1000 seconds. This time was chosen to match the approximate time between the initialization and measurement of the experimental system discussed later, with the average magnetization mavg and local antiferromagnetic parameter qNN tracked during this period. We define qNN as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(1-\langle {S}_{i}{S}_{i+1}\rangle )/2$$\end{document} , where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\langle {S}_{i}{S}_{i+1}\rangle$$\end{document} is the nearest-neighbour correlation, and the mean switching time of a spin j at a time t is given by the Néel-Arrhenius law^55–57^:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\tau }_{j}(t)={\tau }_{0} \exp \left(\frac{\Delta {E}_{j}(t)}{{k}_{B}T}\right)$$\end{document}

Here t is elapsed time, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\tau }_{0}={10}^{-10}{{\rm{s}}}$$\end{document} is the attempt period and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta {E}_{j}(t)$$\end{document} is the energy barrier

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta {E}_{{\mbox{j}}}\left(t\right)={E}_{{{\rm{EA}}}}\left(1+{\mathrm{rand}}\left(j\right)\right)+{E}_{{{\rm{MC}}}}{\sum}_{i,j}{S}_{i}{\left(t\right){D}_{{ij}}S}_{j}(t),$$\end{document}

defined by a configuration of spins \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${S}_{i}$$\end{document} at every time step, magnitude of the point dipolar interaction \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${D}_{{ij}}$$\end{document} between spins i and j, a random ±5% Gaussian variation of anisotropy \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathrm{rand}}(j)$$\end{document} for different spins^24,58^, as well as a time-independent effective anisotropy energy EEA and energy associated with the dipolar coupling between the nearest neighbour spins EMC (see Supplementary Note 1 for details).

The initial switching of individual nanomagnets is promoted by dipolar interactions with other nanomagnets in the lattice given by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${E}_{{{\rm{MC}}}}{\sum}_{i,j}{S}_{i}{\left(t\right){D}_{{ij}}S}_{j}(t)$$\end{document} but, as the system assumes a more ordered state, this term decreases in value. This results in an exponential increase of the switching time due to the exponential term in Eq. 1. Simulations carried out for several sets of EEA and EMC confirm the universally logarithmic dependence of mavg and qNN on the elapsed time (Fig. 1b). An increase in EMC (see columns in Fig. 1b) extends the relaxation dynamics over a longer timescale, while changing EEA (see rows in Fig. 1b) does not change the slope of the curves and only shifts them along the time axis. Irrespective of the EEA and EMC, the spin lattices achieve mavg ~ 0 (blue traces in Fig. 1b) earlier than perfect checkerboard ordering with qNN = 1 (orange traces in Fig. 1b). Furthermore, as the system nears mavg ~ 0, the ordering rate decreases, which can be seen by the decrease in the slope of the orange lines as they enter the green-shaded regions in Fig. 1b. Snapshots of the magnetic state at t = 1000 s (Fig. 1c) provide a clue to why this two-stage ordering process occurs (reflected by the two different slopes in qNN). Panels 1 and 2 in Fig. 1c are snapshots of the system during the rapid change of both mavg and qNN in the first phase of the ordering. Here, individual nanomagnets throughout the system switch (Panel 1) followed by domains of antiferromagnetic ordering starting to form (Panel 2). These doubly degenerate domains expand to completely fill the system (Panel 3 and 4) with boundaries forming between the domains with nanomagnets of opposite polarity (Panels 5 and 6).

From this point on, a further increase in qNN can only be achieved by switching the nanomagnets in the domain boundaries. As these magnets are now stabilized through the dipolar interactions with their neighbours that have already switched, the energy barrier for them to switch is considerably higher than at the beginning of the relaxation process. This increases the switching time exponentially (see Eq. 1) and slows down further ordering. The details of this process are elaborated in Supplementary Note 2. A consequence of the two-stage ordering is that, while a stronger coupling EMC promotes faster demagnetization mavg → 0 (e.g., see blue lines in Panels 7–9 of Fig. 1b, where the time to get to mavg ~ 0 decreases from ~10^-2^ to ~10^-3^), it does not necessarily result in faster ordering qNN → 1 (orange lines in Panels 7–9 of Fig. 1b) as it increases the stabilization of the magnetic states in the domain boundaries. Therefore, EMC must not exceed a certain upper limit for the system to order on a given timescale and to not dwell in a metastable state with mavg ~ 0 and qNN < 1.

The dependence of mavg and qNN on EEA and EMC at t = 1000 s is summarized in the phase diagrams of Fig. 2. In Region I, neither mavg nor qNN experience a notable change compared to the initial state due to the strong anisotropy and weak magnetostatic coupling (corresponding to Panel 1 in Fig. 1c). Increase of the EMC/EEA ratio shifts the system through the transient Region II, characterised by decreasing average magnetization and increasing ordering (corresponding to Panel 2 in Fig. 1c). Further increase of EMC/EEA shifts the system to Region III with mavg ~ 0 and qNN still not 1, characteristic for the second ordering stage (corresponding to Panel 6 of Fig. 1c). Any further increase of EMC/EEA does not affect the average magnetization or ordering at t = 1000 s. Achieving the high ordering of Region IV, given by the dotted black line that envelops a region with qNN ≥ 0.95 in Fig. 2 (exemplified by Panels 8 and 9 in Fig. 1c), requires lowering both EMC and EEA below a certain threshold. The narrow Region V (enclosed by vertical box in Fig. 2) is the regime in which the spins are superparamagnetic at the chosen timescale. The border between Regions V and I at EMC = 0 is therefore a transition between the superparamagnetic and frozen regimes. This border is located at EEA ~ 30 kBT because nanomagnets with such an effective anisotropy have an average switching time of ~1000 s (the time point of these diagrams). Region IV with qNN ~ 1 is, therefore, constrained from the left by the superparamagnetic regime, from the top by the frozen regime, and from the right by states in which switching spins in domain boundaries is too energetically expensive.Fig. 2. Phase diagrams of the average magnetization mavg and the local antiferromagnetic parameter qNN as a function of effective anisotropy EEA and magnetostatic coupling EMC at t = 1000 s.Five regions can be distinguished: I – frozen state, II – transient states, III – state with domains and domain boundaries, IV – checkerboard ordering, V – superparamagnetic state. The white (black) dotted line corresponds to mavg = 0.05 (qNN = 0.95). Labels 3c-3f correspond to the energy landscapes shown in Fig. 3c–f. The three symbols correspond to the experimental lattices with nanomagnet diameter DNM = 170 nm, Co thickness tCo = 1.45 nm and nanomagnet separation SASI = 20 nm (circle), SASI = 25 nm (triangle), SASI = 30 nm (cross). The experimental magnetic configurations of these systems are shown in Fig. 4.

The phase diagrams in Fig. 2 are snapshots at t = 1000 s and, with increasing time, the boundaries between the different regions will shift. Because of this, if the observation time is long enough and EMC is non-negligible, a system located in the frozen region at t = 1000 s may eventually find itself in the more ordered Regions II, III or IV (see phase diagrams for t ~ 27 months in Supplementary Note 3). The state of the nanomagnet array is therefore defined by EEA, EMC and t. Practical applications exploiting systems with specific dynamics defined by EEA and EMC require careful adjustment of the lattice and nanomagnet dimensions, as well as the stack materials and layer thicknesses. We show how this can be achieved in the next section.

Designing the energy landscape of coupled Co/Pt nanomagnets

Artificial spin ices with perpendicular magnetic anisotropy that are thermally-active at room temperature on experimentally measurable timescales have not been realized before. The first reason for this is that achieving large EMC for this geometry is harder due to the confinement of demagnetizing fields in the vicinity of the nanomagnet (central magnet in Fig. 3a) compared to the in-plane case (leftmost magnet in Fig. 3a). The second reason is that balancing the energy contributions to EEA, required to lower the energy barrier to switching, is equally challenging in Co/x (x = Pt, Pd, Ni) multilayers. To address these challenges, we focused on precisely controlling EEA in perpendicular nanomagnets and maximizing EMC between them.Fig. 3. Dependence of the energy landscape of multilayered nanomagnets with perpendicular magnetic anisotropy on the material stack and lateral dimensions.a Magnitude of the magnetic flux density for a permalloy nanomagnet (dimensions: 100 × 300 × 20 nm^3^) with in-plane anisotropy (left) and 200 nm-diameter nanomagnets with perpendicular magnetic anisotropy consisting of 3 Co layers of 0.5 nm thickness (centre) and 8 Co layers of 1.45 nm thickness (right). The Co layers are separated by 0.8 nm of vacuum. b Schematic showing the geometrical parameters used in energy landscape calculations. c–g Energy landscapes calculated for the system shown in b for DNM = 200 nm, SASI = 20 nm and different Co layer thicknesses. The vertical scale is energy, where kB is the Boltzmann constant and *T *= 300 K.

We fabricated the nanomagnets from Co/Pt ferromagnetic multilayers, which are widely used because of their strong perpendicular magnetic anisotropy^59^. The primarily interfacial origin of the anisotropy allowed us to vary EEA and EMC almost independently by changing the number of Co/Pt interfaces and the thickness of the Co layers, respectively. We have calculated the energy landscape of a pair of nanomagnets to guide the selection of nanomagnet diameter DNM, separation SASI, number of Co layer repetitions NCo, and thicknesses of Co and Pt layers, tCo and tPt, respectively (Fig. 3b). We allowed the magnetization to rotate coherently in one nanomagnet ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\vec{M}}_{{{\rm{\theta }}}}$$\end{document} in Fig. 3b) while keeping the other fixed ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\vec{M}}_{{{\rm{fixed}}}}$$\end{document} in Fig. 3b). The energy was then calculated for polar angles θ from 0° to 180° and included four energy terms associated with the uniaxial interfacial anisotropy, demagnetization, magnetostatic interaction between the layers in the nanomagnet and the dipolar coupling between the nanomagnets (see Supplementary Note 4 for details). The resulting energy landscape provides an estimate of the energy barrier that the system needs to overcome to go from the higher-energy magnetic state with parallel moments to the lower-energy antiparallel state.

In Fig. 3c, we show the energy landscape for a pair of coupled nanomagnets made of (Co [0.5] / Pt [1.0])2 / Co [0.5] layers with NCo = 3, DNM = 200 nm and SASI = 20 nm, where the numbers in square brackets are thicknesses in nm. Such stacks have a high anisotropy and are widely used for spintronics applications^60,61^ due to the high thermal stability, which provides more than 10 years of retention time (a measure of how long a device can store information reliably without the nanomagnet switching) even in sub-20 nm diameter nanomagnets^60^. The small number of ferromagnetic Co layers, as well as their low thickness, results in a localized demagnetizing field (central magnet in Fig. 3a), minimizing crosstalk between the nanomagnets. Both properties–high anisotropy and low crosstalk–while useful for information storage applications, contradict the high-EMC, low-EEA requirements of a lattice to give fast spontaneous ordering. Arrays of nanomagnets made of such a stack are located in the “frozen” Region I of Fig. 2 (“3c” label).

To devise a high-EMC, low-EEA stack, we varied NCo and tCo while fixing the nanomagnet diameter DNM and separation SASI. SASI should be minimized to maximize the dipolar coupling, and was set to 20 nm, which was the smallest nanomagnet separation we were able to obtain with confidence when fabricating nanomagnet arrays with electron beam lithography. DNM was set to 200 nm, which was the largest possible nanomagnet diameter that did not result in the formation of multidomain states. This upper threshold value for DNM was determined experimentally by observing the magnetic states in multiple lattices of nanomagnets fabricated with different diameters.

Increasing the number of Co layer repetitions NCo from 3 to 8 produces a stack similar to those used in previous works on artificial spin ices with out-of-plane nanomagnets^25–29^, with the energy landscape shown in Fig. 3d. EMC is approximately quadratically proportional to the total Co thickness and increases by a factor of ~7. However, the increased number of Co/Pt interfaces results in a higher anisotropy, as seen in Fig. 3c, d, and the system remains in Region I of Fig. 2 (“3d” label). Achieving an EMC-dominated lattice requires the lowering of EEA without decreasing EMC. For this, one needs to consider the two additional components of EEA beyond the interfacial anisotropy: the demagnetization energy and dipolar coupling between layers of the stack. The dipolar coupling between the layers is rather weak and increases the perpendicular anisotropy, and thus cannot be used to lower EEA. In contrast, the demagnetization energy is considerable and decreases EEA. The demagnetization energy increases with Co thickness, just like EMC, which has the same physical origin. Therefore, increasing tCo to 1.25 nm results in a reduction in EEA and an increase in EMC (Fig. 3e), shifting the system towards Region II of Fig. 2 (“3e” label). At tCo = 1.45 nm (Fig. 3f), EMC ~ EEA, which results in a highly asymmetric landscape with the lowest EEA and highest EMC among the considered stacks. Such an energy landscape should facilitate spontaneous ordering and, depending on the absolute values of EEA and EMC, an array of such nanomagnets belongs to Region III or IV of Fig. 2 (“3f” label). The demagnetizing field of these nanomagnets (rightmost nanomagnet in Fig. 3a) has a similar extent to that of the Permalloy nanomagnets typically implemented for in-plane artificial spin ices (leftmost nanomagnet in Fig. 3a). Further increase of tCo to 1.65 nm results in the shift of the energy minimum to θ ~ 110°, indicating a loss of perpendicular anisotropy (Fig. 3g).

From calculations, we have shown the effect on the energy landscape of changing the experimental system in terms of layer thicknesses, as well as the nanomagnet diameters and separations. By comparing this to the results from our Monte Carlo simulations, we are then able to predict how changing the system parameters will influence the relaxation timescale. This has enabled us to experimentally realize square lattices of interacting out-of-plane nanomagnets with the ability to control the ordering timescales from sub-seconds to years, as we show in the following section.

Control of relaxation timescales in experimental dipolar-coupled 2D Ising systems

We now turn to experimental systems where limitations in fabrication impose additional constraints on the system design. For example, achieving small SASI—a key parameter influencing EMC—becomes more challenging as the thickness of the stack increases, so we limit NCo to 7. Another effect is that the Pt spacer layers may become discontinuous on decreasing their thickness below ∼0.4 nm and, to be sure that we have a continuous layer we choose tPt = 0.8 nm for the experimental systems, unless otherwise mentioned. These parameters permitted an SASI down to 20 nm, and we varied SASI, DNM and tCo in the experiment.

The multilayers of Ta[3] / Pt[4] / (Co[tCo] / Pt[0.8])6 / Co[tCo] / Ru[2] were deposited by magnetron sputtering on Si substrates with a natural oxide layer. They were patterned into square lattices of 11×11 nanomagnets by electron beam lithography and Ar ion milling. For transport measurements, the Ta/Pt bottom layer was patterned into Hall bars. Magnetic force microscopy (MFM) was used to measure the magnetic state after applying an out-of-plane magnetic field of 0.9 T to initialize all nanomagnets to an “up state”. For the MFM measurements, we covered the samples with a 60 nm-thick Poly(methyl methacrylate) resist to minimize the influence of the stray magnetic field of the MFM probe on the sample.

We first looked at arrays of nanomagnets with tCo = 1.45 nm, DNM = 170 nm and different SASI. The state of the arrays at time t ~ 1000 s after initialization is shown in Fig. 4a. Black contrast indicates that the nanomagnets are in the state initialized by magnetic field, while white contrast indicates nanomagnets that have spontaneously switched. The number of switched nanomagnets gradually decreases with increase in SASI. Since SASI has no effect on EEA and only changes EMC, this result is easy to interpret. The decrease of EMC at a constant EEA lowers the asymmetry of the energy landscape by making the initial energy well deeper and therefore harder to leave (θ = 0° in Fig. 4d) while making the lower-energy well shallower (θ = 180° in Fig. 4d).Fig. 4. Observing the dependence of the ordering timescale on nanomagnet separation and nanomagnet shape with MFM.All systems were initialized to a uniform magnetic state with an out-of-plane magnetic field before the measurements. a MFM images of lattices with tCo = 1.45 nm, DNM = 170 nm and increasing nanomagnet separation SASI observed at t ~ 1000 s after initialization. Black nanomagnets are in the initial magnetic state while white nanomagnets have switched. Image frame colours match those of the corresponding curves in d and frames in e. b Array of 170 nm-wide square nanomagnets separated by 30 nm observed at t ~ 1000 s. tCo = 1.45 nm as in a. The nanomagnet shape is shown in the top left inset. All other results are for circular nanomagnets as shown in the top left inset of a. c The same nanomagnet arrays as in a observed after ~27 months at room temperature. The magnets that have switched during this time are highlighted with yellow circles. d Energy landscapes calculated for the lattices shown in a. e Evolution of average magnetization mavg and local antiferromagnetic parameters qNN, q2NN, q3NN. The fine lines are calculated using Monte Carlo simulations, the bold lines are their averages, and the points are experimental data. f Average magnetization mavg and local antiferromagnetic parameter qNN as a function of SASI determined from the MFM images. Data for a–c are represented as filled circles, open squares and semi-filled circles, respectively. MFM images for SASI = 35 and 40 nm can be found in Supplementary Note 5.

From the MFM images, we calculated the average magnetization mavg and local antiferromagnetic parameter qNN as defined above. With decreasing SASI, the monotonic increase in qNN and decrease in mavg (filled circles in Fig. 4f; MFM images for 35 and 40 nm are shown in Supplementary Note 5) highlight the important role that the magnetostatic coupling plays in the spontaneous ordering. We do not reach qNN = 1, mavg = 0, characteristic of highly ordered states since a further decrease of SASI is challenging in terms of the fabrication. Instead, we changed the shape of the nanomagnets from circular to square with rounded corners (compare insets in Fig. 4a, b) without modifying the stack. The resulting closer proximity of the magnets produces higher EMC and a substantial increase in ordering (MFM image in Fig. 4b; mavg and qNN shown as a square data point in Fig. 4f).

The mavg and qNN of a 2D Ising system with any EEA and EMC will evolve with time (see, for example, Fig. 1b). Therefore, one can achieve a higher degree of ordering simply by waiting long enough. To confirm this, we kept the nanomagnet arrays at room temperature in ambient atmosphere and no magnetic field for ~27 months, which is equivalent to ~7 × 10^7 ^s. We then remeasured the samples without applying any magnetic field and observed the states shown in Fig. 4c (additional intermediate magnetic states are shown in Supplementary Note 5). The nanomagnets that have additionally switched are highlighted by yellow circles. There was no switching of nanomagnets from white to black (i.e. back to the state initialized 27 months prior). The new mavg and qNN are shown in Fig. 4f with half-filled circles. The decrease in mavg and increase in qNN indicate the additional ordering of the lattice, and one can make use of the evolution of these parameters to pinpoint location of the lattices on the phase diagrams of Fig. 2. For this, we calculated mavg, qNN as well as local antiferromagnetic parameters q2NN and q3NN for the 2^nd^ and 3^rd^ nearest neighbours from the MFM images for t = 1000 s and t = 7 × 10^7 ^s. We then fitted these values with mavg(t), qNN(t), q2NN(t) and q3NN(t) calculated using Monte Carlo simulations with the “Hotspice” package^54^, and using EEA, EMC and the ferromagnetic exchange coupling between the magnets J as variables (see Supplementary Note 5 for details). The fitted results are shown in Fig. 4e. The non-zero J in the case of SASI = 20 nm (J ~ 6.5 mT) may indicate incomplete separation of the nanomagnets in the bottom Co layer. The near-zero J for SASI = 25 and 30 nm suggests that the magnets are fully separated. The location of the fitted results on the phase diagrams of Fig. 2 (circles, triangles and crosses for SASI = 20, 25 and 30 nm, respectively) indicates that the achieved levels of EMC are sufficient for a complete ordering (Region V) in t ~ 1000 s but EEA is too high. Note that these phase diagrams were calculated assuming that there is no exchange coupling between the nanomagnets, so J = 0. Therefore, since the SASI = 20 nm “circle” experimental point is for the lattice with J ~ 6.5 mT, its location is approximate.

Having looked at the effect of SASI and time t on ordering in the system, we then looked at the influence of tCo and DNM. MFM images taken at t ~ 1000 s of the nanomagnet arrays with DNM = 200 nm, SASI = 20 nm and varying tCo from 1.1 nm to 1.45 nm are shown in Fig. 5a. As we have seen, tCo has a profound effect on EEA and EMC (Fig. 3e–g), and consequently on the ordering timescale. Indeed, as tCo increases, the corresponding energy landscapes become more asymmetric (Fig. 5b) and, from the MFM images, we find that mavg decreases while qNN increases (Fig. 5c). Increasing DNM to 200 nm provides a lower mavg and higher qNN (0.06 and 0.67, respectively; see data for tCo = 1.45 nm in Fig. 5c) than for the same lattice with DNM = 170 nm (0.47 and 0.44, respectively; see data for SASI = 20 nm in Fig. 4f). However, further increase of DNM to 230 nm in the same stack results in formation of multidomain states within the nanomagnets, providing an upper limit to DNM. Conversely, a decrease of DNM to 140 nm slows down the ordering process so that there is hardly any switching at t ~ 1000 s (top images in Fig. 5d) and only a few reversed magnets at t ~ 7 × 10^7 ^s (bottom images in Fig. 5d). We have therefore experimentally confirmed that DNM and tCo have a profound effect on the ordering dynamics of the lattice with tCo requiring a precise adjustment to achieve a low enough EEA (as shown in Fig. 3d-g). The out-of-plane anisotropy of the nanomagnets means that the state of the nanomagnets can be directly accessed with electrical readout using the anomalous Hall effect. To demonstrate this, we performed electrical measurements on lattices with DNM = 140 and 170 nm and SASI = 30 and 40 nm on timescales of tens of seconds, which showed similar degrees of thermal activity to the lattices measured with MFM (Supplementary Note 6).Fig. 5. Dependence of ordering on the Co layer thickness and the nanomagnet size observed with MFM.All systems were initialized to a uniform magnetic state (black contrast) with an out-of-plane magnetic field before the measurements. a MFM images of lattices with DNM = 200 nm, SASI = 20 nm and increasing Co thickness tCo observed at t ~ 1000 s. Image frame colours correspond to those of the curves in b. White contrast indicates that the nanomagnets have switched. b Energy landscapes calculated for the lattices shown in a. c Average magnetization mavg and local antiferromagnetic parameter qNN as a function of tCo determined from MFM images in a. d MFM images of lattices with DNM = 140 nm and tCo = 1.45 nm observed at t ~ 1000 s and t ~ 7×10^7^ s. The nanomagnets that have switched between the two times are indicated with yellow circles.

We have now demonstrated experimentally that the relaxation dynamics of out-of-plane spin lattices can be engineered with careful tuning of the anisotropy of and magnetostatic coupling between the nanomagnets by altering their dimensions, separation and shape, as well as the stack layer thicknesses. We have provided numerical evidence that the timescale associated with the collective dynamics can be tuned all the way from sub-seconds to years. This opens up unexplored avenues for optimising the performance of reservoir computing systems incorporating these lattices as we show in the next section.

Tuneable-frequency reservoir computing with 2D Ising systems

The 2D Ising system is an appealing platform for a variety of computational approaches. Short-term memory and reset of the system after processing an input play a central role in brain-inspired computing^62,63^, and both of these properties can be achieved by exploiting the dynamics of mavg and qNN that we presented above. This ability to tune the dynamics is particularly important in reservoir and probabilistic computing where, in contrast to conventional von Neumann computing, higher frequencies are not inherently beneficial. Instead, the priority lies in matching the dynamics of the physical system to the timescale of incoming data patterns since this is the only way that the physics of the system can be effectively harnessed for data processing. To demonstrate the effect of frequency matching on the computational performance, we have simulated a signal transformation task using our 2D Ising system for reservoir computing.

Reservoir computing is a framework for neural networks that makes use of a system with non-linear behaviour called a “reservoir” to map input data into a higher-dimensional space, in which the inputs can be separated by a linear transformation. The reservoir does not need to be trained (so it is not modified itself) but does need to be a system with short-term memory and high dimensionality^62,64^, properties that are met by many physical systems^20,21,24,64–69^. Here we employ Monte Carlo simulations using the “Hotspice” package^54^ to simulate the performance of a reservoir comprising the 2D Ising system on transforming sine waves of different frequencies into a sawtooth signal, a task that is often used to test the nonlinearity of a reservoir^69,70^.

For the reservoir, we employed the 11×11 lattice with EEA = 20 kBT ± 5% and EMC = 2.5kBT (Panel 8 of Fig. 1b). We have experimentally demonstrated electrical readout of mavg in the 11×11 lattices using the anomalous Hall effect (see Supplementary Note 6), and therefore we implemented this magnetic state readout method in the simulations. The input for the simulations was applied in the form of an out-of-plane magnetic field B(t) acting on the entire system (Fig. 6a) and was scaled such that its magnitude extends from B0 (minimum) to B1 (maximum) as shown in the top panel of Fig. 6c. Accordingly, the energy barrier term of Eq. 2 was extended to include the Zeeman term associated with the applied magnetic field as follows:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta {E}_{{\mbox{j}}}\left(t\right)={E}_{{{\rm{EA}}}}\left(1+{{\rm{rand}}}\left(j\right)\right)+{E}_{{{\rm{MC}}}}{\sum}_{i,j}{S}_{i}{\left(t\right){D}_{{ij}}S}_{j}(t)+{E}_{{{\rm{Zeeman}}}}\left(B\left(t\right),j\right)$$\end{document}

Fig. 6. Monte Carlo simulations of the performance of an 11 × 11 lattice reservoir for a sine to sawtooth signal transformation.The parameters used for the simulations are EMC = 2.5 kBT, EEA = 20 kBT ± 5% and B0 = −0.35 mT. a Schematic of the lattice used as a reservoir. Blue readout lines indicate the row-by-row averaged magnetic state readout. b Inverse mean squared error 1/MSE as a function of the input magnetic field frequency and amplitude B1. Higher (lower) values in yellow (purple) indicate better (worse) performance. The black contours correspond to 1/MSE ~ 30, highlighting regions where the reservoir performs better than the linear transformation of the original input signal. c Temporal view of the transformation. The upper panel shows the input sinewave signal (black dotted curve) and the target sawtooth signal (black dashed line). Shown in the lower three panels are the prediction with and without the artificial spin ice (ASI) reservoir (blue and red curve, respectively) for different input frequencies with the target sawtooth signal given for comparison (black dashed line).

The average magnetization of each lattice column was used for the readout yi(t) (Fig. 6a). A discussion about the experimental feasibility of such a grid of local readouts is given in Supplementary Note 7. The readout was used to perform a linear regression \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$o\left(t\right)={\sum }_{i=1}^{n}{w}_{i}{y}_{i}\left(t\right)$$\end{document} and transformation^71,72^, where i is the nanomagnet array column number, wi is the associated weight, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$o\left(t\right)$$\end{document} is the predicted signal and n = 11 is the array dimension. The inverse mean squared error (1/MSE) between \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$o\left(t\right)$$\end{document} and the desired result (the sawtooth) was used as a performance metric.

The 1/MSE for a range of frequencies and amplitudes of the input magnetic fields B1 is shown in Fig. 6b. The black contour indicates a threshold of 1/MSE ~ 30, which can be achieved using only a linear transformation without a reservoir. The best results (ii) are achieved at the input frequency of ∼100 Hz. For other frequencies (e.g. at (i) 0.05 Hz and (iii) 100 kHz), the performance is noticeably worse and cannot be improved to the same level by altering the input field magnitude B1 (Fig. 6b). Interestingly, the best performance is achieved at frequencies around the transition between the first (mavg > 0) and the second (mavg ~ 0, qNN→1) ordering stages of the lattice (Panel 8 of Fig. 1b). A low-frequency input corresponds to the longer relaxation time (green-shaded region of Panel 8 in Fig. 1b) and allows the system to reach the same state with mavg ~ 0 regardless of the input data. This results in a square rather than a sawtooth signal after the transformation (blue trace in panel (i) of Fig. 6c). For a high-frequency input (panel (iii) in Fig. 6c), the lattice cannot respond fast enough to the changes in the signal and cannot reproduce the rise or the sudden drop of the sawtooth. At an optimal input frequency (ii), the system achieves an MSE comparable to other magnetic reservoirs^69^. The performance can be further improved by increasing the size of the lattice or introducing a gradient in the effective anisotropy (see Supplementary Note 8).

In summary, we have demonstrated that the performance of the 2D Ising reservoir depends on how well the ordering dynamics matches the frequency of the input. This is characteristic for this type of computation and, therefore, the ease with which one can engineer the relaxation timescale in the 2D Ising lattices is highly advantageous. This avoids the otherwise necessary and computationally expensive preprocessing of the input data^65,68,70,73^.

Discussion

We have shown with Monte Carlo simulations that the timescale of the relaxation dynamics in 2D Ising lattices can be tuned by varying the effective anisotropy and magnetostatic coupling, and that this can be exploited to achieve an enhanced performance of these lattices for reservoir computing. This information allowed us to fabricate experimental arrays of magnetostatically-coupled perpendicular nanomagnets with the desired properties by engineering the energy landscape through careful selection of materials and geometries.

We have shown that we can achieve a high degree of spontaneous ordering in arrays of out-of-plane nanomagnets. Nevertheless, achieving such a high level of ordering on an even shorter timescale as well as reaching qNN = 1 would be beneficial. Placing our experimental results on the ordering phase diagrams in Fig. 2 reveals a way to do so. Indeed, Region IV with qNN ~ 1 is located below the experimental points (circle, triangle and cross in Fig. 2) indicating sufficient EMC and excessive EEA. Therefore, in order to achieve a perfect ordering at a given timescale, we need a method to tune EEA with more precision. While we have shown that EEA can be reduced by increasing the magnetic layer thickness tCo, going beyond the nominal sub-0.1 nm precision of tCo used in this work is challenging. Therefore, to control the perpendicular magnetic anisotropy in Co/Pt multilayers more precisely in the future, one can use thermal annealing^74–77^. Furthermore, laser annealing^78,79^ would provide a means to locally modify EEA and create spatially complex designs.

In addition to making permanent modifications to the magnetic properties, EEA can be temporarily changed by applying current. The generated spin-orbit torque does not lift the degeneracy of “up” and “down” magnetic states in the out-of-plane easy axis geometry^80,81^. Therefore, on application of a current, the Joule heating and spin-orbit torque (requiring the lattice to be placed on an appropriate layer) lower the energy barrier to switching without interfering with the EMC-driven ordering. Experimental confirmation of this is detailed in Supplementary Note 9. Using voltage and current to reduce EEA would provide a way to dynamically change the relaxation timescale of the lattice, opening the way to create systems with an adaptive temporal response, which is a cornerstone of information processing in the brain^82,83^ and brain-like computing^84,85^. In addition, by engineering the current density distribution (as discussed in Supplementary Note 7) or by applying a voltage to individual nanomagnets^86^, a more local electrical control could be achieved.

We have shown that, by careful choice of the layer number and thicknesses in out-of-plane artificial spin lattice, a magnetostatic coupling energy EMC of tens of kBT can be obtained. The ability to create systems with such high coupling provides an alternative, more energy efficient method to link the nanomagnets compared with the external electronic connections often used in probabilistic computing schemes^49,87^ or for time multiplexing in reservoir computing^65,68^. For such applications, a further increase in EMC could be beneficial and can be accomplished by adding a soft magnetic layer with in-plane anisotropy underneath the out-of-plane lattice^28^. We tested this approach on square arrays of nanomagnets with DNM = 140 nm and observed a substantial increase in the spontaneous ordering (Supplementary Note 10). Another approach to increase EMC is to optimise the shape of the nanomagnets (Supplementary Note 11). Such control of the coupling by altering the shape is uniquely suited to out-of-plane nanomagnets since, for in-plane magnets, a change of shape tends to give a change in the distribution of the demagnetization field.

This freedom to control both EEA and EMC locally will provide a means to create lattices with novel emergent properties. In terms of fundamental science, this means that the thermodynamics can now be studied with different lattice geometries to reveal new phases and phase transitions. In addition to creating thermally-active out-of-plane systems, it will now be possible to pattern mixed lattices incorporating both thermally-active in-plane and out-of-plane nanomagnets. These do not necessarily have an equivalent in bulk crystal systems and are likely to display unusual collective phenomena.

In terms of applications, these time-dependent artificial spin lattices of perpendicular nanomagnets offer an exciting platform for reservoir or probabilistic computing in square or other more complex geometries. The fact that these systems can be electrically interfaced, adapted to the input frequency and have short-term memory based on thermal activity, means that they have all the properties required for next generation computing.

Our work therefore opens vistas across the fields of spintronics, reservoir computing and artificial spin ice, providing precise control of the cooperative behaviour and a flexibility in the design that can be finely tuned for different computing applications.

Methods

Sample fabrication

The films were deposited at room temperature onto high-resistivity Si wafers with a natural oxidation layer using DC magnetron sputtering. A base pressure of ≤1 × 10^−8^ Torr and Ar gas pressure of 3 mTorr were used for the sputtering. The deposited films were processed into devices with electron beam lithography and Ar ion milling. 50 nm-thick 2% hydrogen silsesquioxane was used as an electron beam resist to achieve high-resolution patterning of the nanomagnets. A milling current of 60 mA and angle of 5° were used. Poly(methyl methacrylate) 4% 950k was used as an electron beam resist for patterning the Hall bars and electrodes.

MFM measurements

All MFM measurements were performed at room temperature. Low-moment Bruker MESP-V2-LM probes were used. The samples were initialized by applying three 5-second-long pulses of a 0.9 T out-of-plane magnetic field, sufficient to saturate them and give a uniformly magnetized state on removing the magnetic field before the MFM measurements were carried out. Care was taken to make sure that there was no lateral magnetic field component or remanent magnetization associated with the magnetic field source that would affect the state of the nanomagnet arrays. All MFM scans were performed at zero field.

Energy landscape calculations

The energy landscapes were calculated in Mathematica using a custom code. The profiles were calculated for the polar angle of magnetic moment of the nanomagnet going from 0° to 180° in 1° steps. To obtain the saturation magnetization of 1063 kA/m and interfacial anisotropy of 1.46 mJ m^−2^ used for the energy landscape calculations, the thin films were measured using superconducting quantum interference device vibrating sample magnetometry (SQUID VSM). Details of the energy landscape calculations are given in Supplementary Note 4.

“Hotspice” Monte Carlo simulations

The relaxation of the 2D Ising lattices and reservoir capability were simulated in Python using a custom “Hotspice” Monte Carlo code and employing a magnetic moment of 2.37 × 10^−16 ^A m^2^ for each spin and a temperature of 300 K. 600 input periods were used for the linear regression. 400 input periods were used to test the performance of the signal transformation. Details of the simulations are given in Supplementary Note 1.

Supplementary information

Peer review file Supplementary Material