Continuous matrix product states for one-dimensional optical lattices

(1)

Faculty of Science, Department of Physics and Astronomy

Academic year 2019–2020

Continuous matrix product states

for one-dimensional optical lattices

Toon Pillaert

Promotor: Prof. dr.Jutho Haegeman Supervisor: Benoˆıt Tuybens

Thesis submitted in partial fulfilment of the requirements for the degree of Master of Science in Physics and Astronomy

(2)

(3)

Admission for Usage

The author and promoter give the permission to use this thesis for consultation and to copy parts of it for personal use. Every other use is subject to the copyright laws, more specifically the source must be extensively specified when using from this thesis.

Toon Pillaert, June 2020

(4)

(5)

Acknowledgements

First, I would like to thank my promotor professor Jutho Haegeman and supervisor Benoˆıt Tuybens for their support and guidance. Professor Haegeman in particular for putting time and effort in this thesis despite the hectic and unusual times of the last months and Benoˆıt because I could always turn to him for feedback and questions.

I also would very much like to thank my teachers at the conservatory of Ghent who had to cope with my occasionally but sometimes unavoidable absent. They were not only understanding but often took interest in what I was doing besides writing music.

It was quite a disappointment when it turned out that the calculations took multiple days of runtime on my laptop, even for low bond dimensions. I’m therefore grateful that the Pileus computer, owned by UGent, could be used to run my simulations. If not, the data in this thesis would have been very sparse.

My parents can of course not be forgotten in this list. Their interest and my occasional visits to them were always something to look forward to.

Finally, Nathalie and Alma, you made my days in quarantine writing this thesis a lot more enjoyable.

Toon Pillaert Ghent, June 2020

(6)

(7)

Samenvatting

Het doel van deze thesis was het toepassen van de continuous matrix product state ansatz, kortweg de cMPS ansatz, op een één dimensionaal systeem van hardcore bosonen met een peri-odieke potentiaal in de thermodynamische limiet. Toen de cMPS ansatz werd ge¨ıntroduceerd in 2010 als de continue limiet van een specifieke klasse van matrix product states, of MPS, werd ze enkel gebruikt in het geval van uniforme systemen. Hier wordt dit uitgebreid naar peri-odieke systemen door middel van het gebruik van Fourierreeksen voor de relevante grootheden. Daar de bosonen een oneindig grote contact interactie hebben, i.e. hardcore bosonen zijn, zullen ze verschillende eigenschappen krijgen die overeenkomen met die van een fermion gas. In de afwezigheid van een externe potentiaal heeft de Hamiltoniaan van het systeem een U(1) symmetrie die overeenkomt met het behoud van het aantal deeltjes. Men zou verwachten dat deze symmetrie spontaan gebroken wordt en de grondtoestand een Bose-Einstein condensaat is. Omdat het systeem echter één dimensionaal is, verhindert het Mermin-Wagner theorema dit. Het resultaat is dat de grondtoestand zich in een toestand bevindt die gekenmerkt wordt door zogenaamde quasi-long-range order, waarbij de correlatiefuncties asymptotisch een algebra¨ısch verval kennen. Wanneer men een periodieke potentiaal aanlegt (hier met de eigenschap na = 1 waar bij n de gemiddelde deeltjesdichtheid is en a de periodiciteit van de potentiaal) zal er een band gap ontstaan wat gepaard gaat met een exponentieel verval van correlatiefuncties. De bosonische cMPS benadering zal vergeleken worden met het fermion geval en de local density approximation, of LDA. Vervolgens zal ook gesimuleerd worden in welke fase het systeem zich bevindt wanneer de contactinteractie eindige waarden aanneemt in aanwezigheid van een potentiaal.

(8)

(9)

Summary

The purpose of this thesis was to apply the continuous matrix product state ansatz, or cMPS, to a one-dimensional system of hardcore bosons with a periodic potential in the thermody-namic limit. When the cMPS ansatz was introduced in 2010 as the continuous limit of a specific class of matrix product states or MPS, it was only used in the case of uniform sys-tems. Here this is extended to periodic systems through the use of Fourier series for the relevant quantities. Since the bosons have an infinitely large contact interaction, i.e. are hardcore bosons, they will have properties similar to those of a fermion gas. In the absence of an external potential, the Hamiltonian of the system has a U(1) symmetry corresponding to particle number conservation. One would expect this symmetry to be broken spontaneously by the ground state, forming a Bose-Einstein condensate. However, because the system is one dimensional, the Mermin-Wagner theorem prevents this. The result is that the ground state is in a state characterized by so-called quasi-long-range order, in which the correlation functions asymptotically have an algebraic decay. When a commensurate periodic potential (with the property na = 1 where at n is the average particle density and a is the periodicity of the potential) is applied to the system a band gap will occur, indicating an exponential decay of correlation functions. The bosonic cMPS approach will be compared to the fermionic case and the local density approximation, or LDA. Finally, it will be examined which phases occur when a potential is applied to a bosonic system subject to finite contact interactions.

(10)

Conventions

• Natural units are used where ~ = c = kB= 1.

• Curly brackets { · , · } denote anticommutators, square brackets [ · , · ] denote commuta-tors.

(11)

Introducing Tensor Networks

1.1 Why tensor networks?

Upon examining the existing literature on numerical methods for strongly correlated quantum systems, one soon starts to wonder about the wide variety of existing methods. One can for example exactly diagonalize the hamiltonian of the system but this is of course restricted to small systems and far away from the thermodynamic limit, where the subject of this thesis is situated. Another popular method is the mean field approximation, but because for every particle mean field theory averages out the interaction with all other particles quantum correlations are completely lost. Quantum Monte Carlo methods are a very popular approach but it is not applicable to every system because of the sign problem, though it is often a benchmark tool for the systems where it does. Tensor networks are not a perfect ansatz either but it extends the range of accessible systems considerably by making use of the entanglement structures of the concerned systems.

1.1.1 Introduction

Consider the generic many-body wave function of a system consisting of N particles |Ψi =X

{i}

Ci1i2···iN|i1i2· · · iNi (1.1)

where each |ini denotes the local basis of the particle denoted by n. Assume then that the

considered system consists of particles of the same species, each living in its local Hilbert space of dimension d, H = Cd. In order to describe the full phase space of this system, one then needs to keep track of the dN coefficients Ci1i2···iN. This exponential growth in the

number of coefficients with respect to the number of particles renders the exact calculations of macroscopic quantities virtually impossible, considering that one needs particle numbers of the order of the constant of Avogadro, 6.02214076 × 1023, for macroscopic systems. However, in the early nineties of the previous century Steven White, a former PhD student of Kenneth

(14)

Chapter 1. Introducing Tensor Networks 2

Wilson, developed the so-called density matrix renormalisation group method (DMRG), an algorithm constructed to find the ground state(s) of one-dimensional quantum systems. White had understood that most of phase space was irrelevant because the entanglement degrees of freedom were the relevant degrees of freedom of the wave function. It took some years before it was recognized that the DMRG algorithm was in essence a variational optimization algo-rithm over matrix product states, MPS (they will be described in detail in the next section) and that the results that came out of the algorithm were finitely correlated states1. Tensor networks, the generalisation of MPS, has thus provided an invaluable tool by reducing the phase space exponentially so that it now has become much easier to describe ground states and the entanglement structures of various quantum systems2 because the complexity growth now scales polynomialy in N .

Figure 1.1: Various examples of description of Ci1i2···1N. (a) is the visualisation of a matrix product

state (MPS) with periodic boundary conditions, (b) that of the two-dimensional exten-sion of MPS, namely PEPS and (c) that of the multi-scale entanglement renormalization ansatz (MERA).

One can always divide the system described by Eq. (1.1) in two arbitrary subsystems A and

1

Note however that critical phenomena usually fall outside this scope because they are often characterized by scale invariance and thus infinite correlation lengths.

2

It has has been proven[7] that MPS capture the whole subspace of states satisfying an area law if the maximal bond dimension is allowed to be

D ≥ max n min dimH_L(n) L , dimHL (n) R (1.2) Where HL(n) L and HL (n) R

are the left and right Hilbert spaces associated to a ’cut’ of the MPS at site n. In the case of systems with more spatial dimensions, it was proven that tensor networks only describe an exponential small subspace of this the already exponentially small subspace with the property that the correlation functions vanish. When this condition is tightened to exponentially decaying correlation functions it remains an open question if tensor networks completely cover the concerned subspace[5].

(15)

B. If the physical state is chosen at random from the dN-dimensional phase space, and one would calculate the Von Neumann entanglement entropy[2] between those two regions, i.e.

S = −Tr[ρAlnρA] (1.3)

