Benaderende Bethe ansatz voor een spin ladder

(1)

Ghent University Faculty of Sciences Physics and Astronomy

Approximate Bethe ansatz for a spin ladder

Jacob Lamers

Promotor: dr. L. Vanderstraeten

Copromotors: dr. J. De Nardis & Prof. dr. J. Haegeman

Thesis supervisor: M. Van Damme

In order to obtain the degree of Master of Sciences: Physics and Astronomy 2019-2020

(2)

(3)

Faculty of Sciences

Physics and Astronomy

2019–2020

Approximate Bethe ansatz for a spin ladder

Jacob Lamers

Promotor: dr. L. Vanderstraeten

Copromotors: dr. J. De Nardis & Prof. dr. J. Haegeman Thesis supervisor: M. Van Damme

In order to obtain the degree of Master of Sciences Physics and Astronomy

(4)

Preface

Acknowledgements

First of all, the author would like to thank M. Van Damme for navigating me through this thesis and his input of both ideas and data and L. Van-derstraeten for his suggestions and providing me with the necessary tools to conduct this research. Furthermore, we would like to thank J. De Nardis for his help and expertise and of course the VVN for ensuring this thesis reaches a broader public even in these times of Corona.

De auteur geeft de toelating deze masterproef voor consultatie beschik-baar te stellen en delen van de masterproef te kopiëren voor persoonlijk gebruik. Elk ander gebruik valt onder de beperkingen van het auteursrecht, in het bijzonder met betrekking tot de verplichting de bron uitdrukkelijk te vermelden bij het aanhalen van resultaten uit deze masterproef.

15 juni 2020

The author gives permission to make this master dissertation available for consultation and to copy parts of this master dissertation for personal use. In the case of any other use, the limitations of the copyright have to be respected, in particular with regard to the obligation to state expressly the source when quoting results from this master dissertation.

(5)

1 Introduction

Classical integrability is a very well understood field of physics with a lot of useful applications. It has received a lot of attention in its days and it comes thus as no surprise that there have been several mathematical tools developed to solve integrable models.

With the rise of quantum mechanics, attempts were made to find a quan-tum version of integrability. This turned out to be very difficult however and remains a research topic until today. Many models are therefore just assumed to be integrable more than they are proven to be. Although there is no generally accepted definition of quantum integrability as of yet, many definitions have been proposed. (5) In this dissertation, we will use one of these definitions: A model is said to be integrable if the scattering it supports is nondiffractive. (21)

This definition of quantum integrability will be elaborated upon in section 2. After this, in the same section, we introduce the thermodynamic Bethe ansatz and show how it can be used to numerically solve integrable models. In section 3, the method introduced in section 2 is applied to the inte-grable Lieb-Liniger model to test the code.

After this, we apply the thermodynamic Bethe ansatz to the Ising model with transverse field in section 4. Here, we also study the obtained effective velocity of the quasiparticles and show that it corresponds to the observed lightcone dynamics in quenched quantum systems.

This is followed by some attempts at applying this method to the non-integrable spin-1/2 Heisenberg antiferromagnetic ladder without magnetic field in section 5. After this, we apply the thermodynamic Bethe ansatz to this ladder model with a magnetic field instead. We show that the thermo-dynamic Bethe ansatz can indeed be used to solve non-integrable quantum systems in the low temperature limit, but that its usefulness needs some clarification.

In appendix A, we give a brief summary of this dissertation in both En-glish and Dutch. Because of considerations of readability of the main text, we put the lengthy calculation of the Ising model in appendix B. In appendix C is another attempt at solving the problems encountered in section 5 that was not pursued further in this dissertation. Finally, in appendix D, there is a text that is meant to communicate this dissertation to a larger audience without a background in physics. This text was written in Dutch, since it is going to be shared on the website of the VVN student organisation, which main audience is Dutch speaking.

(7)

2 RUDIMENTARY NOTIONS

2 Rudimentary notions

We start this section with a brief introduction of integrability. The main focus will be on the definition of quantum integrability of nondiffractive scattering (21) and its connection to classical integrability. After this, we develop the thermodynamic Bethe ansatz method and clarify its relation to integrability.

2.1 Integrability

2.1.1 Classical integrability

In order to discuss classical integrability, let us assume a system, described by a Hamiltonian H(q, p) where q = (q1, ..., qN) are the coordinates and

p = (p1, ..., pN) are the momenta of the particles. If we are working in one

dimension, N represents the total number of particles in the system1_{. The}

equations of motions are then given by ˙ q = ∂H ∂p ˙ p = −∂H ∂q .

The time evolution of any function of these coordinates and momenta L(q, p) that does not depend explicitly on time, can then be written down, even without solving these coupled differential equations.(15)

dL

dt = {L, H}. Here we introduced the Poisson bracket as

{f, g} = N X i=1 ∂f ∂qi ∂g ∂pi − ∂f ∂pi ∂g ∂qi .

Conservation of these quantities can thus be written down as {L, H} = 0. This system is said to be integrable if there exist N independent conserved quantities L = (L1, ..., LN), with {Li, Lj} = 0, ∀i, j ∈ [1, ..., N]. When this

is the case, one can perform a canonical tranformation from (q, p) to the so-called action-angle coordinates (J, Θ), such that the action coordinates are functions of just these conserved quantities

Ji ≡ Ji(L1, ..., LN).

1_{In general} N

d represents to total number of particles, where d is the number of

(8)

2.1 Integrability 2 RUDIMENTARY NOTIONS

The new Hamiltonian H0 _{equations become}

˙ J = ∂H 0 ∂Θ ˙ Θ = −∂H 0 ∂J .

The first equation equals zero, since the Ji only depend on conserved

quanti-ties and are thus conserved quantiquanti-ties themselves. This means that ∂H0

∂Θ = 0

and the Hamiltonian can not depend on the angle coordinate Θ. The second equation then becomes ˙Θ = −∂H0

∂J := ω(J ) = constant. (22)

It is clear that the equations of motion for an integrable system become very simple. It is then not surprising that integrable classical systems have general mathematical methods to get the exact solution.

Another property of integrable systems is that chaotic behaviour cannot occur. One of the main properties of chaotic systems is their sensitive de-pendence on initial conditions. This means that trajectories starting from almost overlapping initial points will turn out to be very different.(20) 2.1.2 Quantum integrability

A plausible attempt for a definition for quantum integrability would be to replace the Poisson brackets with commutators. In order to prove quan-tum integrability, one would then have to find N independent operators L = (L1, ..., LN) that commute with each other and the Hamiltonian. This

however turns out not to be the right way to do things (5; 21). The main problem is that according to this definition, every quantum system with a finite-dimensional Hilbert space is to be called integrable. Indeed, if one di-agonalises the Hamiltonian to get N eigenvectors |Ψii, the set of projectors

|Ψii hΨi| already satisfies this definition. Indeed:

h |Ψii hΨi| , |Ψji hΨj| i = |Ψii hΨi|Ψji hΨj| − |Ψji hΨj|Ψii hΨi| = 0 and h H, |Ψii hΨi| i = H |Ψii hΨi| − |Ψii hΨi| H = E |Ψii hΨi| − E |Ψii hΨi| = 0.

Instead, one might ask what one would expect from an integrable system. First of all, assume that the total particle number N is conserved. This way the thermodynamic behaviour of the system can be found by letting N → +∞ and L → +∞, while holding d = N/L constant. Here L is the length of the system. Secondly assume that there is only a finite amount

(9)

of particles in this infinite volume, making the density practically zero. In this asymptotic region particles can propagate freely without interacting with other particles. We will also assume that the system supports scattering. By this we mean that if the particles scattered, they will fly away far enough for them to be represented by a plane wave.

Assume a system of N particles with the general Hamiltonian

H = 1 2 N X j=1 p2_j + N X 1=j<k v(xk− xj),

where pj is the momentum operator and v is the interaction between two

particles. We now consider different cases: • N = 1

Since there is only one particle, the system is automatically in the asymptotic region. the wavefunction of the particle is therefore a plane wave ψ = eikx_{, where k is the asymptotic momentum of the particle.}

Using Hψ = ψ, results in = k2_/2_.

• N = 2

If the particles are far apart, they are both described by plane waves and we get for the total momentum and energy:

P = k1 + k2

E = (k1) + (k2).

When the particles do not interact, the wavefunction is Ψ = ei(k1x1+k2x2).

When the particles do interact, the conservation of energy and mo-mentum require that k0

1 = k2 and k02 = k1, where k0 is the asymptotic

momentum after scattering. When after scattering, the particles move far away from each other, they can still be described as a plane wave, since the system supports scattering. The wavefunction then becomes Ψ = ei(k2x1+k1x2)_e−iθ(k1−k2), where θ is the phase shift due to the

in-teraction. The first factor looks like the wavefunction for two non-interacting particles, where the two particles have switched positions (permutation). This was to be expected, since the particles are not distinguishable. One can thus write the total wavefunction as follows

Ψ = ei(k1x1+k2x2)_{− e}i(k2x1+k1x2)_e−iθ(k1−k2) (1)

=X

P

Ψ(P)ei(kP1x1+k_P2x2) (2)

(10)

where the sum goes over all permutations of the particles and Ψ(21)

Ψ(12) = −e

−iθ(k1−k2)_. (4)

• N = 3

Analogously to equation 2, one can write the wavefunction of three particles as follows Ψ =X P Ψ(P)ei(kP1x1+k_P2x2+k_P3x3) + Z Z Z k0 1<k02<k30 dk0₁dk₂0dk₃0S(k₁0, k₂0, k₃0)ei(k10x1+k20x2+k03x3)_, (5)

where the triple integral is such that it conserves total momentum and energy. The first term is the wavefunction if none or only two of the three particles scatter and is of the form that is used in the Bethe ansatz (see later), while the second term is the wavefunction when all three particles interact simultaneously. In the last term, the conservation of energy and momentum is not enough to ensure that the outgoing asymptotic momenta can be written as permutations of the incoming asymptotic momenta. Therefore the system allows diffractive scattering. As will be explained later, the Bethe ansatz method assumes only two-body scattering and can thus only be used when the second term can be written in the form of the first term i.e. when a three body interaction can be written as a sequence of two-body interactions where the sequencing does not matter.

This means that for a non-diffractive system, the Bethe ansatz can provide an exact solution. We could imagine that if there would be a third local conserved quantity next to energy and momentum, we could have enough equations to ensure that the outgoing asymptotic momenta can be written as a permutation of the incoming asymptotic momenta. So whether or not the system allows diffractive scattering depends on the amount of conserved quantities. This link between the existence of an exact solution and the amount of conserved quantities reminds us of the classical integrability case. A proposed definition of quantum integrability is then that the system has to support non-diffractive scattering. (5; 21)

One can make this a little more rigorous by introducing the so called Faddeev-Zamolodchikov operators Z (23; 16). When this operator works on the

(11)

ground state of the interacting theory |ψgsi, it creates the exact elementary

excited state:

|κi_α = Z_α†(κ) |ψgsi ,

where |κiα is an elementary excitation of type α with momentum κ. Since a

permutation of two particles can be seen as an interaction, we can write Z_α†(κ)Z_β†(κ0) |ψgsi = Sα 0_β0 αβ (κ, κ 0 )Z_β†0(κ0)Z † α0(κ) |ψgsi .

A three particle interaction can then be written as Z_α†(κ1)Z † β(κ2)Zγ†(κ3) |ψgsi = Sα 0_β0_γ0 αβγ (κ1, κ2, κ3)Z † γ0(κ₃)Z † β0(κ₂)Z † α0(κ₁) |ψ_gsi .

For an integrable system this S matrix can be written as a sequence of two-body S matrices. This can be done in two ways (see figure 1):

S_αβγα0β0γ0(κ1, κ2, κ3) = Sαβij (κ1, κ2)Sα 0_k iγ (κ1, κ3)Sβ 0_γ0 jk (κ2, κ3) (6) = S_βγjk(κ2, κ3)Siγ 0 αk(κ1, κ3)Sα 0_β0 ij (κ1, κ2). (7)

Since the sequencing itself cannot matter, these two sequences must be equal to each other. This is known as the Yang Baxter equation and is a hallmark of integrable systems.

Figure 1: The three-body scattering matrix can be decomposed in a sequence of two-body scattering matrices in two ways. These two ways must be equiv-alent. This leads to the celebrated Yang Baxter equation. (23)

(12)

Consider now a state of N particles with rapidities λ1, λ2, ..., λN. The

transfer matrix Ta0

a is then defined as (27)

Ta0

a (u) ≡ T a0

a (u|λ1, ..., λN)_bc1₁,...,c_,...,bN_N = S_aba1₁c1(u−λ1)S_aa₁2_bc₂2(u−λ2)...Sa

0_c N

aN −1bN(u−λN),

where a summation over all repeated indices is implied. The indices represent the spin labels of the particles.

It represents the scattering of an auxiliary particle with rapidity u with all N other particles (see figure 2).2

Since the S-matrices should obey the Yang-Baxter equation 7, one can write X a00_,b00 Ta00 a (u)T b00 b (v)S a0b0 a00_b00(u − v) = X a00_,b00 S_aba00b00(u − v)T_bb000(v)Ta 0 a00(u). (8)

Figure 2: Graphical representation of the transfer matrix T : An auxiliary particle with spin a and rapidity u scatters from all other particles b1, ..., bN

with rapidities λ1, ..., λN and ends with spin a0. The state |b1, ..., bNichanges

to |c1, ..., cNi as a result of these scattering processes.

Equation 8 can now be used to prove the commutation relation of these transfer matrices. For now, we assume spin-1/2 particles and thus a, b, c ∈ {−, +}. The commutation relation for spin-1 particles also holds, this is proven in section 5.

T₊+(u)T₊+(v)S₊₊++(u − v) = S₊₊++(u − v)T₊+(v)T₊+(u) T₊+(u)T₊+(v) = T₊+(v)T₊+(u),

2_{It is called a transfer matrix since this object coincides with the transfer matrix after}

(13)

2.1 Integrability 2 RUDIMENTARY NOTIONS T₋−(u)T₋−(v)S₋₋−−(u − v) = S₋₋−−(u − v)T₋−(v)T₋−(u) T−−(u)T−−(v) = T−−(v)T−−(u), T+ +(u)T − −(v)S+−+−(u − v) + T − +(u)T−+(v)S−++−(u − v) = S₊₋+−(u − v)T₋−(v)T₊+(u) + S₊₋−+(u − v)T₊−(v)T₋+(u) T+ + (u)T − −(v) = T − −(v)T++(u)

and in a similar way

T₋−(u)T₊+(v) = T₊+(v)T₋−(u) So that, for T ≡ PaTaa

[T (u), T (v)] = 0.

From this, one can see that T (u) generates conserved quantities. It turns out to be more useful to expand the natural logarithm of T (u) in a Taylor series, since this allows us to construct quasilocal operators. (8) So one can rewrite this commutation relation to

ln T (u), ln T (v) = 0 and expand ln T (u) in a series around u0:

ln T (u) =

∞

X

n

Jn(u − u0)n.

Inserting this into the commutation relation yields [Jn, Jm] = 0.

If one of these Jn is chosen to be the hamiltonian (19), these coefficients

commute with both each other and the hamiltonian. This is exactly the defi-nition of a conserved quantity. So we can see that the Yang-Baxter equation results in an infinite amount of quasilocal conserved quantities. (12; 13) This link between integrability and having a lot of conservation laws resembles the classical case.

One could now ask what this new notion of quantum integrability means in the classical limit. To answer this question, we look at a simulation conducted in Ref (21). There they argue that since quantum mechanics deals with waves (wave mechanics), we obtain the classical limit by looking at beams. By this

(14)

2.2 Bethe ansatz 2 RUDIMENTARY NOTIONS

we mean a number of trajectories with slightly different initial conditions but with the same initial asymptotic momenta. In this simulation, they integrate the equation of motion numerically ¨xj = −_∂x∂H_j, where H is the Hamiltonian

of the system. They found that for an integrable system, the beam stays collimated during the simulation i.e. the wave front remains a straight line. They also found that for non-integrable systems this is not the case and the wave front gets distorted. We thus see a classical version of non-diffractive scattering for integrable potentials. This also reminds us of the absence of chaos in classical integrable models. Indeed, in integrable models similar initial conditions produce similar trajectories.

