Simulation of one-dimensional quantum systems with a global SU(2) symmetry

In this paper, we describe a refined matrix product representation for many-body states that are invariant under SU(2) transformations and use it to extend the time-evolving block decimation (TEBD) algorithm to the simulation of time evolution in the presence of an SU(2) symmetry. The resulting algorithm, when tested in a critical quantum spin chain, proved to be more efficient than the standard TEBD.


Introduction
Quantum many-body systems are characterized by a large Hilbert space, one whose dimension grows exponentially with the system's size. This makes the numerical study of generic quantum many-body phenomena computationally difficult. However, quantum systems are governed by Hamiltonians made of local interactions, that is, by highly non-generic operators. As a result, physically relevant states are atypical vectors in the Hilbert space and can sometimes be described efficiently. Systems in one spatial dimension are a prominent example. Here, the geometry of local interactions induces an anomalously small amount of bipartite correlations and an effective representation is often possible in terms of a trial wave function known as the matrix product state (MPS) [1,2]. This, in turn, underlies the success of the density matrix renormalization group (DMRG) [3], which is an algorithm to compute ground states, and of several recent extensions of it [4,5], including the time-evolving block decimation (TEBD) algorithm to simulate time evolution [4].
Symmetries, of fundamental importance in physics, require special treatment in numerical studies. Unless explicitly preserved at the algorithmic level, they are bound to be destroyed by the accumulation of small errors, in which case significant features of the system might be concealed. On the other hand, when properly handled, the presence of a symmetry can be exploited to reduce simulation costs. The latter has long been realized in the context of DMRG [3,6]. More recently, the Abelian U(1) symmetry has been incorporated into the TEBD algorithm [7].
In this paper, we study how to enhance the TEBD algorithm in systems that are invariant under the action of a more general (possibly non-Abelian) symmetry, as given by a compact, completely reducible group G. (This includes finite groups as well as Lie groups such as O(n), SO(n), U(n) and SU(n).) We present an explicit theoretical construction of a refined MPS representation with built-in symmetry and put forward a significantly faster TEBD algorithm that both preserves and exploits the symmetry. For simplicity and concreteness, we analyze the smallest non-Abelian case, the SU(2) group, which is relevant in the context of isotropic quantum spin systems. In contrast to the much simpler case of an Abelian U(1) symmetry [7], the analysis of the SU(2) symmetry already contains the major ingredients of an arbitrary compact, completely reducible group G. In addition, it can be cast in the language of spin operators, which is often more familiar to physicists when compared with group representation theory. As a test, we have computed the ground state of the spin-1/2 antiferromagnetic Heisenberg chain, obtaining remarkably precise two-point correlators for both short and long distances.
In preparation for describing the SU(2) MPS, we start by introducing a convenient vector basis and discuss the bipartite decomposition of states that are invariant under SU(2).

Total spin basis
Let V be a vector space on which SU(2) acts unitarily by means of transformations e i v· S , where matrices S x , S y and S z close the Lie algebra SU(2), namely [S α , S β ] = i αβγ S γ , and v ∈ R 3 . A total spin basis (TSB) | [V ] jtm ∈ V satisfies the eigenvalue relations 3 and is associated with the direct sum decomposition of V into irreducible representations (irreps) of SU(2) [8], Here,Ṽ ( j) is a d j -dimensional space that accounts for the degeneracy of the spin-j irrep and has basis where m is the projection of the spin in the z-direction, m = − j, . . . , j. Each vector of the TSB factorizes into degeneracy and irrep parts as

Bipartite decomposition
Consider an operator : A → B that acts from vector space A to B and that commutes with the action of SU (2). Let S [A] and S [B] generate the action of the group on A and B. Then we have Describing A and B with a TSB and employing Schur's lemma [8], splits into degeneracy and irrep parts as The operator is dual to a pure state | of a bipartite system with vector space A ⊗ B by taking [B] where ω is completely determined in terms of the Clebsch-Gordan coefficients j 1 j 2 m 1 m 2 | j 1 j 2 ; jm (see, for instance, [10]), namely and where | is an SU(2) singlet, that is, invariant under transformations acting simultaneously on A and B, or Under the transformation in equation (5), equation (4) becomes New Journal of Physics 12 (2010) 033029 (http://www.njp.org/)

