Dipolar interactions of trapped dipolar molecular ions
Anna Filipsson
Abstract
The topic of this thesis is the development of a numerical model of trapped dipolar molecular ions forming a Coulomb crystal, such as in a linear ion trap.
The model regards the system as an ensemble of classical dipoles which form a linear chain in one dimension, and a square lattice in two dimensions. The model was simulated using Metropolis Monte Carlo and the classical Heisenberg model with dipole-dipole interaction. The possibility of phase transitions in the system was investigated by studying the thermodynamic properties at different temperatures, finite system sizes, and boundary conditions.
Previous studies of isotropic classical Heisenberg models in 1D and 2D have shown the absence of of any phase transitions. The 1D dipolar system modeled in this thesis showed no evidence of a critical temperature related to a phase transition, but rather exhibited a slow transition between an ordered and a dis- ordered state for temperatures T > 0. Similar results were obtained for the 2D system, but in this case a method of finite size analysis suggested the existence of a critical temperature. In both systems the dipoles showed a preferential ordering along the z-axis, corresponding to the trap-axis of the ion trap, at low temperature
The results of the simulations are mostly inconclusive due to difficulties in properly evaluating the critical phenomena, and further improvements and extensions to the preliminary model are needed to get more reliable results.
f¨or simulering av f˚angade dipol¨ara molekyl¨ara joner i Coulombkristaller, vilka bland annat bildas i linj¨ara jonf¨allor.
Modellen betraktar systemet som en ensemble av klassiska dipoler, vilka bildar en linj¨ar kedja i en dimension och ett kvadratiskt system i tv˚a dimensio- ner. Modellen simulerades med Metropolis Monte Carlo och den klassiska Hei- senbergmodellen med dipol-dipolinteraktioner. M¨ojligheten till fas¨overg˚angar i systemet unders¨oktes genom att studera termodynamiska egenskaper vid olika temperaturer, finita storlekar p˚a systemet, och randvillkor.
Tidigare studier av isotropiska klassiska Heisenbergmodeller i 1D och 2D har visat p˚a avsaknaden av fas¨overg˚angar. Det endimensionella dipolsystemet som modellerats i denna avhandling visar ingen kritisk temperatur kopplad till en fas¨overg˚ang, utan en l˚angsam ¨overg˚ang mellan ett ordnat och oordnat sy- stem f¨or temperaturer T > 0. Liknande resultat erh¨olls f¨or det tv˚adimensionella systemet, men i detta fall visade en analys av de finita storlekarna p˚a syste- met m¨ojligheten till en kritisk temperatur. I b˚ada systemen visade dipolerna en f¨oredragen ordning l¨angs systemets z-axel, vilket motsvarar den l˚angtg˚aende axeln i jonf¨allan, vid l˚aga temperaturer.
Resultatet fr˚an simuleringarna ¨ar otydliga, till stor del p˚a grund av sv˚arighet- er att utv¨ardera de kritiska fenomenen, och f¨orb¨attringar och ut¨okningar av den prelimin¨ara modellen beh¨ovs f¨or mer p˚alitliga resultat.
Anna Filipsson January 21, 2017
1 Introduction 3
1.1 Ion traps and Coulomb crystals . . . 3
1.2 Thesis outline . . . 4
2 Theory 5 2.1 Coulomb crystal properties . . . 5
2.2 Dipole interaction . . . 5
2.3 Magnetic models . . . 6
2.3.1 Ising model . . . 6
2.3.2 Heisenberg model . . . 7
2.4 Critical phenomena and phase transitions . . . 8
2.4.1 The partition function . . . 8
2.4.2 Thermodynamic quantities . . . 9
2.4.2.1 Energy . . . 9
2.4.2.2 Magnetization . . . 10
2.4.2.3 Heat capacity and susceptibility . . . 10
2.4.3 Phase transitions . . . 11
2.4.4 Spin-spin correlation function . . . 11
2.4.5 Critical behavior . . . 11
2.4.6 Finite size effects and the Binder cumulant . . . 12
2.5 Monte Carlo simulations . . . 12
2.5.1 Detailed balance . . . 12
2.5.2 Ergodicity . . . 13
2.5.3 Markov chains . . . 13
2.5.4 Metropolis algorithm . . . 13
2.5.5 Estimation of expectation values . . . 14
3 Computational model 15 3.1 Heisenberg model of dipolar interactions . . . 15
3.2 Initial configuration . . . 16
3.2.1 Thermalization . . . 16
3.3 Generating spin orientations . . . 16
3.4 Boundary conditions . . . 18
3.5 Thermodynamic quantities . . . 19
Anna Filipsson January 21, 2017
3.5.1 Energy . . . 19
3.5.2 Staggered magnetization . . . 19
3.5.3 Heat capacity and susceptibility . . . 19
3.5.4 Spin-spin correlation . . . 20
3.5.5 Binder cumulant . . . 20
3.6 Metropolis algorithm . . . 20
3.7 Random number generation . . . 21
3.8 Parallelization . . . 21
4 Results and discussion 23 4.1 The 1D system . . . 23
4.1.1 Thermodynamic quantities . . . 24
4.1.2 Correlation function . . . 25
4.1.3 Binder cumulant . . . 28
4.2 The 2D system . . . 28
4.2.1 Thermodynamic quantities . . . 29
4.2.2 Spin-spin correlation function . . . 30
4.2.3 Binder cumulant . . . 31
5 Summary and outlook 33
Chapter 1
Introduction
1.1 Ion traps and Coulomb crystals
The concept of regularly structured, spatially localized crystals of charged par- ticles was first introduced by Wigner in 1934 [1]. He predicted that given low enough electron densities, a condition at which the Coulomb energy exceeds the kinetic energy, the electrons in an electron gas would localize in space and form a crystalline structure. In order to minimize the potential energy in such a structure, the electrons form a string (1D), a triangular lattice (2D), or a body centered cubic lattice (3D) [2].
Coulomb crystals are a class of ordered structures of charged particles, such as atomic or molecular ions, analogous to classical Wigner crystals. While Wigner had no way of realizing such a crystal (the first experimental Wigner crystal would be achieved in the late 1970’s as a monolayer of electrons on the surface of liquid helium [3]), the development of ion traps and laser cooling in the following years made laboratory production of Coulomb crystals possible [4].
Currently, Coulomb crystals of atomic ions are especially relevant in the fields of quantum information processing and quantum simulation [4], while the properties of Coulomb crystals of molecular ions make them ideal for studies of chemical reactions at ultralow collision energies, so-called cold chemistry, and high-resolution precision spectroscopy [5].
The formation of translationally cold molecular ion Coulomb crystals can be done by trapping the ions using a linear-quadrupole trap, known as a linear Paul trap. The ions are symmetrically and dynamically trapped in the xy-plane perpendicular to the trap axis (z-axis), and confined along the trap axis by a static trapping potential. The structure of the resulting Coulomb crystal can differ from the structures of a Wigner crystal, depending on the trap parame- ters [6], but a strong trapping in the xy-plane will generally create a linear chain of ions along the trap axis.
1.2 Thesis outline
While many studies have been done on the topic of atomic ions in trapped Coulomb crystals, theoretical studies of similar systems of molecular ions have not received as much attention [4, 5]. The possibilities for future experiments and applications make Coulomb crystals of molecular ions an interesting field of study.
This thesis aims to study the dipolar interactions in an ensemble of dipolar molecular ions in a trapped Coulomb crystal, in particular the effect on the rotations of the ions due to dipole-dipole interactions, in order to investigate possible critical phenomena or phase transitions. To accomplish this, a classical model of the system is created and numerically simulated using Monte Carlo methods.
Chapter 2 will present the necessary theory behind the model, the methods of studying critical phenomena and phase transitions in statistical mechanics, and will give a general outline of the chosen Monte Carlo simulation technique chosen for the problem. The computational model developed to study the system will be explained in Chapter 3, and the results of the numerical simulations will be presented and discussed in Chapter 4. Chapter 5 will give a summary of the thesis and the obtained results, as well as discuss possible improvements to the model.
Chapter 2
Theory
2.1 Coulomb crystal properties
The properties and structure of a Coulomb crystal were briefly mentioned in the introduction. For a crystal formed in a linear Paul trap, the shape depends on the confining potential:
The one-dimensional crystal structure for strong enough trapping in the xy- plane is an equidistant string of ions [4].
For two-dimensional systems, the crystals have been observed to form a combination of triangular and circular shell configurations, with the inside of a large crystal being an almost perfect triangular lattice and the outside displaying a circular shell structure [6].
2.2 Dipole interaction
The interaction between two molecules with permanent dipole moments µi and µj is illustrated in figure 2.1.
The potential energy Epot can be found from multipole expansion [7], and is given by
Epot(µi, µj, R) = −µi· E(µj)
= − 1
4π0R3[3µiµjcos(ϑi) cos(ϑj) − µi· µj]
= − µiµj
4π0R3[2 cos(ϑi) cos(ϑj) − sin(ϑi) sin(ϑj) cos(ϕi− ϕj)], (2.1) where E(µj) is the electric field due to µj at the point of µi, R is the distance between the dipoles, and 0is the dielectric constant. The angle ϑ is taken from the z-axis and the angle ϕ is taken from the x-axis.
The minimum energy is obtained when the two dipoles are pointing in the same direction along the z-axis, while the maximum energy is obtained with the two dipoles pointing in opposite directions along the z-axis.
R
x
ϑ
j
ϑi
φj
φi
μi
μ
Figure 2.1: Permanent dipoles.
2.3 Magnetic models
Materials which exhibit spontaneous magnetization below a critical temperature Tc, even in the absence of any external fields, are commonly referred to as fer- romagnetic. A ferromagnetic system of dipoles consists of aligned dipoles which all contribute positively to the net magnetization of the system. In contrast, the adjacent dipoles in an antiferromagnetic system are instead anti-aligned causing the net magnetization to be zero. Similar to the antiferromagnet, a ferrimagnet has anti-aligned dipoles of unequal magnetic moment, causing the material to have a non-zero net magnetization
An important tool in statistical physics are the various magnetic lattice mod- els used to explain magnetic behavior in materials. These models are simplified representation of magnetic systems wherein a finite number of magnetic dipole moments are arranged in a lattice and each moment is allowed to interact with its neighbors.
Although the system considered in this thesis is one of classical dipoles, the terminology associated with these magnetic models will be adopted for simplic- ity and because of their predominance.
2.3.1 Ising model
The Ising model [8] represents the simplest and perhaps most well-known mag- netic model. In this model each spin (dipole) is given the discrete value ±1. A ferromagnetic system consists of either all positive or all negative spins of mag- nitude 1. The spins in an antiferromagnetic system instead alternate between +1 and −1.
Anna Filipsson January 21, 2017
The energy of the Ising model is given by the Hamiltonian H = −Jex
X
i,j
σiσj− HX
i
σ2i (2.2)
where the first term is the exchange interaction, with Jexbeing the site-indepen- dent interaction constant, and σi = | ± 1| being the spin at lattice site i. The second term is the energy due to an external field H. For a ferromagnet, Jex> 0, while for an antiferromagnet Jex< 0.
The exact solution for the 1D Ising model was provided in the original article, and shows the absence of a phase transition except at temperature T = 0. The exact solution of the 2D zero-field Ising model on a square lattice does however show the existence of a phase transition [9].
2.3.2 Heisenberg model
The interacting dipoles in section 2.2 are not restricted to point in discrete direction (up or down), instead they are permitted to point anywhere on the unit sphere. In simulating such a system, the simple Ising model proves inadequate and we instead turn to one of the continuous magnetic models.
In the Heisenberg model [10] the unit length spin vectors have components in the x, y, and z-directions, the spin at lattice site i is thus given by Si = (Six, Siy, Siz), where |Si| = 1.
Typically, the Hamiltonian used in the Heisenberg model consist of a number of one- and two-spin terms, e.g.
H = Hex+ HH+ HD+ Hdip, (2.3) where the different terms are explained below.
Exchange interactions
The Heisenberg exchange interaction arises from the Coulomb repulsion of elec- trons and the Pauli exclusion principle. The interaction energy between nearest neighbor pairs hi, ji of spins is given by
Hex= −Jex
X
hi,ji
Si· Sj, (2.4)
where the three-dimensional spin vectors Si are coupled by the interaction con- stant Jex. Jex> 0 leads to ferromagnetic ordering and Jex< 0 to antiferromag- netic ordering, similar to the Ising model. This interaction is isotropic, meaning that in a ferromagnetic material the spins tend to align in parallel, but do so with no preferential direction of alignment.
External field
In the presence of an external field the Heisenberg Hamiltonian contains the term
HH= −H ·X
i
Si, (2.5)
where the field H is coupled to single spins.
Single-ion anisotropy
The single-ion anisotropy, also called the magnetocrystalline anisotropy, is an- other short-range interaction which arises from coupling between the spin-orbit of the ion and the surrounding crystalline lattice. For a system with a preferred orientation along the z-axis, the energy contribution is given by
HD= −DX
i
Si,z2 , (2.6)
where D is the anisotropy constant.
Dipole interaction
The dipole interaction arises from the long-range interaction between magnetic moments in the system. It is analogous to the interaction energy of the two classical dipoles discussed in section 2.2. The Hamiltonian for the dipolar inter- action between spins is given by
Hdip= −JdipX
i6=j
"
3(Rij· Si)(Rij· Sj)
R5ij −Si· Sj
R3ij
#
, (2.7)
where |Rij| = Rij is the distance vector between the spins Si and Sj, and Jdip
is the dipolar interaction constant.
The exact solution for the classical nearest-neighbor isotropic Heisenberg model in one dimension, with H = Hex+ HD, was presented by Fisher in 1964 [11].
The article derives thermodynamic properties for an open chain of (N +1) spins.
In another article [12], the results are extended to a closed ring of N spins in the thermodynamic limit N → ∞. No exact solutions exist for the two- and three dimensional Heisenberg models [13].
Studies of the isotropic 1D and 2D Heisenberg model (systems with only exchange interaction) have shown the absence of any phase transition at non- zero temperatures [14]. However, studies of two-dimensional systems where the Hamiltonian includes anisotropic terms, such as the single-ion anisotropy and dipolar interactions, show the existence of long-range magnetic ordering. In particular, the preferred direction of spins in the presence of long-range dipole interaction induced anisotropy is in-plane [15].
2.4 Critical phenomena and phase transitions
2.4.1 The partition function
In order to study the thermodynamic behavior of a canonical ensemble, a system in thermal equilibrium with a heat bath at a constant temperature T , one makes
Anna Filipsson January 21, 2017
use of the partition function Z. For a classical system this has the form
Z = X
all states
e−H(s)/kBT, (2.8)
where H(s) is the Hamiltonian of the system in state s and kB is Boltzmann’s constant. The partition function is a sum over all possible states of the system, and with it one can calculate the thermodynamic quantities of the system.
The probability that the system occupies a particular state s is given by the Boltzmann distribution
P (s) = 1
Ze−H(s)/kBT. (2.9)
Borrowing terminology from quantum physics, the expectation value of an ob- servable A in the system at a certain temperature and state can be calculated as
hAi =X
s
A(s)P (s) = 1 Z
X
s
A(s)e−H(s)/kBT, (2.10) this calculation can be regarded as a time average over many measurement of A in the system, and is an approximation of the experimental average.
The mean square deviation of the observable A is
h(A − hAi)2i = hA2i − hAi2. (2.11) The square root of equation 2.11 represents the sample standard deviation, and gives a measure of the variation over time of A and the difference between the predicted and observed values.
2.4.2 Thermodynamic quantities
The thermodynamic quantities commonly used in studying magnetic systems can be calculated directly from the partition function 2.10 [16]. The following sections list the quantities of particular interest.
2.4.2.1 Energy
The expectation value of the internal energy is given by hEi = 1
Z X
s
H(s)e−H(s)β, (2.12)
where β = 1/kBT .
In terms of the derivative of the partition function the energy can be calcu- lated as
hEi = −∂ log Z
∂β . (2.13)
By taking the derivative of log Z a second time, the squared energy can be calculated as
hE2i = 1 Z
X
s
H2(s)e−H(s)β =∂2log Z
∂β2 . (2.14)
2.4.2.2 Magnetization
The expectation value of the magnetization is given by hM i = 1
Z X
s
M (s)e−H(s)β. (2.15)
The magnetization M of a system can be regarded as the conjugate variable to an applied magnetic field H and is a response to the variations of this field.
In the Hamiltonian the effect of an applied field appears as −M H, and the magnetization can thus be calculated in terms of the logarithmic derivative of the partition function with respect to H
hM i = 1 β
∂ log Z
∂H , (2.16)
and by taking the derivative of log Z a second time, the squared magnetization can be calculated as
hM2i = 1 Z
X
s
M2(s)e−H(s)β = 1 β
∂2log Z
∂H2 . (2.17)
When using one of the magnetic lattice models, the magnetization can be evaluated from the average over the states s. For the Heisenberg model with S(Sx, Sy, Sz) the magnetization is given by
hM i = hX
s
S(s)i. (2.18)
2.4.2.3 Heat capacity and susceptibility
Fluctuations of observable quantities describe their variation over time and can provide additional important information of the system properties near a phase transition. These are calculated using the mean square deviation (equation 2.11) of the observable.
The heat capacity CV is calculated from the derivative of the internal energy CV = − 1
kBT2
∂E
∂β, (2.19)
or, by combination of equations 2.13-2.14, from the energy fluctuation relation CV = hE2i − hEi2
kBT2 . (2.20)
The magnetic susceptibility χ is calculated from the fluctuation of the mag- netization
χ = hM2i − hM i2
kBT . (2.21)
Anna Filipsson January 21, 2017
2.4.3 Phase transitions
To describe a system in thermal equilibrium which transitions between a dis- ordered and ordered phase, one considers an order parameter: a macroscopic quantity with the property that its value changes from zero to non-zero when undergoing a phase transition [17].
Depending on the kind of physical system studied, the order parameter is defined differently. For instance, the order parameter for a ferromagnetic system is the spontaneous magnetization. Below a critical temperature the system is in its ordered phase and exhibits a non-zero spontaneous magnetization, while above it the system is disordered and the magnetization is identically zero.
For the antiferromagnetic system, the order parameter is instead the stag- gered magnetization [18].
The phase change between paramagnetism and ferromagnetism occurs at the Curie temperature of the material, and the system is only ferromagnetic below this critical temperature. Likewise, the change between paramagnetism and antiferromagnetism occurs at the corresponding N´eel temperature.
If the order parameter exhibits a discontinuity in the first derivative of the free energy at the critical temperature, and thus ”jumps” from one value to another, this is known as a first-order phase transition. In contrast, the fer- romagnetic or antiferromagnetic phase transition in magnetic systems such as the Ising or Heisenberg model, discussed in section 2.3, are examples of con- tinuous phase transition, in which the order parameter changes continuously.
In this type of transition the first derivative of the free energy is continuous at the critical temperature, while the second derivative of the free energy displays singularities.
2.4.4 Spin-spin correlation function
The spin-spin correlation function is a measure of how correlated a spin is with other spins in the system, and is given by
G(R) = hSi· Si+Ri, (2.22)
where R is the distance between the two spins. If the spins are completely correlated G(R) = 1, while if no correlation exists G(R) = 0. The correlation is largest for the nearest neighbors of Si.
In a completely disordered system one expects to find no correlation at all. In general, some small ordering is present in systems such as the Ising or Heisenberg model for finite temperatures, giving rise to a non-zero correlation [13].
The correlation function follows an exponential decay G(R) ∝ e−R/ξ. The correlation length, ξ, is a measure of the typical size of ordered clusters in the system [16].
2.4.5 Critical behavior
By examining the thermodynamic quantities and their variations one finds in- formation regarding the location of a phase transition.
As mentioned in section 2.4.3, the order parameter, assumed here to be the magnetization, will go from non-zero to zero when undergoing a continuous phase transition. In this region, the fluctuations of both the energy and magne- tization will be large, causing the heat capacity CV and magnetic susceptibility χ to diverge at the critical temperature Tc. Similarly, the correlation length ξ divergence to infinity at Tc.
2.4.6 Finite size effects and the Binder cumulant
The theoretical infinite divergence mentioned in the previous section is a prop- erty of an infinite system. When modeling a system using one of the magnetic lattice models, the finite size of the lattice introduces effects which smooth out transitions and makes it difficult to determine critical behavior.
Methods exist for obtaining the properties of an infinite system by observing how the thermodynamic quantities vary with different finite lattice sizes L, this is known as finite size scaling [13].
Additional information on the critical behavior of a system can be found by considering the fourth order Binder cumulant of the order parameter [19], given by
UL= 1 −1 3
hM4i
hM2i2. (2.23)
At large enough lattice sizes, curves of ULas a function of temperature intersect at a fixed point located at the critical temperature Tc.
2.5 Monte Carlo simulations
Monte Carlo (MC) simulations are a collection of numerical tools for studying systems and problems which are difficult to deal with analytically [13]. The interacting many-particle system with partition function and thermodynamic properties mentioned in section 2.4 is an example of such a system. By using MC simulations, one can determine the expectation values of observables (see equation 2.10) without the explicit calculation of the partition function. Instead, repeated random sampling is used to estimate results while taking account of statistical fluctuations.
The general idea behind so-called importance sampling MC simulations is to generate a set of possible states such that the probability of the system occupying a certain state equals the probability given by a chosen distribution, Pn. Two conditions are placed on the process to guarantee that the correct distribution is reached: it must fulfill detailed balance and ergodicity.
2.5.1 Detailed balance
The probability of a system being in a state νmis given by Pm, and the prob- ability that it will be in a state νn a time dt later is given by the transition
Anna Filipsson January 21, 2017
rate W (νm → νn). Detailed balance states that transition rate from state νm
to state νn must equal the transition rate from νn to νm, that is
PmW (νm→ νn) = PnW (νn→ νm). (2.24)
2.5.2 Ergodicity
A process is said to be ergodic if all possible configurations of the system are attainable with a non-zero probability. Thus, if an ergodic Monte Carlo simu- lation were to run for an infinite amount of time, all states would be ”visited”
by the process.
2.5.3 Markov chains
A Markov chains is a time-dependent stochastic process in which successive states are generated from their immediately preceding states, the process retains no memory of any other past states. If the process fulfills the requirements on detailed balance, the Markov chain is said to be reversible. The process is ergodic if there is a number N such that any state in the chain can be reached in N steps [13].
2.5.4 Metropolis algorithm
The Metropolis algorithm [20] is a reversible and ergodic process commonly used to simulate magnetic models. The method avoids the problem of an unknown partition function by generating a Markov chain of states.
The probability to be in state νn is given by the Boltzmann distribution (equation 2.9). The requirement on detailed balance thus tells us that the relative probability between two successive states in the Markov chain is
Pm
Pn = e−(E(νn)−E(νm))/kBT, (2.25) where νn is generated from νm.
Setting ∆E ≡ E(νn)−E(νm), the transition rate of the Metropolis algorithm is
W (νm→ νn) =
( e−∆E/kBT, if ∆E > 0 (2.26a)
1, if ∆E < 0. (2.26b)
If the new energy of the system E(νn) is lower than the previous energy E(νm), the transition to state νn is automatically accepted. If the new energy is higher, the transition is accepted with the probability given by equation 2.26a.
In the case of ∆E > 0, the algorithm looks as follows:
1. Start with an initial state, νm 2. Generate a new state, νn
3. Calculate the energy change ∆E in going from νm→ νn
4. Generate a random number r uniformly distributed in [0, 1]
5. Accept the transition νm→ νn if r < e−∆E/kBT
In simulating a magnetic model one performs steps 1-5 a suitable number of times before collecting information of the thermodynamic quantities of the system and adding them to a total statistical average. This corresponds to one Monte Carlo (MC) step.
2.5.5 Estimation of expectation values
The estimation of expectation values of observables in MC Metropolis simula- tions is done by averaging the observable (the collected total statistical average) over n MC steps [13]. The expectation value of the observable A is then calcu- lated as
hAi = 1 n
n
X
j=1
Aj. (2.27)
The expectation value of A to a power p can likewise be estimated by
hApi = 1 n
n
X
j=1
Apj. (2.28)
Chapter 3
Computational model
3.1 Heisenberg model of dipolar interactions
The system considered in this thesis is an ensemble of dipolar molecular ions in a trapped Coulomb crystal. As mentioned in section 2.1, the ions in such a structure are confined to lattice sites, making the magnetic lattice models discussed in the theory a good choice for simulating the system. The model has continuous degrees of freedom and is thus best simulated using the Heisenberg model.
By assuming a system consisting of single-species classical dipoles, with dipole moment µ = |µ| and potential energy given by equation 2.1, the Hamil- tonian of the system is
H = − µ2 4π0R3
X
hi,ji
[2 cos(ϑi) cos(ϑj) − sin(ϑi) sin(ϑj) cos(ϕi− ϕj)]. (3.1)
This is analogous to the Hamiltonian for Heisenberg dipolar interactions given by equation 2.7 when setting the interaction constant Jdip= µ2/4π0R3, from this point forward referred to simply as J . Because J > 0, the system is similar to an antiferromagnetic material.
Since J ∝ R−3, the dipole-dipole interactions are considered only for nearest neighbor pairs hi, ji, in contrast with the long-range dipole interaction Hamil- tonian presented in the theory.
In the simulations, J = 1. With this in mind, the inter-ion distance is set to be fixed and equal between all dipoles in the system, and the lattice spacing is thus R = 1 in units of [(µ2/4π0)−3].
For consistency, the dipoles in the model system will be referred to as ”spins”
in the following sections.
-1 1 0 3 2 4 -0.4
-0.2 0 0.2 0.4
6 5
(a) 1D system
-0.4 4
3 -0.2
0
0 0.2 0.4
2 1
1
2 0
3 -1
(b) 2D system
Figure 3.1: Initial antiferromagnetic spin configurations in 1D and 2D.
3.2 Initial configuration
The initial configuration is set to be antiferromagnetic, with spins pointing in alternating directions (up or down). The initial angles (see figure 2.1) are (ϑ = π/2, ϕ = 0) for spins pointing up, and (ϑ = π/2, ϕ = π) for spins pointing down.
In one dimension, the model is a linear chain of spins. The number of spins N in the system is equal to the ”length” L of the lattice. The initial configuration of a one dimensional system with N = 6 spins can be seen in figure 3.1a.
In two dimensions, the system is modeled on a square lattice. The number of spins in the system is N = L × L. The initial configuration of a two dimensional system with N = 16 spins can be seen in figure 3.1b.
3.2.1 Thermalization
After initialization, the system is thermalized for nthermMC steps. The thermal- ization steps are performed to ensure that the system approaches an equilibrium state, as well as the desired probability distribution. These steps are then dis- carded, as no recording of statistical data is performed while thermalizing the system.
3.3 Generating spin orientations
The simplest way to generate new spin orientations is to chose a random direc- tion from the total solid angle spectrum, using random spin angles 0 ≤ ϑ < π and 0 ≤ ϕ < 2π. In order for orientations to be homogeneously distributed on the unit sphere, the new angles are chosen as
ϑ = cos−1(2r1− 1) (3.2a)
ϕ = 2πr2, (3.2b)
where r1 and r2 are uniformly distributed random numbers, r ∈ [0, 1).
Anna Filipsson January 21, 2017
(a) The size of the cone is defined by the opening angle δ. The po- sition of the new spin Spis deter- mined by angles α and β.
(b) The coordinate system of Spis rotated such that the center of the cone coincides with the direction of the old spin S0.
Figure 3.2: Generation of new spin direction within a spherical cone about the old spin direction.
Efficiency problems are likely to arise when generating new spin directions from the whole unit sphere. Large changes in direction have a high probability of rejection, especially at low temperatures, since the resulting energy change ∆E is large compared to kBT . One solution to this problem is to instead generate new directions within a cone about the previous direction [18, 21].
The cone generation method used in this study is illustrated in figure 3.2 and works according to the following steps:
1. Define a maximum cone opening angle δ, which restricts the size of the cone.
2. Generate a new direction Sp within the spherical cone with angles 0 < α < δ from equation 3.2a, and 0 < β < 2π from equation 3.2b.
3. Rotate the Sp-coordinate system so that the ˆx-axis coincides with the direction of S0.
4. The rotated Spbecomes the new spin direction Sn with angles ϑ0and ϕ0.
When the cone angle δ is too small, a high number of the suggested direction changes are accepted and the evolution of the system will become too slow. In order to avoid this, the simulations use a method of adjusting the cone angle according to the acceptance rate of the previous MC step. The acceptance rate of a step is defined as
a = number of accepted direction changes total number of suggested direction changes,
where the total number of suggested direction changes in one MC step is equal to the number of spins of the system, N . If a > 0.5, δ is increased by about 2% for the next MC step, if a < 0.5, δ is decreased by about 2% for the next MC step. Hence, the percentage of accepted changes is kept around 50%. The initial size of the opening angle is chosen to be δ = 0.2π.
By running simulations using the different methods of spin orientation (choos- ing α, β on the whole unit sphere or on a spherical cone with or without adjust- ment of δ) the methods were determined to produce very similar results. The cone method with adjustment of δ was the only generation technique which gave acceptance rates a > 5% at the lowest temperatures and was therefore the chosen method for the main simulations.
3.4 Boundary conditions
For systems with free boundary conditions, the spins on the edges of the lattice will interact with fewer neighbors than those far from the boundaries. These edge spins will introduce boundary effects to the system [13].
The boundary effects of a finite system can be eliminated by using periodic boundary conditions in which the lattice sites on one edge of the system are considered neighbors to the lattice sites on the opposite edge. An example of a 2D system with periodic boundary conditions is shown in figure 3.3. Care must be taken when simulating an antiferromagnetic system with periodic boundary conditions, since odd numbered lattice side lengths will break the anti-alignment at the ”seam” of two edges. Thus, all system simulations presented in this thesis have even lattice side lengths (L = 4, 8, 10, etc). Finite size effects will still be present in a system with periodic boundary conditions.
7 8 9
1
4
7
1 2 3
3
6
9
1 2
4 5
3
6
7 8 9
Figure 3.3: Illustration of periodic boundary conditions in 2D system. Lattice sites considered to be neighbors of sites on the system boundary are marked in gray.
Anna Filipsson January 21, 2017
3.5 Thermodynamic quantities
The thermodynamic quantities presented in section 2.4.2 are used to determine the properties of the simulated systems.
The problem has been scaled in terms of the dimensionless parameter (re- duced temperature) T∗= kBT /J . For simplicity, the notation change T∗→ T is made for the remainder of the thesis unless otherwise stated.
3.5.1 Energy
The average energy hEi is calculated using the method in section 2.5.5. The average energy per spin for a specific temperature is calculated as
hei = 1
NhEi, (3.3)
The squared energy hE2i is collected using E in equation 2.28.
3.5.2 Staggered magnetization
The staggered magnetization in the k-direction is calculated as
Mk=
N
X
i
(−1)iSi,k, (3.4)
where i are sites on the lattice. The magnitude of the total magnetization is M =
q
Mx2+ My2+ Mz2. (3.5) The average of the total magnetization, hM i, is calculated using equation 2.18 and the method in section 2.5.5. The average magnetization per spin for a specific temperature is then calculated as
hmi = 1
NhM i. (3.6)
The squared magnetization, hM2i, and hM4i are calculated using equa- tion 2.28.
3.5.3 Heat capacity and susceptibility
The heat capacity per spin is calculated from the fluctuation of the energy CV = Nhe2i − hei2
T2 , (3.7)
and the susceptibility per spin is calculated from the fluctuation of the total magnetization
χ = Nhm2i − hmi2
T . (3.8)
3.5.4 Spin-spin correlation
The spin-spin correlation function G(x) is calculated using equation 2.22, with i = 0 and x corresponding to the discrete lattice sites 0 ≤ x < L. The function becomes
G(x) = hS0· Sxi. (3.9)
In the 2D system, the correlation function is calculated along one row in the center of the lattice, i.e. at L/2.
The whole lattice length cannot be considered when extracting the correla- tion length ξ from the correlation function in systems with periodic boundary condition. The maximum correlation length is then bounded by L/2.
3.5.5 Binder cumulant
The fourth order Binder cumulant associated with the order parameter M is calculated using equation 2.23.
As mentioned in section 3.4, the finite size effects are present in systems with both free and periodic boundary conditions. The cumulant is thus calculated for both these cases.
3.6 Metropolis algorithm
A general form of the metropolis algorithm was mention in section 2.5.4. A de- tailed description of the simulation specific Metropolis algorithm is given below:
For each temperature T , perform n MC steps given by the following scheme:
1. Start with an initial spin configuration (according to section 3.2) 2. Pick a random lattice site, 0 ≤ i < N , containing spin Si(ϑ, ϕ) 3. Generate a new trial spin, S0i(ϑ0, ϕ0), using the method described in section 3.3
4. Calculate the energy change ∆E = E(Si0) − E(Si) using the energy in equation 3.1
(a) If ∆E < 0, accept the direction change and set Si ← Si0
(b) Otherwise, generate a random number r uniformly distributed in [0, 1) and accept the direction change if r < e−∆E/T
(c) Else, reject the direction change 5. Repeat steps 2-4 N times
After each MC step, the values of E, E2, M , M2, and M4 are collected and added to their respective total average. After n MC steps, the estimations of the expectation values, together with CV, χ, UL, and G(x), are calculated according to section 3.5.
Anna Filipsson January 21, 2017
3.7 Random number generation
The importance of random number generation in the Metropolis algorithm can- not be understated, as low quality random numbers can introduce correlations to the generated states and systematic errors to the system [13]. A common choice of pseudo-random number generator (PRNG) is the Mersenne Twister algorithm [22]. The period of this algorithm, the number of output bits after which the generator starts to repeat itself, is 219937− 1,
A faster version of the Mersenne Twister, and the PRNG chosen for use in this thesis, is the SIMD-oriented Fast Mersenne Twister (SFMT) algorithm [23].
The algorithm has support for various periods, with a maximum of 2216091− 1, but for the purpose of this thesis has been set to use the same period as the Mersenne Twister.
3.8 Parallelization
The C-program written to run the simulations was parallelized using Open- MPI [24] in order to speed up the simulations and provide better statistical results.
In the one-dimensional system, the total number of temperatures Tj = j∆T , where ∆T is the size of the temperature step, are evenly divided among the avail- able processors. To make sure that the total simulation time of all processors are approximately equal, the mth processor (among a total of mtot) will calculate the temperatures at j satisfying
j mod mtot= m.
In the two-dimensional system, each available processor runs its own in- stance of the simulation for all temperatures T using different initial states of the PRNG and thus different initial spin configurations after thermaliza- tion. Individual processes calculate ”local” MC averages of the thermodynamic quantities. At the end of a temperature simulation, the local averages are then averaged between all processors to calculate a ”global” average. In doing this, it’s expected that the number of MC steps run by each processor can be reduced while still getting statistically reliable results.
Chapter 4
Results and discussion
4.1 The 1D system
The one-dimensional system was simulated for four different system sizes with free and periodic boundary conditions, each for 40 temperature steps between T = 0.001 and T = 2.0, using ntherm= 104thermalization steps and nMC= 106 MC steps.
Below a temperature of T ≈ 0.5 the spins start to align antiferromagneti- cally along the z-axis of the system. The difference between spin orientations at temperatures above or below this point can be seen in figure 4.1, where fig- ure 4.1a shows the final spin orientations at T ≈ 0.975 and figure 4.1b shows the final spin orientations at T = 0.001.
10 8
6 4
-0.4 2 -0.20.20.40
0 0
(a) Temperature T = 0.975
10
8
6
4
2
0 0
(b) Temperature T = 0.001 Figure 4.1: Spin orientations in 1D system of N = 10 spins, after 104 thermal- ization steps and 106 MC steps.
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 -2
-1.8 -1.6 -1.4 -1.2 -1 -0.8 -0.6 -0.4 -0.2
(a) Energy per spin
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
(b) Magnetization
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
(c) Heat capacity
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5
(d) Susceptibility
Figure 4.2: Thermodynamic quantities and fluctuations of the 1D system with free boundary conditions, for system sizes N = 10, 16, 20, and 32.
4.1.1 Thermodynamic quantities
The energy per spin and staggered magnetization, as well as the corresponding fluctuations (heat capacity and susceptibility) for the different system sizes can be seen in figure 4.2 for a system with free boundary conditions, and in figure 4.3 for a system with periodic boundary conditions.
At T ≈ 0, all systems show a total magnetization M ≈ 1, suggesting an almost perfectly ordered system, the spin orientation in figure 4.1b tells us that this ordering occurs for the z-component of the magnetization, Mz. For tem- peratures T > 0, the systems slowly move to a disordered state. The transition seem to occur over a wide temperature range approximately centered about T = 0.55 for free boundary conditions, and T = 0.65 for periodic boundary conditions.
Anna Filipsson January 21, 2017
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
-2 -1.8 -1.6 -1.4 -1.2 -1 -0.8 -0.6 -0.4 -0.2
(a) Energy per spin
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
(b) Magnetization
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
0 0.5 1 1.5 2 2.5
(c) Heat capacity
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5
(d) Susceptibility
Figure 4.3: Thermodynamic quantities and fluctuations of the 1D system with periodic boundary conditions, for system sizes N = 10, 16, 20, and 32.
4.1.2 Correlation function
The spin-spin correlation function, and the attempted exponential fit, of a sys- tem of N = 10 spins and free boundary conditions is shown for two different temperatures in figure 4.4.
The correlation between spin S0and anti-aligned spins (i.e. at i = 1, 3, 5 . . .) exhibits an oscillating behavior. In particular, the correlation between spins S0 and S1 deviates from the exponential decay, causing the fit to fail and thus rendering the study of correlation length ineffective. This pattern is present for all system sizes and boundary conditions and seems to be a phenomenon related to the dipolar interactions, which will be further commented on in the discussion of the correlation function of the 2D system.
The same results are obtained when starting with a spin Si at a lattice site away from the boundary, i.e. at i > 0, suggesting this deviation is not caused by boundary effects.
In figures 4.5 and 4.6 are shown the correlation function of the x- and z- components of the spins, hS0,xSi,xi and hS0,zSi,zi, respectively. Figure 4.5 shows
0 1 2 3 4 5 6 7 8 9 10 0.9986
0.9988 0.999 0.9992 0.9994 0.9996 0.9998 1
data fitted curve
(a) T = 0.001
0 1 2 3 4 5 6 7 8 9 10
-0.2 0 0.2 0.4 0.6 0.8 1
data fitted curve
(b) T = 2.0
Figure 4.4: Spin-spin correlation function in the 1D system with N = 10 and free boundary conditions at different temperatures T .
0 1 2 3 4 5 6 7 8 9 10
-2 -1 0 1 2 3 4 5 6 10-4
data fitted curve
(a) x-components at T = 0.001
0 1 2 3 4 5 6 7 8 9 10
-0.05 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35
data fitted curve
(b) x-components at T = 2.0 Figure 4.5: Spin-spin correlation function of x-components in the 1D system with N = 10 and free boundary conditions at different temperatures T .
that the x-components display a similar shape to that seen in figure 4.4. The same shape is present in the correlation function of the y-components. Fig- ure 4.6b shows the z-component at a high temperature. The function shows an exponential decay which slows down for decreasing temperatures, until at approximately T = 0.5 where the boundary spins start to deviate from the ex- ponential function. The lower correlation of the boundary spins are believed to be due to boundary effect causing these spins to be slightly less aligned with the z-axis of the system compared to the non-boundary spins.
Anna Filipsson January 21, 2017
0 2 4 6 8 10
0.9988 0.99885 0.9989 0.99895 0.999 0.99905 0.9991 0.99915
data fitted curve
(a) z-components at T = 0.001
0 1 2 3 4 5 6 7 8 9 10
-0.05 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35
data fitted curve
(b) z-components at T = 2.0
0 1 2 3 4 5 6 7 8 9 10
0.25 0.3 0.35 0.4 0.45 0.5 0.55
data fitted curve
(c) z-components at T = 0.5
0 1 2 3 4 5 6 7 8 9 10
0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45
data fitted curve
(d) z-components at T = 0.8 Figure 4.6: Spin-spin correlation function of z-components in the 1D system with N = 10 and free boundary conditions at different temperatures T .
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
0.25 0.3 0.35 0.4 0.45 0.5 0.55 0.6 0.65 0.7
(a) Free boundary conditions
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.55 0.6 0.65 0.7
(b) Periodic boundary conditions Figure 4.7: Binder cumulant of the different 1D system sizes.
4.1.3 Binder cumulant
The fourth order Binder cumulant of the order parameter for the different sys- tem sizes and boundary conditions are shown in figure 4.7. Neither boundary condition show any visible crossing point of the different system sizes. This ap- parent absence of a phase transitions at temperatures T > 0 agrees with previous results of the corresponding 1D Heisenberg model with only exchange interac- tion energy [14], suggesting that the behavior of the purely dipolar Heisenberg model is similar to that of the isotropic system.
4.2 The 2D system
The two-dimensional system was simulated for four different system sizes with free and periodic boundary conditions, each for 40 temperature steps between T = 0.001 and T = 3.0, using ntherm= 103thermalization steps and nMC= 105 MC steps.
The final spin configuration at T = 0.001 of the free boundary condition system of N = 16 spins can be seen in figure 4.8. Similar to the 1D system, the spins are aligned in an antiferromagnetic manner along the z-axis. The spins can also be seen to anti-align in the y-direction.
4 0 3
2 1
1
2 0
3 -1
Figure 4.8: Spin orientations in 2D system of N = 16 spins, at T = 0.001, after 103 thermalization steps and 105 MC steps.
Anna Filipsson January 21, 2017
0 0.5 1 1.5 2 2.5 3
-4 -3.5 -3 -2.5 -2 -1.5 -1 -0.5 0
(a) Energy per spin
0 0.5 1 1.5 2 2.5 3
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
(b) Magnetization
0 0.5 1 1.5 2 2.5 3
0 0.5 1 1.5 2 2.5 3
(c) Heat capacity
0 0.5 1 1.5 2 2.5 3
0 0.5 1 1.5 2 2.5 3 3.5
(d) Susceptibility
Figure 4.9: Thermodynamic quantities and fluctuations of the 2D system with free boundary conditions, for system sizes N = 16, 64, 100, and 256.
4.2.1 Thermodynamic quantities
The results of the 2D system with free boundary conditions can be seen in figure 4.9 and with periodic boundary conditions in figure 4.10. The order parameter in figures 4.9b and 4.10b show a slow transition from an ordered system at T ≈ 0 to a disordered one at T > 0. This is also evident from the heat capacity and magnetic susceptibility. Peaks in these quantities are present near T ≈ 1.77 in both systems.
The thermodynamic quantities in the smallest system, N = 16, is seen to deviate from the otherwise similar quantities of the larger systems, especially in the system with free boundary conditions. This is believed to arise from the finite size effects of the system and from boundary effects, both of which become most prominent for small lattice sizes. The larger lattice sizes can be seen to produce more similar results in terms of energy and magnetization.
The heat capacity and susceptibility show some large deviations at low tem- peratures for systems with a number of spins N ≥ 100. The exact cause of
0 0.5 1 1.5 2 2.5 3 -4
-3.5 -3 -2.5 -2 -1.5 -1 -0.5 0
(a) Energy per spin
0 0.5 1 1.5 2 2.5 3
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
(b) Magnetization
0 0.5 1 1.5 2 2.5 3
0 0.5 1 1.5 2 2.5 3 3.5
(c) Heat capacity
0 0.5 1 1.5 2 2.5 3
0 0.5 1 1.5 2 2.5 3 3.5 4
(d) Susceptibility
Figure 4.10: Thermodynamic quantities and fluctuations of the 2D system with periodic boundary conditions, for system sizes N = 16, 64, 100, and 256.
this problem is unknown, but it is suspected to depend on certain initial states producing systems where spins become ”locked” in aligned positions, altering the staggered magnetization of the system. This phenomenon is shown in fig- ure 4.11. When ”locking” occurs, the system is split in two regions of anti- aligned spins. In figure 4.11, this corresponds to the region of spins at positions z ≤ 3 and the region of spins at positions z ≥ 4. The large energy change
∆E required for a spin to break out of a locked position causes those suggested states to be rejected.
4.2.2 Spin-spin correlation function
The spin-spin correlation function, and the attempted exponential fit, of a sys- tem of N = 100 spins with periodic boundary conditions is shown for T = 3.0 in figure 4.12a.
Similar to the 1D system, the correlation between spins S0 and S1 deviates from the exponential decay. For comparison, the correlation function of a system
Anna Filipsson January 21, 2017
with exchange interaction (H = Hex) is shown in figure 4.12b. This systems exhibits the expected exponential decay, and a calculation of the correlation length ξ shows it to grow larger for T → 0, as discussed in the theory.
Because of the deviating points of the spin-spin correlation function of a system with only dipolar interaction energy, the correlation length cannot be determined. Since no method was found for accurately calculating the correla- tion length of a system exhibiting this sort of oscillating spin-spin correlation function, the study of phase transitions using the correlation length has been ineffective in this project.
4.2.3 Binder cumulant
The fourth order Binder cumulant of the order parameter for the different system sizes and boundary conditions can be seen in figure 4.13. Unlike the 1D system, the 2D system shows a well defined crossing point for the larger systems (N ≥ 64) at the same temperature, T ≈ 1.77, as the peaks in CV and χ seen in figures 4.9 and 4.10.
For comparison, a simulation was run for a 2D system using the exchange in- teraction energy instead of the dipolar interaction energy. The Binder cumulant of the exchange interaction simulations for the different boundary conditions is shown in figures 4.13c and 4.13d. This system shows no crossing point, in agreement with theory [14].
0 1 2 3 4 5 6 7 8 9
0 2 4 6 8 10
Figure 4.11: ”Locked” spins (marked with red) in 2D system of N = 100 spins, at T = 0.001, after 103thermalization steps and 105MC steps.
0 5 10 15 -0.2
0 0.2 0.4 0.6 0.8 1
data fitted curve
(a) Dipole interactions, T = 3.0
0 5 10 15
-0.2 0 0.2 0.4 0.6 0.8 1 1.2
data fitted curve
(b) Exchange interactions, T = 3.0 Figure 4.12: Correlation function for 2D systems with N = 16 and periodic boundary conditions.
0 0.5 1 1.5 2 2.5 3
0.1 0.2 0.3 0.4 0.5 0.6 0.7
(a) Dipole interaction, FBC
0 0.5 1 1.5 2 2.5 3
0.1 0.2 0.3 0.4 0.5 0.6 0.7
(b) Dipole interaction, PBC
0 0.5 1 1.5 2 2.5 3
0.4 0.45 0.5 0.55 0.6 0.65 0.7
(c) Exchange interaction, FBC
0 0.5 1 1.5 2 2.5 3
0.4 0.45 0.5 0.55 0.6 0.65 0.7
(d) Exchange interaction, PBC Figure 4.13: Binder cumulant of the different 2D system sizes, using free bound- ary conditions (FBC) or periodic boundary conditions (PBC), and different interaction energies.