2.2 Bethe ansatz

In this section, we introduce the thermodynamic Bethe ansatz in a similar way to (21). By the end of this subsection, it should be clear that the Bethe ansatz assumes a wavefunction in the Bethe form and thus only takes into account two-body scattering. This means that the Bethe ansatz is only exact for integrable systems and provides an exact solution of these systems. 2.2.1 Zero temperature

As was explained in the previous section, when two indistinguishable particles scatter, the interaction can be seen as a permutation of the two particles with the addition of a phase shift −e−iθ(k1−k2). Consider now a system of N

particles positioned on a circle of circumference L. Note that the choice of the circle will not matter in the thermodynamic limit (N, L → +∞, while holding N/L constant), where we will use the Bethe Ansatz method.

When a particle i moves around the circle, it scatters of all other N − 1 particles and its wavefunction will gain an extra factor

eikiL Y

kj6=ki

(−1)e−iθ(ki−kj)_.

When the particle has moved around the circle and scattered of all other particles, it will be in its original position and there should be no difference with its original wavefunction. Thus

(−1)N −1= ei2πIi _{= e}ikiL Y

kj6=ki

e−iθ(ki−kj)_.

When one is considering an even amount of bosons, then the Ii have to be

(15)

sign, so for fermions or an odd amount of bosons, the Ii have to be integers.

The previous equation can be rewritten as

ki = 2πIi L + 1 L N X j=1 θ(ki− kj). (9)

Although Ii are just numbers, they can be linked to the asymptotic momenta

ki. Indeed, when the phase shift becomes zero, equation 9 becomes

ki =

2πIi

L .

Since Ii can be used to calculate ki, the set {Ii|i = 1, ..., N } defines the

state of the N particles and Ii can be seen as a quantum number (26). The

wavefunction of the system is then Ψ =X

P

Ψ(P)ei(kP1x1+...+k_PNxN)_.

Since only fermions or hard core bosons will be considered here, the Ii must

all be different (every particle needs to have different quantum numbers). Sometimes, the asymptotic momenta will be parametrized with the ra-pidities as follows

ki = m sinh(λi) E = m cosh(λi),

where m is the gap, i.e. the minimum of the dispersion relation and the λi are

the rapidities. One can see that this paramtrization results in E2 _{= k}2_{+ m}2_.

Note that this derivation was based on a two particle interaction. The Bethe ansatz solution can thus only be applied exactly to integrable systems. Since then, the multiple particle interactions can be written as a sequence of two body interactions, where the sequencing does not matter3_.

In the thermodynamic limit, the k’s become continuous and will distribute themselves according to a density ρ(k), between the limits −q and q. The value of q is not known in advance and has to be obtained from

D = Z q

−q

ρ(k)dk; (10)

(16)

where D is the total density D = N/L.

So, the number of k’s in [k, k + dk] is given by Lρ(k)dk and one can write a sum as follows X k = L Z q −q ρ(k)dk. Applying this to the fundamental equation 9 gives

ki = 2πIi L + 1 L N X j=1 θ(ki− kj) 1 = 2π L dIi dki + 1 L X kj θ0(ki− kj) 1 = 2πρ(ki) + Z q −q ρ(kj)θ0(ki− kj)dkj,

where θ0_{is the derivative of θ. In the last line, we used the fact that the}

num-ber of quantumnumnum-bers in [k, k + dk] can be linked to the density according to dIi/dki = Lρ(ki).

2.2.2 Finite temperature

When the temperature is zero, the system will be in the ground state. We will define this ground state as being empty. This means that one can see an excitation on that ground state as being a particle. These particles will arrange themselves according to a density ρ(k). A hole is then defined as an unoccupied but allowed quantum number. These holes will then be arranged according to the hole density ρh(k). According to these definitions, one can

write

ζ(k) = ρ(k) + ρh(k), (11)

where ζ(k) is the density of allowed quantum numbers.

In the thermodynamic limit, the discrete modes k will become continuous. This means that there are many nearly degenerate states and one can switch particles and holes with k ∈ [k, k + dk] around without changing the macro state. This results in an entropy.

(17)

In [k, k + dk], one can choose Lρ(k)dk particles out of Lζ(k)dk possibilities. (Lζ(k)dk)! (Lρ(k)dk)!(Lρh(k)dk)! = exp ln(Lζ(k)dk)! ln(Lρ(k)dk)! ln(Lρh(k)dk)! = exp Ldkζ(k) ln ζ(k) − ρ(k) ln ρ(k) − ρh(k) ln ρh(k). So S = Z +∞ −∞ dkζ(k) ln ζ(k) − ρ(k) ln ρ(k) − ρh(k) ln ρh(k).

With the entropy, the pressure can be found according to p = 1 LT S − E + µN = Z +∞ −∞ dkT ζ(k) ln ζ(k) − T ρ(k) ln ρ(k) − T ρh(k) ln ρh(k) + ρ(k) µ − Υ(k),

where Υ(k) is the dispersion relation (E = R+∞

−∞ dkρ(k)Υ(k)) and µ is the

chemical potential.

One then finds the equilibrium distribution ρ(k) by maximizing p with respect to ρ. δp δρ = 0 (12) = Z +∞ −∞ dk(δρ + δρh) ln (ρ(k) + ρh(k)) (13) − δρ ln ρ(k) − δρhln ρh(k) + δρ µ − Υ(k)/T (14) = Z +∞ −∞ dkδρhln ζ ρh + δρ lnζ ρ + δρ µ − Υ(k)/T . (15) Now, in order to link δρ with δρh, we use again the fundamental equation 9.

kL = 2πI + L Z +∞

−∞

dk0ρ(k0)θ(k − k0). When the derivative is taken with respect to k, one gets

ζ(k) = ρ(k) + ρh(k) = 1 2π 1 − Z +∞ −∞ dk0ρ(k0)θ0(k − k0) ! ≡ 1 2π − Kρ,

(18)

where we used dI/dk = Lζ(k). So

δρh = −δρ I + K.

Substituting this in equation 15 0 = Z +∞ −∞ dkδρ " − K ln1 + ρ ρh − lnρ + ρh ρh + lnρh+ ρ ρ + (µ − Υ(k))/T # 0 = Z +∞ −∞ dkδρ " lnρh ρ − 1 2π Z +∞ −∞ dk0θ0(k − k0) ln1 + ρ(k 0₎ ρh(k0) ! + µ − Υ(k)/T # 0 = lnρh ρ − 1 2π Z +∞ −∞ dk0θ0(k − k0) ln 1 + ρ(k 0₎ ρh(k0) ! + µ − Υ(k)/T. This equation becomes more simple when defining ρh(k)/ρ(k) = exp{(k)/T }:

(k) = Υ(k) − µ + T 2π

Z +∞

−∞

dk0θ0(k − k0) ln1 + e−(k0)/T. (16) This (k) can be interpreted as the dressed energy. Indeed, the first term describes the additional energy the system gets when a new particle is added with momentum k. The second term then describes the extra energy coming from the interaction with the other particles already in the system.

One can now solve this equation with an iteration loop to find (k). This can then in turn be used to find the density. Indeed:

ζ(k) ≡ ρ(k) + ρh(k) (17) ζ(k) ρ(k) = 1 + ρh(k) ρ(k) = 1 + e (k)/T ₍₁₈₎ ρ(k) = ζ(k) 1 + e(k)/T. (19)

This density can now be used to calculate properties of the system, such as the total energy and the effective velocity. The total energy is obtained by

E = Z +∞

−∞

dk

2πρ(k)Υ(k) (20)

and the effective velocity (4) is defined by veff(k) =

E0dr(k) p0dr_(k) =

dEdr

(19)

where the dressing operator is

hdr(λ) = h(λ) − Z +∞ −∞ dλ0 2πθ 0 (λ − λ0)ρ0)hdr(λ0), where λ is the rapidity.

This effective velocity can be interpreted as a generalisation of the group velocity that takes scattering from quasiparticles into account. This is done by the dressing operator.(3; 7). In order to get a better understanding, we look at the following classical system. Follow a tagged particle as it moves through the gas of quasiparticles. its displacement is given by

