Magnetic Bloch oscillations in CoNb

Magnetic Bloch oscillations in CoNb

O

Aleksander Seland

Master of Science Thesis Department of Physics

The Faculty of Mathematics and Natural Sciences UNIVERSITY OF OSLO

Abstract

The material CoNb2O6 has in recent times attracted much attention due to having rich physics in applied magnetic fields. This thesis explores whether or not it is possible to observe magnetic Bloch oscillations in CoNb2O6. We do so by first establishing theoretically CoNb2O6

as a candidate material for these oscillations. Afterwards the magnetic and thermodynamical properties are explored through Monte Carlo simulations, verifying the code as we go, before searching for parameters that make detecting magnetic Bloch oscillation probable. We end up with suggesting a temperature of T = 2.8K and two magnetic field strengths,B = 0.31T and B= 0.38T, giving rise to two detectable Bloch frequencies for future neutron diffraction experiments.

Takk alle sammen!

Dette er den klassiske takketalen for når man er ferdig med masteroppgaven, men noen mennesker fortjener virkelig en stor takk fra meg. Først og fremst, en stor takk til veilederen min, Olav F. Syljuåsen. Takk for alle svar, og for at du i det hele tatt klarer å henge når jeg prøver å forklare spørsmålene mine. Setter også stor pris på ordspillene (slemt felt teori f. eks) og hyggelige samtaler. Neste på lista forloveden min. Hadde det ikke vært for henne, så hadde jeg glemt hva det vil si å leve oppi alt sammen. Jeg hadde mest sannsynlig hatt 1 meter langt skjegg tatt fyr i sollys uten henne, tusen takk! Neste på lista er selvsagt far og mor. Moren min lærte meg å holde ut med hardt arbeid, og faren min tente fysikk interessen med lange turer under stjernehimmelen samt fortellinger om noen rare greier kalt atomer, takk skal dere ha! Til syvende og sist, men kanskje den aller største grunnen til at jeg endte opp med å studere fysikk, er fysikklæreren min på videregående. Etter de første timene med fysikk, syntes jeg faget virket kjedelig, men det klarte han å snu raskt på. Så her sitter jeg nå med en ferdigskrevet masteroppgave og skal snart ut i læreryrke som fysikklærer. Du har hatt en stor innvirkning på det, tusen takk!

Contents

1 Introduction - the physical problem 1

2 The material 1

2.1 Theory of magnetic Bloch oscillations in CoNb2O6 . . . 1 2.1.1 Background . . . 1 2.1.2 1D Hamiltonian and magnetic Bloch oscillations . . . 3 2.1.3 Interchain interactions and the cost of moving a domain wall 5 2.2 The magnetic properties of Co²⁺ . . . 7 2.3 Structure and magnetic phases of CoNb2O6 . . . 8 2.4 A summary of existing models for CoNb2O6 . . . 12

3 The simulation 15

3.1 Deciding on a model . . . 15 3.2 The algorithm . . . 17 3.2.1 A brief summary of the Metropolis Monte Carlo algorithm . . 17 3.2.2 Verifying the algorithm . . . 19 3.2.3 Looking for domain walls and the Wannier-Zeeman ladder. . . 26 3.3 Simulation results . . . 28

4 Conclusion 31

A Direct exchange 31

B Electric Bloch oscillations 34

B.1 The Bloch theorem . . . 34 B.1.1 Semiclassical approach . . . 34 B.1.2 Quantum mechanical approach . . . 35

1 Introduction - the physical problem

Seven years ago, Coldea et al. (2010) performed neutron diffraction experiments on a crystal of CoNb2O6 in a transverse magnetic field and found several close lying dispersive modes. These modes can be viewed as a ladder of energies called the Wannier-Zeeman ladder (WZL) and is the magnetic analog to the Wannier-Stark ladder (WSL) in electrical systems. What inspired this thesis was the fact that the WZL found in CoNb2O6 is very well defined for the lowest lying levels and resemble the conditions for electric Bloch Oscillations already observed by Feldmann et al.

(1992). Since Kyriakidis and Loss (1998) proposed the magnetic analog to electric Bloch oscillations, many scientists have been searching to find a material that has suitable conditions for observing these oscillations. Inspired by the apparent WZL of CoNb2O6, the current work hopes to find a regime of parameters that allows for the observation of magnetic Bloch oscillations in CoNb2O6. This will be achieved by finding a suitable temperature and longitudinal magnetic field such that the cost,

| E|, of moving a domain wall becomes as uniform as possible for all domain walls in the material. This implies a high enough temperature to create domain walls, and a strong enough magnetic field in order for the three dimensional neighborhood at each site to be as uniform as possible. If such a parameter set can be found, neutron diffraction experiments might be able to show a significant peak for energies equal to the Bloch frequency ~! = | E|.

2 The material

2.1 Theory of magnetic Bloch oscillations in CoNb2O6

2.1.1 Background

When studying the results from Bloch (1929), Zener (1934) found that the Schrödinger equation, when given a static electric potential in addition to the periodic poten- tial, had oscillatory solutions when not considering transitions to other bands (Zener tunneling). Even stranger, the amplitude of the oscillations decreased while the fre- quency increased with stronger electric fields. This contradicts classical thinking;

one would assume a particle to accelerate under a constant force, and there is no good classical reason for the force changing the amplitude as well as the frequency of

decades Nenciu (1991). The reason for much of the debate concerned the semiclassi- cal approach of Zener, and it was some decades later that Wannier (1960) solved the system with both the periodic and electrical potential quantum mechanically. He showed (appendix B) that the wave functions of the electrons become localized and that the energy spectrum became discrete within each band. This lifting of the de- generacy in the crystal due to the electric field was dubbed the Wannier-Stark ladder.

It was not until 30 years later that the experimental discovery of these oscillations was done by Feldmann et al. (1992) in semiconductor superlattices.

