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). Our magnetic switching criteria enables ASI dynamics to be captured in a dipole model. Through GPU acceleration, flatspin can simulate realistic dynamics of millions of magnets within practical time frames. We demonstrate flatspin's versatility through the reproduction of a diverse set of established experimental results from the literature. In particular, magnetization details of"pinwheel"ASI during field-driven reversal have been reproduced, for the first time, by a dipole model. The simulation framework enables quick exploration and investigation of new ASI geometries and properties at unprecedented speeds.


I. INTRODUCTION
An artificial spin ice (ASI) is an ensemble of nanomagnets arranged on a lattice, coupled through magnetic dipole-dipole interactions. The vast variety of emergent collective behaviors found in these systems have generated considerable research interest over the last decade [10,26]. Using modern nanofabrication techniques, emergent phenomena can be facilitated through direct control of the ASI geometry, e.g., collective ferromagnetic/antiferromagnetic ordering [27], Dirac strings [19], and phase transitions [1,15]. ASIs offer a unique model system for exploring fundamental physics, since magnetic microscopy enables direct observation of their internal state. There is also a growing interest in ASIs as building blocks for novel devices [11].
Micromagnetic simulations of ASI have been limited to a handful of nanomagnets due to excessive computational cost. Although physically accurate, such high fidelity simulations are unable to capture large-scale emergent phenomena, such as the size of magnetically ordered domains and long-range order. To simulate large ASI systems, an established approach is to sacrifice fidelity for speed by employing a dipole model, i.e., treating each nanomagnet as a single macro spin approximated by a point dipole [3]. Traditionally, Monte Carlo methods have been used in conjunction with the dipole approximation to search for low energy configurations [12,32] or study statistical measures such as vertex populations [3]. However, Monte Carlo methods are inherently stochastic and better suited for ensemble statistics rather than detailed dynamics [2].
flatspin is a simulator for large ASI systems based on a dipole approximation with the ability to capture realistic dynamics, inspired by the work of Budrikis [3]. The flexibility of flatspin enables quick exploration of ASI parameters and geometries. Through GPU acceleration, flatspin can capture realistic dynamics of millions of magnets within practical time frames. * johannes.jensen@ntnu.no In this paper, we present the motivation and design of flatspin. We demonstrate good agreement between flatspin and a variety of published experimental results. We show that flatspin can capture dynamic behaviors observed experimentally, which have previously eluded modeling [16].

II. THE FLATSPIN MAGNETIC MODEL
In this section, we describe the dipole model and the underlying physical assumptions of flatspin. The model is designed to simulate the ensemble state-by-state evolution, i.e., dynamics, of two-dimensional ASI. In short, magnets are modeled as point dipoles (section II A), and each dipole can be affected by three types of external influence: magnetic dipole-dipole coupling (section II C), an applied external magnetic field (section II D), and thermal fluctuations (section II E). The switching of spins is determined using a generalized Stoner-Wohlfarth model (section II F). Imperfections in the ASI are introduced as different coercive fields, set per spin (section II G). Dynamics are modeled using a deterministic single spin flip strategy (section II H).
A. Magnets as dipoles ASI systems are physically realized as elongated islands of a ferromagnetic material, arranged on a twodimensional lattice. The magnets are made small enough to exhibit a single ferromagnetic domain, i.e., coherent magnetization throughout the magnet. The single domain is stable as the energy cost associated with domain walls exceeds the cost associated with the demagnetization energy [7,14]. Since a magnet has coherent magnetization, it can be approximated by a single mesoscopic spin and the magnetic state can be represented by a single vector m.
The magnets will exhibit an in-plane shape anisotropy due to their small thickness, as well as a magnetization direction defined by their elongated shape. Hence indi- Note that the magnetization of spin i is given by its spin, si, and orientation, θi.
vidual magnets can be approximated by classical macro spins with a twofold degenerate ground state defined by the elongated shape of the individual elements. Due to the two degenerate ground state configurations, we approximate each magnet as a magnetic dipole with binary magnetization, i.e., a macro spin, s i ∈ {−1, +1}.
As illustrated in Fig. 1, each magnetic dipole is modelled with a position r i and rotation θ i , which together define the ASI geometry. Furthermore, each magnet is assigned a coercive field, h (i) c , describing its resistance to switching (see section II F). Using reduced units, the magnetization vector of a single magnet can be expressed as wherem i is the unit vector along m i .