it would scale with the volume of the smallest of these two subregions. Physicists have the incredible luck that nature has ensured that the physical relevant Hamiltonians all are lo-cal. For the ground states and low energy excitations, the entanglement entropy of these local Hamiltonians tends to scale only with the surface area of the boundary between the two subsystems[3], S v ∂A. The low-energy states of a local Hamiltonian are thus heavily constrained in the sense that they must obey the area law. It now happens that the states that satisfy this criterion cover only an exponentially small subspace of the unmanageably large many-body Hilbert space.

Figure 1.2: The small corner of the many-body Hilbert space where the area law applies. This is also the region where virtual all physical relevant states reside[1]

Furthermore, it has also been proven quite recently, that when a random state of H⊗N is evolved over time in a polynomial order of N , the manifold of states that can be reached is exponentially small[4]. Given that the physical ground states of local Hamiltonians already live in a tiny subspace of the vast Hilbert space of all possible states, this indicates that it is almost virtually impossible to leave this subspace of finitely-correlated states. Another consequence is that the vast majority of the many-body Hilbert space seems to be nonphysical as it can not be reached in a reasonable amount of time. This begs the question how physically relevant it is to describe physical states as vectors in the Hilbert space of all possible states[4].

(16)

1.2 Diagrammatic notation

A tensor, as far as this thesis is concerned, is a multidimensional array of in general complex numbers. They are usually denoted by a single capital letter which is labeled by indices i, j, ..., the number of which is the dimension of the tensor, called the rank or the order of the tensor. In contractions, the Einstein summation rule is used, meaning that repeated indices indicate a summation there-over.

Cµν =

X

λ

AµλBλν:= AµλBλν (1.4)

When doing calculations with these mathematical objects, one soon runs into problems be-cause of the large notational overhead since there are many indices and contractions. This issue however can be bypassed by using the diagrammatic notation developed by Roger Pen-rose in the early seventies, which provides a clear visual picture of TN computations.

The basic rules of the diagrammatic notation are very simple.

· A tensor is represented by some solid shape where every index is represented by an outgoing line.

· A connection between two index lines implies a contraction, or summation over the con-nected indices

Figure 1.3: Some basic examples of the diagrammatic notation

When using this framework in a physical context, a distinction will be made between two kinds of indices. Namely the physical ones, which are the external indices, and the virtual ones, which are those one contracts over when calculating observables. The dimension of these virtual indices is denoted by the bond dimension D. This is in general not related to the dimension of the physical indices (e.g. a spin 1 particle would have 3 dimensions). There exist several more rules for these schematic notations but these are often author dependent or mere conventions in order to simplify some ubiquitous special tensors such as the metric tensor or the Levi-Civita tensor. However, we will do just fine however with this basic variant.

(17)

1.3 MPS

1.3.1 Construction

Consider a one-dimensional system of N spins living on a lattice L. The individual spin states are represented by |si_i where i denotes the site. Pure states of this system are denoted as |s1s2· · · sNi and any arbitrary state can be described as a superposition of these pure states.

|Ψi =X

{s}

Cs1s2···sN|s1s2· · · sNi (1.5)

Depending on which Hamiltonian is used, the spins will interact with each other in one or the other way. Because this is a quantum system, entanglement is introduced in the system. The reason why matrix product states became a popular technique to tackle such systems is because they are suited to produce long range correlations, even when the system has only nearest neighbour interactions.

In order to introduce entanglement in an efficient manner, virtual degrees of freedom ai, bi

are introduced to every site i. These are called ancillas and live in a generally simple Hilbert space H with dimension D, which is called the bond dimension. The next step is to build maximally entangled states between these virtual degrees of freedom of the neighbouring sites. One can choose for example the canonical maximally entangled state and denote this as

|Kiji = D X α=1 |αi_a i⊗ |αibj (1.6) where |αi_a

i with α = 1, · · · , D is the D-dimensional basis of the ancilla ai. This entangled

structure is projected upon the physical spin states through the linear map A(i) given by

A(i) = As_αβ···(i) |si_i⊗ hα|_a

i ⊗ hβ|bi ⊗ · · · (1.7)

The physical state is of course still a superposition of pure states and is constructed as a direct product of the above elements

|Ψi = N O i=1 A(i)O hi,ji |K_iji (1.8)

where the second tensor product runs over all the nearest neighbour sites hi, ji. One can thus imagine this construction as the physical spins having two virtual subsystems, the ancillas, which are entangled with the corresponding subsystems of its nearest neighbours. It is evi-dent that this construction can be extended to 2D systems, where they are called projected entangled-pair states or PEPS. Extensions to higher dimensional systems also exist but cal-culations in these systems proof to be very cumbersome.

(18)

For a one-dimensional periodic system Eq. (1.8) can be written as

|Ψi =X {si} X {αi} N Y i=1 Asi αi−1αi(i) |s1s2· · · sNi (1.9)

The periodic boundary conditions impose the condition α0 = αN. Also, the bond dimension

Di does not need to be the same for every entangled pair. The name matrix product states

then comes from the interpretation of the indices αi−1αi of Asαii−1αi(i) as the indices of a

Di−1× Di-dimensional matrix Asi. This allows to write the above equation in the simpler

form for a system with periodic boundary conditions

|Ψi =X {si} Tr " _N Y i=1 Asi_(i) # |s1s2· · · sNi (1.10)

which shows that the coefficient of every term |s1s2· · · sNi is given by a product of matrices.

For a system with open boundary conditions D0, DN = 1, so one can leave out the trace.

Ai−2 Ai−1 Ai Ai+1

si−2 si−1 si si+1

|Ψ{A}i =

Figure 1.4: General representation of an MPS in the bulk.

Interpreting the set of MPS with a set of bond dimensions {Di} on a lattice L as the manifold

M_{M P S{D}

i}

M_{M P S{D}_i_}= {|Ψ[A]i , ∀Asi_{(i) ∈ C}Di−1×Di_|∀s

i = 1, ..., qi, ∀i = 1, ..., N } (1.11)

it is clear that this is not a vector space because the sum |Ψ[A]i + |Ψ[B]i = |Ψ[C]i is not contained in MM P S{Di}. In the most general case |Ψ[C]i ∈ MM P S{2Di} because the matrices

Cs(i) are constructed as As(i) ⊗ Bs(i). Multiplication by a constant λ ∈ C however is a valid operation in M_{M P S{D}_i}, meaning that rays of states are included in the manifold. This is

simply obtained by multiplying a random matrix As(i) with λ. The fact that it doesn’t matter which site i is chosen already implicates some freedom in choice in the parameterization of the MPS. Defining now the matrix product parameter space AM P S as

AM P S = N

O

i=1

CDi−1×si×Di _(1.12)

In summary, the matrix product state thus represents the map

(19)

1.3.2 Entanglement entropy

In the introduction it was stated that tensor networks, and matrix product states in specific, obey an area law for the entanglement entropy between two arbitrarily cut subsystems. The proof thereof for MPS is fairly simple because MPS have two convenient properties. The first one is that dividing a 1-dimensional system in two parts simply constitutes of creating an edge between two neighbouring sites. The second is that the ancillas of an MPS are constructed as maximally entangled states, and because of the so-called monogamy of entanglement [17], it is therefore impossible that the constituents of such a maximally entangled state are entangled with another ancilla. The consequence is that for both subregions the upper bound of the entanglement entropy between itself and the other subsystem can be derived assuming that the correlations between the two are exclusively carried by the ancillas. Recalling Eq. (1.3), the Von Nuemann entanglement entropy becomes

S(ρM E) = −TrρM Eln(ρM E) = − D X i=1 1 Dln 1 D = ln(D) (1.14)

Where the well-known formula for the reduced density matrix for maximally entangled states was used: ρM E = D X i=1 1 D|ii hi| (1.15)

It is clear that the entanglement entropy for MPS is independent of the system size. This is however not the case for higher dimensional tensor networks like PEPS. For a 2-dimensional system the amount of maximally entangled pairs is given by the boundary with length L, so the entanglement entropy between a subregion and the environment is bounded by

(20)

Figure 1.5: The edge between a subregion and the environment for a 1-dimensional and a 2-dimensional system where L denotes the number of maximally entangled pairs on the boundary.

1.3.3 Transfer function

One usually starts with computing the overlap of 2 wave functions, which simplifies to the norm if the two MPS are the same. In doing so, an important tool in MPS-calculations, the transfer function, will appear.