∆x = E0(k)∆t −X

i

θ(k − ki),

where E’(k) = dE/dk is the group velocity of the tagged particle when it is freely propagating. Every time it comes across a quasiparticle i, they scatter and the tagged particle accumulates an extra phaseshift and it makes a jump, this ad hoc addition is made to incoorporate the quantum nature of these scattering processes. Since the phaseshift can be both positive or negative (depending on the value of k −ki) the influence of the quasiparticles

can make the effective distance of a particle longer or shorter. When these contributions are averaged the interpretation of the effective velocity becomes clear.(7)

hxi = veff∆t.

2.2.3 Not a low density approximation

Before we continue, we would like to mention that this result is not a low density approximation. This might seem strange at first glance, since we assumed a finite amount of particles in an infinite volume i.e. an almost van-ishing density. This way, we got the asymptotic regime, where the asymptotic Bethe ansatz takes place. The key to the solution of this problem is the virial expansion of the thermodynamics i.e. an expansion in density. The n’th or-der in this expansion represents a n-particle interaction. But for an integrable model, this interaction can be written as a series of two-body interactions. This means that the Bethe ansatz captures all orders of this virial expansion and is thus not a low density approximation. (21)

(20)

3 LIEB-LINIGER MODEL

3 Lieb-Liniger model

In this section, we apply the thermodynamic Bethe ansatz to the Lieb-Liniger model. This model describes a one-dimensional system of interacting bosons and is the non-relativistic limit of many field theories. This model is also integrable and has been used to study integrability. This model has thus been solved before (4; 8; 17), which makes it perfect for testing our method and code.

3.1 Thermodynamic Bethe ansatz

The Lieb-Liniger model is defined by following Hamiltonian (8)

H = − 1 2m N X j=1 ∂2 ∂x2 j + cX j<l δ(xj − xl),

where c is the interaction strenght and ~ = 1. When we set m = 1, we can rescale this Hamiltonian to

H = − N X j=1 ∂2 ∂x2 j + 2cX j<l δ(xj− xl).

Since the Lieb-Liniger model is non-relativistic, one can use k = mλ = λ, where λ is the rapidity and the dispersion relation is Υ(k) = k2_.

The main thing we will need in order to use the equations derived in previous section is the phase shift. In order to calculate this, we look at the two-body interaction. H = − ∂ 2 ∂x2 1 − ∂ 2 ∂x2 2 + 2cδ(x1− x2).

We can write a generic eigenstate as follows

Ψ(x1, x2) = f (x1, x2)Θ(x2− x1) + f (x2, x1)Θ(x1− x2),

where Θ(x − y) is the Heaviside step function defined as

