Dynamics and hysteresis in square lattice artificial spin-ice

Dynamical effects under geometrical frustration are considered in a model for artificial spin ice on a square lattice in two dimensions. Each island of the spin ice has a three-component Heisenberg-like dipole moment subject to shape anisotropies that influence its direction. The model has real dynamics, including rotation of the magnetic degrees of freedom, going beyond the Ising-type models of spin ice. The dynamics is studied using a Langevin equation solved via a second order Heun algorithm. Thermodynamic properties such as the specific heat are presented for different couplings. A peak in specific heat is related to a type of melting-like phase transition present in the model. Hysteresis in an applied magnetic field is calculated for model parameters where the system is able to reach thermodynamic equilibrium.

Artificial spin ices are systems in two dimensions that mimic the usual three-dimensional spin ice materials that exhibit geometrical frustration effects: not all the pairwise spin interactions can be satisfied simultaneously [1][2][3][4][5] . Indeed, the name spin ice comes from the fact that these systems obey the ice rule: for a square lattice, in each vertex four spins meet but two must point inward while the other two must point outward. In general, the artificial compounds are built from magnetic nanoislands (generally made of permalloy) which can be organized in several different geometries in such a way that the frustration is manifested 6-12 . Here, the focus is on the artificial square lattice spin ice array, first fashioned and studied by Wang et al. 6 in 2006. This artificial square ice consists of magnetic nanoislands (with a shape that looks like a "cigar") arranged as shown, for instance, in Ref. 6. Each nanoisland contains a net magnetic moment that tends to point out along the longest axis. When the interactions between neighboring islands are increased, the system of Ref. 6 was increasingly filled with vertices that obeyed the ice rule (twoin/two-out, which is non-degenerate in two dimensions). Despite this, the predicted ground state of square ice was not observed experimentally until the work by Morgan et al. 13 in 2011. Using more advanced techniques, those authors could observe large regions of their samples which were able to adopt the square ice ground state. With this progress, they could also observe the predicted excitations above the ground state, which resemble magnetic monopoles connected by energetic strings [13][14][15][16] (similar to Nambu monopoles 17 ). These elementary excitations are different from those of natural three-dimensional spin ices, which are magnetic monopoles connected by observable but non-energetic strings 4,5 . Therefore, these artificial compounds have attracted large interest in recent years.
The manipulation and control of the excitations could lead to several technological applications. However, there is a primary problem with these artificial systems, particularly with the square ice: the lack of thermalization. Thermal activation of the island's magnetic moment ("spin") configurations is very weak or nonexistent. Lack of thermalization is an important topic for the experimental artificial spin ices and it can be partially alleviated by applying varying external magnetic fields 7,18 . Moreover, reductions in island volume and magnetic moment through state-of-the-art nanofabrication can bring energy scales closer to room temperature, leading to thermally driven slow dynamics. Alternatively, the use of materials with an ordering temperature near room temperature seems to be another important possibility. Really, by using such a material, a recent experimental work on a square lattice in an external magnetic field confirms a dynamical "pre-melting" of the artificial spin ice structure at a temperature well below the intrinsic ordering temperature of the island material, creating a spin ice array that has real thermal dynamics of the artificial spins over an extended temperature range 19 . Another expectation for better understanding of these compounds is to study similar systems such as colloidal systems, which have an advantage over the usual magnetic arrays because thermal activation of the effective spin degrees of freedom is possible 11 . So, a more detailed analysis of the effects of thermal fluctuations and the spin dynamics in a twodimensional spin ice material should be of great interest for a better understanding of these interesting frustrated systems.
Particularly, thermal effects in artificial square ice were studied recently by some of us 20 . The focus was to examine the roles of elementary excitations in the thermodynamic properties of these systems. We have found that the specific heat and average separation between monopoles with opposite charges exhibit a sharp peak and a local maximum, respectively, at the same temperature 20 , T p ≈ 7.2D/k B , where D is the strength of the dipolar interaction and k B is the Boltzmann constant. These features increase logarithmically with the lattice size and, therefore, the possibility of string breaking and free monopoles were discussed. However, all calculations were done by assuming an Ising-like behavior for the net magnetic moment of the nanoislands. The Ising behavior of the islands seems to be the most realistic situation for the typical artificial magnetic ices made of permalloy and therefore, these systems could not display real dynamics; it is one of the causes concerning the problem of lack of thermalization. In this article, we study possibilities beyond the Ising behavior for the nanoislands and consequently, our attention shifts to magnetic ices with real dynamics and the extra features that such dynamics may produce. For this study, the internal structure and shape of the magnetic nanoislands have to be taken into account in details. The theoretical study of the net magnetic moment (with more degrees of freedom) of individual magnetic nanoislands with different sizes and shapes is the initial point. Such an investigation could also predict better circumstances for novel ideas and dictate protocols for the realization of artificial magnetic spin ices with real thermal activation and spin dynamics.
A detailed study concerning non-Ising spins for individual islands can be found in Ref. 21 . By applying the techniques and ideas of our previous paper 21 , we now consider a more intricate situation concerning artificial square spin ices. Here we study the dynamics of the array by assuming that a nanoisland's spin is free to point in any possible direction, but with strong uniaxial anisotropy energies that favor preferred directions. Then, differently from all previous articles published until now on this topic, which consider basically only the dipolar interactions among the islands, here, additional terms concerning these anisotropies must be included in the Hamiltonian. We would like to answer questions like these: 1) Is there any possibility of building artificial magnetic ices (for instance, by using permalloy islands) exhibiting real dynamics, only by decreasing the island sizes? 2) If there is this chance, what are better parameters for the islands and for the arrays with spin dynamics? 3) If there is not such a chance, what should the parameters be for other potential magnetic materials and their ice arrays, so they exhibit real dynamics?
The article is organized as follows: In section II, the model is explained in detail. We also define two order parameters, which are important to identify if the ground state of the models analyzed can be accessed at low temperature. Five models (denoted by A, B, C, D and E) with different lattice and island parameters are proposed to see the possibility of a spin ice dynamics. The numerical calculations are presented in Section III; The dynamics is investigated by using a Langevin approach for the spins. In Section IV , some thermal equilibrium properties of the models A, B, C, D and E are calculated. In Section V , we present the hysteresis calculations and, finally, the conclusions and some discussions are given in Section V I.

