Magnetization dynamics of weakly interacting sub-100 nm square artificial spin ices

Artificial Spin Ice (ASI), consisting of a two dimensional array of nanoscale magnetic elements, provides a fascinating opportunity to observe the physics of out-of-equilibrium systems. Initial studies concentrated on the static, frozen state, whilst more recent studies have accessed the out-of-equilibrium dynamic, fluctuating state. This opens up exciting possibilities such as the observation of systems exploring their energy landscape through monopole quasiparticle creation, potentially leading to ASI magnetricity, and to directly observe unconventional phase transitions. In this work we have measured and analysed the magnetic relaxation of thermally active ASI systems by means of SQUID magnetometry. We have investigated the effect of the interaction strength on the magnetization dynamics at different temperatures in the range where the nanomagnets are thermally active. We have observed that they follow an Arrhenius-type Néel-Brown behaviour. An unexpected negative correlation of the average blocking temperature with the interaction strength is also observed, which is supported by Monte Carlo simulations. The magnetization relaxation measurements show faster relaxation for more strongly coupled nanoelements with similar dimensions. The analysis of the stretching exponents obtained from the measurements suggest 1-D chain-like magnetization dynamics. This indicates that the nature of the interactions between nanoelements lowers the dimensionality of the ASI from 2-D to 1-D. Finally, we present a way to quantify the effective interaction energy of a square ASI system, and compare it to the interaction energy computed with micromagnetic simulations.

bulk spin-ices is the possibility to directly access experimentally, in real space, the microstates through a variety of techniques, such as magnetic force microscopy (MFM), photoemission electron microscopy (PEEM), and resonant transmission X-ray microscopy (TXM). Possible applications of ASIs range from their use in devices as return-point memories 9 , magnetic cellular automata devices 10 , or magnetic metamaterials 11 , due to the possibility of creating and displacing magnetic monopoles (of Nambu type) 12 in the ASIs.
Until recently, studies on ASIs were performed on athermal systems, as the thermal energy needed to flip the magnetic macrospin of the nanomagnets forming the arrays was out of the experimentally accessible range. The studies on these athermal systems reported on effective thermodynamics, frozen excitations, and field demagnetization protocols in an attempt to access the ground-state ordering [13][14][15][16][17][18][19] . Recent reports on thermal ASIs have opened the door to the study of magnetization dynamics in these systems. These systems include thermal annealing processes taking place during fabrication 20 , and systems where the anisotropy barrier of the nanomagnets has been tuned to be in a thermally accessible regime by judicious choice of a magnetic material with a lowered Curie temperature (T C ) 21,22 and by carefully heating the sample above its blocking temperature (T B ) 23 . These reports were shortly followed by studies of thermally fluctuating ASIs which have been imaged via PEEM [24][25][26][27] and TXM 28,29 in real time in a variety of geometries, and recently via MFM imaging of the intermediate thermally stable states after a temperature quenching process 30 . Nonetheless, it has only been very recently that experiments where a phase transition from the superparamagnetic regime ( > > T T T C B ) to the ASI regime have been demonstrated. These report glassy freezings of the magnetization dynamics of square ASI systems measured by X-ray photon correlation spectroscopy 31 and magnetometry 32 , in both cases following a Vogel-Fulcher-Tammann law [33][34][35] , a phenomenological law used to explain, among other systems, the behaviour of spin-glasses.
Building on the seminal MFM measurements of ASI, the huge advances in our knowledge of ASI have typically required access to intense x-ray synchrotron sources. In this context, we present a study to investigate the magnetization relaxation dynamics of square ASI systems by means of SQUID magnetometry. With this technique it is possible to explore the collective dynamics of the whole array of nanoelements composing the ASI, in comparison to previously mentioned techniques where only small portions of the sample are inspected with each measurement. The ASIs are formed by nanomagnets made of Permalloy (Ni80Fe20) with lateral sizes of 68 nm × 22 nm, with two different thicknesses: 5 nm and 6 nm; and three different lattice spacings for each thickness: 138 nm, 175 nm and 208 nm (Fig. 1), making a total of six samples studied. The symmetric disposition of the nanoelements avoid non-uniform magnetization states in them 36 , ensuring their Ising-like behaviour. We observe that the relaxation dynamics of the studied square artificial spin-ices follow an Arrhenius-type Néel-Brown behaviour, contrary to what is reported previously in similar square artificial spin-ice studies 32 . Zero field cooling and field cooling measurements have been performed for all of the samples, together with magnetization relaxation measurements at fixed temperatures. The analysis of the data extracted from the measurements suggest 1-D Figure 1. Square artificial spin ice and its magnetization dynamics process. Scanning electron microscopy images of square ASIs with three different lattice spacings: (a) 138 nm, (b) 175 nm and (c) 208 nm, made of Permalloy nanomagnets with lateral dimensions of 68 nm × 22 nm. Panel (a) shows the two directions [10] and [11] along which the field was applied whilst measurements were performed. The bottom panels (d-f) show schemes of the magnetization states in the pattern when a saturating field is applied right before starting the measurement (d); an intermediate state, at a certain time after starting the measurement, showing individual flips of nanomagnets that create monopolar charges (positive: yellow, and negative: green dots) connected by Dirac strings (encircled by a dashed ellipse in panel (e)) 24 (e); and the final magnetization state, showing ground-state ordering of the square artificial spin-ice, after a full relaxation of the magnetization (f). magnetization dynamics processes. They also show a negative correlation of the average blocking temperatures of the samples with their interaction strength, indicating that the ASIs are in a weakly interacting regime, which is further corroborated following the theory of Shtrikman and Wohlfarth 37 . The experimental results are supported by Monte Carlo simulations of the magnetization processes in the samples studied. In systems of weakly interacting nanoparticles embedded in a polymeric matrix it has been observed that the blocking temperature of the system increases with the interparticle distances 38,39 ; however, this is the first evidence for such a behaviour in a system where the interacting elements are not spherical nanoparticles, but elongated nanomagnets that are arranged in a pre-defined ordered geometry. In addition to the quantification of the interaction strength between the nanoelements in each ASI, the methodology followed here also gives information about the dimensionality of the system. Thus, the present work provides a systematic study of the effect of frustration in the dimensionality of artificial spin-ice systems with different geometries, and opens the door to the design and analysis of desired exotic states and emergent behaviours 40 .

