Collective magnetic dynamics in artificial spin ice probed by AC susceptibility

Merlin Pohlit, ∗ Giuseppe Muscas, † Ioan-Augustin Chioar, Henry Stopfel, ‡ Agne Ciuciulkaite, Erik Östman, Spyridon D. Pappas, § Aaron Stein, Björgvin Hjörvarsson, Petra E. Jönsson, and Vassilios Kapaklis ¶ Department of Physics and Astronomy, Uppsala University, Box 516, SE-75120 Uppsala, Sweden Center of Functional Nanomaterials, Brookhaven National Laboratory, P.O. Box 5000, Upton, New York 11973, USA (Dated: January 15, 2020)


I. INTRODUCTION
Artificial Spin Ice (ASI), i.e. arrays of magnetostatically coupled ferromagnetic islands -mesospins 1 -fabricated by nanolithography 2,3 , exhibit collective phenomena, and importantly, their interaction strength and geometry can be tailored almost at will [4][5][6][7][8] . Properly designed to support thermal fluctuations, ASI systems can serve as a platform for the investigations of thermal magnetization dynamics and freezing transitions in tailored nanostructures, which can also me used to mimic the dynamical properties of frustrated, naturally-occurring magnetic spin systems 4,5,[9][10][11][12] . Insights into the freezing transition and the nature of the frozen low temperature states were obtained by investigations using magnetometry 13 and synchrotron-based scattering-and microscopy-techniques [14][15][16] . With the exception of early work based on temperature dependent magneto-optical measurements 9 and more recent works using synchrotron-based magnetic microscopy 4 and muon relaxation 17,18 , experimental studies of thermally induced transitions are scarce. Furthermore, experimental techniques based on synchrotron radiation and muons impose limitations on availability and accessible time-scales. To this end, AC susceptibility is a well established and accessible technique for probing magnetization dynamics, giving access to a wide frequency range 19,20 .
In this work, we report on AC susceptibility measurements of thermally active extended square ASI arrays measured in a wide frequency and temperature range. We study arrays that are composed of close to identical mesospins, with different gaps between the elements, in order to explore the influence of coupling strength on the collective dynamics. Exploring the frequency dependence of the AC susceptibility signal we employ the VFT law, that can be used for describing the low-field magnetic relaxation of weakly interacting nanoparticle systems 21,22 but has recently been applied also to ASI systems 13,14 , attempting to extract parameters that can be directly related to the magnetostatic energies of the ASI arrays.
We discuss the validity of this simplified approach and address the limitations of such models, in the framework of thermal ASI arrays.

II. EXPERIMENTAL DETAILS
The extended square ASI structures were produced by post-patterning of δ-doped Pd(Fe) thin films 23 , employing electron-beam lithography (EBL) followed by argonion milling. The films, consisting of 40 nm Palladium, 2.2 monolayers of Iron, and a 2 nm Palladium capping layer, were all grown on a Vanadium seeding layer on top of Magnesium oxide (MgO(001)) substrates by DC magnetron sputtering. The total effective thickness of the magnetic layer (Fe and magnetically polarized Pd) was previously estimated to be 2 nm 23,24 . Vibrating Sample Magnetometetry (VSM) revealed that the temperature dependence of the in-field volume magnetization can be described by: M s (T ) = M 0 (1 − T /T 0 ) 0.5 , with T 0 = 410 K 13 . The size of the islands, the lattice parameter and the distance between the islands were determined after the EBL process using Scanning Electron Microscopy (SEM). Typical SEM image is shown as inset of Fig. 1. All the islands have a length of 310 ± 15 nm and a width of 140 ± 15 nm. Arrays with different interisland lattice spacing, d, of 380, 420 and 460 nm were prepared which yields different gaps between the islands (g = 70, 110, and 150 nm). The difference in distances between the islands results in different magnetostatic coupling strengths of the mesopins.
The magnetic moment of a mesospin at T = 5 K was determined to be m 0 = M 0 V = 6.5 × 10 6 µ B , where V is the volume of the magnetic material in the island. 25 Micromagnetic calculations using the MuMax 326 revealed non-collinearities of the moment within the islands. A reduced effective moment of approximately m 0,eff = 0.65 × m 0 = 4.2 × 10 6 µ B at T = 0 K, is used to compensate for the dynamic non-collinear internal magnetic structure of the elements (see Andersson et al. 13 , Bessarab et al. 27 and Bessarab et al. 28 ). Furthermore, the intrinsic ef- The dynamic response of the three extended square ASI arrays was investigated by AC susceptibility employing a magneto-optical Kerr effect (MOKE) magnetometer in longitudinal mode. For this purpose the samples were mounted in a cryostat with optical access (4 K ≤ T ≤ 300 K). All experiments were performed with a 20 mW laser with a wavelength of 660 nm. A pair of Helmholtz coils was used to generate a small sinusoidal magnetic field with a given frequency, aligned along the [10]-direction of the nanostructured array (see inset of Fig. 1, upper panel). For frequencies between f = 0.1 − 333 Hz an amplitude of µ o H ac = 0.1 mT was used, while for the higher frequencies f = 1111 Hz and f = 3333 Hz the amplitude was reduced to µ o H ac = 0.01 mT. A lock-in amplifier (Stanford Research SR830) was used as a voltage source for generating the magnetic field and measuring the AC susceptibility at the corresponding frequency. The sample was shielded from the earth's magnetic field by a double-wall mu-metal cylinder, and demagnetization was performed prior to each cool-down.
Applying the magnetic field along the [10]-direction of the samples (see Fig. 1 ) in the longitudinal MOKE configuration probes primarily the mesospins with long axes parallel to the direction of the magnetic field. Both inand out-of-phase components, χ and χ , were recorded during warming up to the maximum temperature of T = 300 K, starting either at T ≈ 80 K (f = 33 Hz) or T = 200 K (other frequencies) and by stabilizing the cryostat in discrete 1 K temperature steps while keeping the amplitude and the frequency of the magnetic field fixed. In every step, sufficient time (60 s) was taken to allow the measurement system (cryostat) to stabilise before starting a measurement 29 . The time for data acquisition at each temperature step was increased for lower frequencies to ensure an acceptable signal to noise ratio. In this respect it should be noted that the magnetic signal of the samples is extremely small, as it originates from a 2.2 ML magnetic Fe layer in Pd, with a coverage of about ≈ 37 − 54% (depending on the geometry of the pattern). For the weakest coupled sample (d = 460 nm) this equals an effective magnetic coverage of a sub-monolayer thick continuous Fe film in Pd.