B. Magnetic Fields and Temperature
External fields and temperature are modeled as a combination of effective magnetic fields. The total field, h i , affecting each magnet i is the sum of three components: where h (i) dip is the local magnetic field from neighboring magnets (magnetic dipole-dipole interactions), h (i) ext is a global or local external field, and h (i) th is a random magnetic field representing thermal fluctuations in each magnetic element. The contributions from each field is described in the following sections.

C. Magnetic dipole-dipole interactions
The individual magnets, or spins, are coupled solely through dipole-dipole interactions. Each spin, i, is subject to a magnetic field from all neighboring spins, j = i, given by where r ij = r i − r j is the distance vector from spin i to j, and α scales the dipolar coupling strength between spins. The coupling strength α is given by α = µ0M 4πa 3 , where a is the lattice spacing, M is the absolute magnetic moment of a single magnet, and µ 0 is the vacuum permeability. The distance r ij is given in reduced units of the lattice spacing.
The dipole field present at each spin's location is calculated by summing the dipole field contributions from spins in its neighborhood. The size of the neighborhood is user-configurable and defined in units of the lattice spacing. In some geometries, such as square ASI, short range interactions dominate the contributions to h dip [4,31], in which case the neighborhood size can be relatively small, for a benefit of increased efficiency. For geometries where long range interactions are significant, a larger neighborhood is required, e.g., pinwheel ASI [17].

D. External field
Applying an external magnetic field is the primary mechanism for altering the state of an ASI in a controlled manner. The external field can either be set locally on a per-spin basis, h (i) ext , globally for the entire system, h ext , or as a spatial vector field, h ext (r).
Time-dependent external fields are supported, i.e., h ext is a discrete time series of either local, global or spatial fields. A variety of time-dependent external fields are provided, including sinusoidal, sawtooth and rotational fields. More complex field-protocols can be generated, e.g., for annealing purposes or probing dynamic response.

E. Thermal field
Thermal fluctuations are modeled as an additional local field, h (i) th , applied to each magnet individually. Two orthogonal components of the field are independently drawn from the Normal distribution N (0, σ 2 th ). The simulated temperature T is closely related to the value of σ th , where a large σ th corresponds to higher temperatures.
When the material and geometric properties of the magnetic islands are known, it is possible to choose a σ th to match a desired thermal behavior. In some cases, such as for no thermal field and for thermal protocols with no absolute reference point, it is useful to note that σ th = 0 implies T = 0, and that T is a monotonically increasing function of σ th .

F. Switching
Magnetization reversal, or switching, may take place when a magnet is subjected to a magnetic field or high temperature. If the field is sufficiently strong (stronger than some critical field) and directed so that the projection onto m i is in the opposite direction to m i , the magnetization m i will switch direction.
The critical field strength is referred to as the coercive field h c . For elongated magnets, h c depends on the angle between the applied field h i and m i . As illustrated in Fig. 2a, the easy axis, where the magnetization favors alignment, lies along the long axis of the magnet, whereas the hard axis is perpendicular to the long axis. The external field can be decomposed into two components, h and h ⊥ , corresponding to the field component parallel and perpendicular to the easy axis, respectively. We denote the coercive field strength along the hard axis as h k .
A switching astroid is a polar plot of h c at different angles, with h ⊥ on the horizontal axis and h on the vertical axis. It is a description of h c , or h c (φ). For any applied field h i that is outside the switching astroid, the magnet will switch (given that the projection of h i onto m i is oppositely aligned with respect to m i ). Fig. 2c shows the normalized switching astroid for an elliptical magnet (Fig. 2a) as obtained from micromagnetic simulations (red dots). Notice how h c is the smallest at a 45 • angle, indicating that a field directed at 45 • to a magnet's principal axes will require the least field strength in order to switch its magnetization.
The Stoner-Wohlfarth (SW) model captures the angle dependent switching characteristic of single-domain elliptical magnets [29]. The characteristic SW astroid is shown in Fig. 2c (blue line) and is described by the equation In the SW model, switching may occur when the left hand side of Eq. (4) is greater than one. The astroid obtained from micromagnetic simulations and the SW astroid ( Fig. 2c) are nearly identical, indicating that the SW model is a simple and valid description of switching in elliptical nanomagnets.
However, the SW model is only accurate for elliptical magnets, other magnet shapes typically have quite different switching characteristics. Fig. 2d shows the switching astroid for rectangular stadium-shaped magnets (red dots), which is the shape commonly used in most fabricated ASIs (Fig. 2b). Notice how the astroid is asymmetric: rectangular magnets switch more easily along the easy axis than the hard axis.
To capture the asymmetric switching characteristics of non-elliptical magnets, we have generalized the SW switching model to allow an asymmetry between easy and hard axes. Additionally, the model has been extended to support tuning of the curvature of the extrema. In the generalized model, the switching threshold is given by where b, c, β and γ are parameters which adjust the shape of the astroid: b and c define the height and width, respectively, while β and γ adjust the curvature of the astroid at the easy and hard axis, respectively. Introducing these new parameters allows for tuning of the switching astroid to fit with the shape of nanomagnets used in ASIs. With b = c = 1 and β = γ = 3, Eq. (5) reduces to Eq. (4), i.e., the classical Stoner-Wohlfarth astroid is obtained (valid for elliptical magnets).
By tuning the parameters of the generalized SW model, we can obtain the asymmetric switching astroid shown in Fig. 2d (blue line). The astroid is in good agreement with results obtained from micromagnetic simulations (red dots).
In flatspin, the generalized SW model is used as the switching criteria, i.e., a spin may flip if the left hand side of Eq. (5) is greater than one. Additionaly, the projection of h i onto m i must be in the opposite direction of m i :

