Decorated tensor network renormalization for lattice gauge theories and spin foam models

Tensor network techniques have proved to be powerful tools that can be employed to explore the large scale dynamics of lattice systems. Nonetheless, the redundancy of degrees of freedom in lattice gauge theories (and related models) poses a challenge for standard tensor network algorithms. We accommodate for such systems by introducing an additional structure decorating the tensor network. This allows to explicitly preserve the gauge symmetry of the system under coarse graining and straightforwardly interpret the fixed point tensors. We propose and test (for models with finite Abelian groups) a coarse graining algorithm for lattice gauge theories based on decorated tensor networks. We also point out that decorated tensor networks are applicable to other models as well, where they provide the advantage to give immediate access to certain expectation values and correlation functions.


INTRODUCTION
Tensor network algorithms [1][2][3][4][5] have become a well used method for the coarse graining of classical partition functions.In particular for models in which Monte Carlo methods are not available (due to sign problems for fermionic systems or due to a priori complex amplitudes as in spin foams1 ) tensor network algorithms may be the only way to explore the phase structure of a given system.We are in particular interested in the phase structure of spin foams, candidate models for quantum gravity.For these models the key open question is the recovery of a dynamics describing continuum manifolds in the limit involving many building blocks.This work will address the question of efficient algorithms for lattice gauge theories, including generalized models such as spin foams.Here the problem is to deal with the gauge symmetry, that is a large redundancy in the description.This issue has been addressed for a framework where tensor networks are used as a variational ansatz for the low energy wave function ( for quantum lattice models) [11,12], however, to our knowledge not for algorithms in which the tensor network represents the partition function.
A tensor network representation for lattice gauge theories and spin foam models appeared in [13,14] and later also in [15].As we will explain later, these representations are not economic in terms of bond dimension of the initial tensor.As the aim is to go to rather 'large' structure group, that is SU (2) k with k large (or as in [15] even SU (2) which would involve infinite bond dimension), this issue is rather important.
The main aim of this work is therefore to develop methods for coarse graining (generalized) lattice gauge theories, that are sufficiently efficient to be applied to higher dimensions (here we will consider 3D, eventually the methods have to be applied to 4D models) and for 'large' structure groups.
To this end we will go beyond tensor networks, allowing more general indexed objects, which we will term "decorated tensor networks".This allows a more efficient representation of lattice gauge theories and spin foam models and furthermore a comparably simple algorithm that emulates the 2D algorithms of [3,4].Although this generalization of tensor networks and the related algorithm were designed for (generalized) lattice gauge theories the framework is also applicable to e.g.Ising like models in two and higher dimension.Here the decorated tensor networks lead to an interesting mixture of keeping part of the variables of the original model, i.e. the Ising spins, and integrating out and transforming the other part of the variables (by field redefinitions), which is the essence of tensor network models.The field redefinitions inherent in tensor network methods turn however the interpretation of the phases and fixed points into a difficult task, see also the discussion in [16].For this reason the algorithm presented here might be also useful for other models than lattice gauge theories.
Spin foams: In the following we will outline the relevance of our aim for the coarse graining of spin foams.Spin foams [6][7][8] are generalized versions of lattice gauge theories, but feature a more complicated structure.Roughly speaking, the models are parametrized by weights for intertwiners.This means that the weights depend on a combination of several representation labels, instead of being parametrized by face weights, which are functions of one representation label only and describe standard lattice gauge models.The phase structure of spin foams, beyond those phases that also appear in lattice gauge theory, is basically unknown, except for the recent developments in [10,17,18].Tensor network algorithms are indeed so far the only tool that is available for the real space coarse graining of these systems (as arXiv:1409.2407v1[gr-qc] 8 Sep 2014 even the Migdal Kadanoff relations [19,20] do not apply to spin foams).Thus a development of these techniques is worthwhile, although there are still some general issues to address [16,21,22].
The full spin foam models are defined in four dimensions and involve a structure group SU (2) × SU (2) (this is for models describing Euclidean signature metrics, which face less technical challenges than those describing Lorentzian metrics, based on SO(3, 1)).However, these models are barely understood beyond a very few building blocks, additionally the issue of divergences arises [23][24][25].[9,13,14] therefore proposed a program in which one considers a reduction of the structure group to a finite group, and a dimensional reduction to two dimensional spin net models.Nevertheless a key input of spin foam construction, the simplicity constraints [6][7][8]10], and with it, the emphasis on the behaviour of intertwiner degrees of freedom [17,18], was kept in these models.In [18] spin net models based on a quantum group SU (2) k are considered which, at least in 3D, describe gravity with a cosmological constant.In this sense the work [18] lifted the reduction of the structure group and revealed a very rich phase space structure for quantum group spin nets, related to phases in anyonic spin chains [26][27][28].Spin net models can be in fact interpreted as spin foams with a large number of (dual) edges and two vertices [18].The different phases describe the effective coupling between the two spin foam vertices, where this coupling is determined by the (initial) choice of simplicity constraints.
It is now crucial to investigate the phases for spin foams in order to test the conjecture that spin nets could indeed possess a phase structure similar to spin foams (this is similar to the relation between lattice gauge theories and Ising like models defined for the same group [29]).To this end it is useful to incorporate techniques that have been developed for spin nets [17,18] in (a) dealing with (global or gauge) symmetries and (b) allowing for a monitoring of the weights for the intertwiner degrees of freedom during the coarse graining.This was accomplished by a symmetry protecting algorithm, in which magnetic indices are pre-contracted, see also [30] and the tensor decomposes in blocks labelled by intertwiner channels.
Outline: The paper is organized as follows: First we review several equivalent prescriptions of (generalized) lattice gauge theories, which can be represented as tensor networks.We introduce a new type of tensor networks, called decorated tensor networks, which for lattice gauge theories preserve explicitly the Gauss constraints, and in this sense the gauge symmetry, of the theory, under coarse graining.We test different versions of decorated tensor network algorithm for the 3D Ising gauge theory, i.e. lattice gauge theory with Z 2 gauge group. 2 We also 2 This theory is dual to the (non-gauge) Ising model [31], which briefly comment on extensions to 4D and non-Abelian gauge theories.Afterwards we present a similar algorithm applicable to Ising-like systems which analogously preserve the original spin variables of the system, such that correlation functions can be straightforwardly calculated.We close with a discussion of the algorithms and results.