II. THE MODEL SYSTEM
The open square ice system with N c = L 1 × L 2 unit cells can be set up as follows. One can define the sites of a square lattice in the usual way: the k th lattice site (a monopole charge center) is at a point r k = (x k , y k ), where x k = m k a and y k = n k a are integer multiples of the chosen lattice constant a. The points are chosen to fit inside the desired L 1 × L 2 area. For each unit cell of size a × a there are two nano-islands that act as a two-atom basis, with the locations, where the "1" and "2" refer to the sublattices. The nanoisland on the 1 st sublattice has its long axis in the x direction; the other nano-island, on the 2 nd sublattice, has its long axis in the y direction. There are N = 2N c islands in the whole system. A 3D vector magnetic moment µ i , i = 1, 2, 3...N is associated with each island, whose defined center position is some r i . Each r i is selected from the set of r kσ , where σ = 1, 2 denotes the sublattice. We use indeces k, l for locations of the unit cells, and indeces i, j for locations of individual islands or their dipoles.
The dipoles are assumed to have fixed magnitude µ, while their direction is represented by a unit vectorμ i . The magnetic moments interact via long-range dipole forces, and are also affected by two forms of local shape anisotropy. First, there is a uniaxial anisotropy that impedes free rotation in the xy-plane, associated with some energy constant K 1 , and oriented along x for the first sublattice and along y for the second. Dependent on its sublattice, each moment has an axisû i (equal tox orŷ) for this anisotropy. Second, because the nano-islands are thin in the z-direction, the z-direction is a hard axis, and there is a hard-axis anisotropy whose energy scale is determined by a constant K 3 , the same for all the islands. The Hamiltonian is then Here µ 0 is the magnetic permeability of space, andr ij is the unit vector pointing from the position of µ j towards the position of µ i . The first sum is the dipole-dipole interactions, the second sum contains the anisotropy energies and an applied external magnetic induction B ext = µ 0 H ext . A constant is included in the K 1 anisotropy energy so that that energy is zero when a dipole points along its local anisotropy axisû i . It is interesting to note that if a dipole moves in the xy plane, it only pays the cost of the K 1 anisotropy term, but motion up out of the xy plane (say, in the xz plane) involves an energy proportional to the sum of both anisotropies, The motion out of the xy plane is also impeded by the dipolar interactions. With the dipole pair distances scaled by the lattice constant, the effective strength of nearest neighbor dipolar interactions is determined by the dipole energy factor, Depending on the island geometry, which is discussed further below, the anisotropy constants K 1 and K 3 would typically be of similar order of magnitude. Thus, there are three important energy scales: dipolar energy, anisotropy energy, and the thermal energy k B T . The anisotropy constants are proportional to the volume of the islands, as is µ = M s V , where M s is the saturation magnetization of the magnetic material. But then, this dipolar constant D increases as the squared island volume. Thus, it is possible to adjust the island size (and spacing a) to get these energy scales as desired in relation to each other. Typically, the interesting case must have the thermal energy less than both the effective dipolar energy (per site) and anisotropy energy. But note, the effective dipolar energy can be quite a lot larger than that indicated by D, which only measures the energy in a nearest neighbor pair. When the dipolar interactions are summed, the net dipolar energy associated with any island could be much larger than D.
A. Spin-ice ground state and order parameters For the square lattice spin-ice, the ground state is twofold degenerate, and involves alternating dipoles on each of the two sublattices. The ground state fully satisfies the two-in/two-out rule in each monopole charge cell (junction of four islands at the site r k of each unit cell). The unit cell positions are expressed r k = (m k , n k )a, where a is the lattice constant and m k and n k are integers. Then one of the ground states can be constructed by setting the dipole directions as: This formula is arranged so that at a chosen unit cell at position r k , the dipole on one sublattice points inward and the dipole on the other sublattice points outward, thereby globally enforcing the two-in/two-out rule. By reversing the sign on all the dipoles, the other ground state is obtained. With the ground state determined, we can construct a measure of the proximity (in phase space) of any arbitrary state to one of the ground states. This order parameter Z is simply the overlap with this ground state: The index σ labels the sublattice. If the system happens to be found in the ground state defined in Eq. (4), then Z = 1; if the system is in the inverted ground state, then Z = −1. Thus it is possible to show that the range of Z is from -1 to +1. This order parameter is useful for indicating the degree of thermodynamic excitation in the system, by the deviation of |Z| from unity. Further, its sign then gives an indication of processes which involve the transformation from one ground state to the other. Indeed, considered even as a local variable (calculating only near a single charge cell), we can track when the system has different regions close to either of the ground states, possibly with regions separated by domain walls. Fig. 1 shows an example, where the system has a region near one of the ground states, with Z ≈ +1, separated by a domain wall from another region that is near the other ground state, with Z ≈ −1. The net averaged value of Z for the entire system, however, acquires an intermediate value, Z ≈ 0.521, indicating considerable separation from a uniform ground state. The other obvious order parameter to be measured is the areal density of monopole charges, ρ m . We make a simple discrete definition, to connect to Ising spin ice models, and, a more generalized continuous definition that accounts for the greater freedom of the continuous dipoles in the model described here. The discrete definition of a monompole charge involves counting the net number of dipoles that point outward at a chosen charge site r k , and dividing that result by two. There are four dipoles µ i k , i k = 1, 2, 3, 4, surrounding any charge cell center r k . Then the possible monopole charge values are q k = 0, ±1, ±2, although the last values may typically be of low probability. For the discrete charge definition, whether a dipole points outward or inward is determined with a Heaviside step function H(x): The unit vectorsv i k , i k = 1, 2, 3, 4, point outward from charge site r k to each of the four nearest islands. Because this discrete definition can show sudden change when a dipole rotates 90 • from the radially outward direction, we also considered a continuous definition. In the continuous definition, the step function is removed, and only a scalar product is needed, In contrast to the discrete definition, this charge definition varies continuously from q k = −2 to q k = +2. One can also note that for either the discrete or the continuous definition, total monopole charge is conserved. A positive contribution produced by some dipole at one charge cell is accompanied by an equal negative contribution at a neighboring charge cell (each dipole contributes to two charge cells). Then, the total algebraic monopole charge in the system takes the conserved value, zero. In order to get a measure of the monopole charges present, regardless of their sign, we define a density for the system as a whole, by applying absolute value. Thus, the "monopole density" measured in the simulations here is defined as This is normalized by the number of charge cells. By using absolute value, the definition does not allow the cancellation of charges of opposite signs. Note that in either of the ground states, there are no charges at any sites, and ρ m = 0. Charges appear as the system moves away from the ground state. Thus, this is another measure of excitation in the system. (This is true for the spin ice on the square lattice, but not on the Kagomé lattice, whose ground state contains charges, due to there being three dipoles for each charge cell.)

