Exploring phase transitions by finite-entanglement scaling of MPS in the 1D ANNNI model

We use the finite-entanglement scaling of infinite matrix product states (iMPS) to explore supposedly infinite order transitions. This universal method may have lower computational costs than finite-size scaling. To this end, we study possible MPS-based algorithms to find the ground states of the transverse axial next-nearest-neighbor Ising (ANNNI) model in a spin chain with first and second neighbor interactions and frustration. The ground state has four distinct phases with transitions of second order and one of supposedly infinite order, the Kosterlitz-Thouless transition. To explore phase transitions in the model, we study general quantities such as the correlation length, entanglement entropy and the second derivative of the energy with respect to the external field, and test the finite-entanglement scaling. We propose a scaling ansatz for the correlation length of a non-critical system in order to explore infinite order transitions. This method provides considerably less computational costs compared to the finite-size scaling method in [8], and quantities obtained by applying fixed boundary conditions (such as domain wall energy in [8]) are omitted. The results show good agreement with previous studies of finite-size scaling using DMRG.


I. INTRODUCTION
Recently, numerical simulations of many-body quantum systems based on matrix product states (MPS) and their generalizations have called increasing attention as they offer an efficient way to study properties of their ground states and even thermal states. One of the main reasons of their popularity is their applicability to-in principle-any Hamiltonians describing only local interactions on a spin chain and their possible generalizations to higher dimensions. The use of the so-called projected entangled pair states (PEPS) to approximate ground states of such kind of Hamiltonians defined in higher dimensions (e.g. on a square lattice) is considered very promising. (See e.g. Ref. [1] and references therein for a review of the topic.) These methods can provide efficient simulation of systems when other methods such as quantum Monte Carlo break down, e.g. in case of certain frustrated systems. It means for example that one can circumvent the so-called sign problem present in the Monte Carlo simulation. (Note that there exist systems for which there is no other method available, not considering exact diagonalization which works only for small system sizes. One such system is e.g. the frustrated XX model. Ref. [12]) Obviously one of the main goals of investigating ground states is to explore phase transitions and characterize them as precisely as possible. Despite the fact that by definition, the MPS can only describe exponentially decaying correlations, one can obtain accurate results for certain quantities of phase transitions like in the Ising model as it has been demonstrated in nearly all of the papers presenting an MPS based algorithm. (See e.g. Ref. [2] where it has also been used for an Ising model defined on a binary tree.) Models with not only nearest neighbor (NN) interactions and higher than second order transitions have not been studied extensively by means of MPS. This paper has two goals: Find a suitable efficient algorithm based on MPS to simulate Hamiltonians with second neighbor interactions and frustration, as well as demonstrate the use of universal tools, especially finite-entanglement scaling to study different types of phase transitions. Its validity arises from general facts about critical points and the definition of the MPS, hence it can be applied quite universally. As opposed to the usual finite-size scaling it may also have a better scaling for the computational costs if using the same maximal bond dimensions for both cases (see below). Here the 1-dimensional version of the transverse axial next nearest neighbor Ising (ANNNI) model has been chosen which is non-integrable and possesses a quite complex phase diagram, despite being the simplest model with NNN interactions. This model is notable due to several facts. As earlier studies have shown, unlike the Ising model, it has a critical region, i.e. a 2D parameter region where it is critical. A probably infinite order phase transition also takes place (besides second order phase transitions). As mentioned earlier, the detection of an infinite (e.g. Kosterlitz-Thouless) transition is a notoriously hard computational challenge, that's why a computationally beneficial method is needed. The same holds for exploring the physics of a critical region.
In order to study such a model one can choose a finite chain of length L with periodic (PBC) or open boundary conditions (OBC) and then look at the scaling of some physical quantities depending on L and the bond dimension D. It is known that the overall computational costs scales O(D 3 ) in the best case. The latter happens for example when the correlations are very small ξ ≪ L as it is discussed in Ref. [9], however when ξ ≈ L e.g. it can be as large as O(D 5 ) and the precise behaviour of the scaling as a function of D and L is unclear.
Note that in this paper, iMPS and imaginary time evolution are used, hence one works directly on the infinite lattice and the only scaling parameter is D. In this case the computational cost in every small time step ∆t scales as O(D 3 ). Roughly speaking, convergence can be reached if: exp(t∆E) = exp(n∆t∆E) ≪ 1 with the ∆E energy gap and n the number of time steps. Approaching to a critical point, ∆E vanishes as a function of the effective correlation length ξ D , and at the same time we also have a scaling relation like ξ D ∝ D k (see below). With this in mind, it can be estimated how the overall computational time scales as a function of D and one sees why it is slower than O(D 3 ) near to a critical point. For instance, if ∆E ∝ ξ −α D ∝ D −αk , α > 0 then n should scale as O(D αk ) that makes the overall scaling O(D 3+αk ). Previous work about finite-D scaling can be found e.g. in [7] and [10]. Now the idea is to study the scaling of physical quantities not as a function of the system size, but simply as a function of the bond dimension D. Instead of finite-size scaling this can be called as 'finite-D' or 'finite-entanglement' scaling as a MPS with a given D can contain a finite amount of entanglement. As it will be argued in Section IV. the overall computation cost of finite-D scaling can be considerably reduced compared to that of the usual finite-size scaling. In fact, this turns out to be the case for the presented model as well as for the Ising and Heisenberg models, and essentially this depends on the central charge. The key point to this computational cost reduction is our off critical scaling ansatz (4) for the correlation length as a function of D which turns out to be computationally more beneficial than that of used in Ref. [8] for the entanglement entropy.
Clearly, our methods can be applied straightforwardly to any other 1D model with NN and NNN interactions, such as the J 1 − J 2 Heisenberg model, but can also be adapted for the 2D version using PEPS. Despite the polynomial computational time scaling, unfortunately its exponent is so high in the 2D case that in practice one can apply only very small bond dimensions providing only a few and ambiguous results that are not yet sufficient for numerical analysis.