REPRESENTATIONS FOR LATTICE GAUGE THEORIES
Lattice gauge models feature a partition function where g e denotes group elements attached to the (oriented) edges e of a lattice, f denotes an (oriented) face or plaquette of the lattice and h f = e⊂f g e the face holonomy, that is the oriented product of group elements associated to the edges belonging to this face.The face weights ω have to be class functions, i.e. ω(ghg −1 ) = ω(h), for the partition function to be invariant under gauge transformations g e → h −1 s(e) g e h t(e) .Here h v denote group elements, acting as gauge parameter, associated to the vertices of the lattice, s(e) and t(e) the source and target vertex of the edge e respectively.
We assumed a finite structure group in (1), for application to Lie groups we have to replace the sum in (1) by an integral.Later we will perform a duality transformation, in which the partition function appears as a finite/infinite sum for finite/ compact Lie groups.Moreover in this representation we can also easily deal with the quantum group case.
Although the representation (1) features variables based on edges, as in tensor networks, the weights are actually associated to the faces and not to the vertices.Additionally the gauge symmetry leads to an overparametrization of the system, which should be avoided.A duality transformation (or strong coupling / high temperature expansion) rewrites the partition function using a group Fourier transform for the weights: with χ ρ denoting the character of the representation, see for instance [9,31].Here the sum is over the irreducible one can see by solving the Gauss constraints.However such a duality transform is not available for (generalized) non-Abelian gauge theories, to which we hope to apply decorated tensor networks eventually.We therefore keep the Gauss constraints explicit.
Here P e are the Haar projectors, that is representation theoretical objects associated to the edges of the lattice.These Haar projectors map onto the subspace associated to the trivial representation in the tensor product of representations ρ f associated to the faces adjacent to the edge e, that is the invariant subspace in the representation space ⊗ f ⊃e V ρ f .Thus we have magnetic indices m e , n e that are contracted in a particular pattern, as shown in figure 1.
Spin foam models implement simplicity constraints by replacing this tensor P e by a projector P e onto a smaller subspace inside the trivial representation space, thus (3) can serve also as the partition function for spin foam models (for further explanation see for instance [9,13,14]).The particular choice of simplicity constraints is then encoded into a particular choice of projector P e .
For Abelian structure groups, for instance Z K , the partition function (3) simplifies to (we will denote the representations for Z K by k ∈ 0, . . ., K − 1) where o(e, f ) = −1 if face and edge orientation disagree and o(e, f ) = +1 if these agree.Thus the Haar projectors reduce to Gauss constraints, imposing that the representations around the faces meeting in one edge sum up to zero.
Although we have now variables associated to faces, this form of the partition function (3) can be brought into the form of a tensor network.To this end one associates tensors (P e ) to the midpoints of each edge and introduces auxiliary tensors to the midpoints of the faces, that make sure, that each tensor (P e ) around a given face f , sees the same representation label ρ f , see also figure 2. This representation appeared in [13,14] and also later in [15].
However, this representation (3) is rather inefficient if it comes to coarse graining algorithms.The dependence of the tensor in the magnetic indices is prescribed by the gauge invariance of the model, and thus one would rather wish to pre-contract these indices, as in fact happens for the symmetry protecting algorithms developed in [17,18] (for models with global symmetry).Thus one would only have to deal with representation and possibly intertwiner labels (which for multiplicity free3 groups are again representation labels), as in [17,18].This would also allow to understand the coarse graining in terms of weights on intertwiners.
Indeed, a standard technique in spin foams, is to split the projector P e -by using a orthonormal basis |ι e of the invariant subspace of ⊗ f ⊃e V ρ f -into two parts (referred FIG. 3. Construction of the vertex amplitude: The projectors Pe on the edges are expanded in the intertwiner basis {ιe} according to (5), the opposing triangles pictorially representing |ιe ιe|.Together with the contraction of magnetic indices shown in figure 1, the projectors can be split and associated to the vertex of the lattice.
to as intertwiners) 4(P e ) see figure 3.These two parts can be associated to the source vertex and target vertex of the edge e respectively.The sum is over intertwiner labels ι e (which are again representation labels assuming multiplicity free groups) appearing in the tensor product ⊗ f ⊃e V ρ f .Now every half edge, labelled by an edge-vertex pair (e, v), carries an object |ι e .The structure of the magnetic index contraction is such, that the sum over the magnetic indices factorizes over the vertices, i.e. the set of objects ι e associated to a given vertex v can be contracted to a so called vertex amplitude A v ((ρ f ) f ⊃v , (ι e ) e⊃v ) that only depends on representation labels associated to the faces and edges adjacent to the vertex v.This is shown in figure 4. Thus we have now a form of the partition function which almost is a tensor network, if it would not be for the representation labels on the faces.(The face weights can be absorbed into the vertex weights for a regular lattice.)This could again be dealt with by introducing auxiliary tensors put into the middle of the faces.One can then FIG. 4. The construction of the vertex amplitude Av for a 3D cubical lattice with the explicit contraction of magnetic indices.FIG. 5. Definition of a homogeneous, yet anisotropic tensor network from the vertex amplitude representation.By absorbing the auxiliary tensors in the positive quadrants into the vertex, one obtains a cubical tensor network, with tensors T , with additional diagonal edges (some of the auxiliary edges can be combined with the original cubical edges, increasing the bond dimension).The regular cubical edges get equipped with additional data, namely the representation on the face.
attempt to group tensors of this inhomogeneous network into a homogeneous one, resulting in tensors with more edges and also higher initial bond dimension (see [15] where this has been applied to the representation (3)).For a cubical 3D lattice a possible grouping is depicted in figure 5 and would result in 12-valent tensors, which are also connected via diagonal edges.
Alternatively 5 , one can actually copy a given face index ρ f to all edges adjacent to this face.To enforce equality between the copied indices we have to multiply the vertex amplitudes (i.e. the tensors) with Kronecker deltas.For a cubical lattice each edge would then carry four rep-resentation labels coming from the adjacent faces and an additional (intertwiner) label ι e , giving a priori 6 × 5 representation indices for a tensor describing a cubical lattice in three dimensions.These labels would be restricted by the condition that the four representations per edge actually couple to the intertwiner ι e .For Abelian groups the intertwiners are trivial, and this condition turns into the Gauss constraint that the four representation labels sum to the trivial representation.This would then still give for a structure group Z K a bond dimension K 3 .This again is a rather inefficient way to proceed.