B. The undamped dynamics
The zero-temperature, undamped dynamics of each magnetic dipole is determined by a torque equation, where B i is the local magnetic induction acting on the i th dipole, and γ is the electronic gyromagnetic ratio. The local magnetic induction is derived from the Hamiltonian by assuming an energy − µ i · B i for each dipole, i.e., It will be convenient to choose some standard units for the time, the applied field, and so on, to simplify and scale the numerical calculations. The dipole terms are simplified by selection of the lattice constant a as the unit of length. A natural unit to measure field H ext is the saturation magnetization M s from which the particles are made. For example, for Permalloy, with M s = 860 kA/m, this unit as a magnetic induction is close to one tesla: µ 0 M s ≈ 1.08 T. Using this quantity to scale the magnetic field and hence the magnetic induction defines their dimensionless field, When h ext approaches 1.0 the applied field should have a strong tendency to saturate the magnetization of the system (if the dipolar interactions do not impede that). This then indicates how to scale the dipole and anisotropy fields, i.e., by writing the dimensionless local magnetic fields from (10), This suggests the introduction of dimensionless coupling constants that indicate the relative strength of each contribution, These definitions involve the different energy scales divided by an energy unit, that depends on the size of the magnetic islands. The dimensionless dipole parameter d can be seen to be proportional to the volume fraction of the system occupied by magnetic islands, since µ = M s V for each island. Obviously a higher packing of magnetic material into the lattice leads to stronger dipolar effects, and d indicates their effective strength.
The dynamic equation can be scaled in the same way, so that the dimenionless field appears on the RHS. Then the dynamics for the unit vector dipoles is described using a rescaled time τ , With the above scaling of the fields, the unit of time is (γµ 0 M s ) −1 . For the case of Permalloy and using the gyromagnetic ratio as γ = e/m e ≈ 1.76 × 10 11 T −1 s −1 , this unit is (γµ 0 M s ) −1 ≈ 5.26 ps.

