On Tensor Network Methods for One-Dimensional Open Quantum Systems

MSc Dissertation

On Tensor Network Methods for

One-Dimensional Open

Quantum Systems

by

Karel D. Temmink

Submitted in partial fulfilment of the requirements for the degree

of Master of Science

University of Amsterdam

Primary supervisor:

Dr Philippe Corboz

IOP

Secondary supervisor:

Dr Vladimir Gritsev

IOP

A B S T R A C T

Presently, Tensor Network (TN) methods have firmly established themselves as re-liable, efficient, and extremely powerful tools for quantum calculations.TNs have proven specially successful in regimes where the time-evolution is unitary and/or en-tropy obeys an area law, such as ground state calculations in closed quantum systems.

However, less has been accomplished for open quantum systems, where the time-evolution generated by the Lindblad master equation is no longer unitary. These systems, which are analytically relatively poorly understood, are also notorious in computational physics, as they tend to cause all sorts of numerical issues, with the most well-known being that simulated density operators often lose positivity and therefore cease being physical.

In this work, we compare the three most general currently available methods for non-equilibrium steady state (NESS) calculations: the straightforward extension of Matrix Product Density Operator (MPDO)dynamical methods to open quantum sys-tem dynamics, and two more recently proposed methods; the DMRG-like variational method, and the local purification method. This last method is the most recently published of all three, and preserves operator positivity by construction. We first provide proof-of-concept calculations in a well-controlled setting of a 6-site dissipa-tive Heisenberg chain. We then extend our calculations to larger-size chains of 32 sites and compare the methods in each regime.

At low maximum bond dimension size D, we find that the variational method per-forms best in terms of accuracy when the spin profile of an edge-dissipative Heisenberg chain is considered. At higher bond dimension, it becomes impossible to distinguish between the local purification method and the MPDO variational method. It must, however, be noted that the variational method does appear to perform best in approxi-mating the complete NESS. The least accurate method is the naive MPDO dynamical method. This method also causes the simulated density operator elements to acquire a large imaginary part, hence spoiling the physicality of the operators. The compu-tational costs of the methods scale as: ∝ O(d5_D2_{) + O(d}5_D3_{) (MPDO-variational),}

O(d2D3D2_w) + O(d4D2D_w3) (MPDO-variational, with D_w the maximum bond dimen-sion of the MPO representation of the Lindlbladian) and O(d4_KD2

) + O(d3K2D3) (local purification, with K the maximum Kraus dimension) . This means that, in

practise, the local purification method is most efficient in terms of accuracy/cost.

“

In der Beschränkung zeigt sich erst der

Meister.

A C K N O W L E D G E M E N T S

First and foremost, I would like to sincerely thank my primary thesis advisor prof.dr. Philippe Corboz of the Institute for Theoretical Physics Amsterdam at the University of Amsterdam. After a bit of a rough time choosing a project, he generously offered me a project as one of his master’s students. During this project, his door was always open, whenever I had questions or wanted to discuss intermediate results. Even when I doubted myself, he showed continuous support and inspired me to acknowledge and develop my own abilities. For all this, I am forever indebted and grateful.

Besides my primary advisor, I wish to express my gratitude to prof.dr. Vladimir Gritsev, also of the Institute for Theoretical Physics Amsterdam at the University of Amsterdam. At first my primary supervisor, he kindly arranged for a different project for me when our interests no longer aligned. As my secondary advisor, he continued to support me and provide me with useful discussions and suggestions.

I would also like to thank the other members of the tensor network group at the UvA. They were always available for discussion and made my masters’ project a truly memorable experience.

Finally, I must express my very profound gratitude to my parents, to my brother and girlfriend for providing me with unfailing support and continuous encouragement throughout my years of study and through the process of researching and writing this thesis. This accomplishment would not have been possible without them.

Karel D. Temmink

C O N T E N T S

1 i n t ro d u c t i o n 8

i t e n s o r n e t wo r k t h e o ry 10

2 i n t ro d u c t i o n 12

2.1 Why Tensor Networks? . . . 12

2.2 Diagrammatic Notation . . . 13

3 m at r i x p ro d u c t s tat e s 15 3.1 Formalism . . . 15

3.2 Singular value decomposition . . . 17

3.3 Decomposition of a state into an MPS . . . 19

3.4 Norms, overlaps and expectation values with MPSs . . . 23

3.5 Compression of MPSs . . . 26

4 m at r i x p ro d u c t o p e r at o r s 30 4.1 Formalism . . . 30

4.2 Decomposition of an Operator into a Matrix Product Operator . . . . 31

4.3 Operations with Matrix Product Operators . . . 34

4.3.1 Application of a MPO onto an MPS . . . 34

4.3.2 Addition and Multiplication . . . 35

4.4 Norms and Expectation Values with MPDOs . . . 37

5 p u r i f i c at i o n o f m p d o s 42 5.1 Introduction . . . 42 5.2 Formalism . . . 42 ii c l o s e d q u a n t u m s y s t e m s 45 6 i n t ro d u c t i o n 46 7 c a l c u l at i o n s w i t h m at r i x p ro d u c t s tat e s 48 7.1 Closed Quantum System Dynamics Simulations (TEBD/tDMRG) . . 48

7.2 Dynamical Ground State Calculations . . . 51

7.3 Variational Ground State calculations . . . 53

8 c a l c u l at i o n s w i t h m at r i x p ro d u c t d e n s i t y o p e r at o r s 60 8.1 Time Evolution (Real and Imaginary) . . . 60

9 c a l c u l at i o n s w i t h i n t h e l o c a l ly p u r i f i e d f o r m 65 9.1 Time Evolution (Real and Imaginary) . . . 65

iii o p e n q u a n t u m s y s t e m s 67

10 i n t ro d u c t i o n 68

11 c a l c u l at i o n s w i t h m at r i x p ro d u c t d e n s i t y o p e r at o r s 73

Contents 6

11.1 Introduction . . . 73

11.2 Dynamical stationary state calculation . . . 73

11.2.1 Example Calculation . . . 75

11.3 Variational stationary state calculation . . . 82

11.3.1 Example Calculation . . . 86

12 c a l c u l at i o n s w i t h i n t h e l o c a l ly p u r i f i e d f o r m 91 12.1 Dynamical Non-Equilibrium Steady State Calculation . . . 91

12.2 Example Calculation . . . 94 12.2.1 Small System . . . 94 12.2.2 Larger System . . . 96 13 c o m pa r i s o n o f d i f f e r e n t a p p roac h e s a n d g e n e r a l d i s -c u s s i o n 98 13.1 Small Systems . . . 98 13.2 Larger Systems . . . 100 13.2.1 Convergence Monitoring . . . 100

13.2.2 Time-Evolution: Flow and Cost . . . 100

13.2.3 Varying Parameters and Summary . . . 102

13.3 Small vs Larger Systems . . . 104

14 c o n c l u s i o n 106 14.1 Summary . . . 106

14.2 Conclusion . . . 107

1

I N T R O D U C T I O N

Theoretical condensed matter physics involves the use of theoretical models to under-stand and mathematically describe properties of quantum states of matter. Many of these models, specially those in a non-equilibrium setting, cannot be solved analyt-ically, and therefore the demand for stable and efficient numerical methods is ever increasing.

