Entanglement renormalization for disordered systems
Andrew M. Goldsborough Glen Evenbly






[]andrew.goldsborough@mpq.mpg.de
[]www.goldsborough.org.uk


Max-Planck-Institut für Quantenoptik, Hans-Kopfermann-Str. 1, D-85748 Garching, Germany


[]glen.evenbly@usherbrooke.ca
[]www.glenevenbly.com
Département de Physique and Institut Quantique, Université de Sherbrooke, Québec, Canada




We propose a tensor network method for investigating strongly disordered systems that is based on an adaptation of entanglement renormalization [G. Vidal, Phys. Rev. Lett. 99, 220405 (2007)].
This method makes use of the strong disorder renormalization group to determine the order in which lattice sites are coarse-grained, which sets the overall structure of the corresponding tensor network ansatz, before optimization using variational energy minimization. Benchmark results from the disordered XXZ model demonstrates that this approach accurately captures ground state entanglement in disordered systems, even at long distances. This approach leads to a new class of efficiently contractible tensor network ansatz for 1D systems, which may be understood as a generalization of the multi-scale entanglement renormalization ansatz (MERA) for disordered systems.


75.10.Jm, 05.30.-d, 02.70.-c

Entanglement renormalization for disordered systems
    Glen Evenbly
    December 30, 2023
===================================================






§ INTRODUCTION


Disordered systems are naturally of interest in quantum many-body physics as real world materials are rarely perfect. Furthermore it is known that the presence of disorder can dramatically alter the properties of a system, as evidenced in the ground-breaking work of Anderson <cit.> which that showed that the amount of disorder in a non-interacting three dimensional electron system dictates whether the states are localized or extended. More recently there has been great excitement around the idea of many-body localization, which extends Anderson localization to interacting many-body systems <cit.>.

Unfortunately, disordered systems can also be challenging to study using standard approaches. For instance the density matrix renormalization group (DMRG) <cit.>, which is otherwise regarded as an extremely reliable and robust method for 1D quantum systems, may have difficulty even for a relatively simple example like a 1D disordered Heisenberg model due to the presence of long range singlets in the ground state <cit.>. On the other hand, specialist methods developed for disordered systems such as the strong disorder renormalization group (SDRG) method of Ma, Dasgupta and Hu <cit.> may capture the qualitative features of the ground state entanglement, but generally only produce quantitatively precise results when the disorder is infinitely strong.


In this manuscript we propose and benchmark a new tensor network method designed for disordered systems. In most previous tensor network approaches the structure of the network, i.e. the pattern of how the tensors in the network are connected, is fixed independent of the problem under consideration whilst the parameters within the tensors are optimized variationally. Conversely, the method we introduce here is based on allowing both the structure of the tensor network as well as the parameters within the tensors to be adapted to the problem under consideration. This idea is similar to previous work <cit.>, which adapted a coarse-graining transformation, based on blocking together sites using isometries, to specific instances of a disordered system, and thus produced an approximation to the quantum ground state as a disordered tree tensor network (TTN). 
However, in the present work the coarse-graining transformation is based on entanglement renormalization <cit.>.
This differs by including unitary disentanglers, such that the ground state is approximated as a generalized form of the multi-scale entanglement renormalization ansatz (MERA) <cit.>.
Here we use SDRG to adapt the overall structure of the network to the problem under consideration, which helps ensure the qualitative features of the ground state entanglement are captured, before then applying variational sweeping to optimize the tensors and improve the quantitative accuracy significantly.

The manuscript is organized as follows: Section <ref> provides an overview of some of the previous approaches to numerical strong disorder renormalization, and Sect. <ref> discusses the problems of trying to combine SDRG with a variational tensor network algorithm. Section <ref> describes the newly proposed tensor network method for disordered systems, also outlining the variational optimization algorithm that is employed. Finally, Sect. <ref> presents benchmark results for the disordered XX and Heisenberg (XXX) models, which are compared to various other numerical methods.




§ STRONG DISORDER RENORMALIZATION






In this section we review some of previous approaches proposed for strongly disordered systems. Let us consider the disordered XXZ model on a 1D lattice of spin-12 sites, with periodic boundary conditions (PBCs). The Hamiltonian is defined as

    H = ∑_i J_i{1/2[ s^+_i s^-_i+1 + s^-_i s^+_i+1] + Δ_z s^z_i s^z_i+1}

where s^x, s^y and s^z are spin-12 matrices, with s^± = s^x ± i s^y. Notice that the Heisenberg model when is recovered when the anisotropy is fixed at Δ_z = 1, while the XX model is when Δ_z = 0. Disorder is introduced in the couplings J_i which are allowed to vary in strength with position; here we choose these couplings randomly between 0<J_i<J_max according to some probability distribution P(J). 

Many previous studies of the disordered Heisenberg model have used variants of the SDRG <cit.>, the archetypal method for the antiferromagnetic case was set out by Ma, Dasgupta and Hu <cit.>. In this procedure the pair of spins with the strongest coupling J_i are approximated as a singlet, and then second order perturbation theory is used to find an effective coupling J̃_i between the spins either side, as illustrated in Fig. <ref>(a). This step is repeated to produce a random singlet phase, an example is given in Fig. <ref>(a), where singlets created later on in the renormalization scheme can potentially span large distance scales. When averaged over many realizations of randomly chosen couplings J_i, the spin-spin correlation function is known to obey an inverse-square power law decay and the entanglement entropy scales logarithmically in a manner consistent with an effective conformal field theory (CFT) of central charge c=1 <cit.>.





More recent developments have included higher order components to increase the accuracy of observables such as the spin-spin correlation functions <cit.>, which is effectively an application of Wilson's numerical renormalization group (NRG) <cit.> to a system with many impurities. Later it was shown in Ref. GolR14 that this form of SDRG has an efficient tensor network representation in terms of an inhomogeneous TTN, where the structure is set by the renormalization order of the couplings, which we refer to as the tSDRG method. An example of a wavefunction produced from tSDRG is given in Fig. <ref>(b). The basic building block of these wavefunctions are isometries v, which annihilate to an identity with their conjugates v v^† = as also depicted in Fig. <ref>(a), and are responsible for coarse-graining two sites into a single effective site. In the tSDRG approach the isometries are found by diagonalizing the block Hamiltonian supported on the pair of sites with the strongest coupling. Whilst this method was shown to reproduce qualitative features of the correlation functions and entanglement entropy, the ground state energy was not given with high precision <cit.>. Although the accuracy of the energy could be improved through the subsequent use of a variational sweep to minimize the energy of the corresponding TTN, this was found to be at the expense of accuracy in the other features of the state. In particular, the energy minimization update favoured short-range correlations and ignored the long-range singlets, which are key to the proper characterization of the disordered ground states.








§ TENSOR NETWORKS AND DISORDER


The goal of this manuscript is to design a tensor network that can properly capture the structure of entanglement in a disordered ground state, in particular the pattern of long range singlets. 
Whereas the tSDRG approach constructs the TTN and sets the values in the tensors once, in this section we propose a more sophisticated network that allows a variational update whilst still capturing the long range entanglement.
It follows the ideas of entanglement renormalization, employing unitary disentanglers in such a way as to modify the random singlet phase of SDRG to be relevant for finite systems with finite disorder.


Each step of the SDRG can be interpreted as a coarse-graining transformation that maps a lattice Ł to a new lattice Ł' of two fewer sites by fixing two neighbouring spins, specifically those which possess the strongest coupling J_i, in a singlet state as depicted in Fig. <ref>(a). Here we propose two generalizations of this transformation, each a mapping from Ł to Ł'. The first, presented in Fig. <ref>(b), enacts a unitary disentangler u across the coupling J_i to be renormalized before employing isometries v_L and v_R that each coarse-grain two sites into one. The second, presented in Fig. <ref>(c), also enacts a unitary disentangler u across the coupling J_i before employing additional unitary gates w_L and w_R with the boundary spins. The sites at the location of coupling J_i are then fixed in an entangled singlet state, similar to SDRG. 
One should notice that if the spins associated to J_i in the second scheme, Fig. <ref>(c), were instead fixed as a product state then the first scheme of Fig. <ref>(b) is recovered.
This follows as unitary gates w_L and w_R with a fixed (product state) input index can be simplified to isometries v_L and v_R, a property that allows the MERA to be interpreted both as a tensor network and as a quantum circuit <cit.>.
As with other tensor network methods, we set an upper bound to the index dimensions of the tensors in order to truncate of the Hilbert space dimension for the coarse grained sites.
This is known as setting the bond dimension χ.

These coarse-graining schemes satisfy two important properties. Firstly, it is easy to see that there exist choices of tensors in both schemes for which the singlet state of Fig. <ref>(a) is reproduced, thus they should at least match the effectiveness of the standard SDRG. Secondly it should be noted that these schemes preserve locality: any operator supported on two contiguous sites of the initial lattice Ł is mapped to an operator supported on two contiguous sites of the coarser lattice Ł'. Hence it also follows that a nearest-neighbour Hamiltonian on Ł is mapped to a nearest-neighbour Hamiltonian on Ł'. Preservation of locality is important in order to produce a computationally viable numerical algorithm as, in general, the cost of representing a local operator grows exponentially with the size of its support.

The proposed lattice transformations of Fig. <ref>(b-c) can be seen to yield a class of tensor network ansatz, just as entanglement renormalization yields the class of tensor network known as the MERA. For a lattice Ł of N sites, we consider the class of tensor network that encompasses all possible ways of applying N/2 of the transformations from either Fig. <ref>(b) or Fig. <ref>(c) in order to reach a coarse-grained lattice of O(1) sites. An example of a tensor network from this class is depicted in Fig. <ref>(c). Importantly, one should notice that a standard MERA arises as a specific instance in this more general class. More precisely, the modified binary MERA, as introduced in <cit.>, can be understood as a particular (uniform) arrangement of the unit cell of Fig. <ref>(b). Accordingly, we refer to class of tensor network that generalizes the MERA for disordered systems as a disordered MERA or dMERA. Notice that, as with the standard MERA, all instances of dMERA possess a bounded causal structure, due to the preservation of locality in the underlying coarse-graining transformation, and are thus efficiently contractible for expectation values of local observables and for correlation functions. 




§  ALGORITHM


In this section we discuss an algorithm for how the proposed coarse-graining transformations can be implemented for the study of disordered systems. This algorithm consists of a procedure for first determining the overall geometry of the dMERA tensor network, before then applying variational energy minimization to optimize the tensors. 
Here we use the original SDRG algorithm <cit.> to determine the network geometry, then build the tensor network by patterning the unit cell of isometric and unitary gates, as depicted in Fig. <ref>(b-c), to match the order of singlets. 
An example of a random singlet state produced by SDRG and the corresponding dMERA is presented in Fig. <ref>(a,c). 
In general, we employ the transformation of Fig. <ref>(c) throughout, apart from when we wish to increase the bond dimension χ of the network (usually during the earlier RG steps).
The reason is that, despite the block in Fig. <ref>(b) reducing the geometric distance between any two sites connected by a singlet <cit.>, the update process still tends to omit long range correlation.
When the singlets are an explicit part of the coarse-graining block as in Fig. <ref>(c), the long range properties are preserved.

Once the structure of the dMERA has been set, we then adapt the standard MERA energy minimization algorithm <cit.> in order to optimize the tensors. Here we provide only an outline of the algorithm, the full details of which are provided in appendix <ref>. 
In a similar manner to standard MERA, we optimize the tensors by means of a variational sweep over the layers of the network.
Each layer is comprised of a coarse-graining block as depicted by Fig. <ref>(b-c), which transforms between an initial lattice Ł and a coarser lattice Ł'. Following the strategy proposed in Ref. EveV09 the tensors within the cell are updated by first computing their corresponding linearized environments from the coarse-grained Hamiltonian on lattice Ł and the two-site local reduced density matrices on lattice Ł', themselves obtained using the appropriate ascending and descending superoperators respectively. Each of the environments is then decomposed using a singular value decomposition (SVD) in order to find the new isometry or unitary that minimizes the (linearized) energy. The update is performed for all tensors within each layer, and sequentially over all layers, and this variational sweep is iterated until the tensors converge. All the manipulations required to optimize a dMERA scale as 𝒪(χ^7) in the bond dimension χ of the network.
 