hΨ0|Ψi = ht₁t2· · · tN| Tr " X {ti} N Y i=1 ¯ Bti_(i) # Tr " X {si} N Y i=1 Asi_(i) # |s₁s2· · · sNi = ht1t2· · · tN| Tr " X {ti} N Y i=1 ¯ Bti_(i) O X {si} N Y i=1 Asi_(i) # |s1s2· · · sNi = Tr " _N Y i=1 X si Asi_{(i) ⊗ ¯}_Bsi_(i) # = Tr " _N Y i=1 EM(i) # (1.17)

where the mixed transfer matrices EM(i) =

P

s As(i) ⊗ ¯Bs(i) are Dn−12 × Dn2-dimensional

matrices. In going from the first to the second line the property TrXTrY = TrX ⊗ Y was used. In order to obtain the third line, the property X1X2⊗ Y1Y2= X1⊗ Y1

X2⊗ Y2

(21)

Figure 1.6: Diagrammatic notation of the overlap of two MPS.

It is the matrix E(i), when A = B, which (especially in the thermodynamic limit) is called the transfer matrix3 For expectation values of operators the same strategy applies. Because any physical operator can be written as a sum over product operators ˆO1 ⊗ · · · ⊗ ˆON, one

just has to calculate the various single site operators ˆOi. In order to do this, just replace the

product overP

s At(i) ⊗ ¯As(i) for the norm with EO(i) =Ps,t ht| ˆOi|si At(i) ⊗ ¯As(i).

It was previously mentioned that in the case of open boundary conditions, one can even leave out the trace operation in Eq. (1.17). If one then wants to calculate the expectation value of an operator ˆOi which only acts on a single site i, one notices that the product in Eq. (1.17)

can be written as hΨ| ˆOi|Ψi = i−1 Y j=1 E(j) · EO(i)· N Y k=i+1 E(k) (1.19)

It is then natural to define the right product of transfer functions in the above equation as |ρi (i + 1) = E(i + 1) · · · E(N) where

ρR(i) =

X

s

As(i)ρR(i + 1)As(i)† (1.20)

with ρR(N ) defined as ρR(N ) = 1. The difference between |ρRi (i) and ρR(i) is that |ρRi (i)

is a D_i−12 × D2

N = D2i−1-dimensional vector and ρR(i) is a Di−1× Di−1-dimensional matrix.

Likewise, one can also define the left product as hρL| (i − 1) = E(1) · · · E(i − 1) where ρL is

similarly defined as

3_{This name is due to the analogy with the 1-dimensional Ising model. For every nearest neighbour} interac-tion, one has 4 possible spin configurations. Expressing this as a 2 × 2-matrix, i.e. the transfer matrix T , one can for example express the partition function as

(22)

ρL(i) =

X

s

As(i)†ρL(i − 1)As(i) (1.21)

with ρL(0) defined as ρL(0) = 1. These quantities are of course not limited to calculating the

specific expectation values given above. The norm is for example hΨ|Ψi = hρL(i − 1)|ρR(i)i =

TrρL(i − 1)ρR(i)] for every i = 0, ..., N . The merit of calculating the norm and expectations

values using these matrices is that the computational cost is significantly lower in comparison to using the transfer matrices E. Whereas in Eq. (1.17) the computational cost for computing the norm scales as (D2_bulk)3 = D6_bulk, it scales only as D3_bulk when using Eq. (1.20) and Eq. (1.21) because one only has to multiply matrices with the same dimensions as the original MPS matrices As. It is for this reason that open boundary conditions are preferred in MPS calculations. The notations ρL(i) and ρR(i) are of course chosen because they are no other

than the left and right reduced density matrices for a cut at site i.

1.3.4 Canonical forms

Upon inspecting the definition of MPS Eq. (1.9), one can conclude that the physical state is uniquely defined by the tensors Asi_{(i). The inverse, however, is not true. One can also choose}

a gauge transform with a new set of tensors ˜Asi_{(i) as}

˜ Asi n,m(i) = D X k,l=1 Xn,k(i)As_k,li (i)X_l,m−1(i + 1) (1.22)

such that the physical state |Ψ{ ˜A}i = |Ψ{A}i remains unchanged. Furthermore, it has been proven that this is the only freedom one has in reparametrizating the physical state[8][9]. A graphical representation is given in Fig. (1.7).

Figure 1.7: Gauge invariance in MPS

One can exploit this gauge freedom in constructing tensors Asi_{(i) in order to greatly}

sim-plify computations. A convenient choice for this gauge transform is for example when one constructs the tensor AL, in the above manner, as

˜ Asi Li,j(i) = D X n,m=1

Li,n(i)Asn,mi (i)L−1m,j(i + 1) (1.23)

(23)

Chapter 1. Introducing Tensor Networks 11 D X n=1 d X s=1 ˜ Asi Li,n(i) ˜A si Ln,j(i) = 1, ∀i (1.24)

This matrix L can be found by decomposing ρL(i) as ρL(i) = L†(i)L(i), or

ρLn,m =

D

X

k

Lk,n(i) ¯Lk,m(i) (1.25)

One can of course proof this by inserting Eq. (1.23) in Eq. (1.24), rearranging everything in order to obtain Eq. (1.25) and working everything out from there. Resorting to the diagram-matic representation however, prevents this overflow of brackets and indices. The so-called canonical left gauge thus boils down to Fig. (1.8)

Figure 1.8: The canonical left gauge for MPS, the index i is left out for convenience.

One can of course define a right gauge in a similar manner. The main advantage of choosing for example the left canonical gauge is that when calculating some expectation value at some site i one can contract over the whole left environment, resulting simply in a unity tensor because now ˜ρL(i) = 1. There exists also a mixed gauge where one choses a special site i and

uses Eq. (1.22) to construct a left gauge left of i and a right gauge to the right of i. Because this will not be used in this thesis, this will not be further elaborated on.

1.3.5 Thermodynamic limit

Until now, the MPS was described on a finite lattice. The leap to the thermodynamic limit will be made because in that limit the simulations are performed that are the core of this thesis. Because of the macroscopic size of interesting extended quantum systems, one is often interested in the bulk properties of these systems, far away from the boundaries. This also means that the effect of the boundaries on the bulk should vanish in the thermodynamic limit so that only these bulk properties remain.

If one also wants the MPS to be translation invariant, one has to choose the same tensor As _{for every site i.} _{This is computationally a very efficient way because the number of}

variational parameters remains constant so it is possible to go to the thermodynamic limit without increasing the computational cost. Although the continuous matrix product states discussed in the next chapter do not have to be translation invariant, it is chosen to use

(24)

translation invariance here mainly for convenience, as these results can be extended to the general non-translational invariant case with little effort.

Summarized, going to the thermodynamic limit with translation invariance boils down to the following: |Ψ[A]i = lim N →+∞ X {s} v†_L " _N Y i As # vR|{s}i (1.26)

where A is the same D × d × D-dimensional tensor for every site i with d the physical (spin) dimension of every site. v_L† and vR encode the boundary conditions. Because of the infinite

size of the system and the exponential decay of correlation functions, they will not have a physical meaning. Calculating the norm amounts to

hΨ|Ψi = lim N →+∞ vLv † L EN vRv†_R (1.27) This indicates that the transfer matrix E = P

sAs × ¯As should have a largest eigenvalue

equal to one if the MPS is properly normalized. For the overlap of two uniform MPS in the thermodynamic limit hΨ0|Ψi = lim_{N →+∞}TrENM indicates that the overlap is either 0 or 1,

depending on the largest eigenvalue of EM. This result is called the orthogonality catastrophe.

In the case of a largest eigenvalue µ0 6= 1, one can easily rescale the MPS |Ψ(A)i by replacing

A ⇒ A/√µ0. The fact that the largest eigenvalue of E is always a positive number follows

from the fact that E represents a completely positive map. It should be noted that E is not necessarily a hermitian matrix which means that the eigenvalues µi are in general4 complex

with |µi| < µ0, ∀i > 0. Another consequence is that the left and right eigenvectors of E are not

equal. Inspecting Eq. (1.20) and Eq. (1.21), one readily sees that in the bulk where every Asis equal, the two maps ρRand ρLeach convergence to a fixed point, which are the left and right

eigenvectors of E associated with eigenvalue µ0. For the case of open boundary conditions,

the computation of the norm is now simply hρL|ρRi. The calculation of the expectation value

of single site operators then becomes

hΨ| ˆOi|Ψi = hρL|EO|ρRi (1.28)

with EO=Ps,t ht| ˆOi|si At⊗ ¯As defined as before.

In order to be able to say that the bulk properties are really independent of the boundaries, one has to prove that all correlation functions disappear or become arbitrarily small over a substantial distance. It turns out that in uniform MPS they fall of exponentially with the