It was shown [1,2,3,4,5] that many physically interesting states (such as ground states of local, gapped Hamiltonians, low-energy states, or steady states) are entan-gled only (nearly) locally, with an entanglement entropy that scales as an area law. This means that one does not need to consider the entire Hilbert space of the system in order to accurately simulate a system, even when studying a system’s dynamics [6]. Tensor networks provide a powerful and efficient way of performing such calcu-lations, because several key properties such as an entanglement area law naturally emerge in the tensor network representation. Furthermore, the mathematical prop-erties of such networks allow one to systematically and efficiently control the amount of truncation and accuracy of a calculation.

Whilst a plethora of stable techniques and schemes has been developed for (ground state) calculations pertaining to closed quantum systems throughout the years (see e.g. [7,8,9]), studying the non-equilibrium physics of open quantum systems has re-mained largely out of reach. This is because only very few methods can be straightfor-wardly generalised in a reliable way. Although extensions of closed system dynamical algorithms to open systems are relatively straightforward, the main source of compli-cations remains that truncation of the system dimensions, which is a necessary step in keeping algorithms computationally tractable, tend to render simulated density op-erators nonphysical by destroying their positivity. Variational methods have proven more succesful in overcoming this issue [10,11], but their tremendous computational cost renders them inpractical for true many-body systems. Recently, [12] proposed a new dynamical method that preserves operator positivity by construction, which we will study in this thesis.

The goal of this thesis is to study in detail and gain understanding of the most general algorithms available for open quantum systems steady state calculations. We set out to discover to what extent matrix product operator positivity is broken (e.g.

i n t ro d u c t i o n 8

whether it occurs at all system sizes and/or is dependent on quantities as maximum bond dimension size and simulation time), and how different algorithms compare in accuracy, efficiency and stability.

To this end, we provide the reader with a comprehensive overview of currently available tensor network methods for open quantum systems, including example re-sults for small and intermediately sized systems, and in-depth reviews. Furthermore, we provide for the first time detailed studies and critical reviews of these different methods, and discuss their advantages and disadvantages in a comparative setting.

In parti, we introduce the reader to the concept of tensor networks, diagrammatic notation and several (mathematical) tools that are the basis of all methods discussed in this work. We also summarize numerous distinctive representations of and opera-tions with quantum states and operators.

In partii, we briefly summarize the theory of time evolution in a closed quantum system setting. We demonstrate how previously introduced tools can be implemented to perform ground state calculations. Both dynamical and variational methods are discussed and compared in efficiency and performance. This chapter mainly serves as a means for the reader to acquaint themselves with different tensor network schemes, and as a stepping stone towards the final chapter, where slightly more involved meth-ods are discussed.

Part iiiconcludes this thesis with an overview of the mathematical description of the evolution of open quantum systems. We discuss how the schemes and ansätze from the previous chapters translate to the realm of non-equilibrium physics. We specifically focus on a recently proposed approach to Non-Equilibrium Steady State (NESS) calculations that preserves operator positivity by construction [12], and study for the first time in detail how it compares to already established methods that do not explicitly preserve this quantity [10,11].

PART I

TENSOR NETWORK THEORY

“

Everything starts somewhere, though many

physicists disagree.

But people have

al-ways been dimly aware of the problem with

the start of things. They wonder how the

snowplough driver gets to work, or how the

makers of dictionaries look up the spelling

of words.

2

I N T R O D U C T I O N

2.1

Why Tensor Networks?

Consider, for example, a system of N sites, each described by a local Hilbert space Hi

of local dimension d (e.g. for spin-1/2 systems, d = 2). The total Hilbert space of the system is then simply given by the tensor products of the subspaces: Htot = ⊗iH_i, which has a dimension dN _{that obviously grows exponentially with the number of}

sites. A general pure quantum state (or ensemble of pure states) of any system can be written in terms of some defined basis states:

∣ψ⟩ = ∑

{i}

C_{µ1,µ2,...,µN}∣µ₁⊗µ₂⊗. . . ⊗ µ_N⟩ (1) As the Hilbert space grows exponentially with N, specifying a quantum state by providing the coefficients of the wave function Cµ1,µ2,...,µN in some basis is an

ex-tremely inefficient representation.

However, as it turns out, we need not concern ourselves with the entire Hilbert space. It has been shown [1,2, 3,4,5] that low-energy eigenstates of gapped, local Hamiltonians (i.e. Hamiltonians with nearest or next-to-nearest neighbour interac-tions) obey an area law for their entanglement entropy. This is quite remarkable, as by far most quantum states in the Hilbert space obey a volume law. The above reasoning also holds in reverse: only states satisfying an area law can be low-energy states. However, the manifold consisting of these states is an exponentially small subspace of the Hilbert space.

Our luck does not end here, however. One can also show that, after time-evolving a quantum many-body state for a finite time, one still ends up with a state within an exponentially small manifold of states [6]. This implies that the majority of the Hilbert space is only reachable after exponentially (O (eN_{)) long time evolution. In}

practice, this means that, given some initial quantum many-body state, most of the Hilbert space will be unreachable on time-scales much larger than the age of the uni-verse. Hence, if one picks as initial state a quantum state with area-law behaviour, one only needs to consider an exponentially small manifold of the full Hilbert space. This is precisely where the strength of Tensor Network methods lies. It turns out that Tensor Networks offer an extremely efficient ansatz for representing these

2.2 diagrammatic notation 12

entanglement states. As we will see in the sections to come, we can often reduce the exponentially large number of parameters to a number that scales only weakly polynomially with the number of sites.

2.2

Diagrammatic Notation

We now introduce the diagrammatic notation used in tensor network analyses. Gen-erally, a tensor is represented by a shape (below we chose a sphere). The indices of a tensor are visualised as the ’legs’ of the shape [13]. In this notation, a scalar is simply only a shape and a vector has one leg. Likewise, a matrix has 2 legs, and a rank-n tensor has n legs, as can be seen in figure1.

j

(a) The scalar A

i j (b) The vector ⃗A, with compo-nents Ai i j j (c) The matrix A, with compo-nents Aij i j k l m n

(d) The the order-6 tensor A, with components A_ijklmn Figure 1: Diagrammatic notation of various tensors

With this notation comes a very intuitive representation of tensor operations. Con-nected legs represent indices that are summed over. Open (disconCon-nected) legs thus represent ’free’ indices. Below we list a few examples of common tensor operations in diagrammatic notation: v = ∑ i AiBi A B inner product vj= ∑ i AiBij A B _{vector-matrix} multiplica-tion vik= ∑ j AiBij i A j B k matrix-matrix multiplica-tion

2.2 diagrammatic notation 13 viln= ∑ jkm AijkBjlmCkmn A B C i j l k n

m tensor network

contrac-tion

This is a very nice and intuitive way to think about tensors. For example, cyclic properties of the trace operation are immediately visually apparent:

Tr (ABCDEF ) = A B C D E F (2)

3

M AT R I X P R O D U C T S TAT E S

3.1

Formalism

For a system of N lattices with local dimension d, the most general form of a pure state is given by:

∣ψ⟩ = ∑ µ1,µ2,...,µN C_{µ1,µ2,...,iN}∣µ₁⟩ ⊗ ∣µ₂⟩ ⊗. . . ⊗ ∣µ_N⟩ = ∑ µ1,µ2,...,µN C_{µ1,µ2,...,µN}∣µ₁µ₂. . . µ_N⟩ (3) As mentioned before, this exact representation requires an exponentially large num-ber of coefficients, because Cµ1,µ2,...,µN are d

N _{numbers. Connecting to the previous}

section, all the Cµ1,µ2,...,µN can also be understood to be entries of the tensor C, which

has N indices that can take d values each. Specifying all these elements is, however, still computationally inefficient. A solution to this problem is given by so-called Ma-trix Product States (MPSs).