Once the tensor network state is converged we then calculate expectation values and other properties of interest. Calculating expectation values of one- and two-site operators is straightforward; one can simply use the ascending superoperators to coarse-grain the operator through the network as discussed in appendix <ref>. As with optimization, this evaluation has the cost 𝒪(χ^7) in general but can be simplified to 𝒪(χ^6) for operators located on particular sites. Two-point correlation functions can also be efficiently computed by initially coarse-graining each of the operators separately and then fusing them into a single operator at the scale they meet. As with standard MERA, this may incur a higher computational cost as the fusion of the two operators may result in an operator with a larger support. As shown in appendix <ref>, the full cost of calculating the correlation function can vary between 𝒪(χ^6) and 𝒪(χ^11) depending on the particular sites involves and the geometry of the network. Block entanglement entropy is one of the key quantities that will be analysed in sec. <ref>, as one of the aims of this work is to achieve the correct entanglement structure of the wavefunction. However, the calculation of entanglement entropy is generally computationally difficult. Fortunately, the form of the tensor network alleviates this to some degree, as discussed in appendix <ref>, such that the entanglement entropy can be computed from tensor network states with small χ.




§ RESULTS


In order to benchmark the proposed method we focus on the disordered anti-ferromagnetic XX and Heisenberg models, where Δ_z = 0, 1 in Eq. (<ref>) respectively. We select coupling strengths 0 < J_i < 1 from the strong disorder distribution <cit.>, which has the probability density function 

    P(J) = 1/Δ_J J^-1 + Δ_J^-1,

where Δ_J≥ 0 is the disorder strength.
When Δ_J = 1 the distribution is uniform and when Δ_J→∞ the distribution matches that of the infinite randomness fixed point. 
For the remainder of this paper we will concentrate on the stronger disorder regime, when Δ_J≥ 1, as this tends to be more challenging for numerical tenchniques.
Throughout this section we compare the efficacy of the dMERA with various other numerical methods, including the tSDRG algorithm <cit.> and DMRG, which is performed using the ITensor libraries <cit.>. 
It may be possible to obtain more accurate DMRG results by, for example, changing the algorithm to take into consideration inhomogeneities of the system <cit.> or performing measurements at multiple χ and extrapolating to the χ→∞ limit.
Here we concentrate on a standard DMRG implementation without extrapolating so we can compare the finite χ performance of the algorithms.
Finally we note that, as the XX model can be mapped to a free fermion model via a Jordan-Wigner transformation <cit.>, we can obtain exact results to benchmark against, even for large system sizes. 





 §.§ The disordered XX model


We begin by applying a dMERA, built from composition of the unit cell depicted in Fig. <ref>(c), of bond dimension χ = 2 to find the ground state of the (disordered) XX model.
We observe that after optimization with any disorder strength, the unitaries w_L,R, when viewed as a matrix between incoming and outgoing indices, may be written as

    w_L,R = [         1         0         0         0;         0  cosθ_L,R -sinθ_L,R         0;         0  sinθ_L,R  cosθ_L,R         0;         0         0         0         1 ].

with θ_L,R a free variational parameter. 
The disentanglers u converge to a fixed matrix,

    u = [  1  0  0  0;  0  0 -1  0;  0 -1  0  0;  0  0  0 -1 ],

that does not contain any free variational parameters.
Therefore, in this setting, the χ=2 dMERA is fully specified by L - 3 parameters, with L the system size. The optimization of this network can also be simplified from the general algorithm introduced in Sect. <ref>, as discussed in appendix <ref>. 
Furthermore, one should notice that if all the angles are set to θ_L = θ_R = 0, then the w_L,R become identity matrices and the the random singlet phase of SDRG is recovered.
The reason that this parameterization arises is not clear to us, and a proof using analytic methods is beyond the scope of this paper.
In the absence of a proof, the numerical results in the rest of the section are performed using a full optimization, but the simplified update reproduces the results in all test cases.





Results comparing how the average ground state entanglement entropy ⟨ S_A⟩ scales with block size L_A for disorder strengths Δ_J = 1, 2, 4 are presented in Figure <ref>. Note that we average both over all possible blocks of size L_A and 50 different disorder configurations. The results demonstrate that dMERA, even when restricted to bond dimension χ=2, reproduces the ground state entanglement accurately over the range of disorder strengths considered. In comparison the tSDRG algorithm has a tendency to over estimate the entanglement entropy, particularly when the disorder is small, yet still reproduces the expected logarithmic scaling. 
For small disorder, Δ_J = 1, we see that DMRG requires bond dimension χ = 50 in order to rival the accuracy of the χ=2 dMERA, while for stronger disorder DMRG is substantially worse and the entanglement quickly saturates to a constant even for large bond dimension χ. 
Notice also that the results from standard SDRG become more accurate as the disorder increases, which is to be expected as the approach becomes exact in the limit of infinite disorder <cit.>.





The staggered two-point correlation functions ⟨⟨s⃗_x_1 . s⃗_x_2⟩⟩, computed from 50 site systems, are plotted in Fig. <ref>. These have also been averaged both over all positions and over 100 disorder configurations. For length scales smaller than half the system size, the exact results demonstrate the characteristic power law decay of the random singlet phase, with an exponent that depends on the disorder strength Δ_J. While the χ=2 dMERA does reproduce this power law decay, it tends to underestimate the correlation especially for weak disorder Δ_J=1. However, the results are still significantly more accurate than those of the standard SDRG approach or from χ=20 DMRG, where the correlations decay exponentially at larger distances.

