flatspin: A large-scale artificial spin ice simulator

We present flatspin, a novel simulator for systems of interacting mesoscopic spins on a lattice, also known as artificial spin ice (ASI). A generalization of the Stoner-Wohlfarth model is introduced, and combined with a well-defined switching protocol to capture realistic ASI dynamics using a point-dipole approximation. Temperature is modelled as an effective thermal field, based on the Arrhenius-Néel equation. Through GPU acceleration, flatspin can simulate the dynamics of millions of magnets within practical time frames, enabling exploration of large-scale emergent phenomena at unprecedented speeds. We demonstrate flatspin's versatility through the reproduction of a diverse set of established experimental results from literature. In particular, the field-driven magnetization reversal of “pinwheel” ASI is reproduced, for the first time, in a dipole model. Finally, we use flatspin to explore aspects of “square” ASI by introducing dilution defects and measuring the effect on the vertex population.

Figure 1

The representation of nanomagnets as spins $s_i$ and associated quantities: position ${\mathbf{r}}_i$, angle ${\theta }_i$, and distance to neighbor $j$, ${\mathbf{r}}_{ij}$. Note that the magnetization of spin $i$ is given by its spin $s_i$, and orientation ${\theta }_i$.

Figure 2

Probability density functions (PDF) for $h_{\text{th}}$ at different temperatures $T$. The example plots are for a magnetic particle with saturation magnetization $M_{\text{S}}=860\mathrm{kA}{\mathrm{m}}^{-1}$ and volume $V=220\mathrm{nm}\times 80\mathrm{nm}\times 3\mathrm{nm}$, for a time interval $\mathrm{\Delta }t=1\mathrm{ms}$.

Figure 3

Top: Schematic showing hard and easy axes of (a) an elliptical magnet and (b) a rectangular stadium-shaped magnet, as well as the total field acting on the magnet ${\mathbf{h}}_i$ with its parallel and perpendicular components, ${\mathbf{h}}_{\parallel }$ and ${\mathbf{h}}_{\perp }$, respectively. Bottom: Switching astroid for (c) an elliptical magnet and (d) a rectangular stadium-shaped magnet. Red dots show the coercive field obtained from micromagnetic simulations. The blue line in (c) shows the Stoner-Wohlfarth astroid. The blue line in (d) shows the generalized Stoner-Wohlfarth astroid with parameters $b=0.42$, $c=1$, $\beta =1.7$, and $\gamma =3.4$ in Eq. (11). The astroids have been normalized with respect to $h_k$.

Figure 4

flatspin includes the most common ASI geometries: (a) Square (closed edges), (b) square (open edges), (c) kagome, (d) pinwheel “diamond”, (e) pinwheel “lucky-knot”, and (f) Ising.

Figure 5

Overview of the flatspin architecture, with arrows indicating data flow.

Figure 6

M(H) curves of ensembles of noninteracting magnets at different temperatures, simulated with different approaches, and for different coercivities. Also indicated is the analytical $tanh\left(A{\mu }_{0}H\right)$ where $A=M_{\text{S}}V/k_{\text{B}}T$. Note that the MuMax3 results are for binned cell magnetization.

Figure 7

Left: Snapshots of the evolution of a kagome ASI at selected field values. the images have been cropped to show the middle $80\%$ of the total ASI to improve clarity. Right: Comparison of the hysteresis curve of the simulated ensemble (blue line) against a sketch of the hysteresis curve from the experimental results of Mengotti et al. [34] (red-dashed line). The labeled points indicate the points at which the snapshots are sampled from.

Figure 8

(a) Maps showing the net magnetization of the ASI vertices, for an annealed $50\times 50$ square ASI with the given lattice spacing $a$. The white regions have zero net magnetization, and thus correspond to a coherent domain of type I vertices. Colored regions have a nonzero net magnetization, direction indicated by the color wheel, and correspond to type II or III vertices. (b) The absolute value of spin-spin correlation at a given separation for square ASIs of different lattice spacings, $a$. Dashed curves show exponential best fits for data from the original paper [10]. Also indicated is a $1/e$ threshold of correlation.

Figure 9

Hysteresis loop and snapshots of the pinwheel units for various angles $\theta$ of the applied field. Figures (a), (c), and (e) show experimental results, adapted from Li et al. [13], Copyright ©2018 American Chemical Society, CC-BY-4.0 . Figures (b), (d), and (f) show results from flatspin simulation.

Figure 10

(a) Spearman's rank correlation coefficients $\rho$ averaged over 32 different square ASIs, evaluated for different lattice spacings in flatspin and MuMax3. The red line shows the approximate maximum ridge line through the heatmap. (b) The red line shows the true maximum $\rho$ for the lattice spacing pairs. The blue and violet lines show projections of the top row and rightmost column of (a), respectively.

Figure 11

The throughput (number of field calculations per second) as a function of number of spins. Throughput is averaged over 100 simulations of each size. The test was performed on an NVIDIA Tesla V100 GPU with 32 GB of RAM. Note the logarithmic scale of the axes.

Figure 12

A snapshot from flatspin simulations of a pinwheel ASI system with more than one million magnets, as it undergoes reversal by an external field. The angle of the external field is $\theta =0^{\circ }$.

Figure 13

(a) Average complete vertex fraction, by type, as a function of dilution by removal of spins in a $50\times 50$ square ASI. Only complete vertices are included. The averages are taken over 10 different instances for each dilution fraction. (b) Example end states of square ASI of different dilutions. Each pixel correspond to the net magnetization of 4 spins (vertices), as in Fig. 8. Note that white pixels correspond to vertices of no net magnetization (type I vertices if the vertices are complete, as there are no type IV vertices).