The main idea of MPSs is that one writes the coefficients of the wave function not as one large tensor, but rather as a network of N smaller tensors (or, equivalently, sets of matrices) M[i], one for each site in the lattice [14]. The contracted indices

have dimension D, often referred to as the bond dimension. For periodic boundary conditions (BCs), M[i]

∈CD×d×D, ∀i. For closed boundary conditions, the same applies to the tensors in the bulk. For the edge tensors, we have

M[1] ∈Cd×D and M[N ] ∈CD×d. In this work, we will frequently adopt the notation

Mµj

ajaj+1 =M

[j] to avoid overly lengthy notation, as we did above. Below is written

explicitly the decomposition of the coefficients Cµ1,µ2,...,µN for different situations and

the diagrammatic tensor network representation of each.

3.1 formalism 15 General C_{µ1,µ2,...,µN} C

. . .

This representation immediately illustrates one major advantage over the ’regular’ representation: Cµ1,µ2,...,µN requires d

N _{terms, while the tensor network}

representa-tions require at most O(NdD2_{)parameters. This is a large gain: instead of}

exponen-tially many coefficients, we can now describe quantum states as a tensor network with a number of parameters that scales only weakly polynomially with N. Contracting the tensor network will then yield dN _{coefficients back.}

This efficient representation does not come for free, however. Representing the co-efficients as contracted tensors introduces an additional degree of freedom, the bond dimension size D. While at first glance there does not seem to be any intuitive physi-cal meaning behind D, some thought should convince the reader that the connecting indices in the network actually represent the structure of the many-body entangle-ment in the state. In a MPS, the value of D is simply the size of the Hilbert space on the left/right side of the corresponding bond. Hence, a larger amount of quantum correlations requires a larger value of D in order to be represented accurately. As we will see shortly, it will always be possible to construct an MPS representation in an exact way, though it comes at a cost: in order to represent a general state as a MPS exactly, D could need to grow exponentially with the system size.

3.2 singular value decomposition 16

3.2

Singular value decomposition

We will start by reviewing an important tool from linear algebra that is not only widely used, but also essential to the success of tensor network representations: the Singular Value Decomposition (SVD). This tool, as we will see shortly, also has a strong connection with the Schmidt decomposition of quantum states in a bipartite system.

An SVD of a matrix M with dimensions dM = (n × m)is given by [9] :

M = U sV† (4)

Here, the matrix U ( dU= (n ×min(n, m)) ) contains as its columns the left singu-lar vectors. The columns are orthonormal, i.e. U†_{U = 1. In the case where n < m,}

also UU†_{=1 holds true.}

Analogously, the matrix V† _{( d}

V = (min(n, m)) × m ) contains the right singular vectors as its rows, which, too, are orthonormal (V†_{V = 1, and if n > m also V V}†_=1).

The matrix s ( ds= (min(n, m) × min(n, m)) ) is diagonal, with non-negative entries

sii=s_i, the so-called singular values of M. In this work, we assume these entries are sorted in descending order, i.e. s1 >s₂>. . . s_r>0. The number of non-zero singular values r is called the (Schmidt) rank of M.

To make the connection with the Schmidt decomposition, we note that we can write any pure quantum state ∣ψ⟩ on a given bipartite system AB in the following way:

∣ψ⟩ = ∑

i,j

C_ij∣i⟩_A∣j⟩_B (5)

where {∣i⟩A}and {∣j⟩B}form orthonormal bases of A, B with dimensions n, m

respec-tively. If we view the coefficients Cij as entries of the matrix C, we can perform an

SVD on it: ∣ψ⟩ = ∑ i,j min(n,m) ∑ k=1 UikskV_kj† ∣i⟩_A∣j⟩_B (6) = min(n,m) ∑ k=1 (∑ i Uik∣i⟩_A)s_k ⎛ ⎝ ∑ j V_kj† ∣j⟩_B ⎞ ⎠ (7) = min(n,m) ∑ k=1 sk∣k⟩_A∣k⟩_B (8) where the sets {∣k⟩A}and {∣k⟩B}again are orthonormal and can be extended to form

3.2 singular value decomposition 17

Now the connection with the Schmidt decomposition is evident. If one restricts the sum to run over the largest r < min(n, m) singular values, the Schmidt decom-position emerges. A value of r = 1 corresponds to pure quantum states, while r > 1 corresponds to entangled states.

In fact, restricting the values the above sum runs over turns out to be essential to efficient tensor network representations of quantum systems. As an example, consider the optimal approximation of a matrix M (rank r) by a matrix ˜M of rank ˜r < r

obtained by performing an SVD, but only keeping the ˜r largest singular values, setting the rest to zero:

˜

M = U˜sV† (9)

˜s =diag(s1, s2, . . . , s˜r, 0, . . . , 0) (10)

This approximation is optimal in the Frobenius norm ∥M∥2

F = ∑i,j∣Mij∣2 and yields a matrix with exactly the same dimensions as M.

This has a natural extension to states, because for a state ∣ψ⟩, the two-norm is identical to the Frobenius norm of the matrix C in equation (5)

∥∣ψ⟩∥ = ∑

i,j

∣C_ij∣2≡ ∥C∥2_F (11)

iff the sets {∣i⟩} and {∣j⟩} are orthonormal, which is the case here. We can therefore approximate C (originally rank r) by a matrix ˜C of rank ˜r as before, and obtain as an optimal approximation to the state ∣ψ⟩:

˜ ∣ψ⟩ = ˜r ∑ k=1 sk∣k⟩_A∣k⟩_B (12)

3.3 decomposition of a state into an mps 18

3.3

Decomposition of a state into an MPS

As we will see in more detail shortly, the maximum size of the bond dimension between the tensors, D, plays a vital role in the success of the tensor network representation. It is possible to represent any wavefunction as an MPS in an exact way, provided that D is large enough. This, however, means that D could have to grow exponen-tially, depending on the state under consideration, and we will not have gained any advantage by using this representation. It does not end here, however: the mathe-matical tools outlined in the above can also be applied to many-body wavefunctions. To demonstrate this, we introduce here mathematical tools that are of much use, following Schollwoecks work [9,15].

Starting from a general state, given by equation (3), we reshape the coefficient tensor Ci1,i2,...,iN (d × d × . . . × d = d

N _{coefficients) into a (d × d}N −1_{) matrix Ψ with}

entries Ψi1,i2...iN =Ci1,i2,...,iN. Performing an SVD on Ψ gives:

Ci1,i2,...,iN = r1 ∑ a1 Ui1,a1sa1,a1(V † ) a1,i2...iN (13) = r1 ∑ a1 Ai1 a1Ψ(a1i2),i3...iN (14)

where, in going from the first equality to the second, the matrix U was reshaped into a (d × r1) matrix A with entries Aia11 =Ui1,a1. Furthermore, we carried out the

multiplication sV† _{and reshaped the result into a new (r}

1⋅d × dN −2) matrix Ψ with entries Ψ(a1i2),i3...iN = sa1,a1(V

†

)

a1,i2...iN. Upon performing another SVD of the Ψ

matrix, we obtain: Ci1,i2,...,iN = r1 ∑ a1 Ai1 a1 r2 ∑ a2 Ua1i2,a2sa2,a2(V † ) a2,i3...iN (15) = r1,r2 ∑ a1,a2 Ai1 a1A i2 a1,a2Ψ(a2i3),i4...iN (16)

After another N − 3 SVDs, the final expression is:

Ci1,i2,...,iN = r1,...,rN ∑ a1,...,aN Ai1 a1A i2 a1,a2. . . A iN aN −1 (17)

Comparing this to the OBC TN from section3.1, we see that the two expressions are identical. The MPS obtained this way by sweeping from the left to the right (i.e. starting on site 1 and working your way to site N) is also known as a left canonical

3.3 decomposition of a state into an mps 19

MPS. The name comes from the fact that all tensors on sites left of site N inherit the orthogonality properties (i.e. are left normalised)from the SVD matrices:

U†U = 1 → δak,˜ak = ∑ ak−1µk (U†) ak,(ak−1µk)U(ak−1µk),˜ak (18) = ∑ ak−1µk (Aµk†) ak,ak−1A µk ak−1,˜ak (19) = ∑ µk (Aµk†Aµk) ak,˜ak (20)

From which follows that ∑µkA

µk†_Aµk _{= 1. Of course, this property is violated on}

the last site, as AµN was not constructed from a U matrix, but rather from the

product sV , which is generally not normalized in any way. It is instructive to ex-amine the A matrices in a bit more detail. If the decomposition is exact (i.e. all non-zero singular values are kept), we see that we maximally have bond dimensions (1 × d), (d × d2), . . . , (dN /2−1×dN /2), (dN /2×dN /2−1), . . . , (d2×d), (d × 1), which indeed grow exponentially with N, as foreshadowed in section3.1.

Cµ1,µ2,...,µN

↓

↓

⋮

Figure 2: Decomposing a quantum state into a left canonical MPS represen-tation

3.3 decomposition of a state into an mps 20

Similar to the above decomposition, it is also possible to start from the right (site

N), and work one’s way to the left (site 1), performing SVDs to split the C tensor: Cµ1,µ2,...,µN =Ψµ1µ2...µN −1,µN (21) = ∑ bN −1 Uµ1...µN −1,bN −1sbN −1,bN −1(V † ) bN −1,µN (22) = ∑ bN −1 Ψµ1...µN −2,µN −1bN −1B µN bN −1 (23) = ∑ bN −2,bN −1 Uµ1...µN −2,bN −2sbN −2,bN −2(V † ) bN −2,µN −1bN −1B µN bN −1 (24) = ∑ bN −2,bN −1 Ψµ1...µN −3,µN −2,bN −2B µN −1 bN −2,bN −1B µN bN −1 (25) =. . . (26) = ∑ b Bµ1 b1 . . . B µN −1 bN −2,bN −1B µN bN −1 (27) Cµ1,µ2,...,µN

↓

↓

⋮

Figure 3: Decomposing a quantum state into a right canonical MPS repre-sentation

The tensors in this type of MPS (except the one on the first site) are right