II. THE 1D TRANSVERSE ANNNI MODEL
Consider the following Hamiltonian on an infinite chain of spin-1/2 particles This can be considered as the 1 dimensional version of the lattice model with axial next nearest neighbor interactions (ANNNI). If the next nearest neighbor (NNN) interaction term describes antiferromagnetic interactions (J 2 < 0) and J 1 > 0 then the system is frustrated. Note that basically two parameters determines the physical properties of model e.g. κ = J 2 /J 1 and h/J 1 . Changing the basis of every second spin appropriately by rotating around the axis x leads to Z 2i → −Z 2i , so one can switch from an antiferromagnetic NN interaction to a ferromagnetic one and vice-versa, leaving the NNN interaction unchanged. However, there is no rotation switching all couplings to ferromagnetic, thus the sign of the NNN interactions is crucial and besides the magnetic field of h/|J 1 | the ratio κ = J 2 /|J 1 | plays the key role in determining the physical properties. From now on let us set J 1 = 1 > 0 and J 2 < 0 yielding competing ferro-and antiferromagnetic NN and NNN couplings.

III. MPS STATES AND ALGORITHMS
The MPS representation of an infinite, translational invariant chain of d-dimensional systems looks like The normalization constraint requires that the largest eigenvalue of the transfer matrix E = s A s ⊗ A s has to be unique and equals to 1. If it's degenerate (as it happens e.g. for the state Ψ = c 0 ↑↑↑ ... + c 1 ↓↓↓ ...) then it can cause troubles in some types of the algorithms that work with unique leading eigenvectors (Ref. [5], Ref. [6]). Nevertheless, in order to overcome this one can apply tricks like appropriate rotations of some basis vectors, but due to the lots of computations like leading eigenvectors, inverses and squares of matrices these algorithms turn out to be costly. Another representation Ref. [4] uses the Schmidt coefficients in the diagonal matrices λ i explicitly as follows with the normalization constraints in the general case when the dimensions of λ (i−1) and λ (i) can be different. The infinite time evolving block decimation (iTEBD) Ref. [4] uses this description and in case of NN interactions two neighboring sites are updated after the action of e −δtHi,i+1 . One builds up a larger matrix occupying the MPS matrices of the two sites i and i + 1 and decomposing them through a singular value decomposition to obtain the new matrices. A normalization step can be interpreted as a zero time update: we just build a larger matrix and decompose it. The details can be found in Ref. [2]. Several ways exist to do the updates in the case at hand. One option is the use of 'superspins' ( made of 2 neighboring spins). To this end let us join together every second neighboring pairs of spins to form block spins with 2 × 2 = 4 dimensional basis as |t i , t i+1 |t i+2 , t i+3 = |s i |s i+2 with t i ∈ {0, 1} and s i ∈ {0, 1, 2, 3}. In this way we get a translational invariant Hamiltonian with only NN interactions acting between these 4 dimensional block spins so we can use the simplest algorithms designed for this case but with d = 4 and a bit more complicated NN term.
In our scheme, for example the operator σ z i ⊗ σ z i+1 ∈ M 4 (C) acting on the sites i, i + 1 in the original lattice is interpreted as an operator acting on the single superspin labelled by i, describes an interaction between the superspins at the position i and i + 2.
is a NN interaction between the superspins, however it describes a NNN coupling between the original spins. With this in mind one can easily construct the operators in the new basis. We make use of the 'superspins' and the update method described in Ref. [2]. The computation time in every step scales as O(D 3 ) but as opposed to the previous method like in Ref. [6], one does not need to compute neither leading eigenvectors nor inverses or squares of non diagonal matrices hence it is faster and there is no problem with degenerate eigenvalues.