4_{The properties that the transfer matrix has an unique largest eigenvalue and that the corresponding left} and right eigenvectors are full-rank is called the injectivity of MPS. There are MPS were these conditions are not fulfilled but these do not correspond to physically relevant states, i.e. they correspond to so-called cat states.

(25)

distance. Consider a general correlation function hΨ| ˆC_n1Cˆ_m2|Ψi with n < m. Making use of the eigenvalue decomposition of the k-th power of the transfer matrix En,

Ek= |ρRi hρL| +

X

i

µk_i |µii hµi| (1.29)

one calculates this, by making use of the eigenvalue decomposition, Eq. (1.28) and the trans-lation invariance, as

hΨ| ˆC_n1Cˆ_m2|Ψi = hρ_L| ˆC_n1|ρ_Ri hρ_L| ˆC_m2|ρ_Ri +X

i

µm−n−1_i hρ_L| ˆC_n1|µ_ii hµ_i| ˆC_m2 |ρ_Ri (1.30)

The first part is just the product of the expectation values and is called the disconnected part of the correlation function. The rest consists solely of exponential decaying parts. This shows that connected correlations of every MPS decay exponentially. The consequence is that MPS are not suited in describing critical phenomena where correlation functions tend to decay according to power laws. This long range behaviour in MPS will be determined by the second largest eigenvalue µ1, assuming the eigenvalues are sorted according to decreasing

absolute value. The correlation length associated with this exponential decay is then given by ξ = −1/ln(|µ1|).

1.4 Tangent space

Having discussed some basic properties of matrix product states, the last thing left to intro-duce is the MPS tangent space TMM P S[A], which is defined for every point {A(n)|n ∈ Z}.

It was already shown that the MPS manifold for fixed {Di}, i ∈ Z is not a linear subspace.

It therefore makes sense to define a tangent space at every point |Ψ[A]i ∈ M_{M P S{D}_i}, which

is crucial in allowing one to optimize the MPS ansatz for a given Hamiltonian or to perform some time dependent variable principle method.

The definition of a general tangent vector (here in the thermodynamic limit) is |Ψ[B; A]i ,X

i∈Z

X

α

Bα(i) ∂

∂Aα_(i)|Ψ[A(i)]i (1.31)

where α is a collective index that represents both the physical and the virtual indices. It is clear that all the terms in the first sum are constructed by replacing Asi_{(i) with B}si_{(i) at site}

i in |Ψ{A(i)}i. A further discussion about this very useful construct will be postponed to the end of the next chapter were the continuous matrix product state variant will be discussed in somewhat more depth.

(26)

Chapter 2

cMPS