C. Island geometry and energetics
The shape anisotropy constants K 1 and K 3 can be estimated based on the magnetic properties for Permalloy (or other material) and micromagnetics simulations for the choice of island geometries and island volume V . In this work we consider thin elliptical islands. Here L x denotes the major-diameter of the ellipse and L y is the minor-diameter, while L z is the height of the island or its thickness. The semi-major axis is A = L x /2, the semiminor axis is B = L y /2. It is well-known that the an elliptically shaped magnetic particle will have a nonuniform anisotropy 22 within the plane of the island. In Ref. 21, the anisotropy constants (as energies per unit volume) were estimated based on a calculational approach for thin elliptical islands, for a range of thicknesses L z ≪ L x , characterized by an aspect ratio g 3 = L x /L z , and various lateral aspects ratio g 1 = L x /L y . Here we consider some different sizes and shapes for the islands and discuss the advantages and disadvantages, as far as the expectations for their dynamics and and the relative importance of the different energy scales when placed in a square spin-ice array.
Model A. In Wang et al. 6 , experiments on square spin-ice were carried out for (quasi-rectangular) particles with dimensions 220 nm × 80 nm × 25 nm, where the last number is the vertical thickness. In those experiments, the particle sizes were kept fixed, but different lattice parameters a from 320 nm to 880 nm were used. In this first model to consider, we use these numbers to describe elliptical particles: L x = 220 nm, L y = 80 nm, L z = 25 nm. Then the particle volume is V = πABL z = 3.46 × 10 5 nm 3 , and using a saturation magnetization M s = 860 kA/m for Py, the magnetic dipole moment per particle is µ = 2.97 × 10 −16 A·m 2 , the equivalent of about 3.2×10 7 Bohr magneton. For its aspect ratio parameters g 1 = 2.75 and g 3 = 8.8, the anisotropy energy densities can be found by interpolation of the simulation results in Ref. 21 , as K 1 /V = 0.0064 A ex /nm 2 and K 3 /V = 0.0143 A ex /nm 2 , where A ex ≈ 13 pJ/m is the exchange stiffness for Py. The easy-axis anisotropy then works out to K 1 = 2.9×10 −17 joules, while the hard-axis anisotropy is estimated as K 3 = 6.4 × 10 −17 J. These are considerably larger than room-temperature (300 K) thermal energy k B T ≈ 4.1 × 10 −21 J, as needed for stable magnetic moments. The energy unit is ε = µ 0 µM s = 3.21 × 10 −16 J. Then the dimensionless anisotropies are The nearest neighbor dipolar energy scale might be estimated first at lattice constant a = 880 nm, for which it is D = 1.29 × 10 −20 J, about 2000 times smaller than K 1 . If instead the lattice constant a = 320 nm is used, this will scale up by a factor of (880/320) 3 , leading to D = 2.68 × 10 −19 J, or still 100 times smaller than K 1 . The dimensionless dipolar coupling for a = 320 nm is Obviously, values of k 1 , k 3 , and d similar to these are needed to get a spin-ice system, however, dynamics simulations could be difficult with these parameters because the anisotropy is so dominant and Ising-like. Over the time scales that can be accessed in numerical simulations, one would not expect to see much dynamical flipping of the island dipoles, except in the presence of a strong applied magnetic field. Thus it may be interesting instead to consider some other particle sizes where the dynamics can be expected to be more active.
A thinner or smaller island will result in a smaller magnetic dipole moment µ, which leads linearly to weaker anisotropy, but quadratically to weaker dipolar energy. Both energy scales become closer to the thermal energy. Thus we can try to change the particle size in such a way so that room temperature thermal energy is closer to K 1 and perhaps even larger than D.
Model B. Consider smaller elliptical particles, with L x = 100 nm, L y = 12.5 nm, L z = 5.0 nm, on a lattice with constant a = 200 nm. The aspect ratio parameters are g 1 = 8.0, g 3 = 20. In this case, the smaller particle volume V = 4910 nm 3 leads to a smaller magnetic dipole moment µ = 4.22 × 10 −18 A·m 2 , and the anisotropy energy constants are found to be K 1 = 6, 66 × 10 −19 J and K 3 = 7.48 × 10 −19 J. The energy unit is ε = µ 0 µM s = 4.56 × 10 −18 J. This leads to scaled energy parameters k 1 = K 1 /ε = 0.146, k 3 = K 3 /ε = 0.164, and for lattice constant a = 200 nm, one has d = D/ε = 4.87 × 10 −5 . The effective dipolar coupling is very small in this case. The scaled room temperature thermal energy (300 K) is now T = k B T /ε = 0.908 × 10 −3 , stronger than the dipolar coupling d. Thus this model might exhibit observable changes in its properties for easily accessible tempera-tures.
Model C. One can also consider even smaller particles, with L x = 40 nm, L y = 8.0 nm, L z = 4.0 nm, to try and get weaker anisotropy energy scales (still using the Py parameters). The particle volume is now only V = 1005 nm 3 , and the dipole moment is µ = 8.64 × 10 −19 A·m 2 . At aspect ratio parameters g 1 = 5.0, g 3 = 10, the anisotropy energies are found to be somewhat smaller: K 1 = 1.38 × 10 −19 J, and K 3 = 1.12×10 −19 J. The energy unit, though, is now also smaller: ε = 9.34×10 −19 J, leading to dimensionless couplings k 1 = 0.148 and k 3 = 0.120, only slightly different from those in Model B. However, the smaller energy unit means that room temperature effects may be more accessible. The scaled thermal energy at 300 K is increased: It is clear from the above examples, and also from the numerical results presented below, that the strong Isinglike anisotropy for real spin-ice particles dominates over thermal energy at (and below) room temperature. That being the case, we find it interesting also to study some models with fictitious parameters, which might be possible to achieve in other materials with different values of Model D. Rather than assuming a particular particle size and using Py parameters, suppose some particles are arranged so that D = K 1 = K 3 = 1 10 ε. Obviously nature may not easily produce such a system with all equal energy scales, but it may be possible by appropriate materials engineering. We use a fraction of ε, which is required by the definition of d, see Eq. (13) (the volume fraction of dipoles on the lattice cannot be more than unity). Then the scaled energy parameters are all equal: d = k 1 = k 3 = 0.1 . A physical value of ε is needed, based on values of µ and M s for some real particles, to locate room temperature on the temperature scale.
Model E. This is similar to Model D, but with a stronger value for the hard-axis anisotropy: D = K 1 = 1 10 ε, and K 3 = 1 2 ε. This would account for the typically stronger anisotropy for out-of-plane dipolar motion, due to the assumption of thin nano-islands. Only some limited simulations were done for this case, because the results are found to be very similar to those for Model D.

III. NUMERICAL CALCULATION: LANGEVIN DYNAMICS
The dynamics is investigated here using a Langevin approach for the "spins"μ i . This requires including the combination of a damping term and a rapidly fluctuating stochastic torque in the dynamics. The size of the stochastic torques is related to the temperature and the damping constant, such that the system tends towards thermal equilibrium for the chosen temperature. The approach also gives the dynamics in the limit of zero temperature, if needed: the stochastic torques are set to zero, but the damping can still be included.
In fact, it is more physical to think that the dynamics includes random magnetic fields, rather than torques. The dynamical equation for some selected unit dipole exposed to a deterministic field h and a stochastic field h s is written in the dimensionless quantities as The first term is the free motin and the second term is the Landau-Gilbert damping, with dimensionless damping constant α. The dynamics evolves according to the superposition of the deterministic and random fields. One can also split out their contributions separately, which shows the effective dynamical fields (in the brackets) acting on this dipole. Then in order for the stochastic fields to establish thermal equilibrium, they are assumed to have time correlations determined by the fluctuationdissipation (FD) theorem, The indices i, j refer to any of the Cartesian coordinates.
The dimensionless temperature T is the thermal energy scaled by the energy unit, The fluctuation-dissipation theorem indicates that the power in the thermal fluctuations is carried equivalently in the random magnetic fields. For reference, with the physical units replaced the relation can be expressed The Langevin equation in (18) is a first-order differential equation where the noise is multiplicative -the second term involves a product ofμ with a net stochastic field. To discuss the solution method, it is simplest to let y = y(τ ) be a vector that represents the entire set of spins, y = {μ i (τ )}. Then symbolically y obeys a differential equation in the general form, The function f represents the deterministic time derivative on the RHS of (18) and the function f s represents the stochastic part of the dynamics. Each is defined indirectly by comparing this with the Langevin equation.
The fields f , f s and h s are vectors of 3N components, where N is the number of dipoles in the array.

A. Second order Heun integrator
An efficient method for integrating this magnetic dynamics type of equation forward in time is the second order Heun method 23,24 . That is in the family of predictorcorrector schemes and is rather stable. Another method that could be applied with similar precision would be second order Runge-Kutta, however, it would be slightly less efficient for including the stochastic fields.
The predictor stage for the second order Heun algorithm is an Euler step, which is followed by a corrector stage that is equivalent to the trapezoid rule. Each involves moving forward in time over some time step ∆τ , with the needed results obtained by integrating Eq. 22 from an initial time τ n to a final time τ n+1 = τ n + ∆τ , during which the stochastic fields are acting. With notation y n ≡ y(τ n ), the predictor stage produces an initial solution estimateỹ n+1 at the end of one time step, y n+1 = y n + f (τ n , y n )∆τ + f s (τ n , y n ) · (σ s w n ). (23) The effect of the random fields is contained in the last term. The factor σ s w n replaces the time-integral of the stochastic magnetic fields. For each site l of the array, there is a triple of unit variance, zero mean random numbers (w x ln , w y ln , w z ln ) produced by a random number generator. The physical variance σ s needed in the stochastic fields is defined by an equilibrium average over the time step. For an individual component at one site, that is When the FD theorem is applied to one dipole, this gives the variance of the random fields, Thus, the individual stochastic field components, integrated over the time step, are replaced by random numbers of zero mean with the variance σ s . Then, in the corrector stage, the points y n andỹ n+1 are used to get better estimates of the slope of the solution. Their average effect is used in the trapezoid correction stage: The error varies as O((∆τ ) 3 ), hence it is a second order scheme. It is important to note that the same random numbers w n are used in this corrector stage as those applied in the predictor stage, for this individual time step.
In the actual implementation of this method, the functions f and f s do not need to be specifically identified.
Instead, it is seen that the change in any spin over a time step, ∆μ = dτ d dτμ , depends linearly on h∆τ from the deterministic part of the differential equation, and linearly on dτ h s (τ ) from the stochastic part. But the stochastic contribution is replaced by random numbers of the correct variance, that is, Thus, the Euler predictor step is efficiently carried out by evaluating the combined deterministic plus stochastic field contributions, for an individual site, like μ =μ + ∆μ, where the needed effective field that updates this site is the deterministic/stochastic combination, This latter relation generates an effective field, based on the deterministic local magnetic field from expression (12), combined with a stochastic fluctuating impulse. The same type of combination applies in the trapezoid corrector step. The updating field at the end of the time step is calculated, using the predicted position, together with the same random field, That leads to a different estimate for the spin change, Then the corrector stage gives the updated spin according to their averagê µ(τ + ∆τ ) =μ(τ ) + 1 2 ∆μ + ∆μ .
This algorithm does not ensure the conservation of spin length. Thus, the length ofμ can be rescaled to unity after the step. However, in practice, when that becomes necessary, the scheme is already becoming unstable, which can happen if the net field | g| is greater than 1. The best way to avoid numerical instability is to make sure the time step is sufficiently small, so that | g| ≪ 1, or even better, |∆μ| ≪ 1, to the precision desired. The integration requires a sequence of quasi-random numbers (the w n stochastic fields ) with a long period, so that the simulation time does not surpass the period of the random numbers. We have used the gnerator mzran13 due to Marsaglia and Zaman 25 , implemented in the C-language for long integers. This generator is very simple and fast and has a period of about 2 125 , due to the combination of two separate generators with periods of 2 32 and 2 95 .