DECORATED TENSOR NETWORK RENORMALIZATION
Here we propose to generalize the concept of tensor networks in order to keep the gauge constraints satisfied explicitly, and hence allow for more efficient and accurate study of the phase space of a given lattice gauge model.
For the following discussion we found it simpler to transform the representation to the dual lattice (which for spin foams is actually the triangulation).Thus 'tensors based on vertices' are now replaced by 'amplitudes associated to building blocks'.To some degree this is just a change of perspective, see also [32], but it allows for more flexibility in designing (and imagining) algorithms.For a cubical 3D lattice this means that we now associate amplitudes to cubes.Indices, which have been associated to faces, are now identified with the edges ρ f → ρ e , and the intertwiner labels are associated to the faces ι e → ι f , see figure 6.The coupling condition is now also based on the faces, i.e. the representations associated to the edges of a face have to couple to the intertwiner label associated to this face.The coarse-graining / blocking algorithms now proceed by gluing cubes to larger cubes, by summing over the labels associated to bulk edges and bulk faces.

Abelian gauge models in 3D
Let us specialize to Abelian gauge models in 3D.These models are dual to Abelian 'edge' models with global symmetry [31,33,34], as we will see shortly.However, as we want to test algorithms also applicable to non-Abelian gauge models, where such a duality is not available, we will keep the lattice gauge (spin foam) representation of the model.
For Abelian models, for instance with structure group Z K , we denote the representations associated to the (oriented) edges by k e with k e = 0, . . ., K − 1.We have Gauss constraints on the faces, prescribing that the representations of the edges associated to a given face, have to couple to the trivial representation, that is all face FIG. 6. Example of a basic building block of the decorated tensor network.The amplitude tensor Ac is associated to the cube, and contains the information about both edge and face variables, ρe and ι f .Bold edges represent a choice of maximal tree that represents the free variables for Abelian groups.We depict face variables as legs of the central tensor piercing through respective faces.One can think of it as one element of a tensor network decorated with additional edge variables.Later-on we will introduce additional indices associated to edges Ie and/or faces I f .labels are k f = 0.This gives the restrictions where o(e, f ) = −1 if face and edge orientation disagree and o(e, f ) = +1 if these agree.These Gauss constraints can be solved (but as mentioned we are not going to do so) by introducing variables (gauge potentials) k v on the vertices.Setting k e = k s(e) − k t(e) gives a solution to the Gauss constraints (7).As mentioned, this leads back to an edge model in which the gauge group serves now as global symmetry group (e.g.Ising model if we started with Ising gauge theory).Such an edge model could be turned into a tensor network by standard methods.(The easiest way is to use a dual representation -which again will feature Gauss constraints now based on the vertices of the lattice.)Due to the Gauss constraints the edge labels k e are not independent of each other.In fact, for a cube there are five independent constraints of the form (7).An independent set of variables k e can be easily found by choosing a maximal tree (a connected subgraph without loops) among the edges of a cube.Thus by solving the Gauss constraints for the 5 edges which are not part of the tree, the cube amplitude does depend only on 7 variables instead of 12 variables.The initial cube amplitudes are given as where ω are the (Fourier transformed) face weights in (2).Later we will extend this to also contain additional variables, see figure 6.Note that the amplitude can describe not only the cubic models: It is for example possible to embed a tessellation of the cube into 6 tetrahedra into (8), since the total number of independent variables is the same as that of the cube.

Algorithm in the leading order approximation
Let us first convey the spirit of the coarse graining algorithm for the leading order approximation.In this case, throughout the algorithm we will keep the same index structure for the cube amplitudes A c .This amplitude will however flow and loose its initial factorizing form.
The basic algorithm is a (decorated) 3D generalization of the Gu-Wen algorithm [4].We divide each cube into two prisms using a singular value decomposition (SVD).Then four of these prisms are glued to a new cube, see figure 7. (The new building blocks are actually parallelepipeds or rhomboids, but we will refer to them as cubes.)We changed the lattice length, say in the plane (xy), by a factor of √ 2, see figure 7. We then rotate the lattice, so that in the next iteration we coarse grain in an orthogonal plane.
To divide each cube into two prisms, we have to introduce diagonals for two faces (opposite to each other) of the cube.For Abelian models the edges of the subdivided faces determine the representation label for the new edge, as we also impose the Gauss constraints for the new triangular faces.As we introduce two new variables and two additional constraints, the counting of independent variables does not change.
To split the cube we apply singular value decomposition.In general this provides a splitting of a matrix M AB , into two matrices.Cutting off the singular values at the highest one, λ 1 , we obtain a factorizing amplitude We thus would summarize the edge labels of one prism to the index A and the edge labels of the second prism to B. But it seems that we need also to split four of the edges, which are shared by the two prism halves.Again because of the Gauss constraints (now for the new cut faces), these are actually just three independent variables.
We could copy each of the independent labels k e to a pair k e , k e , multiply the cube amplitude with Kronecker delta, δ(k e , k e ) and then proceed with a singular value decomposition into two matrices, where k e and k e will be associated to the two different prisms.FIG. 7. Full iteration of the algorithm.The initial cube is split in two ways into prisms, which are then glued back together into a rhomboid again.The following iteration is performed after a rotation into a different plane.In order to complete the iteration cycle we keep an auxiliary half-edge on two faces, which is summed over in the last step.The bold edges indicate the choice of a maximal tree, the coloured edges highlight the edges shared by the split prisms.For details refer to the text.
However, the Kronecker delta δ(k e , k e ) leads to matrix M AB in block form, where the blocks are labelled by the three independent shared edge variables k e = k e .Thus instead of proceeding with an SVD for the full matrix at once (and just keeping one of the SV's) we perform a singular value decomposition for each of the blocks and keep one SV per block (thus K 3 SV's in total).
That is, we actually perform SVD's of matrices M {ke} A B which are labelled by three independent shared edge labels {k e }.The index A summarizes all the other (initially two independent) edge indices of one prism, and B the other (initially two independent) edge indices of the other prism.
This will result in amplitudes associated to the prisms A p , with four possible prisms resulting from cutting the cube along the two possible face diagonals.The prisms are glued to cubes as depicted in figure 7.These cubes come with a further decoration, namely four half diagonal edges on two opposite sides.Due to Gauss constraints these (eight) additional edges just amount to two further independent variables for the new cubes.These variables can be associated to the faces.
We now rotate the cubes and start the splitting into prisms again.It turns out that we only have to subdivide undecorated faces.During the gluing of the prism the only summation is actually over the decorated faces -thus for gluing a cube we just have to sum over two variables, as the edge shared by all four prisms is Gauss constrained by outer edges.This summation can be done in steps, by gluing first two prisms to half cubes and then just multiplying these to find the new cubes amplitude.Hence the computational complexity scales with K 9 for the last and the most expensive step.It takes negligible time on a laptop for the Ising gauge model in this lowest approximation (less than a second for 50 iterations).
Apart from the simplicity of this algorithm compared to the other options, let us point out another advantage: this is in performing the SVD in block form, where the blocks are labelled by variables, which are straightforward to interpret.In fact, for spin foams the representation labels carry geometric information and represent in 3D indeed the length of the edges.The same strategy applied to the spin net models provided invaluable information about the phase space structure [18] and in general offers an easy way to monitor the outcome of the coarse graining algorithm.For non-Abelian models, the block form SVD is in parallel to the block form in spin net models and the block labels can be interpreted as intertwiner channels.This enforces the hope that spin nets and spin foams do indeed share statistical properties [13,14,17].