III. RESULTS AND DISCUSSION
For each of the three studied arrays, a single peak in both χ and χ is observed, as illustrated in Fig. 1. The shape of the peaks is similar for the three arrays, while a shift in the peak position T m towards higher temperatures is observed with decreasing inter-island distance, i.e. increasing inter-island coupling strength. We attribute the maximum in the AC susceptibility to the condition when the average relaxation time is equal to the observation time window τ m (T m ) = 1/(2πf ) 30 , where f is the frequency of the applied magnetic field and T m is the corresponding temperature commonly referred to as the blocking temperature in superparamagnetic samples. Consequently, by determining the peak positions T m for different observation time windows (frequencies), the average relaxation time of the system as a function of temperature can be extracted from the AC susceptibility data.
In order to investigate both the temperature dependence of the relaxation time and the effect of inter-island interactions, T m was determined for different frequencies, effectively resulting in different observation time windows. Typical results are illustrated for the array with the strongest interactions (d = 380 nm) in Fig. 2. T m was determined by fitting a parabola to a region of interest around the maximum of each curve. For each of the three arrays, nine T m values are extracted, corresponding to the employed excitation frequencies f , as illustrated in Fig. 2. The peaks are found to shift towards higher temperatures with increasing frequency.
A. Fitting the experimental data using the Vogel-Fulcher-Tammann law We start the analysis by employing the empirical VFT law for the relaxation time τ of weakly interacting magnetic particles 31 , an approach which has previously been used to describe the relaxation in artificial spin ice structures 13,14 . Within this approach the relaxation time can be calculated by where τ 0 , E K , k B , T and T F correspond the inverse attempt frequency, the intrinsic energy barrier, the Boltzmann constant, the temperature and the freezing or Fulcher temperature, respectively. While the mesospin's energy barrier E K depends on the shape and magnetisation of the mesospins, the Fulcher temperature T F is indicative of the interaction strength of the elements. The intrinsic energy barrier E K in Eq. (1) is attributed to the shape anisotropy and its temperature dependence can be captured by: where µ 0 is the vacuum magnetic permeability, ∆N is the island differential demagnetizing factor according to Osborn 32 and V the moment of the island. Taking into account the temperature dependence of the island saturation magnetization M s in the VFT law, the temperature dependence of the energy barrier becomes is the energy barrier at zero temperature.  Table I. Inset: Systematic reduction of the extracted freezing Temperatures TF with increased island spacing i.e. reduced coupling strength. Fig. 3 presents the extracted T m values along with their corresponding VFT fits. The fitting was performed by reversing Eq. (2), i.e. solving for T , and using a Levenberg-Marquardt algorithm to find the best fit for E 0 K and T F . This procedure facilitates weighting the data points with their respective uncertainties on the maxima positions T m , while the uncertainty in the magnetic field frequency given by the lock-in amplifier was considered negligible. 33 The uncertainties of T m were taken from the fits used for determining the maximum positions, while the uncertainty in selecting the individual temperature ranges used for peak finding were considered negligible. This approach qualitatively captures the effect that the uncertainty is larger for the highest and lowest frequencies measured. Since all arrays are composed of the same size elements, the flipping time was assumed to be constant, τ 0 = 10 −11 s, in accordance with a previous relaxation study 13 . Same energy barrier, E 0 K was used for all three gap sizes studied. The summary of the results from these analysis is found in Table I. For comparison, an independent magnetostatic estimation of the intrinsic energy barrier was made, calculating ∆N using the Osborn methodology 32 . The effective magnetic moment of an island m 0,eff was used in order to take the non-collinearities of the moment within the islands. By using this approach and the temperature dependence of m eff , the mean value of the energy barrier was determined to be E K /k B ≈ 3200 K at 250 K. Using the temperature-scaling for the fitted value of E 0 K , we obtain E K (T = 250 K) = 2510 K ± 120 K, a value which, for the current case, is not so far from the magnetostatic estimation.