G. Imperfections and disorder
Due to manufacturing imperfections there will always be a degree of variation in the shape and edge roughness of nanomagnets. This variation can be thought of as a disorder in the magnets' inherent properties. Rough edges and a slightly distorted geometry can affect the magnets' switching mechanisms, with defects pinning magnetization and altering the coercive field for each magnet.
In flatspin we model this variation as disorder in the coercive fields. The coercive field is defined individually for each magnet, and a distribution of values can be used to introduce variation. A user-defined parameter, k disorder , defines the distribution of coercive fields, i.e., h

H. Dynamics
flatspin employs deterministic single spin flip dynamics. At each simulation step, we calculate the total magnetic field, h i , acting on each spin. Next, we determine which spins may flip according to the switching criteria Eqs. (5) and (6). Finally we flip the spin where h i is the furthest outside its switching astroid, i.e., where the left hand side of Eq. (5) is greatest. Ties are broken in a deterministic, arbitrary manner, although with non-zero disorder such occurrences are rare. The above process is repeated until there are no more flippable spins.
This relaxation process is performed with constant external and thermal fields. To advance the simulation, the fields are updated and relaxation is performed again. Hence a simulation run consists of a sequence of field updates and relaxation processes.
The dynamical process makes three main assumptions: 1. The external field is quasi-static compared to the time scale of magnet switching.
3. The magnet experiencing the highest effective field compared to its switching threshold is the first to switch.
Assumption 1 means the model holds for low frequency external fields, i.e., where switching settles under a static field. The switching mechanics of nanomagnets are typically in the sub nanosecond range [9,13], and experimental setups often employ external magnetic fields which can be considered static at this time scale. At high applied field frequencies, more complex physical phenomena such as spin waves will have a non-negligible effect on switching dynamics, which are not modeled in flatspin.
Assumption 2 holds if the coercive fields h (i) c , and total field h i , of the magnets are sufficiently non-uniform, so that there will always be a single magnet which will flip first. It is assumed to be unlikely that two magnets will have the same h (i) c and h i simultaneously. However, in those rare cases where two magnets are equally close to switching, overlapping switching events may occur in a physical system. Assumption 3 relies on the fact that all changes in the magnetic fields are effectively continuous, and the change is unidirectional within a simulated time step, i.e., a quasi-static field. Since evaluation happens in discrete time, there will be cases where several magnets are above their corresponding switching thresholds simultaneously. In those cases, the magnet furthest above its switching threshold will have been the first to have crossed the threshold under a quasi-static field. Furthermore, if the angle of the external field is constant, the switching order is invariant to the time resolution of the external field.

