Bilinear-biquadratic spin-1 rings: an SU(2)-symmetric MPS algorithm for periodic boundary conditions

An efficient algorithm for SU(2) symmetric matrix product states (MPS) with periodic boundary conditions (PBC) is proposed and implemented. It is applied to a study of the spectrum and correlation properties of the spin-1 bilinear-biquadratic Heisenberg model. We characterize the various phases of this model by the lowest states of the spectrum with angular momentum J = 0, 1, 2 for systems of up to 100 spins. Furthermore, we provide precision results for the dimerization correlator as well as the string correlator.


I. INTRODUCTION
About 25 years ago the density-matrix renormalization group (DMRG) emerged as a precision tool for the numerical description of one-dimensional quantum many-body systems [1]. Starting with Ref. [2] it was applied successfully to many spin and strongly interacting electron systems. Later it was realized that DMRG can be reformulated in terms of matrix product states (MPS) [3,4]. A comprehensive review of MPS algorithms and their relation to DMRG is presented in Ref. [5]. In this review most algorithms are formulated without regard to symmetries such as U(1) or SU (2).
However, alreadyÖstlund and Rommer [3,4] in an attempt to understand the success of DMRG used SU (2) symmetric MPS. They realized that the local matrices (tensors) decompose into a structural part given by the symmetry and a degeneracy part. The structural part consists of a Clebsch-Gordan coefficient, and the (variational) parameters of the model reside in the degeneracy part only. The number of parameters to be determined is therefore significantly reduced with respect to a non-symmetric theory. This approach, first suggested byÖstlund and Rommer, was generalized to higher order tensors in Ref. [6] and then applied to MERA tensor network calculations [7][8][9].
However, in addition to this basic implementation of symmetry into a tensor network with the purpose of reducing the number of independent parameters, it is possible to entirely eliminate the structural tensors and develop algorithms in terms of the degeneracy tensors only [10][11][12][13]. This reduces the requirements for computational resources significantly, and in turn enables significant improvements to the accuracy of the results that can be obtained. The elimination of structural tensors from the algorithms requires changes to standard (nonsymmetric) tensor network implementations, e.g. via precomputation schemes as suggested in Ref. [9].
It is the purpose of the present paper to develop an algorithm for SU (2) symmetric MPS with periodic boundary conditions (PBC) in terms of reduced tensors only. We stress that we express symmetric tensors (using the Wigner-Eckart theorem) in terms of reduced tensors (degeneracy tensors) and structural tensors consisting of products of Clebsch-Gordan coefficients. To this end we strictly follow the methods and conventions of Edmonds [14]. We do not use the tree decompositions advocated in Ref. [9] for the representation of symmetric tensors. The linear maps relating different tree decompositions of symmetric tensors derived in Ref. [9] directly correspond to expressions relating different coupling schemes [14] in terms of Racah 6j, Wigner 9j, or more general symbols. Consequently, the precomputations for the linear maps suggested in Ref. [9] may be expressed in terms of such symbols.
Furthermore, we apply the proposed algorithm to a physically rather complex one-dimensional model: the spin-1 bilinear-biquadratic Heisenberg (BBH) model on a ring, cos θ s i ⊗ s i+1 + sin θ ( s i ⊗ s i+1 ) 2 (1) with N + 1 set to 1 and s i the spin-s SU(2) matrix representations. In doing so we will reproduce a number of well-known results, e.g. a precise calculation of the Haldane gap in order to check the capabilities of the proposed algorithm. Moreover, we will study the spectrum of the BBH Hamiltonian for PBC and various θ, which have not yet been addressed in the literature. This will shed some light on old questions concerning the phase structure to be discussed below. The BBH model describes the behavior of atomic spinor condensates in optical lattices [15][16][17][18][19] under certain conditions. It also models the physics of some quasi-one-dimensional crystals, e.g., in LiVGe 2 O 6 [20] or Ni(C 2 H 8 N 2 ) 2 NO 2 ClO 4 (NENP) [21]. Furthermore, it was extensively used as a test bed for new tensor network algorithms [22][23][24]. However, these algorithms do not explicitly implement SU(2) symmetry, and the methods developed here (and for infinite systems in Ref. [25]) could serve well to implement SU(2) symmetry into those approaches.
The continuous SU (2) symmetry of the BBH model cannot be broken due to the Mermin-Wagner-Coleman theorem [26,27]. As a consequence the eigenstates of this Hamiltonian can be characterized by the total angular momentum quantum number J. Due to the (discrete) translational symmetry, the eigenstates can also be labelled by the quasi-momentum quantum number p. However, unlike SU(2) symmetry, translational symmetry will not be built into our MPS ansatz explicitly in the present paper. The BBH model has other symmetries not explicitly built into the MPS here, e.g. SU(3) symmetry at θ = 1 4 π and θ = − 3 4 π, a symmetry not easily uncovered in the spin representation of the model [28].
As a function of the control parameter θ the infinitesize spin-1 bilinear-biquadratic Heisenberg model exhibits the rich phase structure shown in Fig. 1 [29]. We briefly discuss the various phases moving around the circle of the phase diagram in a clockwise direction.
For π 4 > θ > − π 4 the system is in the gapped Haldane phase with hidden topological order [30,31]. The ground state has spin J = 0 and is non-degenerate. The first excited state is a triplet (spin-1). At θ = 0 the model corresponds to the simple Heisenberg antiferromagnet and at θ = arctan 1 3 the ground state is the AKLT state [32]. At θ = − π 4 the gap vanishes and there is a second order phase transition to the dimerized phase [33,34].
According to [39,40], the dimer order parameter behaves non-typically in the parameter region θ c ≈ −0.67π > θ > − 3π 4 , and quadrupolar spin correlations increase dramatically there. The conjecture of Chubukov [41] that there is a gapped nematic phase in this parameter region was debated in Refs. [16,42] and rejected in Ref. [43]. However, the existence of the nematic gapless phase in the narrow vicinity of θ = −3π/4 is not ruled out, since the values to be calculated are comparable with the precision of the numerical calculations in that parameter region [40,44].
Using the proposed SU(2) symmetric MPS algorithm for PBC, we calculate energy spectra and characteristic correlation functions in different phases of a spin ring at selected parameter values θ. The calculated energy gaps between these states enable already an elucidation of the phase structure from finite system results. In particular, the existence of a fifth (nematic) phase can be addressed. For the Haldane phase we calculate the string correlator of the ground state and in the dimerized phase the dimerization from the lowest two J = 0 states. These two states form a degenerate doublet in the thermodynamic limit.
The paper is organized as follows. In Sec. II we briefly review the MPS formalism for periodic boundary conditions which is based on the algorithm proposed by Verstraete, Porras, and Cirac [48]. In Sec. III this algorithm is rewritten using SU(2) symmetric tensors only. As already emphasized, it is a major objective of the present paper to eliminate all structural tensors from the PBC algorithm. This is more complicated than for OBC since higher order tensors must be considered. These technical developments are relegated to several appendices, which form an important part of the present article and should enable a straightforward implementation of the algorithm.
As a consequence we only need to handle the degeneracy parts of the tensors explicitly. This enables significant improvements to the efficiency of the implementation. In fact, the MPS virtual dimensions we are able to use are much larger than in various recently proposed implementations for PBC without SU(2) symmetry, e.g. [22,23]. We provide the reduced MPO representation for the bilinear-biquadratic spin-s Hamiltonian as well as the corresponding reduced representation for H 2 and other operators. This enables the calculation of the variance H 2 − H 2 of the various eigen-energies.
We apply the proposed SU(2) symmetric algorithm to the bilinear-biquadratic spin-1 Heisenberg model in Sec. IV and discuss the low lying spectrum in all phases except the ferromagnetic phase. We also briefly address in Appendix F the spin-1/2 Heisenberg model and compare our numerical results with Bethe Ansatz calculations. The results of our work are summarized in section V.