Results
Zero field cooling/field cooling measurements. In order to extract the characteristic relaxation times of the magnetization dynamics of our square ASIs we need to identify the temperature regions in which the samples are thermally active. This region of interest is readily identified using SQUID magnetometry. Details of the measurements can be found in the methods section. Zero field cooling (ZFC)/field cooling (FC) curves have been measured for all of the samples and are shown in panels (a) and (b) of Fig. 2. Upon cooling down in the absence of any external field, from a temperature above the average T B of the system, the nanomagnets will undergo slowing down of the magnetization dynamics until they reach a certain temperature below which the system will freeze into an ordered low energy ground state, where regions of nanomagnets arranged similarly to panel (f) of Fig. 1 are separated by higher energy vertex chains that form domain walls between ordered ground state regions 21 . The range of temperatures at which each sample will be thermally active is identified in the ZFC/FC measurements. The lower bound is given by the temperature at which the magnetization starts to increase in the ZFC (where we expect slow dynamics and long relaxation times), and the upper bound by the average T B . It is expected that for temperatures slightly below the average T B fast dynamics and short relaxation times will be observed. From the lower branch (ZFC curves) of the plots in Fig. 2, the temperature range where the 5 nm thick samples are thermally active and with dynamics observable on laboratory timescales lies between 190 K and 270 K, as observed in panel (a), whereas for the 6 nm thick square ASIs it lies between 300 K and 380 K, as www.nature.com/scientificreports www.nature.com/scientificreports/ identified from panel (b). Note that above the upper bound temperatures specified here the samples are also thermally active, but their dynamics are so fast that we cannot experimentally probe them.
In order to observe the trend of the average T B for the samples studied we have plotted the difference between the The results obtained for the MC simulations show an inverse correlation between the average T B of the samples and the strength of the dipolar interactions between the nanoelements, which is the same effect observed in the samples measured. Nonetheless, a constant temperature difference between the measured and simulated average T B exists for each set of samples: for the 5 nm thick set of samples the measured average T B are 22 ± 8 K lower than the calculated ones, while for the 6 nm thick set of samples they are 30 ± 10 K higher. These differences can be attributed to discrepancies between the nominal and real thicknesses of our samples, with the nominal 5 nm thick set of samples thinner than 5 nm, and the nominal 6 nm thick set of samples thicker than 6 nm. The differences between blocking temperatures of samples with different lattice spacings have been quantified, finding that there is a 2:1 ratio of T B differences between the 138 nm and 175 nm lattice spacing samples, and the 175 nm and 208 nm ones, both for the 5 nm and 6 nm thick set of samples. This 2:1 ratio of T B between different lattice spacings is also maintained in the simulations.

