Entanglement Continuous Unitary Transformations
Abstract
Continuous unitary transformations are a powerful tool to extract valuable information out of quantum manybody Hamiltonians, in which the socalled flow equation transforms the Hamiltonian to a diagonal or blockdiagonal form in second quantization. Yet, one of their main challenges is how to approximate the infinitelymany coupled differential equations that are produced throughout this flow. Here we show that tensor networks offer a natural and nonperturbative truncation scheme in terms of entanglement. The corresponding scheme is called “entanglementCUT” or eCUT. It can be used to extract the lowenergy physics of quantum manybody Hamiltonians, including quasiparticle energy gaps. We provide the general idea behind eCUT and explain its implementation for finite 1d systems using the formalism of matrix product operators. We also present proofofprinciple results for the spin1/2 1d quantum Ising model and the 3state quantum Potts model in a transverse field. EntanglementCUTs can also be generalized to higher dimensions and to the thermodynamic limit.
Introduction. The study of strongly correlated systems entails some of the most important challenges in modern physics. As an example, it is still not fully clear whether the fermionic Hubbard model captures high superconductivity in cuprates [1, 2], or whether some frustrated Heisenberg antiferromagnets can stabilize a quantum spin liquid ground state or not [3]. It is mostly because of challenges like these that many numerical simulation methods were proposed over the years. Still, each one of these methods comes with its own limitations. For instance, exact diagonalization is restricted to systems of relatively small size, and quantum Monte Carlo simulations are hampered by the infamous signproblem [4].
In this context, continuous unitary transformations (CUTs) have become one poweful tool to assess quantum manybody systems [5]. In CUTbased methods the Hamiltonian is continuously transformed by some unitary operator that depends on a continuous parameter . This unitary is constructed in such a way that, when , the effective Hamiltonian becomes diagonal or blockdiagonal in second quantization allowing to extract the important physical properties of complex quantum manybody systems more easily. Nowadays, there are different ways of constructing unitary operators satisfying this property [5, 6, 7, 8]. One of the limitations of this technique, however, is that the evolution of the matrix elements of the Hamiltonian is described by a system of coupled differential equations, and the number of these equations increases with the flow in , quickly becoming intractable. As a consequence, several truncation schemes have been developed over the years which are for example based on perturbation theory [7, 9], realspace approaches using nonperturbative linkedcluster expansions [10] or the spatial extension of operators [11], or in momentum space [12] using the scaling dimension of operators [13].
Another family of methods that has become quite popular in recent years is based on tensor networks (TNs) [14]. Here the wavefunction of the quantum manybody system is represented by a set of interconnected tensors that make explicit the natural pattern of entanglement and correlations in the system. Famous examples of TN methods include the density matrix renormalization group (DMRG) [15] and the timeevolving block decimation (TEBD) [16] for 1d systems, as well as algorithms based on projected entangled pair states (PEPS) [17] for 2d systems, and the multiscale entanglement renormalization ansatz (MERA) [18] for critical systems. Using these methods it is possible to approximate with good accuracy lowenergy as well as timedependent properties. TN algorithms can be applied to finite and infinitesize systems [19, 20, 21], as well as to fermions [22]. Their only limitation is the amount and structure of entanglement in the wavefunction to be described.
Seeing the above in perspective, one might wonder whether the CUT can be formulated in terms of a TN structure that can be used efficiently to truncate the flow nonperturbatively in terms of entanglement. The main point of this paper is that this is indeed the case. More precisely, we show that TNs offer a very natural approximation of the CUT equations, by truncating in the operatorentanglement content throughout the flow. We call this new scheme entanglementCUT (eCUT). This idea can in principle be applied to any dimensionality and system size (including infinite) as long as a faithful TN description of the flow can be set up.
For the sake of concreteness, here we provide an explicit prescription for eCUT applied to finite 1d systems, using the formalism of matrix product operators (MPOs) [14] and we give proofofprinciple results for the spin1/2 1d quantum Ising chain and 3state quantum Potts model in a transverse field and open boundary conditions up to sites. We show that this allows to determine with good accuracy the groundstate energy as well as the 1st excited state – amounting usually to the onequasiparticle (1QP) energy gap –. Moreover, we implement a ground state generator for the unitary transformation, which we shall also discuss as a particular case of a block generator.
CUT and the flow equation. Here we review briefly some basic properties of CUT that will be needed later (the interested reader is referred to, e.g., Ref.[5] and references therein for more information). Let us start by considering an Hamiltonian describing a quantum manybody system that we wish to diagonalize. It is wellknown [5] that this can be achieved by a unitary transformation that depends on a continuous parameter , i.e.,
(0) 
In this equation, . Ideally, the unitary transformation yields a diagonal effective Hamiltonian . In practice, though, one targets a blockdiagonal Hamiltonian. The flowing Hamiltonian Eq. (Entanglement Continuous Unitary Transformations) can be cast in terms of a generator , such that
(0) 
The above is called the flow equation. This generator can be written in terms of the unitary transformation throughout the flow as
(0) 
In the same way, the unitary operator can be written in terms of the generator as
(0) 
where means “ordering”. It is easy to check that the generator is antihermitian. This generator must also be chosen such that the flow drives the Hamiltonian towards some blockdiagonal form, from which the lowenergy properties should be easy to extract. In this paper we choose to work with what we call the blockgenerator, which we define through its matrix elements as
(0) 
In the above equation, are the matrix elements of the Hamiltonian in some relevant quasiparticle basis which is decided beforehand. In practice, the goodness of the choice of basis is checked a posteriori by the goodness of the results (intuitively one chooses a quasiparticle picture that is known to be exact in some limiting regime, e.g., low or large magnetic fields). Also, is the size of the (lowenergy) block being targeted, e.g., for one isolates “just” the ground state, while for one separates the ground state and the first excitation from all other states. For , this generator corresponds to the socalled “variational groundstate generator” [8], which decouples the ground state. For arbitrary , the block generator decouples the dimensional lowenergy subspace, see Appendix.
MPOs and operatorentanglement. Operators with an inner 1d structure admit a TN representation in terms of socalled MPOs, see Fig. (1) for a diagrammatic representation. These are the natural operators acting on the matrix product states (MPS), which are a family of wavefunctions suitable to account for the properties of 1d gapped Hamiltonians with local interactions [23]. At every site of an MPO there is a tensor with two physical indices of size each, accounting for the degrees of freedom of the local Hilbert space at every site (e.g., for spin1/2). Tensors in the MPO are also connected by bond indices with bond dimensions , which in principle can depend on the site . These bond dimensions measure the operatorentanglement content of the MPO ^{1}^{1}1More specifically, the singular values for the bipartition , if normalized such that , satisfy then the constraints , which are nothing but the MPO “operatorversion” of the constraints for the entanglement entropies of an MPS [14].. In practice, we simply call “bond dimension of an MPO” the parameter such that .
eCUT. The main idea of eCUT is to represent the operators for , and in Eqs. (Entanglement Continuous Unitary TransformationsEntanglement Continuous Unitary Transformations) by TNs. As the flow in proceeds, the operatorentanglement content needed to represent these operators may grow, and so will do the bond dimensions. Therefore, the bond indices need to be truncated throughout the flow using standard TN procedures. To achieve this, the flow is discretized in “flow steps” in (similar to the discretization of time in TEBD [16]), and the truncations are done at each discrete step. Generically, such truncations amount to keeping the degrees of freedom in the TN accounting for the majority of the operatorentanglement, which naturally (and indirectly) close the system of coupled differential equations being generated by CUT.
Let us be more specific with how the algorithm works for a finite 1d system. The core of the method is to consider the ordered integral in Eq. (Entanglement Continuous Unitary Transformations) and break it into smaller steps of size , that is,
(0) 
with and the unitary transformation at step being defined as
(0) 
In this aproximation, the error at every step is . In the above equation, is the generator at step , which in turn depends on , the Hamiltonian at step , through, e.g., Eq. (Entanglement Continuous Unitary Transformations). In what follows we show how to approximate and at every step using MPOs for a system in 1d. The procedure reads as follows:
1) Initialization: at step zero, , we write the original Hamiltonian as an MPO with bond dimension , as shown in Fig. (1.a). We also write the original generator as another MPO with bond dimension , see Fig. (1.a). Usually it is easy to get an exact analytical expression for , but if this were unknown, one could then even approximate its MPO representation numerically by using Eq. (Entanglement Continuous Unitary Transformations) and the MPO for . Set also and . In what follows, operators with prime will correspond to truncated MPO approximations to operators without prime.
2) Main loop: for , iterate the following steps:
(i) Approximate : evaluate using a Taylor expansion, e.g.,
(0) 
where we stopped at first order, but higher orders can also be considered. This equation is the addition of two MPOs: one for , and one for , see Fig. (1.b). This, in turn, can be approximated by an MPO of bond dimension for an operator , which is our approximation to .
(ii) Approximate : approximately evolve the Hamiltonian by a step , namely
(0) 
The TN diagrams representing this equation are in Fig. (1.c). The approximation on the right hand side means that the resulting tensor network for is approximated by an MPO of bond dimension for an operator .
(iii) Approximate : compute the generator assuming, e.g., the block generator in Eq. (Entanglement Continuous Unitary Transformations) (but others should also be possible),
(0)  
where the last approximation means that we are approximating the resulting operator by an MPO with bond dimension for the generator operator , as shown in Fig. (1.d).
The whole algorithm follows then by iterating the main loop, that is
(0) 
In this way we obtain a series of approximated Hamiltonians . For sufficiently large , and sufficientlyaccurate MPO approximations, the Hamiltonian will converge to a blockdiagonal form. For the case of the generator in Eq. (Entanglement Continuous Unitary Transformations), the lowenergy rdimensional block will be decoupled from the rest in the chosen quasiparticle basis, from which it is now immediate to extract the ground state energy, as well as (potentially) the lowenergy excitations. Other observables can be computed by flowing the corresponding operator in parallel with the Hamiltonian. In practice, approximations and truncations for MPOs are done using standard TN techniques for 1d finitesize systems, see the appendix material for more information. Notice that our method amounts to finding a quantum circuit that blockdiagonalizes the Hamiltonian which, when run in reverse starting from appropriate product states, will produce the eigenstates of the original Hamiltonian . Therefore, this method is based on an evolution in the Heisenberg picture, i.e., we evolve operators rather than states. This strategy is similar to the one used in some techniques to compute the socalled MERA [24] and, in the continuumflow limit, the socalled continuousMERA [25]. The overall leading computational cost of our algorithm is easily checked to be .
To get the 1st excited state one can consider different strategies. A possibility is to directly consider the block generator for . This is straightforward, but it involves larger bond dimensions in the algorithm, because multiple eigenstates are being simultaneously targeted in superposition. An initial numerical test showed that this approach works well for systems of moderate size, but not too large. To overcome this problem, another possibility is to use an “inverted spectrum technique”, flowing with Hamiltonian , and scanning to hit the desired excitation. Moreover, for the case of a translationallyinvariant system (e.g, with periodic boundary conditions) the 1st excited state is protected by the momentum quantum number and could therefore be targeted as the ground state in a given sector. Our plan is to use this strategy in forthcoming implementations of the method. For the purpose of this paper, however, we choose yet a different option to show proofofprinciple results for the excited state, namely, we use the same strategy as for the ground state but with the “shifted” Hamiltonian , with an energy shift larger than the 1QP gap, and an MPS approximation of the ground state of (which can be computed with, e.g., TEBD or DMRG). This approach has several advantages: first, the required bond dimensions are not too large since only one quantum state is targeted. And second, the stability is well controlled as in the case of the ground state.
Finally, a good practical criteria to determine the best point in the approximated flow is to choose the minimum of the socalled reduced offdiagonality (rod), , where are the Hamiltonian matrix elements and is the subspace that we wish to decouple. If truncations were exact, this quantity would tend to be exactly zero, and therefore its global minimum over the flow defines the optimal point at which we read the energy.
Results. To prove the validity of our method we present here benchmarking calculations for the 1d spin1/2 quantum Ising modelÊ and the 3state quantum Potts model in a transverse magnetic field and with open boundary conditions. Their Hamiltonians are given by
(0) 
with the th Pauli matrix at site , the number of sites, and the Potts matrices are given by
(0) 
In our simulations, we considered a chain of up to spins (though we could easily consider larger sizes), with flow step , Taylor expansions up to third order, and MPO bond dimensions up to . The quasiparticle basis are either the eigenbasis of (basis) or (basis) for Ising,Ê and of for Potts (basis). A priori, the basis is expected to produce good results for low magnetic fields, whereas the basis should provide better results at large fields. We run the flow for up to steps, and target the ground state of the system and its 1st excitation using the method mentioned above.
For quantum Ising, our results for the ground state and are shown in Fig. (2)(a,b). We can see how the basis works better for low fields and the basis for large fields, as expected. The error decreases rapidly with the Taylor order (seemingly exponentially fast), as well as with the bond dimension. Also, we observe that the spectrums of singular values in the MPO decompositions decay exponentially fast (even close to criticality), which validates the precision of the truncation. Moreover, in Fig. (2)(c,d) we show proofofprinciple results for the 1st excited state and , using and previously computed with TEBD. Surprisingly, in this case the basis produces already reasonably good results for all field values, even better than the basis for large fields, indicating that the operator entanglement associated with the  and flows behave differently when it comes to targeting the excitation. Our results for Potts are shown in Fig. (3)(a,b), for , in the basis, and as compared to TEBD results. Again we see that the approach works remarkably well for large fields, since the quasiparticle basis becomes more exact as the field increases.
Conclusions and outlook. Here we have introduced a scheme to solve the flow equation in CUTs using TNs, based on truncating the operatorentanglement content generated throughout the flow. We have provided proofofprinciple results for 1d systems of size up to for the spin1/2 quantum Ising and 3Potts models in a transverse field. Our method allows to extract groundstate energies and, with small modifications, lowenergy excitations as well. We believe that this technique has lots of potential. For instance, it is a natural candidate to study excitations in 2d systems [26]. This technique will also be fruitful in the study of manybody localized phases [27] where, e.g., in 1d, Hamiltonians can be diagonalized by a finitedepth quantum circuit [28, 29]. Moreover, our method may be useful in the context of functionalRG as an alternative truncation scheme of the Wetterich RG flow equation [30]. Finally, by running the flow backwards one would be able to study how quasiparticles get “dressed” by quantum correlations in a quantum manybody system.
Discussions with J. Eisert, J. J. GarcíaRipoll, S. Kehrein, D. G. Olivares, F. Pollmann and M. Rizzi are acknowledged. R. O. and S. S. acknowledge financial support from JGU.
References
 [1] Ã J. Hubbard, Proc. Roy. Soc. (London), Ser. A 276, 238 (1963)
 [2] P. W. Anderson, Science 235, 1196 (1987).
 [3] V. Elser, Phys. Rev. Lett. 62, 2405 (1989); J. B. Marston and C. Zeng, J. Appl. Phys. 69, 5962 (1991); P. Nikolic and T. Senthil, Phys. Rev. B 68, 214415 (2003); R. R. P. Singh and D. A. Huse, Phys. Rev. B 76, 180407 (2007); ibid, 77, 144415 (2008); G. Evenbly and G. Vidal, Phys. Rev. Lett. 104, 187203 (2010); S. Yan, D. A. Huse and S. R. White, Science 332, 11731176 (2011); S. Depenbrock, I. P. McCulloch and U. Schollwoeck, Phys. Rev. Lett. 109, 067201 (2012); Z. Y. Xie, J. Chen, J. F. Yu, X. Kong, B. Normand and T. Xiang, Phys. Rev. X 4, 011025 (2014); D. Poilblanc and N. Schuch, Phys. Rev. B 87, 140407 (2013); D. Poilblanc, N. Schuch, D. PérezGarcía and J. I. Cirac, Phys. Rev. B 86, 014404 (2012); T. Picot, M. Ziegler, R. Orús and D. Poilblanc, Phys. Rev. B 93, 060407 (2016).
 [4] M. Troyer and U.Jens Wiese, Phys. Rev. Lett. 94 170201 (2005).
 [5] F. Wegner, Ann. Phys. (Leipzig) 3, 77 (1994); S. D. Glazek and K. G. Wilson, Physical Review D 49, 4214 (1994); S. D. Glazek and K. G. Wilson, Physical Review D 48, 5863 (1993); S. Kehrein, The Flow Equation Approach to ManyParticle Systems (Springer Science & Business Media, 2006).
 [6] M. Toda, âTheory of nonlinear latticesâ, Springer Verlag Berlin (1989); A. Mielke, Eur. Phys. J. B 5, 605 (1998); S.R. White, J. Chem. Phys. 117, 7472 (2002).
 [7] C. Knetter and G. S. Uhrig, Eur. Phys. J. B 13, 209 (2000).
 [8] C. M. Dawson, J. Eisert, and T. J. Osborne, Phys. Rev. Lett. 100, 130501 (2008).
 [9] J. Stein, J. Stat. Phys. 88, 487 (1997); C. Knetter, K. P. Schmidt and G. S. Uhrig, J. Phys. A 36, 7889 (2003); A. Hackl and S. Kehrein, Phys. Rev. B 78, 092303 (2008); H. Krull, N.A. Drescher, and G.S. Uhrig, Phys. Rev. B 86, 125113 (2012).
 [10] H.Y. Yang and K. P. Schmidt, Eur. Phys. Lett. 94, 17004 (2011); K. Coester, S. Clever, F. Herbst, S. Capponi, and K. P. Schmidt, Eur. Phys. Lett. 110, 20006 (2015).
 [11] T. Fischer, S. Duffe, and G.S. Uhrig, New J. Phys. 12, 033048 (2010); N.A. Drescher, T. Fischer, and G.S. Uhrig, Eur. Phys. J. B 79, 225 (2011).
 [12] C.P. Heidbrink and G.S. Uhrig, Phys. Rev. Lett. 88, 14601 (2002).
 [13] S. Kehrein, Phys. Rev. Lett. 83, 4914 (1999); S. Kehrein, Nucl. Phys. B 592, 512 (2001); M. Powalski, G.S. Uhrig, and K.P. Schmidt, Phys. Rev. Lett. 115, 207202 (2015).
 [14] R. Orús, Annals of Physics 349 117158 (2014); J. Eisert, Modelling and Sim. 3, 520 (2013); N. Schuch, QIP, Lecture Notes of the 44th IFF Spring School 2013; J. I Cirac and F. Verstraete, J. Phys. A: Math. Theor. 42, 504004 (2009); F. Verstraete, J. I. Cirac and V. Murg, Adv. Phys. 57, 143 (2008).
 [15] S. R. White, Phys. Rev. Lett. 69, 2863 (1992); S. R. White, Phys. Rev. B 48, 10345 (1993).
 [16] G. Vidal, Phys. Rev. Lett. 91, 147902 (2003); G. Vidal, Phys. Rev. Lett. 93, 040502 (2004).
 [17] F. Verstraete and J. I. Cirac, condmat/0407066; V. Murg, F. Verstraete and J. I. Cirac, Phys. Rev. A 75, 033605 (2007).
 [18] G. Vidal, Phys. Rev. Lett. 99, 220405 (2007); For an introduction, see, e.g., G. Vidal, chapter of the book Understanding Quantum Phase Transitions, edited by Lincoln D. Carr (Taylor & Francis, Boca Raton, 2010), arXiv:0912.1651v2.
 [19] G. Vidal, Phys. Rev. Lett. 91, 147902 (2003); G. Vidal, Phys. Rev. Lett. 93, 040502 (2004); G. Vidal, Phys. Rev. Lett. 98, 070201 (2007); R. Orús and G. Vidal, Phys. Rev. B 78, 155117 (2008).
 [20] J. Jordan, R. Orús, G. Vidal, F. Verstraete, and J. I. Cirac, Phys. Rev. Lett. 101, 250602 (2008); R. Orús and G. Vidal, Phys. Rev. B 80, 094403 (2009); H. N. Phien, J. A. Bengua, H. D. Tuan, P. Corboz, and R. Orús, Phys. Rev. B 92, 035142 (2015).
 [21] I. P. McCulloch, arXiv:0804.2509; G. M. Crosswhite, A. C. Doherty and G. Vidal, Phys. Rev. B 78, 035116 (2008).
 [22] P. Corboz, G. Evenbly, F. Verstraete and G. Vidal, Phys. Rev. A 81, 010303(R) (2010); P. Corboz and G. Vidal, Phys. Rev. B 80, 165129 (2009); P. Corboz, R. Orús, B. Bauer and G. Vidal, Phys. Rev. B 81, 165104 (2010); C. V. Kraus, N. Schuch, F. Verstraete and J. I. Cirac, Phys. Rev. A 81, 052338 (2010); I. Pizorn and F. Verstraete, Phys. Rev. B 81, 245110 (2010); Q.Q. Shi, S.Hao Li, J.Hui Zhao and H. Qiang Zhou, arXiv:0907.5520; C. Pineda, T. Barthel and J. Eisert, Phys. Rev. A 81, 050303(R) (2010); T. Barthel, C. Pineda and J. Eisert, Phys. Rev. A 80, 042333 (2009).
 [23] F. Verstraete and J. I. Cirac, Phys. Rev. B 73, 094423 (2006).
 [24] C. M. Dawson, J. Eisert, and T. J. Osborne, Phys. Rev. Lett. 100, 130501 (2008).
 [25] J. Haegeman, T. J. Osborne, H. Verschelde, and F. Verstrate, Phys. Rev. Lett. 110, 100402 (2013).
 [26] L. Vanderstraeten, M. Mariën, F. Verstraete and J. Haegeman, Phys. Rev. B 92, 201111 (2015).
 [27] R. Nandkishore, D. H. Huse, Annual Review of Condensed Matter Physics, 6: 1538 (2015); C. Monthus, J. Phys. A: Math. Theor. 49 305002 (2016).
 [28] B. Bauer and C. Nayak, J. Stat. Mech. (2013) P09005.
 [29] F. Pollmann, V. Khemani, J. I. Cirac and S. L. Sondhi, Phys. Rev. B 94, 041116 (2016).
 [30] J. Berges, N. Tetradis, and C. Wetterich, Phys. Rept. 363:223386 (2002).
 [31] T. Fischer, Description of quasiparticle decay by continuous unitary transformations, PhD thesis, Technische Universität Dortmund, 2011.