II. REVIEW: THE MPS FORMALISM FOR PBC
An MPS formalism for PBC was originally proposed in [48] and extended in [49] and [18]. We summarized the algorithm in Refs. [50,51], and therefore we only briefly review here those aspects which are relevant for the present discussion.
The state of a 1D quantum many body system of size N is approximated in terms of a matrix product state Here, the σ i represent the local degrees of freedom at site [i]. The local Hilbert space is assumed to be finite dimensional, and its basis is enumerated by the {σ i }. The elements M Analogously, operators are written as matrix product operators Each W [i] is a rank-4 tensor: W can be expressed in terms of the rank-6 (generalized) transfer tensor E W with the tensor elements The matrices A and B characterize the states |φ and |ψ , respectively. The special transfer tensor E 1 (A, B) represents the matrix element φ|ψ of the identity operator (withĪ = I = 1). Note, that we chose to write the algorithm in terms of higher-rank tensors in order to ease the implementation of SU(2) symmetry in the next section. Parentheses are used to group indices conveniently.
In order no minimize the energy variationally, due to PBC one has to solve at each update step the generalized eigenvalue problem written in terms of the effective Hamiltonian matrix H eff and the effective normalization matrix N eff . The elements of these matrices are obtained from R are contractions of the transfer tensors E W and E 1 (as defined in Eq. (5)) to the left or to the right of the site i, respectively. (In order to define what is left and right we initially label the sites from 1 to N , and keep this labeling throughout the calculation. All sites with label smaller than i are left of site i and all sites with label larger than i are right of site i.) A contraction of a block with a transfer tensor is depicted in Fig. 2. The block and transfer tensors are kept rank-6 throughout the algorithm. We assume summation over repeated indices, however, for clarity we sometimes write out the summations explicitly. Brackets are used in order to indicate fusion of several indices into a combined single index. The update is performed in the standard DMRG 'sweeping' manner as proposed in Ref. [48].
The value ǫ [i] eventually converges to the ground state energy during the iterative update procedure. The updated MPS is obtained by a suitable partitioning of the vector into a tensor: After each update step the local MPS tensor is regauged in order to keep the algorithm stable, i.e. we modify each local MPS tensor M [i] such that one of the following relations is fulfilled, The algorithm presented above not only allows the variational determination of ground states but also the calculation of excited states. They can be constructed by finding the lowest state in the (sub)space orthogonal to the space spanned by the states already found [18]. In order to implement a corresponding update procedure for an MPS in this subspace we need to solve the generalized eigenvalue problem where P [i] projects the effective Hamiltonian (7) and the effective normalization matrix (8) into the relevant subspace. I.e., the local MPS tensor must not only minimize the effective Hamiltonian, but must also be constructed in such a way that the MPS |ψ of a 'new' excited state to be calculated is orthogonal to all MPS |φ k already calculated before: ψ|φ k = 0 for all k. Here, k enumerates these states: k = 0 for the ground state, k = 1 for the first excited state, etc.
We denote the tensors of the states |φ k by Φ  We now introduce SU(2) symmetry into the formalism presented above: we assume an SU(2) symmetric Hamiltonian and SU(2) symmetric states, i.e. the symmetry is not broken. This holds in 1D due to the Mermin-Wagner-Coleman theorem. Technically, all tensor indices can then be chosen to be SU(2) invariants, i.e. they decompose into a total spin (angular momentum) index s, a spin projection index m s as well as a degeneracy index t, which counts the number of times a specific SU(2) representation occurs, e.g. σ = (s, t, m s ). For the physical indices of MPS we have t = 1, however the spins in the bond indices are highly degenerate. In practice, we only need moderate degeneracies t < 10 for our calculations.
First we briefly introduce SU(2) symmetric MPS and discuss the construction of states with given total angular momentum J. We eliminate all structural tensors from the algorithm in subsection III B as well as several appendices.