Magnetization relaxation measurements.
After identifying the interesting region of temperature where each sample shows thermally activated dynamics we can measure the thermal relaxation of the magnetization dynamics in each sample for different fixed temperatures: for the 5 nm thick set of samples between 185 K and 265 K, in steps of 10 K (a total of 9 measurements for each sample); and for the 6 nm thick set of samples between 290 K and 370 K, also in steps of 10 K (again, 9 measurements for each sample). As an example, the recorded measurements of the average magnetization evolution in time, plotted in normalised form, at each fixed temperature for the 175 nm lattice spacing 6 nm thick sample are shown in Fig. 4. Each measurement is normalized to the www.nature.com/scientificreports www.nature.com/scientificreports/ magnetization saturation of each sample measured during the application of the saturating field before starting the magnetization relaxation measurement protocol. The normalized moment m/m S is fitted to a stretched exponential 42 of the form: S r where a characteristic relaxation time, t r , of the magnetization dynamics and the stretching exponent, β, are extracted from each fit. The superimposed lines plotted on top of each measurement (scattered points plots) correspond to the fitted stretched exponential for that measurement. While the meaning of the extracted relaxation times for each sample at each temperature is of essential importance to identify the type of relaxation dynamics followed by the sample, the extracted stretching exponents also possess important information about the relaxation dynamics, as will be discussed. The dependence of the relaxation times and stretching exponents on the temperature, as extracted from the fits of each measurement to Eq. 1 for the six samples studied here, are presented in panels (a) and (b) of Fig. 5, respectively. For the sake of comparison, similar relaxation measurements taken along the [11] direction of the 175 nm lattice spacing 6 nm thick sample (45 degrees from the easy axes of the nanoelements) are also included in Fig. 5.  . Relaxation times and stretching exponents of the samples studied. Temperature dependence of the relaxation times (a) and the stretching exponents (b) extracted from the stretched exponential fits of the measured time evolution of the magnetic moment on the 5 and 6 nm thick samples, for the three different lattice spacings. The measurements were performed along the [10] direction, except for the blue dataset, which corresponds to a measurement of the 6 nm thick 175 nm lattice spacing sample along the [11] direction. (2019) 9:19967 | https://doi.org/10.1038/s41598-019-56219-y www.nature.com/scientificreports www.nature.com/scientificreports/ The temperature dependence of the relaxation times obtained for each of the samples are fitted to an Arrhenius-type Néel-Brown law: where the activation temperature T A , a measure of the energy barrier, depends on the anisotropy energy (E A ) of the island and the interaction energy (E int ) between elements. From this, the activation temperature is given by: is the anisotropy constant, and ΔD is the difference between the in-plane demagnetizing factors of the nanoelements 43 .
For each sample, the fitting of the temperature dependence of the relaxation times to Eq. 2 gives us two parameters: the value for  Table 1 and the ones for the three 6 nm thick samples are presented in Table 2. From the fitted values, and calculating the shape anisotropy energy as previously mentioned, the value of the measured interaction energy is calculated for each of the samples. This interaction energy is compared to the extracted dipolar interaction energies computed via micromagnetic simulations 46 , by subtracting the magnetostatic energies computed for an unfavourable alignment of the neighbouring + NN NN ( 2 ) macrospins from those obtained with a favourable alignment. The discrepancies between the computed and experimentally obtained interaction energies are due to discrepancies between the real volume of the nanomagnets in the samples and the nominal volume used in the micromagnetic calculations.
Even if the relaxation times of the three 5 nm thick samples have been fitted to the Néel-Brown law, a more detailed inspection of Fig. 5 suggests a different behaviour for the 208 nm lattice spacing sample than that of the 138 nm and 175 nm lattice spacing samples. The relaxation time is always shorter for the 138 nm lattice spacing sample than for the 175 nm one for each temperature measured; a trend that is observed, by looking at Fig. 5, not only in the samples with those lattice spacings of the 6 nm thick set of samples, but also in the 208 nm lattice spacing sample. Nevertheless, this trend is not followed by the 5 nm thick 208 nm lattice spacing sample with respect to the other two 5 nm thick samples. Due to the smaller magnetic signal produced by this sample, as it is the one with the least total magnetic material, the magnetic relaxation measurements are more noisy than with any of the other samples measured here. This led to higher correlated fitting parameters, causing the error bars to be relatively small, but the results are less trustworthy due to the higher noise observed in the measurements of this sample. The extracted value for the relaxation and exponent behaviour at = T 205 K is noteworthy: the sample is essentially static at this temperature, and consequently the poor fitting of the stretched exponential to that measurement yield non-realistic values for t r and β. This is reflected in the bigger error bars in t r and β for the 205 K measurement of the 5 nm thick sample shown in Fig. 5.