normal-ized (∑_µ

kB

3.3 decomposition of a state into an mps 21

right canonical.

The final possibility is a combination of both aforementioned decompositions: the

mixed canonical MPS decomposition. To construct from a wave function a MPS that

is mixed canonical with respect to a site j, one makes sites 1 through j − 1 left nor-malised and sites N through j + 1 right nornor-malised. The final ’leftover’ matrices from the canonization are absorbed with the tensor on site j. The formalism is identical to what has already been outlined, so only a diagrammatic representation will be given below. This particular form of a MPS has many advantages, specially when it comes to calculating norms and expectation values, as we will see later.

Note that it is always possible to bring any MPS into a (mixed) canonical form, simply by reshaping the on-site tensors and performing the required amount of SVDs on them, absorbing the U, s, V†_{into the tensors to the left or right, as needed. This is,}

however, not the most efficient or desirable approach, as it requires one to explicitly construct/store the full state in vector form first. It is a common practise to initialise states directly in the MPS representation, which can be done very cheaply (i.e. low

3.4 norms, overlaps and expectation values with mpss 22

3.4

Norms, overlaps and expectation values with MPSs

With the formalism of MPSs in hand, we now review some of the most important tools needed to extract information from these states: Norms and expectation values in tensor network representation. It is most instructive to start with the general overlap between two states. Given two states ∣ψ⟩ = ∑µ,a(∏iMaµiiai+1) ∣µ⟩ and ∣φ⟩ =

∑_µ,a(∏_iNµ_˜aii_˜ai+1) ∣µ⟩, the overlap ⟨φ∣ψ⟩ is given by:

⟨φ∣ψ⟩ = ∑ µ

NµN†⋯Nµ2†Nµ1†Mµ1Mµ2⋯MµN (28) This calculation involves summation over dN _{sets of matrix multiplications, and the}

beauty of MPS seems to be lost. However, we can rewrite the sums as ⟨φ∣ψ⟩ = ∑ µN NµN†⎛ ⎝ ⋯ ⎛ ⎝ ∑ µ2 Nµ2†⎛ ⎝ ∑ µ1 Nµ1†_Mµ1⎞ ⎠ Mµ2⎞ ⎠ ⋯ ⎞ ⎠ MµN (29)

and carry out the sums from the innermost bracket outwards. As can be seen in fig4a, the computational cost can be greatly reduced by contracting only two tensors at once. The optimal contraction sequence thus reads:

1. Mµ1 and Nµ1†are contracted by summing over the physical index µ

1, resulting

in a matrix C with indices a1, ˜a1. This operation is of order O(dD2). 2. The matrix C is then contracted with Mµ2 by summation over a

1 (resulting in

a new tensor C with indices ˜a1, a2, µ2), which scales as O(dD3). 3. Finally, C is contracted with Nµ2†, by summing over ˜a

1, µ2. This gives a new

matrix C with indices a2, ˜a2. This operation is of cost O(d2D2) and essentially leaves us with the same set-up as in step 2, only one site to the right.

4. Steps 2 and 3 are repeated until the final site is reached. The computational cost does not increase any more from this point.

The entire process is illustrated in tensor network notation in figure (4a). Because generally D ≫ d, the maximum cost of the total calculation is of order ∼ O(NdD2

), which is an incredible reduction from exponential complexity to weakly polynomial.

The norm of a MPS can then be calculated in a similar manner, by considering the special case ∣φ⟩ = ∣ψ⟩. But we can do more here: if the state ∣ψ⟩ is in (mixed) canonical form with respect to one of the sites j, every two-tensor contraction of the form Mµi†_Mµi with i ≠ j simply results in the δ line from equation (20). All

contractions of tensors on sites 1 to j thus simply result in a δ line connecting the left auxiliary bonds of Mµj and Mµj†. Of course, the same reasoning holds for tensors

3.4 norms, overlaps and expectation values with mpss 23

to the right of site j: because these are right normalised, contracting them simply gives another δ line, this time connecting the right bonds of Mµj and Mµj†. One is

thus left with the contraction:

⟨ψ∣ψ⟩ = ∑

µj

Tr (Mµj†_Mµj₎ (30)

which is diagrammatically illustrated in figure (4b) below and has computational cost of order O(dD2

), which is a factor ∼ N less than the general calculation.

→

→

(a) Full, most general calculation of the overlap between two MPSs ∣ψ⟩ = ∑_µ,˜a˜(∏iM˜

µi

˜ai˜ai+1) ∣µ⟩˜ and ∣φ⟩ = ∑µ,a(∏iN µi

˜ai˜ai+1) ∣µ⟩. The figure illustrates the most efficient contraction sequence. As a first step Mµ1 and Nµ1† are contracted by summing over the physical

index µ1. In step two, the resulting matrix is then contracted with

Mµ2. The tensor that results from this contraction is contracted

with Nµ2†. These steps are repeated for all remaining sites in the

chain. The computational cost scales weakly polynomially with system size, with maximum cost being of order O(NdD3

).

→

→

→

. . .

→

→

. . .

→

(b) Calculation of the norm of a MPS if the state is in mixed canonical form with respect to site j

Figure 4: Two different approaches to calculating the norm of a MPS. In (4a) , the full and most general overlap between two MPSs ∣ψ⟩ , ∣φ⟩ is shown: the overlap is calculated by contracting the full tensor network that makes up the MPS. The optimal contraction sequence is highlighted at every step. The norm of a MPS can be calculated similarly, by setting ∣ψ⟩ = ∣φ⟩. (4b) shows the tremendous simplification that can be made when the MPS is in mixed canonical form.

Another important quantity of interest is the expectation value of an operator ˆ

3.4 norms, overlaps and expectation values with mpss 24

review) and applied onto MPSs. Given an operator ˆO = ∑µ,µ,b˜ (∏iW µiµ˜i

bibi+1) ∣µ⟩ ⟨µ∣˜

the most general matrix elements of ⟨φ∣ ˆO ∣ψ⟩ are simply: ∑ µN,µ˜N OµN,˜µN_NµN†⎛ ⎝ ⋯ ⎛ ⎝ ∑ µ2 Oµ2,µ˜2_Nµ2†⎛ ⎝ ∑ µ1 Oµ1,µ˜1_Nµ1†_Mµ˜1⎞ ⎠ Mµ˜2⎞ ⎠ ⋯ ⎞ ⎠ Mµ˜N (31)

This calculation is completely analogous to the calculation of the overlap, with the difference that the sums now run over two sets of physical indices and the bond indices of the MPO representation of ˆO. We can, once more, come up with the optimal contraction sequence to reduce the computational cost. Typically, one starts by constructing a tensor C[1]_{, which consists of the contraction of N}[1]† _{with ˆO}[1]

and M[1] _{over the physical indices, leaving a tensor with three free indices. One then}

calculates ⟨φ∣ ˆO ∣ψ⟩ iteratively, by constructing new Cs from the previous ones:

C_a[k]_k,˜ak,bk = ∑ ak−1,µk N_aµ_k−1,akk† ⎛ ⎝ ∑ ˜ µk,bk−1 Oµkµ˜k bk−1,bk ⎛ ⎝ ∑ ˜ak−1 C_a[k−1]_{k−1,˜ak−1,bk−1}Mµ˜k ˜ak−1,˜ak ⎞ ⎠ ⎞ ⎠ (32)

where i runs up to and including N, and C[N ] is once again a scalar, containing the

result. Equation (32) is already written in a form that highlights the optimal contrac-tion sequence, working one’s way from the innermost bracket outwards. For clarity, this contraction is also shown diagrammatically, in figure (5a). The expectation value is then simply the case where ∣φ⟩ = ∣φ⟩. If we take the bond dimension size of the MPO to be Do, the cost of calculating the new C[k] is only 2 × O(dDoD3) + O(d2D_o2D2) with this method, instead of a naive O(d2_D2

oD4).

Of course, not all operators are most conveniently represented as an MPO, such as an operator consisting of a collection of local operators, acting on one site only each ˆO = ⊗i ˆo[i]. In this case, equation (32) still holds and represents the optimal

contraction sequence, with the difference that one must omit the sums over the bond indices {b}, because the operators do not have bond indices in this case. Now, the cost is reduced even further from a naive O(d2_D4

), down to 2 × O(dD3) + O(d2D2). Typically, the operators of interest (e.g. correlators, spin operators, etc) are either (few site) local, or can be written as a sum of local terms. In this case, the double sum reduces to a single sum on several sites, as only the identity acts there. This again allows for significant improvements on the calculation of expectation values by making maximal use of the canonical form of the MPS. If an operator acts only on a certain site k, contracting tensors left and right of k will again result in δ lines, in the same way as before. In the end, one only needs to consider the MPS tensors on site k and the calculation simplifies further to :

⟨ψ∣ ˆO ∣ψ⟩ = ∑

µkµ˜k

Oµk,µ˜kTr (Mµk†_Mµ˜k₎ (33)

3.5 compression of mpss 25

→

→

(a) Full, most general calculation of the matrix elements of the ex-pectation value of an operator ˆO = ∑µ,µ,b˜ (∏iW

µiµ˜i

bibi+1) ∣µ⟩ ⟨µ∣˜ between two MPSs ∣ψ⟩ = ∑µ,a(∏iMµaiiai+1) ∣µ⟩ and ∣φ⟩ = ∑µ,a(∏iN

µi

˜ai˜ai+1) ∣µ⟩˜ . The figure illustrates the most efficient con-traction sequence. As a first step Mµ1 and Wµ1µ˜1 are contracted

by summing over the physical index µ1. Thereafter, the resulting tensor is contracted with Nµ1†. The tensor that results from this

contraction is contracted with Mµ2, after which the process is

re-peated, as in equation (32). The computational cost scales weakly polynomially with system size, with maximum cost being of or-der O(N [2dD3_D

w+d2D2D2_w]), where Dw denotes the maximum

bond dimension size in ˆO.

→

→ . . . →

(b) Calculation of the expectation value (i.e. ∣ψ⟩ = ∣φ⟩) of ˆO from a) for a local operator acting on site j, when the MPS is in mixed canonical form w.r.t. that site.

Figure 5: Two different approaches to calculating matrix elements of an operator ˆO. In (5a) , the full and most general calculation using two MPSs ∣ψ⟩ , ∣φ⟩ is shown. At every step, we have highlighted the computationally opti-mal contraction sequence. The expectation value of an MPDO can be calculated similarly, by setting ∣ψ⟩ = ∣φ⟩. (5b) shows the tremendous sim-plification that can be made when the MPS is in mixed canonical form.

3.5

Compression of MPSs

In tensor network analyses, it is commonplace that the bond dimensions in a MPS become too large for efficient numerical calculations. It is thus worthwhile to review two important methods for approximating a MPS ∣ψ⟩ (built from {Mµi_} matrices)

with maximum bond dimension size D by a MPS ∣ ˜ψ⟩ (built from {Nµi_} matrices)

with maximum bond dimension size D′

<D.

The simplest and computationally cheapest way of accomplishing such a compres-sion is through SVDs. Much like the decomposition of the wave function coefficient

3.5 compression of mpss 26

tensor into a MPS, one can contract the bond between two adjacent MPS tensors and perform a SVD on the 2-site tensor, keeping at most D′ singular values. Depending

on the kind of canonical form one wants for the new MPS, one absorbs the singu-lar values in the left or right tensor. This is then repeated for all bonds in the network.

N[i]

→

⃗

N[i]

Figure 6: In variational methods, every N[i] _{tensor is reshaped into the} ’fat-legged’ vector ⃗N[i] and optimized iteratively.

Whilst computationally cheap, this method is neither the most stable, nor the optimal compression technique. A better, but more costly, alternative is given by variational compression. Mathematically, we want to minimize ∥ ∣ψ⟩ − ∣ ˜ψ⟩ ∥2

F with

respect to ∣ ˜ψ⟩ locally. Hence, at each site i, the equation of interest is:

∂ ∂N[i],† ( ⟨ψ∣ψ⟩ − ⟨ ˜ψ∣ψ⟩ − ⟨ψ∣ ˜ψ⟩ + ⟨ ˜ψ∣ ˜ψ⟩ ) =0 (34) = ∂ ∂N[i],†( ⟨ ˜ψ∣ ˜ψ⟩ − ⟨ ˜ψ∣ψ⟩ ) (35) Ð→ ∂ ∂ ⃗N[i],† ( ⃗N[i],†N ⃗N[i]− ⃗N[i],†⋅ ⃗E ) =0 (36) Where ⃗N[i] is just the tensor N[i] but with all its coefficients arranged as a vector,

as illustrated in fig (6). Likewise, N and ⃗E are the so-called environments of the tensors ⃗N[i],†, N[i] (i.e. the TNS for ⟨ ˜ψ∣ ˜ψ⟩ and ⟨ ˜ψ∣ψ⟩ without ⃗_N[i],†, N[i]) reshaped

to matrices, as illustrated in fig (7) below:

The optimization is then done site-by-site: One sweeps left-to-right and then back right-to-left, updating only one site at a time. From equation (36), we see that we need to solve at each site system of linear equations:

N ⃗N[i]= ⃗E (37)

Which gives a new vector ⃗N[i], which is shaped back into the tensor N[i] and used

as the new MPS entry.

We can, however, improve upon this. We do not need to construct the environ-ments N , E completely from scratch for each site i, but we can recycle information

3.5 compression of mpss 27

−

Figure 7: The energy minimization expression in the parantheses in (36) in diagrammatic notation. Taking the derivative with respect to the tensor N[i],†_{then simply corresponds to removing that tensor} from the network. The parts of the networks outlined incyanand

orangeare reshaped into the matricesN and E, respectively.

from all sites that are not currently being updated: we construct two set of ’blocks’, ˜

S, S, for the left and right term in fig (7), respectively, which are stored in an array.

In this case, the element S[N ] ( ˜S[N ]) is constructed by contracting M[N ] with

N[N ],† (N[N ] with N[N ],†). One then constructs all other elements N ↦ 2 iteratively,

by contracting the previous elements in S, ˜S with the next elements of the environ-ments. Of course, this is completely analogous to the optimal calculation of MPs overlaps (see section3.4for more information). We can thus simply find the optimal contraction sequence from figure (4a):

S_a[k]_k,˜ak = ∑ ˜ak+1,µk Mµ_˜ak,˜ak+1k† ⎛ ⎝ ∑ ak+1 S_a[k+1]_k+1,˜ak+1N_aµ_k,ak+1k ⎞ ⎠ ˜ S_a[k]_k,˜ak = ∑ ˜ak+1,µk Nµ_˜ak,˜ak+1k† ⎛ ⎝ ∑ ak+1 ˜ S_a[k+1]_k+1,˜ak+1N_aµ_k,ak+1k ⎞ ⎠ (38)

Of course, S[1]_{, ˜S}[1]_{need not be constructed for the optimization of the first tensor,}

as the environments simply consist of the entire right part of the chain. After the first tensor is optimized, S[1] _{( ˜S}[1]_{) is constructed by contracting M}[1] _{with N}[1],†

(N[1] _{with N}[1],†_{). As the algorithm sweeps from left to right, after optimizing each}

N[k], the kth elements of S, ˜S are now calculated as:

S_a[k] k+1,˜ak+1 = ∑ ˜ak,µk Mµ_˜ak,˜ak+1k† ⎛ ⎝ ∑ ak+1 S_a[k+1] k+1,˜ak+1N µk ak,ak+1 ⎞ ⎠ ˜ S_a[k] k+1,˜ak+1= ∑ ˜ak,µk Nµ_˜ak,˜ak+1k† ⎛ ⎝ ∑ ak+1 S_a[k+1] k+1,˜ak+1N µk ak,ak+1 ⎞ ⎠ (39)

Note that the indices of the sums and S have shifted with respect to equation (107), and that all open legs are now on the right of the blocks. At any intermediate

3.5 compression of mpss 28 N = ˜S[k−1]⊗ ˜S[k+1] (40) ⃗ E = ∑ ak,ak+1 S[k−1]N[k]S[k+1] (41)

Once a new M[k]has been found, one updates the variational MPS and S, ˜S using

equation (39) (if sweeping to the right) or equation (38) (if sweeping to the left). Another simplification can be made: if one ensures that ∣ ˜ψ⟩ is kept in mixed canonical form with respect to site i throughout the compression sweeps, calculating N is no longer necessary. All contractions between N[k], N[k],† simply result in

δ functions (lines in diagrammatic notation) between the auxiliary bonds. Hence,

N =1, reducing the problem of solving (37) to a trivial:

⃗

N[i]

= ⃗E (42)

Lastly, one can greatly speed up the variational algorithm by providing a proper initial ansatz for the compressed state. Generally, a good initial MPS is given by a SVD compression of the original state.

4

M AT R I X P R O D U C T O P E R AT O R S

4.1

Formalism

Tensor Networks can not only be used to represent pure quantum states, but also mixed states (in the form of a density operator ˆρ) and operators in general. In complete analogy to our consideration of pure states, we note that any operator can be written as: ˆ O = ∑ µ,µ˜ Cµ1,...,µN ˜ µ1,...,µ˜N ∣µ⟩ ⟨µ∣˜ (43)

Where the µ, ˜µ form bases of the outgoing and incoming states of the operator, respectively. Each µi, ˜µi runs over d values, just as was the case in the MPS

rep-resentation. It turns out to be a completely natural generalization to think about coefficients ⟨µ∣ ˆO ∣ ˜µ⟩ of an operator ˆO as a product of sets of matrices, such that:

ˆ

O = ∑

µ,µ˜

Wµ1µ˜1

⋅Wµ2µ˜2⋅. . . ⋅ WµNµ˜N∣µ⟩ ⟨µ∣˜ (44) just as we did for pure states. An operator represented by a multiplication of sets of matrices (or trace over the multiplication in the case of periodic boundary condi-tions) is called a Matrix Product Operator (MPO). If the operator is not just any operator, but the density operator of a quantum system, represented as a MPO, it will be referred to explicitly as a Matrix Product Density Operator (MPDO).

The diagrammatic representation is then as expected: instead of one physical leg on each site, we now have two, one for ingoing states, and one for outgoing states:

General C_µµ_˜1,...,µN 1,...,µ˜N C µ1 . . . µN ˜µ1 . . . ˜µN Periodic BCs a1,a∑2,...,aN Wµ1,µ˜1 a1,a2 W µ2,µ˜2 a2a3 . . . W µN,µ˜N aNa1 =Tr (W_µ[1] 1,µ˜1W [2] µ2,˜µ2. . . W [N ] µN,µ˜N) 29

4.2 decomposition of an operator into a matrix product operator 30 Open BCs ∑ a1,a2,...,aN −1 Wµ1,µ˜1 a1 W µ2,µ˜2 a1a2 . . . W µN,µ˜N aN −1 =W_µ[1] 1,µ˜1W [2] µ2,µ˜2. . . W [N ] µN,µ˜N

Again the same apparent advantage emerges: we were again able to represent exponentially many coefficients by polynomially many in the MPO representation. Once again, we find that this implies that D must grow exponentially.

In the following sections several operations with MPOs are reviewed. Most of these will be completely analogous to what was discussed in chapter3, and we will therefore only briefly review the necessary theory and tools.

4.2

Decomposition of an Operator into a Matrix

Prod-uct Operator

Because any operator can be written as in equation (43), we can always represent it as a tensor network. The big C tensor can be broken down into the individual one-site tensors through successive SVDs. All decomposition methods (left, right and mixed canonical) generalize in a straightforward manner: instead of the single physical index µi, we now consider the two indices µi, ˜µi as one big index (µi˜µi). This then allows for immediate truncation of the MPO, by keeping a maximum of D singular values at each bond.

As mentioned before, it is always possible, but in practise almost never efficient to construct a MPO this way. This is because one would need to construct the full ma-trix form of ˆO first, which has exponentially many elements in N (for Hamiltonians, this is of course dN

×dN elements).

Often, there is a more clever way that does not require a large bond dimension size D. For local operators, such as ˆSz

j, acting only on site j, this is trivial: all

tensors but the one corresponding to that site must be (reshaped) identities, with bond dimension D = 1 between them (of course, one does not even need to explicitly construct the identities in many cases). This can be easily generalised to operators acting only on a few sites. Global operators, such as Hamiltonians, must be treated more carefully. For most systems, one can write the total Hamiltonian as a sum of more local (such as two-local) terms: H = ∑<i,j>hi,j, which can be easily translated

to a set of two-local MPOs, which one can then apply successively onto an MPS or other MPO (See section7.1for more details and an example).

4.2 decomposition of an operator into a matrix product operator 31 Cµ˜1,µ˜2,...,˜µN µ1,µ2,...,µN {µ} {µ}˜

↓

↓

⋮

Figure 8: Decomposing an operator into a MPO by performing successive SVDs. The process itself is a trivial generalisation from an MPS decomposition: the single physical index appearing at each site is replaced by a pair of indices.

Often, one can even write a Hamiltonian as a MPO in an exact way, because of the structure that is present in the terms. As a rather general example, let us consider the isotropic anti-ferromagnetic spin-1/2 Heisenberg chain:

H = N ∑ j=1 S_jzS_j+1z +1 2(S + jS − j+1+S_j−S_j+1+ ) (45) =S₁z⊗S₂z⊗1₃⊗1₄⊗1₅. . . +1₁⊗S₂z⊗S₃z⊗1₄⊗1₅. . . +. . . (46)

It turns out we only need three building blocks (these will be operator-valued matri-ces) to construct its MPO form entirely.

4.2 decomposition of an operator into a matrix product operator 32

The terms in (46) follow a clear pattern, which allows us to define 5 ’states’ of operator strings. At some site in the chain, we can have either one of the following:

1. only units to the right

2. one S+on the right-adjacent site

3. one S−on the right-adjacent site

4. one Sz _{on the right-adjacent site}

5. A completed interaction term to the right

We can then simply list all transitions allowed to form Hamiltonian terms: • 1 ↦ 1, by 1 • 1 ↦ 2, by S+ • 1 ↦ 3, by S− • 1 ↦ 4, by Sz • 2 ↦ 5, by 1 2S − • 3 ↦ 5, by 1 2S + • 4 ↦ 5, by Sz • 5 ↦ 5, by 1

These transitions then define the following building blocks (all other elements are simply the zero operator ):

h1= (_{O 1/2S}− 1/2S+ _Sz 1) , (47) hi= ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ 1 O O O O S+ O O O O S− O O O O Sz O O O O O 1/2S− 1/2S+ _Sz 1 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ i ≠1, N (48) hN = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ 1 S+ S− Sz O ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ (49)

such that we have ˆH = ∏_khk. We have thus found an exact MPO representation

of the Hamiltonian, with bond dimension D = 5. Next-nearest neighbour terms and additional (transverse) field terms can be added trivially by allowing for more states and transitions.

4.3 operations with matrix product operators 33

4.3

Operations with Matrix Product Operators

4.3.1 Application of a MPO onto an MPS

Applying a MPO ( ˆO = ∑ µ,µ˜

(∏_iWµiµ˜i) ∣µ⟩ ⟨µ∣˜ ) onto an MPS (∣ψ⟩ = ∑_µ˜(∏_iMµ˜i) ∣µ⟩˜ ) is done exactly as expected: one simply sums over the indices of the corresponding physical legs to obtain a new MPS:

ˆ O ∣ψ⟩ = ∑ µ,µ˜ (∏ i Wµiµ˜i ) (∏ i Mµ˜i ) ∣µ⟩ (50) = ∑ µ,µ˜ (∏ i Wµiµ˜i_Mµ˜i ) ∣µ⟩ (51) = ∑ µ,µ,˜ a,b (∏ i Wµiµ˜i aiai+1M ˜ µi bibi+1) ∣µ⟩ (52) = ∑ µ,ξ (∏ i Nµi ξiξi+1) ∣µ⟩ (53) = ∑ µ (∏ i Nµi ) ∣µ⟩ = ∑ µ Nµ1_Nµ2_{. . . N}µN_∣µ⟩ (54)

In (52) we reshaped the W and M matrices on the first and last sites such that they have another (dummy) index, a1, b1 for the tensors on site 1, and aN +1, bN +1 for

those on site N, this is purely for notational convenience, as these dummy indices can only take one value: 1. In going from (52) to (53), we carried out the sum-mation over ˜µ, and reshaped the resulting matrices such that they have again only two indices, ξi, ξi+1, with ξi = (a_ib_i), i.e. N_ξµi

iξi+1 =N µi (aibi)(ai+1bi+1)= ∑µ˜W µiµ˜i aiai+1M ˜ µi bibi+1.

This means that although once again tensor networks were able to reduce an expo-nentially complex problem to a polynomially difficult one (i.e. of maximum order O(N d2D2_MPSD_MPO2 )), it comes at the cost of enlarged bond dimensions.

Diagrammatically, the beauty of the simplicity in this operation is restored: one connects the MPS tensors with the corresponding MPO tensors via the incoming physical legs, and contracts over the connected lines, as shown in figure9.

Of course, there exist operators that act more locally then a Hamiltonian, as dis-cussed in section (4.2). For these MPOs, where most tensors will simply consist of the identity with bond dimension D = 1, the application onto an MPS simplifies: Because the sites with just the identity as tensor leave the corresponding sites in the MPS invariant (even the bond dimension does not change size), one can simply ignore them, and only contract the sites where the MPO acts non-trivially.

4.3 operations with matrix product operators 34

ˆ

O

∣

ψ⟩

(a) Full application

DM P S DM P O contraction

ÐÐÐÐÐÐÐ→

Ð

ÐÐÐÐ

→

∣

ψ

′

⟩

(b) Result of the application: a MPS with enlarged bond dimensions Figure 9: a): Application of a MPO onto an MPS (∣ψ′_{⟩ = ˆO ∣ψ⟩) in}

diagram-matic notation. b): The result is again an MPS with the same physical dimensions, but with enlarged bond dimensions that are the product of the bond dimensions of the MPO and the MPS:

D′=DM P O⋅DM P S.

4.3.2 Addition and Multiplication

Adding two MPOs ˆO = and ˆP = with bond dimensions DO and DP respectively can

simply be done by taking the direct sum of the two: the new tensor on every site i will then be Wµiµ˜i_{⊕ Y}µiµ˜i, which has bond dimension D = D

O+D_P.

Equally trivial is the multiplication of two MPOs, ˆO ˆP . Contracting over one set of physcial legs of each MPO, we find that we can simply multiply the MPO dimensions on the tensor level, and construct the product in a straightforward manner. If we write the product ˆO ˆP as a MPO consisting of sets of Zµiµ˜i matrices, the elements of

these are clearly:

Zµiµ˜i

(ai−1bi−1),(aibi)= ∑

µ′ Yµiµ ′ i ai−1,aiW µ′ iµ˜i bi−1,bi (55)

Which is diagrammatically illustrated in fig 10 below to illustrate the similarity with the MPO-MPS product reviewed before.

4.3 operations with matrix product operators 35

ˆ

O

ˆ

P

(a) Full application

DP DP contraction

ÐÐÐÐÐÐÐ→

Ð

ÐÐÐÐ

→

ˆ

O ˆ

P

(b) Result of the application: a MPO with enlarged bond dimensions Figure 10: a): Multiplication of two MPOs ˆO, ˆP in diagrammatic notation.