IV. FINITE-ENTANGLEMENT SCALING TO DETECT PHASE TRANSITIONS
In order to characterize phase transitions using iMPS one needs to study how some special quantities scale as a function of the applied bond dimension D. Let us denote the entropy of entanglement computed from the iMPS by S D and the correlation length ξ D . It is known that in the critical points S D converges logarithmically with corrections of order O(1/ log D) as with c being the central charge and k only depends on c for large enough D as k = 6 The scaling of the correlation length looks like again with corrections of O(1/ log D), see Ref. [7], [9], [10] and [11] for more details.
In practice however, one should go in the other way round: Testing the scaling relations one wants to decide whether the system is critical or not by giving an estimate for the correlation length. The closer we are to the critical point the more delicate the problem due to the high correlation length. One has to decide if the limits S = S ∞ = lim D→∞ S D and ξ = ξ ∞ = lim D→∞ ξ D are finite or not. (S and ξ converge to the exact values of the infinite system as the Trotter errors tend to zero.) It can be especially hard for example in the vicinity of a Kosterlitz-Thouless transition when the inverse of the true correlation length ξ −1 tends to zero very quickly as we are approaching the critical point. As in finite-size scaling let us assume some scaling function for ξ D as Clearly, g(D, ξ ∞ ) ≈ log(ξ ∞ /D k ) for large D if ξ D → ξ ∞ < ∞. It enables us to define f (x) = g(D, ξ ∞ ) with x = D k /ξ ∞ , for which f (x) → 0 if x → 0 and f (x) ∼ − log(x) for large x. Note that besides these limit properties we will not make any further observations about the form of f , just use a simple ansatz what will be justified by the numerics to be a good indicator for the non-critical behaviour.( See below the similar ansatz function for the finite-size case). So for instance, the following simple ansatz meets these requirements Near the critical regime ξ is big but not infinite and for a not big enough D it behaves very similarly to (3) and one cannot distinguish it from the critical scaling. So the natural question arises: For a given and large but finite ξ near the critical regime, at least how big D should be applied in order to rule out the critical scaling with some reasonable accuracy? Essentially it depends on the ratio of log D and f ( D k ξ∞ ). In order to decide whether a given point is critical or not one should test the linearity of log ξ D against log D as a null hypothesis. If it can be rejected then one can fit a trial function according to (4) and (5) to obtain an estimate for ξ ∞ .
Note that(4) resembles of the finite-size scaling of the entropy of the half chain of length L applied in Ref. [8] with the same ansatz for f (x) but now one has to plug x = L/ξ ∞ . Note that this simple ansatz (5) for f (x) has the right limits but for intermediate values it is not correct, however it turns out to be good approximation for S(L) and a good indicator for non-criticality by the numerics. One can obtain (universal) series expansions for f (x) in the small and large limit theoretically, see Ref. [13]. One can also confront this with our ansatz which incorporates the to limits. Of course here, for a given L, the applied bond dimension (or the number of states kept in each DMRG iteration) should be "high enough".
In our case however, we only have the bond dimension as an input parameter which simplifies the picture. Besides, the main advantage is as follows: Having the same ξ ∞ for the infinite system, for the same accuracy one has considerably smaller computational cost if the scaling exponent k > 1. . In order to detect off critical scaling with the same accuracy one expects the necessary maximal length L in the finite-size case must be bigger than the maximal D in the finite-D case, if k > 1. Nevertheless, this can be estimated by studying the scaling anzatz (4) and (6). In terms of D this would result in an overall effective scaling higher than O(D 5 ), in the finite-size case, obtained by the assumption L ∝ D. [1] shows the phase diagram as it has been found in earlier studies by simulations Ref. [8] and perturbative analysis Ref. [3]. We try to explore the phase boundaries by looking at the correlation lengths as a function of the parameters and the second derivative of the energy. A Kosterlitz-Thouless type transition is conjectured between the floating and the paramagnetic(PM) phase where the derivatives of the energy are useless indicators so the expected divergence of the correlation length is checked only. In the critical points we also check the finite entanglement scaling laws and find very good agreement with the theoretical results. Moreover, in some special points we calculate the critical charges as well as the correlation scaling exponents.

Figure
We intend to determine some points of the line of phase transitions between the distinct phases and check these relations as well as the finite-entanglement scaling and compute the central charge. There are three lines of phase transitions, two of them are of second order while an infinite order transition is expected for the third one which needs a special care. We also try to detect the floating phase when κ < 0 and |κ| is large, since its existence for arbitrary large |κ| is still an open question.

A. FM-PM transition
In case of the Ising model, in the vicinity of the critical point as a function of the magnetic field, the correlation length and the magnetization diverges as M ∝ |h − h c | β and the energy and its first derivative are continuous but its second derivative diverges at h = h c . Using D = 40 we get 0.535 < h c < 0.538 by looking at the extrema of the finite derivative ∂ 2 E/∂h 2 . At the transition corresponding to κ = −0.25 we find the same behaviour for the energy as a function of h and a very good agreement with ν = 1 as in the simple Ising case: ν = 1.02 comes from the linear fit on the figure below.

C.
Floating phase-paramagnetic phase transition First, we check the power-law scaling of ξ assumed for critical points as in (3). If this power-law scaling can be rejected statistically, then for a given h and κ, ξ ∞ is computed by the best fit using (4) and finite-D scaling. (See the inset in FIG.4.) We observe that (4) fits well using the ansatz function (5) with weak dependence of the parameter α which can be set to 1.
In fact, the numerics indicate a higher order transition here, with continuous derivatives of the energy. Probably a Thouless-Kosterlitz type transition takes place as suggested by the plot ξ −1 vs h  FIG. 5). Otherwise, where we cannot reject the power-law scaling hypothesis, even though a finite value for ξ ∞ obtained by fitting (4) is rejected. In these cases one cannot rule out the critical scaling statistically. We have got anyway quite ambiguous results for ξ ∞ which -in this case-can be regarded as a consequence of the errors of the simulation. As in Ref. our division for the magnetic field, the critical value can be estimated roughly as 0.42 < h c < 0.44. Note that here, the study of model specific quantities have been avoided.

D. Large frustration
At κ = −10 up to the the applied maximal bond dimension D max = 80 and the division for the magnetic field, the numerics do not seem to indicate the presence of the floating phase because one sees a sharp peak plotting ξ against h. See Fig.6.
Note that in Ref. [8] the existence of the floating phase up to κ = −5 has been indicated but the DMRG method has become imprecise for higher frustrations, however with our new method, the floating phase can be suspected even with D max = 50 as shown in FIG. 7.

VI. SUMMARY AND CONCLUSIONS
The phase diagram of the 1D transverse ANNNI model is studied. We have coupled together neighboring sites of the chain in order to get NN interactions between the 4-dimensional superspins and applied the simple iTEBD with the new NN Hamiltonian. In practice, the iTEBD algorithm seems to be more beneficial than others based on the transfer matrix and its leading eigenvectors. At some typical points of the phase diagram we have studied the phase transitions and estimated the critical quantities and found good agreement with previous DMRG results coming from finite-size analysis and model specific quantities. We have also checked finite-entanglement scaling and its breaking to localize a supposedly infinite order transition between the floating and the paramagnetic phase, as well as pointed out the logarithmic divergence of the entanglement entropy inside the floating phase.  field ∆h = 0.08. However, we can suspect the existence of the floating phase for κ = −5 already with D <= 50. The question of its existence for large frustrations requires a more detailed study with higher bond dimensions, but the applied methods throughout the paper are quite universal and applicable to other frustrated spin systems as well.

VII. ACKNOWLEDGMENT
The author is grateful for the stimulating discussions with Lorenzo Campos Venuti and Marco Roncaglia, as well as their hospitality at the I.S.I. Foundation in Turin, Italy.