That the χ=2 dMERA tends to underestimate correlators can be understood by analysing individual disorder realizations, where it can be seen that the problem occurs when the maximum coupling J_i is not significantly stronger than one of its neighbours, J_i-1 or J_i+1. In such instances the SDRG that is used to determine the structure of the dMERA is also less accurate; although optimization of the tensors allows the wavefunction to be improved the correlations still end up systematically smaller than they should be at larger separations. In disorder instances where the maximum coupling is much greater than its neighbours (J_i≫ J_i-1, J_i+1), which is more likely to occur at stronger disorder Δ_J, we find that the correlation functions are faithfully reproduced by the χ = 2 dMERA.




 §.§ The disordered Heisenberg model



We now address the disordered Heisenberg model, which has Δ_z = 1 in Eq. (<ref>). This model no longer corresponds to free fermions and thus cannot be diagonalized exactly as with the XX model. Also, a dMERA of bond dimension χ = 2 is no longer viable as the set of unitary gates that preserve the SU(2) symmetry would be trivial; thus we increase the bond dimension to χ = 4, 8. 








Figure <ref> shows the average entanglement ⟨ S_A⟩ for a 20 site system when averaged over 50 disorder realizations.
One can infer a rough upper bound to the expected results from tSDRG, which tends to overestimate the entanglement, and a lower bound from SDRG, which tends to underestimate the entanglement. It appears that the approaches based on energy minimization, DMRG and the dMERA, underestimate the entanglement for small bond dimension, approaching from below as the bond dimension is increased. For small disorder, Δ_J = 1 as depicted in Fig. <ref>(a), the dMERA with χ = 8 and 4 appears to be slightly less accurate than DMRG with χ = 50 and 20 respectively. However in the cases with stronger disorder, Δ_J = 2,4 as depicted in Fig. <ref>(b-c), the dMERA clearly reproduces the entanglement more accurately than DMRG, where the entanglement saturates to a constant for larger blocks. As expected, SDRG also becomes more accurate with stronger disorder, almost matching the entanglement from a χ=4 dMERA, although still falling short of the χ=8 case.





The two-point spin correlation functions are plotted in Fig. <ref> for a system of 50 sites, averaged over 200 disorder configurations. The dMERA is seen to reproduce the expected algebraic decay of correlations, producing a pattern similar to that from tSDRG but differing in magnitude at longer distances. However, it is not clear which approach is more accurate as both methods seem quite well converged in bond dimension χ.




§ SUMMARY AND CONCLUSIONS



We have proposed a generalization of entanglement renormalization for disordered systems, which yields a new class of efficiently contractible tensor network ansatz that we call disordered MERA. A key feature of the dMERA is that the geometry of the network can be adjusted to the specific problem under consideration, whereas previous tensor network ansatzes typically use a fixed network geometry only adjust the content of the tensors to the specific problem. This allows much of the entanglement structure of wavefunction to be captured before the variational optimization of the tensors even begins.

Conceptually, the benefit of adjusting the network geometry can easily be understood from an example such as the concentric singlet phase also known as the rainbow state <cit.>. This state, or something close to it, could easily be represented by a dMERA, but would require exponential growth in bond dimension with system size in order to be represented by a  conventional tensor network, such as an MPS or MERA, due to the linear growth of half-chain entanglement entropy.

In practice, the advantage of the this approach was seen in the benchmark results of Sect. <ref>. Here a dMERA of dimension χ=2 yielded remarkably accurate results for the XX model, significantly improving over DMRG calculations with χ=50. 
The dMERA also seemed to capture the entanglement structure in the disordered Heisenberg model, although a larger bond dimension was required to do so adequately. Accordingly, we expect this method may require further development and refinement before it is viable for more challenging models.

There are many possible extensions and modifications to our proposed approach. In the basic unit cell of the coarse-graining transformation we imposed that the truncated sites were either in a product state, as depicted in Fig. <ref>(b), or in a singlet state, as depicted in Fig. <ref>(c). In general this could be chosen as any quantum state; for instance a spin trio type state could be used to describe the baryonic states found in Ref. QuiHM15. It could also be possible to use a different metric to determine the geometry of a dMERA, rather than through use of the SDRG as was considered here, allowing the approach to be extended beyond states based on random singlet phases. For instance, one could try to use the mutual information to set the geometry, similar to the proposal in Ref. HyaGB17. Alternatively, it may be possible to design an algorithm that adjusts the network structure automatically during the variational optimization, which remains an interesting avenue for future work. 


    AMG would like to thank Norbert Schuch and Rudolf A. Römer for fruitful discussions.







§ ALGORITHMIC DETAILS





 §.§ The Hamiltonian


Before the update can begin, the Hamiltonian needs to have its energy shifted so that the entire spectrum is negative. 
This is important to allow the update process to target the ground state, the details of which are discussed in sec. <ref>.
The shift is performed by diagonalizing the two-site Hamiltonian operators and subtracting the identity multiplied by the largest eigenvalue λ_i,i+1^max from each

    h^shifted_i,i+1 = h_i,i+1 - λ_i,i+1^max( _2⊗_2).

The sum of the shifts are stored and added to the final energy after minimization to obtain the true ground state energy once more.




 §.§ The SDRG order


The order in which the coarse-graining blocks are combined to create the full network is given by ordinary SDRG <cit.>.
The algorithm coarse-grains in the order of strongest interaction, renormalizing neighbouring interactions via second order perturbation theory.
For more information see <cit.> and <cit.>.

The SDRG rules for the XXZ model are as follows:
Find the largest coupling J_i in the set of coupling strengths.
The strongest coupled pair form a singlet to first order. 
Second order perturbation theory is then used to find an effective coupling J̃ for the spins either side
    
    J̃ = J_i-1 J_i+1/(1 + Δ_z) J_i.

The location i is stored, then the strongest coupled spins are removed, coarse-graining the chain by two sites.
This is repeated until there are no remaining sites.
The result is a list of the locations of the singlets at each level of coarse-graining as shown in Fig. <ref>(a), which can then be used to set the geometry of the tensor network.




 §.§ The initial tensors


The order of coarse-graining and the maximum allowed χ is used to set the initial state of the tensors.
If χ = 2 the network is in the most simple form; all u and w tensors are 2 × 2 × 2 × 2 and unitary, and each coarse graining block is of the form of Fig. <ref>(c).
Setting these as ( _2⊗_2) means that the initial state is the random singlet state.
Another option is to generate random numbers to fill the tensors.

The top of the network takes a slightly different form due to the small system size and periodic boundaries. 
When the system consists of four sites it is only possible to put one coarse-graining block, then the remaining pair of sites will be connected by a singlet.
As shown in Fig. <ref>, due to the PBCs this second singlet can be considered to be part of the same coarse-graining operation.



The full tensor network for a system of 20 sites with χ = 2 is shown in Fig. <ref>(c).
Notice that this is the network obtained given the order set out in Fig. <ref>(a).

If χ > 2 then the network needs to be altered to allow the bond dimensions to grow.
By looking at Fig. <ref>(c) it is clear that if the tensors at the lower levels of the network have χ = 2, then all of the tensors above must do also.
Thus the coarse-graining block of Fig. <ref>(b) is used to increase the dimension of the effective sites.
Once the bond dimension is above 2 it can grow using the standard coarse-graining block of Fig. <ref>(c).
The dimensions and form of the tensors is set in the initialization stage and will then not change during the algorithm.

There are many choices that can be made relating to how to implement the bond dimension increase.
We choose to only use the χ-increase block [Fig. <ref>(b)] when no information is to be discarded, this is because this form has been shown to be unreliable at retaining long range interactions with lower bond dimensions.
The update is performed without the ρ's so that is is just a coarse-graining of the Hamiltonian.
As such, we only update it on the first sweep as the environment will not change as the algorithm progresses. 
For the rest of the blocks the dimensions of the top two legs of the unitary in the χ-increase block are swapped, e.g. if the bottom two legs have dimensions 2 and 4 then the top legs have 4 and 2.
This is to be consistent with the simplified results of the XX model, which show that the unitary has properties similar to a swap tensor. 
All of these choices have been seen to be more effective as well as being conceptually preferable.




 §.§ The descending superoperators





The two-site local reduced density matrices {ρ} contain all of the tensors that are above the tensor to be updated so that the process minimizes the energy of the network as a whole.
To create the set of ρ's, start at the top and contract down.
The tensors of the coarse-graining block make up the descending superoperators, mapping a density matrix on level Ł' to one at level Ł.
The highest level ρ's consist of the four ways of contracting the top block and its conjugate as shown in Fig. <ref>, where i is the position of maximum coupling from SDRG when the system is four sites.
Here and in the following diagrams the optimal contraction order is found using the netcon function <cit.> and the contraction itself performed using ncon <cit.>.

For the next level down there are six sites on the lattice and the ρ's are created by contracting the coarse-graining block to the correct ρ at the level above (labelled as ρ^' in the figures).
There are three sites (i-1, i and i+1) that are below ρ^'_i-1 shown in Figs. <ref>(b,c,d).
The sites either side of the main coarse-graining block, ρ_i-2 and ρ_i+2 are below ρ^'_i-2 and ρ^'_i respectively as shown in Figs. <ref>(a,e).
Any sites that are not affected by the coarse-graining block have ρ's copied down from the level above. 
This is repeated until the bottom of the network creating the necessary ρ's for the updates of all of the tensors.







 §.§ The update


The process is much the same as for standard MERA <cit.>.
The update of the unitary and isometries is quadratic in the sense that we want to find, for example, the optimal u and u^† at the same time.
Unfortunately there is no exact way to perform a general quadratic update, thus the problem is linearized.
Linearization here means setting u^† as the hermitian conjugate of the current u, so the problem now is just to update u on its own.
Once the new tensor has been found the u^† is replaced with the new u and the process is repeated until convergence.

The update is the process of finding the tensor, here u, that minimizes the energy of the network as a whole.
We call the rest of the network the environment of u, Υ_u.
The fact that the Hamiltonian tensors contain all of the tensors below the coarse-graining block and the ρ's contain all of the network above makes the contraction of Υ_u simply the sum of the three diagrams containing the Hamiltonian components below the u, as shown in Fig. <ref>.
Notice that the linearization is accomplished by contracting u^† into the environment.



Then the optimization of u can be stated as

    E_min = min_u ∈𝒰Tr[ u Υ_u],

where 𝒰 is the set of unitary tensors.
Because the Hamiltonian is shifted, as mentioned in section <ref>, the minimum energy is the negative trace of the singular values of Υ_u <cit.>.
Thus u = -VU^† where U and V are found by performing a singular value decomposition (SVD) of the environment (Υ_u = USV^†)
This is shown in Fig. <ref>, where it is clear that this choice of u results in the energy being the negative trace of S.



The environment is then re-contracted with the new u^† and the process is repeated until convergence.
The update of the isometries is effectively the same but with different environment diagrams to be contracted.
The environment of w_L is given in Fig. <ref> and the environment of w_R constructed similarly.
The main difference to that of u is the i - 2 term, which does not contribute to the environment of u.




The top of the network is again slightly different as there are no ρ's, just the singlets that make up the top, as shown in Fig. <ref>.
This is updated by SVD in the same way as the rest of the tensors.
Diagrams for the environments of w_L and w_R are omitted but are straightforward to recreate.







 §.§ The ascending superoperators


Once the tensors of a block have been updated the Hamiltonian components are raised to the next level of coarse-graining using ascending superoperators <cit.>.
This is accomplished by applying the block to each Hamiltonian tensor individually.
Because the centre three components get mapped onto the same site of the coarse-grained lattice, the new tensor for that site is the sum of these three terms.
The tensor network diagrams for this raising operation are shown in Fig. <ref>.








§ OBSERVABLES





 §.§ Single- and two-site operators


