A new look on the two-dimensional Ising model: thermal artificial spins

We present a direct experimental investigation of the thermal ordering in an artificial analogue of an asymmetric two dimensional Ising system composed of a rectangular array of nano-fabricated magnetostatically interacting islands. During fabrication and below a critical thickness of the magnetic material the islands are thermally fluctuating and thus the system is able to explore its phase space. Above the critical thickness the islands freeze-in resulting in an arrested thermalized state for the array. Determining the magnetic state of the array we demonstrate a genuine artificial two-dimensional Ising system which can be analyzed in the context of nearest neighbour interactions.

We present a direct experimental investigation of the thermal ordering in an artificial analogue of an asymmetric two dimensional Ising system composed of a rectangular array of nano-fabricated magnetostatically interacting islands. During fabrication and below a critical thickness of the magnetic material the islands are thermally fluctuating and thus the system is able to explore its phase space. Above the critical thickness the islands freeze-in resulting in an arrested thermalized state for the array. Determining the magnetic state of the array we demonstrate a genuine artificial twodimensional Ising system which can be analyzed in the context of nearest neighbour interactions. The Ising model invented by Wilhelm Lenz and solved in one dimension by Ernst Ising in 1924 is one of the pillars of statistical mechanics [1,2]. Although built on a simple basis, that of an interacting system composed of a chain of entities with only two discrete states, s = {1, −1}, the Ising model still to this day is used to model magnetic systems and can be applied to a wealth of atomistic and mesoscopic experimental systems ranging from ferromagnetic ordering and atomic-scale antiferromagnets [3] to the ordering of binary colloidal structures [4] and thermal artificial spin systems [5]. In the two dimensional case spins are arranged on a square lattice and each spin interacts with four nearest neighbours, as seen in Fig. 1. The total energy of such a two-dimensional Ising system can be described by the equation where J h and J v correspond to interaction energies in the horizontal [10] and vertical [01] lattice directions of the two dimensional crystal and the sum is taken over all pairs of nearest neighbour spins s i and s j . As opposed to the one dimensional model which shows no spontaneous magnetization at temperatures T > 0 a spontaneous magnetization appears in the two dimensional case [6] with an order parameter given by [7,8] where T corresponds to the temperature and k B is the Boltzmann constant. The model is not complicated by the choice of different values or signs of J h and J v [7]. For J > 0 the interaction between neighbouring spins favors a parallel alignment (see Fig. 1 of J < 0 spins have a preference for an antiferromagnetic alignment. Nano-patterned single-domain magnetic thin film islands have been a prominent candidate in recent years for creating artificial analogues of interacting systems. Using modern lithographic techniques it has become feasible to design directly the shape of such islands and their geometrical arrangement creating two dimensional arrays of interacting artificial spins. Artificial Ising-like spins can be realized by designing elongated islands of thin films materials, for which shape anisotropy confines the magnetization to only two possible orientations. By arranging such islands in different geometries a wealth of interacting systems can be studied including cellular automata [9,10] and frustrated systems such as artificial spin ice [11,12]. The two dimensional nano-scale nature of these systems enables their state to be directly determined by imaging techniques such as mag- netic force microscopy [11,13], photoemission electron microscopy [14][15][16] and Lorentz microscopy [17]. In this paper we investigate a two dimensional array composed of elongated thin film islands in a square lattice (see Fig. 2). During fabrication the array goes through a dynamic phase enabling the system to thermally explore its phase space leading to a low energy ordered state [13]. After this dynamic phase the state of the system becomes frozen-in, generating a snapshot of an arrested thermalized state. The determination of the magnetization of individual islands in the array allows us to demonstrate that the dynamic phase leads to an ordering of the magnetic state of the array, which can be described by the two-dimensional Ising model.
The thickness regime wherein thermal dynamics of the spins are obtained occurs during the deposition of the magnetic material onto prepatterned substrates as shown by Morgan et al. [13,18]. During the deposition, the island thickness (and thereby their volume) becomes gradually larger and the magnetization reversal energy barrier, E r , associated with their shape anisotropy increases. In the initial stages of the island growth this barrier is smaller than the thermal energy enabling the magnetization to spontaneously fluctuate between the two low energy states defined by the shape anisotropy. As the thickness becomes larger the reversal energy increases, eventually overcoming the thermal energy, reaching a threshold where the superparamagnetic behaviour is suppressed as the scale of the combined shape anisotropy reversal energy barrier and the energy landscape of the array are sufficient to lock the magnetization in each of the islands (see Fig. 3). Subsequently the array can be imaged by magnetic force microscopy (MFM) and the magnetization direction of each island can be determined (see Fig.  4).
During the limited time window where fluctuations occur the fast dynamics allow the array to explore its phase space and achieve an equilibrium condition. During this dynamic phase before the magnetization becomes frozen the reversal energy barrier and the interaction energies between neighbouring islands is of the order of the thermal energy k B T . For the island sizes and array parameters used for this study the energy values involved correspond to room temperature values for film thicknesses in the ∼1 nm regime. Uniform thermal dynamics over the entire array therefore require a well defined, stable thickness of each island in order to minimize randomization effects due to variations in the island thickness and film roughness. Amorphous magnetic materials display soft magnetic properties and a high degree of structural uniformity and are thereby suitable for well defined layers below 1 nm in thickness [19]. For this study we therefore choose amorphous Co 68 Fe 24 Zr 8 as the island material previously used for creating ultra-thin magnetic layers [20] and well defined nano-patterned multilayered structures [21]. Furthermore, a field imprinted anisotropy can be induced in Co 68 Fe 24 Zr 8 enhancing the energy barrier for reversal as the magnetization has settled in a fixed direction [22].
Considering magnetostatic interactions in the point dipole approximation between neighbouring islands reveals an interaction scheme which can be mapped to a ferromagnetic interaction in the vertical direction J v > 0 and antiferromagnetic in the horizontal direction J h < 0. The lowest energy state of the array is therefore composed of a staggered arrangement of ferromagnetically aligned chains (see Fig. 2). The lowest energy state of an asymmetric two dimensional Ising system is two degenerate, in our case corresponding to a ferromagnetic ordering in the vertical direction and an antiferromagnetic ordering in the horizontal direction. Excitations above the ground state occur through the reversal of a macrospin and can be viewed in the form of boundary walls separating the two possible ground states as illustrated in Fig. 4(c) and the energy state of individual island can be categorized into 9 different energy states of varying degeneracy shown in Fig. 4(d).
Counting the abundances of excitations composed only of independent vertical or horizontal excitations the relative energy scale between the two directions can be attained. The observed probabilities of these excitations decreases exponentially with increasing energy (inset in Fig. 5) in accordance with a Boltzmann distribution of the states. Determining the ratio of the excitation energies from the inset the energies relating to all energy values of the individual islands can be listed in units of the energy involved with a horizontal excitation |J h |. Within the combined plot (Fig. 5) the observed abundances decrease exponentially establishing the probability for a macrospin to be in an energy state E to be given by a Boltzmann distribution ∼ exp (E/k B T ) and that the system can be sufficiently described by a nearest neighbour interaction model.
Identifying the domain walls separating the two degenerate ground states of the array facilitates the mapping of the magnetic structure of the system as two ordering states, as illustrated in Fig. 3. The mapping of these two states onto two domain colours for the entire array is shown in Fig. 6. Counting the number of islands falling into each of these two domains the order parameter of the array can be obtained. The resulting order parameter for the array, defined by M = (n b − n w )/(n b + n w ) where n b and n w correspond to the domain populations of black and white domains, can then be obtained. Truncating the data array to a square shape the array size is reduced to n = 9828 out of which n b = 5283 islands fall into the black domain while n w = 4545 islands fall into the white domain. The resulting value of M = 0.075 ± 0.015 reveals a slightly higher population of the black domains. Although the statistics of this value are limited by the finite number of observed islands in the array it could be indicative of the array being in a state close to, or above T c , since the lack of global order in the macroscopic arrangement can lead to a finite value of the order parameter. The listed uncertainty of M is determined from the square root of the domain populations, n w and n b . By calculating the pairwise correlation between spins, at a distance r, within the experimental array the correlation function for the system, G(r), can be determined. The resulting correlation array is shown in Fig.  7. The preference for an antiferromagnetic arrangement of neighbouring spins in the [10] direction introduces the possibility of negative values in the correlation and alternating positive and negative values along the [10] direction. Figure 7 therefore shows the absolute value of the pairwise correlation |G(r)| as determined from the array. As can be seen in Fig. 7 the pairwise correlation follows an exponentially decreasing function as expected for an extended array of Ising like spins at T > T c . Furthermore, the asymmetric nature of the array is revealed in the different values of the correlation lengths along the [10] and [01] directions of the array with ξ [10] = 2.87 and ξ [01] = 1.02. The correlation length in the [11] direction, ξ [11] = 0.97, is revealed to be similar to ξ [01] .
The role of the ratio of the interaction strengths on the order parameter and specifically the ordering temperature of the array can be investigated using Onsager's solution, as initially done by Chang [23] or utilizing numerical calculations such as Monte Carlo simulations, see Fig. 8. In the case of an isotropic system, wherein the interaction strength between neighbours, J, is the same in the two main lattice directions ([01] and [10]), the ordering temperature is given by T c = 2J/kB ln(1+ √ 2) ≈ 2.269J/k B . As J v decreases with respect to J h the relative T c of Macrospin array domain configuration. The macrospin configuration of the investigated arrays mapped to black and white domains corresponding to the two degenerate ground state of the system. the system also decreases. Figure 8 shows results from Monte Carlo simulations for decreasing values of the ratio |J v /J h |. The change in T c corresponds well to what is expected in the two dimensional Ising model as the probability of observing two parallel spins in the [01] direction decreases. As can be seen in Fig. 8 the results of numerical simulations correspond well to the analytical solution with small deviations arising from the finite nature of the simulated array.
Assuming a system in a perfectly ordered ground state at 0 K as the temperature is increased thermal energy is introduced into the system and spins will start to change FIG. 7. Pairwise correlations G for the array. The graph shows the absolute value, |G|, for easier mapping of the pairwise correlations for both the ferromagnetic and antiferromagnetic directions. From the slopes the resulting correlation lengths along the different directions is determined to be ξ [10] = 2.87, ξ [01] = 1.02, and ξ [11] = 0.97 their directions. The perfect initial order is therefore broken and the system moves towards a thermally disordered state with a reduced order parameter. From numerical simulations this transition can, furthermore, be observed through the pairwise correlation G(r) between spins. For the two dimensional Ising model the correlation function can be written as G(r) ∼ exp(− r ξ ) for T < T c and T > T c where ξ is the correlation length affected by temperature and the interaction strength between neighbouring spins. The inset in Fig. 8 shows the temperature dependence of ξ for the three major directions of the array, [10], [01], and [11] for the asymmetry ratio |J v /J h | = 0.3 as determined from Monte Carlo calculations. As J h > J v the corresponding correlation length in the [10] direction is larger than in the other directions while the correlation length in the [01] and [11] directions are similar. Considering the results of the Monte Carlo simulations for an interaction ratio of |J v /J h | = 0.3 one can see that for temperatures larger than T c the correlation lengths ξ [10] and ξ [01] follow the same trend as the experimental results, i.e. ξ [10] > ξ [01] and ξ [01] ∼ ξ [11] . In particular at T = 1.9J h /k B the values obtained from the simulations are ξ [10] (2) (solid lines) fits quite well with the obtained numerical data. The inset shows the correlation length ξ as a function of temperature for the ratio |Jv/J h | = 0.3 for the three major crystallographic directions of the array, [10], [01], and [11]. be Representative snapshots of the numerical simulations at temperatures of 1.0J h /kB, 1.34J h /kB, 1.5J h /kB, and 2.0J h /kB, respectively, for an asymmetry ratio of |Jv/J h | = 0.3. Due to the asymmetry the correlation length in the [10] direction is larger than in the [01] direction stretching out the domains in the [10] direction, similar to what is observed in the experimental data.
completely excluded if one desires a full description of the system in any non-arrested state. Further advances in materials science and magnetic imagining techniques will undoubtedly allow us to follow the thermal evolution of a multitude of different artificial spin systems [24] such as artificial spin ice, one and two dimensional Ising arrays and different frustrated arrangements [25,26] as they explore their phase space. Such investigations do not merely offer a direct real space probe of the microstate of thermal systems at a local scale but also a direct determination of the dynamics involved. The possibility of determining directly the state of these systems enables the experimental probing of the effect of external variables such as temperature and applied magnetic field. This introduces the possibility of investigating e.g. the two-dimensional Ising model under applied external field directly for which an analytical solution does not exist.

Methods
The island array was prepared by magnetron sputtering thin film growth onto a pre-patterned fused silica substrate. The 2×2 mm 2 pre-patterned array was fabricated by surface conformal nano-imprint lithography [27,28]. The array was composed of 470 nm × 160 nm islands arranged in a square lattice with a periodicity of 200 nm along the short axis of the islands ([10] direction) and a periodicity of 500 nm along the [10] direction. The amorphous film was deposited by magnetron sputtering from a compound target at room temperature. Argon was used as the sputtering gas at a pressure of 3.0 mTorr. Before growth the base pressure of the sputtering system was below 2 × 10 −7 Pa. Initially a 2 nm thick layer of Al 70 Zr 30 was grown onto the pre-patterned substrate followed by the growth of a 7 nm thick Co 68 Fe 24 Zr 8 and finally a 2 nm thick Al 70 Zr 30 layer to prevent degradation. After patterning the structural quality of the films was investigated by AFM carried out in contact mode using a Nanosurf Mobile S instrument. The thermally ordered state of the array was investigated by magnetic force microscopy using a digital instrument dimension 3100. Topography was imaged by height contrast in tapping mode and magnetic micrographs were scanned in lift mode with a lift scan height of 70 nm. In order to reduce the risk of the stray field from the MFM tip altering individual island states the images were recorded using a low magnetic moment Co alloy MFM tip (PPP-LM-MFMR). Repeated scanning of the same area did not alter the state (see supplementary material). The system was modeled as an Ising system, in which every nano-island is considered as a single macrospin with only two possible orientations (+1 and −1), the simulation cell consists of 150×150 spins. Simulations with both periodic and open boundary conditions were performed. The interactions were only considered between nearest neighbours. The properties of the system were studied via Monte Carlo simulations using the Metropolis-Hastings algorithm. To simulate the process of the spins freezing at a given configuration and then evolving towards the most stable configuration before any measurement is performed the system is considered to be completely disordered at a temperature higher than its critical temperature. T c . During this thermalization phase no measurements are performed, the system is allowed to evolve in accordance to the usual Metropolis-Hastings algorithm. After a certain number of Monte Carlo steps the system is cooled down and the process is repeated once more. This procedure is then repeated until the desired measurement temperature is reached.