The Ising Model of Spin Interactions as an Oracle of Self-Organized Criticality,
Fractal Mode-Locking and Power Law Statistics in Neurodynamics

Chris King – Mathematics Department, University of Auckland

PDF

 

Abstract: This short report highlights properties of fractality and self-organized criticality in the Ising model of ferromagnetism and how these ideas can be applied using wavelet transforms to comparisons with the study of self-organized criticality in neurodynamics. Matlab programs are provided to freely replicate the results.

 

Blog Update on the Quantum Picture: Exceptionally Simple? E8, Golden Ratio and Ising Spins

 

Introduction: The Ising model [1] investigates the phase transition between ferromagnetism and paramagnetism through the Metropolis-Hastings algorithm [2] run inside a Monte Carlo loop [3]. For negative interaction strengths between spin pairs, an anti-ferromagnetic spin-glass results in which adjacent spins have lowest energy when their spins alternate up and down.

 

Fig 1: The Ising model of ferromagnetism as an example of phase transition criticality.  In the lower figure a cellular simulation based on a given cell element interacting with the four nearest neighbours (above and below and to either side) in a rectangular array. The array is iterated according to the Hamiltonian, in this case for 3000 steps. In blue magnetization is plotted as a function of interaction strength of individual spin elements, with the theoretical curve outlined in red. For low interaction strengths (or equivalently for high thermal energies by comparison with the interaction strength) we have paramagnetism in which individual domains of cells of aligned spin are randomly distributed and net magnetization will occur only in the presence of an external magnetic field. As the interaction strength increases, there is a critical phase transition to the ferromagnetic state after which the spins become polarized into large domains with aligned spins and in the asymptotic limit dominance of either spin up or spin down results in bifurcation towards a fully saturated state of net magnetization 1 or -1. At the critical transition value of J = ~0.44 magnetization domains form a fractal distribution similar to states of self-organized criticality, such as earthquakes and sand-piles, in which the critical state is maintained by fractal avalanches. In the top row interaction strength is plotted for both positive and negative interaction strengths on a larger array that shown in the lower figure.  When the interaction strength is negative spins have a lower energy when adjacent spins are opposite, so we then have the anti-ferromagnetism of a spin-glass. Negative values also have a critical transition, but this does not result in net magnetization as the lowest energy states form a chequer-board having zero net polarization. However, as there are two complementary chequer-board arrangements, large domains still develop, with frustration along their boundaries, where the energy cannot be minimized. In the insets are shown corresponding charts of energy versus magnetization and interaction strength.

 

Methods: A set of Matlab m-files performing computational Monte Carlo simulations, developed from [4] following [5] are provided in the link in appendix 2. The equations supporting the algorithm and some of the theory is displayed in appendix 1.

 

Generally the computational simulation is performed using only nearest neighbours directly above and below and to either side of an element of a rectangular array, (although more elaborate neighbourhoods are also used in figure 2). This closest neighbour computational process coincides with the Bethe-Peierls approximation [6] to Ising spin states in statistical mechanics.

 

An array is initiated in a random configuration and is iterated cell by cell, flipping a set proportion of the spins in each iteration if a random variable exceeds the energy difference between the cell and its flipped state in relation to its four neighbours. This provides a thermodynamic model in which spins will flip to a lower energy state but may also, with exponentially diminishing probability, become flipped to a higher energy one.

 

Figures 1 and 2 show the results of this investigation both for the fractal dynamics of the critical phase transition to ferromagnetism with positive interaction energies and for the capacity to form fractal mode-locked states in an anti-ferromagnetic spin-glass state with negative interaction energies. Figure 3 also shows this fractal critical behavior by analyzing the connected regions within the final state of the array.

 

Fig 2: The effect of increasing polarization on an anti-ferromagnetic spin-glass in which the lowest energy state has neighbouring cells having oppositely aligned spins. In the lower figure (a) is shown the magnetization for a 2-D spin glass in which only the four closest neighbours interact, as the external field B is increased.  There is a strong plateau in which a chequer-board lowest energy state with adjacent spins opposite is maintained. (b) There is also a hint of plateuing at a magnetization of ~0.5 in which ¼ of the spins would be aligned one way and ¾ the other. The top row arrays show states corresponding to various values of the external field. (c) The theoretical distribution for a 1-D spin-glass with negative interaction strength (Schr?der) is a ?devil?s staircase? in which mode-locking occurs for an interval neighbourhood of each rational magnetization state in which a regular periodic arrangement, such as uduud, can result in rational periods of any order. (d) The 2-D simulation failed to show any further model-locking even when the interacting neighbourhood was enlarged to include more distant interactions (with an inverse square 2-D dipole interaction), however in the corresponding 1-D simulation (with an inverse linear dipole interaction) a series of plateaus formed at rational magnetization values.

 