B. The dipole fields
The calculation of the dipole term in the local magnetic field, (12), easily consumes most of the calculational effort. In order to have correct results, the fields at all dipole sites must be calculated both in the predictor and corrector stages of a step. We consider a system with open boundaries, rather than periodic boundary conditions. In this way, the dipole calculation is simpler; there are N (N − 1)/2 dipole field contributions to be found at any time. Some gain in speed can be achieved by computing tables of the pair displacements r ij and the powers r 3 ij that are being re-used during the calculations. One of the best ways to speed up the calculation of the dipole fields for larger systems is to write their calculation as a convolution of a Green's function with the source dipoles, and calculate that convolution in reciprocal space, transforming between real and reciprocal space with a fast Fourier transform (FFT). This is simplest if the locations of the dipoles form a simple square grid. Although that is true here, that grid is rotated at 45 • to the x-axis, which brings some complications and only works for the square lattice spin ice.
Instead, we consider that the spin ice involves unit cells on a square lattice, where each cell has a two-atom basis (i.e., two-dipole basis). For the cell whose lower left corner is at position r k = (x k , y k ) = (m k , n k )a, define the two dipoles present. On the "1" sublattice, and on the "2" sublattice, If there is an arbitrary source dipole µ at the origin, then the dipolar field h d it creates at position r = (x, y, z), according to the first term in (12), is well-known, (with all distances measured in lattice constants). This can be used to get the field produced from either sublattice. Summing over source dipoles, the 3 × 3 matrix is a Green's operator G( r) acting on the dipoles at discrete lattice sites. (Here the tilde is used only to indicate a 3 × 3 matrix quantity.) However, to account for the twoatom basis, the Green's matrix is expanded to have an extra pair of indices that refer to the sublattice, one for the field point (α) and one for the source point (β). The Green's matrix for the field produced at point r kα due to the source dipole at point r lβ is It will be important to keep in mind that G αβ actually depends only on the differences of the unit cell positions, r kl ≡ r k − r l . Now the dipole field on the α sublattice, for the cell at r k , is given by a discrete convolution Using fairly obvious notation, µ β ( r l ) is the dipole on the β sublattice for the unit cell at r l . The dot operation represents the matrix multiplication, i.e., an implicit sum over the Cartesian components of G αβ and µ β . Written this way, the same formula could apply to other lattices of interest, such as honeycomb, Kagomé, etc. Note that even at r k − r l = 0, there are contributions that must be included, corresponding to the interactions between the sublattices within an individual unit cell. For a specific example using the square ice sites, one can see that one particular interaction involving the two different sublat-tices (α = 1, β = 2) has a Cartesian element, This is nonzero when r k = r l . Also, the element G xx 21 ( r k , r l ) with source and observer sublattices interchanged can be obtained by changing the sign of a; they are not the same. A similar term with source and observer on the same sublattice is This is equal to G xx 22 ( r k , r l ). It is divergent at r k = r l , however, that would be an interaction of a dipole with itself, a term that must be excluded by definition.
G αβ depends only on the displacements, r kl ≡ r k − r l , which form another square lattice. Then one can find its Fourier transform, using a fast Fourier transform (FFT), setting the arbitrary source point to the origin. The Fourier transform of µ β is also determined. The convolution in real space becomes a simple product of G αβ and µ β in Fourier space, which can then be transformed back to real space by an inverse FFT to obtain h d . Although there is considerable overhead, for larger systems the speedup is tremendous (N ln N operations) when compared to doing the N sums with N terms to get the local dipole fields.
To apply the simplest FFT method, the size of the grid of primitive cells must be 2 Nx × 2 Ny with integers N x and N y . To avoid the wraparound problem, so that the system being simulated is really a single copy of the desired L x × L y size, one needs to choose N x and N y large enough so that 2 Nx > 2L x , and 2 Ny > 2L y . This ensures that the periodic copies of the system, inherent in the application of the Fourier transform, do not "see" or interfere with each other in the convolution.
It can also be kept in mind that there are some symmetries that reduce the calculational overhead. Displacements only on the 1-sublattice or only on the 2-sublattice are the same, so for any of its Cartesian components, Also, the matrix is symmetric in the Cartesian indices, for any sublattice indices, Furthermore, the interactions with both source and observer on the same sublattice are symmetrical in their interchange, Then this implies their Fourier transforms are pure real, leading to some reduction in the computations needed. On the other hand, there is no symmetry generally between different sublattices on different unit cells, so G 12 ( r kl ) = G 21 ( r kl ), see the discussion after Eq. (39). But the interactions within a cell must avoid selfinteractions, so we do define G 11 (0) = G 22 (0) = 0. The interactions between sublattices on the same cell depend only on squared displacements, so G 12 (0) = G 21 (0) = 0. For r kl = 0, these symmetries result in 12 independent elements in G( r kl ) (each of G 11 , G 12 , and G 21 have 4 independent elements), in contrast to the 4 independent elements needed for a single sublattice, see the matrix in Eq. (36).