Discussion
Zero field cooling/field cooling. When studying superparamagnetic relaxation it is generally observed that T B increases as the average interaction strength between particles is increased 47,48 . In the ASIs studied here, where the nanomagnets have an Arrhenius-type Néel-Brown behaviour of the magnetization dynamics, the energy barrier between the two stable magnetization states of the nanoelements (which defines the individual T B of the nanoelement) is given by the sum of the shape anisotropy of the nanoelement (independent of the lattice spacing and common to the three samples on each set) and the interaction energy due to the dipolar magnetic interactions  www.nature.com/scientificreports www.nature.com/scientificreports/ between neighbouring nanoelements. From the ZFC-FC measurements performed in the samples studied here, and supported by Monte Carlo simulations, we can clearly observe that their average T B decreases with the increase of the interaction strength between the nanoelements forming each sample, both for the 5 nm and 6 nm thick set of samples. This suggests that the samples are in a weakly interacting regime 38 , as it is discused hereafter. The Monte Carlo simulations of the ZFC/FC protocols show that the average T B of the ASIs studied are very sensitive to small changes in the thickness of the nanoelements, as derived from the results exposed in panel (c) of Fig. 3.
In a study on systems of interacting ferromagnetic nanoparticles by Mørup and Tronc 38 , they observe an inverse correlation between the average T B and the interaction strength, and they developed a model to explain this effect. The key ingredient of this model is that particles with uniaxial magnetic anisotropy are exposed to time-dependent dipolar fields coming from the neighbouring particles. At a certain time each particle is exposed to a dipolar field, while its magnetization fluctuates between the two stable states. Those fluctuations have frequencies of the order of 10 10 -10 12 s −1 , while occasionally the magnetization vector will access the energy barrier for some value of the angle φ defined by its magnetic moment and the average dipolar field sensed by the particle.
For some values of that angle, the energy barrier is lowered, leading to a decrease of the relaxation time due to the dipolar interactions. They derived an expression for the average T B , which yields an inverse correlation of the average T B with the interaction strength. In our system, we have uniaxial particles (nanomagnets) that are subjected to the dipolar fields coming from the neighbouring nanomagnets. Furthermore, due to the nature of the lithography process in our system, there are very small random deviations of the alignments of the easy axes of the nanomagnets, leading to a similar effect to that accounted for in the model with φ and the dipolar fields.

