MSc Physics
Track: Theoretical Physics
M
ASTERT
HESISCombining Variational Optimization with
Entanglement Renormalization in a Tensor
Network framework
Boris Ponsioen
10279806August 2017
60 ECTS Supervisor: Dr. Philippe Corboz Examiner: Dr. Edan LernerAbstract
In tensor network algorithms for obtaining ground states of one dimensional infinite critical quantum spin systems, two distinctive approaches can be identified. Methods using variational optimization of Matrix Product States (MPS), such as infinite Density Matrix Renormalization Group (iDMRG), tend to produce accurate results in terms of short range observables, for example energy, but are never able to reproduce polynomially decaying correlation functions at large distances, which are typical for ground states of critical systems. In gapped systems, however, the area law for the entanglement implies that an MPS is especially suitable for representing ground states. On the other hand, the Multiscale Entanglement Renormalization Ansatz (MERA), is able to show accurate long range correlations, though it comes at the cost of slow convergence. Our scheme combines elements from both approaches, using the ground state obtained through variational optimization to generate an approximate coarse graining transformation, and aims to provide a controlled tradeoff between long range accuracy and performance, in (near) critical systems.
Samenvatting
In het veld van gecondenseerde materie, een subveld van theoretische natuurkunde, bestudeert men vaak microscopische modellen van interacterende deeltjes die fenomenen kunnen verklaring die wij op macroscopische schaal waarnemen, terwijl zij zo eenvoudig mogelijk worden gehouden. Aangezien het voor veel van deze modellen niet mogelijk is om een exacte beschrijving te vin-den voor hun gedrag, bievin-den numerieke simulaties met behulp van computers uitkomst. Echter bestaan systemen op macroscopische schaal uit enorme hoeveelheden interacterende deeltjes, waardoor het onmogelijk is om de systemen in hun volledigheid te simuleren.
Een bekende techniek uit de theoretische fysica die het mogelijk maakt om uitspraken te doen over het gedrag van modellen in de thermodynamische limiet, een andere veel gebruikte term voor macroscopische schaal, is de Renormalisatie Groep. Door middel van ’uitzoomen’, waarbij informatie over clusters van deeltjes wordt samengevat en de fijne details — welke geen rol spelen op grote schaal — worden weggegooid, is het mogelijk om een effectieve versie van het model te creëren op grotere schaal.
In deze scriptie wordt uitvoerig gebruik gemaakt van zogenaamde Tensor Networks, objecten die op efficiënte wijze toestanden van kwantumsystemen kunnen beschrijven. Deze netwerken vormen de basis voor een scala aan algoritmen voor numerieke simulaties van deze systemen. Wij hebben ons in dit onderzoek gericht op de ontwikkeling van een methode die zulke algorit-men, die vallen onder de noemer van ’variational optimization’, combineert met renormalisatie. Deze methode heeft als doel om op een efficiënte manier simulaties van kwantumsystemen op grote schaal mogelijk te maken, waarbij ook het gedrag van modellen over grote afstanden cor-rect kan worden gereproduceerd. Het effect van de methode is onderzocht en vergeleken met conventionele algoritmen, door middel van het testen op een eenvoudig eendimensionaal model waarvan een exacte beschrijving bestaat.
De resultaten die zijn verkregen banen de weg voor simulaties van ingewikkeldere modellen, waarvan er geen exacte beschrijvingen bestaan en het gedrag op grote schaal nog onbekend is. Bovendien is het mogelijk om de methode uit te breiden naar een tweedimensionale vari-ant. Simulaties van tweedimensionale modellen op basis van Tensor Networks zijn zwaarder dan in één dimensie en onze methode is een potentieel krachtige aanvulling op bestaande algoritmen.
Contents
1 Introduction 1
2 Introduction Tensor Networks 3
2.1 Diagrammatic notation . . . 4
2.2 Matrix Product States . . . 6
2.3 Operators . . . 9
2.4 Entanglement . . . 10
2.5 Multiscale Entanglement Renormalization Ansatz (MERA) . . . 11
2.6 Application - Transverse Ising model . . . 14
3 Methods 16 3.1 Introduction . . . 16
3.2 Energy minimization (infinite DMRG) . . . 17
3.3 Alternative ground state algorithms . . . 23
3.4 Coarse graining and disentangling . . . 28
3.5 Operators . . . 33
4 Results 36 4.1 Ground state simulations . . . 36
4.2 Entanglement renormalization operators . . . 42
4.3 Combined scheme . . . 43
5 Discussion and outlook 55 5.1 Outlook . . . 57
Bibliography 61
1
Introduction
Many problems in theoretical physics cannot be solved analyti-cally, and therefore require numerical methods in order to find approximations. Such methods have become an important part of modern physics and the field of computational physics is con-sidered its own subfield, positioned close to theoretical physics. One of the most commonly studied but also most complex prob-lems is finding the ground state for a given quantum mechanical system. In the following chapters, we take the Ising model as an example for our methods, which consists of an evenly-spaced lattice where each sites houses a spin particle and these spins are allowed to interact with their nearest neighbours [1]. In this thesis, a new technique for finding ground state solu-tions for quantum physical systems is investigated. The scheme consists of two main parts: a conventional ground state algo-rithm, any one of the widely used tensor network methods avail-able, combined with a coarse graining procedure, in close resem-blance to the well-known Renormalization Group [2]. Usually, the ground state algorithms — without any renormalization — can be very effective at describing the behaviour of the ground state for small system sizes. In some cases, also large and even infinite systems can be accurately simulated, but this is highly dependent on the model under consideration.
Tensor network methods, which we will describe in detail in the following chapters, have become popular in recent years due to their efficient representation of quantum states, their concep-tual simplicity and the natural basis they provide for powerful algorithms. The various ground state algorithms [3–5], which we will use as one part of the scheme in this thesis, are able to accurately reproduce the ground state physics in the thermody-namic limit in non-critical systems. If the system is critical, or close to criticality, the algorithms tend to be accurate at short range, but the long-range correlations are no longer accurately reproduced [6, 7].
introduction 2
The second part of the combined scheme is a coarse graining procedure, in which we map the result of the ground state algo-rithm onto a new lattice that is generated by ‘zooming’ out and combining sites into blocks that form the sites on the new lat-tice. Importantly, this transformation is not exact, as it would result in an exponential scaling in the dimensionality of the ten-sors that represent the blocks, but is instead designed to capture the most relevant information and dispose of the rest. By such a transformation, we aim to derive an effective version of the system at larger scale [8].
Our combined scheme then consists of alternating between both parts: first a ground state calculation that identifies the most relevant information on the current scale of the lattice, on which then a coarse graining step is based, forming a new lattice at larger scale. On the new lattice, the ground state calculation can be performed again, continuing the cycle. After a number of these steps the coarse grained lattice effectively represents a large system due to the exponential scaling of the blocks. The goal is then for this representation to be efficient, by keeping only little information in each coarse graining step, while still being fairly accurate. Especially the simulation of long-range physics should then be achievable, up to higher accuracy than the conventional ground state algorithms are able to provide.
2
Introduction Tensor
Networks
The field of tensor networks has been of great interest, especially in recent years. Tensor networks can come in many different forms, but are characterized by a connected network of
multi-dimensional1 arrays. The specific configuration of these tensors 1. The dimensionality of the individ-ual tensors can vary, depending on the application.
and their individual properties can be adapted to the applica-tion, allowing for a high amount of flexibility.
Its success is in part due to the development of the Density Ma-trix Renormalization Group (DMRG) algorithm, by S. White in 1992 [3], which still stands as one of the most powerful numerical methods in one dimensional quantum lattices today. However, the original formulation of the algorithm was quite different from the one used in current literature, including this thesis, as it did not involve tensor networks directly.
One of the main motivations for this new type of numerical meth-ods is based on the fact that the Hilbert space containing the states of quantum systems grows exponentially with the sys-tem size. Therefore, when a given quantum state is evolved inO (poly (N)) time, it can be shown [9] that the subregion of the total Hilbert space that can be reached is also exponentially small. It would take an O(exp(N)) amount of time to reach most of the Hilbert space. In a realistic scenario, where a typi-cal system is ofO(1023) sites, it is clear that a complete (exact)
simulation is infeasible.
A solution came from the realization that amount of entan-glement in the low-energy states of systems with local gapped Hamiltonians is bounded from above. This is formulated in the so-called area law of entanglement, which states that the entan-glement entropy, a measure of the entanentan-glement, scales with the boundary of the subsystem. In a random state, however, the entanglement typically scales with the volume of the subsystem. This special property of the low-energy states can be exploited
introduction tensor networks 4
elegantly in the language of tensor networks. A more detailed description of entanglement and its relation to tensor networks can be found in section 2.2.
Unrelated to the development of DMRG, a certain subclass of quantum states, Matrix Poduct States (MPS) were subject of interest in analytical studies. This class of possesses a particu-larly simple structure and are the most common example of a one dimensional tensor network. However, these MPS were con-sidered more of a conceptual tool rather than a powerful base for new algorithms until it was found that the DMRG algorithm could be elegantly restated in the context of MPS by Ostlund and Rommer [10].
2.1 Diagrammatic notation
At this point, it is best to start with the introduction of the tensors themselves. A tensor is a n-dimensional array of com-plex numbers. The order of the tensor describes the number of dimensions of the tensor: a order-0 tensor is simply a scalar, a order-1 tensor is a vector and a order-2 tensor is a matrix. Although these examples are easily visualized, the tensors can be of any order. A diagrammatic notation is then convenient to represent the tensors.
The pictorial notation that is used commonly in tensor network literature is defined in Figure 2.1. A tensor is denoted by a rounded blue shape and the number of legs correspond to the order. For example, a vector will have a single leg, a matrix has two and the state tensor of a system of N sites will have N legs.
vi1 i1 (a) Mi1,i2 i1 i2 (b) Ψi1,i2,i3,i4,i5 i1i2i3i4i5 (c)
Figure 2.1: Examples of tensor network diagrammatic notation. (a) A vector vi1 with a single index.
(b) A matrix Mi1,i2 with two
in-dices. (c) A tensor Ψi1, i2, i3, i4, i5 with five indices. Usually the explicit labels for the indices are omitted in large tensors.
By using this diagrammatic notation, complicated calculations can be efficiently represented, especially when the networks be-come larger. Therefore, from now on in this thesis, figures are used instead of explicit mathematical notation whenever pos-sible. In this way, the figures serve a similar purpose to the well-known Feynman diagrams in quantum field theory.
introduction tensor networks 5
Multiplication
Whenever two tensors are connected through their legs, multi-plication is implied. This operation for tensors can be cast in the form of matrix-matrix multiplication. An example is shown in Figure 2.2. The multiplication here can be written as:
Mi,j =
∑
k
Ai,k Bk,j (2.1)
which is simply the product of two matrices over one shared index.
Ai,k k Bk,j
i j
= Mi,j
i j
Figure 2.2: Example of a multipli-cation of two matrices in the MPS notation. The result of this operation is again a matrix.
This concept can be extended to arbitrary large networks of such tensors. In general, the connected legs of tensors represent multiplication, while the open legs determine the order of the resulting tensors (after all multiplications have been performed). The process of performing the multiplication over the connected legs in a network is referred to as contracting the network. One important remark is that there is a freedom in the order of contraction, when the calculation is represented by a diagram. This order could have an influence on the computational cost of contracting a network. Generally, determining the optimal contraction order of arbitraty networks is a difficult problem, especially in larger networks, although it can usually be deduced intuitively from a careful consideration of the corresponding fig-ure. An example is shown in Figure 2.3
A B C i l j k m n O(χ5) O(χ4) O(χ6) = i A (m× l) M j k
O(χ4) Figure 2.3: Example showing dif-ferent contraction options with cor-responding computational cost, as-suming all indices have dim∝ χ. In
this figure, contracting over index n first is the best choice. If the indices are of different sizes, another optimal contraction may be found.
The leading computational cost of a given contraction over an index between two tensors can be determined by multiplying the dimensions of all legs connected to both tensors, counting
introduction tensor networks 6
the connecting leg(s) once. In Figure 2.3, the best choice for the first contraction is to contract tensors B and C. Tensor B has indices k, l, n and tensor C has indices m, n. Therefore, the computational cost will be
cost∝ dim(k) · dim(l) · dim(n) · dim(m) ∝ χ4 (2.2) if all indices have dim∝ χ.
Reshaping
The tensors in a tensor network can be freely reshaped in their indices. When a tensor is reshaped, either two or more of its legs are grouped into a combined index, or a single leg is split into multiple legs. In Figure 2.4 this operation is pictured.
i
j k
m
= i
m
Figure 2.4: A tensor can be re-shaped without changing any infor-mation it contains. The combined leg is drawn with a thicker line to ex-press the fact that it is the result of grouping multiple legs together. This operation is also easily reversible, as long as dim(j)× dim(k) = dim(m).
The dimension of a combined leg j is the dimensions of its com-ponents’ legs iα combined:
dim(j) =∏
α
dim(iα) (2.3)
2.2
Matrix Product States
Since this thesis is focused on algorithms in one dimensional sys-tems, the most relevant example of a tensor network will be the
matrix product states2. In this section, an example of a ma- 2. In 2D, a natural extension are the Projected Entangled Pair States (PEPS) [11, 12]. These and other higher dimensional tensor networks lie outside the scope of this thesis.
trix product state is constructed from a general one dimensional quantum state.
A one-dimensional quantum state can be written in its most general form as
|Ψ⟩ = ∑
i1,...,in
Ψi1,...,1n|i1 . . . in⟩ (2.4)
where Ψi1,...,in represents all coefficients for the basis states,
writ-ten as configuration of spins. The summation is over all possible configurations of the n spins in the system. In this way, each possible state can be represented.
introduction tensor networks 7
Ψi1,...,in
i1 . . . in
Figure 2.5: The coefficients of a general one-dimension quantum state represented as a tensor. Note that the downpointing legs represent each site in the system.
The matrix product state is a way of representing the coefficients Ψi1,...,in in a different way. These coefficients can be viewed a
n-dimensional tensor, drawn in Figure 2.5.
Through a series of decompositions, detailed in the following paragraphs, the tensor can be put in the form of a matrix prod-uct state, shown in Figure 2.6. Each site now corresponds to a single tensor, and the full state is represented by the product over their combined indices. Note that an important distinc-tion between Figure 2.5 and Figure 2.6 is that the number of elements in the former scales exponentially with the number of sites, while the latter has a polynomial scaling. However, in ma-trix product states the exponential scaling can simply be hidden in the dimension of the combined indices, known as the bond dimensions, if these are allowed to grow to arbitrary size.
U1 U2 U3 U4 U5
Figure 2.6: An example of a Matrix Product State (MPS). Note that each site, corresponding to a physical index ik, has a separate tensor. The tensors are connected in a chain, where over the connecting indices multitplication is implied.
Singular Value Decomposition
An important operation involved in the process of finding an MPS representation of a given state is the Singular Value De-composition (SVD), which will be used throughout this thesis in multiple applications. An SVD is a common operation in lin-ear algebra which decomposes a given matrix into two unitary matrices U , V , and a diagonal matrix s, which has the singular values of the matrix on the diagonal:
M = U s V (2.5)
Such a factorization always exists and is related to an eigenvalue decomposition. However, the SVD is more general, since it can be applied to any m× n matrix.
introduction tensor networks 8 Ψi1,...,in i1i2 . . . in l = i1 M l = SVD U1 s V i1 m n l ˜ Ψi2,...,in
Figure 2.7: The first step of writing a state Ψ in MPS form. The left-most side is ’split’ off from the rest of the system using a Singular Value Decomposition (SVD). The bold indices are combined indices after reshaping a tensor.
Decomposing the state
To start the process of decomposing the state tensor, an SVD can be applied to the the tensor, where the tensor is first reshaped to a matrix. There is freedom in choosing which of the legs of the tensor to combine into which leg of the matrix, leading to different decompositions. Here we start with the left-most site, so we take the left-most leg of Ψ to be the left leg of the matrix, and group all other legs into the matrix’s right leg. Then, an SVD can be performed on the resulting matrix, resulting in two unitary matrices U and V , as well as a singular value matrix s. This first step of the process is depicted in Figure 2.7.
If we now absorb the singular values, contained in s into the right matrix V , the resulting matrix can be reshaped back into a tensor with n− 1 outgoing legs and one leg connecting it to the left matrix U . This tensor ˜Ψi2,...,in represents the rest of
the system, separately from site 1. In this way, the information about the state of the left-most spin is split from the rest of the system.
Entanglement
Naturally the question of how to deal with possible entanglement arises: it is in general not possible to write the information on site 1 separately from the information on the rest of the state. In the case where this is possible, the state is called separable or
a product state3. 3. not to be confused with the general
matrix product state
The entanglement is here contained in the dimensionality of the connecting leg: it represents how many combinations of the ba-sis states of the left-most spin and the baba-sis states of the rest of the system make up the complete state of the system. If this leg has dimension 1, there is no entanglement in the state. A product state consists of a chain of tensors connected with one-dimensional (trivial) legs. Clearly, the thickness of the leg has
introduction tensor networks 9 ˜ Ψi2,...,in i2 i3 . . . in k m l = SVD U2 s V m l ˜ Ψi2,...,in = U2 Ψ˜i3,...,in i3i4 . . . in k i2
Figure 2.8: The second step in decomposing the state into an MPS. Further steps are analogous. Note that the tensor on site 1 is not drawn.
a direct correlation with the amount of entanglement that can be represented by an MPS. Usually, a maximum is imposed on the dimension of a connecting leg, termed the bond dimension.
In this thesis, the bond dimension is denoted4 by D. The re- 4. Some authors denote the bond di-mension with χ. However, in the al-gorithm presented in this thesis, χ is used for a different parameter.
lation between bond dimension and entanglement is detailed in section 2.4.
Further decomposition
For the second site and further, the previous step can be re-peated. The difference is that now the large tensor ˜Ψi2,...,in has
an extra leg on the left side, corresponding to the right leg of the tensor that belongs to the first site. This leg is grouped with the i2 leg on the left, while indices i3, . . . , inare grouped on the
right side, forming another matrix. After performing a SVD, the i2 index is separated again from the auxiliary index, forming an
order-3 tensor on the second site. Again the singular values are absorbed on the right side, forming a new ˜Ψi3,...,in. This step is
pictured in Figure 2.8.
The process can be repeated until each site is represented by a single tensor, forming a chain of n tensors, connected by auxil-iary bonds. Such an configuration of tensors is termed a matrix product state and is shown in Figure 2.6.
2.3
Operators
In the language of matrix product states, also operators can be represented. An operator is usually drawn as a rectangular ob-ject, with legs pointing upwards as well as downwards. This illustrates the functional nature of operators: after applying an operator on a state, another state is returned. Naturally, oper-ators can have any number of pairs of legs. This thesis mostly
introduction tensor networks 10
deals with one- and two-site operators, e.g. spin and Hamilto-nian (resp.), or combinations of these. An example of the dia-grammatic representation of a one-site operator and a two-site operator is shown in Figure 2.9.
(a) (b)
Figure 2.9: Examples of operators that act on matrix product states.
(a) A one-site operator. Note that the two legs can be attached to one site of either an MPS ket, bra or both.
(b) A two-site operator, which can be placed on two contiguous sites of an MPS.
Expectation values
To compute the expectation value of a state written as an MPS, an operator is placed in between the state and its hermitian conjugate. The hermitian conjugate of a tensor is usually drawn in opposite orientation of the tensor itself. Using this convention, an example of an expectation value is shown in Figure 2.10. Legs that connect tensors of the MPS and the hermitian conjugate MPS directly, without passing through an operator can be seen as the identity operator. Mathematically, this corresponds to implicitly assuming a tensor product of the operator with the identity operators in the spaces it does not cover itself.
Figure 2.10: Computation of an expectation value of an operator (represented by the red square). On legs where the operator is not present, the identity operator can be read.
Besides the operators that correspond to physical observables, which can be placed between an MPS bra and ket to compute its expectation value, other types of operators can be found in tensor networks. These operators can have arbitrary shapes and number of legs. Examples of such operators are introduced in section 2.5.
2.4
Entanglement
The matrix product state is an efficient ansatz for quantum states that adhere to the area law of entanglement. The max-imal entanglement that can be represented in an MPS of bond dimension D can be calculated by considering a block of L sites being separated from its environment on either side. If both
ex-introduction tensor networks 11
ternal indices have dim(ilef t) = dim(iright) = D, the combined
index ¯i ={ilef t, iright} of this block must have dim (¯i) = D2.
If we write the complete state|Ψ⟩ as
|Ψ⟩ =
D2
∑
¯i
|block⟩ ⊗ |environment⟩ (2.6)
the rank of the reduced density matrix ρblock is at most D2,
depending on the entanglement in the state5 [13]. 5. Note: the rank of the reduced den-sity matrix ρenvironmentis also at most D2.
The Von Neumann entanglement entropy, one measure of the entanglement in the system, is given by
S(ρblock) =−tr (ρblock log(ρblock) ) (2.7)
Since ρblock has at most rank D2, it has at most D2 eigenvalues.
If the state is maximally entangled, these eigenvalues will be equal: λi = D−2. The entropy for ρblock written in diagonal
form takes the simplified form, assuming the state is maximally entangled: S(ρblock) =− D2 ∑ i 1 D2log 1 D2 = logD 2= 2logD (2.8)
Note that this value is independent of the block size L. Thus, the entanglement of a block of any size L in any matrix product state of bond dimension D is bounded from above by
S(L)≤ 2logD (2.9)
which satisfies the area law for states it can represent: S(L)≤ αLd−1 - a constant in one dimension.
2.5
Multiscale Entanglement Renormalization
Ansatz (MERA)
Besides matrix product states, a second type of tensor network is of great importance to this thesis. The Multiscale Entangle-ment Renormalization Ansatz (MERA) was introduced by [14– 17] and has been very successful in representing the long-range correlations in critical quantum states. In these states, the MPS ansatz suffers from an inherent limitation in encoding long-range correlations, since it can only represent exponentially decaying
introduction tensor networks 12
correlation functions on larger scale [6, 7], due to its inability to reproduce the entropy scaling S ∝ log(L) that is typical for critical states. Being fundamentally different from MPS, MERA can display the correct polynomial decay of correlation functions by construction.
MERA is closely related to the Renormalization Group tech-niques [2, 18, 19], which have been extremely influential in mod-ern condensed matter physics, and produces a proper renormal-ization flow. Also scale invariance (for critical states) can be imposed, making it possible for MERA to reproduce critical (un-stable) fixed points.
Coarse grainers
At the basis of MERA lies the coarse grainer. This operator has the task of projecting the vector spaces of multiple sites in an MPS onto a combined smaller space. A coarse grainer is defined as an isometry
w :V′ −→ V⊗m (2.10)
with dim(V′)≤ dim(V⊗m) and m corresponds to the number of sites that is coarse grained into a single supersite. Various types of coarse-graining schemes, such as two-to-one or three-to-one can be defined in this way. Since w is an isometry, the following properties must hold true:
w†w =1, ww†= P, P2= P (2.11)
A diagrammatic representation of a three-to-one coarse grainer, which will be the type used in this thesis, is shown in Figure 2.11.
(a) (b)
Figure 2.11: Notation of a three-to-one coarse grainer. (a) A coarse grainer by itself, separate from any connecting tensors. (b) The coarse grainer attached to a three-site MPS. Note that the single upward leg rep-resents the (compressed) combined vector space of the three sites of the MPS.
Disentanglers
A fundamental problem arises when a MERA is constructed us-ing only coarse grainers. Although the coarse grainers effectively
introduction tensor networks 13
Figure 2.12: Example of a MERA using only (three-to-one) coarse grainers (grey triangles). On the central bond, short-range entanglement (if present), illustrated by the red wavy line, is left untouched by coarse graining. Any configuration of the coarse grainers will always leave some short-range entanglement in the system.
remove short-range entanglement in the block of sites it is ap-plied to, it leaves entanglement across the borders of the blocks intact. Some of the entanglement present in the microscopic lat-tice will therefore remain in the system, even after many coarse graining steps. This situation is illustrated in Figure 2.12. To solve this problem, the concept of disentanglers was intro-duced [14]. These are unitary operators that can be placed on a block of sites, usually two sites wide, to effectively remove the short-range entanglement inside the block. Since the opera-tor is unitary, by construction, the transformation is reversible and leaves any expectation values invariant. The disentangler is defined as
u :V⊗m−→ V⊗m (2.12)
with the property
u†u = uu†=1 (2.13)
In this thesis, the disentanglers that are used are exclusively two-site (m = 2 in Equation 2.12) operators. An important remark is that, since the operator has a limited coverage of sites, not all entanglement can be removed. This would require a disentangler of the same size as the entire system. In practice, large disentanglers are computationally expensive to construct, as well as calculations involving them. It is usually more efficient to combine several smaller disentanglers, at the cost of leaving a small amount of entanglement in the system [20].
The construction of the disentanglers is usually achieved by opti-mizing the operator against a measure of entanglement, such as
introduction tensor networks 14
the von Neumann entanglement entropy, under the constraint of unitarity. An explicit method for finding optimal disentanglers is described in section 3.4.
An illustration of a disentangler, as used throughout the rest of this thesis, can be found in Figure 2.13.
(a) (b)
Figure 2.13: Pictorial representa-tion of disentanglers. (a) An isolated disentangler, as a two-site opera-tor. (b) Placement of disentanglers in MERA (network cut off at both sides). Note that each bond on the microscopic lattice is either within a coarse grained block or under a disentangler, ensuring a proper renormalization flow in terms of en-tanglement.
2.6 Application - Transverse Ising model
Tensor networks can represent a broad class of quantum and classical states, making them applicable to a wide range of phys-ical models. As mentioned earlier in this chapter, tensor net-works are especially powerful for states that obey the area law of entanglement. As this thesis focuses on methods to simulate models with ground states for which this does not (always) hold, such as critical models, one obvious choice for a model to serve as an example for applications would be the well-know quan-tum transverse Ising model in one dimension. Its Hamiltonian is given by
H =−J ∑
i,i+1
σxiσxi+1− h σz (2.14)
where J ∈ {+1, −1} toggles between the ferromagnetic and anti-ferromagnetic variants, h ≥ 0 represents the strength of the transverse field and σk are the Pauli spin matrices. This model,
though seemingly simple, has been studied widely for its surpris-ingly rich physics. With h = 0, it reduces to the simple Ising model, which is considered to be the most basic model to ex-hibit a phase transition in D > 1 dimensions. In one dimension, there is an absence of a classical phase transition, however at zero temperature it does show a quantum phase transition, at the critical point h = 1.
Precisely at this critical point, the ground state of the model is a highly entangled state, which has correlations between spins at all length scales, classified by a diverging correlation length.
introduction tensor networks 15
Although this type of state is not ideally represented by a typi-cal one-dimensional tensor network, specifitypi-cally a matrix prod-uct state (MPS), it has the interesting property of being scale-invariant [1, 21]. This can be exploited by a MERA ansatz, described in the previous sections, so it serves as an excellent example for applying the scheme proposed in this thesis, which aims to overcome the limitation in matrix product states by renormalizing the state by techniques from MERA.
3
Methods
3.1 Introduction
BackgroundThe focus of this thesis lies on a new algorithm that consists of a blend of several concepts in tensor networks. The field of tensor networks has sparked the development of many different ground state algorithms, which usually target one aspect specifically yet are less effective in other aspects. For example, the well-known infinite version of the DMRG algorithm, which is used in one part of this thesis, has been very effective at quickly converging to a state with an energy that is very close to the exact ground state energy. As described in detail later in this chapter, infinite DMRG produces a state of an infinite system by assuming a structure of repeating unit cells. Such unit cells form the basis of infinite MPS in many algorithms and can simulate a quantum state up to very high accuracy.
However, while infinite DMRG is very accurate in terms of en-ergy, the true long range physics are in some cases beyond its reach. This is not only an effect of the DMRG itself, since the problem lies in the underlying MPS ansatz and its intrinsic in-ability to reproduce highly entangled states. Although this is by design and actually the basis to its success, due to the area
law of entanglement in most ground states1, difficulties arise in 1. The ability of MPS to represent such states is described in section 2.4
critical systems. The ground state of the transverse Ising model that serves as an example throughout this thesis is, when the system is at criticality, entangled over all length scales. This inherent weakness to critical systems of the DMRG algorithm is one of the main motivations for both MERA and the combined algorithm presented in this thesis.
The MERA anzats follows an entirely different approach than MPS by aiming to effectively renormalize the entanglement in the system. An important property of the MERA scheme is that the distance between sites decreases exponentially in the
methods 17
ber of coarse graining steps: in the three-to-one scheme used in this thesis, two sites that are L sites apart become nearest neighbours after only log3(L) coarse graining steps. Because of this property, MERA is able to accurately reproduce poly-nomially decaying correlation functions [8], present in critical ground states. By comparison, any MPS (and therefore any re-sult of DMRG) is only capable of displaying exponential decay in correlations [7] at large length scales. Despite this powerful advantage, however, a MERA is constructed by iteratively op-timizing each tensor in its usually rather large network, making it quite slow to converge.
Evidently, an ideal approach would seem to be a combination of these ideas: an algorithm that seeks to accurately reproduce long range physics, while being faster than conventional MERA. Such an approach is tested in this thesis, which does not involve optimization of a large MERA network. Instead it builds the network layer by layer using the ground state, produced by a version of infinite DMRG, to compute the coarse grainers along the way. The details of this procedure are described in the re-mainder of this chapter.
Algorithm
The algorithm consists of two steps: an energy minimization step and a coarse graining step. For the minimization step, a form of the infinite DMRG method is used to find the ground state of the matrix product state. In the coarse graining step, following the concepts of MERA, a combination of coarse graining opera-tors and disentanglers are utilized to “lift” the Hamiltonian, the state and any other relevant operators onto the next lattice. On this lattice, the energy minimization step is again performed, and the cycle continues. This leads to an effective version of the Hamiltonian, which corresponds to a certain coarse grained lattice.
3.2
Energy minimization (infinite DMRG)
In order to find the approximate ground state to a Hamiltonian on a given lattice, a version of the infinite DMRG algorithm
is implemented. In this algorithm, a small starting MPS2 is 2. The starting n-site state can be constructed by finding the lowest eigenvector of the n-site Hamiltonian. This is possible only for small sizes, such as four or six sites.
grown by iteratively inserting new sites in the center. Each new added block of sites is constructed in such a way as to minimize the energy of the full system. Since each new contribution is
methods 18
optimized in this way, the resulting growing block of sites is made to be an approximate ground state.
The energy mimization is derived from the optimization of a state|Ψ⟩:
∂
∂Ψ†(⟨Ψ| H |Ψ⟩ − λ ⟨Ψ|Ψ⟩) = 0 (3.1) making use of a Lagrange multiplier to ensure proper normal-ization of the state. DMRG and most other tensor network algorithms rely on performing the energy minimization on only small parts of the tensor network at a time. In each update step, one or more tensors are selected and optimized with respect to their environment, consisting of all other tensors in the network, which is treated as constant during that step. The new tensors are then used in the environment for the selected tensors in the next step. This process is repeated until the full state converges to minimum.
Initialization
The first step in starting up the infinite DMRG algorithm is to compute a ground state of a finite system. This involves constructing a total Hamiltonian for this system, by combining multiple two-site Hamiltonians. The choice of starting system size is free, although small in most cases due to the exponential scaling of the corresponding eigenvalue problem. For a four-site starting system, the Hamiltonian is constructed as:
H4= Hi1,i2⊗ 1i3⊗ 1i4 +1i1 ⊗ Hi2,i3 ⊗ 1i4
+1i1 ⊗ 1i2 ⊗ Hi3,i4
(3.2)
Once a lowest eigenvector has been found, using any iterative method, it can be reshaped into a state tensor. Using the pro-cedure described in section 2.2, the state is transformed into a matrix product state.
Two main components of the algorithm are the left and right environment blocks. These blocks, defined in Figure 3.1, contain a part of the environment needed in the following steps. First, the ground state MPS of the finite system is cut in the center, forming two chains of equal length. In the initialization step, the blocks are built from the expectation value computation of the left and right halves of the MPS. However, the trace over the
methods 19
bond that was cut is not performed, leaving a matrix for each half. L = . . . . . . + . . . . . . + . . . . . .
Figure 3.1: The left environment block used in the infinite DMRG algorithm. The two-site Hamilto-nian operator is applied to each pair of consecutive sites in the chain. Note that the resulting block has the shape of a matrix.
The sites replaced by ’. . . ’ can be understood to have the iden-tity operator applied. The length of these omitted blocks of sites depends on the duration of the simulation thusfar, growing with each DMRG step. A right environment block is constructed analogously, forming another matrix.
Main steps
The main part of the algorithm consists of inserting new blocks of sites, optimized to retain a ground state in a growing system. 1. (Figure 3.4) The full environment E of the new tensor is calculated, using the left and right environment blocks. This environment, which is a O(χ4D4) tensor, can be
re-shaped into a χ2D2× χ2D2 matrix.
2. This step involves finding the two-site block, denoted by e
T, which minimizes the expectation value of the energy. This expectation value is given by
E = eT†· E · eT (3.3)
From Equation 3.1 follows that finding this optimal two-site block is equivalent to solving the eigenvalue problem
E · eT = λ eT (3.4)
in which the environmentE is represented as a matrix and the two-site block as a vector. If the growing MPS is to retain its ground state properties, λ is required to be the lowest eigenvalue.
3. When the lowest eigenvalue λ and corresponding eigenvec-tor eT is found, it can be split into two site tensors U , V
methods 20
using an SVD (Figure 3.2). The new tensors are absorbed into the left and right environment blocks (Figure 3.3), to be used in the next step.
The weights λi can be discarded3, and due to the unitar- 3. They will be replaced by a new
two-site block tensor in the next step. The weights of consecutive steps can be useful to store for monitoring con-vergence.
ity of the tensors Ui and Vi, the MPS will remain in mixed
canonical form. A convenience of this canonical form, aside from the simplification in the computation of the terms de-picted in Figure 3.4, is that the singular values λi represent
the Schmidt coefficients of the full system. Therefore, the truncation in the dimension of the center tensors down to the predefined bond dimension D is guaranteed to be optimal.
T =
SVD
Ui Λi Vi
Figure 3.2: The optimal new two-site block eT is split into a left and
right tensor (by SVD), keeping the weights in the center.
Li = Li−1 Ui Ui† + Ui−1 Ui†−1 Ui Ui† H
Figure 3.3: Update of the left and right environment blocks with the new tensors Uiand Vi, at step i.
Convergence
After a sufficient number of repetitions of steps 1-3 above, the algorithm will converge, with eTi+1 very close to eTi. Effectively,
this means that the center of the current state, which is now 2i sites long, is close to the center of the infinite ground state. Its energy will be close to the ground state energy of an infinite system. Several measures for the convergence have been sug-gested [3, 22], which is a delicate affair. In our case, the overlap of the states of consecutive steps was used.
Infinite ground state
After repeating above steps N times, the system will start to resemble an infinite system. One important difference between the infinite DMRG algorithm as implemented in this thesis and other versions is the fact that the Hamiltonian operator has an even-odd structure. This structure, in which the two-site Hamil-tonian operator on even bonds in the MPS is slightly different
methods 21 L T T† + (1) T T† + (2) T T† + (3) T T† + (4) T T† R (5) ≡ T T† E
Figure 3.4: Computation of the full environment in a DMRG step. The tensor denoted by T is the newly added two-site block, while the dotted line around T† signifies the empty slot for the Hermitian conjugate of T: taking the derivative with respect to T† corresponds to removing the tensor from the network. The result of the addition is the full environment E, which can be reshaped into a matrix. The lowest
eigenvector T∗ is the optimal choice for the newly added sites.
Note that considerable use has been made of the canonical form of the state: in (1), the right side is reduced to an identity; similarly the left side in (5) and both sides in (2), (3), (4).
. . . UA UB UA UB Λ VA VB VA VB . . .
Figure 3.5: Center section of the MPS generated by the infinite DMRG algorithm, with A and B repre-senting the even-odd structure of the state. The arrows symbolize the left- or right-normalization. Evidently the MPS is in a mixed-canonical form with respect to the central Λ.
from the one on odd bonds, a consequence of the coarse graining
element of the full algorithm4. To account for this this even-odd 4. An explanation of this fact is de-tailed in section 3.4.
structure of the Hamiltonian, DMRG simply alternates between the two versions in step 1 of the repeating segment of the algo-rithm.
The resulting MPS will then have a similar structure of even and odd sites, illustrated in Figure 3.5. Yet the formed MPS has no clearly discernible unit cell, since the central Λ element only appears once. In [7] and [22], a simple solution is described. This solution can be easily derived by restating the MPS in Γλ-notation, devised for the iTEBD algorithm [4]. Alternatively, and perhaps more intuitively, it can be derived explicitely from the notation used in this chapter.
At a certain step i, the center block of the MPS is given by
methods 22
where the U and V tensors are respectively left- and right nor-malized. To speed up the eigensolver for the next step, it is possible to make a guess for the new center block.
Inserting a new two-site block in the center, which will then be decomposed by SVD, will result in
. . . UiA λi ViB. . . → . . . UiA
{
Ui+1B λi+1 Vi+1A
}
ViB. . . (3.6) In an MPS with a two-site unit cell, the expectation values of any operator is translationally invariant up to even or odd sites. Due to the canonical form of the MPS at each step in the algorithm, it is very easy to compute the reduced density matrix of a single site. By requiring the reduced density matrix of site A of step i + 1 to be close to that of site A of step i (equivalently for site B), the following guess can be made:
e
Ui+1B eλi+1= λi ViB
eλi+1 Vei+1A = UiA λi
(3.7)
The full guess for step i + 1 will then read, inserting an identity in the second step,
. . . UiA {
e
Ui+1B eλi+1 Vei+1A
}
ViB. . . = . . . UiA
{ e
Ui+1B eλi+1
( eλi+1
)−1
eλi+1 Vei+1A
} ViB. . . = . . . UiA { λi ViB ( eλi+1 )−1 UiA λi } ViB. . . (3.8)
which still requires an explicit (eλi+1)−1 in the center. As the algorithm converges, the matrices λk will converge as well. Due to (small) differences in the even and odd sites in the MPS, a good choice will be to set eλi+1= λi−1, since both were formed as
the center matrix with respect to the same (even or odd) bond. Additionally, the identifications in Equation 3.7, which will ap-proach exactness as the algorithm converges at some step k, allows us to properly identify a unit cell for the state in the thermodynamic limit. After convergence, the tensors can be as-sumed to be equal, up to even/odd structure and left- or right normalization, resulting in a state as follows:
methods 23
. . . UA UB UA λ VB VA VB. . . (3.9) where the center block does not repeat itself. However, we can rewrite the two-site blocks on the left and right side as
UAUB = UA UB λk−1 (λk−1)−1= UAλk VB (λk−1)−1
VA VB= (λk−1)−1 λk−1 VA VB= (λk−1)−1 UA λk VB
(3.10) with λk, λk−1 the weight matrices of the last two steps of the
algorithm. Either of the blocks on the right hand side of Equa-tion 3.10 can be chosen as a unit cell for the ground state in the thermodynamic limit. A simple absorption of the weight matrices into the tensors results in an MPS in the familiar A, B form:
A B ≡ (UAλk) (VB (λk−1)−1) (3.11)
Note that the left- and right normalization properties of the A and B tensors are not made explicit anymore, since the unit cell in its entirety does not have such a property. The state can be brought into a canonical form by usual orthogonalization methods [7, 23].
3.3
Alternative ground state algorithms
Imaginary time evolution (iTEBD)
One popular alternative to the infinite DMRG method was
con-ceived in 2007 [4], named infinite5 Time Evolving Block Deci- 5. This infinite version was based on the original finite TEBD algorithm.
mation, or iTEBD. Instead of optimizing the center tensors in a growing lattice with respect to the approximate environment through finding the lowest eigenvector, iTEBD takes a quite dif-ferent approach.
Any trial wavefunction|Ψtrial⟩ that has a finite overlap with the
ground state6, can be written in terms of the eigenstates of the 6. This should be satisfied by a ran-dom trial wavefunction, although a better guess leads to faster conver-gence. Hamiltonian: |Ψtrial⟩ = ∑ i λi|i⟩ (3.12)
methods 24
A time evolution of such a pure state can be facilitated by a unitary time evolution operator, which depends on the Hamilto-nian:
U (t) = e−iHt (3.13)
This operator in real time can be transformed into an operator that evolves the state in imaginary time, by setting τ = it, τ ∈ R. Note that the imaginary time evolution operator, written as
e
U (τ ) = e−Hτ (3.14)
The application of this operator on the state |Ψtrial⟩ can be
described in terms of the energy eigenstates:
e U (τ )|Ψtrial⟩ = ∑ i λi e−Hτ|i⟩ =∑ i λi e−Eiτ|i⟩ = λ0 e−E0τ|0⟩ + ∑ i̸=0 λi e−Eiτ|i⟩ (3.15)
If we multiply both sides with a factor e+E0τ and take the limit
τ → ∞, we get eE0τ U (τ )e |Ψ trial⟩ = λ0|0⟩ + ∑ i̸=0 λi e <0 z }| { (E0− Ei)τ|i⟩ τ→∞ −−−→ λ0|0⟩ (3.16)
The iTEBD algorithm is constructed by performing a Suzuki-Trotter decomposition on the imaginary time evolution operator
e
U (τ ) into operators acting on two sites [4, 24] for a very small time step δτ ≪ 1: e Uk,k+1(τ )≈ ( e−hk,k+1 δτ )N ≡(Uek,k+1 )N , τ = N· δτ (3.17)
with hk,k+1the Hamiltonian term acting on sites k, k + 1. As for general Hamiltonians the local terms do not commute, this de-composition is not exact, but will have an error ofO(δτ). Higher order schemes are possible to further suppress this Trotter error:
methods 25
in this thesis, a second order decomposition is used, resulting in an error ofO(δτ2).
In its original formulation, the iTEBD algorithm works directly in the thermodynamic limit. In the case of a local two-body Hamiltonian, such as the Ising Hamiltonian, the terms can be separated into even and odd terms, working on even and odd pairs of sites. Although these terms generally do not commute, all terms within the even/odd class commute. This results in a scheme in which a two-site unit cell is updated on alternatively the even and odd bonds.
Recently McCulloch [22] noted that the iTEBD algorithm can be rewritten in the form of the iDMRG scheme, by replacing the eigenvalue problem in Equation 3.4 by a time evolution step:
e
T = eUk,k+1 Tprev (3.18)
with Tprev generated from the previous step (Equation 3.8). As
in iDMRG, this new two-site tensor eT is split by SVD, result-ing in the new Ui, Vi tensors. Besides this minor modification,
the rest of the algorithm can remain unchanged. However, since the update step in iTEBD does not depend on the Hamiltonian terms of the environmentE, computed as in Figure 3.4, its com-putation can be omitted, greatly simplifying and speeding up the algorithm.
Uniform MPS
The iDMRG algorithm has been studied and used extensively, due to its high accuracy and relatively fast convergence. In its most popular implementations, also in this thesis, it is based on a two-site unit cell. The advantage over a one-site implementa-tion is two fold: unlike the one-site scheme, the multi-site scheme allows for systematic truncation of the bond dimension at each update step, as well as a variable bond dimension size, and sec-ondly the one-site scheme has been found to converge toward local minima in some cases. This problem is usually resolved when using a multi-site unit cell.
In some models, however, a single site unit cell may be preferred over larger sizes. One recent approach, named Variational Uni-form Matrix Product States (VUMPS) [5], attempts to provide a more robust alternative to the one-site iDMRG scheme. It is similar to iDMRG, but has one significant difference: instead of growing the finite lattice by inserting one unit cell at a time,
methods 26
while computing the environment termsE by reusing the terms of the previous step, VUMPS works in the thermodynamic limit from the start. In this algorithm, the environment terms are computed by explicitly performing the infinite summation until convergence.
Symmetries
One additional advantage of the one-site VUMPS algorithm is that it allows for a very straightforward implementation of sym-metries. In the case of the quantum transverse Ising model, the behaviour of the ground state is divided into a symmetric and disordered phase. For a transverse field strength h < 1, the ground state is two-fold degenerate, but the symmetry is spon-taneously broken in favor of one the options. At h ≥ 1, the ground state becomes disordered.
In order to conserve this symmetry in the state, the site tensor in a uniform MPS can be split into even and odd sectors [25, 26]. If we adopt a convention in which we assign
|↑⟩ : + (even), |↓⟩ : − (odd) (3.19) the tensors should have the following structure:
A+,m,n = (+ − + U ∅ − ∅ V ) , A−,m,n= (+ − + ∅ X − Y ∅ ) (3.20)
Where + (−) denote the even (odd) sectors of the indices. This arrangement of indices is possible, since the order of auxiliary (bond) indices is arbitrary. At the start of the VUMPS algo-rithm, the tensor can be initiated in this form, fixing the
struc-ture of the auxiliary index (up to the even/odd blocks)7. 7. Special care must be taken at the start, when the initial state is put into a canonical form.
Modifications to infinite sums
At the base of the VUMPS algorithm lies the evaluation of infi-nite summations of Hamiltonian terms, to compute an effective Hamiltonian environment for the center tensors. These infinite sums also arise in the iDMRG scheme, however they are con-structed in an iterative way as the lattice grows. The compu-tation of these terms are shown in Figure 3.3. In the VUMPS scheme, these terms are computed directly in the infinite lattice.
methods 27
The full computation can be written in terms of the transfer matrix [7]T : Lk= L0 · { 1 + T + T2+ . . .} = L0 · k ∑ n=0 Tn = L0 · { N + k ∑ n=0 e Tn } (3.21)
with the convention thatT0≡ 1 and k → ∞. In the third line a
diverging partN of the infinite summation is substracted. Usu-ally, the transfer matrix will have a dominant eigenvector with eigenvalue 1, assuming a normalized state. The divergent part will then simply be the projection onto the dominant eigenspace:
T =∑ i λiPi = P0+ ∑ i>0 λiPi (3.22)
However, in the symmetry-preserving case, simulations in the ordered phase h < 1 of the Ising model will have a two-fold de-generacy in their ground state. Therefore, a slight modification is necessary in the scheme to accomodate for the extra diverg-ing part in the transfer matrix. Assumdiverg-ing a double degeneracy, although this derivation can be easily extended to higher degen-erate cases, the transfer matrix has now two dominant left- and right eigenvectors with eigenvalue 1. Its spectral decomposition will then be T =∑ i λiPi= P0+ P1+ ∑ i>1 λiPi (3.23)
Following the steps from [5], the infinite sum takes the form
∞ ∑ n=0 Tn= ∞ ∑ n=0 { P0+ P1+ ∑ i>1 (λi)nPi } =N +∑ i>1 (1− λi)−1Pi (3.24)
Where we denote the diverging term coming from P0 and P1 by
N and recognize the non-diverging terms as geometric series in λi. The rest of the steps are equivalent to the non-degenerate
methods 28
lim
k→∞Lk· {1 − T + P0+ P1} = L0· {1 − P0− P1} (3.25)
Which is of the form x· A = b, and can be solved using various iterative methods. The computation of the Hamiltonian terms on the right side of the infinite chain is equivalent to the above procedure.
3.4 Coarse graining and disentangling
Once the ground state of the infinite system has been obtained to a sufficient accuracy, the second step can be initiated. This consists of constructing the necessary coarse grainers and disen-tanglers, and applying these to compose an effective Hamiltonian on a larger scale lattice.
In this section, a three-to-one coarse graining scheme is assumed. This scheme has the convenient property that any two-site op-erator remains as such after coarse graining and disentangling, regardless of placement. A careful consideration should be given to the even-odd structure of the ground state MPS provided by
the DMRG algorithm8. In Figure 3.6 the procedure is detailed. 8. In our case, the DMRG algorithm is based on a two-site unit cell. Other unit cell sizes are possible and have been studied [22, 27, 28]. The follow-ing procedure should be adapted to ac-commodate other unit cell choices.
Note that the disentanglers are placed on the boundaries of the three-site blocks that are to be coarse grained. Evidently, two separate contexts can be identified in which the disentanglers are placed: either an AB block or a BA block. This requires two different disentanglers, that are to be calculated separately. Similarly, it follows that also two different coarse grainers are needed.
Figure 3.6: Full coarse graining and disentangling of the Hamiltonian operator. All three possible place-ments of the operator are taken into account.
Note the even-odd structure off all tensors, symbolized by different shades.
Disentanglers
A disentangler is a unitary operator which serves to remove as much entanglement as possible across a given boundary. On this boundary, the system can be cut in two parts. In our case, the system is infinitely long with a two-site unit cell. This leaves two choices for placing the boundary: on an even or an odd bond.
methods 29
On one of the two parts the reduced density matrix can be cal-culated by tracing out the other part of the system. If the MPS is in a canonical form, this partial trace reduces to an identity operation. The result of this procedure is shown in Figure 3.7.
. . .
Figure 3.7: Reduced density matrix of the left half of an infinite MPS. If the right side of the boundary is right-normalized, the partial trace of this entire part of the system reduces to an identity.
Several methods can be employed in order to find the optimal unitary operator to use as a disentangler, . Most of these meth-ods rely on a variational optimization of the operator against a cost function, which should be some measure of entanglement. In our case, this will be the second Rényi entropy, defined as
H2(x) =−log tr(ρ2) (3.26)
which involves the squared reduced density matrix. Although the trace of ρ itself cannot depend on any unitary transforma-tion of the state, the trace of ρ2 will change. The right
dis-entangler should aim to minimize this value. However, when applying any unitary transformation|Ψ⟩ → U |Ψ⟩, the operator U appears four times in the calculation of the Rényi entropy. A variational optimization of all these at once would be com-putationally expensive, therefore a simpler iterative method is utilized.
The iterative method for optimizing a disentangler can be
initi-ated from either a random unitary tensor9 or a specific choice, 9. Such a tensor can be obtained from the singular value decomposition of a random matrix. One of the resulting unitary matrices can then be reshaped into a tensor of the required size.
such as the identity. Then the following steps are repeated from the initialU until convergence:
1. (Figure 3.8) Update the environment: compute ρ2with the
current U applied except for one location in the network. This will serve as the environment, after reshaping into a matrix M.
2. The choice of the newU that minimizes the cost function
methods 30
can be deduced from the singular value decomposition M = VsW†, by using the cyclical property of the trace:
tr(UM) = tr(UVsW†) = tr(W†UVs) ≡ tr(As) (3.28) where A is unitary by construction. Since s is a strictly positive diagonal matrix, the choice of a unitary matrix A that minimizes the cost function is
A =−1 =⇒ U = −WV† (3.29)
leading to tr(UM) = −∑isi, with si the singular values
that appear on the diagonal of s.
= M
Figure 3.8: Full computation of
tr(ρ2), needed for the Rényi entan-glement entropy. The location of the unitary tensorU that is calculated in
step 2 of the method is designated by the dotted shape.
The environment can be reshaped into a matrixM, shown on the right.
Coarse grainers
After the disentanglers have been constructed, the coarse grain-ing operators can be computed. These operators map a given state, in this case the ground state obtained by DMRG, into a coarse grained state. Additionally, they map any operators into the corresponding coarse grained operators. On this new effec-tive lattice, each site represents a block of sites on the previous lattice.
An important part in the role of the coarse grainers is the selec-tion of the subspace onto which it projects. This subspace should be small enough to make DMRG feasible, yet large enough to re-tain the most relevant information of the state. In our case, the
methods 31
= =
Figure 3.9: Full calculation of a coarse grainer. The three central sites form the relevant block, while the sites on the far left and right are traced out, resulting in the reduced density matrix of the three-site block,
after disentangling.
The block can be reshaped into a matrix and the coarse grainer can be extracted through an eigendecom-position.
most important coarse graining operation involves the nian operator and is detailed in the next section. This Hamilto-nian should become an effective operator on a large scale lattice, containing only the most relevant information.
The optimal choice of subspace is the eigenspace of the ground state density matrix [3], since the coarse grainer should retain the support of the ground state reduced density matrix. Figure 3.6 shows that in our case the relevant block for a coarse grainer involves one center site and two disentanglers, spanning five sites
in total10. The reduced density matrix for the three-site block 10. The actual coarse grainer handles only three site, but the disentanglers on the boundaries involve the neigh-bouring sites as well.
that the coarse grainer should cover is obtained by a partial trace over the far left and right sites, drawn in Figure 3.9.
The resulting block can be interpreted as the reduced density matrix of the three sites, after applying disentanglers on the boundaries. From this density matrix, which is Hermitian, the coarse grainer can be constructed from its spectral decomposi-tion: ρ = χ ∑ i λi|Ψi⟩ ⟨Ψi| =⇒ ww†= χ ∑ i |Ψi⟩ ⟨Ψi| (3.30)
where λi are the χ nonzero eigenvalues11. The matrix w can 11. Note that in practical applica-tions, the number of eigenstates that are kept is fixed in such a way that the remaining λi>χare small, so that the support of the ground state den-sity matrix is retained approximately.
then be reshaped into the order-4 tensor shown in Figure 3.9.
Uniform states
When the chosen ground state algorithm produces a uniform ma-trix product state (one-site unit cell), such as with the VUMPS scheme earlier in this chaptersection 3.3, the coarse graining scheme simplifies. Since there is no even/odd variation in the
methods 32
tensors, both disentanglers and both coarse grainers will be equal.
Conserving symmetry
If symmetry is enforced in the tensors during the uniform (VUMPS) algorithm, care must be taken when computing coarse grainers. As detailed in section 2.5, the coarse graining operators are based on the ground state reduced density matrix. The implementa-tion of spin-flip symmetry in the x direcimplementa-tion, corresponding to the Ising Hamiltonian, as described in section 3.3, the reduced density matrix of a three-site block will have a certain structure. A complication arises when we perform the eigendecomposition on the reduced density matrix, to form the coarse grainer, since the order of the eigenvalues can be chosen freely and therefore the order of the rows in the matrix S containing the left eigen-vectors in the decomposition M = SDS−1 can vary. One way around this problem would be to reorder the reduced density matrix based on the parities of the physical indices, resulting in a block-diagonal structure. This is always possible, since the symmetry conservation implies that half of the elements in the reduced density matrix will be zero.
If we define a vector Plfor each site, belonging to a coarse grained
lattice l, which contains the relevant parities corresponding to the physical indices, the parities of the combined physical indices of the three site block are
Pl+1= Pl⊗ Pl⊗ Pl (3.31)
where⊗ denotes the Kronecker product. In the convention that on the microscopic lattice we define P0 =
(
+1−1), the parities
on the first coarse grained lattice l = 1 will read
P1=
(
+1 −1 −1 +1 −1 +1 +1 −1
)
(3.32) We can now compute the permutations required to sort Plin the
order Pl = (+1 ... +1−1 ... −1). Note that the order within the
even/odd sectors is not important. Since the symmetry enforce-ment ensures that in the reduced density matrix all eleenforce-ments within odd blocks are zero, we can permute both indices of the reduced density matrix in the same order as Plto obtain a block-diagonal form:
