Hand-waving and interpretive dance: an introductory course on tensor networks

Tensors are a generalisation of vectors and matrices. A d-dimensional vector can be considered an element of $\mathbb{C}^d$, and a $n\times m$-dimensional matrix an element of $\mathbb{C}^{n\times m}$. Correspondingly a rank-r tensor of dimensions $d_1\times\cdots\times d_r$ is an element of $\mathbb{C}^{d_1\times\cdots\times d_r}$. We can clearly see that scalars, vectors and matrices are all therefore rank 0, 1 and 2 tensors respectively.

In tensor network notation a single tensor is simply represented by a geometric shape with legs sticking out of it, each corresponding to an index, analogous to the indices of Einstein notation. For example a rank-four tensor R can be represented as

In some contexts the shape used and direction of the legs can imply certain properties of the tensor or index—for a general network however, neither carry any special significance. When representing quantum states, it is often convenient to use the direction of legs to denote whether the corresponding vectors live in the Hilbert space ('kets') or its dual ('bras'). By adhering to this convention, certain prohibited contractions can be easily disallowed, such as a contraction between two kets. This is notationally analogous to the convention of upper and lower denoting co- and contra-variant indices in the Einstein or Penrose notation (a specialised form of TNN) employed in the study of general relativity or quantum field theory.

Because quantum mechanics, in contrast to general relativity, is complex, care has to be taken with complex conjugation. This is usually indicated either by explicitly labelling the tensor or adopting some index convention, such as flipping a network (upward and downward legs being exchanged) carrying an implicit conjugation.

The main advantage in TNN comes in representing tensors that are themselves composed of several other tensors. The two main operations we will consider are those of the tensor product and trace, typically used in the joint operation of contraction. As well as these two operations, the rank of a tensor can be altered by grouping/splitting indices.

2.2.1. Tensor product.

The first operation we will consider is the tensor product, a generalisation of the outer product of vectors. The value of the tensor product on a given set of indices is the element-wise product of the values of each constituent tensor. Explicitly written out in index notation, the binary tensor product has the form:

$$\begin{align} [A\otimes B]_{i_1,\ldots,i_r,j_1,\ldots,j_s}:= A_{i_1,\ldots,i_r}\cdot B_{j_1,\ldots,j_s}. \nonumber \end{align} \tag{ 2.2 }$$

Diagrammatically the tensor product is simply represented by two tensors being placed next to each other. As such the value of a network containing disjoint tensors is simply the product of the constituent values.

2.2.2. Trace.

The next operation is that of the (partial) trace. Given a tensor A, for which the xth and yth indices have identical dimensions ($d_x=d_y$), the partial trace over these two dimensions is simply a joint summation over that index:

$$\begin{align} [{\rm Tr}_{x,y}A]_{i_1,\ldots,i_{x-1},i_{x+1},\ldots,i_{y-1},i_{y+1},\ldots,i_r}=\sum_{\alpha=1}^{d_x}A_{i_1,\ldots,i_{x-1},\alpha,i_{x+1},\ldots,i_{y-1},\alpha,i_{y+1},\ldots,i_r}. \nonumber \end{align} \tag{ 2.4 }$$

Similar to Einstein notation, this summation is implicit in TNN, indicated by the corresponding legs being joined. An advantage over Einstein notation is that these summed-over indices need not be named, making the notation less clunky for large networks. For example, consider tracing over the two indices of a rank-3 tensor:

One property of the trace we can trivially see from this notation is that of its cyclic property. By simply sliding one of the matrices around—which only changes the placement of the tensors in the network, and therefore not the value—we can cycle the matrices around (being careful of transpositions), proving ${\rm Tr}(AB)={\rm Tr}(BA)$.

Whilst this serves as a trivial example, the higher rank equivalents of this statement are not always so obvious, and the fact that these properties hold 'more obviously' in TNN is often useful.

2.2.3. Contraction.

The most common tensor operation used is contraction, corresponding to a tensor product followed by a trace between indices of the two tensors. An example would be the contraction between two pairs of indices of two rank-3 tensors, which is drawn as:

Familiar examples of contraction are vector inner products, matrix-vector multiplication, matrix-matrix multiplication, and the trace of a matrix:

Conventional	Einstein	TNN
$\langle \vec{x}, \vec{y}\rangle$	$x_\alpha y^\alpha$
$M\vec{v}$	$M^{\alpha}_{~\beta}v^\beta$
AB	$A^{\alpha}_{~\beta}B^{\beta}_{~\gamma}$
${\rm Tr}(X)$	$X^{\alpha}_{~\alpha}$

2.2.4. Grouping and splitting.

Rank is a rather fluid concept in the study of tensor networks. The space of tensors $\mathbb{C}^{a_1\times\cdots\times a_n}$ and $\mathbb{C}^{b_1\times\cdots\times b_{m}}$ are isomorphic as vector spaces whenever the overall dimensions match ($\prod_i a_i=\prod_i b_i$). Using this we can extend concepts and techniques only previously defined for vectors and matrices to all tensors. To do this, we can group or split indices to lower or raise the rank of a given tensor respectively.

Consider the case of contracting two arbitrary tensors. If we group together the indices which are and are not involved in this contraction, this procedure simply reduces to matrix multiplication:

It should be noted that not only is this reduction to matrix multiplication pedagogically handy, but this is precisely the manner in which numerical tensor packages perform contraction, allowing them to leverage highly optimised matrix multiplication code.

Owing to the freedom in choice of basis, the precise details of grouping and splitting are not unique. One specific choice of convention is the tensor product basis, defining a basis on the product space simply given by the product of the respective bases. The canonical use of tensor product bases in quantum information allows for the grouping and splitting described above to be dealt with implicitly. Statements such as $\left\vert 0\right\rangle\otimes \left\vert 1\right\rangle\equiv \left\vert 01\right\rangle$ omit precisely this grouping: notice that the tensor product on the left is a $2\times 2$ dimensional matrix, whilst the right hand-side is a 4-dimensional vector. The 'tensor product' used in quantum information is often in fact a Kronecker product, given by a true tensor product followed by just such a grouping.

More concretely, suppose we use an index convention that can be considered a higher-dimensional generalisation of column-major ordering. If we take a rank n  +  m tensor, and group its first n indices and last m indices together to form a matrix

$$\begin{align} T_{I,J}:=T_{i_1,\dots,i_n;j_{1},\dots,j_{m}} \nonumber \end{align} \tag{ 2.9 }$$

where we have defined our grouped indices as

$$\begin{align} I:=i_1+d_1^{(i)}\cdot i_2+d_1^{(i)}d_2^{(i)}\cdot i_3+\dots+d_1^{(i)}\dots d_{n-1}^{(i)}\cdot i_{n}, \nonumber \end{align} \tag{ 2.10 }$$

$$\begin{align} J:=j_1+d_1^{(\,j)}\cdot j_2+d_1^{(\,j)}d_2^{(\,j)}\cdot j_3+\dots+d_1^{(\,j)}\dots d_{m-1}^{(\,j)}\cdot j_{m}, \nonumber \end{align} \tag{ 2.11 }$$

where $d_{x}^{(i)}(d_{x}^{(\,j)})$ is the dimension of the xth index of type i(j). When such a grouping is given, we can now treat this tensor as a matrix, performing standard matrix operations.

An important example is the singular value decomposition (SVD), given by $T_{I, J}=\sum_{\alpha}U_{I, \alpha}S_{\alpha, \alpha}\bar{V}_{J, \alpha}$. By performing the above grouping, followed by the SVD, and then splitting the indices back out, we get a higher dimensional version of the SVD

$$\begin{align*} T_{i_1,\dots,i_n;j_{1},\dots,j_{m}}=\sum_{\alpha}U_{i_1,\dots,i_n,\alpha}S_{\alpha,\alpha}\bar{V}_{j_1,\dots,j_m,\alpha}. \nonumber \end{align*}$$

So long as we choose them to be consistent, the precise method by which we group and split is immaterial in this overall operation. As a result we will keep this grouping purely implicit, as in the first equality equation (2.8). This will be especially useful for employing notions defined for matrices and vectors to higher rank objects, implicitly grouping then splitting. Graphically the above SVD will simply be denoted

where U and V are isometric ($U^\dagger U=V^\dagger V=\mathbb{1}$) across the indicated partitioning, and where the conjugation in $V^\dagger$ is included for consistency with conventional notation and also taken with respect to this partitioning. We will refer to such a partitioning of the indices in to two disjoint sets as a bisection of the tensor.

Aside 1 (Why do we care so much about the singular value decomposition?). One of the main uses of tensor networks in quantum information is representing states which belong to small but physically relevant corners of an otherwise prohibitively large Hilbert space, such as low-entanglement states. The central backbone of this idea is that of low matrix-rank approximations. Suppose we have some matrix, and we want the ideal low matrix-rank approximation thereof. Eckart and Young [1] showed that if we measure error in the Frobenius norm, then trimming the singular value decomposition is an ideal approximation. Specifically take $X=USV^\dagger$ to be the SVD of X, then the trimmed version of X is given by

$$\begin{align*} X^{(k)}=US^{(k)}V^\dagger \nonumber \end{align*}$$

where S^(k) has had all but the largest k singular values set to zero (i.e. has matrix-rank k), then Eckart-Young theorem says that $\vert{X-X^{(k)}}\vert_F\leqslant \vert{X-Y}\vert_F$ for all Y of matrix-rank k. Mirsky further generalised this result in [2] to show optimality in all unitarily invariant norms. Whenever we use the term trim, we are referring to this very method of low-rank approximation.

Combining the above tensor operations, we can now give a single definition of a tensor network. A tensor network is a diagram which tells us how to combine several tensors into a single composite tensor. The rank of this overall tensor is given by the number of unmatched legs in the diagram. The value for a given configuration of external indices, is given by the product of the values of the constituent tensors, summed over all internal index labellings consistent with the contractions. A generic example of this is given below:

Whilst tensor networks are defined in such a way that their values are independent of the order in which the constituent tensors are contracted, such considerations do influence the complexity and practicality of such computations. Tensor networks can be contracted by beginning with a single tensor and repeatedly contracting it against tensors one-at-a-time.

The order in which tensors are introduced and contracted is known as a bubbling. As the bubbling is performed the network is swallowed into the stored tensor, until only the result remains.