I. Geometries
The particular spatial arrangement of the magnets is referred to as the geometry. A wide range of ASI geometries have been proposed in the literature. Fig. 3 depicts the geometries included in flatspin, which are the most commonly used ASI geometries: square [31], kagome [22,28], pinwheel [16,17] and ising [24].
Geometries are often decomposed into two or more "sublattices", where the magnets within one sublattice are all aligned (have the same rotation). In Fig. 3, the sublattice a magnet belongs to is indicated by its color. As can be seen, both square and pinwheel ASIs have two perpendicular sublattices, whereas kagome has three sublattices.
flatspin can be used to model any two-dimensional ASI comprised of identical elements. New geometries can easily be added by extending the model with a new set of positions r i and rotations θ i .

III. SIMULATION FRAMEWORK
In addition to a magnetic model, flatspin provides a flexible framework for running simulations, storing re-sults, and performing analysis. Fig. 4 illustrates the overall architecture of flatspin. The ASI model has been described in detail in section II. Conceptually, the ASI model describes the physical system under study. The rest of the components are tools used by the experimenter to interact with the ASI and observe the results. In this section we briefly describe each of these components.
The input encoder translates a set of input values to a series of external fields. Encoders provide a flexible way to define field protocols, and have been designed with neuromorphic computing in mind. A range of encoders are included, e.g., sinusodial, sawtooth and rotational fields.
The responsibility of the runner component is to perturb the ASI model according to the field protocol, and save the results. The model, which is fully parametric, receives parameters from the runner, enabling automated parameter sweeps. In addition, there is support for distributed running of simulations on a compute cluster.
Results are stored in a well-defined dataset format which makes the analysis of a large numbers of simulations straightforward. A suite of analysis tools are included, e.g., for plotting results, visualizing ensemble dynamics and analysis of vertex populations. flatspin is open-source software and released under a GNU GPL license. For more information see the website [6] and User Manual [5].

IV. VALIDATION OF FLATSPIN
To evaluate the ASI model, flatspin simulations were compared to established experimental results from the literature, as well as micromagnetic simulations. In particular we investigate phenomena such as Dirac strings in kagome ASI, the size of crystallite domains in square ASI, and superferromagnetism in pinwheel ASI. Finally, we compare the switching order from flatspin simulations with that of micromagnetic simulations, and investigate the effect of varying lattice spacings.

A. Dirac strings in kagome ASI
To assess the ability of flatspin to reproduce fine-scale patterns, we consider the emergence of Dirac strings in a kagome ASI (Fig. 3c). Applying a reversal field to a polarized kagome ASI results in the formation of monopoleantimonopole pairs [18]. These pairs are joined by a "string" of nanomagnets which have flipped due to the reversal field. As the strength of the reversal field increases, the strings elongate until they fill the array.
We closely follow the methodology set out in an experimental study of Dirac stings in kagome ASI [18], in which a room temperature kagome ASI undergoes magnetiza-tion reversal. We start with an array of 2638 magnets (29 × 29 hexagons) polarized to the left and apply a reversal field H to the right with a slight, downward offset of 3.6 • . This offset breaks the symmetry, such that one of the sublattices is now least aligned with the field, resulting in an increased coercive field on this "unfavored" sublattice.
The time evolution snapshots of Fig. 5 demonstrate a strong, qualitative similarity to the results of Mengotti et al. [18]. We see Dirac strings developing with a preference to lie along the two sublattices most aligned with the field angle. Furthermore, in the final image, we see the vast majority of unflipped magnets (excluding the edges) are on the unfavored sublattice, in accordance with both experimental and simulated results from the literature. Also in Fig. 5, we see the hysteresis of the simulated ensemble (solid line) is similar to that of Mengotti et al. [18] (dashed line) in some sections, but differs near the extrema. The hysteresis can be understood in two phases. The first phase, at roughly M/M S ∈ [−0.6, 0.6], is dominated by the lengthening of the Dirac strings, with almost no activity occurring on the unfavored sublattice. At M/M S < −0.6 and M/M S > 0.6, the ensemble enters a second phase in which the Dirac strings have fully covered the array, and change in magnetization is dominated by switching on the unfavored sublattice. Clearly we see good agreement, within phase one, between our simulated hysteresis and the experimental results. Furthermore, there is a clear phase transition (characterized by a sharp decrease in gradient) in our hysteresis very close to the transition in the experimental hysteresis. Notably however, although the phase transitions occur at a similar time, the change in gradient is less pronounced in our simulated hysteresis. This disparity indicates that, in the second phase, the magnets on the unfavored sublattice flip more easily in our simulation than in the experimental data.
The accuracy of the point dipole approximation is known to suffer when considering kagome ASI. Specifically, it has been shown to underestimate the coupling coefficient of the nearest neighbors by approximately a factor of 5 [23]. Despite this, we observe flatspin accurately reproduces snapshots of the time evolution behavior, while also capturing salient features of the ensemble hysteresis curve.