IV. THERMAL EQUILIBRIUM PROPERTIES
Mostly the magnetic properties of spin-ice materials are investigated in an approximation of zero temperature, because the fundamental interaction strengths of the anisotropy energies and the dipolar energies are much greater than k B T at room temperature (Models A,B,C).
Even so, there could be energetic thermal fluctuations in a magnetic system even at low temperatures, in any situation where the magnetic fluctuations are large, such as near a reversal point in a hysteresis loop. This might possibly lead to enhancement of specific heat in such a situation, and of course, thermal rounding of the reversal paths in magnetization hysteresis loops. Thus it could be interesting to have some calculations of the energy, specific heat, and also of magnetic susceptibilities in a situation of thermal equilibrium.
The time evolution from the Langevin dynamics can be used to get thermal averages, as an alternative to Monte Carlo calculations. As long as the simulation time is long compared to any physical relaxation time, a sequence of energy samples E n and total system magnetization samples M n can be averaged, and their fluctuations can be used to estimate the specific heat and magnetic susceptibility. Suppose there are N s samples taken from the time evolution. The average energy for all N islands is estimated as with a "measurement error" estimated from its standard deviation σ E and the number of samples, The total heat capacity of the system is determined from the fluctuations in the energy, where β = (k B T ) −1 , and from that we obtain the specific heat per island, c = C N /N . The error in the heat capacity really needs to be estimated to be sure that the "measurement" is reliable. It can be calculated by finding the standard deviation of the quantity z ≡ (E − E ) 2 from which C N was obtained, which turns out to be found from the following averages, and the error in specific heat per island is ∆c = ∆C/N . Similar relations were used for the thermally averaged total magnetic moment of the system and the components of magnetic susceptibility. For instance, the susceptibility per island, χ xx , is found from fluctuations in total magnetic moment of the system, M x = l µ x l , and its error ∆χ xx is found from relations similar to those in (47) and (48) used for the heat capacity error, but with E replaced by M x .
Due to the fluctuations caused by the temperature in the simulations, the calculations of C and χ are generally not as precise as those of E and M , without making very long runs. Especially as mentioned above, these calculations are difficult in any physical situation where the magnetization is on the verge of reversal, where the fluctuations are greatest.
For these thermodynamics simulations, the system was started in a random state, with the temperature initially set at the highest value in the range of interest. For the chosen temperature, data samples were taken at some appropriate time interval, that depends somewhat on the energy couplings. For coupling parameters k 1 , k 3 on the order of 0.2 or less, and d several orders smaller (Models A,B,C), we found that a Heun time step ∆τ = 0.01 was sufficient to insure proper energy conservation at zero temperature. Using this time step for finite temperature together with α = 0.1, we averaged over N s = 4000 data samples separated by sampling time interval ∆τ s = 10 3 ∆τ = 10.0 . An initial time interval corresponding to 100 samples was allowed for relaxation before samples were taken. Simulations would be left to run even longer than 4000 samples, if necessary, until the percent error in the magnitude of the total system magnetization was found to be less than 0.1%. The final state at one temperature was then used as the initial state for the next lower temperature in the calculation.
For models D and E, the dipolar coupling is much stronger, and this requires a smaller Heun time step, ∆τ = 0.001, to insure proper dynamics at finite temperature and energy conservation at zero temperature. Besides this, the calculation parameters for averages were the same as for models A,B,C, e.g., N s = 4000 and taking samples at sampling time interval ∆τ s = 10 3 ∆τ , while waiting for 0.1% or better precision in the system magnetization.
Some thermodynamic results for the Wang et al. particles (Model A) are shown in Fig. 2, versus scaled temperature T . A 16 × 16 grid of cells was used (N = 2 × 16 2 = 512). As the dimensionless energy coupling constants are small numbers, the only interesting effects are observed for T < 0.1 . Near T ≈ 0.02 there are peaks in specific heat and in the in-plane components of magnetic susceptibility. As mentioned earlier, though, these features would appear only greatly above room temperature, which is marked with arrows. At these "high" temperatures, one would expect other modifications would take place (besides magnetic effects) and the model would not be applicable. Note that the monopole charge density in Fig. 2c does not go to zero at very low temperature here. It does, however, make a transition to a lower value. This is an indication that the system did not find a state close to a ground state. There is frozen-in disorder at lower temperatures. It is also an indication that the time scale for thermal relaxation to an equilibrium configuration was longer than the time interval used for averaging. Both the discrete and continuous definitions of ρ m exhibit similar behaviors. The order parameter Z (not shown) stayed close to zero for the whole temperature range shown. That is further indication of the system staying far from a ground state, where it would have reached one of the values ±1. Similar results for Models B and C, for smaller particles, are shown in Figures 3 and 4. These suggest that for the typical square lattice spin-ice using Py as the material, the room-temperature thermodynamics is nearly the same as that at zero temperature. There would be some limiting specific heat C ≈ k B and non-zero value for the in-plane susceptibility. However, the monopole density has extreme difficulty to go to zero while scanning from high to low temperature. Then, in fact, the dynamics is a low temperature dynamics in a disordered non-ground state (and non-equilibrium) configuration.
The thermodynamic results for 16 × 16 theoretical Models D and E are shown in Figures 5 and 6. Their general structure is similar: Model D has a strong peak in specific heat near T ≈ 0.22 and a more rounded peak in χ xx ≈ χ yy at a slightly higher temperature. The peak in specific heat in Model E is closer to T ≈ 0.23, and again the susceptibility peak is centered a little higher in temperature. For all of the models studied, the out-of-plane magnetic susceptibility χ zz is considerably smaller than χ xx , and there is only a weak temperature dependence. Models D and E, notably, do reach thermodynamic equilibrium in the simulations. This is seen clearly in the plots of the order parameters Z and ρ m . Now ρ m , for both discrete and continuous forms, tends towards zero at low temperature, as expected for the system moving towards a ground state. Further, the ground state overlap order parameter, Z, tends to go towards unity as T → 0; this is the strongest indication of approaching one of the ground states. At higher temperatures above the peaks in C and χ we see that Z becomes quite close to zero; the system is more random and far from a ground state. In the same high-temperature region, the monopole density tends towards some positive constant limiting value. The discrete definition for ρ m apparently undergoes a stronger change in value as the system makes a transition from its low to high temperature behavior.

V. HYSTERESIS CALCULATIONS
A simple experiment to investigate the magnetic properties of spin ice is the response in an applied external field (hysteresis calculation). To get a general impression of the physical response in any square spin ice, hysteresis calculations were carried out for Model D at fixed scaled temperature T = 0.1 . These were calculated the same as for the thermodynamics, except that it was adequate to average over shorter sequences, N s = 1000, at each step of the applied field. The system was initially set in a random configuration, but with the maximum positive applied field. The field was scanned to lower and negative values along some axis (eitherx or at 45 • to +x) and then allowed to come back to the starting value. In order to interpret the results, it was also important to calculate the order parameter Z and the monopole charge density ρ m during the hysteresis scan.
The results for field applied along thex axis are shown in Figures 7 and 8. In fact, at this temperature, the model does not exhibit any hysteresis: the magnetization per island is the same in backward and forward scans of h ext . However, the magnetization shows regions with distinctly different slopes. In Fig. 7b, one sees that the order parameter Z, however, tends to take on either values close to zero, at strong applied field, or, values near Z ≈ ±1, at weaker field. We note that this temperature T = 0.1 is on the low side of the specific heat peak for zerro magnetic field. Then this shows that in the central region of the MH graph, the system falls into states that are close to the ground states. These states, however, are slightly modified due to tilting of some of the dipoles according to the field strength. Hence, there is close to a linear response with h ext , as the dipoles on the 2 nd sublattice, which are nearly perpendicular to the field, get tilted by it.
The variation of monopole charge density with applied field, for the simulation in Fig. 7, is shown in Fig. 8. Both the discrete and continuous definitions are displayed. The discrete definition has more dramatic changes (although, its errors, not shown, may be larger). Especially, ρ m tends to zero (or a small value for the continuous version) over the same applied field range where Z ≈ ±1. This confirms clearly that the central region of the MH graph corresponds to the system being in states close to the ground states.
For applied field along an axis at 45 • to the +x-axis, the situation is similar, see Figures 9 and 10. The MH and ZH graphs are nearly the same as for those for applied field alongx. In this case, however, the applied field must be causing both sublattice dipoles to tilt at stronger fields. There is a difference, then, in the charge density plot, see Figure 10. Again, in the central region near weak h ext , the monopole density tends to zero, as expected for the system being close to one of the ground states. At strong field, however, the discrete charge density also tends towards zero (as does Z). In this case, though, as the dipoles all tend to align at 45 • to ±x, the net number of dipoles pointing into any charge site is then forced to be zero. This clearly forces the monopole density found by the discrete definition to zero. One sees that the monopole density by the continuous definition, on the other hand, does not fall to zero at high field, and instead behaves the same as it does for field applied alonĝ x. 10 ε, (a) the averaged magnetization per site versus external magnetic field hext applied at 45 • above the xaxis; (b) the order parameter Z versus hext. The field strength was initially set at hext = 0.8, and scanned to hext = −0.8, then back to the starting value as in a hysteresis loop calculation.