The continuous matrix product state |Ψ[Q, R1, ..., RNi with bond dimension D was originally

proposed by Verstraete and Cirac in 2010 as an extension of the MPS formalism[10]. Matrix product states were designed to deal with quantum theories on lattice systems, so the natural extension was to develop a similar method in order to describe quantum field theories. First, they were, as often happens, constructed to describe those theories in 1 spatial dimension. In 2019 the extension to d ≥ 2 spatial dimensionsal field theories has been made by Antoine Tilloy and J. Ignacio Cirac[11]. Meaningful calculations however are only possible for so-called Gaussian continuous tensor networks, where the multi-dimensional variant equivalent of Eq. (2.6) are Gaussian integrals. In general cases is even the definition of a continuous tensor network divergent. This is luckily not so much the case with cMPS. Correlation functions as well as various energies (kinetic, interaction, pairing,...) and densities are all calculable. Those calculations will be described in the next sections and will be applied to the infinite, periodical system of bosons which is the main subject of this thesis.

2.1 Fock Space

Consider a one-dimensional, continuous quantum system of length R = [−L/2, +L/2] on which q particle species labeled by α = 1, ..., q live, to which the creation operators ˆψα†

correspond. The particles can either be bosonic or fermionic so they have to satisfy the (anti)commutation relations ˆ ψα(x) ˆψβ(y) − ηα,βψˆβ(y) ˆψα(x) = 0 ˆ ψα(x) ˆψ_β†(y) − ηα,βψˆ_β†(y) ˆψα(x) = δα,βδ(x − y) (2.1)

where ηα,β = −1 if both particle species are fermionic and ηα,β = 1 if one of the particle

species is bosonic.

(27)

Chapter 2. cMPS 15

In order to construct a rigorous definition of the Hilbert space for the cMPS one needs to resort to the Fock construction. In the case of a single particle (of a random species α) the quantum mechanical wave function on R is a vector in the Hilbert space of square integrable functions on the region R, H(1)R = L2(R). The N -particle Hilbert space likewise

is the set of square integrable functions on the Cartesian product of N copies of R, H(N )R =

L2

R(N )_S,A, where S(A) denotes the symmetrical(anti-symmetrical) subspace of R(N ) in the case of bosons(fermions). The Fock space containing 0, ..., ∞ particles of a single species is then constructed as H(F )R = +∞ M N =0 H(N )R . (2.2)

The extension to a system with q particle species, denoted by α, is then readily constructed as H(F )R = +∞ M N0=0 ... +∞ M Nq=0 H{NRα}α=1,...,q (2.3) where H{NRα}α=1,...,q = L2 q Y α=1 RNα (S,A)α ! . (2.4)

is the Hilbert space of N particles of each species α. The unique vacuum state of this Fock space is given by |Ωi ∈ H{NRα=0}α=1,...,q.

2.2 cMPS manifold and definition

The Fock construction from the previous section allows a definition of the variational manifold of continuous matrix product states McM P S(D)∈ H

(F ) R as

M_{cM P S(D)}= {|Ψ[Q, R1, ..., RN]i |∀Q : R → CD×D; ∀Rα: R → CD×D, ∀α = 1, ..., N } (2.5)

where the cMPS |Ψ[Q, R1, ..., RN]i with bond dimension D, defined on the continuum R =

[−L/2, +L/2], is given by |Ψ[Q, R₁, ..., RN]i = Traux BPexp " Z +L/2 −L/2 dxQ(x) ⊗ ˆ1 + N X α=1 Rα(x) ⊗ ˆψα†(x) #! |Ωi (2.6)

With the trace running on an auxiliary space CD, i.e. the ancilla space. Q(x) and Rα(x),

which are D × D-dimensional matrices, are elements of the space of linear operators acting on the ancilla space L(CD) , CD×D. PexpR_−L/2+L/2dx(...) is the path ordered exponential, which

(28)

Chapter 2. cMPS 16

is completely analogous to the time ordered exponential (used for example in the calculation of the Dyson series), and orders the terms in the series expansion of the exponential from left to right for increasing values of x. |Ωi is the vacuum that is annihilated by the annihilation operators ˆψα(x). Finally, B is the D × D-dimensional matrix that encodes the boundary

conditions, which thus also lives in the space of linear operators acting on the ancilla space B ∈ L(CD). This parameter will not be of much use to us because the system that is used in this thesis is periodic and infinite so the boundary conditions are irrelevant. However, one can choose B = 1D when dealing with a finite system with periodic boundary conditions. In

case of open boundary conditions, one can choose B = vRv†_Lwith vR and vL D-dimensional

boundary vectors. The matrices Q(x) and Rα(x) are themselves also subject to boundary

conditions but because the system to which the cMPS ansatz will be applied in this thesis is infinite, they are not relevant here and will therefore not be discussed.

2.3 From MPS to cMPS

Before exploring continuous matrix product states some more, an explicit derivation of the connection between cMPS and MPS will be given. As mentioned before, cMPS is the continuum limit of a special subset of MPS. First, one can approximate the continuum R = [−L/2, +L/2] by a lattice L with spacing a, which thus contains N = L/a sites. The equivalent of the creation and annihilation fields ˆψ†α(x) and ˆψα(x) in quantum field theory

are the creation and annihilation operators ˆc†α(n) and ˆcα(n) that live on the lattice. The

re-lations that connect the former and the latter are ˆψ†α(x) = ˆc†α(n)/

√

a and ˆψα(x) = ˆcα(n)/

an MPS |Ψ{A}i with matrices As(n) with bond dimension D where s can take the values 0, α, β, (α, β), .... For an infinite-dimensional local basis, this MPS is only a formal construct.

2.3.1 Correspondence between As_{, Q and R}

Take now the continuum limit and send a → 0 so the number of lattice sites L/a in L becomes infinite. This is however not trivial. In general an MPS is not limited to a single Fock space in the thermodynamic limit and two different MPS are in general orthogonal due to the infrared orthogonality catastrophe. Because of this we need to restrict to a specific subset of MPS, namely those where the total average particle number h ˆN i is finite. Put differently, when taking the thermodynamic limit, we restrict our MPS to a single Fock space. The trick is then of course to find the correct description of the matrices As_{(n) in function of the cMPS}

matrices Q(x) and Rα(x). When looking at the definition of the cMPS it is clear that Q(x)

is connected to A0(n), with x = na, and that the latter has to be the dominant matrix in the MPS description. This is because of the finite number of particles, most sites will contain

(29)

Chapter 2. cMPS 17

zero particles. Likewise, R(x)α=1,··· ,q corresponds to As, ∀s 6= 0. The recipe Verstraete and

Cirac found was that the cMPS |Ψ[Q, R1, · · · , RN]i could be obtained from the continuum

limit a → 0 of the MPS |Ψ[A]i by identifying

A0(n) = 1D+ aQ(na)

Aα(n) =√aRα(na)

α 6= β ⇒ Aα,β(n) = a 2 h

Rα(na)Rβ(na) + ηα,βRβ(na)Rα(na)

i α = β ⇒ Aα,β(n) = a 2Rα(na) 2 ... (2.7)

where the vacuum |Ωi of the cMPS limit can be identified with the vacuum |0i = ⊗n∈L|0i_n,

∀n = −L/2a, −L/2a + 1, ..., +L/2a − 1 and ˆψ†α(x) = ˆc†α(n)/

√

a. This seems to look sensible but one wants to show the correspondence in an exact manner. The first logical step that immediately comes to mind to prove this connection is to expand the path ordered exponential

|Ψ[Q, R1, ..., RN]i = +∞ X N =0 Z −L/26x16···6xN6L/2 dx1· · · dxN Traux " B Q(x1) ⊗ 1 + q X α1=1 Rα1(x1) ⊗ ˆψ † α1(x1) ! × · · · × Q(xN) ⊗ 1 + q X αN=1 RαN(xN) ⊗ ˆψ † αN(xN) !# |Ωi (2.8)

As a side note, the rearrangement of the above equation by grouping the subsequent occur-rences of the Q term readily shows that cMPS can be viewed as a superposition over different particle number sectors.

|Ψ[Q, R1, ..., RN]i = +∞ X N =0 q X α1,...,αN=1 Z −L/26x16···6xN6L/2 dx1· · · dxN Traux " BMQ(−L/2, x1)Rα1(x1)MQ(x1, x2) · · · RαN(xN)MQ(xN, L/2) # ˆ ψ†_α₁(x1) ˆψ†α2(x2) · · · ˆψ † αN(xN) |Ωi (2.9) where MQ(x, y) = +∞ X k=0 Z x6z16···6zk6y dz1· · · dzkQ(z1) · · · Q(zk) = Pexp Z y x Q(z)dz ! (2.10)

(30)

Chapter 2. cMPS 18

A quick comparison with Eq. (2.3) reveals that this superposition is not equivalent to the different particle number sectors H{Nα}α=1,...,q

R in the Fock space H (F )

R because in the path

ordered integral only the total number of particles N is fixed whereas in the Fock space con-struction the number of particles Nα of every individual species is fixed.

In order to continue and to avoid too much overhead in the notation, the special case of a system with periodic boundary conditions and with a single particle species will be used, which boils down to leaving out the index α and the summations there-over. The bosonic version (but with open boundary conditions) of this system is also the system we will be using in our simulations. Eq. (2.8) then becomes

|Ψ[Q, R₁, ..., RN]i = +∞ X N =0 Z −L/26x16···6xN6L/2 dx1· · · dxN Traux " Q(x1) ⊗ 1 + R(x1) ⊗ ˆψ†(x1) ! × · · · × Q(xN) ⊗ 1 + R(xN) ⊗ ˆψ†(xN) !# |Ωi (2.11)

Now compare this with the MPS |Ψ{A}i where again n = −L/2a, −L/2a + 1, ..., +L/2a − 1, |mi_n= ˆc†(n)m/√m! and where Am(n) stands for the presence of m particles at site n

|Ψ[A]i =X {sn} Tr "+L/2a−1 Y n=−L/2a Asn_(n) # |· · · s_n−1snsn+1· · ·i =Traux "+L/2a−1 Y n=−L/2a A0(n) ⊗ |0i_n+ A1(n) ⊗ |1i_n+ · · · # =Traux "+L/2a−1 Y n=−L/2a 1 + aQ(na) ⊗ |0in+ √

aR(na) ⊗√a ˆψ†(na) |0i_n+ · · ·

#

(2.12)

where Eq. (2.7) is used in going from the second to the last line. It is now clear that the continuum limit a → 0 of this MPS is equal to the cMPS of Eq. (2.11) by again identifying |Ωi as |0i = ⊗n∈L|0in and by ordering the terms in the last line in increasing order in a,

which takes on the role of dx, to obtain the summation

+∞

P

n=0

.

This correspondence between MPS and cMPS is fruitful on the side of cMPS because one can derive for example the expectation values of various observables by taking the continuum limit of the MPS case. It is also useful for concluding that the entanglement entropy when working with open boundary conditions of one half of the chain with the other half is bounded in the upper limit by log(D).

(31)

Chapter 2. cMPS 19

2.3.2 Continuum limit of the transfer matrix

The various calculations of MPS derived in the previous chapter relied heavily on the transfer matrix E. If one wants to do calculations with cMPS one will encounter a similar object which will also be called a transfer matrix, but now denoted with T in order to avoid confusion with the MPS transfer matrix E. In Section 2.5 this transfer matrix T will appear naturally when calculating the norm of a cMPS, just as E was found when calculating the norm of a MPS. Getting ahead of things, it will first be derived that the transfer matrix T is the coefficient of the first order term of the expansion of the transfer matrix E in the lattice spacing a.

E(n) = X s As(n) ⊗ ¯As(n) = 1D+ aQ(na) ⊗ 1D+ a ¯Q(na) + √ aR(na) ⊗ √a ¯R(na) + · · · = 1D⊗D+ a Q(na) ⊗ 1D+ 1D⊗ ¯Q(na) + R(na) ⊗ ¯R(na) + O(a2)

(2.13)

The MPS transfer function can thus be written as E = 1D⊗D + aT + O(a2) where T thus

denotes the cMPS tranfer function. This expression also shows that where the largest eigen-value of E has a norm equal to one, the largest eigeneigen-value of T is zero. The consequence is that all the other eigenvalues of T are thus negative.

2.4 Regularization

Having defined the cMPS ansatz, one needs to check if the physical N -particle wave function is differentiable in each of its arguments, otherwise cMPS states will not have a finite kinetic energy ˆ T = Z +L/2 −L/2 ˆ t(x)dx (2.14)

where the kinetic energy density is given by

ˆ t(x) = N X α 1 2mα d ˆψα† dx (x) d ˆψα dx (x) (2.15)

Differently put, one requires that the norm of d ˆψα

dx (x) |Ψ[Q, R1, ..., RN]i is finite. In Section (2.7.1) the kinetic energy will be computed by making use of a generating functional but here the tactic will be to drag the field operator through the path-ordered exponential, so it can act directly on the vacuum. Using the results of Appendix C, one obtains

ˆ ψα(x) |Ψ[Q, R1, ..., RN]i = Tr h B ˆUα(−L/2, x)Rα(x) ˆU (x, +L/2) i |Ωi (2.16)

(32)

Chapter 2. cMPS 20 ˆ U (x, y) = P exp " Z y x dz Q(z) ⊗ ˆ1 + q X α=1 Rα(z) ⊗ ˆψα(x) !# ˆ Uα(x, y) = P exp " Z y x dz Q(z) ⊗ ˆ1 + q X β=1 ηα,βRβ(z) ⊗ ˆψβ(x) !# (2.17)

Where ˆU (x, y), ˆUα(x, y) ∈ L(H ⊗ CD) with CD the ancilla space. In differentiating Eq. (2.16)

with respect to x one encounters the derivatives

d Ûα dx (−L/2, x) = Ûα(−L/2, x) Q(x) ⊗ ˆ1 + q X β=1 ηα,βRβ(x) ⊗ ˆψ†_β(x) d Û dx(x, +L/2) = − Q(x) ⊗ ˆ1 + q X α=1 Rα(x) ⊗ ˆψα†(x) ˆ U (x, +L/2) (2.18)

Putting everything together results in

d ˆψα dx (x) |Ψ[Q, R1, ..., RN]i = Tr " B ˆUα(−L/2, x) [Q(x), Rα(x)] + dRα dx (x) ˆ U (x, +L/2) # |Ωi +Tr " B ˆUα(−L/2, x) q X β=1 [ηα,βRβRα− RαRβ] ⊗ ˆψ_β†(x) ˆ U (x, +L/2) # |Ωi (2.19) The norm of the term on the first line will be computed when calculating the expectation value of the kinetic energy in Section (2.7.1) and is thus finite, on the condition that the derivative of Rα(x) is well behaved. The term on the second line however creates particles

of any species at the fixed position x, indicating that it is not normalizable because of the divergent δ(0) contribution in the norm of Eq. (2.19). This forces one to impose the regularity condition

ηα,βRβRα− RαRβ = 0, ∀x ∈ R (2.20)

Differently put, the matrices Rα should satisfy the same statistics as the creation operators

they couple to. The system which forms the core of this thesis consists of a single species of bosons, so the condition Eq. (2.20) is automatically fulfilled. In systems with multiple species of bosons, the matrices Rα(x), Rβ(x) should thus commute at the same point x.

For a system consisting of a single fermionic species, the regularity condition imposes that R2(x) = 0, ∀x ∈ R. Anticipating Section (2.7.1), this factor will also be present in the N -particle wave function thus imposing the Pauli exclusion principle.

(33)

Chapter 2. cMPS 21

2.5 Overlap of two cMPS and the transfer matrix

Having introduced continuous matrix product states as the extension to the continuum of a certain class of matrix product states, calculating expectation values is the next step. In MPS one makes use of the transfer function E = P

sA ⊗ ¯A in calculating both the norm of the

wave function and the expectation values of operators. The same goes for cMPS. Because the transfer function E first appeared when calculating the overlap of two MPS, here the same strategy will be used. From now on only a single particle species will be used in calculations because this is also the system that is exploited in this thesis. The extension to multiple species α = 1, ..., q is straightforward.

In calculating the overlap between |Ψ[Q0, R0]i and |Ψ[Q, R]i the difficulty lies in the path ordered exponential. Therefore it is convenient to use Eq. (2.9) as a starting point and try to write the overlap back again as a path ordered exponential of some quantity.

The overlap hΨ[Q0, R0]|Ψ[Q, R]i then becomes

hΨ[Q0,R0]|Ψ[Q, R]i = +∞ X N =0 Z −L/26x16···6xN6L/2 dx1· · · dxN Traux " BPexp Z x1 −L/2 Q(z)dz R(x1)Pexp Z x2 x1 Q(z)dz · · · R(x_N)Pexp Z L/2 xN Q(z)dz # × Tr_aux " BPexp Z x1 −L/2 Q0_(z)dz R(x1)Pexp Z x2 x1 Q0_(z)dz · · · R(x_N)Pexp Z L/2 xN Q0_(z)dz # (2.21) Where the following orthogonality relation was used

hΩ| ˆψ(x0_N) · · · ˆψ(x0₂) ˆψ(x0₁) ˆψ†(x1) ˆψ†(x2) · · · ˆψ†(xN) |Ωi = δ(x01− x1)δ(x02− x2) · · · δ(x0N − xN)

(2.22) In order to write this as a path-ordered exponential expansion, one needs everything in one Tr[· · · ], therefore use the identities for D ×D-dimensional matrices already used in calculating the overlap of a MPS in addition to eX ⊗ eY _{= e}(X⊗1+1⊗Y ) _{to obtain}

(34)

Chapter 2. cMPS 22 hΨ[Q0, R0]|Ψ[Q, R]i = +∞ X N =0 Z −L/26x16···6xN6L/2 dx1· · · dxN Traux " B ⊗ BPexp Z x1 −L/2Q(z) ⊗ 1D+ 1D ⊗ Q0_(z)dz R(x1) ⊗ R(x1) · · · R(xN) ⊗ R(xN)Pexp Z L/2 xN Q(z) ⊗ 1D+ 1D⊗ Q0(z)dz # (2.23)

Comparing this to Eq. (2.8) one readily notices that this is indeed the path-ordered exponen-tial of the quantity T, the transfer matrix for cMPS derived in Eq. (2.13) when Q = Q0 and R = R0 hΨ[Q0, R0]|Ψ[Q, R]i = Traux " B⊗BPexp Z L/2 −L/2Q(x)⊗1 D+1⊗Q0(x)+R(x)⊗R0(x)dx # (2.24) For the norm, this equation can thus be written in a formula similar to the formula of the norm of a general MPS hΨ[A]|Ψ[A]i = Tr

" N Q n=1 E(n) # : hΨ[Q, R]|Ψ[Q, R]i = Tr_aux " B ⊗ BPexp Z L/2 −L/2T(x)dx # (2.25)

2.6 Reduced density matrices

It is often useful when thinking about cMPS to make the connection with MPS and reason first for the latter to understand how to do something for the former. In Eq. (2.24) a formula for the overlap of two cMPS was formulated. This is however of little practical use because one still has to calculate a path-ordered exponential integral. With MPS, in particular in the thermodynamic limit, doing the various calculations was greatly simplified by making use of the left and right reduced density matrices. A very similar simplification exists for cMPS. One can also define a left reduced density matrix ρL and a right reduced density matrix ρR

for a cMPS. Before giving the correct definitions, it proves useful to visualize cMPS some more.

If one wants to calculate the expectation value of ˆO =P

i

ˆ

Oi for a MPS one can visualize the

(35)

Chapter 2. cMPS 23

Figure 2.1: The contraction over a single site operator ˆOi[13]

Having this in mind, one can then define two linear maps T , ˜T : L(CD_{) 7→ L(C}D_{) that map}

virtual operators f , which are D × D matrices. T and ˜T are the mathematical formulation of Fig. (2.1) for cMPS, depending if the open legs of the MPS contraction are either on the right or the left side

T(x)_{(f ) = Q(x)f + f Q(x)}†_{+ R(x)f R(x)}†_, _(2.26)

˜

T(x)(f ) = f Q(x) + Q(x)†f + R(x)†f R(x) (2.27) As the system considered is using open boundary conditions, one can thus define a left and a right reduced density matrix through the initial conditions ρL(−L/2) = vLv†_Land ρR(L/2) =

vRv†_R and the first order differential equations

d dxρL(x) = ˜T (x) _ρ L(x) (2.28) d dxρR(x) = −T (x) _ρ R(x) (2.29) where ρL(x), ρR(x) ∈ L(CD). Just as in the case with the reduced density matrices of the

MPS, one can associate vectors |ρLi , |ρRi ∈ CD⊗ CD with these D × D matrices. The above

equations are formally solved as

(ρL(x)| = (ρL(−L/2)|Pexp Z x −L/2 T(z)dz |ρ_R(x)) = Pexp Z L/2 x T(z)dz |ρ_R(L/2)) (2.30)

Comparing these equations with Eq. (2.25) it is easy to see that the norm can thus be written as

(36)

Chapter 2. cMPS 24

2.7 Practical calculations

Before the calculation of expectation values is described, one may argue that Eq. (2.31) is not of any more practical use than Eq. (2.25). The former does not contain a path-ordered exponential but the definitions of the reduced density matrices however, do. In order to do practical calculations one then must find a way to bypass those path-ordered exponentials. One can take the differential equations Eq. (2.28) as a starting point and make numerical approximations of ρL(x) and ρR(x) given the matrices Q(x) and R(x). The approach in this

thesis is simpler and exploits the structure of Eq. (2.28).

Because of the periodic potential that will be applied on our infinite system of bosons, a decomposition of the relevant quantities in Fourier modes is chosen. The sinuses and the cosines then fix the x-dependencies of the system and the only task left to do is to find the correct matrix coefficients for the series expansions.

A similar problem was encountered by the Swiss-American Physicist Felix Bloch when he wanted to describe the conduction of electrons in crystalline[15]. He found that when one wants to describe the one-dimensional wave function of an electron in a periodic potential one could describe it as

ψ = exp(ikx)u(x) (2.32)

Bloch was however not the first to solve this problem. The French mathematician Gaston Floquet had almost fifty years earlier proven a similar theorem that was later named after him. It stated that solutions of differential equations of the form, where C(x) is a periodic function with period P ,

d

dxF (x) = C(x)F (x) (2.33)

always had a solution where F could be written as the product of a complex exponential and a function ˜F (x) with periodicity P [16]. This was the same result as Bloch would derive some time later1.

Applying Floquet’s theorem on Eq. (2.28) (ρRis completely analogous), where thus Q(x) and

R(x) are expressed as finite Fourier series and after eliminating the common factor exp(λx), one obtains

λ ˜ρL(x) = ˜ρL(x)Q(x) + Q(x)†ρ˜L(x) + R(x)†ρ˜L(x)R(x) − ˜ρL(x)0 (2.34)

with λ a complex number. After expanding ˜ρL(x) in a finite Fourier series as well, one obtains

a system of coupled differential equations for the different Fourier modes of ˜ρL(x) for given

1_{Floquet was however also not the first to prove this theorem, despite it being named after him. The first} one to prove this was George William Hill in 1877, followed by Floquet in 1883 and Alexander Lyapunov in 1892. Bloch obtained his result only in 1928.

(37)

Chapter 2. cMPS 25

Q(x), R(x). By comparing the first two terms with the third term at the right-hand side it is clear that one has to choose twice the number of modes for Q(x) one uses for R(x). For the sake of completeness, the formula for the right reduced density matrix is given by

− λ˜ρR(x) = Q(x) ˜ρR(x) + ˜ρL(x)Q(x)†+ R(x) ˜ρL(x)R(x)†+ ˜ρL(x)0 (2.35)

2.7.1 Expectation values

In order to calculate the expectation values of normally ordered operators ˆO =: O[{ ˆψ†}{ ˆψ}] : one can introduce a generating functional Z[ ¯J , J ]

Z[ ¯J , J ] = Traux B ⊗ BPexp Z L/2 −L/2 T(x) + J (x)[R(x) ⊗ 1] + J (x)[1 ⊗ R(x)]dx (2.36)

Working in a system with open boundaries allows to write the Traux

B ⊗ B · · · from the above expression as (l(−L/2)| · · · |r(L/2)). The advantage of this approach is that it can be used in the same manner as in quantum field theory to express the expectation value of every possible measurable quantity of the system. The recipe goes as follows. First, make sure that all the operators are normally ordered; next, take for every ˆψ†in the operator ˆO the functional derivative ∂

∂ ¯J of Z[ ¯J , J ]. Likewise, taking the functional derivative ∂

∂J of Z[ ¯J , J ] corresponds with a ˆψ in the operator O. These functional derivatives are the natural extension of ordinary derivatives of functions to these so-called functional derivatives of functionals i.e. functions of functions. They obey the axiom

∂

∂J (x)J (y) = δ(x − y) (2.37)

and follow the same rules for composite functions as ordinary derivatives.

An example will provide some more clarity. Take for example the expectation value

hΨ[Q, R]| ˆψ†(x) ˆψ†(y) ˆψ(y) ˆψ(x)|Ψ[Q, R]i (2.38) which is the integrandum of the expectation value of the generic two-particle interaction

ˆ W = 1 2 RL/2 −L/2dx RL/2 −L/2dy w(x, y) ˆψ†(x) ˆψ†(y) ˆψ(y) ˆψ(x).

Applying the recipe just described means that one can express this as

hΨ[Q, R]| ˆψ†(x) ˆψ†(y) ˆψ(y) ˆψ(x)|Ψ[Q, R]i = ∂

4_{Z[ ¯}_{J , J ]} ∂ ¯J (x)∂ ¯J (y)∂J (y)∂J (x) _¯ J ,J =0 (2.39) Taking these functional derivatives then amounts to

(38)

Chapter 2. cMPS 26

hΨ[Q, R]| ˆψ†(x) ˆψ†(y) ˆψ(y) ˆψ(x)|Ψ[Q, R]i = θ(y − x)(l(x)|R(x) ⊗ R(x)Pexp Z y x T(z)dz R(y) ⊗ R(y)|r(y)) + [x ↔ y] (2.40)

If one wants to find a closed expression for ˆW , one can express (l(x)|R(x)⊗R(x)Pexp Ry

x T(z)dz

as a first order differential equation similar to Eq. (2.34) and solve R(l)x (x) = (l(x)|R(x)⊗R(x)

for every x ∈ L and y > x as Rx(l)(y).

The expectation value of the kinetic energy hΨ| ˆT |Ψi, which will prove to be a bit more cumbersome due to the limits associated to the derivatives, is then

hΨ[Q, R]| d dxψˆ †_(x) d dxψ(x)ˆ |Ψ[Q, R]i = lim x→y d2 dxdy ∂2Z[ ¯J , J ] ∂ ¯J (x)∂J (x) _¯ J ,J =0 = lim x→y d2 dxdy θ(y − x)(l(x)|(1D⊗ R(x))Pexp Z y x T(z)dz (R(y) ⊗ 1D)|r(y)) + θ(x − y)(l(y)|(R(y) ⊗ 1D)Pexp Z x y T(z)dz (1D⊗ R(x))|r(x)) (2.41) In order to evaluate this expression, first act with d

dy on the terms in the big brackets. The results of the action of d

dy upon the different y-dependent factors in the first term are listed below. The second term is analogous, just as the derivatives with respect to x.

d dyθ(y − x) = δ(y − x) d dyPexp Z y x T(z)dz = T(y)Pexp Z y x T(z)dz d dy(R(y) ⊗ 1D) = dR(y) dy ⊗ 1D d dy|r(y)) = −T|r(y)) (2.42)

The derivatives of the Heaviside functions nicely cancel for both derivatives with respect to y and x because the matrices R(x) obey the same statistics as the creation operators. One can obtain the derivative of the path-ordered exponential on the second line by writing it as an infinite sum Pexp Ry

x T(z)dz = 1 + R y x T(z1)dz1 + Ry x Rz1 x T(z1)T(z2)dz2dz1+ · · · . The

(39)

Chapter 2. cMPS 27

reduced density matrix. Filling the above expressions for the various derivatives listed in Eq. (2.41) results in hΨ[Q, R]| d dxψˆ † (x) d dxψ(x)ˆ |Ψ[Q, R]i = lim x→y " θ(y − x)(l(x)| h T(x), 1D⊗ R(x) i + 1D ⊗ dR(x) dx Pexp Z y x T(z)dz × h T(y), R(y) ⊗ 1D i + dR(y) dy ⊗ 1D |r(y)) +θ(x − y)(l(y)| h T(y), R(y) ⊗ 1D i +dR(y) dy ⊗ 1D Pexp Z x y T(z)dz × _h T(x), 1D⊗ R(x) i + 1D⊗ dR(x) dx |r(x)) # (2.43) The four terms with the commutators, e.g.

h

T(x), 1D⊗ R(x)

i

, can be rewritten using T(x) = Q(x) ⊗ 1 + 1 ⊗ Q(x) + R(x) ⊗ R(x) to obtain for the given example

h T(x), 1D ⊗ R(x) i = 1 ⊗ h Q(x), R(x) i

. After applying this to all four terms and taking lim

x→y, which causes the

Heaviside functions and the path-ordered exponentials to disappear so the tensor products can be worked out, one finally obtains

hΨ[Q, R]| d dx ˆ ψ†(x) d dx ˆ ψ(x) |Ψ[Q, R]i = (l(x)| [Q(x), R(x)] + dR(x) dx ⊗ [Q(x), R(x)] + dR(x) dx |r(x)) (2.44)

2.7.2 Gauge transformations and canonical forms

The gauge freedom of a cMPS is easiest derived from the correspondence between the matrices As and Q and R, according to Eq. (2.7). In the previous chapter it was shown that the state |Ψ{A}i remains unchanged when G(n − 1)G−1(n) is placed between As(n − 1) and As(n), with G(n) a position dependent, invertible matrix. Using the correspondence Eq. (2.7), and defining g(na) , G(n) one can easily compute, using some Taylor expansions

(40)

Chapter 2. cMPS 28

A0(n) = 1D+ aQ(na) ⇒ ˜A0(na) = g((n − 1)a)−1A0(n)g(na)

= g((n − 1)a)−1g(na) + ag((n − 1)a)−1Q(na)g(na) = 1D + a −dg −1 dx (na)g(na) + g(na) −1 Q(na)g(na) + O(a2) = 1D + a ˜Q(na) + O(a2)

A1(n) =√aR(na) ⇒ ˜A1(na) = g((n − 1)a)−1A1(n)g(na) =√ag((n − 1)a)−1R(na)g(na) =√ag(na)−1R(na)g(na) + O(a32)

=√a ˜R(na) + O(a32) · · · (2.45) And using dg −1 dx (na)g(na) = −g(na) −1dg

dx(na), one obtains the gauge transformation rules for cMPS:

˜

Q(x) = g(x)−1Q(x)g(x) + g(x)−1dg

dx, R(x) = g(x)˜

−1_R(x)g(x) _(2.46)

where one requires the matrix-valued function g(x) to be twice differentiable in order to have a defined first-order derivative for ˜Q(x) and ˜R(x). The extra term g(x)−1dg

dx in the expression of the gauge transformation of Q(x) can be recognized as the infinitesimal parallel transport encountered in Yang-Mills gauge theories where it appears for example in the gauge transfor-mations of the W bosons.

The caveat with these transformations is that when one applies a gauge transform on the state |Ψ[Q, R]i, one also has to apply these transformations on the boundaries; i.e. B =˜ g−1(L/2)Bg(−L/2). This poses however an extra condition on the possible gauge transfor-mations, which is dependent on the chosen boundary conditions. In the case of PBC, i.e. B = 1D, this results in g(L/2) = g(−L/2). OBC are a bit more cumbersome because they

result in the coupled equations

v_L†g(−L/2) = v†_L g(L/2)−1vR= vR (2.47)

One then has to choose to either solve these equations for a function g(x) which satisfies both2 or only solve one, e.g. for g(−L/2), and include the other boundary vector, e.g. vR,

in the variational parameters. Although one has to bear in mind these restrictions when constructing for example the canonical forms discussed in the following paragraph, they are not of much importance to this thesis because of the infinite size of the system that will be

2

(41)

Chapter 2. cMPS 29

considered.

As is the case with MPS, one can use this gauge freedom to define a left and a right canonical form in order to simplify calculations. Applying the gauge transformations Eq. (2.46) on the definition of the reduced density matrices Eq. (2.28), one easily derives the following transformation rules

˜

ρL(x) = (g−1(x))†ρL(x)g−1(x) ρ˜R(x) = g(x)ρR(x)g(x)† (2.48)

Using this, one can define ρL(x) := 1D. Eq. (2.28) then becomes ˜Q(x)†+ ˜Q(x) + ˜R(x)†R(x) =˜

0. If one would then write this again in function of Q(x), R(x) and g(x) using Eq. (2.46), one would find that g(x) is fixed up to a unitary prefactor, an element of U(D), which can then be used to diagonalize r(x). Hence obtaining the left-canonical form. The right-canonical form is analogously defined by choosing ρR(x) := 1D and using the remaining U(D) freedom

of g(x) to diagonalize ρL(x) for every x.

2.8 Tangent space

As for MPS, a short description will be given about the tangent space of a cMPS, but this will not be elaborated too much upon as only a basic notion of the tangent space will be needed. For MPS it has been proven that it can be described as having the structure of a variational manifold with a well-defined tangent space[14]. The infinite dimensionality of the variational parameter space and Hilbert space for cMPS however, complicates such a rigorous description.

In the following, this remark is ignored because of the usefulness of tangent space vectors in e.g. finding the ground state and in their use in the time-dependent-variational-principle (TDVP) scheme. The construction of the tangent space is analogous to that for MPS as described in Section (1.3.1). Every open subset of cMPS with fixed bond dimension D has the structure of a manifold, denoted by McM P S. For every point |Ψ[Q, R]i ∈ McM P S one

can construct thus a tangent space T|Ψ[Q,R]iMcM P S ⊂ H. The definition of a tangent vector

|Φ[V, W ; Q, R]i in T|Ψ[Q,R]iMcM P S is given, in analogy with Eq. (1.31), by

|Φ[V, W ; Q, R]i = Z L/2 −L/2 dx D X α,β Vαβ(x) ∂ ∂Qαβ(x) + Wαβ(x) ∂ ∂Rαβ(x) |Ψ[Q, R]i (2.49)

(42)

Chapter 2. cMPS 30

Without going into detail, the description of |Φ[V, W ; Q, R]i in terms of V (x) and W (x) is redundant because of the gauge freedom described in the previous section. This implies that the representation of cMPS tangent vectors has a gauge invariance and it is given by

V (x) ← V (x) + [Q(x), h(x)] +dh dx(x) W (x) ← W (x) + [R(x), h(x)]

(2.50)

with h(x) ∈ gl(C, D) ≡ CD×D. In the above formulas, the matrix B which encodes the boundary conditions, has been considered unchanged. This forces conditions on the gauge transformation g(x) introduced in the previous section and on the set of allowed functions h(x):

g(+L/2)Bg(−L/2)−1= B

h(+L/2)B − Bh(−L/2) = 0 (2.51)

Further ignoring the implications on h(x) and B, in constructing the tangent space for some base point |Ψ[Q, R]i one would also like to impose the orthogonality of the base point to its own tangent space

hΨ[Q, R]|Φ[V, W ; Q, R]i = Z

P

dx(l(x)|V (x) ⊗ 1 + W (x) ⊗ R(x)|r(x)) + B.T. (2.52) In order to do this one can choose either a left gauge or a right gauge (not to be confused with the canonical left and right gauges ρL= 1 and ρR= 1) where e.g. the left gauge amounts to

setting

(l(x)|hV (x) ⊗ 1D+ W (x) ⊗ R(x)

i

= ρL(x)V (x) + R(x)†ρL(x)W (x) = 0 (2.53)

Which amounts to the following identity in the canonical left gauge

(43)

Chapter 3

Methods

The main goal of this thesis is to describe a one-dimensional system of spinless, non-relativistic bosons with a contact interaction in a periodic potential on the continuum. The starting point isthe Hamiltonian: ˆ H = Z +∞ −∞ dx d ˆψ † dx (x) d ˆψ dx(x) + c ˆψ†(x)2 ψ(x)ˆ 2+ ˆV (x) ˆψ†(x) ˆψ(x) (3.1) where the bosons interact with each other through a delta-interaction with strength denoted by c. When considering the above Hamiltonian without a potential term, it is equal to the Lieb-Liniger model, which is an example of an exactly solvable quantum system. Lieb and Liniger determined the ground state[18], and Lieb alone found the whole excitation spectrum[19]. The Lieb-Liniger model contains a continuous U(1) symmetry corresponding to particle number conservation, so one would expect that the ground state would be de-generate because of spontaneous symmetry breaking. This is however not the case because of the Mermin-Wagner-Hohenberg-Coleman theorem[20][21][22]. This theorem states that for spatial dimensions d ≤ 2 no spontaneous symmetry breaking can occur at finite temperatures. It has thus an unique groundstate.

3.1 Lieb-Liniger model: benchmark 1

In order to study the Lieb-Liniger model, a chemical potential will be added to the Hamilto-nian of the system in order to produce a non-trivial ground state

ˆ ˜ H = Z +∞ −∞ dx d ˆψ † dx (x) d ˆψ dx(x) + c ˆψ†(x)2 ψ(x)ˆ 2− µ ˆψ†(x) ˆψ(x) (3.2) which determines the expectation value of the particle density. Here, the units m = 1/2, ~ = 1 have been used. Although this model seems to have two parameters, the particle density ρ and the strength of the delta interaction c, rescaling the length x → ρx shows that it is

(44)

Chapter 3. Methods 32

essentially a one-parameter model, here denoted by γ = c/ρ. This implies that dimensionless quantities like e(c, ρ)/ρ3, with e(c, ρ) the energy density, can be expressed as functions of only γ by rescaling everything with the density

e(c, ρ) = ρ3e(γ, 1) , ρ3e(γ) (3.3)

This is translated into the matrices Q(x) and R(x) by the mapping

Q(x) 7→ ρQ(x) , R(x) 7→√ρR(x) (3.4)

It is thus natural to choose ρ = 1 when performing simulations. There are two ways in varying γ. The first is keeping the interaction strength c constant and varying the density, and the second is keeping the density constant and varying the interaction strength. In combination with the fact that it is natural to work with ρ = 1 in the cMPS ansatz, the natural choice is to let the interaction strength vary. It is interesting to note that as a consequence the Tonks-Girardeau limit can not only be reached by sending c → +∞, but that it can also be simulated by a very low density bosonic gas. However, the focus will be on the first inter-pretation. In the limit of infinite interaction the particles are not allowed to pass through each other1 anymore. Hence imitating the Pauli exclusion principle of fermions, i.e. fermion-ization. It was proven that such a one-dimensional hardcore boson gas indeed shares some properties with a free fermion gas[26], discussed in the next section. The main advantage of the c → +∞ limit is that the exact value of some quantities like the energy density are exactly known, e.g. e(γ)/ρ3 = π2/3. In the c → 0 limit the boson gas is clearly free.

Computing the energy density with the cMPS ansatz amounts to

e(γ) − µρ = hΨ[Q, R]| d ˆψ † dx (x) d ˆψ dx(x) + c ˆψ†(x)2 ψ(x)ˆ 2− µ ˆψ†(x) ˆψ(x)|Ψ[Q, R]i = hΨ[Q, R]|DR ⊗ DR + cR2⊗ R2− µR ⊗ R|Ψ[Q, R]i (3.5) where DR = [Q, R] + dR/dx. The energy associated with the chemical potential is separated such that e(γ) has the same meaning here as in the Bethe ansatz, which will be discussed in Section (3.4).

1