Self-assembling tensor networks and holography in disordered spin chains

We show that the numerical strong disorder renormalization group algorithm (SDRG) of Hikihara et. al. [Phys. Rev. B 60, 12116 (1999)] for the one-dimensional disordered Heisenberg model naturally describes a tree tensor network (TTN) with an irregular structure defined by the strength of the couplings. Employing the holographic interpretation of the TTN in Hilbert space, we compute expectation values, correlation functions and the entanglement entropy using the geometrical properties of the TTN. We find that the disorder averaged spin-spin correlation scales with the average path length through the tensor network while the entanglement entropy scales with the minimal surface connecting two regions. Furthermore, the entanglement entropy increases with both disorder and system size, resulting in an area-law violation. Our results demonstrate the usefulness of a self-assembling TTN approach to disordered systems and quantitatively validate the connection between holography and quantum many-body systems.


I. INTRODUCTION
There is currently a lot of excitement around the socalled AdS/CFT correspondence and possible applications in condensed matter physics. 1 The AdS/CFT correspondence is most well known in high energy physics where it was noted 2 that there exists a duality between certain theories of gravity on D + 1 dimensional Anti de Sitter (AdS) spacetime and conformal quantum field theories (CFT) living on its D dimensional boundary. In condensed matter systems, the AdS/CFT correspondence can provide a geometric interpretation of renormalization group (RG) techniques since the additional holographic dimension can be interpreted as a scale factor in the RG coarse graining. 1 It has been argued recently 3,4 that certain RG approaches to the Hilbert space of critical many-body interacting system in D dimensions, such as the multi-scale entanglement renormalisation ansatz (MERA) for tensor networks, share many of their geometric properties with D + 1 dimensional AdS. This connection is based on ideas 5 that suggest that the entanglement entropy of a region on the boundary is related to the minimal surface in the holographic bulk that separates the region from the rest of the surface. These ideas were further developed by Evenbly and Vidal 6 to discuss the underlying geometric structure of entanglement and correlation functions in such tensor networks in general.
Tensor network methods provide elegant and powerful tools for the simulation of quantum many-body systems. Their original manifestation, the density-matrix renormalization group (DMRG) 7 is now understood to be based on a variational update of a matrix product state (vMPS), 8,9 and has found applications in a wide range of fields such as quantum chemistry 10 and quantum information 11 as well as condensed matter physics. 12 More recent developments have extended the methods to, e.g., critical systems, 13 two-dimensional lattices [14][15][16] and topologically ordered states. 17 For disordered quantum many-body systems, the strong-disorder renormalization group (SDRG) provides a similarly unifying approach. 18 It was originally devised by Ma,Dasgupta and Hu 19,20 for the random antiferromagnetic (AFM) Heisenberg chain where s i is the spin-1/2 operator and J i is the coupling constant, which takes a random value between 0 < J i < J max according to some probability distribution P (J). The principle behind the SDRG is to eliminate the most strongly coupled pairs of spins and replace them with an effective interaction that couples the spins at either side of the pair, as shown in Fig. 1a. The pair of spins coupled by J max are thought as being frozen into a singlet ground state as the neighbouring interactions are significantly weaker -ultimately leading to the random singlet phase, which is the ground state of the system. 21,22 This freezing of degrees of freedom is remarkably close to an update process in entanglement RG for tensor networks 13 and suggests the possible usefulness of the AdS/CFT correspondence also for disordered spin chains. By analysing the probability of survival through the SDRG algorithm it is possible to predict that mean correlations will have a power-law decay 21 with negative power 2. Similarly, the entanglement entropy can be shown to scale logarithmically with block size, 23 where the amount of entanglement between blocks A and B is quantified by the Von Neumann entropy with ρ A the reduced density matrix obtained by tracing over the B components of the density matrix.
In this work, we have developed a self-assembling tree tensor network (TTN) algorithm based on the previous ideas of SDRG. 24,25 This allows us to calculate properties such as expectation values, correlation functions and entanglement entropy directly and efficiently from the geometry of the TTN. In particular, we find that the distance dependence of the spin-spin correlation function can be studied not only via direct calculation of the correlation functions, but also via the holographic distance dependence along the tree network connecting two sites. In section II we will briefly review the numerical strong disorder renormalization group of Hikihara et. al. 25 and define the states and operators that form the basis for our work. Section III shows how the numerical SDRG on a matrix product operator (MPO) self-assembles the TTN. Finally, in section IV we compute correlation functions and entanglement entropy (i) directly using the TTN as well as (ii) via simply counting the path lengths and connectivities in the holography. We find that both approaches give consistent results.

II. THE MPO IMPLEMENTATION OF THE SDRG
A. The numerical SDRG The SDRG method was extended to both ferromagnetic (FM) and anti-ferromagnetic (AFM) couplings by Westerberg et. al. 24,26 The approach finds the neighboring pair of spins s i , s i+1 with the greatest energy gap ∆ i between the ground state and excited state and combines them into a single effective spinS (Fig. 1b). The effective couplings between the new spin and its neighbours are then recalculated using Clebsch-Gordan coefficients and the new gaps∆ i−1 and∆ i updated. SDRG was once more extended by Hikihara et. al. 25 to include higher states at each decimation, in the spirit of the numerical renormalization group (NRG) 27 and the DMRG. 7 This method therefore decomposes the system into blocks rather than larger spins allowing for more accurate computation of, e.g., the spin-spin correlation functions. The more states that are kept at each decimation the more accurate the description and is exact in the limit of all states kept.
Consider a point in the algorithm where the Hamiltonian is made up of blocks H B i at each site and couplings H C i,i+1 between them as in Fig. 1. The couplings take the form of a two spin Hamiltonian where s R i is the spin operator of the right hand spin of block i and s L i+1 is the left hand spin of block i + 1. In full the Hamiltonian is where N B is the number of blocks.
Let us now define the gap ∆ i as the energy difference between the highest energy SU(2) multiplet that would be kept and the smallest multiplet that would be discarded in a renormalization of block H B i,i+1 . The scheme works by searching for the pair of blocks with the largest gap ∆ im and then combines the coupling and the blocks that it connects into a single block This block and the couplings either side are then renormalized by a matrix (V χ ) of the eigenvectors corresponding to the lowest χ eigenvalues of the block, such that only full SU(2) blocks are kept. The process is repeated until the system is represented by one block. The details of the algorithm are described in appendix A.

B. Numerical SDRG as an MPO process
Hikihara's numerical SDRG can be naturally described as a set of operations on an MPO (see Appendix B for more details). First, we contract the MPO tensors for the pair of sites with the largest gap, sites i m and i m + 1 (Fig. 2a) Here we have σ im = 1, . . . , χ for the physical indices and, for the Heisenberg model (1), the virtual indices are b im = 1, . . . , 5. Next, we perform an eigenvalue decomposition on the on-site components of the new MPO tensor keeping the eigenvectors of the lowest χ eigenvalues (V χ ) As with the Hikihara's algorithm, only the χ eigenvalues that make up full SU(2) multiplets are used. Then we contract V χ and V † χ with the new MPO tensor to perform the renormalization (Fig. 2b). For the moment write the two-site combined MPO W [im,im+1] in terms of an effective site with index τ = 1, . . . , χ, χ+1, . . . , χ 2 , i.e. W τ,τ bi m−1 ,bi m +1 . Similarly, we can write the set of eigenvectors as [V χ ]σ im τ . Then the contraction is explicitly given as whereσ im = 1, . . . , χ is the spin index of the renormalised site i m . Hence we replace sites i m and i m +1 with a single renormalized site and relabel the remaining indices. The contraction makes the on-site component of the new MPO simply a diagonal matrix of the lowest χ eigenvalues Λ χ (Fig. 3a). It also has the effect of renormalizing the coupling spins just as in the Hikihara approach (Fig.  3b) The contraction therefore maps two MPO tensors onto one while preserving the indexing structure of the MPO. As the final step, we diagonalize the neighbouring blocks to update the distribution of gaps. The procedure is then repeated until the system is just one site, and we diagonalise to obtain the ground state energy E g of the system.

III. TREE TENSOR NETWORKS AND SDRG
The MPO description of SDRG given above amounts to a coarse-graining mechanism that acts on the operator. Alternatively, we can view it as a multi-level tensor network wavefunction acting on the original operator. To illustrate this, we can split the τ index of V † χ as in Eq. (8) back to the original spin indices σ im , σ im+1 to create an or ww † = 1 1 = w † w (Fig. 4a). A renormalization in the SDRG algorithm as in Fig. 3b+c can then be rephrased graphically as in Fig. 4b. This makes the notion of mapping two MPO tensors to one immediately explicit. When viewed in terms of isometries, the algorithm can be seen to self-assemble a tensor network based on the positions of largest gaps before each renormalisation. When written in full, it builds an inhomogeneous binary tree tensor network (TTN) as shown in Fig. 5. Tree tensor networks are one of the major areas of tensor network research and TTNs with regular structures have been extensively studied. [29][30][31][32] The isometric nature of the isometries allows for calculations to be performed in a highly efficient manner. 28,30,33 When calculating expectation values, such as the two-point correlation function (Fig. 6a), only those tensors that effect the sites that the operators act on need to be included, this is known as the past causal cone 28 and is drawn as a blue shadow in the holographic bulk. This allows for a reduction in the number of contractions that need to be performed to obtain a result. Calculation of the entanglement (Von Neumann) entropy (2) can also be made more efficient as shown in Fig. 6b. In addition to the reduction due to the isometries, we note that the entanglement entropy is not affected by the isometries acting just on A. 30 However, the entries in the density matrix will change so we label it A .

IV. RESULTS
In the following, we shall compare results for the disordered anti-ferromagnetic Heisenberg model (1) when using a modern DMRG implementation, e.g. variational MPS (vMPS), with those obtained from our TTN SDRG strategy (tSDRG). The set of couplings (J i ) shall always be taken from a box-type distribution, 25 i.e. constant in the range 0 < 1 − ∆J/2 < J i < 1 + ∆J/2 < 2 and zero outside. Unless stated otherwise, we use strong disorder ∆J = 2 − in the following. We assume open (hard wall) boundary conditions throughout.

A. Convergence and ground-state energies
In Fig. 7 (main), we show show the dependence of the disorder averaged ground state energy per site, E g /L, on ∆J for constant L. We find that for both vMPS and tSDRG, the E g /L values decrease for increasing ∆J, i.e. the ground state energy lowers as disorder in the J i couplings allows the system to form particularly energetically favourable spin configurations. We also see that the vMPS for the chosen values of χ and L reaches lower energies. This suggests that it is yet more efficient in finding an approximation to the true ground state energy. However, upon increasing ∆J, the difference between vMPS and tSDRG is getting smaller. This is expected since SDRG is based on the idea that the contribution from the non singlet interactions is small, which is more accurate an assumption the greater the disorder. The figure also shows that increasing χ can considerably improve the results of the tSDRG. 34 In Fig. 7 (inset) we show E g /L as a function of L for various values of χ at the strongest permissible disorder ∆J = 2 − . We find that the values of E g /L have do not vary much anymore for system sizes L ≥ 100. Conversely, E g /L values for L < 100 are clearly dominated by the presence of the open boundary conditions.

B. Correlation functions
The correlation functions for a strongly disordered Heisenberg chain are expected to average out to be a power-law decay 21 where s x1 · s x2 is understood to be the disorder averaged expectation value of the two point spin-spin correlation function. This r −2 scaling of the correlation is The diagram in the right-hand side of (b) has been reduced in the horizontal direction to highlight the reduction in complexity due to the isometries. a feature of the disorder in the system 21 and should be contrasted with the well-known power-law dependence of correlation functions 35 in interaction-driven Luttinger liquids. 36 In loop-free tensor networks, correlations scale as e −αD(x1,x2) , where D(x 1 , x 2 ) is the number of tensors that connect site x 1 to x 2 . 6 DMRG is based on the MPS and as such it has one tensor per site, i.e. D MPS ≈ |x 2 − x 1 |. Therefore correlations in DMRG scale exponentially. This suggests that for long chains it will be necessary to keep large numbers of states to be able to model a power-law correlation of the system. 6 tSDRG on the other hand has a holographic geometry based on a random TTN with path length D TTN ≈ log |x 2 − x 1 |, i.e. scaling logarithmically with distance when averaged. This makes it much more suited to capture the desired power law decay In Fig. 8, we show the behaviour of s x1 · s x2 com- puted directly as well as its holographic estimate based on (13). We find that the behaviour for |x 2 −x 1 | 1 and |x 2 −x 1 | < L/2 is indeed very similar for both approaches. The best fit value for α is 0.62 ± 0.02 where the error is the standard error. 37 We find that in the indicated distance regime, both measures of s x1 · s x2 are consistent with the expected r −2 behaviour. For |x 2 − x 1 | > ∼ L/2 we see that the boundaries lead to an upturn on the behaviour of s x1 · s x2 for both direct and holographic estimates. This upturn is a result of boundary effects and can easily be understood in terms of the holographic TNN: for |x 2 − x 1 | ≥ L/2, the average path length in the tree decreases (cp. Fig. 6). This is also consistent with periodic systems where we expect correlation functions to be equal for |x 2 − x 1 | = r and L − r. In the inset of Fig. 8 we show the distance dependence of D TTN with χ = 10. For |x 2 − x 1 | < L/2, the data can be described by as linear behaviour in log |x 2 − x 1 | with slope 2.94 ± 0.02. Note that this slope along with the value of α = 0.62 ± 0.01 gives an estimate of power-law exponent a = 0.62 × 2.94 = 1.84 ± 0.04 for fixed L = 150. Figure 9 shows that as L increases, the resulting value of the scaling power a also increases towards the expected value of 2 for larger systems upon increasing L. We have also checked that the differences between χ = 10 and 20 remain within the error bars and hence we use χ = 10 for calculations of s x1 · s x2 in Fig. 8.

C. Entanglement entropy
In general, the entanglement entropy S A|B is difficult to compute as the size of the reduced density matrix ρ A scales exponentially with the size of block A. While for special cases, such as the XX model, 38 S A|B can be computed more easily, the general strategy involves finding the eigen-or singular values of ρ A . 9 The TTN representation of tSDRG gives an alternative means of calculating S A|B for any bipartitions A and B of the system. In a similar manner to the correlation functions, the geometry of the tensor network is related to its ability to capture S A|B . Briefly, S A|B is proportional to the minimum number of indices, n A , that one would have to cut to separate a block A of spins from the rest B of the chain (cp. Fig. 6). 6 This dependence is related to the famous area law, which states that for the ground state of a gapped system, the entanglement entropy of a region is proportional to the size of the boundary that separates the two regions. 39,40 The MPS is a simple line of tensors (cp. appendix B) and thus the number of indices that separate one region from another is a constant and independent of the size of the block and its position in the chain. Unlike the MPS, for the TTN the position of the block in the chain alters the number of indices that have to be cut to separate it from the rest of the system. This suggests that there are spatial regions in the chain that are more and less entangled, which is likely to be true for a strongly disordered spin chain. The concept is hence similar to discussing the entanglement in the Ma, Dasgupta, Hu implementation 19 of SDRG, where the entanglement entropy is related to the number of singlets that have to be broken to separate a region from the rest. 23 In Fig. 10 we show that the average value of S A|B remains approximately constant upon increasing the disorder, while the average of the maximal S A|B shows a pronounced increase. This indicates that the full distri- bution of S A|B develops long tails with large S A|B values when increasing ∆J. For strong disorders ∆J > ∼ 1.5 we find that tSDRG is more accurate than vMPS. The vMPS estimates of S A|B are consistently below the values obtained by the tSDRG. Only when increasing χ do we reduce the deviation. This behaviour is most pronounced for the average of the maximal S A|B values. For example, with χ = 20, the S A|B values obtained for vMPS deviate from the tSDRG results around ∆J ≈ 1.2. Hence we see that an increase in S A|B requires a considerable increase in χ for vMPS to accurately capture entanglement. On the other hand, for weak disorders ∆J < ∼ 0.5, vMPS gives consistent results already for small χ = 10. The values obtained for S A|B from tSDRG are much higher in this regime. We believe this to be an overestimation of S A|B by the tSDRG because, as discussed before, tSDRG selects most strongly the singlet pairs in the disordered system, which of course become less prevalent for low disorder. Figure 11 shows that when L is increased for ∆J = 2, both the average and average peak values of S A|B increase logarithmically in L. This again implies that as L is increased, the χ value for vMPS needs to be increased also to be able to capture the entanglement. On the other hand, the holographic nature of the TTN means that the minimal surface in the network increases with system size and thus describes this entanglement without the need to increase χ. Although S A|B is therefore captured well by the network, contracting ρ A for larger L becomes more increasingly difficult, even with the simplifications suggested in section B, since the size of the matrices scales as O(χ nA ). We therefore have to restrict ourselves to smaller χ and L values than in sections IV A and IV B.
where region B is a block of extent L B in the centre of the spin chain. Note that this implies an effective central charge 23 ofc = 1 · log 2. This is different from the bipartition entanglement S A|B that we considered before. We show the resulting S A,B in Figure 12. The figure clearly indicates that finite size effects become prevalent for large L B , so we fit for L B ≤ L/2 only. The resulting scaling behaviour S A,B ≈ (0.22 ± 0.02)log 2 L B is consistent with Eq. (14). We finally also examine the entanglement entropy per bond, S/n A , of a TTN for both bipartitions A|B and blocks A, B with χ = 10 when averaging over 500 disorder configurations with L = 50. Figure 13 shows that away from the boundaries S/n A saturates to the same constant 0.47 ± 0.02 for bipartitions and blocks. 42 This is consistent with Ref. 6 and implies that the entanglement entropy is proportional to the length of the holographic minimal surface that connects the two blocks. Note that for L B ∼ L/2, we find that up to 20% of our samples for χ = 10 lead to calculations of S A,B consuming memory beyond 100GB. This is currently out of reach for us and we disregard the configurations. Nevertheless, we think that this is purely a numerical artefact and does not change the average values of S A,B /n A reported here. Calculations with smaller χ confirm this. 42