A. Construction of SU(2)-symmetric MPSs for PBC
First we discuss how to construct SU(2) invariant states, i.e. states with total spin J = 0. The set of MPS matrices defined in Eq. (2) must be brought into an SU(2) symmetric form at every site i.
The structure of the rank-3 local tensor M [i] at position i is determined by the Wigner-Eckart theorem [9] according to which the tensor decomposes into a degeneracy part M and a structural part C with the structural part fixed by SU(2) symmetry and given in terms of the Wigner 3j symbol [52] The degeneracy part (also often called "reduced tensor") is denoted by a script letter while the corresponding full tensor is denoted by a roman letter. Here, for simplicity, we omit the position index i (this index will be reintroduced whenever necessary). It is often convenient to use The degeneracy part M does not depend on the spin projections and contains the (variational) parameters of the state. Note that the reduced matrix elements of the MPS for which the spins j, s, j ′ do not fulfill the 'triangle rule' can be set to zero, as the corresponding 3j symbol vanishes under these conditions. This makes the reduced MPS matrices rather sparse, in fact, they have a banded block structure.
We will call the set of all (jt) which characterize a reduced matrix its degeneracy set d, e.g. d = {(1/2, 2), (3/2, 4), (5/2, 3)}. In practice, we want to choose the smallest possible set with few different spins j and small degeneracies t. However, it must be noted that at present we have no method in order to choose this set algorithmically. Of course, for a large system we expect that the required degeneracies will be large. On the other hand the size of the degeneracy sets is strictly limited by the available computing resources. We will discuss this issue further in the following sections.
We now turn to the construction of covariant states with total spin J = 0. To achieve this a fictitious noninteracting (local) spin J > 0 [8,11] is inserted into the system. The fictitious spin is inserted at the site N + 1, i.e. between site N and site 1. The tensor at this additional site takes the form and the resulting MPS has total spin 0. For completeness, we define F 0,0 = 1. After optimization of the tensor network with the fictitious spin one needs to obtain the covariant state |JM , which is given by This state must be normalized if required. In practice, we seldom need to reconstruct the state using Eq. (17), since the calculation of SU(2) invariant observables can be expressed in terms of the reduced tensors only as is shown further down in this section.