B. Domain size in square ASI
In order to demonstrate simulation of large-scale behavior, we have reproduced the emergence of large do-  I I I I I I I I I   0.85Hc   II II II II II II II II II   0.92Hc   III III III III III III III III  mains of magnetic order in square ASI, similar to experimental results of Zhang et al. [32]. One of the main advantages of flatspin over typical alternatives is the scalability and high throughput of large systems with many magnets. Some emergent ASI phenomena require large systems in order to be fully quantified and studied with high fidelity, such as the domain size of magnetic charge crystallites. For ASIs with strongly coupled magnets, typical domain sizes can become too large for direct experimental observation. Thus, an accurate estimate of the domain size for ASIs with a small lattice constant is, in part, limited by the number of directly observable magnets.
For a given range of lattice constants covering both strongly coupled ASIs and weakly coupled ASIs, a corresponding range of large to small magnetic order coherence lengths is expected. In this study, we consider square ASI (closed edges, Fig. 3a) with different lattice constants, a, ranging from 320 nm to 880 nm. 50 × 50 square ASIs were annealed in flatspin using a thermal protocol of exponentially decreasing σ th . A switching astroid for 220 nm × 80 nm × 25 nm was obtained through micromagnetic simulations, described by generalized astroid parameters b = 0.4, c = 1.0, β = 3.0, and γ = 3.0. Additionally, h k = 0.186, k disorder = 0.05, and a neighbor distance of 10 magnets was used. The thermal protocol was chosen such that the total dipole interaction energy was not significantly reduced by increasing simulation steps.
In the annealed state, the spin-spin correlation as a function of their lateral separation was calculated across the ensembles. Analysis of the average correlation of annealed states provides insight about the typical coherence length of magnetic order, i.e., magnetic charge crystallite size, or domain size. Here, the correlation of two spins is defined as +1 (-1) if their dipole interaction is minimized (maximized). Averaging correlation across distinct types of spin pairs, in the annealed ASI, gives a measure of how coherent the ASI is at that particular neighbor separation. How quickly the average correlation decreases as a function of separation can be used to estimate the characteristic domain size. In particular, it can be argued that the separation where the correlation falls below 1/e is the characteristic domain radius [20,32]. Typical domain structures and correlation results can be seen in Fig. 6. The domains shown in Fig. 6a and the correlation curves in Fig. 6b are, for the most part, in good agreement with experimental results [32]. A qualitative comparison of the domain sizes and structures in Fig. 6 shows that the domains tend to be larger, with smoother domain boundaries, for smaller a. The analysis of coherence as a function of separation also shows identical trends and similar values, where an increase in a leads to low correlation, even between nearest neighbors.
For cases where a < 400 nm, domains do not significantly increase in size when the lattice constant is further reduced. This discrepancy with the experimental results is not completely unexpected: the point-dipole approximation is known to underestimate nearest-neighbor interaction for magnets placed close together [23]. In ad- dition, a stronger interaction between spins would cause each spin flip to contribute a greater change in the total dipole energy. This makes a gradual descent towards the ground state by random spin flips (the thermal fluctuations as modeled by flatspin) harder to achieve. These issues may be addressed by increasing the coupling parameter α for nearest neighbor spins, and by a longer and slower annealing protocol. A longer and slower annealing protocol will inevitably come at the cost of longer computation times. In a future version of flatspin, it might be beneficial to include other thermal effects such as a declining saturation magnetization of the constituent material.
These results show that flatspin provides sufficient flexibility, fidelity and performance required to reproduce experimentally observed large-scale emergent behavior in ASIs.