Improving the truncation
In order to improve the truncation one would take into account more than one singular value per block in the process of splitting of cubes into prisms 6 .This requires the introduction of additional indices on the faces, here denoted by I f . .., that describe the correlations between the two split parts of the cube.
In the following iterations, faces are split into triangles and hence face indices have also have to be 'divided': We do so by putting the face index I f onto the newly introduced diagonal edge e and denote it by I e (= I f ), the diagonal e now carries an index pair (k e , I e ).Thus if we glue the entire complex, the index on the diagonal edge will ensure that the four triangular faces that are glued in pairs onto each other, carry the same index I e , as is the case before cutting.Eventually, all the faces will carry variables I f , and all the edges will have variable pairs (k e , I e ).We define the initial amplitude as where the introduced indices can be understood as (initially empty) higher-order corrections decorating the cube amplitude, cf.figure 6.
Let us return to the process of splitting the cube into prisms, which is implemented by a SVD.There are now three independent indices k e and four indices I e that are shared by both prisms.These indices could therefore be used as parameters for the SVD, as described before in the case of just having the three k e -indices.One would now have K 3 × χ 4 e blocks where χ e is the lengths of the I e indices.Assuming we choose also χ e singular values per block we would take into account K 3 × χ 5  e singular values in this splitting and thus have an 'effective' bond dimension of χ eff = K 3 × χ 5 e .However we have now a situation where one expects that the singular values from a block labelled by I e = 2 for all four edges are much smaller than from the block labelled by I e = 1.This can be accommodated by comparing the singular values for all block labels among each other and selecting the χ largest ones, where χ is now a choice for the effective bond dimension.The χ largest singular values have to be distributed among the various blocks, resulting into blocks of varying size and also into vanishing blocks if these do not get any singular values assigned. 7The size of an index I e can thus get reduced a posteriori.This makes the algorithm much more complicated to implement.A possibility, to be explored in future work, is to fill entries up with zeros to reach a constant block size and then dealing with sparse matrices or tensors.
more complicated version [13,14] is to compare singular values across blocks and introduce a global cut-off for the number of singular values.Thus the size of the blocks might actually vary.Experience with the models in [13,14] has shown that reducing the size of a block might lead to instabilities, which can be addressed by introducing an algorithm in which the block size can increase, but not decrease. 7See [13,14,38] for an implementation of an algorithm with varying block size in 2D An alternative procedure is to keep just the block form with the k e indices as parameters and to accommodate the I e indices by first multiplying the amplitude with Kronecker delta's δ IeI e and then to split by an SVD with the I e indices assigned to one prism and the I e indices assigned to another prism.This automatically compares the singular values across the I e block structure, with a constant block size χ e (i.e.K 3 × χ e singular values are taken into account with each splitting) , however with the disadvantage that one has to split factors of Kronecker deltas with an SVD.
After the splitting of the cubes into prisms one proceeds as before, gluing four prisms to a cube in two steps.We sum over only the bulk variables, i.e.I e and k e labels on the half diagonals, as well as face indices I f for the faces without half diagonals.The most expensive step is the gluing of the two pairs of prisms into the cube, the cost of which scales with K 9 × χ 27 e .A reduction in costs can be obtained by summarising the variables on the top and bottom faces (with total bond dimension χ 4 e each), T and B respectively, into single face indices, T and B (treating the one independent k e index of one half diagonal as parameter).Such a truncation map, introducing a cut-off χ 4 e → χ e , can be implemented using the SVD: where R denotes all the remaining indices.After applying the truncation map on both sides (here (V † ) T T from the right and V B B from the left) the new truncated amplitude becomes In fact we could have also chosen U as a truncation map.We can apply the technique from [5] to determine which of the maps, U or V , gives a better truncation by comparing the root mean square values of the rejected singular values in both cases.
The simplification reduces the memory usage from O(K 9 × χ 24 e ) to O(K 9 × χ 18 e ), which allows for numerical testing of this higher-order algorithm.For the 3D Ising gauge model we obtain an improved critical temperature, β c = 0.710 for χ e = 2.
The memory and computational costs of this improvement scheme are quite high, also because additional indices are associated to edges and faces.An alternative scheme only introduces additional face indices.Imagine we put a vertex in the middle of the cube and attach a tensor with six edges carrying the I f -type indices.Then this splitting would correspond to splitting of a tensor network edge.To do so we proceed as discussed before: We formally replace the face index I f with a pair (I f , I f ) (for the two triangular faces) and multiply the amplitude of the cube with two Kronecker deltas δ(I, I ) for the two pairs of triangular faces.We then proceed with an SVD, treating the k e indices as parameters but not the two I f indices that are being split.Thus, the number of singular values taken into account at each splitting is K 3 × χ f , that is χ f singular values per block labelled by the three independent indices k e of the newly introduced face.This defines the effective bond dimension, i.e. the number of indices per tensor leg, as K 3 × χ f .
Proceeding as before with gluing the prisms to cubes, the cost of the algorithm goes as K 9 ×χ 14 f .We can reduce it further by applying the embedding maps on the top and bottom faces in two steps: χ 2 f → χ f after gluing two prisms, and once again after two pairs of them are put together.The overall complexity of the algorithm scales with K 9 × χ 10 f .By comparison the (higher order singular value) 3D algorithm of [5] has scaling of χ 11 for the computational costs (for models that can be brought into tensor form with bond dimension χ).For the Ising model we found a steady improvement of the critical temperature with increasing χ f , reaching β c = 0.729 for χ f = 8.

