{'fraction_non_alphanumeric': 0.03594935169562376, 'fraction_numerical': 0.009621961785763144, 'mean_word_length': 4.280783432332638, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 14, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Improved numerical methods for infinite spin chains with long-range interactions

We present several improvements of the infinite matrix product state (iMPS) algorithm for finding ground states of one-dimensional quantum systems with long-range interactions. As a main new ingredient we introduce the superposed multi-optimization (SMO) method, which allows an efficient optimization of exponentially many MPS of different length at different sites all in one step. Hereby the algorithm becomes protected against position dependent effects as caused by spontaneously broken translational invariance. So far, these have been a major obstacle to convergence for the iMPS algorithm if no prior knowledge of the systems translational symmetry was accessible. Further, we investigate some more general methods to speed up calculations and improve convergence, which might be partially interesting in a much broader context, too. As a more special problem, we also look into translational invariant states close to an invariance braking phase transition and show how to avoid convergence into wrong local minima for such systems. Finally, we apply the new methods to polar bosons with long-range interactions. We calculate several detailed Devil's Staircases with the corresponding phase diagrams and investigate some supersolid properties.


I. INTRODUCTION
A numerical method for the simulation of large quantum systems needs to meet two requirements: i) an ansatz suitable for the problem in question, and ii) efficient algorithms to find the (at least nearly) optimal solution within the chosen ansatz. For one-dimensional quantum systems on lattices, the currently most powerful numerical tools are matrix product states (MPS) based algorithms including the density matrix renormalization group (DMRG) [1][2][3][4][5]. Their primary limitation is given by the amount of entanglement they can handle. Several extensions of MPS have been conceived to overcome this restriction as, (e.g. [6][7][8]), but for most practical applications MPS are still the first choice. This is mainly due to the performance of the underlying optimization routines. Although the general task to find ground states is known to be NP hard [9,10], commonly used algorithms seem to have no problem to attain optimal MPS solutions within computer precision for a plenitude of physical relevant systems. Nevertheless, for some physical systems of interest these algorithms still face severe difficulties.
In this paper, we treat such problematic cases given by ground states of infinite spin chains with long-range interactions. An increasing interest in reliable numerical methods for these states is, e.g., triggered by the excellent experimental control of ultracold gases and the possibility to realize systems with long-range dipole-dipole interactions such as Rydberg atoms or polar molecules [11][12][13]. Although these systems are of finite size, one is often interested in the thermodynamical limit (i.e., infinite systems) for a better insight.
Different strategies are known for the numerical study of ground states in the thermodynamic limit. One might try to extrapolate results from a series of increasingly large finite systems [14][15][16] or directly construct the infinite state itself. The latter is, e.g., done by the infinite time-evolving block decimation (iTEBD) algorithm [17], which is based on an explicit translational invariant ansatz. This ansatz is quite elegant for interactions which are restricted to nearest neighbors, while it gets impractical for long-range interactions.
A comfortable way to incorporate long-range interactions is to encode them into a matrix product operator (MPO) [18][19][20][21]. This concept can be integrated in an infinite matrix product state (iMPS) algorithm [22,23].
The basic idea of this iMPS algorithm is to obtain the ground state of an infinite system as the fixed point of a constantly growing finite state by inserting iteratively new sites into its middle until convergence is reached. A major disadvantage of this approach is that the algorithm generally fails to converge if the ground state has a nontrivial translational symmetry.
In this paper, various extensions to the basic iMPS algorithm are presented. As a central element, the superposed multi-optimization (SMO) is introduced (section III B), which provides a remedy for the just mentioned convergence problem. Here, the key idea is to join the optimization of exponentially many MPS in a superposition and solve it efficiently. Due to this superposed optimization, the effective overall problem becomes translational invariant again and poses no longer a hindrance for convergence.
We also introduce several improvements which work independently of the SMO method. As such, we present two different modifications of the MPS optimization routine: On the one hand, we use physical insight for systems close to translational invariance breaking phase transition to suggest a method to reduce the danger of being trapped into a local minimum of the energy (section III D). On the other hand, we provide a more technical discussion how to recycle information from previous optimizations to speed up calculations (section IV C). As a part of these considerations we provide a simple implementation of a Davidson like algorithm [24] based on information from previous optimizations (appendix D). This implementation is not bonded to the MPS framework and hence might be used in much broader context.
All algorithms are described in depth, such that readers who are willing to reproduce our results should find all the information needed for successful programming. As a consequence, readers who just like to understand the crucial ideas might find the amount of algorithmic details far beyond their interest. Having these two types of readers in mind, we partitioned the material according to its level of detail into different chapters, such that entire passages can be omitted.
In section II, we review the basic concepts and the iMPS algorithm as presented in [22]. Readers who feel save to skip this part find in Fig. 1 and Fig. 2 a pictorial description of all symbols used in this section. In Section III, the main new concepts of this paper are presented, above all the SMO method. Changes of the algorithm are kept to a minimum in here in contrary to section IV where several algorithmic improvements and their numerical realization are presented in detail. Each subsections of section IV contains an individual topic, which is presented in the first lines. These subsections can be skipped without danger of losing the ability to understand the rest of the paper. A slight exception might be section IV C, in which the concept of iterative eigenvector solver based on subspace projections is reviewed. Familiarity with this concept is assumed in section IV D, IV E and appendix B 1. In section V, we apply our new algorithm to a system of polar bosons with long-range interactions. Detailed calculations of Devil's Staircases and phase diagrams are shown and a supersolid like phase is investigated. Finally, the paper is complemented by an appendix, into which several details have been outsourced.

II. BASIC ALGORITHM
In this section we review fundamental concepts [1] and the iMPS algorithm as presented in reference [22]. For this algorithm to work not only the Hamiltonian but also the ground state have to be translational invariant. The extension to ground states with broken translational symmetry will be introduced in section III.