Many networks admit both efficient and inefficient bubblings, highlighting the need for prudence when planning out contractions. Take for example a ladder-shaped network (we'll see a few of these in the following lectures). One bubbling we may consider is to contract along the top of the ladder, then back along the bottom. Showing both this bubbling, as well as the partially contracted tensor that is kept in memory (in red), we see this bubbling looks like:

The scaling of this procedure is however quite unfavourable; consider a ladder of length n. At the midpoint of this contraction, when the top has been contracted, the tensor being tracked has rank n, and thus the number of entries is scaling exponentially with n. As such the memory and time footprints of this contraction are also exponential, rendering it infeasible for large n. If however we contract each rung in turn, the tracked tensor has a rank never more than 3, giving constant memory and linear time costs.

The memory footprint at any step during the contraction corresponds to the product of the dimensions of each leg passing through the boundary of the contracted region (see the red legs in equation (2.18)). Whilst the above ladder arrangement possesses both good and bad bubblings, some networks possess an underlying graph structure that does not admit any efficient contraction ordering. A good example of this is the 2D grid; due to the 2D structure of this lattice, it is clear that the contracted region must, somewhere near the middle of the contracting procedure, have a perimeter on the order of $\sqrt{n}$ where n is the number of tensors. As a result such contractions generically take exponential time/memory to perform. An example of a high cost step during such a bubbling is shown below, with the prohibitively large perimeter indicated by the red legs.

Although the bubblings we have depicted here involve picking a single tensor and contracting others into it one-by-one, this will frequently not be the most efficient order; often a multibubbling approach is faster.

Pfeifer et al [3] provides code which allows for finding optimal bubbling order for networks of up to 30–40 tensors. This code interfaces with that provided in [4] and [5], providing a complete tensor network package.

Above we've described that there exist networks which stymie the specific contraction procedures we've outlined. In this section we'll see that there also exist networks for which there are complexity theoretic obstructions which do not allow for any contraction procedure to be efficient.

We will now consider the computational complexity associated with tensor network contractions. Whilst all of the tensor networks we will consider in later lectures constitute memory-efficient representations of objects such as quantum states, not all permit efficient manipulation. This demonstrates that how one wishes to manipulate a tensor network is an important part of considering them as ansätze.

Whilst algorithms which can speed up tensor network contractions by optimising the bubbling used [3–5], as discusssed above, the underlying computational problem is $\mathsf{NP}$-complete [6, 7]

Even ignoring the specific bubbling used, the complexity of the overall contraction procedure can also be shown to be prohibitive in general. Consider a network made from the binary tensors e and n. The value of e is 1 if and only if all indices are identical, and zero otherwise, whilst n has value 1 if and only if all legs differ and 0 otherwise. Take an arbitrary graph, and construct a tensor network with an e tensor at each vertex and n tensor in the middle of each edge, with the connectedness inherited from the graph.

By construction, the non-zero contributions to the above tensor network correspond to an assignment of index values to each vertex (enforced by e) of the original graph, such that no two neighbouring vertices share the same value (enforced by n). If each index is q-dimensional this is a vertex q-colouring of the graph, and the value of the tensor network corresponds to the number of such q-colourings. As determining the existence of a q-colouring is an $\mathsf{NP}$-complete problem [8], contracting this graph is therefore $\mathsf{\#P}$-complete [9]. Indeed similar constructions exist for tensor networks corresponding to $\#\mathsf{SAT}$ and other $\#\mathsf{P}$-complete problems [10]. As we will see later in section 7, there also exists a quantum hardness result which shows approximate contraction to be $\mathsf{Post{\rm{-}}BQP}$-hard, putting it inside a class of problems not believed to be efficiently solvable on even a quantum computer.

Problems 1. Solutions in supplementary material (stacks.iop.org/JPhysA/50/223001/mmedia).

• 1.  
  Consider the following tensors, in which all indices are three-dimensional, indexed from 0:

  

  

  Calculate the value of the following tensor network:

  
• 2.  
  In this question we are going to consider expanding out a contraction sequence, in a manner which would be needed when coding up contractions. Given a network, and an associated bubbling, we wish to write out a table keeping track of the indices of the current object, the tensor currently being contracted in, the indices involved in that contraction, and new indices left uncontracted. For example for the network

  

  where the bubbling is performed in alphabetical order, then the table in question looks like

  Current	Tensor	Contract	New
  —	A	—	$\alpha, \beta$
  $\alpha, \beta$	B	α	γ
  $\beta, \gamma$	C	$\beta, \gamma$	δ

  For the tensor network

  

  construct a corresponding table, where contraction is once again done in alphabetical order.
• 3.  
  • (a)  
    Calculate the contraction of the tensor network in equation (2.19) for bond dimension 3, i.e. calculate the number of three-colourings of the corresponding graph.
  • (b)  
    Using the e and n tensors from section 2.5, come up with a construction for a tensor network which gives the number of edge colourings. For a planar graphs, construct an analogous network to count face colourings.
  • (c)  
    Using tensor networks, determine the minimum number of colours required to vertex and edge colour the below graph (known as the chromatic number and index respectively).

    
• 4.  
  Much like the singular value decomposition, given a bisection of the indices we can consider norms of tensors.

  • (a)  
    Does the operator norm depend on the bisection, i.e. are the operator norms across any two bisections of the same tensor necessarily equal?
  • (b)  
    What about the Frobenius norm? If they can differ, give an example, if not draw a tensor network diagram that shows it to be manifestly independent of bisection.
• 5.  
  Write out the Einstein notation corresponding to the network in equation (8.1).

The Bell basis forms a convenient orthonormal set of two qubit states that exhibit maximal entanglement. The standard notation for this basis is

$$\begin{align*} \vert {\Phi^\pm} \rangle :=\Big(\vert {0} \rangle \otimes\vert {0} \rangle \pm\vert {1} \rangle \otimes\vert {1} \rangle \Big)/\sqrt{2} \qquad{\rm{and}}\qquad \vert {\Psi^\pm} \rangle :=\Big(\vert {0} \rangle \otimes\vert {1} \rangle \pm\vert {1} \rangle \otimes\vert {0} \rangle \Big)/\sqrt{2}. \nonumber \end{align*}$$

The first of this basis, $\vert {\Phi^+} \rangle$, we shall denote $\vert {\Omega} \rangle$ and simply refer to as the Bell state. Thought of as a matrix, $\Omega$ is proportional to the one qubit identity,

$$\begin{align} \vert {\Omega} \rangle {=\frac{1}{\sqrt{2}}}\left( \begin{array}{c} 1\nonumber \\ 0\nonumber \\ \hdashline0\nonumber \\ 1 \end{array} \right) \quad \mathop{\leftrightharpoons}\limits_{\rm Matricise}^{\rm Vectorise}\quad {\frac{1}{\sqrt{2}}}\left(\begin{array}{cc} 1&0\nonumber \\0&1 \end{array}\right)=I{/\sqrt{2}}. \nonumber \end{align} \tag{ 3.1 }$$

In tensor network notation, this is represented simply as a line connecting its two legs.

Next we will define $\Omega(O)$ to be the vectorisation of an operator O, such that $\vert {\Omega(O)} \rangle =(O\otimes I)\vert {\Omega} \rangle$.

Given this definition, we can see that the Bell basis simply corresponds to a vectorisation of the Pauli operators

$$\begin{align*} \vert {\Phi^+} \rangle =\vert {\Omega(I)} \rangle,\quad \vert {\Phi^-} \rangle =\vert {\Omega(Z)} \rangle,\quad \vert {\Psi^+} \rangle =\vert {\Omega(X)} \rangle,\quad \vert {\Psi^-} \rangle \propto\vert {\Omega(Y)} \rangle. \nonumber \end{align*}$$

Thus we see that the Bell basis is intimately linked to the Pauli operators, with the Euclidean inner product on Bell basis states corresponding to the Hilbert–Schmidt inner product on Paulis.

Given this notation for the Bell basis, we can now understand Quantum Teleportation in TNN. The idea here is for two parties (Alice and Bob, say) to share a Bell state. Given this shared resource of entanglement, we then allow Alice to perform local operations on her half of the pair, and an arbitrary fiducial qubit. After transmitting only two classical bits, Bob can then correct his half of the pair such that he recovers the state of the original fiducial qubit, successfully teleporting the data within.

The procedure for teleportation goes as follows. First Alice performs a projective measurement in the Bell basis on both the fiducial qubit and her Bell qubit. The result of this measurement is then (classically) transmitted to Bob, requiring two communication bits. Bob then performs the corresponding Pauli on his Bell qubit, correcting the influence of the measurement. Taking the fiducial state to be $\vert {\psi} \rangle$, and supposing the measurement outcome corresponds to $\vert {\Omega(\,p)} \rangle$, then this procedure gives Bob a final state of $\vert {\phi} \rangle =\vert {\psi} \rangle /2$:

$$\begin{align} \vert {\phi_B} \rangle =\overbrace{\Bigl(\,p_{B}\Bigr)}^{{\rm{Correction}}}\overbrace{\Bigl(\langle{\Omega_{A_1A_2}(\,p)}\vert\Bigr)}^{{\rm{Teleportation}}}\overbrace{\Bigl(\vert {\psi_{A_1}} \rangle \otimes\vert {\Omega_{A_2B}} \rangle \Bigr)}^{{\rm{Setup}}}=\vert {\psi_B} \rangle /2, \nonumber \end{align} \tag{ 3.4 }$$

where A₁ and A₂ correspond to the single qubit registers of Alice, and B to Bob's qubit. In tensor network notation this can be clearly seen:

$$\begin{align} = \vert {\psi} \rangle /2, \nonumber \end{align} \tag{ 3.7 }$$

where the dashed line indicates the physical separation of the two parties.

As such we can see that $\vert {\psi} \rangle$ is correctly transmitted for any measurement outcome p, each of which is seen with probability 1/4. Thus we see that in spite of the non-deterministic intermediary states, the overall procedure is deterministic. Analogous procedures can work for p being elements of any set of operators which are orthonormal with respect to the Hilbert–Schmidt inner product, e.g. higher dimensional Paulis.

3.2.1. Gate teleportation.

The idea behind gate teleportation is similar to regular teleportation, but utilises a general maximally entangled state instead of the Bell state specifically. Suppose we prepare a maximally entangled state $\vert {\Omega(U^T)} \rangle$ corresponding to a unitary U, and post select on a Bell basis measurement of $\vert {\Omega(\,p)} \rangle$, followed by a correcting unitary C_p, then Bob ends up with the state:

$$\begin{align} \vert {\phi} \rangle ={\overbrace{\Bigl(C_p\Bigr)}^{{\rm{Correction}}}\overbrace{\Bigl(\langle{\Omega_{A_1A_2}(\,p)}\vert\Bigr)}^{{\rm{Teleportation}}}\overbrace{\Bigl(\vert {\psi_{A_1}} \rangle \otimes\vert {\Omega_{A_2B}\big(U^T)} \rangle \Bigr)}^{{\rm{Setup}}}} \nonumber \end{align} \tag{ 3.8 }$$

$$\begin{align} = C_pUp^\dagger \vert {\psi} \rangle /2 \nonumber \end{align} \tag{ 3.11 }$$

If we take $C_p:=U pU^\dagger$ then Bob receives $U\vert {\psi} \rangle$ for all measurement outcomes, i.e. $\vert {\phi} \rangle \propto U\vert {\psi} \rangle$. If U is a Clifford operator⁵, this correction is also a Pauli, making the procedure no more resource intensive in terms of the gates used than standard teleportation.

An example of where this is useful is in the case where Paulis can be reliably performed, but Cliffords can only be applied non-deterministically. Gate teleportation allows us to prepare the $\vert {U^T} \rangle$ first, simply retrying the non-deterministic procedure until it succeeds. Once this has succeeded, we can use gate teleportation to apply this unitary on the data state using only Pauli operations. As such we can avoid needing to apply non-deterministic gates directly on our target state, endangering the data stored within.

For a given mixed state ρ, a purification is a pure state $\vert {\psi} \rangle$ which is extended into a larger system (the added subsystem is known as the purification system), such that the reduced density on the original system is ρ. One such purification is given by $\vert {\psi} \rangle \propto(\sqrt{\rho}\otimes I)\vert {\Omega} \rangle =\vert {\Omega(\sqrt{\rho})} \rangle$, which can be simply seen by considering the corresponding tensor networks. The definition of the state is

which gives a reduced density of

By dimension counting, it can be shown that the above purification is unique up to an isometric freedom on the purification system, i.e. all purifications are of the form $\left(\sqrt{\rho}\otimes U\right)\vert {\Omega} \rangle$ where $U^\dagger U=\mathbb{1}$. Equivalently all purifications can be considered to be proportional to $(\sqrt{\rho}\otimes I)\vert {\Omega} \rangle$, where $\vert {\Omega} \rangle$ is some maximally entangled state other than the Bell state.

Stinespring's Theorem says that any quantum channel $\mathcal{E}$—a completely positive trace preserving (CPTP) map—can be expressed as a unitary map V acting on a larger system followed by a partial trace, i.e.

$$\begin{align} \mathcal{E}(\rho)={\rm Tr}_1 \left[V^\dagger\left(\rho\otimes \vert {0} \rangle \langle{0}\vert\right) V\right]. \nonumber \end{align} \tag{ 3.14 }$$

Physically this means that dynamics of an open system is equivalent to those of a subsystem of a larger, closed system—the founding tennet of the Church of the Larger Hilbert Space. Any CPTP map can be represented by a set of Kraus operators K_i such that

$$\begin{align} \mathcal{E}(\rho)=\sum_i K_i^\dagger \rho K_i{\rm{~where~}}\sum_i K_iK_i^\dagger = I. \nonumber \end{align} \tag{ 3.15 }$$

In TNN this looks like

where the transposition in the Hermitian conjugate is done with respect to the horizontal legs, and the upper leg corresponds to the virtual index i.

Next we define the tensor U as

where we can see that U is an isometry ($U^\dagger U=I$), which we can think of as a unitary V with an omitted ancilla

Using this, and partial tracing over the upper index, we get the Stinespring Dilation Theorem as desired:

$$\begin{align} ={\rm Tr}_1 \left[V^\dagger\left(\rho\otimes \vert {0} \rangle \langle{0}\vert\right) V\right] \nonumber \end{align} \tag{ 3.23 }$$

Problems 2. Solutions in supplementary material (stacks.iop.org/JPhysA/50/223001/mmedia).

• 1.  
  Consider the inverse of teleportation. Alice wishes to send classical bits to Bob, and possesses a quantum channel through which she can send Bob qubits. How many bits of information can be communicated in a single qubit? For simplicity consider the case where Bob can only perform projective measurements.
• 2.  
  Suppose Alice and Bob initially shared a Bell pair. Does this pre-shared entanglement resource boost the amount of classical information that can be successfully communicated, and if so by how much? Hint: Notice that the four possible Bell states differ by a Pauli acting on a single qubit.

In addition to the above construction, MPS can (equivalently) be viewed as a special case of the projected entangled pair states (PEPS) construction [18, 23, 24]. This proceeds by laying out entangled pair states $\vert {\phi} \rangle$ on some lattice and applying some linear map $\mathcal{P}$ between pairs

where

is the chosen entangled pair. In section 7, we will generalise this construction to arbitrary dimensions and arbitrary lattices.

It is clear that this construction is equivalent to the tensor network construction by letting $\vert {\phi} \rangle =\sum_{j=0}^{d-1}\vert {dd} \rangle$. We can write the linear map $\mathcal{P}$ as

$$\begin{align} \mathcal{P}=\sum_{i,\alpha,\beta}A_{i;\alpha,\beta}\vert {i} \rangle \langle{\alpha\beta}\vert. \nonumber \end{align} \tag{ 4.13 }$$

The tensor A is exactly the MPS tensor introduced above, and the choice of entangled pair ensures that the A tensor corresponding to a pair of PEPS 'projectors' applied to the Bell state above is exactly the contraction of the corresponding A tensors:

$$\begin{align} \mathcal{P}^{(1)}\otimes\mathcal{P}^{(2)}\vert {\phi} \rangle _{2,3}=\sum_{i_1,i_2;\alpha_1,\beta_1,\alpha_2,\beta_2,j}A^{(1)}_{i_1;\alpha_1,\beta_1}A^{(2)}_{i_2;\alpha_2,\beta_2}\vert {i_1i_2} \rangle \langle{\alpha_1\beta_1\alpha_2\beta_2}\vert(\mathbb{1}\otimes\vert{\,jj} \rangle \otimes\mathbb{1}) \nonumber \end{align} \tag{ 4.14 }$$

$$\begin{align} =\sum_{i_1,i_2;\alpha_1,\beta_1,\beta_2}A^{(1)}_{i_1;\alpha_1,\beta_1}A^{(2)}_{i_2;\beta_1,\beta_2}\vert {i_1i_2} \rangle \langle{\alpha_1\beta_2}\vert. \nonumber \end{align} \tag{ 4.15 }$$

Thus, we see that the two descriptions are equivalent, and interchanged through the applications of local unitaries to the virtual indices of A or equivalently changing the maximally entangled pair in the PEPS.

We note that this should not generally be seen as a practical preparation procedure. Generically the PEPS tensors will map states down into a non-trivial subspace, with the physical implementation of this requiring post-selected measurements. If one of these fails, we need to go back and begin the construction from the start, meaning this procedure is not generally scalable.

4.2.1. Product state.

Let

$$\begin{align} A_0=\left(\begin{array}{cc} 1 \end{array}\right), \qquad \qquad A_1&=\left(\begin{array}{cc} 0 \end{array}\right). \nonumber \end{align} \tag{ 4.16 }$$

This gives the state $\vert {00\ldots0} \rangle$, as does

$$\begin{align} A_0=\left(\begin{array}{cc} 1&0\nonumber \\ 0&0 \end{array}\right), \qquad \qquad A_1&=\left(\begin{array}{cc} 0&0\nonumber \\ 0&0 \end{array}\right). \nonumber \end{align} \tag{ 4.17 }$$

4.2.2. W state.

What state do we get when we set

$$\begin{align} A_0=\left(\begin{array}{cc} 1&0\nonumber \\ 0&1 \end{array}\right), \qquad \qquad A_1&=\left(\begin{array}{cc} 0&1\nonumber \\ 0&0 \end{array}\right),\label{eqn:w-state} \nonumber \end{align} \tag{ 4.18 }$$

and we choose the boundary conditions of the MPS to be

We have $A_0A_0=A_0$, $A_0A_1=A_1$, $A_1^2=0$ and ${\rm Tr}[A_1X]=1$, so we get

$$\begin{align} \vert {W} \rangle =\sum_{j=1}^{N} \vert {000\ldots01_j000\ldots0} \rangle, \nonumber \end{align} \tag{ 4.20 }$$

the W-state [18].

4.2.3. GHZ state.

If we choose $\vert {\phi} \rangle =\vert {00} \rangle +\vert {11} \rangle$ and $\mathcal{P}=\vert {0} \rangle \langle{00}\vert+\vert {1} \rangle \langle{11}\vert$, or the equivalent MPS tensor

$$\begin{align} A_0=\left(\begin{array}{cc} 1&0\nonumber \\ 0&0 \end{array}\right), \qquad \qquad A_1&=\left(\begin{array}{cc} 0&0\nonumber \\ 0&1 \end{array}\right), \nonumber \end{align} \tag{ 4.21 }$$

then we get the Greenberger-Horne-Zeilinger (GHZ) state [18]

$$\begin{align} \vert {GHZ} \rangle =\vert {00\ldots0} \rangle +\vert {11\ldots1} \rangle. \nonumber \end{align} \tag{ 4.22 }$$

4.2.4. AKLT state.

Suppose we wish to construct an SO(3) symmetric spin-1 state [21, 22, 25]. Let $\vert {\phi} \rangle =\vert {01} \rangle -\vert {10} \rangle$ be the SO(3) invariant singlet state. Let $\mathcal{P}:\mathbb{C}^{2\times 2}\to \mathbb{C}^3$ be the projector onto the spin-1 subspace

$$\begin{align} \mathcal{P}=\vert {\tilde{1} \rangle } \langle{00}\vert+\vert {\tilde{0} \rangle } \frac{\langle{01}\vert+\langle{10}\vert}{\sqrt{2}}+\vert {-\tilde{1}} \rangle \langle{11}\vert. \nonumber \end{align} \tag{ 4.23 }$$

The advantage is that the spin operators on the corresponding systems pull through $\mathcal{P}$, meaning it commutes with rotations. Let $(S_x, S_y, S_z)$ be the spin vector on the spin-1 particle, and $(X_i, Y_i, Z_i)/2$ the spin vector on the ith qubit, then this means:

$$\begin{align} S_Z\mathcal{P}=\left(\vert {\tilde{1} \rangle} \langle{\tilde{1}}\vert-\vert {-\tilde{1} \rangle} \langle{-\tilde{1}}\vert\right)\left(\vert {\tilde{1} \rangle } \langle{00}\vert+\vert {\tilde{0} \rangle } \frac{\langle{01}\vert+\langle{10}\vert}{\sqrt{2}}+\vert {-\tilde{1}} \rangle \langle{11}\vert\right) \nonumber \end{align} \tag{ 4.24 }$$

$$\begin{align} =\vert {\tilde{1}} \rangle \langle{00}\vert-\vert {-\tilde{1}} \rangle \langle{11}\vert\nonumber \end{align} \tag{ 4.25 }$$

$$\begin{align} =\mathcal{P}\frac{Z_1+Z_2}{2}\nonumber \end{align} \tag{ 4.26 }$$

$$\begin{align} S_X\mathcal{P}=\frac{\vert {\tilde{0} \rangle } \left(\langle{\tilde{1}}\vert+\langle{-\tilde{1}}\vert\right)+\left(\vert {\tilde{1} \rangle } +\vert {-\tilde{1} \rangle } \right)\langle{\tilde{0}}\vert}{\sqrt{2}}\left(\vert {\tilde{1} \rangle } \langle{00}\vert+\vert {\tilde{0} \rangle } \frac{\langle{01}\vert+\langle{10}\vert}{\sqrt{2}}+\vert {-\tilde{1} } \rangle \langle{11}\vert\right)\nonumber \end{align} \tag{ 4.27 }$$

$$\begin{align} =\left(\frac{\vert {\tilde{0} \rangle } \left(\langle{00}\vert+\langle{11}\vert\right)}{\sqrt{2}}+\frac{\left(\vert {\tilde{1} \rangle } +\vert {-\tilde{1}} \rangle \right)\left(\langle{01}\vert+\langle{10}\vert\right)}{2}\right)\nonumber \end{align} \tag{ 4.28 }$$

$$\begin{align} =\mathcal{P}\frac{X_1+X_2}{2}, \nonumber \end{align} \tag{ 4.29 }$$

with the same holding for S_Y. Thus the state obtained after this projection is fully SO(3) symmetric, but has a nontrivial entanglement structure (which would not be obtained if the state was simply a singlet at each site for example).

This state has many interesting properties. We can write a 2-local Hamiltonian for which this is the ground state. Let $\Pi_2$ be the projector onto the spin-2 subspace of a pair of spin-1 particles. This operator has eigenvalues {0,1}. $\Pi_2$ annihilates an adjacent pair of spin-1 particles, since they are built from two spin-1/2s and a spin-0, so have no overlap with the spin-2 subspace. It is simple to check that on periodic boundary conditions the ground state of $H=\sum \Pi_2$ is unique (and gapped).

If we examine the action of rotations about the three axes of the spin-1, we see that

$$\begin{align} R_{\hat{n}}(\theta)\mathcal{P}=\mathcal{P}R_{\hat{n}}(\theta)\otimes R_{\hat{n}}(\theta). \nonumber \end{align} \tag{ 4.30 }$$

In particular, $R_{\hat{x}}(\pi)\mapsto -XX$, $R_{\hat{y}}(\pi)\mapsto -YY$, $R_{\hat{z}}(\pi)\mapsto -ZZ$. In section 5 we will see that this tells us the AKLT state is in a nontrivial symmetry protected topological (SPT) phase.

4.2.5. Cluster state.

It is convenient to write a bond dimension 2 MPS for this state where a physical site contains a pair of spins. Let

$$\begin{eqnarray} A_{00}=\left(\begin{array}{cc} 1&0\nonumber \\ 1&0 \end{array}\right) \qquad A_{01}=\left(\begin{array}{cc} 0&1\nonumber \\ 0&1 \end{array}\right) \qquad A_{10}=\left(\begin{array}{cc} 1&0\nonumber \\ -1&0 \end{array}\right) \qquad A_{11}=\left(\begin{array}{cc} 0&-1\nonumber \\ 0&1 \end{array}\right), \nonumber \end{eqnarray} \tag{ 4.31 }$$

or equivalently the map from virtual to physical spin-1/2 particles

$$\begin{align} \mathcal{P}&=\left(\begin{array}{cccc} 1&0&1&0\nonumber \\ 0&1&0&1\nonumber \\ 1&0&-1&0\nonumber \\ 0&-1&0&1 \end{array}\right), \nonumber \end{align} \tag{ 4.32 }$$

where the entangled pairs are in the Bell state $\vert {\phi} \rangle =\vert {00} \rangle +\vert {11} \rangle$. The map $\mathcal{P}$ corresponds to the circuit

Notice in this case our PEPS tensor $\mathcal{P}$ simply corresponds to unitary circuit. As such this is one of the exceptional cases in which the PEPS description can be considered a scalable preparation procedure.

Given an explicit MPS description of this state, we can now back out a Hamiltonian for which it is a ground state, allowing us to infer certain properties.

The initial state is constructed from entangled pairs $\prod\vert {\phi} \rangle _{2j, 2j+1}$, and is the unique ground state of the Hamiltonian

$$\begin{align} H=-\sum_j \left(X_{2j}X_{2j+1}+Z_{2j}Z_{2j+1}\right). \nonumber \end{align} \tag{ 4.34 }$$

Applying the circuit (between Bell pairs with first qubit odd and second even), we see that this transforms to

$$\begin{align} H^{\prime}=-\sum_j \left(Z_{2j-1}X_{2j}Z_{2j+1}+Z_{2j}X_{2j+1}Z_{2j+2}\right) \nonumber \end{align} \tag{ 4.35 }$$

$$\begin{align} =-\sum_k Z_{k-1} X_k Z_{k+1}. \nonumber \end{align} \tag{ 4.36 }$$

This is precisely the cluster state Hamiltonian. The physical symmetry of this model is $\mathbb{Z}_2\times\mathbb{Z}_2$, where $S_1=\prod_j X_{2j-1}$ and $S_2=\prod_j X_{2j}$. Pushing this backwards through the circuit, we see that it is equivalent to act on the virtual spins with $S_1=\prod_j Z_{2j} Z_{2j+1}$ and $S_2=\prod_j X_{2j} X_{2j+1}$.

This action tells us that, just like the AKLT state, the cluster state possesses SPT order.

MPS form a vanishingly small corner of the full Hilbert space, and thus we cannot hope to use them to approximate arbitrary states. If physically relevant states correspond to those which can be well approximated by MPS, and MPS manifest the same non-generic properties as these physical states, then they represent an extremely useful tool to study these systems.

4.3.1. Decay of correlations.

We have already seen that MPS have bounded levels of entanglement, manifesting as strict area laws. We will now investigate the type of correlations which can be represented. Let $\mathcal{O}$ be some operator for which we wish to compute the two point correlator

$$\begin{align} \langle{\psi[A]}\vert\mathcal{O}_0\mathcal{O}_{j+1}\vert {\psi[A]} \rangle, \nonumber \end{align} \tag{ 4.37 }$$

where the subscript denotes the site at which the operator $\mathcal{O}$ is applied. Graphically this expectation value is written as:

We refer to the object

as the $\mathcal{O}$-transfer matrix. Note that we usually just refer to $\mathbb{E}_{\mathbb{1}}$ as the transfer matrix and simply denote it $\mathbb{E}$.

The correlator (in the thermodynamic limit) can then be written as

$$\begin{align} \langle{\psi[A]}\vert\mathcal{O}_0\mathcal{O}_{j+1}\vert {\psi[A]} \rangle ={\rm Tr}\left(\mathbb{E}^{\infty}\mathbb{E}_{\mathcal{O}_0}\mathbb{E}^{\,j}\mathbb{E}_{\mathcal{O}_{j+1}}\mathbb{E}^{\infty}\right) \nonumber \end{align} \tag{ 4.40 }$$

$$\begin{align} \propto V_{\rm L}^\dagger \mathbb{E}^{\,j}V_{\rm R}.\label{eqn:MPSCorrelator} \nonumber \end{align} \tag{ 4.41 }$$

where V_L and V_R are the dominant left and right eigenvectors of $\mathbb{E}$ respectively. The only change required when calculating longer range correlators is inserting higher powers of $\mathbb{E}$ in equation (4.41). The decay of correlators is therefore controlled by the eigenvalues of $\mathbb{E}$. We can normalise A so that the dominant eigenvalue of $\mathbb{E}$ is 1, with the rest lying inside the unit disk. Thus any correlator can either decay exponentially with distance or be constant. Thus we see that MPS can only capture states with exponentially decaying correlations [22].

4.3.2. Gauge freedom.

Not all MPS represent different physical states [18]. The set of transformations of the description (i.e. the MPS) which leaves the physical state invariant are known as gauge transformations. In the case of MPS, these correspond to basis transformations on the virtual level:

$$\begin{align} =\vert {\psi[B]} \rangle, \nonumber \end{align} \tag{ 4.44 }$$

where $B_j=MA_jM^{-1}$. Note that M is only required to have a left inverse, so can be rectangular and enlarge the bond dimension.

Another freedom is blocking. We can combine several MPS tensors $A_{i_1}, A_{i_2}, \ldots, A_{i_j}$ into a single effective tensor B_k, on a larger physical region

A number of canonical forms exist which partially gauge fix the MPS description. One of the most common is the left-isometric or left-canonical form (with right-isometric or right-canonical defined analogously). Here the MPS tensors obey

$$\begin{align} \sum_{j=0}^{d-1}A^\dagger_jA_j=\mathbb{1}_{D\times D}, \nonumber \end{align} \tag{ 4.45 }$$

This is most useful on open boundary systems where a simple algorithm exists to put any MPS into this form. It is frequently used in numerical applications, in particular when using variational minimisation to optimise an MPS description of a ground state (DMRG), a mixed left/right isometric form is used.

Putting an MPS into this form is a partial gauge fixing. The remaining freedom is that of a unitary⁷ on the virtual level, rather than general invertible matrix. This technique is heavily used in tensor network algorithms as a method of increasing numerical stability.

When we renormalise a system, we usually think about attempting to write down an effective model at a longer length scale which captures the low energy portion of the original model. This can be achieved by blocking sites together, then discarding degrees of freedom to ensure the description remains useful. In the MPS, blocking can be achieved by simply contracting tensors together. How to discard only high energy degrees of freedom is a challenging question. MPS allows us to avoid having to answer this question all together [26].

Since we care only about expectation values of operators, we can work entirely in the transfer matrix picture. Blocking sites together simply consists of taking products of transfer matrices

$$\begin{align} \tilde{\mathbb{E}}=\mathbb{E}\mathbb{E}\mathbb{E}\mathbb{E}\mathbb{E}\ldots\mathbb{E}, \nonumber \end{align} \tag{ 4.47 }$$

with sandwiched operators $\mathbb{E}_O$ being renormalised similarly. Note that the dimension of $\tilde{\mathbb{E}}$ remains D⁴ at all times, so we never need to worry about discarding degrees of freedom. We can also use transfer matrices formed from different MPS to get off-diagonal terms of the form $\left\langle{\psi}\vert{O}\vert{\phi}\right\rangle$.

As described above, an MPS can be used to represent a pure state. How is a mixed state represented in this language?

Let $\vert {\psi{[A]}} \rangle$ be some (pure) MPS state. We can write the density matrix corresponding to $\vert {\psi{[A]}} \rangle$ as

$$\begin{align} \rho[A]=\vert {\psi{[A]}} \rangle \langle{\psi{[A]}}\vert\nonumber \end{align} \tag{ 4.48 }$$

The reduced density matrix on some subset of spins R will therefore be represented as

$$\begin{align} \rho[A]_{\rm R}=\vert {\psi{[A]}} \rangle \langle{\psi{[A]}}\vert\nonumber \end{align} \tag{ 4.50 }$$

where we have used the left and right normal forms to bring in the boundary terms.

The above network is an example of what is referred to as matrix product operators (MPOs) [21, 27, 28]. The general form of MPOs we will be considering is

In addition to being used to represent density matrices, MPOs can be used to represent a large class of many body operators, including small depth quantum circuits and local Hamiltonians. For example, the transverse field Ising Hamiltonian

$$\begin{align} H=-J\sum X_jX_{j+1}-h\sum Z_j \nonumber \end{align} \tag{ 4.53 }$$

can be represented on a line with the (operator valued) matrix

$$\begin{align} M&= \left(\begin{array}{ccc} \mathbb{1}&0&0\nonumber \\ X&0&0\nonumber \\ -hZ&-JX&\mathbb{1} \end{array}\right)\label{eqn:triangle1} \nonumber \end{align} \tag{ 4.54 }$$

and end vectors

$$\begin{align} v_{\rm L}=\left(\begin{array}{ccc} 0&0&1 \end{array}\right) \qquad{\rm{and}}\qquad v_{\rm R}=\left(\begin{array}{c} 1\nonumber \\0\nonumber \\0 \end{array}\right). \nonumber \end{align} \tag{ 4.55 }$$

The Hamiltonian on N sites is then obtained as

$$\begin{align} H=v_{\rm L} M^N v_{\rm R}. \nonumber \end{align} \tag{ 4.56 }$$

The Heisenberg model

$$\begin{align} H=-J_X\sum X_j X_{j+1}-J_Y\sum Y_j Y_{j+1}-J_Z\sum Z_j Z_{j+1}-h\sum Z_j \nonumber \end{align} \tag{ 4.57 }$$

can be obtained in the same fashion with

$$\begin{align} v_{\rm L}=\left(\begin{array}{ccccc} 0&0&0&0&1 \end{array}\right),\qquad M= \left(\begin{array}{ccccc} \mathbb{1}&0&0&0&0\nonumber \\ X&0&0&0&0\nonumber \\ Y&0&0&0&0\nonumber \\ Z&0&0&0&0\nonumber \\ -hZ&-J_X X&-J_Y Y&-J_Z Z&\mathbb{1} \end{array}\right),\qquad v_{\rm R}=\left(\begin{array}{c} 1\nonumber \\0\nonumber \\0\nonumber \\0\nonumber \\0 \end{array}\right).\label{eqn:triangle2} \nonumber \end{align} \tag{ 4.58 }$$

More generally, an MPO can be used to represent any operator which does not increase the Schmidt rank of any state too much. An existing explicit analytic construction of MPOs for 1D local Hamiltonians, as well as a new generalisation for higher dimensional Hamiltonians, is covered in more detail in appendix A.

Problems 3. Solutions in supplementary material (stacks.iop.org/JPhysA/50/223001/mmedia).

• 1.  
  Describe the state given by an MPS with tensor

  

  where index ordering is as shown and indices 1 and 2 are combined. Boundary conditions require inserting a Pauli Z before closing periodic BCs, similar to equation (4.19).
• 2.  
  Describe the state given by the MPS whose only nonzero components are

  

  where the left and right boundary conditions are $\vert {0} \rangle$.Hint: Writing out the matrices corresponding to fixing the physical index might help!
• 3.  
  Describe the qudit state given by the MPS

  

  where $i, j\in \mathbb{Z}_d$, $\oplus$ denotes addition mod d, the left boundary condition is $\vert {0} \rangle$, and the right boundary is $\vert {q} \rangle$ for some $q\in\mathbb{Z}_d$.
• 4.  
  Let $\mathcal{G}$ be some group. Describe the operator given by the MPO with

  

  where the left boundary condition is $\vert {1} \rangle$, the right boundary is $\vert {q} \rangle$ for some $q\in\mathcal{G}$, and $g\cdot h$ denotes group multiplication.
• 5.  
  Suppose the local basis is labelled by particle number. What is the action of the following operator (bond dimension linearly increasing left to right)?

  

  with left vector $L=\vert {0} \rangle$ and right vector $R=\sum_{i=0}^N i\vert {i} \rangle$.
• 6.  
  Write an MPO for the transverse-field-cluster Hamiltonian

  $$\begin{align} H=-J\sum_{j} Z_{j-1}X_jZ_{j+1}-h\sum_j X_j. \nonumber \end{align} \tag{ 4.64 }$$

  Hint: This can be done with bond dimension 4.
• 7.  
  Use the ideas of MPSs and MPOs to prove that $\log$ depth quantum circuits can be simulated efficiently on a classical computer.

The classical definition of a phase, or more particularly a phase transition, is usually associated to some nonanalytic behaviour of the free energy density

$$\begin{align} f(\beta,{\bf v})=-\frac{\log{\rm tr}~{\rm e}^{-\beta H({\bf v})}}{\beta}, \nonumber \end{align} \tag{ 5.1 }$$

where ${\bf v}$ is some vector of parameters of the model (pressures, masses, coupling strengths, etc) and H the Hamiltonian of our system. Clearly when we take the quantum limit ($\beta\to\infty$), the free energy is simply the ground state energy. A quantum phase transition is thus associated with the ground state [29].

At a classical phase transition, correlations become long ranged

$$\begin{align} \langle \mathcal{O}_0\mathcal{O}_x\rangle-\langle\mathcal{O}_0\rangle\langle\mathcal{O}_x\rangle\sim\vert x\vert ^{-\nu}, \nonumber \end{align} \tag{ 5.2 }$$

where the averages are taken with respect to some thermal distribution. We therefore say that a thermal (classical) phase transition is driven by thermal fluctuations, where the variance measures the increasingly long range of these fluctuations. A quantum phase transition also has divergent correlation length, however there is no thermal average— the statistics are purely quantum in origin [29].

A classical phase corresponds to a range of deformations of H and β which can be made without causing nonanalyticities in the free energy f. Likewise, a quantum phase transition occurs where the ground state energy becomes nonanalytic (in the thermodynamic limit) as a function of some Hamiltonian parameters (not temperature this time!). Suppose we have a continuous family of quantum Hamiltonians $H(\lambda)$. The lowest energy levels generically act in one of the following ways [29]:

On the left, there is no phase transition, whilst on the right a transition occurs when the roles of the ground and first excited states cross.

For our purposes, a phase transition will be associated with a gapless point in the spectrum. Therefore, we will say that two states $\vert {\psi_0} \rangle$ and $\vert {\psi_1} \rangle$ are in the same phase if there is a continuous family of Hamiltonians $H(\lambda)$ such that $\vert {\psi_0} \rangle$ is the ground state of H(0), $\vert {\psi_1} \rangle$ is the ground state of H(1), and the gap remains open for all $\lambda\in[0, 1]$.

An equivalent notion is finite time evolution under a local Hamiltonian [30]. Two states are in the same phase if they can be interconverted by time evolution for a finite period. This is linked to the possibility of one state naturally evolving into the other.

It is simpler, and essentially equivalent, to ask which states can be interconverted by a local quantum circuit of depth constant in the system size [20, 31]. We will work within this framework. One may also ask the more complicated question of how phases change if we impose a symmetry; if we insist that all of the Hamiltonians $H(\lambda)$ commute with some symmetry group $U_g(\lambda)$. In the circuit picture, this corresponds to restricting the gate set to only gates which commute with this symmetry [31–33].

In this lecture, we will restrict ourselves to the case of injective MPS [18, 34]. If we assume the MPS is in left canonical form

then injective MPS are those for which the identity is the unique eigenvalue 1 left eigenvector of the transfer matrix. Moreover this means that there exists a unique full-rank⁸ density matrix ρ which is a 1 right eigenvector, i.e.

$$\begin{align} \sum_{j=0}^{d-1}A_j\rho A^\dagger_j=:\mathcal{E}(\rho)=\rho\label{eqn:MPSchannel} \nonumber \end{align} \tag{ 5.4 }$$

These MPS correspond to unique gapped ground states of local Hamiltonians [35]. The arguments we will present here generalise to non-injective MPS, however they become very technical.

We will refer to states which cannot be connected by any constant depth local circuit as being in distinct topological phases, or having distinct topological order. This is to distinguish them from the symmetric phases we will discuss later in this lecture. In fact, we will see that there are no nontrivial topological phases in 1D [20].

Let A_j define some injective MPS, and construct the transfer matrix $\mathbb{E}$⁹

As discussed in the previous lecture, this can be used to renormalise the MPS. Taking products of this transfer matrix corresponds to blocking sites of the original MPS. Since the MPS is injective, the leading eigenvalue of $\mathbb{E}$ is 1 and all other eigenvalues are strictly smaller. Therefore, by taking the kth power of the transfer matrix, we obtain a new transfer matrix which is

where $\vert \lambda_2\vert <1$ is the second eigenvalue of the transfer matrix and ρ is the fixed point of the channel. This transfer matrix can be decomposed to give a new effective MPS tensor describing the long wavelength physics

On the regions we blocked together, we could have first applied a unitary to the state without changing the blocked transfer matrix. Since we only required a constant number of sites to be blocked to achieve this MPS tensor, this unitary freedom is restricted to a constant depth unitary circuit—precisely the equivalence we wish to allow.

Now, let V be some unitary which acts as $\sum_{j, k} \sqrt{\rho}_{j, k}\vert{\,j, k} \rangle \to\vert {0, 0} \rangle$ on the state given by $\sqrt{\rho}$ and arbitrarily on the rest of the space. We can now use this to apply two circuit layers to the MPS

which completely disentangles the MPS, giving the state $\vert {00\cdots0} \rangle$.

Notice that this was all achieved by simply blocking a constant number of sites together, so we have only used a constant depth quantum circuit. Therefore, all injective MPS are in the same (topological) phase as the product state, and therefore each other.

The proofs in this section are translated into TNN from [34].

Since there are no nontrivial topological phases, we will now examine what happens when a symmetry restriction is imposed on the allowed gates. Let $\mathcal{G}$ be some symmetry group for a state which acts on-site as $U_g:=u_g^{\otimes n}$ for each $g\in\mathcal{G}$, where u_g is a unitary representation of $\mathcal{G}$ acting on a single site. Recall that for u_g to be a representation, we must have

$$\begin{align} u_gu_h=u_{gh} \nonumber \end{align} \tag{ 5.10 }$$

for all $g, h\in\mathcal{G}$.

Let A be an MPS tensor such that $\vert {\psi[A]} \rangle$ is symmetric, meaning that $U_g \vert {\psi[A]} \rangle ={\rm e}^{{\rm i}\phi_g}\vert {\psi[A]} \rangle$ for all $g\in\mathcal{G}$. We will now examine how this symmetry is realised on the MPS tensor itself.

We require an understanding of the action of unitaries on the physical level of an MPS, and when they can be 'pushed through' to act on the virtual level. There, they won't be touched by the action of constant depth symmetric circuits on the physical legs, so any properties associated with the virtual action of the symmetry will be an invariant of the phase.

We require two lemmas.

Lemma 1. Let u be some unitary and A an injective MPS tensor. Then the largest eigenvalue λ of the u-transfer matrix

is contained within the unit disk.

Proof. Let $v^\dagger$ (note that we are not assuming that this is unitary) be a left eigenvector of $\mathbb{E}_u$

We therefore get for some density matrix ρ

Once again let ρ be the (unique) right eigenvector of $\mathbb{E}$ with eigenvalue 1. We can view the above expression as an inner product between two vectors

We can now apply the Cauchy-Schwarz inequality across the dotted line, giving

where the vertical lines indicate absolute value. Thus we have

and so $\vert \lambda\vert \leqslant1$. □

Lemma 2. Equality is achieved in lemma 1 if and only if there exists a unitary v and an angle θ such that

Proof. First we prove the 'if' direction. Assume that equation (5.19) holds. Then

and so we have found a left eigenvector $v^\dagger$ of $\mathbb{E}_u$ with a modulus 1 eigenvalue of $\lambda={\rm e}^{{\rm i}\theta}$.

Now we prove the 'only if' direction. Assume there exists a left eigenvector $v^\dagger$ with eigenvalue of modulus 1, then the Cauchy-Schwarz inequality equation (5.15) must become an equality. Therefore, there is some scalar α such that

Taking the norm of each side as vectors, we have

Therefore, $\vert \alpha\vert =1$, so $\alpha={\rm e}^{{\rm i}\theta}$.

Since ρ is full rank, it is invertible, so

Now, rearranging this and left multiplying by $v^\dagger$, we have

We therefore see that $v^\dagger v$ is a left eigenvector of the transfer matrix $\mathbb{E}$ with norm-1 eigenvalue. By assuming injectivity however we require that the only norm-1 eigenvalue is the non-degenerate  +1 eigenvalue, whose left eigenvector is the identity. Thus we conclude v is, after rescaling, unitary, and that equation (5.19) therefore holds. □

So far, we have established that a unitary u can be 'pushed through' the MPS tensor if and only if the u-transfer matrix has an eigenvalue of unit magnitude. We will now show that u is a local symmetry if and only if it can be pushed through. This will complete our understanding of the action of local symmetries on MPS tensors.

Theorem 1 (Symmetries push through). Let $\mathcal{G}$ be a group. A unitary representation u_g is a local symmetry if and only if

for v_g unitary and $\theta_g\in[0, 2\pi)$.

Proof. If equation (5.31) holds, it is clear that u_g is a symmetry since v_g is simply a gauge transformation on the MPS.

Let

be the reduced density matrix on k sites, where ρ is the right fixed point of $\mathbb{E}$. By construction, ${\rm tr}(\sigma_k)=1$, but $\sigma_k$ will generically be mixed, so ${\rm tr}(\sigma_k^2)\leqslant1$. Recall that the purity of a density matrix is lower bounded by the inverse of the matrix-rank, i.e. ${\rm tr}(\sigma_k^2)\geqslant 1/{\rm rank}(\sigma_k)$. Since our reduced density matrix is obtained from a bond dimension D MPS, it has rank at most D². Therefore

where the second equality holds because u_g is a local symmetry.

Here, the left and right boundary vectors ($\mathbb{1}$ and ρ) are independent of the number of sites upon which $\sigma_k$ is supported, so this inequality holds for all k. This can only be the case if $\mathbb{E}_{u_g}$ has an eigenvalue of magnitude 1, as it would otherwise have to possess exponential decay. From lemma 2, this implies that there exists some unitary v_g and an angle $\theta_g$ such that

which completes the proof. □

We now investigate the properties of the virtual action of the symmetry. As discussed above, if we apply a constant depth circuit with symmetric gates to the MPS (i.e. mapping us to any other state in the phase), we can push the symmetry action first through the circuit and then onto the virtual level. Therefore, any properties it has will be an invariant of the phase.

Aside 2 (Projective representations). Let $\mathcal{G}$ be some group. A (linear) representation u_g obeys

$$\begin{align} u_gu_h=u_{gh}\hspace{5mm}\forall g,h\in \mathcal{G}.\label{eqn:rep} \nonumber \end{align} \tag{ 5.37 }$$

This is not the most general way of acting with a group however. We could also ask for

$$\begin{align} v_gv_h=\omega[g,h]v_{gh}\hspace{5mm}\forall g,h\in \mathcal{G}, \nonumber \end{align} \tag{ 5.38 }$$

where $\omega[g, h]={\rm e}^{{\rm i}\phi[g, h]}$ is a scalar which depends on both g and h independently. This is known as a projective representation. One might ask whether this is simply a more complicated way of writing a linear representation. Maybe we can rephase v_g to obtain equation (5.37). Let $\beta[g]$ be some phase depending only on g then after a rephasing $v_g\mapsto \beta[g]v_g$, we have

$$\begin{align} v_gv_h=\omega[g,h]\frac{\beta[gh]}{\beta[g]\beta[h]}v_{gh}=\omega^\prime[g,h]v_{gh}. \nonumber \end{align} \tag{ 5.39 }$$

We say that ω and $\omega^\prime$ are equivalent if they are related in this way, so

$$\begin{align} \omega\sim\omega^\prime\iff \exists\beta : \omega^\prime[g,h]=\frac{\beta[gh]}{\beta[g]\beta[h]}\omega[g,h]. \nonumber \end{align} \tag{ 5.40 }$$

A projective representation is therefore equivalent to a linear representation if the phases can be completely removed, i.e. there exists a β such that

$$\begin{align} \omega[g,h]= \frac{\beta[g]\beta[h]}{\beta[gh]}. \nonumber \end{align} \tag{ 5.41 }$$

As you will show in problems 4, there are projective representations which are not equivalent to any linear representation.

Suppose we act with u_g followed by u_h on the MPS tensor, then

We could also have combined $u_gu_h=u_{gh}$ before pushing through, which tells us

Therefore

$$\begin{align} (v_g\otimes v_g^\dagger)(v_h\otimes v_h^\dagger)=\frac{{\rm e}^{{\rm i}\theta_{gh}}}{{\rm e}^{{\rm i}\theta_g}{\rm e}^{{\rm i}\theta_h}} v_{gh}\otimes v_{gh}^\dagger, \nonumber \end{align} \tag{ 5.44 }$$

so $(v_g\otimes v_g^\dagger)$ is equivalent to a linear representation. We can split this across the tensor product, telling us that in general

$$\begin{align} v_gv_h=\omega[g,h]v_{gh}, \nonumber \end{align} \tag{ 5.45 }$$

where ω is some phase. We cannot say anything about the phase in this case, since anything would be cancelled by tensoring with the conjugate.

The only freedom we have to change v_g within a phase is local rephasing, therefore the equivalence classes of ω label the different phases of injective MPS with a symmetry restriction. These equivalence classes are indexed by the so-called second group cohomology class of the group $\mathcal{G}$, an object usually written as $\mathcal{H}^2(\mathcal{G}, U(1))$ [24, 30].

Problems 4. Solutions in supplementary material (stacks.iop.org/JPhysA/50/223001/mmedia).

• 1.  
  The group $\mathbb{Z}_2\times \mathbb{Z}_2$ has the presentation $\mathbb{Z}_2\times \mathbb{Z}_2=\langle x, z\vert x^2=z^2=e, xz=zx\rangle$. Show that the Pauli matrices form a projective representation of $\mathbb{Z}_2\times \mathbb{Z}_2$.Hint: let $v_x=X, \, v_z=Z, \, v_{xz=zx}=Y$ and show that $v_gv_h=\omega[g, h]v_{gh}$, where ω is some phase.
• 2.  
  Determine the factor system $\omega[g, h]$ for the Pauli matrices.
• 3.  
  Show that the Pauli projective representation is not equivalent to a linear representation.Hint: xz  =  zx, can we rephase v_x and v_z to make $v_xv_z-v_zv_x=0$?
• 4.  
  Recall from section 4.2 that the symmetry of the cluster state is $\mathbb{Z}_2\times \mathbb{Z}_2$, with the action on the MPS tensor being

  

  What can we conclude about the cluster state?

By far the most studied and successful of the algorithms in the field is DMRG. For clarity we will be restricting ourselves to finite DMRG, though there do exist thermodynamic variants. DMRG is an umbrella term which encompasses several similar algorithms, the algorithm we will discuss here is a simplified but nonetheless effective example. As the introduction of this algorithm in [53] pre-dates TNN, its description has historically been presented in a far more physically motivated and technically complicated manner. Due to the corresponding shift in interpretation, the original acronym now holds little relevance to the modern tensor network interpretation of DMRG, and so for clarity we intentionally omit defining precisely the expansion of DMRG as an acronym. For a full review in pre-TNN notation see [54], and see [21] for a TNN treatment.

Representing the Hamiltonian by an MPO, optimising the Rayleigh quotient over MPS looks like the following:

The difficulty is that as we need the contraction of these MPS tensors; the overall objective function is highly non-linear, but it does however only depend quadratically on each individual tensor. The key heuristic behind DMRG is to exploit the simplicity of these local problems, approximating the multivariate (multi-tensor) optimisation by iterated univariate (single tensor) optimisations.

Note that while the DMRG algorithm we are going to outline only calculates ground states, related generalisations exist which can be used to simulate excited states, dynamics etc.

6.1.1. One-site.

The simplest interpretation of the above sketch of DMRG is known as DMRG1 (or one-site DMRG). For a fixed site i, the sub-step involves fixing all but a single MPS tensor, which is in turn optimised over, i.e.

$$\begin{align} A_i\longleftarrow\mathop{{\rm arg\,min}}_{A_i}\frac{\left\langle{\psi(A_i)}\vert{H}\vert{\psi(A_i)}\right\rangle}{\langle{\psi(A_i)}\vert{\psi(A_i)}\rangle}. \nonumber \end{align} \tag{ 6.3 }$$

In TNN these step look like:

Next we define the environment tensors

which correspond to taking closed tensor networks—the expectation values of H and the I respectively—and removing the objective tensor. Given these environments, the sub-step in equation (6.4) becomes

Vectorising this equation yields

$$\begin{align} A_i\longleftarrow\mathop{{\rm arg\,min}}_{A_i}\frac{\left\langle{A_i}\vert{\mathcal{H}_i}\vert{A_i}\right\rangle}{\left\langle{A_i}\vert{\mathcal{I}_i}\vert{A_i}\right\rangle}. \nonumber \end{align} \tag{ 6.8 }$$

Finally we can simplify the denominator of this objective function by appropriately gauge-fixing our MPS to be in canonical form. By putting the parts of the MPS left of our site in left-canonical form, and those to the right in right-canonical form, then we get that $\mathcal{I}_i$ simply reduces to the identity:

Given this canonicalisation, the problem thus reduces to

$$\begin{align} A_i\longleftarrow\mathop{{\rm arg\,min}}_{A_i}\frac{\left\langle{A_i}\vert{\mathcal{H}_i}\vert{A_i}\right\rangle}{\langle{A_i}\vert{A_i}\rangle}. \nonumber \end{align} \tag{ 6.10 }$$

As $\mathcal{H}_i$ is Hermitian, this optimisation has a closed form solution given by the minimum eigenvector¹¹ of $\mathcal{H}_i$. By sweeping back and forth along the chain, solving this localised eigenvector problem, and then shifting along the canonicalisation as necessary, we complete our description of the algorithm.

The main advantage of DMRG1 is that the state stays within the MPS manifold without the bond dimension growing, meaning that the algorithm is greedy¹². This strict restriction on the bond dimension can however be a double-edged sword; this means that there is no particularly convenient method of gently growing the bond dimension as the algorithm runs¹³, and no information is gained regarding the appropriateness of the choice of bond dimension. Both of these problems are addressed in turn by the improved, albeit slightly more complicated, DMRG2 algorithm.

6.1.2. Two-site.

The idea with DMRG2 is to block two sites together, perform an optimisation in the vein DMRG1, then split the sites back out. This splitting process gives DMRG2 its power, allowing for dynamic control of the bond dimension, as well as providing information about the amount of error caused by trimming, which helps to inform the choice of bond-dimension.

First an optimisation is performed:

which can once again be solved by taking the minimum eigenvector of an environment tensor with respect to two sites, $\mathcal{H}_{i, i+1}$, once again in mixed canonical form. After this the two-site tensor is split apart by performing an SVD¹⁴ and a bond trimming:

This trimmed SVD has two key features. Firstly the bond dimension to which we trim could be higher than that we originally started with, allowing us to gently expand out into the space of higher bond dimension MPS. Secondly we can use the truncated singular values to quantify the error associated with this projection back down into the lower bond dimension space, better informing our choice of bond dimension.

Time-evolving block decimation (TEBD) [58, 59] is a tensor network algorithm that allows the dynamics of 1D spin systems to be simulated. By simulating imaginary-time-evolution low-temperature features such as the ground state may be calculated as well.

To simulate imaginary-time-evolution, we need to approximate the imaginary-time-evolution operator $U(\tau)=\exp(-\tau H)$. The problem here is that whilst we may have an efficient representation of H, any exponential of it will not necessarily have a succinct representation. Take the example of a two-body Hamiltonian with corresponding imaginary-time-evolution operator

$$\begin{align*}U(\tau)={\rm e}^{-\tau \sum_ih_i}\qquad{\rm{where}}\qquad H=\sum_i h_{i}\nonumber \end{align*}$$

and h_i is an interaction term acting on spins i and i  +  1. Whilst H has a constant Schmidt rank, admitting an efficient representation as an MPO, $U(\tau)$ generically has exponential bond dimension for almost all τ.

Let $H_o(H_e)$ denote the sum of terms h_i for odd(even) i. As all the terms within $H_o(H_e)$ are commuting, ${\rm e}^{-\tau H_o}({\rm e}^{-\tau {\rm i}t H_e})$ can be efficiently computed and represented. The problem of approximating $U(\tau)$ can therefore be reduced to the problem of approximating ${\rm e}^{-t(A+B)}$ when only terms of the form ${\rm e}^{-\tau A}$ and ${\rm e}^{-\tau t B}$ can be computed.

The central mathematical tool to TEBD are the exponential product approximations. The first order of these approximation is the Suzuki-Trotter formula, which approximates the total evolution by simply evolving each subsystem:

$$\begin{align*}{\rm e}^{-\tau (A+B)}={\rm e}^{-\tau A} {\rm e}^{-\tau B}+\mathcal{O}(\tau^2).\nonumber \end{align*}$$

It turns out there exist entire families of such approximations [60], though for our purposes we will just illustrate the procedure for Suzuki-Trotter.

The TEBD algorithm works by approximating the imaginary-time-evolution operator by the above exponential product formulae, applying it to a given MPS, and trimming the bond dimension to project back down into the space of MPS.

Our approximation to the imaginary-time-evolution operator is given by a product of layers containing only nearest-neighbour two-site operators, meaning we need only be able to contract these operators into our MPS. Suppose we want to apply an operator U to the spins at sites i and i  +  1. The idea is to apply the operator, contract everything into a single tensor, then once again use an SVD trimming to truncate the bond dimension back down.

The benefits this trimming procedure gave to DMRG2—namely control over bond dimension growth and quantification of trimming errors—are also seen in TEBD. As the above procedure is entirely localised, TEBD also admits a large amount of parallelisation, not typically available to DMRG.

From-scratch implementation of these simple algorithms can be achieved with relative ease, however several high performance libraries exist for research level simulations. We direct the interested reader to investigate ITensor [61] (C+  +), evoMPS [62] (Python), Matrix Product Toolkit [63] (C+  +), uni10 (C+  +) [64], Tensor Operations [65] (Julia) among others. A simple tensor class can also be easily written in MATLAB.

Problems 5. Solutions in supplementary material (stacks.iop.org/JPhysA/50/223001/mmedia).

• 1.  
  Consider the critical transverse Ising model

  $$\begin{align} H=-\sum_{i=1}^{n-1}X_iX_{i+1}-\sum_{i=1}^{n}Z_i. \nonumber \end{align} \tag{ 6.13 }$$

  For open boundary conditions, it is known that the ground state energy as a function of n has the form [66]

  $$\begin{align} E(n)=1-\csc\left(\frac{\pi}{\alpha n+\beta}\right) \nonumber \end{align} \tag{ 6.14 }$$

  for some integers α and β. Using either DMRG or TEBD, estimate the ground state energy for several chain lengths and calculate α and β.
• 2.  
  It is known that the local hamiltonian problem is in $\mathsf{P}$ for 1D gapped Hamiltonians [40, 42–45]. DMRG and TEBD are the most common techniques for numerically finding the ground states of such systems. For a gapped and 1D local Hamiltonian, prove that DMRG or TEBD converge.

We have already seen the PEPS construction in 1D. Let $\vert {\phi} \rangle \in \mathbb{C}^D\otimes \mathbb{C}^D$ be some (usually) entangled pair and $\mathcal{P}:\mathbb{C}^D\otimes\mathbb{C}^D\to\mathbb{C}^d$ some linear map. Then

where

is the chosen entangled pair. As we saw, we have a large choice in the exact description we use. We can transform the local basis of each spin in the entangled pair by any (left) invertible matrix

$$\begin{align} \vert {\phi} \rangle \to (A\otimes B)\vert {\phi} \rangle, \nonumber \end{align} \tag{ 7.3 }$$

since we can modify $\mathcal{P}$ to compensate

$$\begin{align} \mathcal{P}\to \mathcal{P} (B^{-1}\otimes A^{-1}).\label{eqn:gauge} \nonumber \end{align} \tag{ 7.4 }$$

One thing to note is that $\vert {\phi} \rangle$ does not necessarily need to be a valid quantum state. We usually leave it unnormalised for convenience.

In addition to this gauge freedom, we have additional choices in the description. We could use entangled triplets for example. Let $\vert {\psi} \rangle =\vert {000} \rangle +\vert {111} \rangle$, then we could choose our PEPS to be

Clearly this doesn't offer any more descriptive power than using entangled pairs. Suppose we have some PEPS projector $\mathcal{Q}$ acting on pairs, then we can extend this to a $\mathcal{P}$ acting on triplets by

$$\begin{align} \mathcal{P}=\mathcal{Q}\left(\mathbb{1}\otimes(\vert {0} \rangle \langle{00}\vert+\vert {1} \rangle \langle{11}\vert)\right). \nonumber \end{align} \tag{ 7.6 }$$

In the other direction, we can build a product of triplets using a minor modification of the GHZ MPS presented above and then use $\mathcal{Q}$ to build our state of interest.

The extension from one to higher dimensional systems proceeds straightforwardly. We will discuss the simple case of a hypercubic lattice, but the framework can be carried out on any graph. In particular, we will restrict to 2D.

As before, we allow $\vert {\phi} \rangle$ to be some entangled pair. The PEPS is built as the natural generalisation to 2D

where

$$\begin{align} \mathcal{P}:\left(\mathbb{C}^D\right)^{\otimes 4}\to\mathbb{C}^d \nonumber \end{align} \tag{ 7.8 }$$

is some linear operator from the virtual to the physical space.

Clearly there is a large amount of gauge freedom in this description as there was in the 1D case. Any invertible transformation of each virtual spin can be compensated in the definition of the PEPS 'projector' $\mathcal{P}$, analogous to equation (7.4).

As in the MPS, one may ask whether using different entanglement structures leads to greater descriptive power. It is easy to see that this is not the case in general. Suppose we choose to lay down plaquettes in a GHZ state and then act with PEPS projectors between plaquettes.

We can use a standard PEPS to prepare this resource state, so any state which can be prepared from this 'projected entangled plaquette' construction can be prepared from a PEPS at small additional cost.

We will now look at several example PEPs.

7.3.1. Product state.

We have already seen this example in 1D. Exactly the same thing works in 2D, for example take

$$\begin{align} \mathcal{P}=\vert {0} \rangle \left\langle\begin{array}{ccc} &0&\nonumber \\ 0&&0\nonumber \\ &0& \end{array}\right\vert. \nonumber \end{align} \tag{ 7.10 }$$

7.3.2. GHZ state.

Directly generalising the 1D case, we can use

$$\begin{align} \mathcal{P}=\vert {0} \rangle \left\langle\begin{array}{ccc} &0&\nonumber \\ 0&&0\nonumber \\ &0& \end{array}\right\vert +\vert {1} \rangle \left\langle\begin{array}{ccc} &1&\nonumber \\ 1&&1\nonumber \\ &1& \end{array}\right\vert \nonumber \end{align} \tag{ 7.11 }$$

to build the GHZ state.

7.3.3. RVB state.

Let D  =  3 be the bond dimension and let

for $\alpha\in\{1, 2\}$, as well as all rotations on the virtual level, be the only nonzero elements of the PEPS tensor. Suppose we tile these tensors and project the dangling indices onto the $\vert {2} \rangle$ state. What is the resulting physical state?

This state is known as the resonating valence bond state [22, 67, 68] and consists of a superposition of all complete tilings of the lattice with maximally entangled pairs

where

Aside 3 (Kitaev's Toric code). Kitaev's Toric code [69] is a canonical example of a topologically ordered model Here we will construct a Hamiltonian with the code space as the ground space of the model. The ground state of this Hamiltonian is the superposition of all closed loops of flipped spins.

We place qubits on the edges of a square lattice.

We wish to create a Hamiltonian with closed loop states (of flipped spins) as the ground state. Suppose all spins are initially in the $\vert {0} \rangle$ state. Then around every vertex v place an interaction

To be in the ground state of this term, the number of edges flipped to $\vert {1} \rangle$ neighbouring a given vertex must be even. Drawing edges carrying flipped spins in red, we can trace the effect of this on the lattice

We can see that on a square graph, requiring an even number of edges incident on each vertex enforces that all of our loops are closed.

At this point, our ground space contains all states with only closed loops. We want an equal superposition of all closed loop states. This is achieved by placing an interaction around plaquettes or squares on the lattice, which convert between loop states. To be an eigenstate, all loop states reachable from the vacuum state must be in the superposition. At each plaquette p, place an interaction

This has the desired effect. Placing the interaction at the indicated plaquette performs the following transformation of loops

It's not hard to convince yourself that all loop states can be reached from the empty state, so all closed loop patterns must be in the superposition. The final Hamiltonian is

and the ground state is an equal superposition over all closed loop states:

Note that the Toric code Hamiltonian is usually presented in the $\vert {+} \rangle /\vert {-} \rangle$ basis rather than the $\vert {0} \rangle /\vert {1} \rangle$ basis.

7.3.4. Toric code ground state.

The simplest way to construct a PEPS for the toric code uses the structure of the ground state. The PEPS tensor is constructed to ensure the superposition of closed loop patterns is achieved upon contraction. The most natural way to achieve this it to write a single tensor for every second plaquette rather than each site.

We begin by adding new edges to the lattice. These edges will become the bonds in the tensor network.

where the plaquettes are numbered for clarity.

Recall that the ground state is built using loops of $\vert {1} \rangle$ in a background of $\vert {0} \rangle$. We choose the state of the added edges such that the loop pattern is preserved

where indicates a spin in the $\vert {1} \rangle$ state on that edge. We choose the following convention when it is ambiguous

which makes everything consistent.

Interpreting these added edges as bonds in a tensor network, we obtain a PEPS tensor for every second plaquette in the original lattice with four physical indices. The nonzero components are

where $i, j, k, l\in\mathbb{Z}_2$. In this tensor the straight legs indicate virtual indices, and the wavy legs physical indices, specifically the four qubits on the given plaquette. The network looks as below, with the dotted lines representing the original lattice:

This tensor simply ensures that if adjacent physical indices are in the $\vert {1} \rangle$ state, i.e. carrying a loop, then the virtual index between them does not carry a loop which would leave the plaquette. Conversely, if only one is in the $\vert {1} \rangle$ state, the loop must leave the plaquette.

Since an even number of the virtual bonds must be in the $\vert {1} \rangle$ state for the tensor entry to be nonzero, the PEPS tensor has a property called $\mathcal{G}$-injectivity [35]. This means that there is a symmetry on the virtual level

This turns out to be closely related to the topological order present in this model.

Let D  =  2 be the bond dimension and let

be the only nonzero elements of the PEPS tensor. The physical state generated is the 2D cluster state, a universal resource for measurement based quantum computing [70, 71].

If we could efficiently take the inner product between PEPS (i.e. contract a square grid network), then we can clearly classically simulate single qubit post selected measurements by simply contracting rank 1 projectors onto the physical indices of these PEPS tensors. This shows us that we cannot contract even simple PEPS states efficiently, unless post-selected quantum computing can be classically simulated ($\mathsf{Post{\rm{-}}BQP}=\mathsf{BPP}$) [72].

7.4.1. Numerical PEPS.

Above, we saw a number of properties of 1D PEPS or MPS. We will now see which properties hold in two dimensions. One might naïvely expect MPS and more general PEPS to share similar properties. As we will see below, these two tensor network states share qualitatively different properties, both in terms of the physics the corresponding states exhibit, and in the computational power of the tensor networks.

Aside 4 (Tensor network for classical partition function). Let $H[s]=\sum_{\langle i, j\rangle }h[s_i, s_j]$ be some classical Hamiltonian. We frequently want to calculate the partition function $\mathcal{Z}=\sum_{\{s\}}{\rm e}^{-\beta H[s]}$ for such a system at a temperature β. We can use a simple tensor network to help.

Define the two tensors

Placing a D tensor at every classical spin and an M tensor corresponding to each interaction, the following network evaluates to the partition function.

Thermal expectation values can be calculated by inserting local tensors into this network. For example

where

has been inserted at site n.

Notice that by combining D and M tensors, the partition function can be described with a single tensor

Let

$$\begin{align} H=-J\sum_{\langle i,j\rangle }s_i s_j, \nonumber \end{align} \tag{ 7.26 }$$

where $s\in\{\pm1\}$, the classical Ising model. The tensor Q then simplifies to

7.5.1. Algebraic decay of correlations.

As we saw above, MPS can only capture states with exponential decay of correlations (or constant correlations of course). We will now see if this holds in the case of PEPS. We can build a PEPS state corresponding to a classical partition function by modifying the above construction [67]. Let

or equivalently combine these into

This defines a PEPS state

Note this is a pure state, and not a thermal state. It is however not normalised, with $\langle{\psi}\vert{\psi}\rangle=Z$. Correlation functions computed using this state are equal to those computed using classical statistical physics. Suppose we were to consider a classical model with a thermal phase transition (such as the Ising model above). Such a model will exhibit algebraic decay of correlations at the critical temperature, implying that the corresponding PEPS does as well. Thus we can see that unlike MPS, the states described by PEPS can exhibit algebraic decay of correlations.

7.5.2. Gauge freedom.

The gauge freedom of a PEPS tensor is a simple generalisation of the MPS freedom. As before, we can block tensors together without changing the global state. In addition, we can perform the following transformation (on a translationally invariant PEPS):

where N and M are invertible matrices.

Recall that in the MPS case, we could use this freedom to bring the tensors into a canonical form. This cannot be done exactly in the case of PEPS, though there do exist numerical methods to bring PEPS into approximate canonical forms [76].

Problems 6. Solutions in supplementary material (stacks.iop.org/JPhysA/50/223001/mmedia).

• 1.  
  What is the PEPS tensor required to build the GHZ state on the honeycomb lattice where spins reside on vertices?
• 2.  
  Which 2 qubit gate is obtained by contracting the following tensors along the horizontal index?

  
• 3.  
  The cluster state can be prepared from the all $\vert {+} \rangle$ state by applying CZ between all adjacent spins. Show that equation (7.20) indeed gives the cluster state.Hint: Consider the decomposition of a gate given in the above problem.
• 4.  
  Investigate how logical operators on the physical spins of the Toric code can be pulled onto the virtual level of the PEPS. Can you see why $\mathcal{G}$-injectivity is so important for topologically ordered PEPS?
• 5.  
  Convince yourself that evaluating expectation values on the PEPS constructed from a classical partition function indeed reproduces the thermal expectation values.

8.1.1. Logarithmic violation of the area law.

One of the key properties realised in the MERA which cannot be realised in MPS is a scaling of entanglement entropy. This is easily seen by bond counting. Recall that if n bonds must be broken to separate a region from the rest of the network, the maximum entanglement entropy that can be supported is $n\log D$, where D is the bond dimension. Recall that in the case of MPS any reduced state on a contiguous region can be removed by cutting n  =  2 bonds.

By inspecting the diagram

it is straightforward to see that to remove a block of N physical indices from the rest of the network, $\mathcal{O}(\log N)$ bonds must be cut. This shows that the maximum entropy scales as $\log N \log D$ [77, 78].

8.1.2. Power law decay of correlations.

Using the constraints on the tensors (equation (8.2)), we can simplify the evaluation of a two point correlator on a MERA state [79].

Note that the length scale behaviour of the correlator is completely determined by the application of a superoperator

where the w tensor can be viewed as a set of Kraus operators

obtained by grouping the indices indicated.

Thus, $\mathcal{S}$ is a completely positive, unital map and all eigenvalues λ of $\mathcal{S}$ are $\vert \lambda\vert \leqslant 1$. We can bring operators separated by N sites together by applying $\mathcal{S}$ $\sim \log N$ times. Considering eigenoperators of the $\mathcal{S}$ superoperator, the correlator acts as

$$\begin{align} \langle A_jB_k\rangle\sim\frac{\langle A_0B_1\rangle}{\vert\,j-k\vert ^{\Delta_A+\Delta_B}},\label{eqn:decay} \nonumber \end{align} \tag{ 8.10 }$$

where $\Delta_{\phi}=-\log_3\lambda_\phi$, $\Delta_\phi\geqslant0$ are known as scaling dimensions, where $\lambda_\phi$ is the corresponding eigenvalue of $\mathcal{S}$. Therefore, a MERA state can support algebraic decay of correlations. Although this discussion required the operators to be placed at special sites, it can be easily generalised.

8.1.3. Efficient manipulation.

As described in section 2.5, a good tensor network ansatz should fulfil two properties. First, it should be efficiently storable. All of the networks we have discussed thus far have this property, as only a small number of coefficients are required to represent these states. The second property is more subtle; one should be able to extract physical data efficiently. Although this works for the 1D MPS network, it fails for 2D PEPS states; the contractions required to calculate expectation values of local operators is incredibly hard.

It turns out the MERA has both of these properties. One can efficiently store the state data, and, thanks to the constraints in equation (8.2), one can efficiently compute local expectation values and correlators. We have already seen how this works. The isometric constraints ensure that local operators on the physical level of the network are mapped to local operators on the higher levels [86]. Therefore, computing expectation values only requires manipulation of a small number of tensors in the causal cone of the operator

where the shaded region indicates the causal cone of the five site operator on the physical level indicated in yellow. Notice that the number of tensors on each subsequent level does not grow. Indeed, after a single layer of tensors, the operator becomes a three site operator, and the range never grows. Thus, we see that the layers of the MERA act to map local operators to local operators.

Much of the discussion above concerned interpretation of the layers of the MERA as Kraus operators, defining a unital CP map on local operators. Evaluating expectation values can be seen as application of many superoperators followed by the inner product with some state on a smaller number of sites

$$\begin{align} \langle{\psi_{0}}\vert O\vert {\psi_{0} \rangle } =\langle{\psi_{k+1}}\vert\mathcal{A}_k\circ\cdots\circ\mathcal{A}_1\circ\mathcal{A}_0(O)\vert {\psi_{k+1}} \rangle, \nonumber \end{align} \tag{ 8.12 }$$

where $\mathcal{A}_j$ is a map from 3^N−j spins to 3^N−j/3 spins. This can be seen as a renormalisation group or scale transformation. The state $\vert {\psi_{j}} \rangle$ is supported on 3^N−j spins, and contains only the physical data necessary to understand the physics on that length scale. As we saw, if O is a local operator, $\mathcal{A}(O)$ is easy to evaluate. This allows us to understand the effective operator as a function of length scale [77, 79, 80].

The thermodynamic or macroscopic observables can be seen as the operators obtained by applying a formally infinite number of MERA layers to the high energy or microscopic observables. Thus, the macroscopic physics, or phase structure, is determined by fixed points of the maps $\mathcal{A}$. Some particularly interesting states are the scale invariant states. If the MERA tensors are all the same after some layer, the state is scale invariant. For these states, we do not expect the physics to change as a function of length or energy scale. The fixed point observables of these states are particularly simple to understand, and distinct scale invariant states characterise the different phases.

Since there is no characteristic length scale set by either the spectral gap or correlation length, gapless ground states are expected to be scale invariant. The MERA therefore allows us to understand the long range physics of these states incredibly efficiently [79, 86]. Another way to achieve a scale invariant state is to have zero correlation length—these states characterise gapped phases.

In the appropriate limit, the low energy physics of the gapless spin chains considered here is described by a conformal field theory (CFT) [87, 88]. The physics of CFTs is thought to be related to gravitational theories in one additional dimension [89–91].

This duality can be observed in the MERA network [92–94]. Imposing the graph metric on the MERA, we find a discretised anti-de Sitter (AdS) metric [92], whilst the edge theory is a 'discretised' CFT. In addition to being a concrete realisation of the holographic principle, the MERA/CFT duality provides avenues towards designing quantum error correcting codes [95].

We note that the AdS/MERA connection remains an open research question. Limits on the ability of MERA states to replicate physics on scales less than the AdS radius have been shown [94]. Additionally, whether the geometry is best understood as anti-de Sitter [92] or de Sitter [93] is currently unclear. Whatever the status, the connection is intriguing. We encourage the interested reader to explore the rapidly expanding literature on the topic [94–103].

8.4.1. Product state.

Let

$$\begin{align} w&= \begin{array}{cc} & 0 1\nonumber \\ \begin{array}{c} 000\nonumber \\ 100\nonumber \\ 010\nonumber \\ 110\nonumber \\ 001\nonumber \\ 101\nonumber \\ 011\nonumber \\ 111\nonumber \\ \end{array} &\left(\begin{array}{cc} 1/2&0\nonumber \\ 1/2&0\nonumber \\ 0&1/2\nonumber \\ 0&1/2\nonumber \\ 1/2&0\nonumber \\ 1/2&0\nonumber \\ 0&1/2\nonumber \\ 0&1/2\nonumber \\ \end{array} \right) \end{array} \nonumber \end{align} \tag{ 8.13 }$$

and $u=\mathbb{1}$.

If we build $\log_3 N$ layers using these tensors, we end up with a state on N sites. The network still has a free index at the top, so we need to define a one-index 'top tensor' T to obtain the final state. Let $T=\vert {+} \rangle$. The state obtained is $\vert {+} \rangle ^{\otimes N}$.

8.4.2. GHZ state.

Let

and $u=\mathbb{1}$. Let the top tensor be $T=\vert {+} \rangle$. The state obtained is $\frac{\vert {0} \rangle ^{\otimes N}+\vert {1} \rangle ^{\otimes N}}{\sqrt{2}}$.

8.4.3. Cluster state.

It is more convenient to define the cluster state on a binary MERA than a ternary. Place two spins at each site and let

where is a controlled-Z gate and H is the Hadamard. If we pick a top tensor $T=\vert {++} \rangle$, we obtain the cluster state on periodic boundary conditions.

8.4.4. Gapless states.

Recently, a family of analytic MERA for the critical point of the transverse field Ising model was proposed [104]. One can also use numerical techniques to obtain a MERA approximation to the ground state of a local Hamiltonian however. Here, we will present some physical data obtained for a model known as the transverse field cluster model [105]. In particular, we will present the ground state energy and the decay exponents ($\Delta_\phi$ in equation (8.10)).

This model is most straightforwardly defined with a pair of spin half particles at each site. The Hamiltonian for this model is

$$\begin{align}H=&-\sum_j\left(X_{j}^{(1)}+X_{j}^{(2)} + Z_{j-1}^{(2)}X_{j}^{(1)}Z_{j}^{(2)} + Z_{j}^{(1)}X_{j}^{(2)}Z_{j+1}^{(1)}\right)\nonumber\\&-\lambda\sum_j\left(X_{j}^{(1)}X_{j}^{(2)}+Z_{j}^{(1)}Y_{j}^{(2)}Y_{j+1}^{(1)}Z_{j+1}^{(2)}\right). \nonumber \end{align} \tag{ 8.16 }$$

This is the cluster state Hamiltonian with transverse fields and an additional interaction with variable strength. The Hamiltonian remains gapless for a range of values of λ, over which the ground state energy varies continuously as seen in figure 1(a). The decay exponents also vary over this range, meaning that the thermodynamic physics or RG fixed point is dependent on λ. These exponents can easily be extracted from an optimised MERA by finding the eigenvalues of the $\mathcal{S}$ superoperator in equation (8.8). The MERA results are shown in figure 1(b).

Zoom In Zoom Out Reset image size

Figures reproduced from [105].

Problems 7. Solutions in supplementary material (stacks.iop.org/JPhysA/50/223001/mmedia).

• 1.  
  Can you find a MERA for the W state?
• 2.  
  What state is given by the MERA with

  

  $u=\mathbb{1}$ and top tensor $T=\frac{1}{\sqrt{2}}(\vert {00} \rangle +\vert {11} \rangle)$?
• 3.  
  The above state is the ground state of the Hamiltonian

  $$\begin{align} H=-\sum_{j=1}^{N/2}(X_{2j}X_{2j+1}+Z_{2j}Z_{2j+1}) \nonumber \end{align} \tag{ 8.18 }$$

  on periodic boundary conditions. Is that clear? Can you find a unitary U_2j−1,2j which transforms this Hamiltonian into

  $$\begin{align} H=-\sum_{j=1}^{N/2}(Z_{2j-1}X_{2j}Z_{2j+1}+Z_{2j}X_{2j+1}Z_{2j+2})? \nonumber \end{align} \tag{ 8.19 }$$
• 4.  
  Act with the above transformation U on the MERA tensor to obtain another MERA tensor. What is this state?
• 5.  
  What is the maximum range of thermodynamic observables in a ternary MERA scheme?
• 6.  
  What does the reduced density matrix on a few sites of the MERA look like? Notice that it corresponds to the top tensor being passed through a CPTP map several times, this is usually called the descending superoperator.
• 7.  
  Do tree tensor networks (i.e. MERA for $u=\mathbb{1}$) have any area law violation on contiguous regions?

In this notation, if we label our indices $I=\lbrace\cdot, 1, \rightarrow\rbrace$, then the transverse Ising model Hamiltonian given in equation (4.54) is given by

where the boundary terms fix the far left and right indices to $\vert {\rightarrow} \rangle$ and $\vert {\cdot} \rangle$ respectively.

One common interpretation of this construction is in terms of finite-state automata, with the index values corresponding to the automaton states, and the non-zero index values to the transition rules. The automaton moves from left to right¹⁵, with the boundary vectors setting the initial state to $\vert {\rightarrow} \rangle$ and final state to $\vert {\cdot} \rangle$. With only these restrictions, the automaton can transition from $\vert {\rightarrow} \rangle$ to $\vert {\cdot} \rangle$ either directly (giving the field term  −hZ), or via 1 (giving the Ising term  −JXX) at any location.

To make the higher dimensional generalisation clear we will slightly modify this finite-state automata language, to that of particles and their decay. We can think of  →  as a right-moving particle, and · as the vacuum. The first two transition rules (A.1) correspond to both the vacuum and particle being stable states, with the remaining transitions (A.2) to valid decay routes of the particle. Thus we can interpret the value of the overall MPO as being a superposition over all decays, with each corresponding to a term in the Hamiltonian.

A.1.1. Heisenberg model.

Suppose we wish to construct a Hamiltonian containing multiple two-body terms, such as the Heisenberg anti-ferromagnet, which contains the terms  −J_XXX, −J_YYY, −J_ZZZ, as well as a field  −hZ. An MPO of this model is given in standard notation in equation (4.58).

Added Hamiltonian terms can be accommodated in this construction by extra decay chains. Take our index set to be $I=\lbrace\cdot, x, y, z, \rightarrow\rbrace$ and our MPO to have terms:

Again equations (A.3) correspond to stable vacuum and particles, and each of the transition rules equations (A.4)–(A.7) to each term in the Hamiltonian.

A.1.2. Cluster model.

The Cluster Hamiltonian contains three body terms of the form ZXZ. Larger terms such as this can be accommodated by longer decay chains. Take an index set $I=\lbrace \cdot, 1, 2, \rightarrow\rbrace$ and include the standard stable vacuum/particle terms as well as

By combining the above two techniques, we can construct arbitrary local Hamiltonians.

In higher dimensions we can use a similar construction. Suppose we want to construct a 2D field Hamiltonian, consisting of a Z at every site. Take our index set to be $I=\lbrace \rightarrow, \cdot\rbrace$. Our typical stable vacuum/particle terms that we will always include now become

For the field Hamiltonian we need only allow for a simple particle decay of

As for the boundary conditions, along the top, right and bottom boundaries we will once again fix the only non-zero indices to be the vacuum $\vert {\cdot} \rangle$. Along the left edge, the boundary condition is a virtual W-state (see equation (4.18)) on indices $\lbrace\rightarrow, \cdot\rbrace$, i.e. the equal superposition of all single-particle states. As such we can see that all the non-zero contributions to the Hamiltonian are of the form:

As with 1D, by introducing intermediary states and different decay rules, arbitrary local Hamiltonians in any dimension can be similarly constructed. For example suppose we wanted a 9-body plaquette term of the form:

$$\begin{align*}\begin{array}{ccc} J&K&L\nonumber \\ M&N&O\nonumber \\ P&Q&R \end{array} \nonumber \end{align*}$$

Take $I=\lbrace\cdot, 1, 2, \rightarrow\rbrace$ and our non-trivial decay modes to be

then we can see that the non-zero contributions to the Hamiltonian are of the form

Below are several more example of Hamiltonian constructed by the above method.

A.3.1. Toric code (Wen plaquette).

A.3.2. Quantum compass model/Bacon–Shor code.

A.3.3. 2D transverse Ising.

A.3.4. 2D cluster state.