Expectation values are calculated by applying the coarse-graining blocks on operators.
Because the block maps two-site operators on one level to two-site operators on the level above, calculating the expectation values means simply applying the raising operators on the operator in the same way as the Hamiltonian operators during the update.
However due to the unitary or isometric properties of the tensors, usually only a fraction of the tensors need to be contracted and the rest cancel to identities.
The remaining tensors are said to be in the causal cone of the operator.
An example of this is given in Fig. <ref>, which shows an operator acting on sites 11 and 12 in the network from Fig. <ref>(c).
The causal cone is highlighted in blue, the remaining tensors are not involved in the expectation value.
This allows expectation values to be calculated in a highly efficient manner.




The diagrams for the raising operation are effectively the same as those in Fig. <ref> except the summation of terms in (a) is unnecessary as the expectation value will only include one of the three diagrams at each level.
The numerical result of the expectation value is found when the top is contracted, as shown in Fig. <ref>.
The leading order cost is 𝒪(χ^7), but can be 𝒪(χ^6) for operators at some sites in some network geometries.
Single-site operators can be expressed as two-site operators via tensor product with an identity. 
The expectation values of which can then be contracted using the same diagrams as before.
This is the most practical method as a single-site operator will be mapped to two sites by the coarse-graining block anyway.







 §.§ Two-point correlation functions





Correlation functions are calculated by applying the coarse-graining blocks on the two operators independently just as single- and two-site expectation values.
As before, only the tensors within the causal cone need to be contracted. 
However, at some point in the process the two sites will be joined by a single block.
This is the point at which the causal cones of the two operators join, as shown in Fig. <ref>.
The joining means that the calculation of two-point correlation functions can be significantly more costly than a single-site operator.
This additional expense was also observed for standard MERA and was one of the reasons for the development of ternary MERA as opposed to the original binary variant <cit.>.
The reason is that when the two operators meet the outcome is not necessarily another two-site operator.
For the ternary MERA case however, there are choices of sites where the cost is 𝒪(χ^8).
This is when the two sites are at the centre leg of the isometry, i.e. at separations r = 3^q with q = 1,2,3 …. 
Otherwise the cost grows beyond 𝒪(χ^8).
For dMERA the problem can not be so simply worked around.
When performing calculations for disordered systems, discarding individual correlations will result in errors in the averaging and therefore losing the characteristic power law of the random singlet phase.
It is therefore necessary to calculate correlations between all pairs of sites or discard all data for a given separation r if one is too costly to compute.

When the operators in dMERA meet the possible outcomes are two-, three- and four-site operators, as shown in the tensor network diagrams of Figs. <ref>(a-f).

The raising of two-site operators is the same as in sec. <ref>, three- and four-site operators are shown in Figs. <ref>(g-l) and <ref>(m-t) respectively.
The expectation value of the correlation function is finally determined by contraction with the top of the network as shown in Fig. <ref>.
In full, the cost of a calculating a two-point correlation function can be between 𝒪(χ^6) and 𝒪(χ^11) depending on the sites and structure of the network.
As mentioned in sec. <ref>, for algorithmic simplicity a single site-operator is encoded in a two-site operator by tensor product with an identity.
This suggests that there are two ways of encoding the each operator: O_i⊗_i+1 and _i-1⊗ O_i.
Due to the fact that the network is inhomogeneous, this choice of encoding can affect the cost of contracting the correlation function. 
It is therefore important to choose the operator encoding that minimizes the cost.










 §.§ Entanglement entropy





The entanglement entropy S_A of a region A given by

    S_A = -Trρ_Alog_2ρ_A,

where ρ_A is the reduced density matrix, obtained by tracing over region B.
The contraction of ρ_A in general is inefficient as the computational complexity grows exponentially with system size.
The construction of the dMERA wavefunction makes the problem significantly easier as we only need to contract the tensors that join blocks A and B <cit.>, as illustrated by Fig. <ref>. 
This makes it possible to calculate S_A when χ and L are relatively small.
We note that in all cases the entanglement is still captured by the tensor network and it is just the calculation of the entropy that is problematic.




§ ENERGY MINIMIZATION OF THE XX MODEL


As discussed in sec. <ref>, dMERA for the χ = 2 XX model has a simplified form where the u tensors have no variables and the w_L,R have one angle each.
The number of parameters in this tensor network is then L - 3.
This number comes from the fact that there are L/2 - 1 coarse-graining blocks as the top one is trivial, each containing two w tensors.
The final pair of angles are also correlated; the difference between the two being the only free parameter.
The update process can be performed by simply differentiating to find the minimum energy.
For example, we require the minimum energy

    E_min = min_θtr[ Υ w(θ) ].

Writing the environment Υ, and w as matrices as in Eq. (<ref>), the energy can be written as

    E    = tr[ Υ w(θ) ] 
       = Υ_1,1 + Υ_4,4 + ( Υ_2,2 + Υ_3,3) cosθ
            + ( Υ_2,3 - Υ_3,2) sinθ,















where Υ_i,j are elements of the environment.
Differentiation shows that the stationary points are at

    tanθ = Υ_2,3 - Υ_3,2/Υ_2,2 + Υ_3,3,

and the second derivative is used to find if it is a maximum or minimum.
The environment is then rebuilt with the new w and the process is repeated until convergence.



45
fxundefined [1]
 ifx#1

fnum [1]
 #1firstoftwo
 secondoftwo
 
fx [1]
 #1firstoftwo
 secondoftwo
 





noop [0]secondoftwo

ref[1]@startlink#1@href
href[1]#1@endlink
anitize@url [0]`
12`$12`&12`#12`1̂2`_12`%12
startlink[1]
endlink[0]

rl [1]href #1












@bib@innerbibempty

[Anderson(1958)]And58
  author author P. W. Anderson, @noop  journal journal
  Phys. Rev. volume 109, pages 1492
  (year 1958)NoStop
[Basko et al.(2006)Basko,
  Aleiner, and Altshuler]BasAA06
  author author D. Basko, author I. Aleiner,   and author B. Altshuler, http://dx.doi.org/10.1016/j.aop.2005.11.014 journal
  journal Ann. Phys. volume 321,  pages 1126  (year 2006)NoStop