Θ(x − y) = (

1 if x > y 0 if x < y . For f, we take the general wavefunction from equation 3.

(21)

3.1 Thermodynamic Bethe ansatz 3 LIEB-LINIGER MODEL

We now let H work on the wavefunction − ∂ 2 ∂x2 1 Ψ = −Θ(x2− x1)∂x21f (x1, x2) − Θ(x1− x2)∂ 2 x1f (x2, x1) +∂x1f (x1, x2)δ(x2− x1) − ∂x1f (x2, x1)δ(x1− x2). So HΨ = k2 1 + k22Ψ + 2δ(x2− x1)

c(A12+ A21) + i(k1 − k2)(A12− A21)

ei(k1+k2)x1_.

If Ψ is an eigenvector of H, then the second term should be zero: A12 i(k1− k2) + c = A21 i(k1− k2) − c A21 A12 = i(k1− k2) + c i(k1− k2) − c .

One can now use equation 4 in order to find the phaseshift. Rewriting the previous expression: A21 A12 = −c + i(k1− k2) c − i(k1− k2) = −c 2_{− (k} 1− k2)2+ i2c(k1− k2) c2_{+ (k} 1 − k2)2 = − " c2− (k1− k2)2 c2_{+ (k} 1− k2)2 + i 2(k1− k2)c c2_{+ (k} 1− k2)2 # . The modulus of this complex number is

s (k1− k2)4+ c4− 2(k1− k2)2c2 (k1− k2)2+ c2 + 4(k1− k2) 2_c2 (k1− k2)2+ c2 = 1.

This means that the result of the two-body interaction is a pure phase, as was expected. This phase is called the phaseshift.

A21 A12 = −e−iθ(k1−k2)_, where θ = −2 arctan 2(k1−k2)c c2_+(k 1−k2)2 1 + c2−(k1−k2)2 c2_+(k 1−k2)2 ! = −2 arctan k1− k2 c . (21) Equation 16 requires the derivative of the phaseshift:

dθ dk1 = − 2 c 1 + (k1−k2)2 c2 = −2c c2_{+ (k} 1− k2)2 . (22)

(22)

3.2 Results 3 LIEB-LINIGER MODEL

3.2 Results

Together with equation 22, one can solve equation 16 with an iteration loop to find = T lnρh(k) ρ(k) . (k) = k2− µ − T 2π Z +∞ −∞ dk0 4c 4c2_{+ (k − k}0₎2 ln 1 + e−(k0)/T.

We would now like to calculate the density in order to calculate properties of this system. Therefore, we use equation 19.

ρ(k) = ζ(k) 1 + e(k)/T =

1 1 + e(k)/T,

where ζ(k) = 1 because if there are multiple states with asymptotic momen-tum k then the phaseshift θ(k, k) = 0. This means that A12/A21 = −1 and

the wavefunction becomes zero.

A plot of the density of the Lieb-Liniger model can be seen in the figures below (figures 3 and 4).

Before we move on, we want to remark that in the code, we made following approximation: Z +∞ −∞ dkf (k) ≈ Z 10 −10 dkf (k).

This approximation is valid since f(k) → 0 when ρ(k) → 0 and it is clear from figure 3 that ρ(k) ≈ 0 for values of k outside of the interval ] − 10, 10[. In figure 5, the total energy of the system is shown as a function of temperature for different values of the systems parameters.

In figures 6 and 7, the effective velocity is plotted as a function of tem-perature and asymptotic momentum.

(23)

(a) One can see that the density goes to zero when the temperature goes to zero and µ = 0. This agrees with our defini-tion that the ground state of the system is empty (had no particles). When µ 6= 0, it becomes energetically advantageous to produce particles and the density does not go to zero anymore for small tempera-tures. When the temperatures rise, more and more particles are being created and the density rises.

(b) When c increases, the interaction between the particles increases, which makes the creation of particles less favoured even though µ = 1. So with in-creasing c, the density decreases.

Figure 3: The density for the Lieb-Liniger model as a function of temperature for different values of µ (a) (c = 1) and multiple values of c (b) (µ = 1).

(24)

(a) (b)

Figure 4: The density for the Lieb-Liniger model for T = 1 as a function of asymptotic momentum for different values of µ (a) (c = 1) and multiple values of c (b) (µ = 1). One can see that for all values of µ and c, the density is everywhere almost zero, except for k ∈ [−5, 5]. This makes our approximation of the integral boundaries reasonable.

(25)

(a) One can see that the energy goes to zero when the temperature goes to zero and µ = 0. This agrees with our defini-tion that the ground state of the system is empty. When µ 6= 0, it becomes advan-tageous to produce particles and the en-ergy increases. When the temperatures rise, more and more particles are being created and the energy rises.

(b) When the interaction strength in-creases, the creation of particles becomes less favoured. So with increasing c, the energy decreases.

Figure 5: The energy for the Lieb-Liniger model as a function of temperature for different values of µ (a) (c = 1) and multiple values of c (b) (µ = 1). We set the ground state energy to zero.

(26)

(a) (b)

Figure 6: The effective velocity for the Lieb-Liniger model as a function of asymptotic momentum for T = 1 and different values of µ (a) (c = 1) and multiple values of c (b) (µ = 1).

(a) (b)

Figure 7: The effective velocity for the Lieb-Liniger model for k = 5 as a function of temperature for different values of µ (a) (c = 1) and multiple values of c (b) (µ = 1).

(27)

4 ISING MODEL

4 Ising Model

Now we will study the Ising model in a transverse field (18). This system has following Hamiltonian

H = −J

N

X

i

(gˆσx_i + ˆσ_izσˆz_i+1), (23) where J is an exchange constant to set the energy scale, g is the coupling with the transverse field and σ are the usual pauli matrices. It represents a one-dimensional lattice of spins that interact with both the transverse field and their nearest neighbours.

This model has two regimes depending on the coupling constant g. If g ≈ 0, then H ≈ −J P_iσˆz

iσˆi+1z and the ground state is twofold degenerate:

|↑↑ ... ↑i or |↓↓ ... ↓i . If g is very large, H ≈ −J Pigˆσ

x

i and the ground state becomes

non-degenerate

|→i⊗N =h√1 2

|↑i + |↓ii⊗N.

Since the number of ground states is an integer number, there must be a value of g for which the number of ground states goes uncontinuously from two to one. Indeed, there turns out to be a phase transition at g = 1. (10) Here we will focus on the regime g > 1.

4.1 Energy and density

From the calculation done in appendix A, we know that the Ising model has free fermionic excitations with following dispersion relation

Υ(k) = 2Jp1 + g2_{− 2g cos(k).}

Since these free fermions do not interact, the S-matrix becomes very simple and is just -1 for all values of the asymptotic momenta. This means that the derivative of the phaseshift is zero, which makes all equations of the Bethe ansatz much simpler.

With this in mind, the dressed energy becomes for µ = 0 (k) = Υ(k).

(28)

4.1 Energy and density 4 ISING MODEL

Indeed, since the particles don’t interact, the rest of the system does not need to adjust to a particle being created and the dressed energy just becomes the dispersion relation.

The density is then

ρ(k) = 1

eΥ(k)/T _{+ 1}.

This is indeed what we expect for free fermions.

Figure 8: The density of quasiparticles in the Ising model (with J = 1/4, g = 4 and T = 1) as a function of asymptotic momenta (left) and the total density as a funtion of temperature (right).

The total energy is given by E = Z dk 2πρ(k)Υ(k) = Z dk 2π Υ(k) eΥ(k)/T _{+ 1}

and is plotted in figure 9, where it is compared to some simulation based on the free fermion approach, which is exact here. One can see that the results overlap quite nicely and the deviation we do see is coming from finite precision effects in the code.

(29)

4.2 Lightcone dynamics 4 ISING MODEL

Figure 9: The energy of the Ising model (with J = 1/4 and g = 4) as a function of temperature (left), compared with the exact result. These results overlap of course, since the particles are in fact really free fermions, therefore we plotted the difference between these two functions in the right plot. The deviation we do see is coming from finite precision effects in the code.

4.2 Lightcone dynamics

Since the phase shift is zero, the dressing operator defined in equation 2.2.2 becomes the trivial operator

hdr(k) = h(k). Because of this, the effective velocity becomes

veff=

E0dr(k) p0dr_(k) =

dΥ(k) dk

and is thus independent of the temperature. Indeed, the existence of other particles does not matter here, since they do not interact anyway. The effec-tive velocity of the Ising model is plotted in figure 10. When a local quench is applied to the system, it will be out of equilibrium at that point. A lot of quasiparticles will be created there and will transport the energy density away with this effective velocity. Since this is a finite velocity, this results in a lightcone-type behaviour. (1; 2) These lightcones have already been experimentally observed. (14)

In figure 11, we plot the obtained maximum effective velocity on a sim-ulation of this lightcone behaviour based on the MPS methodology (see Ref (11)). We see that the boundaries of this lightcone is indeed determined by the maximal velocity of the quasiparticles.

(30)

4.2 Lightcone dynamics 4 ISING MODEL

Figure 10: The effective velocity of the Ising model (with J = 1/4, g = 4 and T = 1) as a function of asymptotic momenta. Since the particles do not interact, this figure is valid at all temperatures.

Figure 11: This is a plot of the time evolution of the energy density in the Ising model (with J = 1/4 and g = 4) in the ground state (left) and at T = 2 (right). The white lines are determined by x = vefft, where veff is the

maximal effective velocity. We can see that the boundaries of the lightcone corresponds with this maximal effective velocity. The points on the time-axis correspond to the time steps in the simulation (0.005 s).

(31)

5 HEISENBERG SPIN LADDER

5 Heisenberg spin ladder

Up till now, we used the thermodynamic Bethe ansatz on integrable models and checked the exact solution. Now we will use it to solve a non-integrable model, namely the spin-1/2 Heisenberg antiferromagnetic ladder. As was ex-plained in the first section of this dissertation, the Bethe ansatz can only be applied to systems that satisfy the Yang-Baxter equation i.e. integrable sys-tems. Then the real three body interactions could be written as a sequence of two body interactions. It might thus seem strange that we would like to apply this ansatz to non-integrable systems. The idea behind this is that, since we defined the ground state of the system to be empty, we expect that for low temperatures, there will not be enough particles for the three body interactions to be important (Remember that it were the multi-body interac-tions that violated integrability). The system would thus be approximately integrable for low temperatures. This assumption allows us to use the ther-modynamic Bethe ansatz for low temperatures in non-integrable systems. Since the thermodynamic Bethe ansatz does take the two body interactions into account, we expect this approximation to give better results than other approximations used to solve these systems, such as the free fermion approx-imation. We of course expect that this approximation will deviate from the exact solution when the temperature rises, since then the approximation is no longer valid.

5.1 Spin ladder without magnetic field

We can now begin by defining the Hamiltonian of the spin 1/2 Heisenberg antiferromagnetic ladder as follows (23; 24)

H =X i,l ~ Si,l· ~Si+1,l+ γ X i ~ Si,1· ~Si,2− h X i,l S_i,lz , (24)

where l = 1, 2 are the two legs of the ladder and Si,lis the spin operator on the

i’th site of the l’th leg of the ladder (see figure 12). The first term represents the interaction between nearest neighbours along a leg of the ladder, the second term is the nearest neighbour interaction between two spins on the same rung of the ladder. The last term is the coupling of all spins to a field. For now, we set h = 0.

(32)

5.1 Spin ladder without magnetic field 5 HEISENBERG SPIN LADDER

Figure 12: Graphical representation of the Hamiltonian of the spin 1/2 Heisenberg antiferromagnetic ladder. Every point is a site on which a spin 1/2 fermion lives. It can interact with its nearest neighbours on the same rung (J⊥) and along the same leg (Jk). In equation 24, we have set Jk = 1

and γ = J⊥/Jk.

The low energy spectrum can be intuitively derived by considering fol-lowing limits: (23)

• γ → ∞

In this limit, the interaction of nearest neighbours along the legs of the ladder can be ignored. The system now becomes a set of independent couples, each consisting of two spins. In the ground state, these two fermions are in the spin 0 state. As usual, this ground state is defined as the vacuum of the system. Excitations (particles) upon this ground state are created when one of these couples forms a spin 1 state. These excitations are called magnons and have spin 1.

• γ → 0

In this limit, there is no nearest neighbour interaction along the rungs and the system splits into two separate Heisenberg antiferromagnetic chains. The excited states of these chains are spinons (spin 1/2). When the chains are coupled, these spinons confine into the magnons.

Since these magnons have spin 1, there are three types of particles (the three components of the triplet). They all have the same dispersion relation, plotted in figure 13. Their interaction partner can also be one of these three

(33)

Figure 13: The dispersion relation of all three types of magnons in the spin 1/2 Heisenberg antiferromagnetic ladder for γ = 4.

types. This means that we will have a 9 × 9 S-matrix. The SU(2) symmetry of the system will make sure that this S-matrix is constant within each sector: (23; 24) S =   −eiθ0 _{× 1} 1×1 −eiθ1 × 1 3×3 −eiθ2 _{× 1} 5×5  ,

where the θ’s are the phaseshifts of every sector, plotted in figure 14. One can now easily perform the thermodynamic Bethe ansatz on every sector separately. It assumes that all particles scatter in that specific sector. The real system however, consists of particles that can have the three possible O(3) spin labels (−1, 0, 1). The sector in which two particles scatter depends on the combination of both of these spin labels. So a single particle is not a member of a certain sector, since the scattering sector depends on the other particle. So summing all contributions of the three sectors separately is not the way to go. So we will need to find another way to go about this.

5.1.1 Auxiliary particles

One attempt at solving this problem is to perform a derivation similar to Ref (27). There they consider the following S-matrix:

(34)

Figure 14: The phaseshift θ(k1−k2)of every spinsector of the ladder model for

γ = 4 as a function of the asymptotic momenta of both scattering particles. where here a, b, c and d are either − or + and

σT(λ) = λ Sv(λ) λ − iπ σR = −iπ Sv(λ) λ − iπ ; Sv(λ) = Γ 1 2 + λ 2iπ Γ − λ 2iπ Γ 1 2 − λ 2iπ Γ λ 2iπ . The S-matrix has thus the following form

S =     σT + σR σT σR σR σT σT + σR     .

According to the discussion in section 2, the main goal is to diagonalise the trace of the transfer matrix. The way to do this, is by introducing additional degrees of freedom, so called auxiliary or virtual particles. (16; 27) Since every particle has a rapidity, this addition comes down to adding extra rapidities.

From section 2, we know the transfer matrix Ta0

a to be defined as Ta0 a (u) ≡ T a0 a (u|λ1, ..., λN)cb11,...,c,...,bNN = S a1c1 ab1 (u−λ1)S a2c2 a1b2(u−λ2)...S a0cN aN −1bN(u−λN),

(35)

where a summation over all repeated indices is implied.

Since the S-matrices should obey the Yang-Baxter equation 7: X a00_,b00 Ta00 a (u)Tb 00 b (v)Sa 0_b0 a00_b00(u − v) = X a00_,b00 S_aba00b00(u − v)T_bb000(v)Ta 0 a00(u).

One can now start from the ’bare vacuum’ state |0i, where all the spins of the real particles are down |− − ...−i. One can see that this is an eigenstate of T = T+ + + T − −: T+ +(u) |0i = S a1c1 +− (u − λ1)Saa12−c2(u − λ2)...S +cN aN −1−(u − λN) |c1, ...cNi ,

where we used that in |0i, all spins are down: b1 = −, ..., bN = −. Remember

that there is an implicit sum over all repeated indices. This means that aN −1 = ±. We see however that if aN −1 = −, S−−+cn(u − λN) = 0 for both

values of cN due to the δ functions in equation 25. The only contribution

different from zero will be the one with aN −1= + and cN = −.

The second to last term then becomes S+cN −1

aN −2−. This term is again only

different from zero if aN −2 = + and cN −1 = −. One can continue this

reasoning until all S-matrices are fixed: T+ +(u) |0i = S +− +−(u − λ1)S+−+−(u − λ2)...S+−+−(u − λN) |− − ...−i = N Y i σT(u − λi) |0i .

A similar derivation can be followed for T− −:

T₋−(u) |0i = Sa1c1

−− (u − λ1)Saa12−c2(u − λ2)...S

−cN

aN −1−(u − λN) |c1, ...cNi .

Here, the first S-matrix is only different from zero when a1 = c1 = −. This

fixes all S-matrices to the following form

T₋−(u) |0i = S₋₋−−(u − λ1)S−−−−(u − λ2)...S−−−−(u − λN) |− − ...−i

= N Y i σT(u − λi) + σR(u − λi) |0i = N Y i Sv(u − λi) |0i .

(36)

Adding these two contributions results in

T (u) |0i =T₊+(u) + T₋−(u)|0i =

N Y i 2σT(u − λi) + σR(u − λi) |0i T (u) |0i = N Y i σT(u − λi) + Sv(u − λi) |0i , so the bare vacuum state is indeed an eigenstate of T (u).

Now the space of states can be constructed by applying T−

+(v) to this

bare vacuum. A generic state is then given by

|v1, v2, ..., vMi = T+−(v1)T+−(v2)...T+−(vM) |0i .

This comes down to adding M auxiliary particles to the bare vacuum. In general, such a state is not an eigenvector of T (u).

From equation 8, one can see that this state does not depend on the sequencing of these transfer matrices. Indeed, equation 8 with a = b = + and a0 _{= b}0 _{= −} _is X a00_,b00 Ta00 + (u)Tb 00 + (v)S −− a00_b00 = X a00_,b00 S₊₊a00b00T_b−00(v)T_a−00(u).

The only contribution that does not return zero is T₊−(u)T₊−(v)S₋₋−− = S₊₊++T₊−(v)T₊−(u)

[T₊−(u), T₊−(v)] = 0. And in a similar way

[T₋+(u), T₋+(v)] = 0.

Before continuing, let us first take a look at some other commutation relations. We obtain one of them by looking at equation 8 with a = +, b = +, a0 = + and b0 = −: X a00_,b00 Ta00 + (u)Tb 00 + (v)S +− a00_b00(u − v) = X a00_,b00 S₊₊a00b00(u − v)T_b−00(v)T_b+00(u) T+ +(u)T − +(v)S +− +−(u − v) + T − +(u)T + +(v)S +− −+(u − v) = S++++(u − v)T − +(v)T + +(u)

(u − v)T₊+(u)T₊−(v) = iπT₊−(u)T₊+(v) + (u − v − iπ)T₊−(v)T₊+(u). And in a similar way:

(37)

One can now use these commutation relations to determine a constraint when a general state is an eigenvector of T (u) = T+

We can now use the fact that |0i is an eigenstate of T+

+ to write T+ +(u) |v1v2...vMi = M X j=1 −iπ vj − u M Y i6=j vi− vj + iπ vi− vj N Y k=1 σT(vj − λk) |v1, ..., vj, ..., vM, ui + M Y j=1 iπ + vj − u vj − u N Y k=1 σT(u − λk) |v1, ..., vMi ,

where v1, ..., vj, ..., vM is v1, ..., vM where vj is removed.

In a similar way, we find

T− −(u) |v1v2...vMi = M X j=1 −iπ u − vj M Y i6=j vj− vi+ iπ vj− vi N Y k=1 Sv(vj − λk) |v1, ..., vj, ..., vM, ui + M Y j=1 iπ + u − vj u − vj N Y k=1 Sv(u − λk) |v1, ..., vMi .

The condition that |v1, ..., vMi is an eigenstate of T (u) is then that the first

term of both previous equations sum to zero.

M X j=1 −iπ vj− u M Y i6=j 1 vi− vj " (vi− vj + iπ) N Y k=1 σT(vj − λk) − (vi − vj − iπ) N Y k=1 Sv(vj− λk) # = 0 M Y i6=j vi− vj− iπ vi− vj + iπ N Y k=1 σT(vj− λk) Sv(vj − λk) = 1.

(38)

We have now found a condition under which a general state is an eigenvector of T (u). We are thus now able to find all the eigenstates and thus diagonalise T (u).

5.1.2 Application to spin 1/2 Heisenberg antiferromagnetic lad-der

We now try to perform a similar derivation, but for spin-1 excitations. First of all, we will need the S-matrix in the same form as equation 25. Since the low-energy properties of gapped systems are well described by the O(3) nonlinear sigma model (6), one expects the S-matrix to have the following form.

S_abcd(λ) = δabδcdσ1(λ) + δacδbdσ2(λ) + δadδbcσ3(λ),

where the σi are yet to be determined.

Now that we have the S-matrix of the spin 1/2 Heisenberg antiferromag-netic ladder in the same form as equation 25, we can try and perform a similar derivation. Remember that here a, b, c, d ∈ [−, 0, +], since we have spin-1 excitations.

We again consider a state of N particles with rapidities λ1, λ2, ..., λN. One

can then introduce the transfer matrix Ta0

a as Ta0 a (u) ≡ T a0 a (u|λ1, ..., λN)cb11,...,c,...,bNN = S a1c1 ab1 (u−λ1)S a2c2 a1b2(u−λ2)...S a0cN aN −1bN(u−λN),

where a summation over all repeated indices is implied.

Since we assume that this system approximately satisfies the Yang-Baxter equation, we can again write down equation 8 and use it to calculate some commutation relations. Inserting a = b = a0 _{= b}0 _{= +} _gives

T+ +(u)T++(v)S++++(u − v) + T − +(u)T − +(v)S−−++(u − v) + T+0(u)T+0(v)S00++(u − v)

= S₊₊++(u − v)T₊+(v)T₊+(u) + S₊₊00 (u − v)T₀+(v)T₀+(u) + S₊₊−−(u − v)T₋+(v)T₋+(u). Every other term will be zero due to the dirac δ fuctions in the S-matrix. Do-ing the same for for some well chosen combinations of a, b, c, d and summDo-ing everything results in T (u)T (v) =T+ +(u) + T00(u) + T − − (u) T+ + (v) + T00(v) + T − −(v) = T (v)T (u) [T (u), T (v)] = 0.

This is exactly what we had before. In the next step however, we are going to encounter some problems. The next step was to start from the bare vacuum

(39)

state |0i and check if it is an eigenstate of T (u). The problem is that the ground state consists of all the spin pairs being in the singlet state. This singlet state however cannot be represented by one of the O(3) spin labels (+, 0, −). This way the transfer matrix cannot be constructed, since the S-matrix requires a b-label that is one of these O(3) spin labels. This represents the fact that the transfer matrix scatters a particle from all other particles in this bare vacuum state. But there are no particles in this vacuum state.

Using another bare vacuum state doesn’t solve this. Indeed, a vacuum state of the form |0i = |− − ...−i will turn out to be not an eigenstate of Ta

a(u). To see this, we look at

T₊+(u) |0i = Sa1c1

+− (u − λ1)Saa12−c2(u − λ2)...S

+cN

aN −1−(u − λN) |c1, ...cNi .

Due to the extra dirac δ function in the S-matrix compared to the one in the previous section, this is not enough to fix all S-matrices in this equation. Indeed, both the first and the last S-matrix have two nonzero contributions. Writing the sum over a1 explicitly:

T+ +(u) |0i = S +− +−(u − λ1)S+−a2c2(u − λ2)...Sa+cN −1N −(u − λN) |−, c2, ...cNi + S₊₋−+(u − λ1)S−−a2c2(u − λ2)...Sa+cN −1N −(u − λN) |+, c2, ...cNi . Here, Sa2c2

+− (u−λ2)has two nonzero contributions and S−−a2c2(u−λ2)has three.

Now writing the sum over a2 explicitly:

T+ +(u) |0i = S +− +−(u − λ1)S+−+−(u − λ2)S+−a3c3(u − λ3)...Sa+cN −1N −(u − λN) |−−, c3, ...cNi + S₊₋+−(u − λ1)S+−−+(u − λ2)S−−a3c3(u − λ3)...Sa+cN −1N −(u − λN) |−+, c3, ...cNi + S₊₋−+(u − λ1)S−−++(u − λ2)S+−a3c3(u − λ3)...Sa+cN −1N −(u − λN) |++, c3, ...cNi + S₊₋−+(u − λ1)S−−00 (u − λ2)S0−a3c3(u − λ3)...Sa+cN −1N −(u − λN) |+0, c3, ...cNi + S₊₋−+(u − λ1)S−−−−(u − λ2)S−−a3c3(u − λ3)...Sa+cN −1N −(u − λN) |+−, c3, ...cNi .

From this it is already clear that |0i will not be an eigenstate of T+ + (u).

But maybe we can find some conditions under which this does become an eigenstate. For that, we look at T−

−(u): T₋−(u) |0i = Sa1c1 −− (u − λ1)Saa12−c2(u − λ2)...S −cN aN −1−(u − λN) |c1, ...cNi = S₋₋++(u − λ1)S+−a2c2(u − λ2)...S −cN aN −1−(u − λN) |+, c2, ...cNi + S₋₋00 (u − λ1)S0−a2c2(u − λ2)...Sa−cN −1N −(u − λN) |0, c2, ...cNi + S₋₋−−(u − λ1)S−−a2c2(u − λ2)...S −cN aN −1−(u − λN) |−, c2, ...cNi .

Every state coming from T+

+(u) |0i that is different from |0i should sum to

zero with its counterpart coming from T−

−(u) |0i. So there will be a condition

(40)

5.2 Spin ladder approximation 5 HEISENBERG SPIN LADDER

Spin combination Spin sector possibilities

↑↑ spin 2 ↑ 0 1 2 spin 2 + 1 2 spin 1 ↑↓ 1 6 spin 2 + 1 2 pin 1 + 1 3 spin 0 0 ↑ 1₂ spin 2 +1₂ spin 1 00 2₃ spin 2 +1₃ spin 0 0 ↓ 1₂ spin 2 +1₂ spin 1 ↓↑ 1 6 spin 2 + 1 2 spin 1 + 1 3 spin 0 ↓ 0 1 2 spin 2 + 1 2 spin 1 ↓↓ spin 2

Table 1: Spin sector possibilities for all combinations of two spin 1 particles

One of these conditions is remarkable, indeed T+

+(u) |0idoes not give rise

to a state with c1 = 0, while such a state does arise from T−−(u) |0i. This

means we need following condition: X a2,...,aN −1 S₋₋00 (u − λ1)S0−a2c2(u − λ2)...S −cN aN−(u − λN) |0, c2, ...cNi = 0, or S₋₋00 (u − λ1) = 0.

But this condition is not fulfilled.

It would thus seem that it is highly nontrivial to generalise the method used in Ref (27).

5.2 Spin ladder approximation

Since the previous attempt at applying the thermodynamic Bethe ansatz to the ladder model failed, we try the following approximation instead. Remem-ber that the first approximation was that at low temperatures, the ladder model is approximately integrable. So applying the Bethe ansatz is already an approximation in and of inself.

the main additional approximation here is that all particles scatter in a single auxiliary sector. The phaseshift of this auxiliary sector is then a linear combination of the ones plotted in figure 14. In order to determine this linear combination, we take a look at the Glebsch-Gordan coefficients. (9) With these coefficients, one can determine the spin sector possibilities in which two spins will scatter. The results are shown in table 1. Since the dispersion relation is the same for ↑-, 0- and ↓- particles and they will scatter

(41)

with the same phaseshift, we assume that, on average, they are created in equal amount.

When a particle is created, there is thus a 1/3 chance of it being a ↑-particle. When this particle moves around and scatters from another particle, there is again a 1/3 chance that the other particle is also a ↑-particle and the spin combination will be ↑↑. This way, each spin combination has an equal 1/9 chance of occurring. Each spin combination has some chances of occurring in a certain spin sector (see table 1). So has the ↑↓ combination a 1/6 chance of occurring in the spin 2 sector, a 1/2 chance of occurring in the spin 1 sector, and so on. This way the average phase shift can be obtained by θav(k1− k2) = 5 9θ2(k1− k2) + 3 9θ1(k1− k2) + 1 9θ1(k1 − k2) (26) Before moving on, we realize that the thermodynamic Bethe ansatz changes a bit when there are three kinds of quasiparticles. Indeed, equation 11 becomes

ζ(k) = ρ1(k) + ρ0(k) + ρ−1(k) + ρh(k),

where ρ1 is the density of particles with the + spinlabel, ρ0 is the density of

the particles with the 0 spinlabel, etc.

Since we assumed that each kind of quasiparticle is created in equal amounts, one can write

ζ(k) = 3ρ0(k) + ρh(k).

Redoing the derivation following equation 11, one can see that the equation for the dressed energy 16 remains the same, if ρh(k)/ρ0(k) = exp{(k)/T }

and θ0_{(k − k}0_{) = θ}0

av(k − k

0₎_{. The formula for the density changes however.}

ζ(k) = 3ρ0(k) + ρh(k) ζ(k) ρ0_(k) = 3 + ρh(k) ρ0_(k) = 3 + e (k)/T ρ0(k) = ζ(k) 3 + e(k)/T.

Remember that these quasiparticles were the result of two spins entangling to a triplet state. Since there are three ways of forming this triplet, we have three kinds of quasiparticles. It is however impossible for two quasiparticles to be at the same place, since that would require these two spins to be entangled in such a way that they result in two types of triplets. This is impossible. As a result of this, ζ(k) = 1.

(42)

The total density, given by

ρ(k) = ρ1(k) + ρ0(k) + ρ−1(k) = 3ρ0(k),

then becomes

ρ(k) = 3ρ0(k) = 3 3 + e(k)/T.

All quantities, such as the total energy can then be calculated with the usual formulas if the density is changed to the equation above.

When applying the thermodynamic Bethe ansatz to this system, it turns out that the range of asymptotic momenta for which we have the phase shift (plotted in figure 14) is not large enough to capture the whole peak in the density. Therefore we have first extrapolated these phaseshifts to a larger range of asymptotic momenta. We did this by fitting a first order polynomial a + bk1 to the phaseshift at a fixed value of k2. For the spin-1 sector, we fitted

two linear functions, one for the left side of the discontinuity and one for the right side. Afterwards we looked how these fitted parameters a and b changed when k2 was changed. We again fitted a first order polynomial to this and

used this to extrapolate to a larger range of asymptotic momenta. The result can be seen in figure 15.

Now one can apply the thermodynamic Bethe ansatz to this ladder model. the results are plotted in figures 16, 17 and 18. In figure 18, we compare our results to a simulation based on the MPS methodology. (24) It is surprising that the Bethe ansatz produces a worse result than the free fermions. We would expect it to perform better since it takes interactions into account. This makes us assume that the approximation of the average phaseshift is not a good one.

(43)

Figure 15: The extrapolation of the phaseshift for the spin-0 sector (top) and the spin-1 sector (middle) and the spin-2 sector (bottom) of the ladder model.

(44)

Figure 16: The density for the ladder model without magnetic field and for γ = 4 as a function of asymptotic momenta (left) and the total density as a function of temperature (right). On the left plot, it is clear that the initial range of asymptotic momenta ([3.05, 3.25]) was not large enough. The extrapolated phaseshift however does capture the whole peak of excitations.

Figure 17: The effective velocity for the ladder model without magnetic field and for γ = 4 as a function of asymptotic momenta (T = 2) (left) and function of temperature (k = 1) (right).

(45)

5.3 Spin ladder with magnetic field 5 HEISENBERG SPIN LADDER

Figure 18: Plot of the total energy of the spin-1/2 Heisenberg antiferromag-netic ladder model without magantiferromag-netic field and for γ = 4 as a function of temperature. It is surprising that the Bethe ansatz actually gives a worse result than the free fermion method.

5.3 Spin ladder with magnetic field

Since the previous approximation turned out to be a bad one, we try some-thing else: we apply a magnetic field (h = 3). A consequence of this magnetic field is a shift in the dispersion relation. Before, every particle had the same dispersion relation, but that degeneracy gets broken by the magnetic field. The dispersion relation of the particles with spin +1 gets lowered by 3, for spin 0 particles, it remains the same and for spin −1 particles, it gets raised by 3. This ensures that for the spin +1 particles, the gap to the ground state is much smaller than for the other particles due to alignment to the magnetic field. This way, approximately all particles will have spin +1. Since two spin +1 particles will always scatter in the spin 2 sector, one can safely ignore the other two sectors.

(46)

5.3.1 Thermodynamic Bethe ansatz

Now we can apply the thermodynamic Bethe ansatz. The resuls are given in figures 19, 20, 21 and 22. In figure 19, a comparison is made between the Bethe ansatz and the free fermion approach. They are both compared to the results of a simulation based on MPS. (24) It is clear that both methods are very similar to the results given by the simulation for β > 1.5. For higher temperatures, both methods deviate. Although the Bethe ansatz approach lies closer to the results of the simulation than the free fermion method, the improvement is almost negligible. This can be explained however by the fact that the interaction between particles is not that high in the spin-2 sector. Indeed, looking at figure 14, the phaseshift in the spin-2 sector is of the order 0.01 and is thus very close to zero. Of course, when the phaseshift is zero in the Bethe ansatz, the equations become those of the free fermions.

It is clear that for large temperatures the Bethe ansatz result deviates largely from the one obtained with simulations. This is of course expected since then the assumption that the tree-body interactions are negligible is not reasonable anymore. here, the effects of non-integrability come into play. A similar plot of the lightcone dynamics to the one of the Ising model is given in figure 23. We again observe an overlap between the line correspond-ing to the maximal effective velocity and the boundaries of the lightcone.

Figure 21: The density for the ladder model with magnetic field (h = 3) and γ = 4 as a function of asymptotic momenta (left) and the total density as a function of temperature (right).

(47)

Figure 19: Comparison of the total energy for the ladder model with magnetic field (h = 3) and γ = 4 as a function of temperature of the Bethe ansatz approach to the one obtained with a simulation. We also included the result of the free fermion method as a comparison. It turns out that the free fermion approach and the Bethe ansatz give very similar results. This is because in the phaseshift in the spin 2-sector is of the order 0.01 and is thus similar to 0.

(48)

Figure 20: Comparison of the total energy for the ladder model with magnetic field (h = 3) and γ = 4 as a function of temperature of the Bethe ansatz approach to the one obtained with a simulation, but now zoomed in.

Figure 22: The effective velocity for the ladder model with magnetic field (h = 3) and γ = 4 as a function of asymptotic momenta (T = 2) (left) and function of temperature (k = 1) (right).

(49)

Figure 23: This is a plot of the time evolution of the energy density in the ladder model with magnetic field (h = 3) and γ = 4 in the ground state (left) and at T = 2 (right). The white lines are determined by x = vefft, where

veff is the maximal effective velocity. We can see that the boundaries of the

lightcone corresponds with this maximal effective velocity. The points on the time-axis correspond to the time steps in the simulation (0.005 s).

(50)

6 CONCLUSION

6 Conclusion

In this dissertation, we showed that the thermodynamic Bethe ansatz is a method that can be used effectiveley to solve integrable quantum systems. Indeed, both the total energy of the system and the effective velocity pro-duced by this method are very close to what one expects and it can explain the appearance of the observed lightcone-like dynamics.

This method can also be used to solve non-integrable systems in the low temperature limit, since there the three-body interactions and hence the non-integrability remains negligible. The advantage over other methods such as the free fermion method remains inconclusive, since we were only able to perform the Bethe ansatz on non-integrable models with a very low phaseshift. Indeed, it turns out to be nontrivial to generalise the method outlined in section 5.1 and we were forced to work with a magnetic field. We however remain convinced that the thermodynamic Bethe ansatz will produce much better results for systems with a higher phaseshift.

(51)

A SUMMARY

A

Summary

A.1 Summary

Methods of solving nonintegrable systems remain low in number. In this dissertation, we show that the thermodynamic Bethe ansatz can be used to solve the nonintegrable spin-1/2 Heisenberg antiferromagnetic ladder in the low temperature limit. Normally, the Bethe ansatz can only be used for integrable models, i.e. models which satisfy the Yang-Baxter equation. This makes sure that the multiple-body interactions can be written into a series of two-body interactions rendering the wavefunction in the Bethe form. This is not the case for non-integrable models, exept in the low temperature limit. Indeed, in this limit, the density of quasiparticles will be low enough to ignore the three-particle interactions, rendering the model approximately integrable. This means that in the low temperature limit, one can approximately use the thermodynamic Bethe ansatz.

A.2 Samenvatting

Tot op de dag van vandaag blijft het aantal methodes om niet-integreerbare systemen op te lossen beperkt. In dit proefschrift tonen we aan dat de ther-modynamische Bethe ansatz gebruikt kan worden om de niet-integreerbare spin-1/2 Heisenberg antiferromagnetische ladder op te lossen in de lage tem-peratuurslimiet. Normaal gezien kan de Bethe ansatz enkel gebruikt worden in integreerbare modellen, i.e. modellen die aan de Yang-Baxter vergelijking voldoen. Dit zorgt ervoor dat de meerdere-lichaamsinteracties geschreven kunnen worden als een reeks van tweelichaamsinteracties, wat ervoor zorgt dat de golffunctie in de Bethe vorm geschreven kan worden. Deze redenering gaat niet op voor niet-integreerbare systemen, behalve in de lage temperatu-urslimiet. Inderdaad, in deze limiet, zal de densiteit van quasideeltjes klein genoeg zijn om de drielichameninteracties te verwaarlozen. Dit maakt het model benaderend integreerbaar, wat betekent dat men, in de lage temper-atuurslimiet, de thermodynamische Bethe ansatz benaderend kan gebruiken.