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

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, ${\bf{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 neighbors, as seen in figure 1. The total energy of such a two-dimensional Ising system can be described by the equation

$$\begin{eqnarray}&&E=-{J}_{h}\displaystyle \sum _{i,j}^{h}{{\bf{s}}}_{i}\cdot {{\bf{s}}}_{j}-{J}_{v}\displaystyle \sum _{i,j}^{v}{{\bf{s}}}_{i}\cdot {{\bf{s}}}_{j},\end{eqnarray} \tag{ 1 }$$

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 neighbor spins ${{\bf{s}}}_{i}$ and ${{\bf{s}}}_{j}$. As opposed to the one-dimensional model which shows no spontaneous magnetization at temperatures $T\gt 0$K a spontaneous magnetization appears in the two-dimensional case [6] with an order parameter given by [7, 8]

$$\begin{eqnarray}&&M(T)={\left(1-{\left[\mathrm{sinh}\left(\displaystyle \frac{2{J}_{h}}{{k}_{{\rm{B}}}T}\right)\mathrm{sinh}\left(\displaystyle \frac{2{J}_{v}}{{k}_{{\rm{B}}}T}\right)\right]}^{-2}\right)}^{\frac{1}{8}},\end{eqnarray} \tag{ 2 }$$

where T corresponds to the temperature and ${k}_{{\rm{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].

Zoom In Zoom Out Reset image size

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 with in plane [11–13] and perpendicular realizations [14]. The two-dimensional nano-scale nature of these systems enables their state to be directly determined by imaging techniques such as magnetic force microscopy (MFM) [11, 15], photoemission electron microscopy [16–18] and Lorentz microscopy [19].

In this paper we investigate a two-dimensional array composed of elongated thin film islands in a square lattice (see figure 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 [15]. 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.

Zoom In Zoom Out Reset image size

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 [15, 20]. 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 behavior 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 figure 3). Subsequently the array can be imaged by MFM and the magnetization direction of each island can be determined (see figure 4).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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 neighboring islands are of the order of the thermal energy ${k}_{{\rm{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 [21]. For this study we therefore choose amorphous Co₆₈Fe₂₄Zr₈ as the island material previously used for creating ultra-thin magnetic layers [22] and well defined nano-patterned multilayered structures [23]. Furthermore, a field imprinted anisotropy can be induced in Co₆₈Fe₂₄Zr₈ enhancing the energy barrier for reversal as the magnetization has settled in a fixed direction [24].

Considering magnetostatic interactions in the point dipole approximation between neighboring islands reveals an interaction scheme which can be mapped to a ferromagnetic interaction in the vertical direction ${J}_{v}\gt 0$ and antiferromagnetic in the horizontal direction ${J}_{h}\lt 0$. The lowest energy state of the array is therefore composed of a staggered arrangement of ferromagnetically aligned chains (see figure 2). The lowest energy state of an asymmetric two-dimensional Ising system is twofold 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 figure 4(c) and the energy state of individual island can be categorized into nine different energy states of varying degeneracy shown in figure 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 figure 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 (figure 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 $\sim \mathrm{exp}(-E/{k}_{{\rm{B}}}T)$ and that the system can be sufficiently described by a nearest neighbor interaction model.

Zoom In Zoom Out Reset image size

Identifying the domain walls separating the twofold degenerate ground states of the array facilitates the mapping of the magnetic structure of the system as two ordering states, as illustrated in figure 3. The mapping of these two states onto two domain colors for the entire array is shown in figure 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\pm 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.

Zoom In Zoom Out Reset image size

By calculating the pairwise correlation between spins, at a distance ${\bf{r}}$, within the experimental array the correlation function for the system, $G({\bf{r}})$, can be determined. The resulting correlation array is shown in figure 7. The preference for an antiferromagnetic arrangement of neighboring 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({\bf{r}})|$ as determined from the array. As can be seen in figure 7 the pairwise correlation follows an exponentially decreasing function as expected for an extended array of Ising like spins at $T\gt {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 ${\xi }_{[10]}=2.87$ and ${\xi }_{[01]}=1.02$. The correlation length in the [11] direction, ${\xi }_{[11]}=0.97$, is revealed to be similar to ${\xi }_{[01]}$.

Zoom In Zoom Out Reset image size

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 [25] or utilizing numerical calculations such as Monte Carlo simulations, see figure 8. In the case of an isotropic system, wherein the interaction strength between neighbors, J, is the same in the two main lattice directions ([01] and [10]), the ordering temperature is given by ${T}_{c}=\frac{2J/{k}_{{\rm{B}}}}{\mathrm{ln}(1+\sqrt{2})}\approx 2.269J/{k}_{{\rm{B}}}$. As J_v decreases with respect to J_h the relative T_c of 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 figure 8 the results of numerical simulations correspond well to the analytical solution with small deviations arising from the finite nature of the simulated array.

Zoom In Zoom Out Reset image size

Assuming a system in a perfectly ordered ground state at T>0 K as the temperature is increased thermal energy is introduced into the system and spins will start to change 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({\bf{r}})$ between spins. For the two-dimensional Ising model the correlation decays exponentially with distance towards zero for $T\gt {T}_{c}$ while for $T\lt {T}_{c}$ it decays exponentially towards the square of the spontaneous magnetization per spin [26]. The correlation function can be written as $G({\bf{r}})\sim \mathrm{exp}\left(-\frac{{\bf{r}}}{\xi }\right)$ for $T\gt {T}_{c}$, where ξ is the correlation length affected by temperature and the interaction strength between neighboring spins. The inset in figure 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}\gt {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 ${\xi }_{[10]}$ and ${\xi }_{[01]}$ follow the same trend as the experimental results, i.e. ${\xi }_{[10]}\gt {\xi }_{[01]}$ and ${\xi }_{[01]}\sim {\xi }_{[11]}$. In particular at $T=1.9{J}_{h}/{k}_{{\rm{B}}}$ the values obtained from the simulations are ${\xi }_{[10]}=1.95$, and [01], ${\xi }_{[01]}=1.06$, ${\xi }_{[11]}=1.08$, similar to the values obtained experimentally.

The interaction strength between neighboring islands was calculated using micromagnetic modeling using rigid spins and assuming a simple dipole model of the macrospins for neighbors up to three lattice constants (see further details in supplementary material). The resulting interaction energies were then used in conjunction with MC simulations to include long-range interactions. However, when including these long-range interactions the correlation lengths were not observed to behave as it is in the experiment, with the [01] direction showing larger correlation lengths than the [10] direction, in contradiction with the experimental observations. This lack of agreement between the micromagnetic analysis and the empirically analyzed array indicates the need for carrying out a detailed analysis of the stray field landscape for the experimental system during the transition from a thermally active state to an arrested state as any modifications to the magnetization within the islands away from a rigid spin state would introduce higher order magnetostatic interactions. The energy landscape is furthermore modified by the arrangement and density of the islands as flux lines are screened by nearer neighbors. Such modifications would affect the energy landscape of the array favouring an increase in energy for nearest neighbors while reducing longer range interactions. From the data we deduce that the magnetic state of the experimental array corresponds to a temperature above the ordering temperature of the system i.e. that the thickness of the islands when the magnetization becomes frozen in is lower than the thickness where the system would undergo the phase transition. However if these thicknesses are similar, effects due to critical slowing down at the phase transition can not be neglected and might affect the final observed state [27]. The good agreement obtained from the nearest neighbor approach applied in the manuscript supports the case that this system is sufficiently described by a nearest neighbor Ising model of dipolar interactions. These results, therefore, indicate that in the arrested state the experimental observations can be described with an effective nearest neighbor Ising model. Hence, long-range dipolar interactions are not needed to model the present system in the arrested stated. However, it should be noted that long-range dipolar interactions cannot be 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 [28] such as artificial spin ice, one- and two-dimensional Ising arrays and different frustrated arrangements [29, 30] 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.

Acknowledgments