The negative interaction strength regime was also explored using variations on the Matlab algorithm, in both 2-D and 1-D modes, to test for the fractal devil?s staircase of mode-locked states [7] referred to in Schr?der [8]. The lack of additional plateaus was hypothesized to result from energy interaction confinement to nearest neighbours, so variants of the algorithm were developed in which more distant neighbours were also involved in the energy function through dipole interactions in their respective dimensionalities. While the 2-D model failed to demonstrate additional mode-locking states, the 1-D model did then return a series of near stationary intervals surrounding several fractional states of magnetization, consistent with the results in Schr?der.

 

Fig 3: Log-log plot of region size versus number of regions of this size range for J = 0, 0.2, 0.4406868, 0.6, 1 coloured blue, green, black, magenta and red respectively, each for an average of five Monte Carlo runs of 3000 iterations, showing that the critical state has a consistent fractal power law distribution (black) by contrast with stochastic sub-critical J which lack large regions and partitioned super-critical J which have a preponderance of a few large regions, (lower red and magenta sections) containing a small number of small random incursions (left).

 

Recent research in neurodynamics by Bullmore?s group [9] has highlighted the close correspondence between the fractal dynamics of the critical state of the Ising model, other models such as the Kuramoto model [10] and self-organized criticality in brain dynamics. We thus extended our use of the Matlab algorithm to replicate some of the methods used in their work. We thus adapted a wavelet analysis algorithm to portray the time behaviour of the Ising model in terms of wavelet amplitudes and the phase correlations between pairs of channels.

 

Fig 4: Above: Wavelet transforms of the summed outputs of an 8x8 sub-array of a 96x96 Ising simulation at the values J = 0, 0.4406868, 1 for zero external field demonstrate fractal time variation of a combined channel at the critical value. The dominance of low frequency wavelet amplitudes at J=1 reflects the partitioning into large domains. The high frequency contributions at J=0 reflect the small scale randomness. Below: Plots of phase lock of the complex correlation function between superimposed channels for the same values of J as above, showing that the fractal nature of the critical state also extends to a fractal  distribution of temporal phase lock episodes. The black regions for J = 1 result from zero wavelet amplitudes of one or other channel at higher frequencies causing the argument of the complex correlation function to be undefined.

 

In figure 4 is shown an absolute wavelet transform for an exponential series of frequencies for sub-critical, critical and super-critical interaction strengths, demonstrating that the fractal nature of the critical state can be readily detected using wavelet transforms. We defined and used a complex version of the Morlet wavelet  rather than the Hilbert wavelet transform of Bullmore?s group, but have otherwise followed their methods, using an Ising simulation on a 96x96 array and then summed the binary [-1,1] values on each 8x8 sub-array to form 12x12 - effectively 64-bit channels of ?real? amplitudes – iterated over the last 8192 iterations of a 12,192 iteration run.

 

The successive iterates were firstly given a direct wavelet portrait using the absolute amplitude of the real part of the wavelet coefficient array to generate the upper series of profiles in figure 4.  Pairs of channels were then phase correlated and the phase angle of the resulting complex argument taken in the range  to generate a set of two-channel correlation profiles, giving a portrait of the fractal distribution of the critical state in terms of time intervals of phase correlation between a pair of channels.

 

 

Fig 5: Log-log plot of number of intervals on which a pair of channels have phase lock with respect to the length of the phase lock. J = 0, 0.4406868, 1 coloured blue, black, and red respectively. The critical state (black) approximates a power law distribution, while J = 0 lacks longer phase lock intervals and J=1 is dominated by large intervals corresponding to the large domains and very short intervals caused by frustration between competing domains.

 

Finally we pooled phase correlation intervals of several pairs of channels using Bullmore group?s criteria to test for a power law relationship at criticality over 4 orders of magnitude of interval length (exponentially double the range of fig 6). These results were not averaged over short time intervals as in their work, and were performed only for a small number of pairs, but do show an approximate power law distribution in the log-log plot of figure 5 as a proof of principle for the Matlab program which can be compared with their results as shown in figure 6.

 

Fig 6: (Left): Probability distribution of phase lock interval (PLI) between pairs of processes at critical (black line)

and at hot temperature (low interaction strength)  (blue) plotted on a log-log scale. The black dashed line represents a power law with slope a~{1:5. (Right) Probability distribution of lability of global synchronization (the number of channel pairs in phase lock over a given interval) at critical temperature (black line) and at hot temperature (low interaction strength) (blue); the black dashed line represents a power law with slope a~{0:5) (Bullmore et. al.).

 

Appendix 1: Summary of the Theory and Equations

 

The Hamiltonian for the interaction is .

The ratio of probabilities for a flip is . The spins of particles in which r > 1 or r < 1 but greater than a uniformly distributed random number have the potential to be flipped if a second random number exceeds a given threshold regulating the proportion flipping in each iteration.

 

The values for the theoretical thermodynamic curves shown in figure 1 are outlined below:.

 

 

The wavelet coefficients  the phase coherence  and the lability of global synchronization  are given by:

Appendix 2: Matlab Programs Download Link

 

Ising.zip

 

 

References