Further simplifications
The 3D algorithm still requires large computational resources, especially related to memory consumption of the decorated tensor, which largely pertains to dealing with the dimension higher than two.However, another reason is that the algorithm is based on cubes, which are not the most elementary building blocks.A 3D cube has 12 edges and 6 faces, whereas a tetrahedron has only 6 edges and 4 faces.The difficulty is to come up with regular coarse graining schemes involving tetrahedra.One straightforward scheme of coarse graining tetrahedra into tetrahedra right away is based on cubes subdivided into 6 tetrahedra each.In this case 8 smaller tetrahedra are coarse grained into a larger tetrahedron with all edges doubled in length.
Here we propose an algorithm based on prisms, which can be obtained from the cube algorithm described above by cutting the cubes.Instead of starting with cubes and subdividing these into prisms, we start with prisms right away, see figure 8.We divide two of these into pairs of tetrahedra and pyramids, before re-glueing into larger prisms back again.
The advantage is, that in the leading order approximation, the basic amplitude carries fewer indices, decreasing the computation cost to O(K 8 ), as well as the memory consumption from O(K 9 ) to O(K 7 ).While we find no qualitative differences to the cubic algorithm, the obtained critical temperature is near β c = 0.630 (within 18% of the Monte Carlo result [35][36][37]).Although less accurate, the technique might find applications in the models with large structure groups, where the cubic al-FIG.8. Simplified algorithm based on prisms, tetrahedra and pyramids.Basic building blocks are now prisms, which after splitting and re-glueing form the coarse-grained polyhedra.
gorithm would be too expensive.

OUTLOOK: NON-ABELIAN MODELS AND FOUR DIMENSIONAL MODELS
Here the initial amplitude of the cube depends not only on the representation labels ρ e attached to the edges but also on intertwiner labels ι f for the faces.Often, these face intertwiner labels can be chosen to be equal to representation labels ρ f , describing the coupling of a pair of edges via ρ f to the other pair of edges.If this pair of edges is chosen to be adjacent, then ρ f can be attached to a corresponding diagonal as in the Abelian case.Details for this algorithm for Non-Abelian groups will appear elsewhere, together with an application to spin foams with quantum group symmetry.
Let us also point out that the algorithms generalize to four dimensions.We have now hypercubes, the (in 3D) edge indices ρ e are replaced by plaquette indices (i.e.attached to the two-dimensional objects) ρ p and the (in 3D) face indices ρ f are now attached to the 3D cubes which make up the boundary of the hypercube.Thus the splitting and coarse graining can as before be considered to take place in a plane.

DECORATED TENSOR NETWORKS FOR OTHER MODELS
In this section we lay out that the main idea of the decorated tensor network algorithm can also be applied to other models, besides (generalized) lattice gauge theories.
Let us describe the general structure of decorated tensor models for lattice gauge theories: A tensor carrying the indices I, J, M, . . .sits in the centre of each cube with its edges piercing the faces of the cube, equipping them with indices.The tensors in neighbouring cubes are glued together by contracting the index of the edge connecting them, which is equivalent to summing over the index on the common face.In addition, this tensor network comes with supplemental decorations, i.e. the labels attached to the edges of the cube.In fact these are the variables which keep most of the physical information.Indeed one can imagine to contract the 'internal' tensor network by summing over all indices I, J, M, . . .and be left with a model with the original type of variables.
This model would in general be non-local, if the I, J, M, . . .indices have non-trivial range.Thus we can directly understand how non-local couplings, which appear in other real space renormalization schemes, are turned into additional indices and how the cut-off on these indices can also be understood as a truncation of non-local interactions that can occur, see also the discussion in [32,39].
Let us proceed with the application of this idea to Ising-like models in 2D 8 .The algorithm will amount to a sophisticated decimation procedure, i.e. keep a subset of the original Ising spins, separated at larger and larger distances, while the decimated spins are absorbed into a tensor network.These Ising spins can serve as signature observables, for instance one can straightforwardly compute a correlation function from the final tensor (this will be however at e.g. a distance half the system size in the contraction of six tensors to a torus topology).
Again instead of starting with a tensor or a vertex model, we rather think of squares as the basic building blocks.To these squares we assign an initial amplitude, which is a function of the four Ising spins sitting on the four corners of the square.Additionally, we will use a dual lattice for the tensor network, i.e. each square is decorated with a vertex in its centre with four emanating edges that cut the edges of the squares at the midpoints.The indices on these tensor edges will allow to go to a higher order approximation.Thus we again obtain the structure of a decorated tensor network.
The coarse graining algorithm is a straightforward adaptation of the Gu-Wen algorithm [4] taking the dec-oration by the Ising spins into account, see figure 9.The idea is to cut the squares along their diagonals into triangles by a SVD.As in the Gu-Wen algorithm, this requires to split the 4-valent tensor in the center of the square into two 3-valent ones, but additionally one has to 'split' the Ising spins connected by the diagonal.Naively, one can double these spins and identify them by introducing Kronecker Deltas, which allows us to straightforwardly perform the SVD, yet with a significant amount of irrelevant data in case the doubled spins are not equal.Instead it is more efficient to use the pair of split spins as block labels: the matrix that we intend to split decomposes into smaller block matrices, labelled by the pair of spins, on which we apply a SVD separately.Consequently the new index of the 3-valent tensor piercing the edge connecting the pair of split spins is precisely labelled by these spins.
After splitting the square along both diagonals, we obtain four triangular amplitudes which are connected to form the new (rotated) square as in figure 9.This incorporates the summation over four tensor indices, as in the Gu-Wen algorithm, and one Ising spin in the center of the new square.In the lowest possible approximation, i.e. tensor index range one, one only sums over the Ising spin.

Simplification -Triangular algorithm
In fact, the code depicted in figure 9, and also the related Gu-Wen algorithm, can be simplified by shifting the perspective from square building blocks to triangular ones.
If we consider the second part of the iteration, i.e. the glueing of the triangles and contraction of the respective indices, we realize that this can be performed in two steps: First we glue two small triangles into a larger triangle, then we glue the two larger triangles together to form the new square.However, in the next iteration of the algorithm we split this square again into the larger triangles we previously glued.
Hence, instead of performing the summations and subsequent splitting of the square, we propose an algorithm purely based on triangular amplitudes, see figure 11.Therefore we keep track of four different triangles, denoted as S 1 , . . ., S 4 : these amplitudes are pairwise glued together, once vertically, once horizontally, forming new larger triangles, where the new long edge carries two tensor indices and one Ising spin.To arrive back at an amplitude comparable to the initial one, we have to define a truncation map (three-valent tensors) blocking these indices into a new common tensor index, whose index range determines the quality of the approximation.Again these maps are computed via a (Higher Order) SVD as in [5]; note that the truncation maps applied to opposing triangles should be related by hermitian conjugation, such that they form a partition of unity (before implementing FIG. 9.The initial amplitude T associated to the square is split along its two diagonals.The two spins (gi, gj) connected by the diagonal get duplicated and label the blocks to which SVD is applied.This gives four new triangular amplitudes, S1, . . ., S4, which are glued together to a new square amplitude T .The last involves summation over the bulk tensor indices and one spin g.

the truncation on the number of singular values).
This modification of the algorithm is not restricted to decorated tensor networks and can be readily applied also to the original Gu-Wen algorithm, for which it leads to a reduction of computational costs from χ 6 to χ 5 .
The interesting feature of this decorated tensor network algorithm is the mixture of keeping part of the original variables and allowing for the additional tensor indices field redefinitions.An essential ingredient is a SVD per block as extensively used in the algorithms developed in [13,14,17,18].Especially for the non-Abelian models the knowledge of the singular values per block was key to be able to interpret the fixed point tensors.This might help to address some criticism to tensor network methods in [16], which points out the difficulty in extracting physical information from the fixed point tensors and separating physical from unphysical (due to the field redefinitions there is a kind of weak gauge symmetry for the tensor models) information.In terms of computation efficiency the algorithm is rather the same as the original Gu-Wen algorithm (apart from the reduction due to basing the algorithm on triangles instead of squares).
Results: We have applied tensor network methods to the 2D Ising model defined on a square lattice.To lowest approximation, i.e. χ = 1 per configuration of spins on the new edge corresponding to an effective bond dimension χ eff = 4, we observe a phase transition at the inverse temperature β ≈ 0.49664.In comparison to the exact solution, 2) 2 ≈ 0.44068, this is an error of roughly 13 %, yet still an improvement over the Migdal-Kadanoff methods, where one finds β (MK) ≈ 0.3047, an error of roughly 31 %.Additionally, we can compute critical exponents, which essentially confirm the results presented in [16] on the Gu-Wen algorithm on the square lattice for the lowest cut-off.Exploiting the fact that decorated tensor networks preserve (some) of the original Ising spins, we can compute correlation functions.In figure 12 we show the correlation functions (and its first derivative) of two 'neighbouring' FIG.11.
Triangular version: Starting from four triangular amplitudes, S1, . . ., S4, we glue them once vertically and once horizontally into larger triangles.The long edge carries two tensor indices and one spin, which we combine into a new tensor index by the embedding maps U , U † and V , V † respectively.These larger triangles, S 1 , . . ., S 4 , are the new triangular amplitudes for the next iteration.
spins on a 2D lattice with periodic boundary conditions after different numbers of iterations, i.e. different distances between the spins on the original lattice.This correlation function clearly indicates the phase transition between the ordered and disordered phase, where the correlation functions intersect right on the phase transition because of scale invariance of the fixed point tensor.
Increasing the cut-off the phase transition temperature approaches the exact result; the highest cut-off tested is χ eff = 60, for which we observe a phase transition around β ≈ 0.440706 (an error of roughly 0.06 %).
For the 3D Ising model our method gives in the lowest order approximation a phase transition between The correlation functions (top) and their first derivative (bottom) for two neighbouring Ising spins evaluated for amplitudes after different number of iterations, here n = 20, 21, 22, 23, in the lowest order approximation, χ eff = 4, with periodic boundary conditions.Due to the different number of iterations, the correlations functions correspond to different system sizes and different distances of the two spins on the original lattice.The more iterations have been performed the steeper the correlation function becomes on the phase transition and its first derivative diverges.The correlation functions for different number of iterations meet in one point, which marks the phase transition temperature.0.25 − 0.251, whereas Monte Carlo approaches [35][36][37] give around β c ≈ 0.2217.

DISCUSSION
Local gauge symmetry, and thus a redundancy of information, in lattice gauge theories poses a serious challenge to tensor network algorithms aiming to efficiently coarse grain these systems and extract their effective dynamics at larger and larger scales.Yet their capability to also deal with negative and complex amplitudes, as encountered in fermionic systems and spin foams in quantum gravity, adds to their desirability.
Therefore we have introduced and tested a novel tensor network renormalization algorithm, called decorated tensor networks: Instead of encoding all variables into the tensor network and thus potentially blurring the gauge symmetry (e.g.represented in the form of Gauss constraints for the representation labels), in particular after several coarse graining steps and variable redefinitions, we base our algorithms on building blocks of the discretisation, here cubes in the 3D cubical lattice, which are coloured by the original variables of the system and carry a tensor in their center; hence the tensor gets 'decorated' by the original variables of the system.
We tested the method with the Ising gauge model in 3D, which led to acceptable results.(We are not aware of other results with tensor network algorithms applied to gauge theories.)The method seems to be in particular effective for large structure groups, which will appear in spin foams and lattice gauge theories (in fact with infinite initial bond dimension).For spin foams decorated tensor networks provide in fact the first systematic coarse graining scheme.
We presented several possibilities to improve the lowest order approximations.We observed a systematic improvement of the phase transition temperature with increasing cut-off.Again it is difficult to compare with other tensor network methods -general discussions suggest that implementation of entanglement filtering [22,40] might improve very much the effectiveness of 3D algorithms.
The decorated tensor network scheme allows for a clear identification and preservation of (part of) the original variables of the system.As an additional advantage the fixed points of this renormalization procedure are more straightforward to interpret.Note that the tensor network itself plays a distinctly different role in this algorithm: Instead of wholisticly encoding the dynamics of the system, it rather appears as way of locally encoding higher order corrections.To lowest cut-off, i.e. if all tensor indices are trivial, the systems keeps its original form.Higher order corrections can be understood as a higher multiplicity of particular colouring / configurations of the building blocks.
In fact, a slightly different perspective can be taken as well: the network of tensors, located in the center of the building blocks, can be contracted at any point.After this contraction one arrives back at a system depending on the same type of variables as the original ones, yet with generically non-locally interacting building blocks.Hence one can understand the decorated tensor network as encoding non-local interactions between the building blocks, where the cut-off on the tensor indices can be directly translated into a truncation of the non-local interactions, see the discussions in [32,41].This perspective will help to compare tensor network results with other, for instance continuum methods, based on truncations of (non-local) coupling terms in the effective action.
In this vain we hope that decorated tensor networks might eventually facilitate new methods mixing both analytical (for the decoration indices) and numerical tech-niques (for the tensor indices).This can possibly provide new approaches to go beyond finite groups towards Lie groups with infinite initial bond dimension and to allow to deal with (gauge) divergences occurring in spin foam models [23][24][25].
Regarding the application to spin foams, the decorated tensor networks allow to keep the geometric interpretation of the spin foam variables.The coarse graining process might thus allow to conclude whether the blocking of variables (as encoded in so-called embedding or truncation maps [32,41]), determined by the SVD truncation, also proceeds in a geometric manner.As models of quantum gravity spin foams pose additional challenges.One is that the notion of scale, in terms of complexity of boundary data, as well as diffeomorphism symmetry have to be emergent [41][42][43].To this end non-local amplitudes are unavoidable [44,45], which indeed can be taken care off with the decorated tensor networks.As described in [41] tensor networks can be understood as a truncation scheme to compute physical states and a physical vacuum.We hope that the decorated tensor networks allow in particular to retain the geometric interpretation of the initial amplitudes, with the actual tensor network providing an improvement [46,47] over the bare (initial) amplitudes.

FIG. 1 .
FIG. 1. Projector Pe is associated to each edge of the lattice.The magnetic indices are contracted with neighbouring projectors sharing the same face f .In case the orientations of the faces f ⊃ e and the edge e match, i.e. o(e, f ) = 1, ∀f ⊃ e, then Pe : Inv f ⊃e Vρ f → Inv f ⊃e Vρ f , the contragredient representation has to be used for opposing orientations.

FIG. 2 .
FIG.2.Tensor network representation of a spin foam defined on a cubical lattice.The tensors Te on the edges capture the intertwiner degrees of freedom, while the tensors T f on the faces are an auxiliary structure.For a given face f they ensure that all edge tensors Te with e ⊂ f carry the same representation label ρ f assigned to that face.

FIG. 10 .
FIG.10.Schematic illustration of the 3D algorithm: Similar to the gauge theory version, cubes are split along their diagonal faces, here labelled by four spins that are split into two.As a result, one considers generalized cubic amplitudes, where the 'crossed' carry an additional Ising spin.