[Nandkishore and Huse(2015)]NanH15
  author author R. Nandkishore and author D. A. Huse, 10.1146/annurev-conmatphys-031214-014726
  journal journal Annu. Rev. Condens. Matter
  Phys. volume 6, pages 15 (year 2015)NoStop
[Altman and Vosk(2015)]AltV15
  author author E. Altman and author R. Vosk,  10.1146/annurev-conmatphys-031214-014701 journal journal Annu. Rev. Condens. Matter Phys. volume 6, pages 383 (year
  2015)NoStop
[Abanin and Papić(2017)]AbaP17
  author author D. A. Abanin and author Z. Papić, 10.1002/andp.201700169 journal
  journal Ann. Phys. (Berlin) volume
  529, pages 1700169 (year 2017)NoStop
[White(1992)]Whi92
  author author S. R. White, 10.1103/PhysRevLett.69.2863 journal journal Phys. Rev. Lett. volume 69, pages 2863 (year
  1992)NoStop
[Östlund and Rommer(1995)]OstR95
  author author S. Östlund and author S. Rommer, @noop  journal journal Phys.
  Rev. Lett. volume 75, pages 3537
  (year 1995)NoStop
[Schollwöck(2005)]Sch05
  author author U. Schollwöck, 10.1103/RevModPhys.77.259 journal journal Rev. Mod. Phys. volume
  77, pages 259 (year 2005)NoStop
[Schollwöck(2011)]Sch11
  author author U. Schollwöck, 10.1016/j.aop.2010.09.012 journal journal Ann. Phys. volume
  326, pages 96 (year 2011)NoStop
[Juozapavi ččius  et al.(1997)Juozapavi ččius, Caprara, and Rosengren]JuoCR97
  author author A. Juozapavi ččius, author
  S. Caprara,  and author
  A. Rosengren, 10.1103/PhysRevB.56.11097 journal journal Phys.
  Rev. B volume 56, pages 11097
  (year 1997)NoStop
[Juozapavi ččius  et al.(1999)Juozapavi ččius, Urba, Caprara, and  Rosengren]JuoUCR99
  author author A. Juozapavi ččius, author
  L. Urba, author S. Caprara,  and author A. Rosengren, 10.1103/PhysRevB.60.14771
  journal journal Phys. Rev. B volume 60, pages 14771 (year
  1999)NoStop
[Goldsborough and Römer(2014)]GolR14
  author author A. M. Goldsborough and author R. A. Römer, 10.1103/PhysRevB.89.214203 journal journal Phys. Rev. B volume
  89, pages 214203 (year 2014)NoStop
[Ruggiero et al.(2016)Ruggiero, Alba, and Calabrese]RugAC16
  author author P. Ruggiero, author V. Alba,   and author P. Calabrese, 10.1103/PhysRevB.94.035152 journal journal Phys. Rev. B volume 94, pages 035152 (year 2016)NoStop
[Ma et al.(1979)Ma,
  Dasgupta, and Hu]MaDH79
  author author S.-k. Ma, author C. Dasgupta,  and  author C.-k. Hu, http://link.aps.org/doi/10.1103/PhysRevLett.43.1434 journal
  journal Phys. Rev. Lett. volume 43,  pages 1434 (year 1979)NoStop
[Dasgupta and Ma(1980)]DasM80
  author author C. Dasgupta and author S.-k. Ma, 10.1103/PhysRevB.22.1305 journal
  journal Phys. Rev. B volume 22,  pages 1305 (year 1980)NoStop
[Lin et al.(2017)Lin,
  Kao, Chen, and Lin]LinKCL17
  author author Y.-P. Lin, author Y.-J. Kao,
  author P. Chen,  and author Y.-C. Lin, 10.1103/PhysRevB.96.064427 journal journal Phys.
  Rev. B volume 96, pages 064427
  (year 2017)NoStop
[Vidal(2007)]Vid07
  author author G. Vidal, http://link.aps.org/doi/10.1103/PhysRevLett.99.220405
  journal journal Phys. Rev. Lett. volume 99, pages 220405 (year
  2007)NoStop
[Vidal(2008a)]Vid08_2
  author author G. Vidal, 10.1103/PhysRevLett.101.110501 journal journal Phys. Rev. Lett. volume 101, pages 110501 (year
  2008a)NoStop
[Evenbly and Vidal(2009)]EveV09
  author author G. Evenbly and author G. Vidal, 10.1103/PhysRevB.79.144108 journal
  journal Phys. Rev. B volume 79,  pages 144108 (year 2009)NoStop
[Pfeifer et al.(2009)Pfeifer, Evenbly, and Vidal]PfeEV09
  author author R. N. C.  Pfeifer, author G. Evenbly,  and author G. Vidal, 10.1103/PhysRevA.79.040301 journal
  journal Phys. Rev. A volume 79,  pages 040301 (year 2009)NoStop
[Evenbly et al.(2010)Evenbly, Pfeifer, Picó, Iblisdir, Tagliacozzo, McCulloch, and  Vidal]EvePPI10
  author author G. Evenbly, author R. N. C. Pfeifer, author V. Picó,
  author S. Iblisdir, author L. Tagliacozzo, author
  I. P. McCulloch,  and author
  G. Vidal, 10.1103/PhysRevB.82.161107 journal journal Phys.
  Rev. B volume 82, pages 161107
  (year 2010)NoStop
[Evenbly and Vidal(2014)]EveV14
  author author G. Evenbly and author G. Vidal, 10.1007/s10955-014-0983-1 journal
  journal J. Stat. Phys. volume 157,  pages 931 (year 2014)NoStop
[Fisher(1992)]Fis92
  author author D. S. Fisher, 10.1103/PhysRevLett.69.534 journal journal Phys. Rev. Lett. volume 69, pages 534 (year 1992)NoStop