V. CONCLUSION
In this work, we demonstrate the validity and usefulness of a suitably adaptive tensor network approach to locally disordered one-dimensional quantum many-body systems. In contrast to traditional vMPS approaches to disordered systems, where the initial geometry of the MPS ignores the disorder and only takes it into account at the stage of variational sweeps, 43 our approach incorporates the disorder into the fabric of its tensor network. We believe this strategy to be inherently more suited to disordered systems -the results presented here show that the accuracy of tSDRG is already comparable to vMPS without including any additional variational updates. This advantage is particularly evident for longranged correlations and an entanglement entropy that violates the area law.
Our results furthermore show that, when disorder averaged, a random AFM spin 1/2 system is well characterised by an effective CFT on the boundary of a discretized holographic bulk. We believe, to the best of our knowledge, that we have thus shown the quantitative validity of holography for the first time here. In particular, our spin-spin correlation function, Fig. 8, as well as the block and bipartition entanglement entropies, Fig. 13 show excellent qualitative and numerical agreement with their holographic counterparts. Such an agreement also reconfirms that the self-assembly of the TTN produces the necessary tensor network geometry.
Our local RG procedure selects spin pairs based on en-ergy gaps. It is tempting to reformulate this based on the local entanglement content of such pairs. However, it is not straightforward to find such a local measure that captures energies and wave functions well simultaneously.
In particular, we do not find a convenient local entanglement measure that would have a simple relation to the local values of J i . More promising might be the implementation of a variational TTN. 30 Our initial results suggest that this does indeed improve the energy values, but at considerably increased efforts in implementation and computation -every disorder configuration of course necessitating its own variationally updated tree structure.

ACKNOWLEDGMENTS
We are grateful to Andrew Ferris and Nick d'Ambrumenil for valuable discussions. We would like to thank the EPSRC for financial support (EP/J003476/1) and provision of computing resources through the Mid-Plus Regional HPC Centre (EP/K000128/1). AMG would like to thank the organisers and participants of the Networking Tensor Networks 2012 workshop at the Centro de Ciencias de Benasque.
where σ i are the physical indices of the lattice and enumerate the states in the local Hilbert space. The tensor C σ1...σ L can be decomposed into a tensor network, the most common of which is the Matrix product state (MPS) (B2) Here, a i = 1, . . . , χ for site i in the bulk. It is convenient when studying tensor networks, such as the MPS, to give the equations a diagrammatic form (Fig. 14a). Each tensor is drawn as a shape where each line coming out represents an index and connected lines represent tensor contractions.

Matrix product operators
In a similar manner to the matrix product state, operators acting on lattice wavefunctions can be decomposed into a network of more simple tensors. A general operator on a lattice can be written as:  Simply multiplying W b1 W b1,b2 · · · W b L−2 ,b L−1 W b L−1 results in (B5). The top right element of Eq. (B7) and equivalent elements in Eqns. (B6) and (B8) are referred to as the on-site elements. This is where an external magnetic field of the form h i S z i would be introduced. Furthermore, it is possible to include longer range interactions in the elements away from the top and right row and column. 44,45 Another way of describing the contents of an MPO is a matrix product (MP) diagram. 44 This is a pictorial representation of the elements in the tensor, whereby the indices are numbered circles and the corresponding elements are paths that connect any two indices (Fig. 15). Matrix multiplication, or contraction, is then represented by the sum of the unique paths that connect the indices on the far left and right of the diagram when multiple matrices are placed end to end (Fig. 16). For MPOs it is understood that the binary operator between terms is a tensor product. Tracing out the different paths in Fig. 16 results is the standard form (B5) with L = 4. The MP diagrams give a convenient means of visualising the components of the MPO and are particularly useful  when creating operators with long range components or periodic boundary conditions.