Appendix
Proof of validity for the block generator
We want to show that the block generator transforms the subspace of the block in the same way as the quasiparticle generator [6, 7]. To achieve this, we start with the special case , and follow the derivation in Ref.[31]. We denote the vacuum state with and, instead of transforming the Hamiltonian with , we transform the vacuum as (Schrödinger picture). The flow of the vacuum state is then given by
(0) 
We can now introduce an orthonormal basis and obtain
(0) 
with the matrix elements given by
(0) 
The above equations yield
(0) 
We shift the dependency to the vacuum state and obtain , with , and where we defined the dependent projector .
Now we want to generalize the derivation above for . We denote the states with , e.g, and . Consider the flow of any state , given by
(0) 
The matrix elements for the quasiparticle generator are given by , and for the block generator by
(0) 
Hence we have for the quasiparticle generator
(0) 
By adding and subtracting terms and shifting the dependence to the states one gets
(0) 
For the blockgenerator we get (after shifting the dependence to the states)
(0) 
which is a generalization of the equation found previously for . Notice that the transformation of the subspace with for both generators is independent of all the other states with , and it only depends on the initial Hamiltonian and the states themselves. The two generators only differ in the last term of Eq. ( ‣ Proof of validity for the block generator). But this term only includes the offdiagonal elements in the block, which vanish for the quasiparticle generator but do not transform for the block generator. Thus both generators transform the considered subspace in the same way, but the quasiparticle generator also diagonalizes the block. In our case, we find it convenient to diagonalize exactly the residual small block procuded by the flow.
Details on the numerical Tensor Network approach
Truncation. At every step in the algorithm, truncations need to be done in the MPO bond dimension in order to avoid an exponential growth in the number of parameters, which amounts to truncate in the operatorentanglement content. This could be achieved in different ways. For instance, one could do a variational optimization by sweeping throughout the MPO tensors as in DMRG [14, 15]. An alternative strategy, however, is to find the canonical form of the MPO and then truncate in the bond dimension, as in the TEBD algorithm [19]. In our 1d calculations we have seen that both approaches produce similar results, but the second is a bit more efficient.
Ê Stability.ÊWe have found several sources of instability in the code that need to be very taken into account for large sizes. First, the algorithm is very sensitive to the Taylor order, since errors from the nonunitarity of propagate fast. It is thus important to do a precise truncation in the Taylor series for , which in our case is sufficient at order three. The second source of error has to do with the overall norm of the MPOs. In TN algorithms based on the Schrödinger’s picture, where quantum states are evolved (e.g., TEBD), tensors can be normalized almost at will throughout the flow in order to keep them numerically wellconditioned, as long as the final observables are computed dividing by the norm of the whole quantum state. However, this is not the situation for eCUT, since it is based on the Heisenberg picture. To be more precise: multiplying the MPO tensors by constants in order to keep them numerically wellconditioned must be done such that the overall norm of the MPO does not change, because otherwise one introduces big errors in the Taylor expansion, as well as in the eigenenergies. This turns out to be quite delicate when combined with the canonical form of MPOs mentioned before [19], since tensors usually have very small components, while matrices have them very large (and therefore the overall norm of the Hamiltonian MPO is always , with the size of the system). In practice, we found that the following “tricks” produce wellbehaved tensors for large :

After every MPO truncation, normalize s and s as , with a given site. This normalization, which we call relative normalization to the left, produces wellconditioned tensors and preserves the overall MPO norm.

After computing a new MPO, contract s and s as , and operate with tensors until the next truncation (i.e., absorb s to the left).

Compute the canonical form by applying until convergence the identity gate to the MPOs [19], sequentially and always from left to right.

In the truncation, use an adaptive bond dimension, i.e., remove diagonal values of below some cutoff , being this cutoff proportional to the overall norm of , e.g., .
The combination of the above tricks helps greatly to have wellconditioned tensors throughout the flow while keeping the overall norm of the MPOs constant.