TVE 14 028 juni
Examensarbete 15 hp Juni 2014
Study of Magnetic Nanostructures using Micromagnetic Simulations and Monte Carlo Methods
Nils Bäckström
Jonathan Löfgren
Vilhelm Rydén
Teknisk- naturvetenskaplig fakultet UTH-enheten
Besöksadress:
Ångströmlaboratoriet Lägerhyddsvägen 1 Hus 4, Plan 0
Postadress:
Box 536 751 21 Uppsala
Telefon:
018 – 471 30 03
Telefax:
018 – 471 30 00
Hemsida:
http://www.teknat.uu.se/student
Abstract
Study of Magnetic Nanostructures using
Micromagnetic Simulations and Monte Carlo Methods
Nils Bäckström, Jonathan Löfgren, Vilhelm Rydén
We perform micromagnetic simulations in MuMax3 on various magnetic
nanostructures to study their magnetic state and response to external fields. The interaction and ordering of nanomagnetic arrays is investigated by calculating the magnetostatic energies for various configurations. These energies are then used in Monte Carlo simulation to study the thermal behaviour of systems of nanomagnetic arrays. We find that the magnetic state of the nanostructures are related to their shape and size and furthermore affect the emergent properties of the system, giving rise to temperature dependent ordering among the individual structures. Results from both micromagnetic and statistical mechanic simulations agree well with available experimental data, although the Monte Carlo algorithm encounter problems at low simulation temperatures.
Ämnesgranskare: Jonas Lindh
Handledare: Vassilios Kapaklis, Erik Östman
Summary in Swedish
Mikromagnetism ¨ar idag ett viktigt omr˚ade inom forskning med m˚anga anv¨andnings- m¨ojligheter, inte minst inom informationslagring som i exempelvis h˚arddiskar. Mikromag- netisk forskning handlar om att studera hur sm˚a magneter, ofta i storleksordning mikro- eller nanometer, beter sig och interagerar med varandra. P˚a grund av sv˚arigheter i att direkt observera och kontrollera dessa sm˚a strukturer anv¨ands ofta datorsimuleringar f¨or att h¨arma de verkliga processerna.
Liksom alla system i naturen tenderar magneterna att r¨ora sig mot ett tillst˚and med s˚a l˚ag energi som m¨ojligt, det s˚a kallade grundtillst˚andet. Det visar sig att en viss typ av mikromagneter har grundtillst˚and som liknar en vanlig stavmagnet, med en nord- och sydpol. Ett exempel visas i figur 1(a), d¨ar de f¨argade omr˚aderna ¨ar magneter och deras riktning ges av pilarna. Denna egenskap g¨or att magneterna l¨ampar sig v¨al f¨or att skapa modeller av andra fenomen i naturen, som inte n¨odv¨andigtvis har med magnetism att g¨ora.
I figur 1(b) visas hur v¨ate- och syreatomerna i vanlig is strukturerar sig. Isen uppvisar h¨ar ett fenomen som kallas frustration, vilket betyder att systemet inte har ett v¨aldefinierat grundtillst˚and, utan ist¨allet flera olika. Detta leder till en viss slumpfaktor i hur systemet ordnar sig. Frustration ¨ar emellertid sv˚art att studera direkt i naturligt f¨orekommande material som is, eftersom detta kr¨aver att man kan g¨ora observationer p˚a atomniv˚a. En l¨osning ¨ar att bygga modeller med de f¨orh˚allandevis stora mikromagneterna.
(a) (b) (c)
Figure 1: (a) Fyra mikromagneter i en ordnad struktur. (b) Molekylstrukturen i vanlig is.
(c) Ett rutn¨at med upprepning av (a).
Upps¨attningen med fyra magneter i figur 1(a) motsvarar en syreatom (bl˚a sf¨ar) i figur 1(b). Om man upprepar denna upps¨atning magneter i ett rutn¨at enligt figur 1(c) f˚ar man ett frustrerat system som ¨ar f¨orh˚allandevis enkelt att finjustera och g¨ora m¨atningar p˚a. Dessa system ¨ar inte bara intressanta ur en modellsynpunkt, utan ¨ar ¨aven en m¨ojlig kandidat till framtida informationslagring. Systemen kan ¨aven anv¨andas f¨or att studera magnetiska laddningar, som kan m¨ojligg¨ora nya s¨att f¨or informations¨overf¨oring i framtiden.
I detta arbete anv¨ands datorsimuleringar f¨or att simulera mikromagneter och f¨or att studera frustrationen i de st¨orre modellsystemen. V˚ara resultat visar att egenskaper hos de enskilda magneterna, s˚a som dess storlek och avst˚and mellan dem, har stor p˚averkan p˚a hur de frustrerade systemet ordnar sig. ¨Aven tempereraturen hos systemet p˚averkar dess grundtillst˚and.
Contents
1 Introduction 1
1.1 Objective . . . 1
2 Theoretical background 2 2.1 Magnetic islands and their magnetization . . . 2
2.1.1 Magnetostatic interactions and the forming of domains . . . 2
2.1.2 Hysteresis . . . 3
2.2 Spin ice . . . 3
2.3 Vertices and artificial spin ice . . . 5
2.4 Micromagnetic simulations . . . 6
2.4.1 MuMax3 . . . 6
2.5 Equilibrium Monte Carlo methods . . . 7
2.5.1 Markov processes . . . 7
2.5.2 The Metropolis algorithm . . . 8
2.5.3 Equilibrium and sampling . . . 9
3 Method 10 3.1 Simulating micromagnetic islands . . . 10
3.2 Simulating vertices . . . 10
3.3 Monte Carlo simulation . . . 11
3.3.1 Algorithm details . . . 11
4 Results 12 4.1 Domain formation in micromagnetic islands . . . 12
4.1.1 Hysteresis in circular islands . . . 12
4.2 Vertex energies . . . 13
4.3 Monte Carlo simulation . . . 14
4.3.1 Varying temperatures . . . 15
4.3.2 Fixed temperature . . . 15
5 Discussion 18 5.1 Formation of single domain . . . 18
5.2 The effect of neighbouring islands and vertices . . . 18
5.3 Vertex energy with varying lattice constant . . . 19
5.4 Monte Carlo simulation . . . 19
5.5 Suggested further development . . . 20
6 Conclusions 21
1 Introduction
In materials science, recent advances in lithographic patterning techniques have opened up the possibility to nano-engineer new types of magnetic materials. As a consequence, com- pletely new areas of study have emerged which demand deeper understanding of magnetic properties and interactions of materials on the nanometer scale. Measurements performed on these new materials may help in understanding the physics of certain phenomena, such as frustration. Micromagnetic simulations provide an alternate way to further investigate the structure of these materials, and to discover new areas where experiments may be performed in future research.
The phenomena of frustration may for example be found in ordinary water ice, in which the arrangement of the atoms give rise to pairwise interactions that may not all be satisfied simultaneously. As a result, a high number of degenerate ground states exist leading to a built in entropy even at zero Kelvin. This, in turn, give rise to emergent properties and ordering in the material structure. Artificial spin ice is a model system used to study the physics of frustration, in which numerous magnets in the nanometer length scale are used to model the frustrated geometries.
The study of these systems may help developing front-end technology, where magnets of decreasing size are used in technological devices such as hard disk drives. Other possi- bilities exist as well, including the use in magnetic transistors. Furthermore, researchers have discovered that two magnetic poles may be separated through the spin ice structure, effectively creating artificial magnetic monopoles.
1.1 Objective
The scope of this work is limited to simulating micromagnets and larger model systems containing those. We will use micromagnetic simulations to study the static behaviour of the micromagnets, such as magnetic domain formation, stray field minimization and remanent magnetization. We will also study interactions between the individual magnets to estimate energies for the building blocks of the system, and use them in statistical mechanics simulations of large systems using a Monte Carlo (MC) algorithm. Results from these simulations are then compared and discussed in the framework of recent research in the field.
2 Theoretical background
Although the scope of this work is limited to performing simulations, some knowledge of magnetism and magnetic materials is necessary for the reader. In the following section, some introduction on this topic is presented. For further reading, we refer to textbooks on the subject, for example by David Jiles [1]. Following this introduction, micromagnetic structures and systems will be described. Finally, we look at the basics of MC methods and why they are so useful.
2.1 Magnetic islands and their magnetization
A magnetic island is an isolated patch of a magnetic material with a certain shape and size. If made small enough, the minimization of energy may cause the magnetization of such islands to form a single domain. The magnetization may for example be uniform in a certain direction, or in the shape of a vortex, as shown in figure 2(a)-(b). The preferred axis for uniform magnetization is known as the easy axis. If a single domain magnetization is not formed, the magnetic island is said to be multi domain, which can be seen in figure 2(c). Between the internal domains, domain walls are formed.
A more detailed description of the definitions above is given in the section below.
(a) (b) (c)
Figure 2: (a)-(c) Magnetic islands with different shape, size and magnetization: (a) 450 × 30nm circular island with vortex single domain magnetization. (b) 470 × 160 × 12nm elongated rod with uniform single domain magnetization. (c) 2820 × 960 × 12nm elongated rod with multi domain magnetization.
2.1.1 Magnetostatic interactions and the forming of domains
To understand why magnetic domains are formed, one must understand the different in- teractions in magnetic materials. In ferromagnetic materials there are long range magne- tostatic interactions and short range quantum mechanical exchange interactions between atomic spins [2]. The material will always try to minimize its total energy, and the com- petetive effects of the interactions determine the magnetization in the ground state. The minimization of the ferromagnetic exchange energy will align nearby spins parallell to each other. This in turn creates a larger stray field which increases the magnetostatic energy.
If a magnetic island is not symmetrical, the created stray field will depend on the direc- tion on the magnetization. This gives rise to shape anisotropy, which means one axis of magnetization is energetically favourable, the easy axis.
Since the exchange energy is short range it will tend to act on nearby spins, forming magnetic domains in which the magnetization is uniform. The minimization of the long range magnetostatic energy will then try to align neighbouring domains anti-parallel. Be- tween domains a domain wall is formed, also trying to minimize its own energy which results in different types of domain walls, as seen in figure 3. Domains will form when the decrease in magnetostatic energy is larger than the energy cost of forming domain walls, which depends heavily on the geometry and material parameters.
Figure 3: The striped white regions represent the two different types of domain walls in magnetic materials. The magnetization in the domain walls is perpendicular to the neigh- bouring domain magnetizations, and out of plane (Bloch) or in-plane (N´eel) respectively.
2.1.2 Hysteresis
When a ferromagnetic material is exposed to an external magnetic field, the internal spins will prefer to align in the same direction as the field. When the field is removed, some of the alignment will remain resulting in a remanent magnetization, or remanence. To remove the magnetization in the material, an external field in the opposite direction must be applied. The strength of the field required to do so depends on the materials coercivity, where high coercivity corresponds to high resistance in reversing the alignment of the spins.
Coercivity is a complex subject, highly related to the magnetic domains and domain walls discussed in 2.1.1 and dependent on the atomic and molecular structure of the material.
If the external field is varied over a loop and the magnetization is measured, a hysteresis curve can be obtained [1]. An example is shown in figure 4.
2.2 Spin ice
When water freezes, the molecular structure shown in figure 5 is formed. This structure has a built in frustration, caused by pairwise interactions in the substance which cannot all be satisfied at the same time. The hydrogen atoms are arranged on the lines connecting the oxygen atoms in the tetrahedical structure and for each oxygen atom, two hydrogen atoms are arranged closer to it and two further away. This is referred to as the ice rule, and can be satisfied in a number of ways resulting in residual entropy even at very low
−HC
Hs
H B
BR
Figure 4: A hysteresis loop showing how the magnetization of a ferromagnetic material responds to an applied external field. The external field begins at zero, increases until the material has reached its saturation magnetization, meaning that practically all domains are aligned with the field. HS is the the strength of the external field required to reach saturation. The field is then reversed and finally back to the starting direction. BR is the remanence magnetizaion showing the magnetization still left when the external field is removed after saturation. −HC is the coercivity which is the strength of the field needed to remove the magnetization.
temperatures. This was shown theoretically for water ice in 1935 [3] and later verified experimentally [4].
Figure 5: The molecular structure of water ice [5]. For every molecule, two hydrogen atoms are arranged near the oxygen atom and two farther away. This is illustrated by arrows:
two pointing in, two pointing out. The many ways in which this condition can be satisfied give rise to geometric frustration.
A magnetic analogue to water ice is spin ice, found in certain pyrochloric [6] com- pounds. The general structure is the same, although the frustration in spin ice arises from individual ion spins and their attempt to minimize their interaction energy. The anisotropy in the material forces the spins to point in one of the two directions, and since the inter- action energy between all pairs can not be minimized at the same time, frustration and a
high number of ground states exist. If thermal fluctuations are present, individual spins can be flipped in such way that the ice rules are not satisfied. This results in so called defects, where individual tetrahedra show a magnetic net charge. These defects show up in pairs of what can be seen as an monopole and an anti-monopole.
However, water ice and spin ice are difficult to study because of the length scales involved and the difficulty in observing individual spins and atoms. This begs the question:
is it possible to mimic these structures to make it easier to study the physics of frustration?
2.3 Vertices and artificial spin ice
The single domain magnetic islands described in section 2.1 can be seen as macro spins only able to point in one of two directions. By placing four such islands in a vertex configuration, the molecular structure of spin ice can be modeled. The vertices represent a 2D projection of the tetrahedra and each vertex connects four spins, just as in spin ice.
When placed on a lattice, these vertices make up artificial spin ice.
Artificial spin ice vertices can be ordered into four different types (figure 6), depending on how many macro spins are pointing in and out of its center. Since type I and II vertices have equally many macro spins pointing in as out, they have the lowest energy levels.
Following the same reasoning, type III and IV vertices have the highest energy levels [7].
Artificial square ice verticies may be placed on a lattice, as shown in figure 7, where the lattice constant is defined as the distance between the centers of two neighbouring vertices. As a result of this placement, each pair of adjacent vertices share one magnetic island, meaning that changing its magnetization will change the energies of both vertices.
Figure 6: All sixteen artificial spin ice vertex configurations [8], divided into the four different types (I, II, III, IV). All same-typed vertices have degenerate energy levels.
The particular structure mentioned above is known as square artificial spin ice.
This structure is similar to water- and spin ice, but easier to manipulate and analyze.
In practice, these nanoscale magnets are fabricated using lithographic techniques and can be observed in detail using magnetically sensitive imaging techniques. The fabrication gives the possibility to tune the energy scales and geometry of the system by changing the sizes and material of the magnetic islands. To simulate systems like this a combination of micromagnatic simulations and MC methods may be used.
Figure 7: Artificial spin ice vertices placed on a lattice [9].
2.4 Micromagnetic simulations
Micromagnetic simulations provide the means to describe and understand the magnetism of physically small magnetic objects. For many micromagnetic systems, analytical solutions are not available making numerical simulations necessary. Historically, the Object Oriented MicroMagnetic Framework (OOMMF) [10] has been widely used for simulations of this kind. Since OOMMF uses the CPU for calculation, the performance is greatly limited by the speed of the processor used. Many of the computations are highly parallelizable, and even though later versions of the program supports parallellization, GPU-accelerated programs can improve performance due to a much larger number of processing cores.
This increase in performance means simulations can run with finer discretizations, thereby improving the accuracy of the results.
2.4.1 MuMax3
MuMax3 is a GPU-accelerated simulation program that can simulate both magnetostatic systems and dynamic systems in time and space. For large systems, a speedup of a factor 100 and higher compared to CPU-based programs can be obtained [11].
The differential equation used to predict the evolution of the magnetization in ferro- magnets is the Landau-Lifshits Equation (LLE)
∂M
∂t = − γ
1 + α2M × Hef f − αγ
Ms(1 + α2)M × (M × Hef f) (1) where Ms is the saturation magnetization, γ the gyromagnetic ratio and α the damp- ening factor. Hef f is the effective field, which combines several field terms (one of which is the external magnetic field) with additional effects such as the possible anisotropy and thermal effects. MuMux3 solves the LLE numerically by discretizing the magnetic sample into cubiodal cells with a uniform magnetization, and then integrating the equation using a Runge-Kutta time stepping scheme. Since MuMax3 uses a method based on the fast fourier transform to calculate the magnetostatic field term, space has to be divided into
equal cells on a regular grid. Because of this, geometries other than perfect rectangles has to be approximated. However, smaller cell sizes can be used to improve this approximation.
For more information on the theory behind MuMax3 we recommend reading the intro- ductory article by A. Vansteenkiste [11]. Practical information and examples to get started using the program can be found on the official website [12].
2.5 Equilibrium Monte Carlo methods
1Monte Carlo methods are simulation techniques based on random sampling from a certain probability distribution. A model system is built and let to pass from state to state by a chosen dynamic that mimics a real system. For a correct choice of dynamics, the probability of the system being in a certain state at some time is equal to the probability weight of the real system, given by a Boltzmann distribution for a system with a thermal reservoir [13]. It states that the frequency or probability of a state with energy E is given by e−E/kT up to some normalizing constant, where k is Boltzmann’s constant and T is the temperature of the system. By sampling from the states the system passes through during the simulation, it is possible to get average values for different properties of the system.
It turns out this method give accurate estimates even when sampling only a small fraction of the states, if they are chosen in a correct manner. For example, a 10x10 island artificial spin ice model would have an order of 2100 states (not discounting for symmetries), so it is clearly impractical to calculate sums of all these states to get measurements for quantities. This is why MC methods are so useful. The main disadvantage is that by only taking some of the states into account, statistical errors are inevitable in the results.
The MC method can be thought of as a way of sampling states from the often astronom- ical number of states, according to a certain probability distribution. Sampling according to a Boltzmann distribution means the relative frequency of each state will correspond to the time spent in each state by the real system. Instead of averaging over all states, the values are calculated only from states which contribute most to this sum. This may seem clear, but the difficulty lies in choosing the dynamics such that the set of states gener- ated correspond to the Boltzmann distribution. One way of doing this is by using Markov processes.
2.5.1 Markov processes
A Markov process can be seen as a machine which, given a state of a certain system, returns a new state of the system. When this is repeated the Markov process generates a Markov chain of states, containing the states we wish to sample from when making our measurements. This process is necessary as we cannot simply choose a state at random and use the Boltzmann distribution to decide whether we should use it or not: almost all states would be rejected due to the exponentially decreasing probability in the Boltzmann distribution. A Markov process is a stochastic process, since it is characterized by a set of
1This section is based on the textbook Monte Carlo Methods in Statistical Physics [13] by Newman and Barkema, which is recommended for a more thorough discussion and analysis of the subject.
transition probabilities P (µ → ν), giving the probability of generating state µ given state ν. For a true Markov process, these probabilities should be constant in time and depend only on the given states, not on the history of the process.
There are two more important conditions for our Markov process, the condition of ergodicity and the condition of detailed balance. Ergodicity demands a nonzero probability of reaching any state from any other arbitrary state, given enough time. Detailed balance is the condition that guarantees that we have a Boltzmann distribution of states once we have reached equilibrium. Equilibrium means that the total transition rate into and out of any state is equal [13]. From this starting point, we come to the condition that pµP (µ → ν) = pνP (ν → µ), where pµ and pν are the relative frequencies of the states.
Combining this with our expression for the Boltzmann distribution gives the following condition for our transition probabilities
P (µ → ν)
P (ν → µ) = e−β(Eν−Eµ) (2)
where β = 1/kT , giving freedom in our choice of probabilities. Figuring out the mechanism of generating states in the Markov process to give the exact right transition probabilities is very difficult. To increase the freedom in the Markov process, we introduce the acceptance ratio and split our transition probability in two parts
P (µ → ν) = g(µ → ν)A(µ → ν) (3)
where g(µ → ν) is the probability that our Markov process will generate state ν from state µ or the selection probability, and A(µ → ν) is the probability that we will accept this state and update our model. This gives us more freedom in choosing our Markov process, since we always tune the acceptance probabilities to satisfy the condition in equation 2.
However, in an ideal algorithm, we would want the selection probabilities to be equal to the desired transition probabilities. This means all the acceptance ratios are equal to unity and no time is wasted rejecting moves. In fact, this is one the most important considerations when designing a MC algorithm. If the acceptance ratios are very low for most moves, the algorithm will slow down dramatically giving us fewer states to sample from. It may also make amount of steps needed to reach equilibrium so large that no results are generated from the simulation.
2.5.2 The Metropolis algorithm
As indicated in the previous section, acceptance ratios as close to unity as possible are highly desirable. The Metropolis algorithm is a widely used Monte Carlo method that achieves this in a systematic way. Here, we focus on methods with single-spin-flip dynamics, meaning that in a spin model, the Markov process only proposes moves which differ from the previous by one single spin. Hence, the non-zero selection probabilities are given by
g(µ → ν) = 1
N (4)
for a model with N spins. It is clear this single-spin-flip dynamic fulfills the condition of ergodicity discussed in the previous section, since we can move from any state to another by simply flipping each differing spin, one at a time. It also assures that the proposed new state has a small difference in energy from the previous one, which is analogous to a real system that rarely makes high energy transitions. The next step is choosing the acceptance ratios, keeping in mind that we must satisfy equation 2 and also keep them as close to unit as possible. This is the essence of the Metropolis algorithm, which gives the acceptance ratios
A(µ → ν) =
(1 if ∆E ≤ 0
e−β∆E if ∆E > 0 (5)
where ∆E = Eν − Eµ is the energy difference for the proposed change in states. These ratios are easily derived by maximizing them given the conditions of equation 2 and 3.
Since the selection probabilities cancels out in equation 2, we are left with a ratio between the two acceptance ratios. The largest one of these can then be set to unity and the other chosen to satisfy the condition.
The algorithm described above is a Metropolis algorithm for a spin system with single- spin-flip dynamics, because of its use of equation 5 for calculating the acceptance probabil- ities. In this work, the Metropolis algorithm is applied to artificial spin ice model system, which is one case in which the Metropolis algorithm prove useful. The algorithm does, however, have certain limitations.
2.5.3 Equilibrium and sampling
The goal of applying the Metropolis algorithm to our model system is to simulate the evolution of the system at a certain temperature, and eventually reaching equilibrium. As previously mentioned, this means that the final probability distribution of states is given according to their respective Boltzmann weights. In our model system, this is verified by looking at the distribution of the vertex types at the end of the simulation. A more naive way of determining if the simulation has reached equilibrium is to examine some property of the system over time, for example the energy, and check if it seems to have reached stable minimum. The simulation must then be run some time at equilibrium to gather states to sample for our measurements of the system.
As a measure of time in the simulation, the concept of a sweep is often used. In a spin model with N spins this corresponds to N selected new states, regardless of whether the state was accepted or not.
3 Method
The MuMax3 simulations presented in section 2.4 were run for elongated rods, circular islands and vertices to examine the domain formation and energy levels. The energy values obtained were used to simulate square artificial spin ice systems using the MC methods discussed in section 2.5.
3.1 Simulating micromagnetic islands
Firstly, simulations were made to confirm that the micromagnetic islands formed single domain magnetization for sufficiently small islands. For these simulations, island sizes 300 × 100 × 12nm, 470 × 160 × 12nm and 2820 × 960 × 12 were used. The geometry of an island was set up as a rectangle with semicircles conjoined to both ends, creating elongated rods. An external field was applied, either perpendicular to the islands easy axis, or with 45◦ angle from the easy axis. After the systems magnetization had aligned with the external field, the field was removed and the system was allowed to reach equilibrium.
In addition, two simulations with circular islands were made. One with a single circular island and one with four circular islands placed 550nm apart, all with a diameter of 450nm and thickness of 30nm. First, the external magnetic field was set to zero and initial magnetization randomized. The magnetization was allowed to reach equilibrium and a magnetic field was applied along one axis. The field strength was varied between −0.15T to 0.15T in small steps, letting the magnetization reach equilibrium in each step. The data for the magnetization was then used to create hysteresis curves for both systems.
3.2 Simulating vertices
After simulating individual islands, simulations of the vertex types were run for different lattice constants and island sizes. Only one of the degenerate states for each type were simulated. The islands were given an initial uniform magnetization and were then relaxed, in order to acquire vertex energies for the MC simulations. The energy values used are the values for magnetostatic energy.
The simulated island sizes and lattice constants are presented in table 1.
Table 1: Simulated island sizes and lattice constants for the artificial spin ice vertices.
Dimension LxWxH (nm) Lattice constants (nm)
300 × 100 × 12 [380 420 460 540 800 1200 2000]
470 × 160 × 12 [600 660 720 850 1400 2200 3600]
In addition, a simulation with two neighbouring vertices sharing one magnetic island was performed. In this simulation, the island size was 470 × 160 × 12nm and the lattice constant 600nm.
3.3 Monte Carlo simulation
The Metropolis algorithm, described in section 2.5.2, was used to simulate the artificial spin ice systems. Our implementation uses single-spin-flip dynamics where each magnetic island is considered a macro-spin with two possible states of magnetization along its long axis. The islands were placed on an n × n lattice similar to that of figure 7, although rotated 45 degrees. Periodic boundary conditions were used to mimic an infinite lattice.
The total energy of the state was approximated as the sum of all individual vertex energies, with energy values taken from the MuMax3 simulations. Since each island is connected to two vertices, the energy difference associated with a single-flip could be cal- culated as the energy difference of these vertices before and after the flip. This was our
∆E used in equation 5.
The algorithm above was repeated for a given number of sweeps, a fixed temperature and a set of vertex type energies. Simulations were run for a wide range of temperatures, and in order to determine whether the system had reached equilibrium, the evolution of the total system energy was analyzed. Finally, a verification that the type distribution corresponded to a Boltzmann distribution was made.
A more detailed description of the algorithm follows below.
3.3.1 Algorithm details
The algorithm was implemented in MATLAB R2013b. In the algorithm, a flip attempt began with the proposal of a new state, according to the selection probabilities discussed in section 2.5.2. In practice, this was done by taking one island in the grid at random with equal probability 1/n2, and propose to flip this island as our next state in the Markov chain. Thus, we only proposed states which differ by one spin as discussed in section 2.5.2.
We then used the acceptance ratio given in equation 5 to determine if we should accept this state as our next. We used the approximation for total system energy discussed in section 3.3 to calculate the energy change.
4 Results
4.1 Domain formation in micromagnetic islands
Figure 8 show the magnetization for elongated rods, relaxed from different initial states.
The magnetization fluctuates for some time while stabilizing in a single domain formation along the easy axis.
Note also, in figure 8 (d) and (h), that the magnetization at the edges of the rods is not entirely homogeneous.
(a) (b) (c) (d)
(e) (f) (g) (h)
Figure 8: MuMax3 simulation of 470 × 160 × 12nm elongated rods, showing the formation of single domain magnetization. (a)-(d) The magnetization when relaxed from a state with an external field perpendicular to the easy axis. (e)-(h) The magnetization when relaxed from a state with an external field with 45 degrees to the easy axis.
Figure 9 shows the relaxation of a much larger magnetic island. The island creates sev- eral magnetic domains showing that the domain formation in a magnetic material depends on shape and size.
(a) (b) (c) (d)
Figure 9: MuMax3 simulation of 2820 × 960 × 12nm elongated rods, showing the formation of multi domain magnetization. (a)-(d) The magnetization when relaxed from a state with an external field perpendiculat to the easy axis.
4.1.1 Hysteresis in circular islands
Results from the hysteresis simulation of the circular magnetic islands are presented in figure 10. As the external field increases the center of the vortex is pushed outwards perpendicularly to the applied field and the magnetization of the material is aligned with
the external field. As the external field decreases back towards zero, it does not take the same path along the magnetization curve. Since the magnetization curve always intersects the same point when the external field is zero, the magnet will always return to a single domain vortex when no field is applied.
(a) 0 mT (b) 34 mT (c) 92 mT (d) 93 mT (e) 38 mT (f) 4 mT
−0.15−1 −0.1 −0.05 0 0.05 0.1 0.15
−0.5 0 0.5 1
External magnetic field (T)
Relative magnetization
a) b)
c)d) f)
e)
(g)
−0.15−1 −0.1 −0.05 0 0.05 0.1 0.15
−0.5 0 0.5 1
External magnetic field (T)
Relative magnetization
(h)
Figure 10: (g) Hysteresis curve with sample points for (a)-(f). (h) Hysteresis curve for four circular islands, placed in square with a center-to-center distance of 550nm. Here, the islands switch from vortex to uniform and back at a lower field value than the single island.
Neither curve show remanence when the magnetic field is zero.
4.2 Vertex energies
Figure 11 shows MuMax3 simulations of the different types of vertices. The interaction between the islands can be seen at the edges of the islands where the magnetization bend away from the easy axis. In figure 11 (a), (b) and (c) the stray field at the edges align and decrease the energy for the vertex. In figure 11 (d) however, the stray field interacts by repulsion and the energy for the vertex increases. Note also that the magnetization of all individual islands are kept single domain after the relaxation even when interacting with nearby magnetic islands.
Energies for all vertices for varying distances and sizes are shown in figure 12. The energy for type IV vertices decreases with increasing distance, while the energy of all other types increases. Furthermore, if the individual islands are moved far from each other, no
(a) (b) (c) (d)
Figure 11: Magnetization of the four different vertex types, from MuMax3 simulations with an island size of 300 × 100 × 12nm. (a) Type I, (b) type II, (c) type III and (d) type IV vertices.
energy difference exists between the four types.
500 1000 1500 2000
1.9 2.3 2.7
x 10−17
Lattice constant (nm)
Demagnetizing energy (J)
Type I Type II Type III Type IV
(a)
500 1500 2500 3500
3.2 3.6 4 4.4 4.6x 10−17
Lattice constant (nm)
Demagnetizing energy (J)
Type I Type II Type III Type IV
(b)
Figure 12: Vertex type energies as a function of lattice constant, for (a) 300 × 100nm and (b) 470 × 160nm islands, from MuMax3 simulations. As the lattice constant increases, the coupling decreases until no interaction remains.
When two vertices of equal type were simulated together and thus sharing one island the energy was up to 15% lower, depending on the chosen type.
4.3 Monte Carlo simulation
In the following simulation results, temperatures in terms of relative energies of the system are more relevant than in units of Kelvin, because of the different energy scales in the simulations compared to experimental artificial spin ice systems. To generalize the results for other energy scales, temperature (T) is therefore weighted to energy properties of our system: the temperature in Kelvin (K) is multiplied by the Boltzmann constant (J K−1)
10−3 10−2 10−1 100 101 Temperature
Weighted distribution
Type I Type II Type III Type IV
(a)
10−3 10−2 10−1 100 101 6.8
7 7.2 7.4 7.6 7.8x 10−14
Temperature
Total energy (J)
(b)
Figure 13: Monte Carlo simulation of 64 × 64 lattice with island size 470 × 160 × 12nm and lattice constant 600nm. Last 400 step average of (a) vertex type distribution and (b) total energy of the system, as a function of temperature. At a temperature of just below 10−1.5, system dynamics are brought to our simulation method. Because of the weighted distribution in (a), the y-scale is omitted.
and then divided by the energy (J ) of a type I vertex. From this section and on, this dimensionless parameter will be the measure of temperature.
4.3.1 Varying temperatures
The results from MC simulations with varying temperature are shown in figure 13 and figure 14. In the former, we see the type distribution (a) and total energy (b) for the system. The presence of type III vertices start occurring for temperatures just below 10−1.5. Here, the systems energy and type distribution form deterministic curves as opposed to the stochastic behaviour that occur for lower temperatures. As the temperature rises further the type distribution evens out, until finally equally distributed between the types.
The formation of chains for different temperatures is presented in figure 14. In (a), the system is stuck in a local energy minimum. In (b), (c) and (d) the chain formation and number of exited types is related to the temperature. An increased temperature results in more chains and more high energy vertices.
4.3.2 Fixed temperature
The results from MC simulations for T = 3.3e − 2 are shown in figure 15. During the first few sweeps the energy of the system decreases rapidly and the distribution of types quickly favour the lower energy types. At the end of the simulation, when the system is in equilibrium, a Boltzmann distribution is indeed found in figure 15(c). Note also that the higher energy vertices in figure 15(d) does not occur at entirely random positions, but form chains which also run through the boundaries.
1 2 3 4
(a)
1 2 3 4
(b)
1 2 3 4
(c)
1 2 3 4
(d)
Figure 14: The vertex type distribution after a 10 000 sweep MC simulation for temper- atures (a) 3.3 × 10−4, (b) 1.9 × 10−2, (c) 3.3 × 10−2 and (d) 1.2 × 10−1. Brighter colours represent higher energy vertex types. The small number of chains in (b) corresponds to a lower system energy level, whereas the large number of high energy vertices in (c) and (d) correspond to a higher system energy. In (a), no type III vertices are present, and the system is stuck in a local energy minimum.
0 500 1000 1500 2000 6.9
7 7.1 7.2 7.3 7.4 7.5 7.6
x 10−14
Monte Carlo steps
Total energy (J)
(a)
0 500 1000 1500 2000
0%
20%
40%
60%
80%
Monte Carlo steps
Distribution
Type I Type II Type III Type IV
(b)
3.4 3.6 3.8 4 4.2 4.4
x 10−17 10−2
100 102
Type energy (J)
Weighted type count
Type I−IV Linear regression
(c)
1 2 3 4
(d)
Figure 15: Results of Monte Carlo simulation of 64 × 64 lattice with island size 470 × 160 × 12nm, lattice constant 600nm and T = 3.3 × 10−2. (a) Energy of the system and (b) weighted type distribution for first 2000 sweeps of the Monte Carlo simulation. (c) Average type count for last 400 sweeps and (d) graphical display of type distribution after 10 000 sweeps. Brighter colours represent higher energy vertex types.
5 Discussion
The results of the MuMax3 simulations do indeed show interesting properties, only observ- able on the nanometer lengthscale that was used. The simulations give precise energies, which seem to help in providing an accurate behaviour of the larger scale MC simulations.
However, the chosen MC algorithm does differ from experimental results for low simulation temperatures.
5.1 Formation of single domain
The formation of single domain magnetization along the easy axis in the smallest magnetic islands (figure 8) was vital to our work, as it makes sure the islands can be seen as macro spins in the artificial spin ice model. The larger structures showing multi domain magne- tization (see figure 9) is an example of the opposite phenomena, when the total energy is minimized by instead forming several domains. These results show that, as expected from domain formation theory [2], sufficiently small magnetic islands tend to be single domain whereas larger islands are multi-domain.
Our simulation also shows that the magnetization bends off from the easy axis at the edges to minimize the magnetostatic energy. This means that our simulated energies are lower than those with a completely uniform magnetization. This is a desirable result that give us more exact energies for the MC simulations.
The circular islands show single domain magnetization although in the shape of a vortex. Note that the hysteresis curve, in both cases, passes through the same point when the external field approaches zero. This means that the magnet is a stable single domain whose magnetization acquires the shape of a vortex when no external field is present. This is consistant with experimental work, for example by Shinjo et al. [14], which demonstrated the existence of vortices in permalloy dots with diameters between 300 nm and 1 µm.
When four circular islands were placed in the same alternating magnetic field, a similar behaviour was seen except for a somewhat smoother hysteresis curve, as seen in exper- iments [15]. The four islands seem to interact in a way that pushes the magnetization towards a vortex magnetization. This is also seen in the opposite direction where a weaker external field is needed to force the four islands to a uniform state. The increase in total magnetization for the four islands seem to amplify the effect of the external field.
5.2 The effect of neighbouring islands and vertices
For the elongated islands, the non-homogeneous magnetization at the edges becomes even clearer in figure 11, where the four magnetic islands in the vertex minimize their total energy not only with itself, but also with respect to the other three islands in the vertex.
Furthermore, the presence of single domain magnetization, despite the interaction with other islands, strengthen the argument that the single domain magnetization will remain even in larger system where more islands interact, such as those simulated with MC meth- ods.
We believe the 15% lowering in energy when simulating two adjacent vertices to be a combination of two sources. First of all, the system with two vortices sharing one of the islands has one island less that contributes the the total energy. Secondly, the shared island has both edges connected to a vertex with interacting islands. No further simulations were made for larger constellations of vertices in MuMax3, but in a lattice, where every vertex have four neighbours, it is reasonable to assume that the energy would be lowered even more.
5.3 Vertex energy with varying lattice constant
The change in energy that comes with varying the lattice constant (see figure 12), implies that the islands making up the vertices are coupled for short distances. If the islands are moved far away the energy difference disappears completely, suggesting the islands are now uncoupled. Furthermore, for type I-III the energy increases with increasing distance, which means that these vertex configurations are energetically favourable due to the minimization of the stray field. For type IV, however, the vertex configuration leads to a higher energy level. That is, the lower type vertices prefer to be in this formation, whereas the type IV vertices would rather be completely decoupled.
5.4 Monte Carlo simulation
The results of the MC simulations reveal several interesting phenomena in the artificial spin ice systems. They show how the system properties and short range interactions between the nanoscale islands create long range structure and order in the system. As seen in figure 13, this ordering is also highly dependant on the temperature and in agreement with previously performed MC simulations [16].
Because of the many possible states of the system, it was necessary to use MC simula- tions to find equilibrium. Since all individual parts of the system were pairwise coupled, and the way the system changed state depended on the predicted energy change, we needed an algorithm that could excite the system out of local minima while still fulfilling a final Boltzmann distribution related to the temperature. The Metropolis algorithm does both and did very well as long as the acceptance rate was high enough.
The introduced periodic boundary conditions were used in order to mimic an infinite size lattice. This makes sure that there are no special boundary islands with different properties in the system.
Moreover, the total energy of the system was approximated as the sum of all individual vertex energies. In absolute values this is a rough approximation but the relative differences between the energies should remain. Thus, the behaviour of the system is expected to be properly simulated in our algorithm.
The formation of type II chains seen in figure 14 is in agreement with experimental results [17, 18]. The length of the chains seem to be highly related to the minimization of energy, as shown in figure 13. Also, the occurrence of higher energy vertices increase as
the temperature rises. This comes as a result from the increased flip probability that come with higher system energies.
Clearly, a minimum exists around T = 3 × 10−2. For temperatures below, the system is ”freezed in” and stuck at a local minimum, thus being unable to further minimize the system energy as expected by experimental comparisons [19]. The problem arises when the probability for a flip to occur becomes so small, that the required number of simulation loops makes the algorithm highly impractical. To see the dynamics of the system at lower temperatures, the algorithm must be more effective in evaluating the flips that are not accepted.
An immediate consequence of the single-flip algorithm is the following. For a vertex of type I, II and IV, a single-spin-flip will always result in a type III. This is also valid the other way around, meaning that the appearance of type III is required for the system to have dynamics.
The presence of a Boltzmann distribution, as observed in figure 15(c), implies a correct final distribution of vertex types for the MC simulations when the system did not freeze in a local minima.
5.5 Suggested further development
During our MC simulations, we focused solely on the simulation of square artificial spin ice.
There are, however, other geometries as well. For instance, the circular islands simulated in this work could also be placed in a rectangular lattice and geometries such as Shakti [16] may also be simulated using MC methods.
The problem with the system freezing in due to low transition probabilities for low temperatures may be solved using a continuous-time algorithm. By doing so, time is increased in non-equal time steps, thus decreasing the need for long simulation run times.
With the current algorithm the system is only allowed to flip one island at a time, forcing dynamics to depend heavily on the type III vertices. Another approach is to implement a loop algorithm which evaluates several islands (two per vertex the loop passes) at the same time [20]. This algorithm would make the system less dependant on individual vertex excitations.
In order to achieve total system energies which correspond to real experiments, one must further investigate the influence adjacent vertices have on each other. If a correct parametrization can be found, more accurate energies could probably be obtained from the MC simulations.
6 Conclusions
• MuMax3 provides fast and accurate simulations of micromagnetic islands.
• The formation of magnetic domains in micromagnetic structures is highly related to their shape and size. Sufficiently small magnets will be single domain, where elongated rods will get uniform magnetization and circular islands will get vortex magnetization. When the size is increased, over some critical size, the magnets will be multi domain.
• If a circular magnetic island is placed in a varying external field, a hysteresis curve can be produced. However, for the sizes used in our simulations, no remanent mag- netization would remain when the field is turned off.
• For all vertex types, except for those of type IV, the coupling between the islands leads to a lower magnetostatic energy. If the lattice constant is increased, this effect is reduced. If several vertices are placed next to each other, the energy will be lower than the sum of the individual vertex energies.
• When square artificial spin ice is simulated with a metropolis Monte Carlo method, chains of type II vertices will form for a certain temperature interval. These chains will be surrounded mostly by type I vertices. As the simulated temperature increases so does the energy level of the system, resulting in the occurence of type III and IV vertices.
• All type transitions in the algorithm have to pass through a type III vertex. Thus, for low temperatures no transitions can occur due to lower energy levels, and the system freezes in at a local energy minimum. This is believed to be a result of the fixed step algorithm used which, in combination with the low acceptance ratios, leads to problems in reaching equilibrium.
• For sufficiently high temperatures, a final state where the occurence of vertex types corresponds to a Boltzmann distribution was reached. This indicates that the MC algorithm used in this work is valid.
References
[1] David C. Jiles. Introduction to magnetism and magnetic materials. Second edition.
CRC press, 1998.
[2] Alex Hubert and Rudolf Sch¨afer. Magnetic Domains, The Analysis of Magnetic Mi- crostructures. 3rd Printing. Springer-Verlag Berlin Heidelberg, 2009.
[3] L. Pauling. “The structure and entropy of ice and of other crystals with some random- ness of atomic arangement”. In: Journal of the American Chemical Society 57.2680 (1935).
[4] W. F. Gaiuque and J. W. Stout. “The entropy of water and the third law of thermo- dynamics. The heat capacity of ice from 15 to 273 K”. In: Journal of the American Chemical Society 58.1144 (1936).
[5] Steve T. Bramwell. “Condensed-matter physics: Great moments in disorder”. In:
Nature 427 (2006), pp. 273–274.
[6] S.T. Bramwell and M.J. Gingras. “Spin ice state in frustrated magnetic pyrochlore materials”. In: Science 294.5546 (2001).
[7] R. F. Wang et al. “Artificial spin ice in a geometrically frustrated lattice of nanoscale ferromagnetic islands”. In: Nature 439 (Jan. 2006).
[8] Vassilios Kapaklis et al. “Melting artificial spin ice”. In: New J. Phys. 14 (Mar. 2012).
[9] Unnar B. Arnalds. “Magnetic order in artificial structures”. PhD thesis. Uppsala University, 2012.
[10] M. J. Donahue and D. G. Porter. OOMMFF User’s Guide, Version 1.0. Tech. rep.
Gaithersburg, MD: National Institute of Standard and Technology, Sept. 199.
[11] A. Vansteenkiste and B. Van de Wiele. “MuMax: A new high-performance micromag- netic simulation tool”. In: Journal of Magnetism and Magnetic Materials 323.2585 (2011).
[12] DyNaMat group - Ghent University - Belgium. mumax3. May 2014. url: http : //mumax.github.io/.
[13] M. E. J Newman and G. T. Barkema. Monte Carlo Methods in Statistical Physics.
Oxford University Press, 1999.
[14] T. Shinjo et al. “Magnetic vortex core observation in circular dots of permalloy”. In:
Science 289.5481 (2000).
[15] R.P. Cowburn et al. “Single-domain circular Nanomagnets”. In: Phys. Rev. Lett. 83 (1999).
[16] Gia-Wei Chern, Muir J. Morrison, and Cristiano Nisoli. “Degeneracy and Criticality from Emergent Frustration in Artificial Spin Ice”. In: Phys. Rev. Lett. 111 (Oct.
2013).
[17] Z. Budrikis et al. “Domain dynamics and fluctuations in artificial square ice at finite temperatures”. In: New Journal of Physics 14 (2012).
[18] A. Farhan et al. “Direct observation of thermal relaxation in artificial spin ice”. In:
Phys. Rev. Lett. 111 (Aug. 2013).
[19] S. Zhang et al. “Crystallites of magnetic charges in artificial spin ice”. In: Nature 500 (2013), pp. 553–557.
[20] Hiroshi Shinaoka and Yukitoshi Motome. “Loop algorithm for classical Heisenberg models with spin-ice type degeneracy”. In: Physical Review B 82 (2010).