B. On the extraction of characteristic interaction energies
With a single energy barrier fitted for all three datasets, the impact of the different gap sizes is inherently accounted for by the Fulcher temperature, T F . This dependence is represented in the inset of Fig. 3, highlighting a connection between T F and the inter-island magnetostatic interactions. This further raises the question of whether the Fulcher temperature can be systematically and accurately used to extract the characteristic interaction strength between mesospins.
In the framework of a weak-coupling regime in spin glasses, Shtrikman and Wohlfarth 31 offered a recipe for extracting the typical interaction energy between the magnetic components by using the mean-field based formula: where E 0 K represents the average intrinsic energy barrier of a single magnetic element, while E 0 i is the characteristic interaction energy between the magnetic constituents,  Fig. 3, assuming a fixed τ0 = 1 · 10 −11 s, as done in Andersson et al. 13 for the same samples. Right column: Characteristic magnetostatic interaction energies obtained from the VFT fits in the weak coupling formalism using Eq. (4). The given uncertainties correspond to the fitted parameter's standard errors as provided from the weighted non-linear fitting procedure described in the text. in turn defined as E 0 i = µ 0 m 0,eff · H 0 i , with H 0 i representing the characteristic mean interaction field. Note that the latter quantity is determined by considering the ground state manifold of the magnetic system, i.e. the configurations for which the mean interaction field is the strongest.
We have extracted the characteristic interaction energies for the three studied arrays, presented in Table I, and compared them with the corresponding ground-state energies of a square spin ice, i.e. type I tiling, given by conventional interaction models based on the point-dipole approximation, the dumbbell representation, as well as micromagnetic simulations 34 . As shown in Fig. 4, the obtained values deviate considerably from the corresponding estimations of all the models employed, with resulting energies that can be more than 5 times smaller than expected. Furthermore, although limited by the number of experimentally-available points, the gap-dependent scaling of the extracted values appears to be much less pronounced than the model predictions. This approach was previously also employed by Morley et al. 14 -attempting to extract experimental values of the characteristic mesospin interaction energies -which reported a significant underestimation of the extracted energies.
While several factors can contribute to the mismatch between the experimental values and the various models, we mainly attribute these discrepancies to the incompatibility between the frameworks for Eq. (4) and artificial spin systems. The basic assumption of the formalism is that the VFT law can be treated as an interacting Arrhenius law, i.e. an Arrhenius-like expression in which the intrinsic energy barrier is further biased by a local interaction field: where H F i (T ) represents the temperature-dependent mean interaction field. This equality imposes a certain temperature dependence on this variable, parameterized  FIG. 4. (a) Comparison between the interaction energies per island (E 0 i ) extracted for the three experimental gap sizes (black dots) and the type I tilling ground-state interaction energies per island computed for various inverse gap sizes (1/g) using micromagnetics (blue diamonds), the dumbbell model (red line) and the point-dipole model (green line) for the given sample geometries. The shaded areas account for the error originating from a gap uncertainty of ±15 nm, extracted from the analysis of the SEM images. (b) Normalized temperature dependence of the normalized mean interaction field for various thermal dynamics ranges, calculated using the Fulcher-based expression (dotted lines) and the statistical approach (continuous lines). The black lines correspond to the case with temperature-independent energy barriers and island moments, while the coloured lines consider the temperature scaling for different normalized Curie temperatures T0/TF . The vertical dashed lines mark the boundaries of each considered range. The inset gives the normalized temperature dependence of the ratio between the Fulcher-based and statistical-based interaction fields. The experimental values provided in Table I  by E K and T F : Notice that this interaction field presents a divergence around the Fulcher temperature, a physically unrealistic feature. Taking now an at-equilibrium statistical perspective, the mean interaction field should converge at low temperatures towards a finite value corresponding to the ground-state manifold. We shall further consider the choice of Shtrikman and Wohlfarth 31 assuming a hyperbolic tangent behavior, akin to the case of a paramagnetic system in an externally applied field: This brings us to the remaining conditions for the validity of Eq. (4): • T T F , i.e. the temperature should be sufficiently far away from the Fulcher temperature.
• k B T E 0 i , i.e. the energy associated to the thermal bath should be much larger than the characteristic interaction energy.
It should also be noted that the intrinsic energy barrier is also assumed to be much higher than both the thermal bath and the interaction energy, i.e. E K k B T and E K E 0 i . With these in place, the functions describing the two interaction fields, H F i (T ) and H S i (T ), become compatible to a first order expansion. This validity region is illustrated in Fig. 4(b). Notice how the two interaction fields present almost identical temperature scaling once the temperature is about two orders of magnitude higher than T F .
If we now consider the temperature dependence of the intrinsic energy barrier and the magnetic moment, a key feature of mesoscopic spin systems, the possibility of finding a compatibility region for the two interaction fields is severely limited. Fig. 4(b) illustrates the impact of having a range in thermal dynamics bound by the Fulcher and Curie temperatures. Here, the energy barrier, E K , from Eq. (5) is replaced with its temperature dependent form given by Eq. (3), while the characteristic interaction field from Eq. (7), H 0 i , is similarly replaced with a temperature dependent expression, H 0 i · (1 − T /T 0 ) 0.5 , thus accounting for the scaling of the inter-island couplings. As it can be seen, only for a thermal range spanning several orders of magnitude can one achieve a compatibility region that accommodates the aforementioned formalism. This is particularly highlighted by the inset of Fig. 4(b), where the ratio between the two interaction field expressions is plotted as a function of temperature. The red lines correspond to the most weakly interacting sample, with an average gap of 150 nm and characterized by a T 0 /T F ∼ = 4 ratio. Even for this scenario, there is no clear overlap between the two mean interaction fields within the experimental temperature window, which should therefore compromise the matching with the interaction models considered in Fig. 4(a).
A microscopic description of the phenomenological parameter T F was provided in a study of the AC susceptibility of weakly interacting magnetic nanoparticles by Vernay et al. 22 . The VFT law, assuming the case of both weak dipolar interactions and surface anisotropy for the magnetic nanoparticles, can be transformed into a semianalytical expression, linking T F to deviations from uniaxial anisotropy and dipole-dipole interactions 22 . While this analysis assumes weak interactions being valid also in our case, the modelling of the interactions with discrete point dipoles is an imprecise description of our spatially extensive thermal mesospins (see Fig. 4(b)). On the other hand this approach 22 can serve as a stimulation for the development of a new and revised formalism capable of capturing in detail the collective temporal 5,13 and thermal dynamics 4,5,7,8 of artificial spin ice.