VI. DISCUSSION AND CONCLUSIONS
We have studied the possibilities of spin dynamics in artificial spin ice systems. To do this, we have analyzed frustrated arrays consisting of two-dimensional square lattices of elongated magnetic nanoislands. The internal structure of the magnetic nanoislands was taken into account in such a way that the dynamics of all spins inside each magnet were inserted in the interactions. Then, depending on the island shapes, aspect ratios, sizes, elements and organization in the array, we have looked for possible departures from the Ising-like behavior commonly found for these systems. In general, it is not easy to build square ice arrays with dynamics by using nanoislands made of permalloy. As a consequence, we have found that the systems without real dynamics (islands practically with an effective Ising behavior) have great difficulty to achieve the ground state (see, for instance, models A, B and C). Indeed, the order parameter Z (defined in Section II) never reaches the values −1 or +1 (the two degenerate ground states), even for very low temperatures. This result agrees with all experimental studies 6,13 concerning the square lattices. On the other hand, by considering fictitious materials (as described by the constants D, K 1 and K 3 ), we found interesting deviations from the Ising behavior and consequently, more easily thermalized spin dynamics along the array. For these systems (see models D and E), the ground state can be easily obtained for low temperatures.
Some of the results obtained here can be directly compared to those of Ref. 20, where the artificial square spin ice with point-like dipoles with Ising-like behavior was studied by using conventional Monte Carlo simulations. Here, for simplicity, the model which replaces the magnetic moment of the islands by a point-like and Ising-like dipole (see, for instance, Refs. [14][15][16]20 ) will be referred to as model F . In model F , the ground state (Z = ±1) of the square ice can be easily obtained for very low temperatures, i.e., by using conventional Monte Carlo simulations; the ground state naturally appears at low enough temperatures. However, our results for models A, B and C, obtained by using a Langevin dynamics, showed that the ground state does not appear at low temperatures. This may indicate that there is a kind of dynamic constraint (effectively, excessively long relaxation time) in the system's dynamics that prevents it from reaching the ground state over a moderate time of observation. Indeed, a similar result was found when the dynamics of the model F in the presence of external magnetic fields is considered 12,[27][28][29] . Thus, it may indicate that arti-ficial square spin ices made with Permalloy may never reach the ground state, since both external field dynamics and thermally driven dynamics have bottlenecks that prevent the access of the ground state. On the other hand, our results for models D and E, for materials with fictitious constants, indicate that it may be possible to access the ground state using another kind of material. It seems then that, by changing the island's anisotropies and interactions, the dynamical bottleneck is eliminated.
This difference may be associated with the way that the system explores the phase space. First, let us consider systems where the islands have an Ising-like behavior. Since we are dealing with classical particles, each island must pass through an energy barrier to change its magnetization direction (even being Ising-like), either by the creation and propagation of a domain wall or by rotation of a single domain. Either way, there is no option for tunneling and some energy must flow to the island. This internal energy barrier was not taken into account in the Monte Carlo calculations of Refs. 14,15,20 and probably this is why the ground state was obtained. The results of the present study and those from Refs. 12,27-29 include the energy barrier for spin flips and we may expect that the existence of a huge energy barrier is responsible for the difficulty to access the ground state. Moreover, in models D and E, where the energy barrier is smaller, the ground state is accessible. It seems then that the key factor blocking access to the ground state of artificial square spin ice is the energy barrier for spin flips.
For model F , the specific heat exhibits a peak at a temperature around 7.2D/k B , where it is suspected that the string connecting the Nambu monopoles is broken and the system is able to support free monopoles 20 . The specific heat for models A, B, C, D and E also exhibits a characteristic peak. However, for models A, B and C, the peak arises for much higher temperatures (around 20D/k B ); it may seem unreasonable since the magnetic moment of the islands for these models has many more degrees of freedom than that of model F . Nonetheless, we remark that the ground state can not be obtained for models A, B, C even for zero temperature and, therefore, there is not a way of establishing equivalences between the results of these models and results of model F . Alternatively, the ground state of models D and E can be found at zero temperature and then, we can compare the results for the specific heat with the ones obtained for model F . Really, the peak in the specific heat for models D and E arises at a temperature around 2D/k B , which is about three times smaller than the case of model F . It is obviously expected because there are more spin degrees of freedom for models D and E than for model F . Furthermore, the peak for model E occurs at a temperature slightly above that one for model D because the out-of-plane spin degree of freedom is more restricted in E.
To complete this study, we have also calculated the hysteresis for the models that exhibit dynamics. It is an important calculation to get a general impression of the physical response in any square spin ice, since the simple experiments to investigate the magnetic properties of these systems is the response in an applied external field. Finally, we would like to say that we really hope that all the investigations developed here could help experimental advances toward spin ice systems in which the ground state could be achieved and/or the transition rendered by appearance of free monopoles occurs around room temperature. Experimentally, a recent work was already addressed in this direction. Indeed, Kapaklis et al. 19 have proposed an experimental system (in an external magnetic field) where thermal dynamics can be introduced by varying the temperature of the array. On a square lattice, they use a material (based on δ-doped P d(F e)) with an ordering temperature near room temperature to confirm a dynamical "pre-melting" of the artificial spin ice structure at a temperature well below the intrinsic ordering temperature of the island material. Such a procedure is capable of creating a spin ice array that has real thermal dynamics of the artificial spins over an extended temperature range 19 . The possibility of observing emergent monopoles is therefore conceivable, following the general approach that the authors of Ref. 19 described in the design of spin ice arrays. This is a first step towards realization of artificial spin ices as conceived in models D and E.