[Fisher(1995)]Fis95
  author author D. S. Fisher, 10.1103/PhysRevB.51.6411 journal
  journal Phys. Rev. B volume 51,  pages 6411 (year 1995)NoStop
[Westerberg et al.(1995)Westerberg, Furusaki, Sigrist, and  Lee]WesFSL95
  author author E. Westerberg, author A. Furusaki, author M. Sigrist,   and author P. A. Lee,  @noop  journal journal Phys. Rev. Lett.  volume 75, pages 4302 (year
  1995)NoStop
[Altman et al.(2004)Altman,
  Kafri, Polkovnikov, and Refael]AltKPR04
  author author E. Altman, author Y. Kafri,
  author A. Polkovnikov,  and  author G. Refael, 10.1103/PhysRevLett.93.150402 journal journal Phys. Rev. Lett. volume 93, pages 150402 (year 2004)NoStop
[Iglói and Monthus(2005)]IglM05
  author author F. Iglói and author C. Monthus, @noop  journal journal Phys.
  Rep. volume 412, pages 277 (year 2005)NoStop
[Refael and Moore(2004)]RefM04
  author author G. Refael and author J. E. Moore, 10.1103/PhysRevLett.93.260602 journal journal Phys. Rev. Lett. volume 93, pages 260602 (year
  2004)NoStop
[Refael and Moore(2007)]RefM07
  author author G. Refael and author J. E. Moore, 10.1103/PhysRevB.76.024419 journal
  journal Phys. Rev. B volume 76,  pages 024419 (year 2007)NoStop
[Hikihara et al.(1999)Hikihara, Furusaki, and Sigrist]HikFS99
  author author T. Hikihara, author A. Furusaki,
   and author M. Sigrist, 10.1103/PhysRevB.60.12116 journal journal Phys. Rev. B volume 60, pages 12116 (year 1999)NoStop
[Wilson(1975)]Wil75
  author author K. G. Wilson, http://link.aps.org/doi/10.1103/RevModPhys.47.773
  journal journal Rev. Mod. Phys volume 47, pages 773 (year
  1975)NoStop
[Evenbly and Vidal(2013)]EveV13
  author author G. Evenbly and author G. Vidal, in @noop  booktitle Strongly
  Correlated Systems: Numerical Methods, series and number
  Springer Series in Solid-State Sciences, editor edited by  editor A. Avella and editor F. Mancini (publisher
  Springer Berlin Heidelberg, year 2013) Chap. chapter 4NoStop
[Evenbly and Vidal(2011)]EveV11
  author author G. Evenbly and author G. Vidal, 10.1007/s10955-011-0237-4 journal
  journal J. Stat. Phys. volume 145,  pages 891 (year 2011)NoStop
[Laflorencie(2005)]Laf05
  author author N. Laflorencie, 10.1103/PhysRevB.72.140408 journal journal Phys. Rev. B volume
  72, pages 140408 (year 2005)NoStop
[Ite()]Itensor023
  http://itensor.org/ title ITensor
  library, note Version: 0.2.3, http://itensor.orgNoStop
[Henelius and Girvin(1998)]HenG98
  author author P. Henelius and author S. M. Girvin, 10.1103/PhysRevB.57.11457 journal
  journal Phys. Rev. B volume 57,  pages 11457 (year 1998)NoStop
[Fisher(1994)]Fis94
  author author D. S. Fisher, 10.1103/PhysRevB.50.3799 journal
  journal Phys. Rev. B volume 50,  pages 3799 (year 1994)NoStop
[Vitagliano et al.(2010)Vitagliano, Riera, and Latorre]VitRL10
  author author G. Vitagliano, author A. Riera,   and author J. I. Latorre,  http://stacks.iop.org/1367-2630/12/i=11/a=113049 journal journal New J. Phys. volume
  12, pages 113049 (year 2010)NoStop
[Ramírez et al.(2014)Ramírez, Rodríguez-Laguna, and Sierra]RamRS14
  author author G. Ramírez, author J. Rodríguez-Laguna,  and author G. Sierra, http://stacks.iop.org/1742-5468/2014/i=10/a=P10004 journal
  journal J. Stat. Mech. Theor. Exp. volume 2014, pages P10004 (year
  2014)NoStop
[Quito et al.(2015)Quito,
  Hoyos, and Miranda]QuiHM15
  author author V. L. Quito, author J. A. Hoyos,   and author E. Miranda, http://link.aps.org/doi/10.1103/PhysRevLett.115.167201 journal journal Phys. Rev. Lett. volume 115, pages 167201 (year
  2015)NoStop
[Hyatt et al.(2017)Hyatt,
  Garrison, and Bauer]HyaGB17
  author author K. Hyatt, author J. R. Garrison,
   and author B. Bauer, https://arxiv.org/abs/1704.01974 journal journal
  arXiv:1704.01974 [cond-mat.str-el]  (year 2017)NoStop
[Pfeifer et al.(2014a)Pfeifer, Haegeman, and Verstraete]PfeHV14
  author author R. N. C.  Pfeifer, author J. Haegeman,  and author F. Verstraete, 10.1103/PhysRevE.90.033315 journal journal Phys. Rev. E volume
  90, pages 033315 (year
  2014a)NoStop
[Pfeifer et al.(2014b)Pfeifer, Evenbly,
  Singh, and Vidal]PfeESV14
  author author R. N. C.  Pfeifer, author G. Evenbly, author S. Singh,   and author G. Vidal, https://arxiv.org/abs/1402.0939 journal journal
  arXiv:1402.0939 [physics.comp-ph]  (year
  2014b)NoStop
[Vidal(2008b)]Vid08
  author author G. Vidal, @noop  journal journal
  arXiv:0707.1454v2  (year 2008b)NoStop
[Tagliacozzo et al.(2009)Tagliacozzo, Evenbly, and Vidal]TagEV09
  author author L. Tagliacozzo, author G. Evenbly,  and author G. Vidal, 10.1103/PhysRevB.80.235127 journal
  journal Phys. Rev. B volume 80,  pages 235127 (year 2009)NoStop