Magnetization relaxation.
According to previous studies of interacting superparamagnets 49 , higher (lower) interaction energies are expected in samples with higher (lower) average T B . Based on this, and with the behaviour of the average T B observed in our samples, it is expected that the effective interaction energies extracted from the magnetic relaxation measurements for the strongest interacting samples (smaller lattice spacings) will be smaller than the ones for the weaker interacting samples (bigger lattice spacings).
From the data presented in Table 1 it is observed that the extracted interaction energies for the three 5 nm thick samples overlap each other, as we are in a such a weakly interacting regime it is impossible to measure a distinct effect on this set of samples, based only on the interaction energies extracted from the relaxation measurements. This means that, although with static magnetometry measurements (ZFC/FC) we can observe differences in the behaviour of our samples, with magnetization dynamics measurements we cannot infer a different behaviour in these samples. Nonetheless, from the experimentally obtained values for the interaction energies of the 6 nm thick set of samples, presented in Table 2, it is clear that, while the 175 nm and 208 nm lattice spacing samples have similar overlapping interaction energies, the interaction energy of the 138 nm lattice spacing sample is smaller than the other two. This means that the difference in interaction strength between the 138 nm and 175 nm lattice spacing samples marks a threshold, below which the extracted interaction energies are indistinguishable between samples from the magnetization dynamics measurements. Furthermore, the fact that the interaction strength of the 138 nm lattice spacing sample is smaller than those of the 175 nm and 208 nm lattice spacing samples is in good agreement with the inverse correlation of the T B and their interaction strengths, obtained from the ZFC/FC measurements.
Following the theory of Shtrikman and Wohlfarth 37 , we can establish a criterion to determine whether our samples are in a weakly interacting regime or not. This criterion is based on the comparison between the interaction energy of the system and the anisotropy energy of the single element, and states that if  E E int A the samples are in a strongly interacting regime, whereas if  E E int A the samples are in a weakly interacting regime. In our case, the interaction energies are around two orders of magnitude smaller than the anisotropy energy of the single islands composing the ASI arrays, thus meaning that we are in a weakly interacting regime in all the cases studied here. In previous studies of magnetization dynamics in square artificial spin ice systems 31,32 , the lattice spacings and dimensions of the nanomagnets composing the samples studied there indicate that these samples should be in a weakly interacting regime according to the Shtrikman and Wohlfarth theory. However, they did not observe an inverse correlation of the average T B with the interaction strength. However, in those ASIs the activation volumes were much smaller than the nanoelement volumes, which is an effect previously reported in ASIs with larger nanoelements 30 . This means that the nanoelements are able to reverse at much lower temperatures than those expected for the full nanoelement volume barrier. The nanoelements of the study presented here have an activation volume much closer to the full nanoisland volume, so that the Shtrikman-Wohlfarth theory of being in a truly weakly interacting limit applies, hence we expect to observe an inverse correlation between T B and the interaction strength.
Comparing the experimentally obtained interaction energies with the calculated values of the dipolar energies for each sample, it can be seen that although they do not overlap, they are of the same order of magnitude and their discrepancies are attributed to a reduction in the volume of the real sample with respect to the values used for the micromagnetic computations.
It is worth noting that the calculations of the dipolar energies presented in Tables 1 and 2 have been obtained relaxing the magnetization state from a situation in which all the nanomagnets were in a perfect single-domain state, resulting in a canting of the magnetic moments along the edges of the nanoelements, but not involving any dynamic process. Therefore, we do not expect to observe an inverse correlation between the dipolar interaction energies extracted from the calculations and the interaction strengths, as this is a result of the magnetization dynamics processes in the samples, as discussed previously.
From an inspection of panel (a) of Fig. 5 it is evident that the measurements with the samples mounted along the [10] and [11] directions are indistinguishable, as the relaxation times match perfectly, as expected, due to the fact that the underlying magnetization dynamics processes are identical in the two sublattices forming the square ASI arrays.

Scientific RepoRtS |
(2019) 9:19967 | https://doi.org/10.1038/s41598-019-56219-y www.nature.com/scientificreports www.nature.com/scientificreports/ The meaning of the stretching exponents obtained from the magnetization relaxation measurements is related to the dimensionality of the dynamic processes taking place. As a result of the lithography process there is a distribution of the energy barriers between the two stable magnetization states in each of the nanoelements, giving rise to a random distribution of energy potentials in the square array. This maps on to the so-called trapping model 50 , allowing us to extract information about the dynamic processes from it. This model predicts that the stretching exponent obtained from the fits of the dynamics measurements to Eq. 1 takes the form: where d is the dimensionality of the system. The stretching exponents for the 5 nm thick set of samples (panel (b) of Fig. 5) scatter around values of β = 1/3, excluding the irregular behaviour of the 5 nm thick 208 nm lattice spacing sample. This value of β = 1/3 suggests 1-D dynamic processes, similar to those observed in square artificial spin-ice systems studied by PEEM 24,26 . These 1-D processes consist of the formation and propagation of chains of nearest neighbour nanomagnets undergoing reversal processes, being in the so-called string regime. The string and domain regime formation in thermally active ASIs was first shown by Farhan et al. 24 for strongly interacting arrays. More recent work comprised the study of both strongly interacting square ASIs, showing similar results, and weakly interacting ones where vertical string propagations occur 30 . Both cases present 1-D thermal excitations which happen in the initial part of the relaxation process or, compared to the present work, at lower temperatures. Based on this, vertical 1-D excitations are expected for the 5 nm thick samples, whereas 1-D diagonal strings are expected for the 6 nm thick samples, leading in both cases to a convergence of β = 1/3. The more complex temperature dependence for the 6 nm thick set of samples is attributed to a mixed fluctuating state of nanoelements of both string formation and propagation processes, together with domain formations. This leads to a departure from the one dimensional behaviour, increasing the value of β towards 1/2, which corresponds to a 2-D fluctuating system. The temperature dependence of the stretching exponents shows an average decrease in β for higher temperatures in all the samples. This is attributed to the confinement of the fluctuations that become domain walls, as shown in a previous study by Budrikis et al. 51 , leading to a recovery of the 1-D nature of the thermal excitations at higher temperatures and, therefore, to a convergence towards β = 1/3 of the stretched exponents. It is worth noting that the measurements at both edges of the temperature region studied have relaxation times that are in the limits of the detection of the technique, and the corresponding fits to the measurements have a higher χ 2 value.

conclusions
To summarize, we have studied the magnetization dynamics of sub-100 nm square artificial spin ice samples with different thicknesses and lattice spacings by means of SQUID magnetometry. From the measurements we can conclude that the magnetization relaxation times obtained as a function of the temperature follow a simple Arrhenius-type Néel-Brown behaviour. This is expected from interacting superparamagnetic nanoparticles 52 that do not freeze into a glassy state 32 .
The average blocking temperatures have an inverse correlation with the interaction strength of samples formed with nanomagnets with similar dimensions. These results are supported by Monte Carlo simulations. This unexpected negative correlation of the interaction strength with the average T B is a first-time observation in systems of elongated nanomagnets disposed in a non-random ordered geometry, although it has been observed previously 38,39 for interacting spherical randomly distributed ferromagnetic nanoparticles. In a study carried out by Mørup and Tronc 38 this effect is explained with a model that assumes uniaxial magnetic anisotropy in the weakly interacting nanoparticles that are exposed to dipolar fields from the neighbouring nanoparticles, both characteristics being found in the square ASIs studied here.
The magnetization relaxation measurements for each temperature are fitted to a stretched exponential function, from which we obtain a characteristic relaxation time and a stretching exponent. For each set of samples with the same thickness, the relaxation times have a positive correlation with the lattice spacing, the only exception being the least interacting sample (208 nm lattice spacing, 5 nm thick), whose nanoelements are in the limit of being non-interacting to super weakly interacting. The stretching exponent values give information about the dimensionality of the magnetization dynamics processes, and a value of β = 1/3 implies one-dimensional magnetization dynamics processes. Departures from β = 1/3 are attributed to mixed regimes of string formation and propagation processes with domain formation. The reduced dimensionality of the system, which shows 1-D magnetization processes in a 2-D ensemble of nanomagnets, is a direct consequence of the effects of the geometric frustration present in the square ASIs.
The experimental method proposed here is not only able to quantify the interaction energy of ASIs, but can also be applied in general to any ensemble of interacting nanomagnets. This is not restricted to ensembles following an Arrhenius-type Néel-Brown behaviour, but also to those following any other law (e.g. Vogel-Fulcher-Tammann) where an analytic expression for the temperature dependence of the relaxation times exist.

Methods
Growth and structural characterization. The ASI samples studied here have been fabricated by means of electron-beam lithography, following standard procedures. Firstly, a layer of ZEP resist is spin-coated on a Si [100] substrate after cleaning the substrate. Then, standard exposure to the electron beam procedures are followed, to lithographically define the nanoelements with the desired lateral dimensions and lattice spacings, followed by a developing process of the resist after exposure by rinsing in a chemical developer. A thin layer of Permalloy www.nature.com/scientificreports www.nature.com/scientificreports/ (Ni80Fe20) is deposited onto the masked substrate, followed by a 2 nm thick Al cap (to prevent the samples from oxidation) and finally a lift-off process results in the square ASI patterns (Fig. 1). Surfaces of 2 mm × 2 mm were lithographically defined and covered by the patterns, to optimize the signal in the magnetometer.