B. The optimization algorithm for SU(2) symmetric MPS for PBC
In order to make best use of the symmetry, it is advantageous to introduce a reduced tensor structure not only for the MPS but for all other tensors of the algorithm as well. This enables to express the whole algorithm in terms of degeneracy tensors only, and spin projection indices and the structural parts of the MPS will be completely eliminated from the algorithm. As a consequence, the computations are significantly more efficient, and calculations with much larger virtual dimensions become feasible.
In the present section we will explain how the PBC algorithm presented in section II is rewritten in terms of reduced tensors only. Technical details will be relegated to several Appendices. The results contained in these Appendices will enable a rather straightforward implementation of the proposed algorithm. For open boundary conditions a similar strategy was followed by McCulloch [13] but not described in much detail. However, for PBC we face a number of differences, in particular, the 'block' tensors, i.e., the products of transfer tensors, are rank-6 tensors for PBC, while they are rank-3 in the OBC implementation.
In the following, we will provide the general definition of reduced tensors and exemplify the construction of reduced tensors for rank-4 W tensors, which are the building blocks of MPO. We then go on to describe the variational algorithm for the determination of MPS in terms of reduced tensors.
In order to introduce rank-k SU(2) invariant tensors T a1,...,a k each index has to be decomposed into a spin index, a degeneracy index, and a spin projection index, e.g. a 1 = (j 1 , t 1 , m j1 ) in the same way we did for the tensor indices of MPS tensors in Eq. (14). For k ≥ 3 an element of a SU(2)-invariant rank-k tensor may be obtained from the generalized Wigner-Eckart theorem [6,  . Each vertex represents a C factor or 3j symbol. Summation over the indices of internal edges is implied.
A rank-k intertwiner may be decomposed into a product of C factors defined in Eq. (15). Different decompositions are possible and they may be represented as different coupling schemes as explained in more detail in Ref. [14] and Ref. [9]. Here, we will choose coupling schemes which we find convenient for our purposes.
Let us illustrate the construction of SU(2) invariant tensors using an MPO tensor W as example. Such tensors have rank-4, therefore we need one intermediate index e = (j e , m e ). The MPO has two physical indices (s, m s ) and (s, m ′ s ) (assuming that at each site there are particles with the same spin s) and two virtual indices with a rank-3 reduced tensor W je γ,γ ′ (the factor √ 2s + 1 is introduced for convenience in order to free the reduced unity tensor from an s dependence). Our definition of the reduced tensor corresponds to the coupling scheme shown in Fig. 3. We note that t e = 1 because fusion of two equal spins s gives a non-degenerate spin. Due to the restrictions imposed by the structural C factors as well as the high sparseness of the MPO in its full form, many reduced tensor elements can be chosen to be zero, and the reduced tensor W assumes a characteristic sparse structure. This is demonstrated by the reduced tensor representations of the Hamiltonians needed in the present paper given in Appendix A. We also provide a reduced representation for the H 2 operator in Appendix D.
As was already stated, the W tensors of many important operators can be obtained exactly as well as their reduced tensor representation W. This representation is characterized by a degeneracy set d = {(j 0 , t 0 ), . . . (j n , t n )}, e.g. for the BBH Hamiltonian one Similarly, one can expand higher rank tensors like the block tensors H R and H L into reduced tensors. Details are given in Appendix B. As a consequence, the optimization (update) step Eq. (6) can be formulated as a generalized eigenvalue problem in terms of the reduced effective Hamiltonian H eff and the reduced effective normalization matrix N eff (both given explicitly in Eq. (C1)) We stress again that the reduced normalization matrix arises due to PBC and both H eff and N eff change at each iteration step. The updated matrix M is then obtained from a suitable partition of the vector ν into a matrix, as well as zeros for those tensor elements which do not fulfill the triangle rule. Details of the derivation of the optimization step in terms of reduced tensors are given in Appendix C. A reduced form of the projection operator P defined in Eq. (12) is also given there. This enables the calculations of excited states in different spin sectors. The regauging step (Eq. (10)) in terms of reduced tensors is briefly described in Appendix E. A first demonstration that the proposed algorithm works as expected is given in Appendix F.

C. Calculation of observables
If the operator O of an observable is SU(2) invariant, its expectation value Eq. (4) can be calculated using reduced tensors only. In the present paper we do not consider covariant operators. Using Eq. (3) for the operator and Eq. (17) for the wavefunction |JM , one can write The expectation values of the SU(2) invariant operator are equal for different spin projections: JM |O|JM = JM ′ |O|JM ′ . Therefore, we can write if the MPS tensor F at the fictitious site is normalized such that Q R = 1. For simplicity we did not label the MPS matrices with their corresponding local spin index.
In the last line we have represented the expectation value as a reduced left block tensor as defined in Appendix B. It is this represenation which we use for a recursive calculation of expectation values using the formulas provided in Appendix B. An example of an SU(2) invariant operator is the dimerization operator (22). On the other hand, the string correlator (23) is not SU(2) symmetric due to the parity factor i+l−1 j=i+1 exp(iπs z j ). For such operators a full MPS needs to be reconstructed using (17), and the MPO in its full form must be used.

IV. BILINEAR-BIQUADRATIC SPIN-1 HEISENBERG MODEL WITH PBCS
In this section we present a numerical study of the spin-1 bilinear-biquadratic Heisenberg (BBH) model Eq. (1). This model has been studied extensively in recent years, and we have given a number of pointers to the literature in the Introduction. The present study is distinguished by its particular treatment of SU(2) symmetry in periodic systems. With moderate computational effort and relatively small numbers of variational parameters we obtain results of high precision. We present results for the low lying spectrum of the BBH model as well as for selected correlation functions. In order to do calculations using the proposed algorithm one needs to choose an appropriate degeneracy set for the virtual spins. The larger this set the larger is the computational effort. An algorithmic procedure for automatic selection of a suitable degeneracy set is presently under development.
In order to select the degeneracy sets we use the entanglement structure for guidance. Obviously, for highly entangled states one needs large degeneracy sets with many virtual spins. In order to estimate the entanglement of the states at different θ, we use exact diagonalization for a system of N = 10 spins and determined the negativity [53] of the ground state to quantify nearest-neighbor qutrit-qutrit entanglement (see Fig. 4). The negativity N is obtained from the reduced 2-qutrit density matrix ρ, where T 1 denotes the partial transpose of ρ with respect to the first qutrit and . the trace norm. It is obvious from Fig. 4 that bipartite entanglement is largest in the dimerized phase, while in the Haldane phase and the critical phase the states are less strongly bipartite entangled. Therefore, we cautiously expect that the largest computational resources will be needed in the dimerized phase. At the AKLT point θ = arctan 1 3 there is a local entanglement minimum. At this point we just need the trivial virtual spin representation {(1/2, 1)} for a quasi-exact calculation for any system size. The bi-partite entanglement vanishes at the point θ = π/2 (which happens even for large systems). At θ = −3π/4 the negativity vanishes as well, however, at this point the ground state is in fact highly entangled. It resembles a state similar to the Greenberger-Horne-Zeilinger (GHZ) state, which is maximally entangled but measures zero nearest-neighbor entanglement.
Our choice for the virtual spin representations and their degeneracies is partly guided by these entanglement results. The negativity shows characteristic singularities at various points, foreshadowing clearly already for N = 10 the phase structure one observes in the thermodynamic limit. Interestingly, there is a weak singularity at the Heisenberg point θ = 0, where there is no phase transition in the thermodynamic limit.
Together with the bipartite entanglement we also plot in Fig. 4 the quantity J/(2N ) (J is the total spin of the ground state) for N = 10. In the parameter region of the critical phase J = 1, while it is zero in all other phases except the ferromagnetic phase.
The MPSs we construct are eigenstates of the BBH Hamiltonian as well as eigenstates of the total angular momentum operator (therefore we label the states with angular momentum J). However, they are not necessarily eigenstates of quasi-momentum operator as well. Nevertheless, if the constructed eigenstates are non-degenerate (apart from the SU(2) degeneracies), we determine in some cases the quasi-momentum of the constructed states and label the states and their energies with the momentum label p. The operator T n of the translation over n sites acts on an eigenstate with well-defined momentum p as T n |ψ = e −ipn |ψ .
Therefore, for such states p = 2πn p /N (n p = 0, . . . , N − 1) can be determined from the expectation value of the translation operator, ψ|T n |ψ = e −ipn , which is easily calculated in the MPS framework using Eq. (4) and a cyclic shift of the MPS tensors in |ψ .
A. Ground state energy Fig. 5 shows results for the ground state energies per site in the J = 0 and J = N sectors for N = 100 spins as a function of the parameter θ. The 'global' ground state, i.e. the lowest state of all spin sectors, is J = 0 except in the critical phase ( 1 4 < θ/π < 1 2 ), where the 'global' ground state is J = 1 for N = 100, and the ferromagnetic region. Due to the symmetry H(θ) = −H(θ + π) the minimal and maximal eigen-energies are related by E min (θ) = −E max (θ +π) (this maximal energy is also shown in Fig. 5). With few exceptions we used the representation {(1/2, 6), (3/2, 6), (5/2, 5)} in order to produce this plot. Note that results for N = 100 are rather close to the thermodynamic limit.
At θ/π = − 3 4 , − 1 2 , − 1 4 , 1 4 , 1 2 exact ground state energies per site are known from Bethe Ansatz calculations. We compare with these results in Table I and find rather good agreement. It is interesting that our numerical results are obtained using half-integer virtual spins only (we obtain similar values with integer virtual spins, however, at a somewhat larger numerical cost). Precise calculations are most resource intensive for the two SU(3) symmetric points, θ = π/4 and θ = −3π/4. In particular, virtual spins j = 7/2 are necessary in the degeneracy set d in the vicinity of θ = −3π/4. E T are finite-size predictions for θ/π = 0.5, −0.75, −0.5 (e.g., [38]). Infinite-size Bethe ansatz results for θ/π = 0.25 [35,47]   In the ferromagnetic region 1 2 < θ/π < 5 4 the bilinearbiquadratic spin-1 Heisenberg model has a 2N + 1 degenerate spin-N ground state for any system size N . It is possible to work out the ground state energy analytically within the proposed MPS algorithm as the virtual spin configuration is simply {(N/2, 1)}. We demonstrate this in Appendix G.

B. Low lying spectrum
While it is possible to extract much information on the phase structure from an analysis of the ground state properties, in the present paper we will concentrate on the low lying excitation spectrum, and we will discuss the spectrum at selected points in the phase diagram. We will emphasize characteristic differences of the low lying spectrum in different phases which show up already for finite systems.
The spectrum for a system of 50 sites within the range − 3 4 < θ < 1 2 is presented in Fig. 6. We include the lowest two spin-0 states, one or two spin-1 states and one spin-2 state. The states are labelled by their momentum p in some cases. We discuss the spectra in the various phases   shown in Fig. 6 from right to left: Characteristically, in the critical phase ( 1 4 < θ/π < 1 2 ) the ground state for 50 spins is doubly degenerate spin-1. There is small gap to the lowest spin-2 state, and there are slightly larger gaps to two lowest spin-0 states. Consequently, the lowest excitation is expected to be quadrupolar in this phase, which confirms predictions for spin chains in [40]. Our numerical results for the ground state for the systems with 50 or 100 sites at θ = 0.4π agree quantitatively very well with predictions in Ref. [54] obtained by extrapolation of the results for relatively small systems.
The ground state in the entire Haldane phase (− 1 4 < θ/π < 1 4 ) is J = 0, and the lowest excited state is J = 1. The gap between these two states remains finite in the thermodynamic limit (Haldane gap). Numerical data for the J = 0 state and J = 1 state for a purely bi-linear Heisenberg ring with N = 100 sites are given in Table II. From these results one obtains the Haldane gap ∆ = 0.41096 as compared to ∆ = 0.41048 given in Ref. [22,55]. At the given precision N = 100 cannot be distinguished from infinite systems. The calculated energy of the J = 0 state agrees well with the results from previous studies [56,57]. The data presented in Fig. 6 refer to 50 sites. It is observed that at θ = 0 the first excited state in the spin-0 sector lies much above the lowest spin-1 state, and it is nearly degenerate with the lowest spin-2 state (see also [55]).
The dimerized phase (− 3 4 < θ/π < − 1 4 ) has been the subject of many investigations in the history of the BBH model, and we will address the calculation of the dimerization in the next subsection. This phase is characterized by the fact that the first excited spin-0 state is lower than the lowest spin-1 and spin-2 states. The spin-1 state is lower than the spin-2 state for θ/π > − 1 2 but above this state for θ/π < − 1 2 . At θ = −π/2 the lowest spin-1 and spin-2 states form a degenerate pair (see also [36]). In the thermodynamic limit, the two lowest J = 0 states are degenerate, however with a finite gap to the lowest J = 1 or J = 2 states. As shown in Table III the calculated energies of the three lowest states at the biquadratic point θ = −π/2 agree very well with the Bethe Ansatz calculations of Sørensen and Young [38]. According to the heuristic entanglement analysis presented in the introduction to this section we expect calculations to be   [38].
rather "hard" at this point. And indeed, we need large degeneracy sets in order to achieve reasonable precision.
In the parameter region of the debated nematic phase close to θ = −3π/4 PBC rings of N ≤ 16 sites were investigated in [42]. We extend this work for systems of up to 50 sites. We observe that the lowest spin-2 state is significantly lower than the first excited spin-0 state as well as the lowest spin-1 state. As a consequence, quadrupolar correlations increase strongly in this region as pointed out in [40].
One may be tempted to conjecture the existence of the separate phase in this region from our results. However, in order to do so one must study the behaviour of the spectrum in the thermodynamic limit, which we do using finite-size scaling analysis for the gaps ∆ 00 and ∆ 02 . The conjecture of Ref. [41] about the existence of the gapped nematic phase (i.e., ∆ ∞ 00 > 0 and ∆ ∞ 02 = 0) was rejected in earlier works [16,42,43]. In fact, these calculations suggest that ∆ ∞ 00 = 0. However, there exists still the possibility that ∆ ∞ 02 = 0 in a narrow region of − 3 4 < θ/π −0.7. Earlier DMRG calculations for OBC [40] were not able to resolve extremely small values of ∆ ∞ 02 . From a finite size scaling analysis we obtain ∆ ∞ 00 2 · 10 −4 at θ/π = −0.72. This suggests that the gap between the two lowest spin-0 states closes in agreement with earlier calculations. In fact, ∆ 00 (N ) can be fitted very well by a power law ∆ 00 (N ) = BN −α . In line with this, the dimerization remains finite, D ∞ = (4 ± 2)·10 −3 , as discussed in more detail in the following subsection. In contrast, ∆ 02 (N ) shows exponential behavior. From a fit to our data we obtain a gap between the lowest spin-0 and the spin-2 states of size ∆ ∞ 02 ≃ 0.009 at θ = −0.72π. Moreover, we observe a monotonic decrease of the scaled gap N ∆ 00 (N ) and a monotonic increase of the scaled gap N ∆ 02 (N ) with increasing of N , which also indicates the absence of a nematic state at this θ. However, as pointed out in [44], the possibility of the existence of such a state even closer to θ/π = − 3 4 cannot be excluded. Calculations closer to the ferromagnetic region for larger systems are presently under way. for the bilinear-biquadratic Heisenberg ring of N = 50 sites as a function of θ. D(N ) ≃ 0.6015 at θ = −π/2. The theoretical prediction for the infinite system at θ = −π/2 is indicated by a cross [58,59], and our extrapolated value D∞ ≃ 0.568 is in good agreement with theory.

C. Correlation functions
Finally, we also calculate a few physically interesting correlations in the BBH model. In the area −3/4 < θ/π < −1/4 we determine the dimerization correlator shown for a system of N = 50 sites in Fig. 7. In fact, due to translational symmetry, the ground state as well as the excited states of a finite-size PBC ring are not dimerized. In order to calculate the dimerization we take a symmetric/antisymmetric superposition of the lowest two J = 0 states with different momenta. These two states are separated by a very small gap for finite systems, and they should develop into the degenerate doublet in the thermodynamic limit. The dimer order parameter calculated in this way is The dimerization D can be calculated analytically for θ = −π/2 in the thermodynamic limit, where the BBH model can be mapped to a spin-1/2 XXZ model. According to Refs. [58,59] n=1 tanh 2 (n arccosh 3 2 ) ≃ 0.5622, which agrees with the Monte-Carlo result of Ref. [38].
We fitted our results for 30, 40 and 50 sites at θ = −π/2 by a function given in Eq. (6) of [44]: The fit gives D ∞ ≃ 0.568 in good agreement with the theory, and ξ ≃ 20.2. Calculations for larger systems (about 100 sites) are needed to obtain a more precise result.
The dimer order parameter (22) was calculated in [42] for finite chains of up to 48 sites. In such OBC calculations there is no translational symmetry, and the degenerate doublet is mixed automatically. For systems of up to 50 sites the results for PBC are always slightly larger than for OBC. There are also calculations of the dimerization for spin rings in Ref. [18], however for the four-point correlator 0 (0) |D 2 |0 (0) . Our extrapolated result for θ = −0.65π is consistent with the result of [18].
Finally we come back again to the Haldane phase: This phase is characterized by the presence of a nonzero string correlator of the ground spin-0 state [60], that does not decay in the limit N, l → ∞. Fig. 8 shows the string correlator g(l) as a function of θ. The system size is N = 100, and the string length is taken as l = 30. This choice allows sufficient length of both the string and the rest of the system, in order to resemble the infinitesystem properties (numerical results indicate that g(l) is constant to a high degree for 20 ≤ l ≤ 40). Our results in the range −0.2 ≤ θ/π ≤ 0.2 agree quantitatively very well with infinite-size calculations presented in [61].
Our numerical calculations indicate that long-ranged string order is present in the Haldane phase only (however, it is numerically hard to obtain precise values at the boundaries of the Haldane phase). The string correlator for the AKLT point (θ = arctan 1/3) can be calculated exactly [32]: g(∞) = 4/9, which is reproduced extremely well in our calculation for N = 100 sites. The curve has the maximum at this point. The best estimate for a string correlator for the bilinear Heisenberg spin-1 model (θ = 0) is g(∞) ≃ 0.374325 [2]. Our result g(30) = 0.374330 calculated for a system size of N = 100 spins agrees very well with this value.

V. CONCLUSIONS
In the present work we developed an algorithm for SU(2) symmetric matrix product states (MPS) with periodic boundary conditions (PBC) that involves only reduced tensor elements. It was applied to a study of the lowest-lying states of the spectrum of the spin-1 bilinearbiquadratic Heisenberg model of up to 100 sites. The characteristic differences of the spectrum in the various phases of the model were stressed. Dimerization and string order were studied in the dimerized and Haldane phase, respectively. Our results agree rather well with previous studies based on DMRG or Bethe Ansatz calculations, and we could extend previous studies to a more complete coverage of the full parameter range, in particular in the dimerized phase close to the transition to the ferromagnetic phase. We confirm the absence of a nematic phase for θ/π > −0.72 within our numerical precision.
The precision of the results we can achieve with our implementation of the algorithm depends on the degeneracy set d we choose for the representation of the virtual spin indices of the MPS. Due to fact that we have eliminated all spin projection indices, our algorithm achieves rather large 'effective' virtual spin dimensions in comparison to a non-symmetric implementation. E.g. a state with reduced MPS dimension m = 30 described by the degeneracy set d = {(0, 6), (1/2, 6), (1,6), (3/2, 6), (2, 6)} corresponds to an effective MPS virtual dimension m = 90, which is already rather large for a PBC calculation. Moreover, due to the sparse structure of each MPS reduced tensor, such a state is described by only 288 complex parameters per site compared to 16200 complex parameters per site necessary to specify a non-symmetric state. The number of parameters is also roughly one order smaller than for a U(1) symmetric state. For this state about 98% of the parameters used in a nonsymmetric calculation are actually zero due to symmetry. As a consequence, the numerical effort for the solution of the generalized eigenvalue problem at each update step is reduced correspondingly. In fact, we suspect that the states determined without explicit consideration of the symmetry are incorrect even though their energy is determined precisely. This needs to be investigated in more detail in future work.
Of course, the numerical effort to be expected grows with the size of the chosen degeneracy set d. And at this stage this set is fixed from the start of our calculation. It would be desirable that the algorithm dynamically chooses this set according to suitable algorithmic criteria. The implementation of such a procedure is under way using ideas presented in Ref. [56] and discussed for PBC in Ref. [51].

ACKNOWLEDGMENTS
We thank Briiissuurs Braiorr-Orrs for discussions. Mykhailo V. Rakov thanks Physikalisch-Technische Bundesanstalt for financial support during short visits to Braunschweig.

Appendix A: Reduced representation of MPOs
In Eq. (19) we have defined the W tensors of MPO in terms of reduced tensors W. Here, we will present explicit results for the reduced tensors needed in the present study. In fact, all formulas in this Appendix are valid for arbitrary spin-s.
In order to construct the reduced tensor W H for the spin-s BBH Hamiltonian, we rewrite this Hamiltonian in terms of the tensors C s,je,s ms,me,m ′ s , i.e. replace the spin matrices by C tensors. In general, we would expect terms with j e = 0, 1, · · · , 2s− 1, 2s, however one only finds nonzero terms with j e ≤ 2 for this Hamiltonian. The corresponding W H tensors have dimension 11×11×s×s, which can be grouped into three SU(2) singlets, one triplet and one quintet. The reduced tensor W H is therefore labeled by the degeneracy set d = {(0, 3), (1, 1), (2, 1)}. One then finds for the first site i = 1 the W H matrices, And for all other sites i = 2, · · · , N one obtains and Σ 0 = Sgn(sin θ) and Σ 1 = Sgn(sin θ − 2 cos θ). It follows that for s = 1/2 we have w 2 = 0 and W 2 = 0. Due to PBC the tensors at the first site have a different structure than the tensors at the other sites. The reduced tensors of the corresponding unity MPOs are given by i.e. a diagonal matrix constructed from all j's of the degeneracy set d. E.g. for the set d = {(0, 3), (1, 1), (2, 1)} the reduced unity operator is given by This matrix is used at the site of the fictitious spin. For the Heisenberg antiferromagnet (θ = 0) or ferromagnet (θ = π) the MPOs can be simplified. One obtains: w 0 = w 2 = 0, Σ 0 = 0, Σ 1 = −Sgn(cos θ) = ∓1, and W 2 = 0. Therefore, the quintet and one singlet in the basis can be omitted, so that the basis is {(0, 2), (1, 1)}. Then, the reduced tensors take the form for i = 1 with ω 1 = 3s(s + 1). The reduced unity tensor at the site of the fictitious spin shrinks to

Appendix B: Reduced block tensors
The calculation of matrix elements of MPS in MPO Eq. (4) reduces to a finite product of transfer tensors ('blocks'). Similarly, the block tensors H L , H R , N L , N R , O L,k , O R,k which are necessary to build the generalized eigenvalue problem (7), are finite products of transfer operators. The size of such blocks may be increased by multiplication of a transfer tensor on the right or on the left of an existing block. Generically, we will denote these block tensors by B L and B R , and the index distinguishes if the block is produced by left or the right multiplication.
In reduced form these tensors have three intermediate indices e 1 , e 2 , e 3 , and, therefore, a reduced block tensor B L has rank 9. To define these reduced tensors we choose the coupling scheme shown in Fig. 9, which corresponds to the following expression We can take t e1 = t e2 = t e3 = 1 without loss of generality. For N -and O-blocks jĪ = j I = 0 (and γĪ = γ I = 1) and, consequently, j e1 = j e2 = j e3 . For right blocks a similar expression may be written down. Using the orthogonality relations of the 3j symbols [14], one can express the reduced tensor elements of the reduced block B L in terms of the tensor elements of B L .
During the update procedure of the algorithm the length of the blocks has to be increased by right or left multiplication of transfer tensors as illustrated in Fig. 2. However, we only need to determine reduced blocks and want to avoid reconstruction of the full block. In order to do that the summation over intermediate spin projection indices must be carried out analytically using appropriate summation formulas for 3j symbols. Such formulas can be found e.g. in Edmonds [14], and we will strictly follow the conventions of that reference.
An increase of the reduced left block by one transfer matrix may be formulated in terms of the following recursion formula,

(B
[i+1],je 1 ,je 2 ,je 3 L ) (γĪ ,γā,γb),(γI ,γa,γ b ) = (2j e2 + 1) The expressions in curly brackets denote the Wigner 9j and the Racah 6j symbol, respectively. The matrices α and β are the reduced tensors corresponding to the MPS tensors A and B, respectively, and W is the reduced W tensor of an MPO. The local spin s i = s for i ≤ N and J for the fictitious site. In order to avoid unnecessary indices we again do not label this expression with the local spin index, which is fixed in our calculation. The recursion is started from a reduced identity block tensor, which is easily obtained from Eq. (B1): (γĪ ,γā,γb),(γI ,γa,γ b ) = (2j e1 + 1)(2j e3 + 1) δ γĪ ,γI δ γā,γa δ γb,γ b δ je 1 ,je 2 δ(jā, j e1 , jb) δ(jĪ , j e3 , j e1 ) (B5) with δ(j 1 , j 2 , j 3 ) = 1 if j 1 , j 2 , j 3 fulfill the 'triangle rule' and zero otherwise. Superficially it appears that the recursion formula for B L is (appart from a trivial prefactor) independent of j 1 and analogously B R independent of j 2 . However, this is in fact not the case since both j 1 and j 2 enter via the initial condition of the recursion. We remark that the reduced block tensors B are very sparse. Therefore, an implementation of arbitrary rank sparse tensors is needed for their calculation and manipulation.
In the same way one also obtains the reduced vector y k for the calculation of excited states according to Eq. (12), where the index k labels the different excited states, k denotes the reduced tensor of the MPS tensor Φ [i] k in Eq. (13). Analogously, the block tensors B that correspond here to O tensors are labeled by the additional index k. and the two right virtual indices are fused into a single virtual index each, and a summation over intermediate spin projections of the physical spin s is performed.
In the SU(2) symmetric formalism, we need to express the reduced W tensor of the squared operator O 2 in terms of the reduced W tensor of the operator O. To this end we follow section 7.1 of Edmonds [14]. Using Eq. (7.1.1) and Eq. (7.1.5) of this reference one obtains = (2j e + 1) (2j + 1)(2j ′ + 1)(2s + 1) (−1) Here, j e1 and j e2 are the intermediate indices of reduced W O tensors while j e is the intermediate index of reduced W O 2 tensor. The left virtual indices (j 1 t 1 ) and (j 2 t 2 ) of W O tensors are fused into one index (jt), and, analogously, the right indices (j ′ 1 t ′ 1 ) and (j ′ 2 t ′ 2 ) are fused into the index (j ′ t ′ ). One can immediately read off Eq. (D1) that for the bilinear-biquadratic spin-1 Heisenberg model the reduced W H 2 tensors have dimensions 3 × 36 × 36 with the virtual spins described by the degeneracy set {(0, 11), (1,11), (2,10), (3,3), (4, 1)}. These tensors are rather sparse. Note, that the result (D1) holds for any spin s. The matrix H eff depends on the Hamiltonian under consideration. For the spin-1 bilinear-biquadratic Heisenberg model one finds (H eff ) (1,1),(1,1) = N (cos θ+sin θ) Consequently, one obtains the energy of the J = N state as E = 1 N (H eff ) (1,1), (1,1) (N eff ) (1,1),(1,1) = cos θ + sin θ as expected. Despite the analytical simplicity of this result it cannot be obtained numerically within the present implementation. As is obvious from the analytical results the matrix elements of matrices H eff and N eff decrease with increasing N . At about N ≃ 25 the numerical values cannot be represented by machine precision numbers any more, and the calculated values are arbitrary.