The discovery of electric Bloch oscillations brought with it speculations about whether or not magnetic Bloch oscillations could exist. Kyriakidis and Loss (1998) explored this idea for quantum spin chains and modeled the border between two domains (""| ##) of spin directions as particles called “solitons”. In the literature, these solitons also go by the names “domain walls” and “kinks”, from which we will use

“solitons” and “domain walls” interchangeably. Their study, taking into account only one soliton, showed that these solitons become localized when an external magnetic field longitudinal to the chain is applied, and the energy spectrum becomes discrete.

This new spectrum was dubbed the Wannier-Zeeman ladder (WZL) where each step on the ladder is determined from the number of spins opposing the magnetic field (the location of the soliton). They showed that when applying a constant magnetic field, these solitons undergo oscillations, changing the size of the magnetic domain.

The conditions for observing magnetic Bloch oscillations are:

1. The material must have a strong anisotropy 2. There should be no interband tunneling

3. There should be no inelastic scattering by phonons or other particle like phe- nomenons in the crystal

4. Elastic scattering from impurities in the crystal should not occur, or at best, be kept at a minimum.

Condition 1 tells us that there should be a preferred direction of spin through a dominant exchange interaction. If this condition isn’t met, the notion of magnetic domains and their borders becomes hard to define. Condition 2 puts restrictions on the magnetic field by not allowing the chain to becomes fully polarized as this simply removes the soliton. The third condition permits the destruction of the oscillatory

motion through collisions and the last condition implies that impurities may be a problem for large oscillation amplitudes. CoNb2O6 is one such material with a strong anisotropy, and will be explored in the next sections.

2.1.2 1D Hamiltonian and magnetic Bloch oscillations

In order to show that CoNb2O6 is a candidate material for such oscillations, let us consider a general spin-1/2 XXZ Heisenberg Hamiltonian in an applied longitudinal magnetic field hz = gµBB (as proposed by Cabrera et al.). To allow for a tidy notation, we shall only consider a single chain of length 2N + 1 and let each atom on the chain be indexed by i. The hamiltonian is

Hˆ = X

i

⇣ J_zSˆ_i^zSˆ_i+1^z J_xy⇣

Sˆ_i^xSˆ_i+1^x + ˆS_i^ySˆ_i^y⌘

J_zSˆ_i^zSˆ_i+2^z h_zSˆ_i^z⌘

(2.1) By introducing the ladder operators,Sˆ^± = ˆS^x±iSˆ^y, we may rewrite the hamiltonian as:

Hˆ = X

i

✓

J_zSˆ_i^zSˆ_i+1^z J_xy 2

⇣Sˆ_i⁺Sˆ_i+1+ ˆS_i Sˆ_i+1⁺ ⌘

J_z⁰Sˆ_i^zSˆ_i+2^z h_zSˆ_i^z

◆

(2.2) When constructing a basis for one such chain, we need to chose an approximation.

Kyriakidis and Loss (1998) looked at chains containing only one domain wall, and thus neglecting the J_xy term from (2.2) as it produces more domain walls. This approximation was pointed out by Shinkevich and Syljuåsen (2012) to be hardly realizable in materials under application of a longitudinal field. The magnetic field will try to align the spins with the field, and when considering ferromagnetically locked chains as in CoNb2O6, the domains created will give rise to two domain walls.

We will thus consider states |m, l, Qi such that

|m, l, 1i = · · · ""|{z}^m###

l

" · · · +

|m, l,1i = · · · ##|{z}^m"""

l

# · · · +

where m gives the position of the first domain wall between position m 1 and m

application of our hamiltonian (2.2), we see that theJ_xy term will produce additional domain walls in our system except for when i = m 1 and l = 1. If given such a state, it’s action is to move flipped spin as

J_xy 2

NX1

i= N

⇣Sˆ_i⁺Sˆ_i+1+ ˆS_i Sˆ_i+1⁺ ⌘

|· · · ##"## · · ·i = J_xy

2 (|· · · #"### · · ·i+|· · · ###"# · · ·i) This translation of the domain walls by one site does not change the domain size, and can be neglected as it does not contribute to magnetic Bloch oscillations. Let us now consider our chain of length 2N + 1, and look at the action of the different terms in the spin-1/2 Hamiltonian when given a two-soliton state with l 6= 0:

NX1

i= N

J_zSˆ_i^zSˆ_i+1^z |m, l, Qi = J_z|m, l, Qi (2.3)

NX2

i= N

J_z⁰Sˆ_i^zSˆ_i+2^z |m, l, Qi = J_z⁰ (2 _l1)|m, l, Qi (2.4)

h_z XN

i= N

Sˆ_i^z|m, l, Qi = Qlh_z|m, l, Qi (2.5) where the ground energy contribution (l = 0 state) has been set to zero when considering large chains such that boundary conditions are negligible. Let us now in the spirit of Kyriakidis and Loss (1998) define the hermitian operator Qˆ such that Qˆ|m, l, Qi = Q|m, l, Qi and consider l 2 such that our hamiltonian becomes

Hˆ = J_z 2J_z⁰ Qlhˆ _z (2.6)

We see that this hamiltonian does not have any dynamics to move the domain walls and thus cannot produce magnetic Bloch oscillation. To remedy this, we may intro- duce a transverse magnetic field along the x-direction with hamiltonian contribution

Hˆ_x = 1 2h_x

XN

i= N

⇣Sˆ_i⁺+ ˆS_i ⌘

(2.7)

Its action on the basis states, keeping the two domain wall approximation, is Hˆx|m, l, Qi= 1

2hx(|m 1, l+ 2, Qi+|m+ 1, l 2, Qi) (2.8) We see that this terms expands or shrinks the domain by one spin at each end, thus having the dynamics necessary for magnetic Bloch oscillations. We can conclude that the effective one dimensional hamiltonian for magnetic Bloch oscillations in CoNb2O6 is:

Hˆ_1D = J h_x 2

XN

i= N

⇣Sˆ_i⁺ + ˆS_i ⌘

Qlhˆ _z (2.9)

where J = J_z 2J_z⁰ is an effective exchange between two adjacent atoms. This shows that CoNb2O6 indeed is a candidate material for magnetic Bloch oscillations.

The fact that it needs an external transverse field to have the dynamics needed for magnetic Bloch oscillations, means we’re in control of both the amplitude and frequency of the oscillations (Kyriakidis and Loss, 1998).

2.1.3 Interchain interactions and the cost of moving a domain wall

When solitons propagate in a constant longitudinal magnetic field, they change the energy of the system by flipping spins. In order to observe them, we want each move to cost exactly the same energy such that the associated Bloch frequency,

~! = | E|, is a unique value so as not to have several Bloch frequencies in the material. The cost moving a domain wall (l ! l + 1) given Hˆ_1D is E = h_z, but as shown in figure 3 of section 2.3, there exists interchain interactions that add a dependence on neighboring chains to the cost. By including these interactions and using | S| = 1 for S = 1/2, we get

| E| = J₁X

hiji

( S_i)S_j +J₂X

hiki

( S_i)S_k +h_z S_i

= J₁X

hiji

S_j +J₂X

hiki

S_k+h_z (2.10)

where the spin being flipped is at sitei, and hiji,hikidenotes the nearest interchain

values of | E|, and thus Bloch frequencies, we expect to see in CoNb2O6 in various magnetic fields h_z.

The AF ground phase (table 1 in section 2.3), as shown in figure 1, has the following cost of moving a domain wall

| E|AF = J₁X

hiji

S_j +h_z =|±J₁+h_z| =|sign(h_z)J₁+h_z| (2.11) due to the J₂ coupled spins canceling out. The spins coupled by J₁ are aligned and sum up to 1. The ± sign was resolved by noting that J₁ < 0, and that J₁ opposes the alignment imposed by h_z implying they should have opposite signs. We may expect the cost of moving domain walls to be dominated by this value for low fields and low temperatures.

J₁ J₂ a

b

Figure 1: The ground state of CoNb2O6forT <1.97K. The red and blue points represent different directions of spin along thez-direction.

For high fields, the system moves into a saturated (fully polarized) paramagnetic phase P1. In this phase, the cost of moving a domain wall is

| E|P⁰ = J₁X

hiji

S_j +J₂X

hiki

S_k +h_z = |sign(hz) (J₁+ 2J₂) +h_z| (2.12) where no spins cancel out and the factor of 2 comes from four J₂ coupled spins con- tributing. It is expected that the simulations show a transition from being dominated by | E|AF to | E|P M as the field increases.

For intermediate fields, the system is in the SF1 phase shown in figure 2, and has two values of | E|

| E|SF1 =

(| E|P⁰

|h_z| (2.13)

J₁ J₂ a

b

Figure 2: The magnetically induced SF1 phase with propagation vector (0,¹₃,0).

2.2 The magnetic properties of Co²⁺

The magnetism in CoNb2O6 is due to the cobalt atoms. The two outer4selectrons in the cobalt atoms are given to a ligand atom through ionic bonding leaving us with 7 3delectrons in the outermost shell. As known from quantum theory, this corresponds to n = 3 and l = 2 leaving us with 5 orbital states |n, l, m_li = |3,2, m_li available to the electrons. When taking into account spins and the Pauli exclusion principle, we notice that two of the five orbital states become occupied by two electrons with opposite spins. This leaves 3 electrons, each with its separate quantum number m_l and total spin given by the Clebsch-Gordan tables as either 1/2 or 3/2. In order to determine which spin has the lowest energy, we need to consider the spin orbit coupling contribution to the hamiltonianHˆ_SO = ˆL·S. The coupling strength factorˆ in the case of Co²⁺ has been found to be negative (Low, 1958), and as such, the S = 3/2 is the ground state while S = 1/2 is the excited state. These states have a gap of S = 1, allowing us to model these atoms as effective spin-1/2s with a relatively high g-factor to preserve the magnetic moment µ = gµ_BS.

2.3 Structure and magnetic phases of CoNb2O6

CoNb2O6 has an orthorhombic crystal structure (space group P bcn) and consists of weakly coupled chains in the c direction on an isosceles triangular lattice (a-b plane) as shown in figure 3 (Cabrera et al., 2014a; Kobayashi, Mitsuda, et al., 1999;

Lee, Kaul, and Balents, 2010). The atoms are arranged as Co-O-Co-O-Co in the c-direction, bonded ionically with ionization Co²⁺ and O² . The interactions in the a-b plane happen on longer paths with Co-Nb-Nb-Co-Nb-Nb-Co in the a direction and Co-Nb-Co in the b direction, where each ion is connected by O² (Kobayashi, Mitsuda, et al., 1999). The canting angle ✓ along each chain is29.6 (Cabrera et al., 2014a).

J₁ J2

a

b

J_z J_z⁰

✓ c

(a) (b)

Figure 3: (a) Schematic structure of CoNb2O6 with each point corresponding to a Co²⁺ ion. The interactions between nearest neighbors along each chain are ferromagnetic, while the others are antiferromagnetic withJ_z |J_z|>|J₁|>|J₂|. (b) The three dimensional chain structure showing a90 exchange path between adjacent cobalts due to the octahedral oxygen structure. The image is copied from Coldea et al. (2010)

Let us first consider the exchange interaction between adjacent cobalts along the chain. Since the cobalt ions are separated by diamagnetic oxygen ions, the exchange is mediated by the superexchange mechanism (Weitzel, 1976). In the case of CoNb2O6, the superexchange is found to be ferromagnetic (FM) due to an angle of 90 (figure 3) made by two Co²⁺ and a O² (Kobayashi, Mitsuda, et al., 1999;

Weitzel et al., 2000). To see how this angle makes the usually anti-ferromagnetic superexchange ferromagnetic, we need to consider the effects of direct exchange.

It is possible to show (appendix A) that direct exchange between two electrons in

different orbital states on the same atom prefer a ferromagnetic alignment. If there is a significant overlap between two orbitals, as in current case between the d_z² orbitals of Co²⁺ and the p-orbitals of O² , the alignment must be anti-ferromagnetic in order to satisfy the Pauli exclusion principle (Pavarini et al., 2012). When placed on a straight line (figure 4), we see that the overall exchange between the cobalts is anti-ferromagnetic due to the orbitals of the two cobalts overlapping with the same orbital in O² (here called p_z). In the 90 angle case (figure 5), the d_z² orbitals overlap with separate orbitals (here called p_x and p_y) (Pavarini et al., 2012; M.

Suzuki and I. S. Suzuki, 2009). This causes the overall superexchange interaction to be ferromagnetic.

Figure 4: 180 superexchange through an intermediate O² atom.

Figure 5: 90 superexchange through an intermediate O² atom.

The other interactions are the anti-ferromagnetic next-to-nearest neighbor inter- action along the chains, as well as the anti-ferromagnetic interchain interactions in thea-bplane with the weakest along the adirection (Weitzel et al., 2000; Kobayashi, Mitsuda, et al., 1999; Cabrera et al., 2014a). The second interaction along the chains arises from the zigzag structure shown in figure 3 (Kj¨all, Pollmann, and Moore, 2011;

Cabrera et al., 2014a). The interaction is proposed to be antiferromagnetic, and al- though not discussed, we may assume it to be a result of direct exchange between two atoms.

In the zero field limit, the material exhibits 3 magnetic phases at low temper- atures (Scharf et al., 1979). At T_{IN C} = 2.95K, the material transitions from a paramagnetic (PM) phase to an incommensurate spin density wave (INC) phase.

The second phase transition happens at T_AF = 1.97K, where it transitions into a non-colinear, commensurate anti-ferromagnetic (AF) phase. The noncolinearity will not be considered in this thesis due to our Ising approximation which contains only one spin axis. The material also exhibits several magnetically induced phases for applied longitudinal, as well as transverse magnetic fields (Heid et al., 1995; Weitzel et al., 2000). We shall focus on the longitudinal fields H k c, and a phase diagram is supplied in figure 6 with phase descriptions given in table 1. Each of the ordered

phases are characterized by a locking in thecdirection, leaving us with a pronounced 2D isosceles triangular lattice behavior (Kobayashi, Mitsuda, et al., 1999). The spin flip phases SF1 and SF2 are characterized by reversals of whole chains in the c direction, and the SF2 phase is unstable due to the mixing of of two different prop- agation vectors. A way to stabilize this phase is proposed by Heid et al. (1995), which is to apply a magnetic field in the a direction. This suggestion was explored by Weitzel et al. (2000), by applying a magnetic field at different angles from the c direction in the ac-plane. In these experiments, the SF2 phase stabilized with propagation vector 0,¹₂,0 , but a new phase appeared calledSF3 with propagation vector (0,¹₄,0) characterized by a ""## pattern along the b direction.

Figure 6: A sketch of the magnetic phase diagram for magnetic fields applied in the c-direction.

The original diagram can be found in figure 1a of Weitzel et al. (2000)

Phase Propagation vector~k Illustration

AF 0,¹₂,0 "#"#

INC (0, k_y,0), ¹₂ < k_y 0.37 None

SF1 0,¹₃,0 #""#""

SF2 0,¹₂,0 mixed with 0,¹₃,0 None

P’ (0,0,0) """"

P Not ordered None

Table 1: The magnetic phases of CoNb2O6in an applied longitudinal fieldH kcwith propagation vectors and illustrations of the different magnetic phases. Note that the illustration shows the spin alignments along the crystallographic b direction, which is the easy-axis of the material in theab-plane.

2.4 A summary of existing models for CoNb2O6

We shall consider a classical Ising approach to this material where we consider all interactions to be of the Ising type. Before stating the hamiltonian, a number of considerations has to be done. As the material has a zig-zag structure along the chains in the c-direction, one should consider next-to-nearest interactions (Kj¨all, Pollmann, and Moore, 2011; Cabrera et al., 2014a). The interactions in theab-plane are not of equal strength due to the different exchange path lengths, and have a ratio J₂/J₁ = 0.75 Cabrera et al. (2014a) and Kobayashi, Mitsuda, et al. (1999). Given these considerations, a general 3D model to describe CoNb2O6 in a longitudinal magnetic field is

H = J_zX

r

S_r^zS_r+c/2^z J_z⁰ X

r

S_r^zS_r+c^z J₁X

r

S_r^zS_r+b^z

J₂X

r

⇣S_r^zS_r+(a+b)/2^z +S_r^zS_{r+(b a)/2}^z ⌘

h_zX

r

S_r^z (2.14) whereris the position of each Co²⁺ion, and the different exchange constants couples Co²⁺ ions as shown in figure 3 with J_z >0 and J₁, J₂, J_z⁰ < 0.

There exists two articles, Cabrera et al. (2014a) and Kobayashi, Mitsuda, et al. (1999), that list exchange constants for CoNb2O6 while considering the full 3D structure with J₁ 6= J₂ (see table 2a). Other articles use lower dimensional models, perfect triangle models (ab-plane), or include additional contributions such as dipolar energy. The origins of the sets we will use are different as in what they try to accomplish. The sets from Kobayashi et al. tries to reproduce the magnetic phase

diagram for magnetic fields along the c direction, while the sets from Cabrera et al.

tries to model spin waves in a transverse field. There are some peculiarities to these sets which we will address now.

The difference between the two Kobayashi sets is the value of J₂. This comes from two different considerations of how the spins interact along the (a±b)/2 directions. The first set (Kobayashi 1) is the one used in their model of CoNb2O6, but it doesn’t take into account the canting angle ✓ of the zig-zag chains. When considering the canting angle, the spin interaction gets weaker as only a portion of the spins are on the same axis. The value of J₂ in the “Kobayashi 2” set is modified by a factor of cos 2✓ to have the correct “amount” of spin interacting. This value of J₂ also reproduces an observed propagation vector(0,0.37,0) through mean field calculations temperatures right below T = 2.95 K (Kobayashi, Mitsuda, et al., 1999). The second thing to notice is that both Kobayashi sets do not include a next-to-nearest neighbor interaction term, and no specific reason is given for this choice.

The Cabrera et al. sets we’re all produced by a ² best fit from data retrieved at high applied transverse magnetic fields (7 and 9 T), causing a breakdown of the magnetic ordering. The authors themselves list the values as “representative”, only giving the ratio J₂/J₁ with reasonable accuracy. They also discuss various other models including restricted ones (J_z⁰ = 0), but claim they didn’t get a good fit.

Set J_z J_z⁰ J₁ J₂ S

Kobayashi 1 0.6015 K 0 K -0.0508 K -0.0812 K 3/2 Kobayashi 2 0.6015 K 0 K -0.0508 K 0.75J₁ 3/2 Cabrera 1 2.19 meV -0.29 meV -0.031 meV 0.77J₁ 1/2 Cabrera 2 2.4 meV -0.36 meV -0.036 meV 0.75J₁ 1/2

(a) The original sets of coupling constants

Set J_z J_z⁰ J₁ J₂

Kobayashi 1 1 J_z 0 -0.0845 J_z -0.135J_z Kobayashi 2 1 J_z 0 -0.0845 J_z 0.75J₁

Cabrera 1 1J_z -0.13J_z -0.014 J_z 0.77J₁ Cabrera 2 1J_z -0.15J_z -0.015 J_z 0.75J₁

(b) Table (a) written in terms ofJz for each set.

Table 2: Overview of the coupling sets.

in the value of J_z. This difference is also clear when considering the ratio J₁/J_z for each set (table 2b). In order to better compare the models, we should rewrite J_z ! J_z 2J_z⁰ in order to have an effective exchange along the chains. The factor of 2 appears because of two next-to-nearest interactions trying to weaken the interaction between two adjacent spins. We may also multiply the resulting couplings by S² in each set, and end up with the sets written in terms of the same spin variable S = 1.

The resulting values are summarized in table 3a and 3b.

The difference in J_z across the sets has decreased, but the values are still clearly separated. This difference is most likely due to the approach taken by the authors.

Mean field calculations are known to not be accurate in three dimensions (Kopietz, Bartosch, and SchÃŒtz, 2010, p 30), which puts uncertainty on the values provided by Kobayashi et al. . On the other side, the values given by Cabrera et al. are

“representative” and there exists many other sets that also provide good fits for the

2 method (Cabrera et al., 2014a). The ratios are J_z^{Cab 1}

J_z^Kob = 4.67

1.353 ⇡ 3.45, J_z^{Cab 2}

J_z^Kob = 4.87

1.353 ⇡ 3.60 (2.15) One interesting thing to notice, is that the values for J₁ and J₂ seem to be in fair agreement with each other despite different methods used to get them. As seen in table 3b, the difference in J_z causes them to vary a lot from each other when considering ratios. As we will see in the next sections, the ratios determine the separation between T_{IN C} and T_AF, from which we will determine what set to use in our simulations.

Set J_z J_z⁰ J₁ J₂ S Kobayashi 1 1.353 K 0 K -0.1143 K -0.1827 K 1 Kobayashi 2 1.353 K 0 K -0.1143 K 0.75J₁ 1 Cabrera 1 4.67 K 0 K -0.090 K 0.77J₁ 1 Cabrera 2 4.87 K 0 K -0.10 K 0.75J₁ 1

(a) The different sets brought to the same form by absorbing the spin information into the coupling constants and redefining Jz in order to absorb J_z⁰ (see discussion above). The conversion from meV to K was done by using1meV⇡11.6 K

Set Jz J_z⁰ J1 J2

Kobayashi 1 1 0 -0.0845 J_z -0.135J_z Kobayashi 2 1 0 -0.0845 J_z 0.75J₁

Cabrera 1 1 0 -0.019J_z 0.77J₁ Cabrera 2 1 0 -0.021J_z 0.75J₁

(b) Table (b) written in terms ofJz(see discussion above)

Table 3: The sets brought to the same form

3 The simulation

3.1 Deciding on a model

In order to decide on a model, we wish to get the correct statistics for our system.

We have chosen to model the system as straight chains on a isosceles triangular lattice, even though there exists a canting angle in the a-c plane (figure 3). This simplification is justified by having a smaller J_z to compensate for the fact that canting angle reduces the amount of spins interacting along the chains, and it is the reason for the J_xy term we discussed in the previous sections. We shall also consider all the coupling constants in terms of J_z and set J_z = 1. In order to choose a set, we perform an exact numerical calculation of the heat capacity for all sets on a small lattices, and compare them to the experimentally measured ones (Weitzel et al., 2000;

Liang et al., 2015). The set chosen is then used in a Monte Carlo simulation and compared to the exact solution. All calculations are done with periodic boundary conditions in order to best approximate the thermodynamical limit (Newman and G. T. Barkema, 1999).

The partition functions (figure 7) show all the features we need in our model

models are at what temperatures these transitions happens as well as the shape of the heat capacity. The phase transitions discussed in section 2.3 satisfy T_IC/T_AF ⇡ 1.5, and although this ratio cannot be validated from a2⇥2⇥x system, we will at least be able to choose a system size that shows what seems to be the right shape. In addition to that, we can make considerations as to whether or not our single-spin-flip algorithm will work given the location of the phase transitions.

All sets show a translation of the extrema towards lower temperatures as the system size increases along the c direction, but the Cabrera sets translates more than the Kobayashi sets. This is due to competing interactions along the chains make it easier for spins to flip, and thus a lower temperature is needed for the algorithm to lock the chains. Another feature of the Cabrera sets is the very low temperature at which the AF phase locks in, seen as the local extrema to the left.

When increasing the system size, we see that the local extrema does not move and thus, if we assume the sets have the correct ratio for large systems, the global extrema will be translated to very low temperatures. When doing Metropolis Monte Carlo with a single spin flip algorithm on a frustrated triangular lattice, the system usually freezes at low temperatures and extremely long runs are needed to get the correct statistics (Wang, Sterck, and Melko, 2012). One could argue that such a scheme could be solved by cluster algorithms, but they tend to freeze when turning on an external magnetic field, and are thus not applicable to our goal Newman and G. T.

Barkema (1999). The conclusion is the dismissal of the Cabrera sets due to our single spin-flip algorithm.

The two Kobayashi sets both have good temperature regions with the ratios that seem to converge on the right temperatures. From those and from the heat capacity gotten from experiments, we see that the right tendencies seem to happen for systems with dimensions a⇥a⇥ ³₂a, and we will choose such dimensions to run simulations on. Due to the value of J₂ in the first Kobayashi set, we deduce that a |J₂| > |J₁| will lead the system to a ferromagnetic ordering along the b direction which is not compatible with the observed propagation vector (0,¹₂,0). We conclude that the set we should be using is the “Kobayashi 2” set where the canting angle along the chains have been taken into considerations when scaling downJ₂. This is further confirmed as the right choice by a later paper by some of the same authors (Kobayashi, Hosaka, et al., 2014). Another confirmation comes from the ratioJ₂/J₁ = 0.75 for the scaled down version ofJ₂ which reproduces the right propagation vector ⇡ (0,0.37,0)right below T_IC (Cabrera et al., 2014a; Kobayashi, Mitsuda, et al., 1999).

0 0.5 1 1.5 2 2.5 3 Temperature [K/J

z] 0

0.1 0.2 0.3

Heat Capacity [J z]

2x2x2

Kobayashi 1 Kobayashi 2 Cabrera 1 Cabrera 2

0 0.5 1 1.5 2 2.5 3

Temperature [K/J z] 0

0.1 0.2 0.3 0.4 0.5

Heat Capacity [J z]

2x2x3

Kobayashi 1 Kobayashi 2 Cabrera 1 Cabrera 2

0 0.5 1 1.5 2 2.5 3

Temperature [K/J_z] 0

0.2 0.4 0.6

Heat Capacity [J z]

2x2x4

Kobayashi 1 Kobayashi 2 Cabrera 1 Cabrera 2

Exact heat capacity for different coupling sets and system sizes

0 0.5 1 1.5 2 2.5 3

Temperature [K/J_z] 0

0.2 0.4 0.6 0.8

Heat Capacity [J z]

2x2x5

Kobayashi 1 Kobayashi 2 Cabrera 1 Cabrera 2

Figure 7: Numerical calculation of the heat capacity per spin from the partition function for different system sizes

3.2 The algorithm

The current work relies on a single spin flip Metropolis Monte Carlo simulation.

The code used is a direct implementation as can be found in any introductory book on Monte Carlo methods in statistical mechanics. We will therefore present the theory, and considerations made when writing the code (section 3.2.1), and let the verification (section 3.2.2) show that it produces the right magnetic phases and statistics. At the end (section 3.2.3), we make considerations as to how to look for domain walls.

3.2.1 A brief summary of the Metropolis Monte Carlo algorithm

The Metropolis algorithm deals with systems where the correct distribution is known

edge allows the algorithm to do random walks based on probabilities determined from the Boltzmann distribution. The result is that averages can be calculated by an equal weighing of each measurement as the frequency of configurations made are already sampled by the Boltzmann distribution (Metropolis et al., 1953).

As already mentioned briefly in the previous section, there exists many different implementations to spin systems based on multiple updates at once (cluster algo- rithms), single updates, and combinations of the two (Wang, Sterck, and Melko, 2012). Due to the presence of a magnetic field, the present work is done using a single spin flip algorithm, with the core algorithm (Newman and G. T. Barkema, 1999) being:

1. Create a random initial configuration

2. Place yourself at a random location in the system

3. Evaluate the difference in energy E from flipping the spin 4. Compute a random number w

5. If e ^E  w, accept the spin flip. If not, reject the flip.

6. Repeat steps 2-5 until satisfied.

The steps 2-5 is considered one Monte Carlo step. If N is the number of spins in the system, we define a Monte Carlo sweep as N steps. It is then usual to measure the magnetization and energy of the system for each sweep and generate statistics from these measurements.

In order to get the right statistics, one should discard an amount of sweeps before storing data. This is due to the Boltzmann distribution not being flat, and thus we need to be sure our system is in equilibrium for the respective temperature. The equilibration time (number of discarded sweeps) is usually determined by looking at the time evolution of a run with two or more systems. We reach equilibrium when the energy and magnetization become stationary around the same value for both systems. At this point, a number of sweeps equal to the equilibrium time are used to estimate the correlation time (also called autocorrelation) in the magnetization data. The correlation time is a measure of how often we can make uncorrelated measurements. Ideally, one would make very long runs at all temperatures before finding the correlation time. The problem with this approach is the time consumed

when calculating. Some temperatures don’t need to be run for a long time in order to get good statistics, while other temperatures need much more computation time due to rapidly increasing correlation times at phase transitions as well as very low temperatures. We therefore need to estimate the correlation time at a relatively early stage in the run such that we can have a notion of how long the total run should be at each temperature. We may thus take the early calculated correlation time as a measure of how long some temperatures need to run relative to others.

The next problem is how to deal with correlations within the datasets. If the data is very correlated, the entire analysis might be wrong, and thus one should be careful.

Correlation in the dataset means the algorithm spends considerable time locked in states that produce magnetizations such that sign(m(t) hmi) doesn’t change. This is typical for freezing processes, but eventually, some change will happen before the system freezes at some new magnetization. If the run is done long enough, the mean can be well defined.

In order to do a proper error analysis, one should consider some methods to deal with the correlations in the generated set. Some of the most used are the binning method, the jackknife method and the bootstrap method (Newman and G. T. Barkema, 1999). The bootstrap method is the most suitable for large datasets and is the one we will use. Ideally, we would run many separate runs in order to generate many mean values and get the standard error from these, but this is generally too time consuming. By bootstrapping the dataset, we simulate many runs without actually having to do them. This is done by resampling n random measurements from the dataset and allowing for one measurement to be chosen more than once. The number n is the total number of measurements in the dataset, and as such, we can create many different means from the dataset fast. The standard error is then calculated as the standard deviation of the resampled (bootstrapped) means (Newman and G. T. Barkema, 1999).

3.2.2 Verifying the algorithm

In order to verify the algorithm, a run lasting 1500 correlation times was preformed on a 2 ⇥2⇥ 3 system and compared to the exact numerical solution using S = 1.

Units were chosen such that ~ = k = 1. Plots of the heat capacity, magnetic susceptibility as well as the magnetization and it absolute value were plotted (figure

heat capacity down to temperatures of T^⇤ = kT /J_z = 0.4, and a slightly less good match for the magnetic susceptibility. A possible explanation for the less accurate results for the magnetic susceptibility compared to the heat capacity is that 1500 correlation times is too short of a run. This is evident by some temperatures in the range 0 < T^⇤ < 0.5 having the proper variation in magnetization and energy while others have zero. Another thing to consider is the spread of values in the magnetic susceptibility. The error bars show little uncertainty in the values calculated, but the values may differ from the exact solution due to the geometrical frustration in CoNb2O6. This frustration does not impact the heat capacity as much as the the magnetic susceptibility due to the ground state energy being highly degenerate.

Also, some spin flips leave the energy unchanged, while the magnetization changes for every flip.

The temperature range for T^⇤ < 0.4 shows that the algorithm suffers from freeze- ins. This effect is due to exp ( E) !0 for T^⇤ ! 0, which denies the algorithm from flipping a spin in most cases. This effect is very usual for frustrated systems at low temperatures when applying a single spin-flip algorithm Wang, Sterck, and Melko (2012), and we cannot remedy this as the magnetic field introduced later leaves cluster algorithms faulty. As such, ergodicity is violated for these temperatures, and although some results in this temperature range will be presented, they may not be valid.

0 0.5 1 1.5 2 0

0.1 0.2 0.3

0.4 Heat Capacity per spin

Monte Carlo Simulation Exact solution

0 0.5 1 1.5 2

0 0.5 1

1.5 Magnetic Susceptibility per spin

Monte Carlo Simulation Exact solution

0 0.5 1 1.5 2

-0.5 0 0.5

Magnetization

Magnetization per spin

System size = 2x2x3

0 0.5 1 1.5 2

0 0.1 0.2 0.3 0.4

|Magnetization|

|Magnetization| per spin

Monte Carlo Simulation Exact solution

Figure 8: Plots showing the different statistics produced by the simulation in zero magnetic field.

The temperature given isT^⇤=kT /J_zand the run lasted for 1500 correlation times per temperature on a 2⇥2⇥3 lattice.

In order to settle on a value for J_z, a run was done for a larger system with dimensions 8⇥8⇥12 and run length equal to 3000 correlation times (figure 9). This larger system shows properties closer to the measured heat capacity and magnetic susceptibility (Weitzel et al., 2000; Nandi, Prabhakaran, and Mandal, 2016; Liang et al., 2015). Just right below T^⇤ = 1, we see a change in the slope of the magnetic susceptibility (T_{IN C}^⇤ ⇡ 0.85) corresponding to the transition from the PM phase into the INC phase at T = 2.95 K. There is also a change in slope for the absolute value of the magnetization occurring at T_AF^⇤ ⇡ 0.58. This signals the transition into the AF phase. If this model indeed is good, it should reproduce the ratio T_{IN C}/T_AF ⇡ 1.50. We have T_{IN C}^⇤ /T_AF^⇤ ⇡ 1.47, which is close. Given that T_{IN C}^⇤ will move towards lower temperatures when the system size is increased further, we may assume the ratio to converge on the right value. This effect allows us to use

This value is close to J_z = 3.5 K given by Weitzel et al. (2000) who considered a very realistic model including many different types of interactions. We may assume that T_{IN C}^⇤ won’t change dramatically, and T_AF^⇤ even less so, allowing us to take J_z = 3.5 K as the correct value for our model. The difference between the ratios given in table (3b) is the relative strength of the interchain interactions. Due to the ratio T_{IN C}^⇤ /T_AF^⇤ for the Kobayashi set, we may assume that the stronger interchain couplings compensate for the missing terms in our hamiltonian compared to the one by Weitzel et al. (2000). The conclusion is that our model is able to capture the thermodynamical statistics of CoNb2O6 for temperatures T J_z/2 in the zero field limit.

0 0.5 1 1.5 2 0

0.2 0.4 0.6

0.8 Heat Capacity per spin

0 0.5 1 1.5 2

0 2 4

6 Magnetic Susceptibility per spin

0 0.5 1 1.5 2

-0.05 0 0.05

Magnetization

Magnetization per spin

System size = 8x8x12

0 0.5 1 1.5 2

0 0.02 0.04 0.06

|Magnetization|

|Magnetization| per spin

Figure 9: A run done for 3000 correlation times on a 8⇥8⇥12 lattice.

A run was also made for different magnetic fields as shown in figure 10. When looking at this, we must remember that our algorithm suffered from freeze-ins at T^⇤ < 0.5, and as such, the results for T^⇤ 2 {0.33,0.35,0.4} may not be valid for low fields. What we can see, is that for T^⇤ = 0.35, we have behavior at low fields close to what we would expect theoretically. It shows a transition from the AF ground state at h_{z,IN C} = 0.05J_z and a transition into the saturated P’ phase at hz,P⁰ = 0.5Jz, giving a ratio of hz,P⁰/hz,IN C = 100 which is exactly the same as

h_{AF F R} = 315 Oe and h_{F R P M} = 3150 Oe, which corresponds to magnetic fields strengths ofB_{AF F R} = 0.0315 T and B_{F R P}⁰ = 0.315T. This is consistent with the magnetic phase diagram in figure 6. From these values, we may find gµ_B

gµ_B = h_z

B = 0.05·3.5K

0.0315T ⇡ 5.56K

T, µ_B = Bohr magneton (3.1) For low temperatures, we see two stable magnetic phases for 0.1  h_z < 0.25 and 0.25  hz < 0.35, with the first being consistent with the induced INC phase and the second being the ferrimagnetic SF1 phase (figure ). Although a structure of spins as ""# would have a magnetization of 1/3, we notice a value of ⇡ 0.4. This is mainly due to 8 atoms along the b axis is not commensurate with the current spin distribution., The result is an extra frustration in the simulation with the magnetic field tending to align them up, yielding a higher magnetization.

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8

Magnetic field h

z [J

z] 0

0.2 0.4 0.6 0.8 1

|Magnetization|

|Magnetization| vs magnetic field for different temperatures

T* = 0.33 T* = 0.35 T* = 0.4 T* = 0.5 T* = 0.6 T* = 0.7 T* = 1

Figure 10: Magnetization per spin as a function of magnetic fields plotted for different tempera- tures. The lattice size is 8⇥8⇥12 and the run lasted for 2000 correlation times.

Snapshots of the lattice taken for various magnetic fields and temperatures are supplied in figure 11 and 12. The snapshots for T^⇤ = 0.4 and 0.6 shows that the

algorithm succeeds in creating the right magnetic phases. Further, the lattices at T^⇤ = 0.6 and T^⇤ = 0.8 show an increase the number of domain walls with increasing temperature as expected. The ground state for in table 11a is not the same as the one reported with propagation vector(0,¹₂,0), but it has the same total energy. This is because any translation of a basal chain in the b direction, has the same degree of frustration. We also notice that for T^⇤ = 0.8, we see many domain walls created, along with what seems to be a very uniform local neighborhood at each site.

(a) AF ground phase (b) SF1 phase

(c) SF1 !P’ transition (d) P’ Phase

Figure 11: Different magnetically induced phases for T^⇤= 0.4

(a) SF1 phase (b) P’ phase with thermal excitations

(c) INC phase close to the P phase (d) P’ phase with thermal excitations

Figure 12: Snapshots of different phases for temperaturesT^⇤ = 0.6and 0.8.

3.2.3 Looking for domain walls and the Wannier-Zeeman ladder.

With the validation in place, we may turn to the problem of finding suitable condi- tions for the observation of magnetic Bloch oscillations. When looking for domain walls and the associated cost of moving them one site, we will apply the following algorithm (written for a general chain length L_z):

1. Generate a set of uncorrelated lattice configurations

2. Sweep through each lattice isolating one chain at the time

3. Calculate total spin along the chain S_chain =P

iS_i 4. If |S_chain| = L_z, move to the next chain.

5. Else, sweep over the chain looking for domain walls given by S_i + S_i+1 = 0 (periodic boundary conditions implied).

6. When found, calculate the cost | E| associated with flipping S_i and S_i+1. 7. Count the number of domain walls found and add them to a total. Do the same

for the number of times each separate | E| has been detected.

8. Move to the next chain and repeat 2-7

It was noted in section 2.1.3 that if we are to find a unique Bloch frequency, the cost of moving a domain wall should be as uniform as possible. This suggests we calculate the frequency of occurrence for each value of | E|. The calculation is done through dividing the number of times one specific cost was detected by the total amount of costs calculated. In addition to this, we ought to calculate the density of domain walls per chain. This is important because if few domain walls were created, one can end up with a high frequency of occurrence. This would suggest to look for magnetic Bloch oscillations where there are too few domain walls, and would not show up as a significant peak in neutron scattering experiments. Another problem is colliding domain walls when the density is too high. If this happens, the oscillatory motion is destroyed. The density is calculated through

⇢_DW = DW_tot

N ·L_z (3.2)

where DW_tot is the total number of domain walls across N lattice configurations with chain lengthL_z. We hope to find a parameter set that orders the material such that the interchain neighborhoods are as uniform as possible. As seen in section 2.1.3, the two magnetic phases that has only one type of neighborhood, are the AF phase and the saturated P’ phase (figure 6). We will therefore compare our results for fieldsh_z < 0.3J_z to| E|AF (2.11), and with| E|P M (2.12) for fieldsh_z 0.3J_z and hope to see a good match.

3.3 Simulation results

The simulation was set on a lattice with dimensions 8 ⇥8 ⇥ 12 as in the previous sections, due to it offering a nice compromise between calculation time and accuracy.

The run lasted for 5 days on 11 CPU cores, where the algorithm discussed above was applied for T^⇤ between 0.4 and 0.8. In order to compare with the theoretical costs

| E|AF (2.11) and | E|P⁰ (2.12), the cost of moving a domain wall was calculated using S = 1/2 instead of S = 1. The implications of this choice will be discussed later.

The data produced by the simulation is plotted vs h_z in figure 13, 14 and 15. We have chosen to not consider the T^⇤ = 0.4 data as the domain wall density per chain is nearly zero as seen in the density plot. The size of the circles is determined by the frequency of occurrence, and if the frequency is greater than 0.4, its value is plotted inside the circle. The lines correspond to the evolution of the energy cost according to | E|AF for h_z < 0.2, and | E|P⁰ for h_z 0.3 as already mentioned above.

0 0.2 0.4 0.6 0.8 1

0 0.5 1 1.5 2

2.5 Domainwall density per chain

T* = 0.4 T* = 0.5 T* = 0.6 T* = 0.7 T* = 0.8

Figure 13: Domain wall density per chain ⇢_DW vs. magnetic field h_z

0 0.2 0.4 0.6 0.8 1 0

0.2 0.4 0.6 0.8

Energy cost of moving a domainwall one site, T* = 0.5 1

0.67 0.66 0.42 0.65

0.96 1 1 1 11

1 1 11

1 11

1 11

1 11

0 0.2 0.4 0.6 0.8 1

0 0.2 0.4 0.6 0.8

1Energy cost of moving a domainwall one site, T* = 0.6

0.48 0.56

0.59 0.91

0.99 1

1 1

1

Figure 14: The results from moving domain walls. The circle sizes are scaled by the occurrence frequency of calculated energy costs. The lines represent the evolution of the cost given the AF ground phase and the saturated PM phase.

0 0.2 0.4 0.6 0.8 1

0 0.2 0.4 0.6 0.8

1Energy cost of moving a domainwall one site, T* = 0.7

0.48 0.46

0.53 0.84

0.96 0.98

0.99 0.99

1

0 0.2 0.4 0.6 0.8 1

0 0.2 0.4 0.6 0.8

1Energy cost of moving a domainwall one site, T* = 0.8

0.47 0.75

0.9 0.95

0.97 0.98

0.99

Figure 15: Continuation of figure 14

The results all show a plethora of energy costs at lower magnetic fields h_z  0.4, and a higher degree of ordering for h_z 0.4. This implies that in a field interval close to h_z = 0.4, major changes in magnetization take place and may result in many domain walls. By looking at the density plot and the magnetic