b): The result is again an MPO with the same physical dimen-sions, but with enlarged bond dimensions that are the product of the bond dimensions of the MPOs: D′_=D

O⋅DP.

We have seen that the application of MPOs generally leads to an increase in bond dimension size. If the result of a MPO application is an MPS, one can of course compress this using the methods outlined in section3.5. If this is not an option, one can turn to compressing the MPO itself before application, which is done in the exact same way as MPSs are compressed, both by SVDs and iteratively, again with the double index reshaped into a single large index acting as a big MPs index.

It is important to mention that there has also been proposed a manner to compress the MPO-MPS product by SVDs: [16]. The idea is to bring both the MPO and the MPS into the same canonical form, to ensure that the block states are almost orthonormal. One then compresses the product, which is carried out tensor-by-tensor or ’on-the-fly’ for computational efficiency, by successive SVDs. Even if the result is not of great quality, it often forms a great initial ansatz for variational compression methods.

4.4 norms and expectation values with mpdos 36

4.4

Norms and Expectation Values with MPDOs

Because we are now working with MPDOs, and not with MPSs, calculating expec-tation values and norms is done in a slightly different manner than before. Where before we needed multiplications involving ⟨ψ∣ and ∣ψ⟩, now calculations involve trac-ing over (part of) the density operator ˆρ.