IV. CONCLUSIONS
We studied the AC susceptibility of thermally active square ASI arrays of varying interaction strength. The freezing of the mesospin dynamics was measured over a wide range of observation times and a systematic dependence of the freezing temperature on the inter-island coupling strength was found. Extracting magnetostatic energies from the frequency dependence using the VFT law, which was recently applied to ASI 13,14 , revealed significant discrepancies of the obtained interaction energies compared to theoretical estimates, similar to the case of Morley et al. 14 . Besides experimental uncertainties, we attribute this to the violation of the requirements of weak coupling to extract energies in the VFT formalism, along with the inability to obtain measurements at temperatures far above the freezing transition, while staying well below the material's Curie temperature. We note that these requirements are generally difficult to meet in thermally active mesospin systems. Therefore a more advanced model enabling the extraction of microscopic variables and accounting for the details of our mesospin systems, such as the internal magnetic structure, temperature dependence of the energy barriers and interaction energies, is highly desirable and will be developed in future works.
Nonetheless, AC susceptibility using a longitudinal MOKE setup was proven to be a simple yet powerful technique for studying magnetization dynamics of thermally active nanostructures. Its high sensitivity allows investigations of minute changes in the mesospin dynamics in arrays with an intrinsically low magnetic moment. The method lends itself to temperature and frequency dependent studies on a laboratory scale 20 , which facilitates an evaluation of well-established models, or scaling laws for the description of mesospin systems. The method can further serve as an excellent tool for the characterization of thermal dynamics, collective behaviour 35 and thermodynamic phase transitions in magnetic metamaterials 18,36 , such as ASIs 5,17 .
The data that support this study are available via the Zenodo repository 37 .