4
The above decomposition is quite sensible: it shows that a coefficient N j 1 t 1 m 1 j 2 t 2 m 2 may be nonzero only if (i) j 1 = j 2 (only the product of two spin j irreps can give rise to a spin 0 irrep, that is, the singlet | ) and (ii) m 1 = −m 2 , which guarantees that the z-component of the spin vanishes.
In addition, equation (9) embodies the essence of our strategy: to isolate the degrees of freedom that are not determined by the symmetry-in this case the degeneracy tensor T j t 1 t 2 . We now consider the singular value decomposition (SVD) of the tensor T j t 1 t 2 for a fixed j, and define By combining equations (9)- (11), we arrive at our canonical symmetric bipartite decomposition (CSBD) which is related to the Schmidt decomposition by the identifications α → ( jtm), λ α → η j t ω j m and where some of the Schmidt coefficients λ α are negative. (If one prefers to work with nonnegative Schmidt coefficients, negative signs can be absorbed into e.g.
jtm of a bipartite system C ⊗ D can be expressed in terms of TSBs for C and D as [9] where tensor X relates to degeneracy degrees of freedom and tensor C is given by the Clebsch-Gordan coefficients

Matrix product decomposition
We now consider a chain of n quantum spins with spin s, represented by a one-dimensional (1D) lattice where each site, labeled r (r = 1, . . . , n), carries a (2s + 1)-dimensional irrep of SU (2). where m , can be codified as a MPS [1,2], Following the conventions of Vidal's paper [4], here λ [r ] α are the Schmidt coefficients of | according to the bipartition [1 · · · r ] : [r + 1 · · · n] of the spin chain and the tensor [r ]m αβ relates the Schmidt vectors for consecutive bipartitions, When | is a singlet, that is, equations (12) and (15) supersede equations (13) and (19) and each tensor λ and in equation (18) decomposes into degeneracy and irrep parts, see figure 1, whereC is related to the Clebsch-Gordan coefficients C bỹ Here, the term (ω j m ) −1 arises in conjunction with the transformation equation (5). The SU (2) MPS is defined by equations (21) and (22). In this representation, the constraints imposed by the symmetry are used to our advantage. By splitting tensors λ and , we achieve The TEBD algorithm is based on updating the MPS when a gate U acts on two neighboring sites. This diagram generalizes figure 3(i) of [11] after the replacements λ → (η, ω) and → (X,C) of equations (21) and (22) for a SU (2) MPS.
two goals simultaneously. On the one hand, the resulting MPS is guaranteed, by construction, to be invariant under SU (2) transformations. That is, any algorithm based on this representation will preserve the symmetry exactly and permanently. On the other hand, all the degrees of freedom of | are concentrated in smaller tensors η and X (tensors ω andC are specified by the symmetry), and thus the SU (2) MPS is a more economical representation. If | · | denotes the number of coefficients of a tensor, then where λ and λ are the tensors to the left and to the right of and where, following spin composition rules, the last sum is restricted to pairs ( j, j ) such that | j − j | s.

Simulation of time evolution
Our next step is to generalize the TEBD algorithm [4] to the simulation of SU(2)-invariant time evolution. This reduces to explaining how to update the SU (2) MPS when an SU(2)invariant gate U acts between contiguous sites, see figure 2. The update is achieved by following steps analogous to those of the regular TEBD algorithm, see figure 3 of [11], involving tensor multiplications and one SVD (figure 3). Following the procedure described in [11], we first absorb the gate into the tensors of the MPS by contracting the tensor network in figure 3.
Here V Xi and V Ci enact a change of basis from the product basis {| j 1 m 1 t 1 } ⊗ | j p m p t p to the coupled basis {| jmt }. Specifically, V Ci are the Clebsch-Gordan coefficients that describe the change of basis | j 1 m 1 ⊗ | j p m p → | jm , and V Xi denotes the change of basis | j 1 t 1 ⊗ | j p t p → | jt . A straightforward choice for V Xi is an invertible linear map that uniquely associates each pair j, t to a tuple j 1 , t 1 , j p , t p . Define (V Xi ) jt j 1 t 1 j p t p = 1 for each such unique association and as 0 otherwise. The network contraction in figure 3 is a sum of the terms corresponding to various compatible values of j 1 , j 2 and j 3 . For each such term, the top line of tensors comprising the η's, X 's and V X 's are multiplied together. By Schur's lemma, the product of the tensors ω's, C's, V C 's and U on the lower leg is proportional to ω, with a proportionality constant α that depends on j, j 1 , j 2 , j 3 and j p and that can be pre-computed because none of these tensors depend on | . α is then absorbed into the product of the degeneracy matrices, which are then  (2) MPS, analogous to figures 3(i)-(iii) of [11] for a regular MPS. Once U has been applied to two spins, additional tensors V Xi and V Ci implement a unitary transformation required to reabsorb these spins into blocks and obtain an updated representation for the bipartition [1 · · · r ] : [r + 1 · · · n]. Then, for each fixed value of the j indices (discontinuous lines), the η's, X 's and V X 's are multiplied together and the result, with a weight α j 1 j 2 j 3 j coming from the product of the ω's,C's, V C 's and U (that can be pre-computed because none of these tensors depend on | ), is added together to give rise to tensor T .
added together to give rise to tensor T . A SVD of the tensor T jt jt ensues (equation (10)), giving rise to tensors R and S and the updated tensor η [r ] . The coupled basis | jmt is reorganized back into product basis by contracting R and S with (V Xi ) −1 (shown in figure 4(B) for R only). Once again, the contraction is effected by multiplying (V Xi ) −1 with R (or S) for each value of j. Finally, contractR (andS) with (η [r −1] ) −1 (and (η [r +1] ) −1 ) to obtain the updated tensors X [r ] and X [r +1] . Note that all matrix multiplications and SVDs now involve smaller tensors, and only tensors X and η of the SU(2) MPS need to be updated. This results in a substantial reduction in computational space and time and thus an increase in performance. For instance, the SVD of in figure 3 of [11], where | | ≈ (2s + 1) 2 |λ| 2 is now replaced with the SVD of jt jt (see figure 3) for each value of j, where | jt jt | = ( j+2s j j−2s d j ) 2 . The cost c svd (A) of computing the SVD of a matrix A grows roughly as |A| 3/2 and is the most expensive manipulation of the TEBD algorithm. We obtain the following comparative costs:

Example
For illustrative purposes, we consider a quantum spin chain with s = 1/2 and with the Hamiltonian that is, the spin-1/2 antiferromagnetic Heisenberg model, which is SU(2) invariant and quantum critical at zero temperature. Representing the ground state of this model with a MPS is 9 challenging, because both the critical character of the system and the presence of the SU(2) symmetry contribute to its entanglement. As a result, compared to a gapped system or a system with smaller symmetry, an accurate MPS approximation to the ground state of (28) requires a larger value of χ. We have computed a SU(2) MPS approximation to the ground state of H , in the limit n → ∞ of an infinite chain, by simulating imaginary-time evolution [11] starting from a state made of nearest-neighbor singlets ( With the constraint j d j = 600, we have obtained that the following irreps j, with degeneracies d j , contribute to the odd and even bipartitions 4 of the resulting state: that is, with a regular MPS, storing the same state would require about 50 times more computer memory, whereas performing each SVD would be about 300 times slower.
We have computed the two-point correlators 5 z , and C 2 (r ) ≡ S [1] z S [r +1] z , and the average C 2 (r ) ≡ (C 2 (r ) + C 2 (r ))/2. For small r , they read as follows: where, for C 2 (r ), the square brackets show the first digits that differ from the exact solution ( [12] and reference therein), from which we recover e.g. nine significant digits for r = 1. An expression for the correlator C 2 (r ) is also known for large r 6 . There, for r ≈ 4000, 10 000 and 13 000, our results approximate the asymptotic solution with errors of 1, 5 and 10%, respectively, see figure 5. In comparison, with a regular MPS and similar computational resources, we lose three digits of precision for r = 1, whereas a 10% error is already achieved for r ≈ 500 instead of r ≈ 13 000. It is important to note that the appearance of two different correlators C 2 (r ) and C 2 (r ) is an artifact of the present algorithm and that the average C 2 (r ) ≡ (C 2 (r ) + C 2 (r ))/2 is a legitimate procedure to eliminate it, since it corresponds to taking as the ground state of the infinite chain an equally weighted superposition of the MPS and its translation by one site. It is also of interest to note the effect that a finite χ in the MPS has on symmetries. If the infinite Heisenberg chain is 4 For semi-integer spins, e.g. s = 1 2 , the CSBD of equation (12) for a partition [1 · · · r ] : [r + 1 · · · n] of the chain has only integer (semi-integer) values of j when r is even (odd). 5 C w 2 (1) (w = , ∇) are computed by contracting a small tensor network involving the tensors (η, ω) and (X,C) for two contiguous sites. For r > 1, C w 2 (r ) is computed in the same way, but first we simulate r − 1 (SU(2) invariant) swap gates that bring the two relevant sites together. 6 We use c = 1.05 in equation (5.25) of Lukganov and Terras [13].  (24), that is, the rank of an equivalent (regular) MPS. The lowest line, χ 6 = 2200, shows the errors in the data presented in the table. Bottom: numerical results for C 2 (r ) for up to r = 20 000 sites, for different sizes χ i , together with corresponding errors i .
addressed with a standard MPS, translation invariance is largely preserved but SU (2) invariance is clearly broken. If a SU (2) MPS is used instead, SU(2) invariance is exactly preserved but translation invariance is clearly broken. The extent to which either symmetry is broken decreases with increasing χ.

Final remarks
The above test with a critical spin-1/2 chain unambiguously demonstrates the superiority of the SU(2) MPS and TEBD with respect to their non-symmetric versions. Promisingly, these techniques can now be used to address systems that remain otherwise largely inaccessible to numerical analysis due to the large dimension of the local Hilbert space. These systems include a chain made of large spins, say s = 4, or a spin ladder with several legs. We regard the latter as a chain with several spins per site, where each site decomposes into SU(2) irreps as in equation (2) [9]. In addition, the SU(2) MPS is not restricted to the representation of SU(2) singlets. On the one hand, it can be used to represent any SU(2) invariant mixed state ρ of the chain, which decomposes as (see equation (2)) This is achieved by attaching, to the end of the chain, an environment E that duplicates the subspace V of the chain on which ρ is supported and by considering a singlet purification We first build a SU(2) MPS for the purification and then trace out E. The resulting structure is a matrix product representation that retains the advantages of the SU(2) MPS. In particular, note that when ρ corresponds to a single irrep j, the environment is a site with a spin j, and the chain together with the environment is just an extended spin chain, with the purification being of the form On the other hand, the SU (2) MPS can also be modified to represent any pure state | [V ] jm of the chain with well defined j and m. To achieve this, we first consider a mixed state ρ as in equation (31), that is, a symmetrization of | [V ] jm , and then consider a purification | ρ for ρ as in equation (31), for which we can build a SU (2) MPS. Finally, we recall that | [V ] jm = [E] jm | ρ , which leads to a simple, SU(2) MPS-like representation for | [V ] jm in terms of the SU(2) MPS for the purification | ρ . The time-evolution simulation techniques described in this paper can be applied to the above generalized representations.
After nearly completing the original manuscript, we became aware of the related results obtained by I McCulloch, which were derived independently in the context of DMRG [14].