C. Superferromagnetism in pinwheel ASI
In this section, we use flatspin to reproduce the dynamic behavior of pinwheel ASI, which had yet to be demonstrated with a dipole model [16]. We find that our switching criteria plays a key role in replicating magnetization details during the field-driven array reversal.
Pinwheel ASI is obtained by rotating each island in square ASI some angle about its center. A rotation of 45 degrees results in a transition from antiferromagnetic to ferromagnetic order [17]. The dynamics of pinwheel ASI in many ways resemble continuous ferromagnetic thin films, with mesoscopic domain growth originating from nucleation sites, followed by coherent domain propagation and complete magnetization reversal [16].
Here we demonstrate that flatspin is able to replicate the experimental reversal processes presented in Li et al. [16], where pinwheel "diamond" ASI ( Fig. 3d) is subject to an external field at different angles. A key result is that the angle θ of the external field controls the nature of the reversal process. When θ is small (equally aligned to both sublattices), reversal happens in a single avalanche, whereas when θ is large (more aligned to one sublattice), reversal happens in a two-step process where one sublattice switches completely before the other. Previous attempts at capturing this behavior in a dipole model have proven unfruitful [16].
To replicate this process in flatspin, an asymmetric switching astroid is required, i.e., the threshold along the parallel component is reduced by setting b < 1 in Eq. (5). From micromagnetic simulations of a single 470 × 170 × 10 nm magnet we obtain the following characteristic switching parameters: b = 0.28, c = 1.0, β = 4.8 and γ = 3.0. Other simulation parameters include α ≈ 0.00033, h k = 0.098, k disorder = 0.05 and a neighbor distance of 10. Figs. 7a to 7d show hysteresis loops and array snapshots when the field is aligned with the array (θ = 0 • and θ = −6 • ). As can be seen, the results from flatspin (Figs. 7b and 7d) are qualitatively very similar to experimental results (Figs. 7a and 7c). In all cases, the ASI undergoes reversal in a single avalanche. Reversal begins at a small number of nucleation points close to the edge, followed by domain growth and domain wall movement perpendicular to the direction of the field. The simulated system appears to have an anisotropy axis of 0 • as opposed to −6 • observed experimentally. Hence Fig. 7b is most similar to Fig. 7c and Fig. 7d is most similar to Fig. 7a. It should be noted that the non-zero anisotropy axis found experimentally has not yet been explained.
Figs. 7e and 7f show the hysteresis loops and array snapshots when the field is misaligned with the array (θ = 30 • ). Again, flatspin simulations (Fig. 7f) replicate key features observed experimentally (Fig. 7e). Reversal now happens in two steps: the sublattice whose magnets have their easy axis most aligned with the field will switch first, followed later by the other sublattice. This two-step reversal process results in an emergent rotation of the collective magnetization. The magnetization is constrained to follow the orientation of the magnets, resulting in reversal via stripe patterns at 45 • .
Li et al. [16] report they were unable to replicate the magnetization details using a point-dipole Monte Carlo model. One crucial difference between flatspin and their dipole model is the switching criteria. They use the simpler criteria h i · m i < h (i) k , which considers only the parallel field component and will be largely inaccurate for fields which are not aligned with the magnet's easy axis. Indeed, using this simpler switching criteria in flatspin results in a very different reversal process and magnetization details (not shown).

D. Comparison to micromagnetic single-spin switching order
Micromagnetic simulations in MuMax3 [30] have been shown to agree with experiment due to high simulation fidelity. It is, therefore, of interest to study how well flatspin agrees with MuMax3 at the level of detail expressed in flatspin.
Here we evaluate the switching strategy outlined in Section II H, by comparing the switching orders obtained in flatspin and MuMax3, of a square ASI as it undergoes reversal by an external field. As a similarity measure, Spearman's rank correlation coefficient ρ [8] is used, where a value of 1 indicates perfect correlation and 0 indicates no correlation between switching orders.
In the weakly coupled regime, the switching order is dictated by the coercive fields of each individual magnet. In flatspin, the coercive field can be set directly by modifying h (i) k . In MuMax3, we control the coercive field implicitly, by varying the first-order, uniaxial, magnetocrystalline anisotropy, K The system we considered was a 4 × 4 square (closed) ASI, each magnet measuring 220 nm × 80 nm × 25 nm. flatspin was run with parameters b = 0.38, c = 1, β = 1.5, and γ = 3.2. In both simulators, we applied a gradually increasing reversal field at θ = 44 • .
At a certain point the dipolar interactions begin contributing to the switching order. To verify that flatspin still captures switching dynamics, we perform a comparision of the switching orders for all pairs of lattice spacings in both simulators. Fig. 8a shows the correlations for each pair of lattice spacings as an average over 32 different square ASIs. We observe a clear linear relationship between the two simulators, with higher lattice spacings exhibiting higher correlation. The non-zero y-intercept in the heatmap indicates that, as suspected, the coupling strength is slightly underestimated by the dipole approximation employed in flatspin, in particular for lower lattice spacings. For example, flatspin with 300 nm lattice spacing is most similar to MuMax3 with 380 nm.
The red line in Fig. 8b traces the ridge in the heatmap, i.e., the highest ρ, for each flatspin lattice spacing. As can be seen, a near-perfect agreement between the simulators is found in the weakly coupled regime (high lattice spacing). As lattice spacings decrease, the effect of dipole interactions become apparent. Below 450 nm the correlation drops. Since flatspin does not account for the micromagnetic state, complete correlation is not expected.
The particular selection of h (i) k values introduces an inherent bias in the switching order. One might expect that this bias causes the dipole interactions to have a negligible impact on the switching order, leading to an inflated correlation between flatspin and MuMax3, regardless of lattice spacing.
The violet and blue lines, (rightmost column and top row from Fig. 8a), show the deviation from the inherent switching order in MuMax3 and flatspin, respectively. Their rapid decline confirms that the switching order is not dominated by the inherent bias for highly coupled systems. The red line clearly shows a stronger agreement between MuMax3 and flatspin (higher than the blue and violet lines), even at smaller lattice spacings.