A. MPS and MPO
In this paper we deal with spin chains. The quantum state of a spin chain is determined by the inner degrees of freedom of its components |ψ = s1s2···sn A s1s2···sn · |s 1 ⊗ |s 2 ⊗ · · · ⊗ |s n . (1) Since the size of the tensor A s1s2···sn grows exponentially with the number of sites, a more economical representa-tion is needed. A matrix product state (MPS) [1,3,4] consists in the ansatz where we used the Einstein summation convention. For a general exact quantum state exponentially growing bond dimensions α i are needed, but even for infinite systems excellent approximations are possible with moderate bond dimensions if the ground state fulfills an area law for the entanglement entropy [25].
Many relevant Hamiltonians are represented by MPO with relative small bond dimensions. For our purposes, it is important to mention that this is often also true for Hamiltonians with long-range interaction terms (see, e.g., [21] [R] s n sn (5) Further more, all tensors H are independent of the total number of sites. As a consequence, the MPO of a translational invariant Hamiltonian for n sites can easily be augmented to n + 1 sites by just inserting another tensor H.
Apart from an efficient representation of quantum states and operators, MPS and MPO also provide an efficient way to calculate expectation values. For details, we refer to the literature (as, e.g., [1]).
Remark: In this paper we apply the term translational invariant also to finite systems (with open boundary conditions) and their Hamiltonians which are used in the iMPS algorithm to approach the infinite case.

B. Overview of the iMPS algorithm
Any algorithm which deals with infinite MPS could be addressed as iMPS algorithm. In this paper we use this term exclusively for algorithms of the type as presented in reference [22]. This algorithm aims at finding an MPS representation for the ground state of an infinite one-dimensional quantum system. In its plain version the iMPS algorithm takes translational invariance for granted such that all sites behave equally. Thus all we need to construct the entire state is a perfect description of one site and its entanglement features with its environment given by the rest of the system. This environment is dominated by nearby neighbor sites while the influence of sites far away can be neglected in any onedimensional system with an asymptotic decay of correlations faster than r −1 . Therefore the environment built up by an infinite system can be simulated with a finite system. Correspondingly the center site of a sufficiently large but finite system provides a good approximation for its counterpart in the infinite case.
The iMPS algorithm is built upon a finite system, which is iteratively enlarged by inserting new sites into its middle. Since we express quantum states by MPS each of these new sites is represented by an individual tensor A [n] (2). Before a new tensor A [n] is inserted it is optimized such that the resulting energy is minimized. Hereby all previously inserted tensors A [j<n] are left unchanged. Of course, these local optimizations of the new tensors A [n] are generally not sufficient to find the lowest energy state of the entire system. What we are supposed to get is a ground state approximation which might be bad at the outer edges but close to the center, it should become better with each new tensor inserted. This is all we need to obtain an adequate description of the center site's environment, since the influence of the outer sites fades away with distance, anyway. Therefore, we expect the environment of the center site to converge towards its infinite counterpart, and with that the new tensors A [n] inserted into each round should converge, too For ground states which violate translational invariance, it is no longer given that all sites behave equally. At this point, our intuitive argumentation breaks down and the algorithm generally fails to converge. We will study this case in section III. For the time being, we assume that the algorithm ends up with a converged tensor A [converged] . With the help of this one tensor, the entire iMPS can be constructed using the rules we will encounter in II C.

Long-range interactions
An important ingredient for a fast computer code is a clever book-keeping of the interaction terms, which al-lows us to save many calculations due to recycling. In the case of long-range interactions, this task becomes tricky, since the iMPS algorithm permanently splits the MPS and adds new sites. By this means, the distances between sites on the left and right halves change each time and with them all distance-dependent interaction terms. Nonetheless, for translational invariant Hamiltonians recycling can still be done in a well-arranged fashion by encoding the Hamiltonian into an MPO as in equation (5). Every time the MPS is enlarged by a new tensor, the MPO is enlarged by the standard tensor H µ l µr s s , too. This simple procedure automatically corrects all distance-dependent interaction terms.

C. Constructing the MPS
So far, we just mentioned that the iMPS algorithm constructs the MPS by constantly inserting new tensors in its middle. We will now specify on that. First, it is convenient to treat the MPS as divided into two halfs, left and right from the center. Each time the optimization procedure described in II D provides a new tensor A α l αr s , we have to decide into which half A α l αr s is to be absorbed. This is done under the following rules: 1. According to the half in which A α l αr s is to be absorbed, decompose it ( Fig. 1 Q is orthogonalized such that where the asterisk denotes the complex conjugation (see the upcoming equation (10) for details).

With each new tensor
A overwrite the λ of the tensor before.
Each A [n] is optimized in such a fashion that it compensates for the overwritten λ [n −1] . For the orthogonalization (9) we proceed as follows: In a first step, we write the tensor A α l αr s as matrix A m,a where m is a multi-index. If the tensor A α l αr s is to be absorbed into the left MPS half m = (s, α l ) elsewise m = (s, α r ). For the leftmost and rightmost tensors of the MPS, which have only two indices, m = s. The index a corresponds to the leftover index of A α l αr s , which is not in m. Next, the matrix A m,a is decomposed into an orthogonal part Q and a "rest" part λ. Different decompositions would fulfill this task, but for the working of the algorithm it is best to resort to a singular value decomposition Rewritten as tensors we end up with equation (8).
A consequent application of these rules yields an MPS built of orthogonalized tensors Q and one single matrix λ from the very last A [new] in the center ( Fig. 1 (ii)). (11) This MPS standard form is numerical robust and has an easily calculated norm ψ = ψ|ψ . To see this let us multiply Q (the left half of equation (11)) with its complex conjugated and analog for the right half. Thanks to equation (9) all Q · Q * -pairs turn into δ-functions and the norm of the MPS (11) equals the remaining λ (see Fig. 1 (iii)) which also equals the norm of the last inserted tensor A [new] . Thus

D. Tensor optimization
The iMPS algorithm is an iterative procedure. As described in section II B new tensors A [n] are constantly inserted into the MPS, which represents the finite state ψ. Each of these new tensors A [n] is optimized such that the energy E = ψ|H|ψ ψ|ψ is minimized. As we will discuss in more detail below, this can be written as • |A [n] is the vectorized form of the tensor s with the multi-index i = (α l , α r , s). • I [n] is the identity operation thanks to the orthogonalized standard form of the MPS, see (12).
Equation (14) is solved by setting |A [n] equal to the lowest eigenvector of H [n] .
Remark: We will repeatedly use the notation |T or T | for a vectorized tensor T .

Effective operator H
In II C we mentioned that it is convenient to treat the MPS as divided into two halves: left and right from the newest tensor A [n] . For the same reason we decompose and a right half R α r µrαr [n] , which are connected by a single MPO tensor H µ l µr s n sn (5) corresponding to the new site (see Fig. 2 (iii)).
Since the iMPS algorithm is an iterative procedure, L [n] and R [n] are built up iteratively, as well. Suppose we intend to absorb the A [n−1] of the previous optimization step into the left half. First, we use equation (8) to gain the orthogonal tensor Q αlα l   Fig. 1). The boxes indicate the left and right halfs L α l µ l α l (16) and R α r µr αr (17). (iii) Same object as above with contracted inner indices of L α l µ l α l and R α r µr αr . The object in the H-shaped box corresponds to the effective operator H (15).

E. Algorithm
After having presented the decisive ingredients of the iMPS algorithm, we like to emphasize the steps one actually has to perform on the computer. The algorithm consists of an initializing procedure (see II E 2) and an iteration loop, which is repeated until convergence is reached (7). The tensor A [converged] is all we need to calculate expectation values. We do not hold any copy of the MPS we are calculating. The only objects stored (in the purest version of the algorithm) are the actual versions of L α l µ l α l [n] (16), R α r µrαr [n] (17) and A α l αr [n] s .

Initialization
At the beginning we have to initialize the values of L α l µ l α l and R α r µrαr . This can be done with the help of an exact solution ψ for a small system best with an even number of sites n = 2k (we mostly used n = 8 or 10 for two-level sites). The state ψ is split into a left and right part, which allows to calculate L and R. In more details, note the following: 1. Write the Hamiltonian of the small system as matrix with multi-indices H (s 1 ···s n ),(s1···sn) and solve for ψ (s1···sn) .

Use a singular value decomposition
6. Use V analog to U to calculate R α r µrαr .
Comparing with (11) we find that

III. BROKEN TRANSLATIONAL INVARIANCE
In this section, we present some new concepts for the iMPS algorithm which arise from the need to deal with spontaneously broken translational invariance. A principal shortcoming of the basic iMPS algorithm is its failure to converge in such cases. To overcome this deficiency, we introduce the superposed multi-optimization (SMO) method in section III B. Once convergence is restored, we turn our attention in section III C to the question of how to obtain a specific solution out of the ground state manifold degenerate due to broken translational invariance. In addition, in section III D, we treat the special case of a non-degenerate ground state which is separated by a very small energy gap from a state that breaks translational invariance.

A. Preliminary considerations
Just one tensor A [converged] suffices to construct an entire iMPS. At first sight, one might therefore think that such an iMPS is only capable of describing states where all sites behave equally, which is no longer true for states which break translational invariance. But still, also these states can be handled. This is due to the construction rules presented in section II C. These result in an iMPS structure given by equation (11), where the matrix λ marks a special position (and with that breaks translational invariance) if it can not be commuted to its neighbor sites.
The real problem is to find A [converged] . In II B, we already mentioned that our argumentation in favor of the convergence A [n] → A [converged] is no longer valid in the case of broken translational invariance. The iMPS algorithm is grounded on local optimization and therefore it is vulnerable to locally altering states, as they appear on a physical level for states with broken translational invariance. In this case, simple local optimization will not result in a global optimal fixed point A [converged] .

Known solutions
When we write about the breakdown of translational invariance, we mean that the state is no longer invariant under the shift of one site. Still, the state can maintain invariance under the shift of k sites. If k is known, one can introduce new super sites where one super site encompasses k of the old sites. Now, the system is translational invariant with respect to the shift of one super site. This involves that we have to optimize tensors which represent the super sites. Due to the exponential increase of the physical dimension, this method is practical for very small k only. To avoid the scaling problem, Crosswhite [22] suggested to use an MPS ansatz for the super sites. In practice this means, we insert k old sites at once and optimize the corresponding tensors. We are not aware if this was ever tested successfully. Aside from that, one still needs a prior knowledge of the value k.
Alternatively, one can extend the standard MPS structure with auxiliary tensors [26,27], which allow us to introduce a symmetry breaking element for the price of a non-linear optimization. To our knowledge, this has not been tested for long-range interactions, so far.

B. Superposed multi-optimization (SMO)
Our solution of the convergence problem induced by locally altering states does not depend on any prior knowledge and stays within the standard MPS framework. The key idea is to wash out local dependency of the optimization by choosing each new tensor A [n] such that it minimizes the sum of the energy of exponentially many different MPS instead of just one. These MPS are different ground state approximations to the qualitatively same Hamiltonian applied to systems of different sizes. All necessary minimizations can be joined in a superposition and solved by one optimization. The time relevant steps stay the same as in the single MPS optimizing algorithm presented so far, so there is no noticeable loss in speed.
As explained in the following, after each optimization round, the number of MPS joined in the superposition increases by a factor of 4. Thus, in the nth round, the optimization (14) gets formally extended to The MPS representing the |ψ i are of different length and the position of the tensor A [n] is no longer in the center but varies from MPS to MPS. In the basic algorithm, the tensors A [n] experience quite individualized environments and the optimization adapts to these local circumstances. In the modified algorithm, A [n] faces exponentially many different environments averaging out local effects and emphasizing common global features. This enforces heavily the desired convergence (18) is based on 2 2(n−1) = 4 n−1 different MPS. These MPS are not hand picked, but indirectly generated by the algorithm. The only modification needed to create such a superposition concerns the left and right halves L α l µ l α l and R α r µrαr . We still use equation (16) and (17) to perform the iteration steps L

Modification of the algorithm
and R but afterwards we add the new and the old results to achieve superpositions Since this is a simple addition of two tensors, the structure and size of L and R stay the same and do not entail any computational complications.
In the basic algorithm we have to decide at each iteration step whether to absorb the tensor A [n−1] into the left half L α l µ l α l [n] (16) or into the right half R α r µrαr [n] (17). Only the half of choice is modified. Now, we symmetrize the algorithm and modify both halves in each step. With this modification, the operator sum H [n] is calculated with one single use of equation (15). As a further advantage of the symmetrization the modified iMPS algorithm can now take advantage of mirror symmetric Hamiltonians. As this subject is a bit off-topic, we refer the interested reader to Appendix F.
Since each of the tensors A [1] , . . . , A [n−1] is absorbed into the left half L and the right half R the longest MPS encoded in the operator superposition of H [n] contains 2(n − 1) tensors (where we neglect the contribution from the initialization routine explained in section II E 2 and the hole for the new tensor A [n] ). But in the general MPS each of the 2(n−1) tensors only appears with a probability of 50% due to the addition of the new and old L and R in equation (19). The possible tensor combinations give rise to the heralded 2 2(n−1) different MPS from which k 2(n−1) are of length k. More precisely of these MPS contain k [L] tensors left and k [R] tensors right from the hole for the new tensor A [n] .

Comments
The crucial observation is that the iteration steps to perform are always the same independent of the tensor position and the size of the MPS. Therefore all the different MPS can be optimized together combined in a superposition. The reader who is more familiar with finite MPS calculations might wonder about the complete loss of information concerning the single MPS in the superposition, which comes along with equation (19). We have to remind ourselves that the main objective of the iMPS algorithm is to get the tensor A [converged] , which suffices to construct the infinite MPS. The finite MPS are just tools to obtain this tensor. Once we have it the finite MPS are no longer needed.
Another interesting question is whether the operator sum H [n] in equation (18) constructed via equation (19) is really suitable for a variational ansatz. This question has two aspects.
1. Is every optimal iMPS solution of equation (18) also a minimum of the Hamiltionian?
2. Does the algorithm always converge to this optimal solution?
The intuitive argumentation given in section II B also suggests that the answer to the second question should be yes, but actually even for the well-established DMRG algorithm the answer has to be no, since otherwise NP hard problems could be solved. Nonetheless, it is a matter of fact that DMRG converges extremely well for most practical purposes. The same "practical proof" can be given for equation (18) as demonstrated by our applications (section V).
The first question can be answered more formally. Equation (18) represents the sum of exponentially many energy terms. The lowest conceivable value of this sum is reached, if each energy term takes its individual minimal value. If all individual energies are minimized, the obvious answer to the first question is yes. Hence, we have to ask, whether it is possible to minimize all individual energies at once, having in mind that all MPS involved are created indirectly via equation (19). With finite MPS this might only be possible up to a certain relative error. But for the limit of infinite MPS this relative error shrinks to zero and the problem is trivially solved by uniform iMPS, i.e., in the case where all tensors stem from the same A [converged] . Then, all iMPS appearing in equation (18) look alike and either none or all of them minimize their Hamiltonians. Even in the case of broken translational invariance, one can always find at least one translational invariant ground state, which can be written as uniform iMPS and hence optimizes all terms in the sum at once.
Next, we like to further inspect the numerical consequences of equation (19)   Although we never experienced any practical problems, the cases 1 and 4 are, at least from the theoretical point of view, a bit troublesome. In case 4 the tensor A [n−1] is inserted twice. But A [n−1] was never optimized for double insertion. Close to the end, when A [n] → A [converged] is almost achieved, this should pose no problem. Meanwhile, at an early stage the effect should be more severe. On the other hand, even in the basic algorithm, the MPS description is not perfect -especially not at the beginning.
Case 1 might seem trivial, since everything stays the same. Potential difficulties arise in the superposition with the other cases. According to equation (8) the tensors A are decomposed into Q and λ and only Q is absorbed. The matrix λ is overwritten with the next A (respectively Q). In the basic algorithm, this is easy to justify: The next A [n] can compensate for λ [n −1] . But, a perfect compensation can only be achieved for one λ, not for many of them. Here is the problem: In case 1, old λ [n−2] , λ [n−3], . . . of the previous steps are conserved, while in 2., 3. and 4. a new λ [n−1] comes into play. All λ have to be compensated for. The stronger the λ alter, the less adequate is their compensation. The variation of the λ can be reduced by enforcing A [n] −A [n−1] to be small.
Towards the end of the optimization A [n] − A [n−1] is small anyway. At an early stage one might have to resort more strongly to the convergence enforcing method we will discuss in section IV A 2.

C. Selecting a specific ground state
We have seen how to ensure convergence in the case of broken translational symmetry. But so far, we have no control to which of the degenerate ground states the algorithm converges. Some of these ground states might be more favorable for our purposes than others, and we now answer the question as to how to obtain them. Any further degeneration aside from broken translational invariance is excluded from this consideration.
Let us look at two fully converged MPS A and B where B be a representation of the wished-for ground state which fits our purposes best, while A stands for any ground state to which the algorithm actually has converged. According to equation (11), both MPS have the following structure In Appendix C, we show that it suffices to replace the matrix λα 0 α0 in A by the new matrix γα 0α0 to obtain an MPS which represents exactly the same physical state as B In other words, we do not need to take care to which ground state the algorithm converges, since after it has converged, we are able to transform the obtained solution easily into any other. We do not even have to know B, as long as we have a description such as, e.g., "the ground state with the highest expectation value for the operator X". All we have to do is a one-time optimization of the new matrix γα 0α0 under the desired side condition.

One tensor update versus multi tensor update
So far, we focused on uniform MPS which result from an algorithm that inserts one new tensor each round. Crosswhite [22] suggested that one might also insert a certain number of q tensors per round in the form of a small MPS. Although we are so far not aware of any successful practical applications of this ansatz, it is worth having a closer look. For the single-site algorithm to work in the presence of broken translational, invariance we introduced the SMO method, which washes out local variations and thereby fortifies convergence. Still, the SMO is a general method and could also be implemented in an q-site algorithm.
As we have seen in the section above, the single-site iMPS algorithm will come up with a solution that encodes all possible ground states. This abundance has its price. Given the situation that we know the periodicity q of the ground state of interest, we could use an iMPS algorithm which inserts q sites at once. This would enable us to find some specific lowly entangled ground states which could be expressed by a non-uniform MPS with a far smaller bond dimension. Since q translationally shifted copies of such a non-uniform MPS always allow us to construct a uniform MPS, the maximal gain in bond dimension is given by a factor q and the maximal difference in the entanglement entropy of the half chain is ∆S = log 2 (q). We could confirm this difference for the model studied in the applications (section V) varying the matrix λ (20) over the set of ground states, as described in the section above.
From the perspective of the needed bond dimension, the multi-tensor update is clearly superior to the-single tensor update for systems with a periodicity q > 1. Sill, in this paper we favor the single-tensor update, which needs no prior knowledge of the periodicity and results in a well-converging algorithm, which has proofed its reliability in practical tests.

Degenerate tensor
In the case of broken translational, invariance one can jump from one ground state solution to another just by changing the matrix λ, which is part of the bigger tensor A [n] (8 To suppress this effect, we can resort to the convergence enforcing method we will present in IV A 2. In addition, we will describe in III D and IV D a numerical method to enforce a translational invariant solution which eliminates the above-mentioned degeneration.

D. Translational invariant ground states and local minima
In this section, we consider possible convergence problems due to translational invariance breaking states which lie closely above the non-degenerate ground state level. In such cases, the infinite system still provides a translational invariant ground state, while for finite systems even small alterations of the energy spectrum due to boundary effects suffice to favor a ground state with broken translational invariance. Since the iMPS algorithm is based on growing finite systems, it might start out converging into a false minimum and get trapped there. Even if the algorithm escapes out of this trap later, it supposably costs many optimization rounds and significantly slows down convergence.
To avoid these problems we suggest to modify the algorithm such that it only converges to translational invariant states. This is no limitation: In the case of an unique ground state, the sole solution has to be translational invariant, anyway. If the ground state level is degenerate, one of the solutions is translational invariant and according to equation (21) we can still transform it into another type of solution after the algorithm has converged.
Whether the fully converged MPS A, is translational invariant or not depends on its matrix λ.
At this point we should be more precise and write λ [L] or λ [R] , depending on whether λ stems from a left or a right decomposition (8). Actually, as a consequence of decomposition (8) In this case λ can be commuted to any position and hence does no longer mark any specific site of the MPS. This is what we are aiming for. In order to end up with an A [converged] where λ [L] = λ [R] we alter the minimization routine which computes the tensors A such that solutions with small differences λ [L] − λ [R] i.e. big overlap λ [L] |λ [R] are preferred. In the long run this should accumulate to λ [L] = λ [R] .
As a first straight forward way we tried to extend the minimization (18) of A| H|A to with a suitable coupling parameter γ [λ] . This is no longer a simple to solve bilinear problem since one needs to perform the decomposition (10) to get λ [L] and λ [R] . To avoid this complication and restore bilinearity we tried to resort to the easily calculated approximationsλ [L] and λ [R] (E2) derived in the appendix E, but the results we obtained in this way were not very convincing.
In section IV D we introduce a less conventional approach which turned out to work far more satisfyingly for us. Instead of extending the minimization of A| H|A by a new term as suggested in equation (24), we alter the routines of the iterative eigenvector solver we use to solve it. The modus operandi of these solvers is reviewed in section IV C. Until after then, we suspend further explanations.

IV. ENHANCED ALGORITHM
The considerations of the last section were mainly conceptual. The only actual change of the algorithm we performed is given by equation (19), which incorporates the SMO method. In this section, we delve far more into numerical details and extend the algorithm by further routines to make it more efficient. A reader not interested in technical details of the algorithm might proceed directly to section V

A. Enforcing convergence
The goal of the iMPS algorithm is the global convergence A [n] → A [converged] . This property has to emerge over the long term, while it is not part of the evaluation system of the local minimization from which each A [n] is drawn. As a consequence, small local improvements might be purchased with strong fluctuating A [n] counteracting global convergence. In an unstable scenario of overcompensation, these fluctuations might even inflate in a fatal manner. To prevent this from happening we extend the algorithm by two methods. The first method (superposition method) aims at attenuating the influence of problematic A [n] on the ongoing calculations, while the second method (gain function method) directly modifies the optimization routine such that excessive variation of the A [n] are suppressed. Both methods are complementary and worked well together in our calculations.

Superposition method
The first method takes advantage of the fact that the and with that inhibit a chain of overcompensation. We remind the reader that A [n] is absorbed into L [n+1] (16) and R [n+1] (17) before equation (19) is used to build up superpositions. This latter equation is now replaced by The only new ingredient compared to equation (19) is the adjustable weight 1 ξ [n] > 0 calculated as with N = A [refer] [n+1] . The weights ξ [n] used in equation (27) are the same as in equations (25) and (26).
For equation (26) to work we still have to define the weighted average ∆A [n] . Various definitions are possible. As a heuristic choice, we picked the following one: Finally, we remark that we end up in a deadlock if ξ [n] = 0. To prevent this from happening we will introduce the parameter ∆ max in equation (30) of the upcoming subsection.

Gain function method
The idea of the gain function method is to manipulate the minimization procedure of A .
In the appendix B we show how to approximate γ efficiently such that where 0 < c [n] 1 and ∆ max are parameters of our choice. Limiting ∆A [γ] [n] by assigning e.g. ∆ max = 10 · ∆A [n] (28)  2. This formula was found heuristically and worked fine for us, although more adequate choices might exist.
When the iMPS algorithm finally approaches its end γ becomes very small and its effect might be overruled by numerical imprecision. To prevent this, we recommend defining a lower limit for γ above the limit of the numerical precision.

B. Energy overgrow
If the average energy per site of an infinite state does not equal zero, the total energy of the entire state is ±∞. Of course, we never have to deal with an infinite value since our numeric is restricted to finite systems. Nonetheless, a problem remains. In the long run, the numeric value of all the information encoded in the tensor H stays more or less the same except for the energy, which grows with each new site. The tensor H gets more and more ill conditioned since the numeric value of the energy overgrows other information and thereby reduces the achievable precision. To avoid this problem, we advise to subtract each iteration step the energy  (16) and (17). As shown in Appendix A, the MPO tensor H µ l µr s s has a slot which represents a local interaction term. To this local interaction we add −E [n] · I s s .

C. Minimization routine and information recycling
With an increasing number of rounds n the successive minimizations of the different A [n] | H [n] |A [n] become more and more similar, which opens the opportunity to speed up the minimization recycling information from preceding turns. In order to understand these ideas (and also those of section IV D, section IV E, and Appendix B 1), we have to review the principles of the iterative eigenvector solvers we use [28]. In the MPS context, these solvers come with the major advantage that H [n] never has to be constructed explicitly; it suffices to be able to assemble H [n] |A for any given |A . Further, we do not need to perform the minimization to its very end. For the algorithm to work, it suffices to perform a limited amount of iterations, such that the resulting |A might not be optimal but still significantly improved. Due to information recycling, these improvements accumulate, such that the optimal solution emerges in the long run.
Iterative eigenvector solvers are very well suited for the outer eigenvalue spectrum. Already with modest effort we can expect to find a good approximation |e 0 ≈ |E 0 for the lowest eigenvector of H [n] . The central idea is to project the problem defined on a huge space of dimension N onto a much smaller subspace of dimension k N and solve it there. For this to work, we have to build up iteratively a small set {|A 1 , . . . , |A k } of k orthonormal vectors which enables us to express the minimizing eigenvector |A We need to solve for a [min] i , which is obviously the minimizing eigenvector for the k × k matrix H As long as r is too big we have to extend the set {|A 1 , . . . , |A k } iteratively by a further vector |A k+1 . Any form of educated guessing for a suitable new |A k+1 is allowed. The Lanczos [29] and Arnoldi [30] algorithm use a different way of calculation but end up with with |A k+1 ⊥|A 1 j k by construction.
In contrast to the basic iMPS algorithm, which sets |A 1 equal to the lowest eigenvector |e [n−1] 0 = |A [n−1] of the last round [22], we choose are already good approximations for |E 0 (=|E [n] 0 ). In many cases, we could obtain a considerable speed up extending this idea to a few more than just the first vector of the set Further, we observe that all |A [refer] (27) are derived from the best eigenvectors of the previous rounds. As an additional extension we also tried to include the next best eigenvectors |e [n−1] j>0 of the last round The improvements we achieved in this way were relatively poor. A much more promising way to take advantage of the |e [n−1] j is to use them for an efficient approximation of the inverse operator D = e 0 · I − H [n] −1 (D6), which allows a handy implementation resembling the Davidson (or Jacobi-Davidson) [24] method. As a result the update equation (37) is replaced by the more appropriate ansatz (D7). Details are explained in the Appendix D.
At the end of this section, we like to caution the reader that the methods presented here might counteract the methods presented in section IV A. Global convergence and improved local minimization often go hand in hand, but not always. If the algorithms indicates to run unstable, one should consider to partially switch off the improvements just presented. This is likely to happen if |e 0 is degenerate. In this case, the recycled knowledge from the past strongly increases the probability that already a shallow optimization suffices to find alternative solutions, which might result in unwanted fluctuations as, e.g., described in III C 2. Usually, this problem is announced in advance. In the Davidson implementation, one should not resort to eigenvectors with eigenvalues too close to the best. Similarly, once the small set of recycled initial values |A [refer] [n−j] suffices to get a second best eigenvalue very close to the best, it might be wise to abandon this method and only use |A alone.

D. Enforcing translational invariant ground states
In this section, we demonstrate the algorithmic realization of the considerations put forward in section III D. There, we argued that it is beneficial to push the algorithm towards translational invariant iMPS solutions to avoid getting trapped in local minima. We further showed that translational invariance is assured if the decomposition (8) (23). This is what we are aiming for.
The approach we are about to present is not very intuitive. Therefore, we start our explanations with an intermediate step and introduce a less practical but easier to understand procedure which consists of the following steps and has to be performed with each new tensor A after it has been optimized: . Nonetheless, the changes in A might be too pronounced to be acceptable. To soften this approach one can ignore line 4. and just go through 1. to 3. once. After that we generally still have λ [L] = λ [R] but with a reduced distance λ [L] −λ [R] compared to the initial value. This is all we need to achieve λ [L] = λ [R] in the long run. But the new A is still likely not to qualify for the optimizing tensor we are looking for. Now, we come to the procedure we really use. Instead of symmetrizing the tensor A after its optimization we integrate the symmetrization into the optimization routine. As recapitulated in section IV C, the optimization routine expresses the vectorized tensor A [n] = |A as a linear combination of a small set of basis vectors |A i (34). The idea is to alter these basis vectors |A i such that we have a similar effect as the procedure above. At the stage of the optimization the Q The |A i are created iteratively. In each iteration step we first create a new |A i as we used to do (37), (D8) and then alter it. Therefor we introduce |Ā i defined as where we tensorized the vector |A i in lines 2 and 3. With that, we replace |A i by an orthonormal version of |Ā i : For a better understanding we insert the |Ā i in the linear combination (41). As shown in Appendix E, we get which mimics the effect of the procedure presented above. But in contrast to the procedure above the story does not end here. The important point to notice is that the algorithm can still adopt to the alteration (43) of the basis vectors |A i and come up with alternative solutions. More favorable weights a i than those used in equation (44) are presumably to be found. Even the |A i themselves are likely to be different since they are calculated iteratively according to the needs of the minimization. While there are still enough resources to compensate sufficiently for the negative effects of the enforced alteration, the positive effects should survive since the arguments in their favor are largely independent of the a i , |A i chosen by the optimization routine. Still, this alteration is a trade off, but we have good reasons to believe that we gain more than we sacrifice. For practical applications, we only need a few lines of code to implement equation (42), which is also easy to turn off for systems where it is not needed, i.e., when the unaltered algorithm shows no tendency to run the risk of being trapped in a local minimum. In such a case, the alteration is likely to slow down the algorithm slightly. For the applications tested by us the loss in performance was only marginal. On the other hand we also encountered many cases where the altered algorithm clearly outperformed the unaltered one, which was partially even unable to find the correct ground state within the observed run time.
Although we strongly recommend to implement the alteration (43) as presented, one could also use a compromise and only alter the first basis vector |A 1 = A

[refer] [n]
, which has already a strong impact on the outcome of the optimization. This reduced version does not come with the need to program a new eigenvector solver. Each solver which accepts an initial vector |A 1 will do. In any case, the gain function in equation (29) is understood to change accordingly to the alteration of . E. Length of the MPS After n optimization steps, even the longest MPS encoded in the superposition created by the SMO method does not surpass the length l = 2 · n + l 0 (where l 0 is the initial length (see II E 2)). For some systems with long-range correlations this might be too short unless n reaches some considerably high number, which would go along with an extended calculation time. To shorten this calculation time two methods might be of help:

Use a tensor A α l αr
[n] s with small bond dimension χ [small] until a certain MPS length is reached, then increase the bond dimension to its final value χ [big] .
2. Use fast Krylov subspace methods [28] to insert the same tensor many times (e.g. 10 5 ) into the MPS.
A simple and comfortable way to increase the bond dimension from χ [small] to χ [big] is to use an isometric and proceed as follows after A [n] has been optimized but still not been inserted into L [n] (16) and R [n] (17): Next, the new tensor and R β r µrβr [n] as usual but without the superposition building step (19) respectively (25). To avoid trapping into a local minimum one might also consider to add a small amount of noise to A α l αr [n] s before applying equation (46). A possible strategy for the small bond dimension χ [small] is to proceed until convergence has been reached A [n] → A [converged] , but before the bond dimension is increased, many more copies of A [converged] are inserted into the MPS without any further optimization. These insertions can be done in the standard fashion or generally much faster by projecting the problem onto a small subspace, similar to the way the eigenvector problem is solved (see IV C). To formalize this method, let us introduce the operator I which inserts one copy of A [converged] into L [n] (16), i.e., With that we build up the Krylov subspace K r , and similar with R [n] (17). As in IV C we create an orthonormalized system of basis vectors |L k and calculate I ij the subspace projection of I Keeping in mind that the subspace projection of L [n] is simply given by the vector l j = 1 0 0 0 . . . T we find The number of basis vectors |L i should be chosen such that this approximation is perfect within computer precision. Further errors are introduced by an imperfectly converged A α l αr [converged] s and from the energy overgrow effect described in IV B which should rule out attempts to go for p → ∞. Still, a small amount of the last two errors is acceptable since they have a similar effect as the afore-mentioned extra noise to avoid local minima.

V. APPLICATIONS
In the last section we presented various methods to improve the performance of the iMPS algorithm with long-range interactions. The main subject was to ensure convergence, where special attention was paid to broken translational invariance. This so far troublesome case can now be tackled mainly due to the newly introduced method of superposed multi-optimization (SMO).
The modified iMPS algorithm has superior convergence properties compared to the basic version but it does not surpass its precision, which is determined by MPS and MPO inherited limitations. In cases where both versions converge, the quality of the results is identical. Readers who are interested in the achievable precision of the iMPS method in comparison with analytical solutions are therefore referred to the literature [22,23].
We checked our algorithm with different models. For a reliable basic benchmark, we examined, e.g., states with long-range chiral order in the next nearest neighbor Heisenberg model and found the expected agreement with the results given in the references [31,32].
Here, we will present results for a model of polar bosons described by a Bose-Hubbard like Hamiltonian with a long-range interaction term. In the thermodynamic limit the ground state of this model exhibits symmetry breaking crystalline phases as well as incommensurate phases with algebraically decaying long-range correlations.
The long-range interaction of the Hamiltonian we consider decays as ∼ r −3 . To model this interaction with an MPO the decay is approximated as weighted sum of 20 exponential functions (see also equation (A17) in the appendix and reference [19]).  For odd distances r = 2j + 1 (right picture) n0nr is suppressed by a factor of roughly 100 compared to even distances r = 2j (left picture). To demonstrate the algebraic decay of n0nr the function f (r) = α · r −0.5 + n0nr→∞ with an adequate α is included on top of both plots as a red line.

A. Bose-Hubbard model with long-range interaction
We study the thermodynamic limit ground states of polar bosons in an one-dimensional optical lattice described by the following effective Hamiltonian [11,33] whereĉ † j andĉ j are the creation and annihilation operators for a boson on site j andn j =ĉ † j ·ĉ j . This model is characterized by a hopping amplitude t, an onsite interaction energy U , a chemical potential µ, and a long-range dipole-dipole coupling V /r 3 . For µ > 0, the chemical potential favors as many bosons as possible in the ground state, while the dipole-dipole coupling together with the on-site interaction try to avoid two bosons coming too close to each other. For certain parameter regimes, this interplay allows for translational invariance breaking crystalline phases with optimized distances between the bosons where q sites accommodate exactly p bosons. We will refer to them as p/q-phases. The model is known to host an entire Devil's Staircase of crystalline phases for t = 0 if the joined potential of on-site interaction and dipole-dipole coupling is convex [11,34]. splits into two branches, each decaying ∝ r −α + constχ, with constχ < 0 for the lower branch. This constant is a purely numerical effect, as can be seen by the scaling with the bond dimension χ.
In the following, we will investigate two qualitative different regimes of this model: U → ∞ and U = V.

Devils's Staircase for U → ∞
For U → ∞, each site can accommodate at most one boson and the effective dimension d eff of the local Hilbertspaces reduces to d eff = 2. Figure 3 displays the average ground state densities of the bosons and the localization of the corresponding p/q-phases. To determine the periodicity q of the phases we counted the number of eigenvectors of the transfer matrix T ij with an absolute eigenvalue of one. Once q is known p follows from the average density. Figure 3 (iii) shows a magnification of the area between the 1/4-phase and the 1/5-phase. The biggest phase between these two phases is the 2/9-phase, which can be understood as primary compromise (2/9 = [1 + 1] / [4 + 5]). In the same fashion we find e.g. that the biggest phase between the phase 2/9 and 1/5 is the 3/14-phase. The maximal detectable value of q is given by the bond dimension χ of the MPS, which is 64 in case of Figure 3 (iii). However, since the range of µ covered by the different phases diminishes with growing value of q most phases beyond q = 30 escaped our resolution. The highest value we hit was p/q = 11/52.

Devils's Staircase for
For sufficient small U the ground states of the Hamiltionian (53) might accommodate more than one boson per site, which allows for new types of Devil's Staircases. An example is given by Figure 4, which shows the densities and phases for U = V = 1. Here, simple translational invariance is broken by an underlying occupation pattern given by . . . , x i , 0, x i+2 , 0, x i+4 , 0, . . . with x j = 1 or 2. Of course, for any non-zero hopping amplitude t > 0 we expect fluctuations around this pattern such that a more accurate description might be given by . . . , x i , ε i+1 , x i+2 , ε i+3 , . . . , which we need in the next subsection V A 3. At a certain point, these fluctuations will become so strong that the underlying pattern is destroyed, but this is not the case for the entire region of Figure 4.
In the lobes of the new Devil's Staircase the sublattice . . . , x i , x i+2 , x i+4 , . . . crystallizes in regular pattern of single and double occupied sites. These lobes exhibit an approximate symmetry under the exchange of single and double occupied sites. This is a nontrivial symmetry in contrast to the exact particle hole symmetry of Figure 3 (i).
Outside the crystalline phases Burnell [11] predicted a supersolid like phase. In the following we present numerical evidence which supports this claim.

Supersolids for
A supersolid is characterized as spatially ordered phase which also exhibits superfluid properties. We already mentioned the spatial order, belonging to Fig- ure 4, which is given by the occupation pattern . . . , x i , ε i+1 , x i+2 , ε i+3 , . . . . In our numerical studies we consider translational invariant superpositions of the ground states, where the occupation pattern is still visible in the two-point correlation functions as n 0nr and ĉ † 0ĉ r shown in Figure 5 and Figure 6. Due to the spacial order both correlation functions are split into two branches, where correlations belonging to odd distances r = 2j + 1 are strongly suppressed compared to correlations belonging to even distances r = 2j. Furthermore, both branches exhibit an algebraic decay of the same power. For a superfluide, the power of the decay of n 0nr is supposed to be reciprocal to the power of ĉ † 0ĉ r . Luttinger liquid theory predicts [35] n 0nr = n 0nr→∞ + const 1 · cos (2πρ 0 r) · r −2K + . . .
with ρ 0 = n i . These relations are confirmed by our numerical results, as displayed in Figure 7. Both correlation functions give rise to the same Luttinger liquid parameter K within a small error range.
In the thermodynamical limit algebraically decaying correlation functions go hand in hand with an infinite entanglement entropy of the half chain S, which can not be represented by any MPS with finite bond dimension χ. Nonetheless, it was shown [36,37] that the numerically obtained entropy S χ for such critical phases shows a predictable scaling as function of the MPS bond dimension where c represents the central charge. A demonstration of the scaling behavior is given in Figure 8. From this numerical sample one obtains S χ=256 − S χ=16 = (0.218 ± 0.003) × 4, which matches S χ=256 − S χ=16 ≈ 0.224 × 4 drawn from equation (56) for c = 1.

VI. CONCLUSION
We have presented several extensions to the basic iMPS algorithm for systems with long-range interactions, some of them with the potential to be useful in a much broader context. A special focus was set on problems arising from broken translational invariance. Here, convergence was ensured by various means, but mainly due to the SMO method which irons out local variations by optimizing exponentially many MPS simultaneously. The new algorithm was successfully applied to calculate detailed Devil's Staircases and phase diagrams for polar bosons (53) and was also suitable to verify supersolid properties. Theoretical restraints as considered in the comments III B 2 seem to have negligible influence on practical applications such that the new version of the iMPS algorithm is a genuine improvement in the sense that it can all the old version could plus more. the Austrian Ministry of Science BMWF as part of the UniInfrastrukturprogramm of the Research Platform Scientific Computing at the University of Innsbruck. We further like to thank Andrew J. Daley, Marcello Dalmonte, Johannes Schachenmayer, Ulrich Schollwöck and especially Andreas Läuchli for helpful discussions.

Appendix A: MPO representation for Hamiltonians
For self-consistency we give a short account based on some examples how to construct an MPO representation for a given Hamiltonian (see also [18][19][20][21] First, we need a neat way to write down the explicit form of the fourth order tensors H µ l µr s s . We write them as matrices whose entries are matrices, too H µ l µr s s = H s s µ l µr . (A2) As an example we consider the Ising Hamiltonian state level due to broken translational invariance [R] s1 · q α1α2 [R] s2 · · · (C1) If degenerations due to further symmetries are involved, A and B are supposed to have the same characteristic values for these symmetries. This allows us to operate as if no further symmetries exist, since all operations we are about to use leave these characteristic values unchanged. Under this condition we will prove the existence of a matrix γα 0α0 such that the ground state B can be expressed using the tensors Q stemming from the iMPS A where γα 0 α0 = ωα 0β · ξ βδ · ν δα0 . Let us start by surveying the elements of the proof. In order to show the claimed equation (C2) we will actually prove the gauge transformation [R] s2 · · · (C3) This gauge transformation will be proven for the case that the two underlying MPS represent the same physical state -which is not given for A and B (C1). In order to use the gauge proof for our purpose we need to find one physical state described by two different MPS where the first MPS be constructed using the Q of A and the second by the q of B. This one state which allows us to complete our proof is the translational invariant ground state. According to our preliminary remarks this state is unique. Hence, if we succeed to construct two different MPS which represent a translational invariant ground state, we know that they represent the same physical state as demanded. We will prove the following construction for this state with new tensors τα 0α0 and θα 0α0 . We start by showing the gauge transformation (C3). In order to increase clarity we define Although not an identity, the operator − → Q α k s1···s k · − → Q * α k s1···s k acts like such if it is applied from the left on any MPS with the structure represents an arbitrary right side of the MPS. Using the identity (C6) we get In the following we assume that the two MPS represent the same physical state. Now, let us apply the operator − → Q β k s1···s k · − → Q * β k s1···s k on this equation. Due to equation (C7) the left side stays unchanged. Hence the physical state represented by the right side does not change either with ..sn) is smaller than the bond dimension α k , the value of α k can be reduced, since in that case it turns out to be unnecessarily high. Applying ← − M −1 from the right on equation (C9) we end up with For k → ∞ this covers the first part of the heralded gauge equation (C3). The second part of equation (C3) is proved by a straight forward application of the arguments used above on the right side of the MPS. Now, we have to show the MPS construction of the translational invariant ground state (C4). As above, we assume that the ground state level of the Hamiltonian under consideration is g times degenerate due to broken translational invariance. Let T L be the operator which shifts all sites of an MPS by one position to the left and T R = (T L ) −1 . For the MPS A representing any of the possible ground states we get where T represents an unnormalized version of the translational invariant ground state, which we like to construct. As an intermediate step we like to prove equation (C15) below. Therefor we have to look at the effect T L has on A (using equation (C1) for A ) [R] s2 · · · (C12) Next, we look at the MPS (T R ) g · T L · A, which has the form · · · Q α−2α−1 Since the two MPS (T R ) g · T L · A = T L · A describe the same physical state we are allowed to apply the gauge transformation (C10) and identify · · · Q α−2α−1 [L] s0 ·ω γ α0 .
(C14) Inserting this expression into equation (C12) we arrive at the following description for T L · A T L · · · Q α−2α−1 As we see, applying the operator T L on A has the same effect as the replacement of the matrix λα 0α0 by ωα 0α0 . Since higher powers of T L can be created by an iteration of the arguments just presented, it follows that the effect of any (T L ) j on A can be accounted for by an accordingly calculated ωα 0α0 [j] . Following that construction the only difference between the various MPS (T L ) j · A is their tensor ωα 0 α0 [j] and hence the task of equation (C11) to sum up these MPS reduces to a summation of the ωα 0α0 Let us review our arguments: By virtue of equation (C16) we can transform the MPS A and B given in equation (C1) such that we end up with the translational invariant ground state T Q and T q as claimed in equation (C4). Since we work under the condition that T Q = T q we are allowed to use the gauge transformation (C10) to replace the q of T q by the Q of T Q . The same replacement is possible in B because the q in B and T q are identical (as are the Q in A and T Q ). This concludes the proof of equation (C2) we aimed for.

Appendix D: Davidson implementation
In this subsection we introduce a practical implementation resembling the Davidson [24] method based on recycled information of the previous round, which allows to improve the update equation (37) of the iterative eigenvector solver presented in section IV C. We adopt the same notation as in that section but mostly drop the index [n] to keep the formulae clean.
The best possible new vector|A k+1 the iterative eigenvector solver could come up with to replace equation (37) is an orthonormalized version of |∆ 0 = |E 0 − |e 0 . Hence we use the accumulated data at the end of round n−1 for an efficient one time estimation of D [n−1] , which we will apply in round n.
In order to calculate D we need a simplified form of H which allows easy inversion. We know k approximated eigenvectors |e i<k of H. In order to have an orthonormal basis for H we imagine N − k further |ê k≤i<N , where dim H = N × N . With that we approximate H as H|e with an average eigenvalue α = const. for the unknown eigenvectors. One might be tempted to simplify equation (D3) using H|e i ≈ e i · |e i -but we do not, since the eigenvectors are not very well approximated except for |e 0 . To be able to perform the inversion in equation (D6) it suffices to resort to the exact result e i | H|e j = e i · δ ij without i , which is a consequence of the construction (35). Further, with the results gathered during the optimization (IV C) the H|e i are as quickly calculated as the |e i . Next, we take the trace of equation (D3) and set α such that both sides are equal tr H = e i | H|e i + α · (N − k) The trace of the exact H is efficiently calculated by already tracing over its components L α l µ l α l and R α r µrαr before assembling them (15). Now we insert the approximated H (D3) in D The inversion is solved by as can be verified inserting the result in D · D −1 = I. The final question we have to answer is which value we assign to the unknown exact eigenvalue E 0 . The best approximation (which we already used in equation (D1)) is E 0 ≈ e 0 , but this produces a singularity in D. There are two ways out. First, we can always pick E 0 a little bit lower E 0 := e 0 − ε. Second, we should discard the troublesome term ∼ |e 0 e 0 | in D anyway for the following reason: We replace equation (37) by and with that |e 0 e 0 |r ≈ 0·|e 0 . Afterwards, residual parts of the |e 0 e 0 | term are exfiltrated again because |A k+1 has to be orthogonalized (and normalized) These considerations are also part of the more elaborated Jacobi-Davidson [24] method to which this implementation can be extended. We might further consider to omit terms with e i e 0 since their influence shrinks with (e 0 − e i ) −1 . At the end of the day, the effort to construct D as well as the effort for each application scale with N times the number of |e i e i | terms used in D.