Magnetic characterization. The magnetic characterization was performed using a commercially available
Quantum Design SQUID magnetometer. Zero field cooling curves have been measured by heating the samples to 400 K (above TB) and bringing them to 10 K in a field less than 10 Oe (remanent field when not applying any field by the magnetometer), to ensure that the samples are in the lowest energy state (ground state depicted in panel (f) of Fig. 1), followed by a measurement of the moment on the sample, as a function of the temperature, from 10 K to 400 K, in a probe field of 30 Oe. Field cooling curves have been measured by heating the samples to 400 K, and then measuring their magnetic moment from 400 K to 10 K under a probe field of 30 Oe. The zero field cooling and field cooling curves presented in Fig. 2 have been measured with the fields applied and the samples mounted along the [10] direction, and with a heating/cooling rate of 2 K/min. The data shown in panels (c) and (d) of Fig. 2 have been smoothed following an adjacent averaging smoothing method, where the four neighbours of each data point are averaged, replacing the value of the data point by the new average value. Thermal relaxation measurements have been performed following this procedure: firstly, a saturating field of 5000 Oe was applied, forcing the magnetic moments to align with the applied field, resulting in a magnetic configuration of the array similar to the one depicted in panel (d) of Fig. 1. Then, the magnetic field is removed and the time evolution of the magnetization (in the absence of any external field) is measured for 600 s, resulting in graphs similar to the individual scattered plots shown in Fig. 4. Note that the magnetization plots shown in that figure are normalised to the magnetization saturation measured during the application of the saturating field, in order to perform the fitting of the stretched exponential function to extract the characteristic relaxation time. All the measurements performed in this study have been done with 1 s averaging time. The square patterns have been measured mounting the samples in such a way that the field was applied along the [10] (parallel to one of the sublattices of the square array) and the [11] (45 degrees from the [10]) directions, showing similar results (see blue and red datasets in Fig. 5).

Micromagnetic simulations.
The micromagnetic simulations of the magnetostatic energies have been performed by means of OOMMF 46 , assuming the nominal island sizes, with cell sizes of 2 × 2 × 1 nm 3 , well below the exchange length of Permalloy in every dimension, and considering nearest-neigbour and second nearest-neighbour interactions between elements. The material properties used are the ones defined for Permalloy in the OOMMF package by default. The use of cell sizes 2 × 2 × 1 nm 3 is motivated by the fact that similar results are obtained with 1 × 1 × 1 nm 3 cell sizes, and the former reduces drastically the computing time.
Monte carlo simulations. The energetics and magnetisation processes of the system described in Fig. 1 were also theoretically investigated using a standard Metropolis Monte Carlo algorithm, assuming that the nanoelements forming the system are identical. Here, the spin Hamiltonian H has the form 53 : denoting terms for the dipolar interaction and externally applied field, respectively. We consider the magnetic nanoislands to be well-separated enough so that they can be considered as Ising-like spins and can be taken as point dipoles. In this case, the interaction between the magnetic moments is given by the expression 54 where l is the lattice spacing and V is the volume of the naoelements 55 .
The effect of an external field, such as the one applied to the real system, can be calculated by evaluating s i i app The method followed in order to perform the Monte Carlo simulations is similar to those employed in literature in the context of nanoparticles and fine magnetic nanostructures [56][57][58][59] and, more in particular, in Artificial Spin-Ices 51,60,61 . The method is described as follows: for the field cooling simulation, the initial state at a temperature T well above TC of the ensemble consists of all the magnetic moments saturated and pointing towards the easy axis of one of the sublattice (the [10] direction described in Fig. 1). The standard Metropolis Monte Carlo method then calculates the equilibrium states for each temperature, where 10 5 steps are sufficient for convergence, and once equilibrium is reached the temperature is decreased to the next simulation temperature. The temperature step used for the simulations is 7 K. A field of 30 Oe is applied during the simulation of the cooling process along the [10] direction. For the zero field cooling protocol, the initial state is the one depicted in panel (f) of Fig. 1 at a temperature of 10 K, and the temperature is increased in steps of 7 K up to 400 K under a probe field of 30 Oe. The components of the magnetization that contribute to M are those along the [10] direction. The stopping criterion for the simulation is reached when there are no significant spin fluctuations within each Monte Carlo step, corresponding to spin fluctuations of less than 2% in the averaged magnetization. The sensitivity of the Monte