The trace norm Tr( ˆρ) is thus a quantity of central interest. Let ˆρ be a density operator with operator basis ∣µ⟩ ⟨ ˜µ∣ as in equation (43). The mathematical descrip-tion of the trace Tr( ˆρ) = ∑µ,µ˜(C

µ1,⋯,µN

˜

µ1,⋯,µ˜N ∏kδ

µk ˜

µk) has a natural and intuitive tensor

network representation:

Tr( ˆρ) =

. . .

. . .

We see that taking the trace simply corresponds to connecting both sets of physical legs and contracting over all connected legs, which is typically done tensor-by-tensor. For physical states, we have the obvious requirement that Tr( ˆρ) = 1. This is, how-ever, not a strict condition in numerical analyses, because any MPDO can always be normalised by dividing it by its norm.

Expectation values of operators are then given by ⟨ ˆO⟩ = Tr( ˆρ ˆO), or Tr( ˆρ ˆO)/ Tr( ˆρ) if the density operator was not properly normalised. If ˆO = ∑

µ,µ′ (∏_iWµiµ ′ i_{) ∣}µ⟩ ⟨µ′∣and ˆ O = ∑ µ,µ˜

(∏_iWµiµ˜i) ∣µ⟩ ⟨µ∣˜ , the corresponding tensor network expression is:

Tr( ˆρ ˆO) = ∑ µ,µ,µ˜ ′ (∏ k Wµkµ′kYµ ′ kµ˜k_δµk ˜ µk) (57) =

. . .

. . .

ˆ

O

ˆ

ρ

Which is again best executed on a site-by-site basis, sweeping from the left to the right.

4.4 norms and expectation values with mpdos 37

It is, however, not uncommon that operators can be written as sums of more local terms e.g. Hamiltonians: ˆH = ∑⟨i,j⟩ˆhij or total-spin operators: ˆS

z

tot = ∑kSˆkz. Below

we provide an example for a six-sites chain:

E =Tr( ˆρ ˆH)/Tr( ˆρ)

=

(

+

+

+

+

)/

(59) This can be straightforwardly extended to the completely local case. In these cases, we can devise a clever contraction scheme such that as much information as possible is recycled. The idea is very similar to the scheme described in section 3.4. One constructs two ’blocks’, C and N , which are iteratively grown by contraction with tensors on the N sites. We initialize with:

C[1]_a 1a3 = ∑ µ1,µ˜1,µ′1 µ2,µ˜2,µ′2 hµ˜1µ˜2 µ′ 1µ′2(∑_a 2 Wµ1µ˜1 a1a2W µ2µ˜2 a2a3)δ µ1 µ′ 1δ µ2 µ′ 2 (60) = a2 a1 a3 ˜µ1 ˜µ2 µ1 µ′ 1 µ2 µ′ 2 (61) =

a

a

4.4 norms and expectation values with mpdos 38 And for N[1]_: N[1]_a 1a2 = ∑ µ1µ˜1 Wµ1µ˜1 a1a2δ µ1 ˜ µ1 (63) = a₁ a₂ µ1 ˜µ1 (64) =

a

a

All successive blocks are then constructed as follows:

C[_ak] 1ak+2= ∑ ak+1 C[_ak−1] 1ak+1 ⎛ ⎝ ∑ µk+1,˜µk+1 Wµk+1µ˜k+1 ak+1ak+2δ µk+1 ˜ µk+1 ⎞ ⎠ + ∑ ak N[_ak−1] 1ak ⎛ ⎜ ⎜ ⎜ ⎜ ⎝ ∑ µk,µ˜k,µ′ k µk+1,µ˜k+1,µ′ k+1 hµ˜kµ˜k+1 µ′ kµ ′ k+1 ⎛ ⎝ ∑ ak+1 Wµkµ˜k akak+1W µk+1µ˜k+1 ak+1ak+2 ⎞ ⎠ δµk µ′ k δµk+1 µ′ k+1 ⎞ ⎟ ⎟ ⎟ ⎟ ⎠ (66) = a1 ak+2 µk+1 ˜µk+1 ak+1 + ak+1 ak a1 ak+2 ˜µk˜µk+1 µk µ′ k µk+1 µ′ k+1 (67) = a₁ a_k+2 (68)

4.4 norms and expectation values with mpdos 39 And N[_ak] 1ak+1= ∑ ak N[_ak−1] 1ak ⎛ ⎝ ∑ µkµ˜k Wµkµ˜k akak+1δ µk ˜ µk ⎞ ⎠ (69) = a1 ak+1 µk ˜µk a_{k (70)} = _{a1 ak+1} (71)

The site index k runs from k = 2 to k = N − 1 for the C blocks, while k ∈ [2, N] for the N blocks. The expectation value is then given by:

⟨ ˆO⟩ = C[N −1]/N[N ] (72)

For the case of local operators, we only need to change the C’s:

C[1]_a 1a2 = ∑ µ1,µ˜1,µ′1 O[1]_µ˜ 1µ′1W µ1µ˜1_δµ1 µ′ 1 = a1 a2 ˜µ1 µ1 µ′ 1 =

a

a

4.4 norms and expectation values with mpdos 40 C[_ak] 1ak+1= ∑ ak C[_ak−1] 1ak ⎛ ⎝ ∑ µk,˜µk Wµkµ˜k akak+1δ µk ˜ µk ⎞ ⎠ + ∑ ak N[_ak−1] 1ak ⎛ ⎜ ⎝ ∑ µk,µ˜k,µ′ k O[k] ˜ µkµ′_kW µkµ˜k akak+1δ µk µ′ k ⎞ ⎟ ⎠ (74) = a1 ak+1 µk ˜µk ak ₊ a1 ak+1 ˜µk ak µk µ′ k (75) = _{a1 ak+1} (76)