V. PERFORMANCE
Although the total simulation time will depend on many factors, it is of interest to measure how simulation time scales with the number of spins. As the number of spins are increased, simulation time will be largely dominated by the calculation of the effective field, h i , acting on each of the N spins in the lattice. Computing time for h  edge magnets is negligible (in the common ASI geometries). Computing h i for all spins will take O(N ) time, i.e., computation time grows no faster than linear in N . Fig. 9 shows the throughput (number of field calculations per second) as a function of the number of spins. Here a field calculation is defined as the computation of h i for a single spin i, hence for N spins there will be N such field calculations. The geometry used was square ASI (open edges) using a standard 8 spin neighborhood for calculating h (i) dip . The throughput was averaged over 100 simulations of each size. The test was performed on an NVIDIA Tesla V100 GPU with 32GB of RAM.
At around 200 000 spins, the throughput saturates at 10 8 field calculations per second. On our test setup, computing h i for one million spins takes approximately 10 ms. Above 200 000 spins we are able to fully utilize the GPU resources.
To simulate the reversal of an ASI by a gradually increasing external field, at least one field calculation per spin flip is required, i.e., at least N field calculations. If the external field gradually changes with a resolution of K values, the worst case will be when all spins flip during a single field value. In this case the number of field calculations required will be N + K − 1 since there will be K − 1 field calculations which results in no spin flips.
The total simulation time depends largely on the particular experimental setup, parameters and other system characteristics. Time will be spent on things other than field calculations, e.g., organizing and writing results to storage. Hence the total simulation time will be longer than predicted by field calculations alone. As an example, the simulations from section IV C of 25×25 pinwheel ASI with 1250 magnets took approximately 6 seconds with K = 2500, for one reversal. Fig. 10 shows a snapshot from flatspin simulations of 708 × 708 pinwheel ASI as it undergoes reversal by an external field. The ability to simulate such large systems allows an experimenter to explore phenomena at much larger scales than can be directly observed experimentally. With more than one million magnets, the simulation of array reversal took several days to complete. A video of the full reversal is available as Supplemental Material [Note1].

VI. CONCLUSION
flatspin is a highly effective simulator for ASI ensemble dynamics. At its heart lies a robust magnetic model based on dipole-dipole interactions with a switching crite- The flatspin ASI model has been verified against micromagnetic simulations and experimental results from the literature. On a detailed level, we found good agreement between micromagnetic simulations and flatspin in terms of magnet switching order. Emergent fine-scale patterns in kagome ASI were replicated successfully, where the formation of Dirac strings matched experimental results. Large-scale domain sizes in square ASI were reproduced, and good agreement was found between flatspin and experimental results. Finally, using flatspin, the experimental magnetization reversal of pinwheel ASI was reproduced for the first time in a dipole model.
Through GPU acceleration, flatspin scales to large ASI systems with millions of magnets. High speed, parallel computation allows for a large number of ASI simulations to be executed, enabling quick exploration of parameters and novel geometries. The flexibility and performance offered by flatspin opens for unprecedented possibilities in ASI research.