# Matrix Product States with dynamically emerging global symmetries

###### Abstract

Quantum many body physics simulations with Matrix Product States can often be accelerated if the quantum symmetries present in the system are explicitly taken into account. Conventionally, quantum symmetries have to be determined before hand when constructing the tensors for the Matrix Product States algorithm. In this work, we present a Matrix Product States algorithm with a dynamical symmetry. This algorithm can take into account of, or benefit from, or symmetries when they are present, or analyze the non-symmetric scenario when the symmetries are broken without any external alteration of the code. To give some concrete examples we consider an XYZ model and show the insight that can be gained by (i) searching the ground state and (ii) evolving in time after a symmetry-changing quench. To show the generality of the method, we also consider an interacting bosonic system under the effect of a symmetry-breaking dissipation.

###### pacs:

## I introduction

Matrix Product States (MPS) algorithms have become one of the best tools to study one-dimensional quantum manybody systems FannesWerner1992 (); OstlundRommer1995 (); Perez-GarciaCirac2007 (); Schollwock2011 (). It has been shown that they can represent the ground state of a wide class of local Hamiltonians Hastings2007 (); Vidal2008 (). They can also be applied to study the time evolution of both unitary and open quantum manybody systems Vidal2003 (); Vidal2004 (); DaleyVidal2004 (); WhiteFeiguin2004 (); VerstraeteCirac2004 (); Daley2014 (); DeVegaBanuls2015 (); GuoPoletti2018a (); XuPoletti2019 (); ProsenZnidaric2009 (); PalmeroPoletti2019 (), provided the entanglement entropy, or the operator space entanglement entropy for dissipative systems PizornProsen2008 (), do not grow too fast.

Symmetries play a very important role in physics as they imply the presence of conserved quantities. These make it possible to write the Hamiltonian of a system in a block-diagonal form, where each block has an associated quantum number, thus significantly reducing the storage requirements and the computational cost of diagonalizing it. In a number of systems and models the total number of particles, or the total magnetization, are conserved: in these scenarios the systems have a symmetry. In other cases a system could have a discrete symmetry, like a symmetry, which is a subgroup of . In such situations, it is often advantageous to explicitly encode the quantum symmetry into the tensors of the MPS DelayVidal2004 (); Perez-GarciaVidal2008 (); SinghVidal2011 (), which are also called symmetric Matrix Product States (). would usually result in more efficient numerical algorithms since a much smaller number of states have to be considered. are also particularly important when searching for the ground state in a certain symmetry sector of an Hamiltonian, otherwise the ground state search algorithm would return the lowest ground state of the global Hilbert space, i.e. of all the symmetry sectors. Typically, one would thus use a different MPS depending on whether the Hamiltonian is symmetric or not, or maybe use a lesser efficient option. This is particularly true in the study of dynamical systems in which the equations of motion change the symmetry properties of the system during the evolution. For example, one could have prepared a quantum state in the ground state of an chain, which has symmetry, and then suddenly tune the parameter so that the system is modeled by a symmetric chain. With conventional , one would have to manually convert the ground state from to symmetry for the time evolution.

In this work we propose MPS with dynamical global symmetry to overcome this difficulty, which we will refer to as dynamically symmetric Matrix Product States (). The central idea of our approach is to generalize so that instead of a fixed total quantum number, it allows a superposition of different total quantum numbers. As a result, a is able to treat systems with symmetry or systems whose symmetry group is a subgroup of , or even no symmetry at all, on the same footing.

This paper is divided into the following sections: In Sec.II, we introduce in detail how to implement , including the ground state search and time evolution algorithms; In Sec.III, we show two concrete examples in which we can demonstrate the use of , one for a unitary system and the other for a dissipative system; In Sec.IV we draw our conclusions.

## Ii method

Similarly to , the building blocks of are symmetry protected tensors. The definition of a symmetry protected tensor, as well as the basic tensor operations based on it, has been presented in detail in the literature. One particularly detailed example is SinghVidal2011 (), although the use of symmetry protected tensors has been presented earlier, for instance in McCulloch2007 (). In Sec.II.1, we will first give a minimal introduction to symmetry protected tensors, then in Secs.II.2 and II.3 we, respectively, describe the dynamically symmetric MPS and Matrix Product Operator (MPO). In Sec.II.4 we describe how to search for the ground state and in Sec.II.5 how to implement time evolution.

### ii.1 Symmetry protected tensors

A symmetry protected tensor can be represented as a list of dense tensors labelled by quantum numbers, the dimension of the space corresponding to that quantum number and, very importantly, a direction for each quantum number stating whether it is incoming or outgoing from the tensor. This direction is key to preserve the symmetry for the overall state. For example, a dimensional symmetry protected tensor could be written as

Here is a dimensional dense tensor with dimensions , and for labels the four quantum number (note that here we don’t differentiate between superscripts and subscripts). The quantum numbers have two directions, and , depending on whether the quantum numbers flows out or in respectively. The symmetry can be explicitly enforced by the fusion rule

(1) |

where and are subtracted because their quantum numbers flow outside of the tensor , while and are added because their quantum numbers flow inside of the tensor . For the implementation of algorithms with symmetric Matrix Product States, it is useful to use an appropriate library that takes into account which tensors can be added, multiplied, decomposed and merged. Of particular importance would be a structured tensor which takes into account of the possible number conserving combinations (see, for example, SinghVidal2011 ()). To simplify the notation, in the following we will denote simply as and as , implicitly indicating the size ; the tensor in Eq.(II.1) would then be rewritten as .

### ii.2 Dynamically symmetric Matrix Product States

An MPS representing a quantum state on sites is a chain with dimensional tensors. The tensor on site of the MPS is labelled by a physical quantum number , an auxiliary quantum number which connects it to its left hand site, and which connects it to its right hand site. As a result, can be written as

and the fusion rule Eq.(1) becomes

(2) |

We encode the flow of quantum numbers as in Fig.1 (bottom). Conventionally, if an MPS has a fixed total quantum number , then one would have for the right boundary, and for the left boundary, and we could denote it as . The central observation is that one can generalize this approach such that still holds, while could have multiple values depending on how many possible total quantum numbers the quantum state has. As a result, different will correspond to a different total quantum number. Simply speaking, a is a collection of many s which we write symbolically as

(3) |

As for typical , the overall bond dimension (defined as the size of the auxiliary dimensions) of a is the sum of the sizes of all the blocks corresponding to different auxiliary quantum numbers. In this way, while keeping each block to a manageable size, it is possible to effectively use a very large overall bond dimension for the .

### ii.3 Dynamically symmetric Matrix Product Operator

Dynamically symmetric Matrix Product Operator () can be treated in a similar manner as they form a chain of dimensional symmetry protected tensors. The tensor on the -th site of a is label by two physical quantum numbers and , an auxiliary quantum number which connects it to its left hand site, and which connects it to its right hand site. Thus can be written as

where the dimensions of the tensors are implicit in this notation, and the quantum numbers satisfy

(4) |

The flow of quantum numbers is represented in Fig.1 (top). For symmetric systems, the MPO does not change the total quantum number of the quantum state. If the quantum state had total quantum number before being multiplied by the MPO, after the multiplication between MPO and MPS, the total quantum number for the quantum state is still . As a result we can denote a symmetric MPO as , with the particular number preserving condition that . In order to generalize this concept, we allow the quantum number at the left boundary to take multiple values while still keeping . Therefore the dynamically symmetric Matrix Product Operator can be written as

(5) |

where if and only if the system has symmetry. In case the system has only a discrete symmetry, for example parity, a symmetry, then and satisfy where is an even number. If the system has no symmetry at all, then and can be arbitrary. However, in this case, it may be more efficient to use non-symmetric and since the number of states will be the same for both approaches and the non-symmetric version has a much simpler data structure.

In the next sections we show how to adapt the ground state search and time evolution algorithms using and .

### ii.4 Ground state search algorithm

Here we present the ground state search algorithm that we have implemented based on . The core structure retraces that of two-sites variational Matrix Product States algorithm which we briefly present in the following. In this algorithm it is convenient to keep the MPS in the “mixed-canonical” form

(6) |

where the are called “left-canonical” because

and the are “right-canonical”, i.e.

In both Eqs.(II.4,II.4) represents the identity matrix for the two indices , . Here means to take the element wise conjugate of the tensor , during which the direction of each quantum number also gets reversed. The two tenors with a commom quantum number which is summed over should have opposite directions. No restrictions is made on the tensors . A two-site ground state search variational MPS algorithm works by iteratively minimizing the energy of an effective Hamiltonian on each neighbouring pair of sites. This minimization is done iteratively, from left to right and then back from right to left, in what is called a sweep. Thus the central step for the two-site ground state search algorithm is to build the local effective Hamiltonian . To compute on sites and (also referred to as bond ), one needs to trace out the left environment consisting of the sites to the left of the -th site to obtain the tensor

and the right environment consisting of sites to the right of the -th site to obtain the tensor

With , and , , one can construct as

Grouping the tensor indexes and , one could treat as a matrix and find its lowest eigenstate and eigenvalue. is then decomposed into two dimensional tensors by singular value decomposition (SVD)

and the results are used to update or depending on whether the sweep is from left to right or from right to left.

To adapt the above algorithm to , one notices that the only difference between , and , is that for the latter there is a much larger variety of possible values for the left-most index ( or ) of the tensors at the left boundary. Conventionally, (auxiliary index on the ) can only take a single value and (auxiliary index on the ) has to be , therefore one can straightforwardly construct as

where 1 means a trivial dense tensor with a single element and where we have used the label to remind the reader that this concerns symmetric Matrix Product States. For and , instead, must be able to include all possible quantum numbers that could appear after the has been applied on the . Thus a way to construct is to look at the boundary indexes and and find all the relevant combinations , namely

With such a construction of , the variational ground state search algorithm would be able to automatically detect and preserve the symmetry in the system, whether it is or, for example, .

### ii.5 Time evolution algorithm

Time evolution with can be implemented in various ways. In the following we describe two of them.

The first approach is based on multiplication between MPOs and MPSs, as in time evolution with Tchebychev polynomials Garcia-Ripoll2006 (); HalimehMcCulloch2015 (); ZaletelPollmann2014 (). These type of algorithms only depend on MPO and MPS arithmetic, namely MPO multiplication with MPS and the addition and subtraction of MPSs. Hence these approaches can be readily implemented with . For a generic construction of MPO from arbitrary Hamiltonian (including symmetric ones), one can refer to ClaudiusSchollwock2017 ().

While the approach just described can be readily implemented, it is in general slower than Matrix Product States algorithms which use Suzuki-Trotter decomposition Trotter1959 (); Suzuki1976 (). We thus describe here one way to use Suzuki-Trotter decomposition for . In these algorithms, one splits the time evolution operator into many local operations which only affect the MPS locally. A typical example is shown in the central portion of Fig.2, with orange background, in which the evolution for a certain time is divided in evolution on even and odd bonds (respectively with the Hamiltonians and ), and which can be parallelized (for more details on this MPS algorithm see, for instance, the review Schollwock2011 ()). In the following, for clarity, we focus on a unitary evolution with a Hamiltonian , however in the example discussed in Sec.III.2 we will show how this approach can also be readily applied to dissipative evolutions. For , special attentions should be paid to the terms of the Hamiltonian which change the total quantum number. Considering, for instance, that there is a local term acting on the -th site of a initially in a single number sector , the resulting state will be a superposition of three number sectors . This means that even if only acts locally on site , the left boundary tensor at site of the has to be updated and thus the effect on is non-local. As a result, we should not simply absorb non-number conserving terms into a local operator and perform the usual Suzuki-Trotter based algorithms. One way to overcome this complication is to separate an Hamiltonian into a symmetric, number-conserving part, , and an asymmetric, non-number-conserving part, :

(7) |

It is now possible to use a second order time evolution operator as (setting )

(8) |

as shown in Fig.2. For the symmetric portion, , one can perform a Suzuki-Trotter based algorithm, while for , one can treat it as an MPO and perform an MPO based time evolution. This hybrid time evolution algorithm would be efficient if the bond dimension of is small which, for instance, is the case if the non-number-conserving term is local. To give a more concrete idea before we discuss examples in the following section, local non-number conserving components of an Hamiltonian are typical in coupled resonator arrays, see for instance the review NohAngelakis2016 ().

## Iii Examples

In the following we present two exemplary applications of dynamically symmetric MPS. In Sec.III.1 we discuss a unitary evolution while in Sec.III.2 we present a boundary driven dissipative systems.

### iii.1 Ground state and time evolution of an XYZ chain

The Hamiltonian of an spin chain of size can be written as

(9) |

where , and are the operators corresponding to Pauli matrices, is the tunneling, is the interaction strength and denotes the anisotropy ranging from to , is the strength of magnetization. For the special case , this model reduces to an chain with longitudinal field which conserves the total spin and has symmetry, while for , the Hamiltonian commutes with the parity operator defined as

(10) |

thus having a discrete symmetry. In the following we work in units for which . We first apply the ground state search algorithm to with . In the following we keep the interaction strength and magnetization strength , and we compute the ground state for different values of , starting the search of the ground state from a state with zero total magnetization . For the case , the resulting ground state will be the ground state in the same number sector as the trial starting state, namely in the number sector with total spin. However, for , the resulting ground state will automatically be in the even number sector, which is a superposition of states from different number sectors that all have an even number for the total spin. If one uses a symmetric MPS for the case, one will obtain the same ground state but only in the number sector labelled by (corresponding to , i.e. total spin modulo is ). A major difference between using a symmetric MPS and our dynamically symmetric MPS is that with our method it is straightforward to extract the information about the occupation of states with different total spins from the ground state. In fact, one can simply split a dynamically symmetric MPS into many symmetric MPSs with a fixed total spin as in Eq.(5), by just splitting the tensor on the left boundary. Therefore, with dynamically symmetric MPS, one can easily compute the distribution of the ground state in different number sectors , which is shown in Fig.3(a). One can see clearly for , the ground state has a single total spin (a kronecker delta at represented by a dashed line), while the distribution becomes broader as increases to .

We then do a quench, starting from the ground state corresponding to , , and then evolving it with the Hamiltonian (9) for , as if a sudden quench occurred. This evolution can be represented as

(11) |

We chose to perform the evolution using a fourth order Runge-Kutta method based on MPOMPS multiplication. To ensure convergence of our results, we have done simulations with bond dimensions and checked various observables, for instance, for the average spin occupation we get a difference of the order of . With the approach it is straightforward to see that the quantum state of the system gradually expands into different number sectors with an even number of total spin , as shown in Fig.3(b). We note that despite the initial state is computed for an Hamiltonian with a different symmetry from that of the evolution Hamiltonian, we use the same structure of Matrix Product States both for searching the ground state and the evolution, making the passage between the two different symmetry regimes seamless.

### iii.2 Time evolution of a dissipatively boundary driven Bose-Hubbard chain

The method to evolve in time the dynamically symmetric Matrix Product States can be readily extended from unitary to dissipative dynamics. We show an example in the following. The dynamics of the system is described by a master equation of Gorini-Kossakowski-Sudarshan-Lindblad form GoriniSudarshan1976 (); Lindblad1976 ()

(12) |

where we refer to as the Lindbladian, as dissipator, and where the Hamiltonian can be written as

(13) | ||||

The dissipator can be written as

(14) |

In Eqs.(13,14) we have used the notations for the bosonic creation and annihilation operators and . is the tunneling amplitude and is the local interaction strength. denotes the strength of the dissipation acting at the edges or . If , the dissipator will impose a local thermal distribution at the edges with average occupation respectively. Similarly as for the case of the chain, in the following we work in units for which .

In general, to simulate the time evolution of this system with Matrix Product States, it is first useful to stress that, while at a given site it would be possible to have a very large (even infinite) number of bosonic particles, the first step is to limit the local Hilbert space size to a finite number, which we refer to as . For unitary systems, for an interaction strength , it is usually sufficient to keep . However, to faithfully represent a local bosonic thermal state, has to be much larger due to the long tail of the thermal distribution. In our simulations we have used and to ensure the accuracy of the results.

A standard way to deal with density matrices with Matrix Product States is to consider them as vectors and, as a result, the local Hilbert space will be of size . If one uses non-symmetric Matrix Product States for the time evolution, for instance splitting the operator into many two-sites operators, then each two-body operator would be a matrix. For , one such operator will contain (possibly complex) numbers, thus consuming a large amount memory. Furthermore, the most numeric expensive part of the Matrix Product States algorithm would be a two-site singular value decomposition performed on a matrix. For a reasonable bond dimension and local Hilbert space , this tensor would have a size of , for which it is prohibitive to run a singular value decomposition on most personal computers. On this point we remark that in GuoPoletti2015 () we kept only the diagonal and a few off-diagonal elements of the single-site reduce density matrix thus reducing the local space from to with . This was possible for the dissipator used (same as the one in Eq.(14), and the system analyzed), as the local density matrix is mainly diagonal. Such approximation is however not sufficient to study larger systems. Using it is possible, as explained in the following, to study larger systems.

We note that the unitary part of is number conserving, while the dissipator in Eq.(14) is local, thus the MPO for the non-number-conserving part has a bond dimension equal to , which is the smallest possible size. In addition, the dissipator in Eq.(14) has an ulterior property: if we write the density operator using the basis where is the total number of bosons in the ket ( in the bra), and is the vector detailing how the bosons are distributed over the sites (respectively describes how the bosons are distributed between sites), then it appears clear that the dissipator couples the element corresponding to with . This implies that if the initial condition belongs to a single number sector (e.g. at initial time ), then the density operator will have only non-zero terms for blocks in which . This block-diagonal structure of the density operator has been discussed in GuoPoletti2017b () and exploited in the context of studying transport with exact diagonalization tools. This symmetry of the density operator will automatically be preserved with and , making our hybrid time evolution algorithm introduced in Sec.II.5 an ideal tool to study this type of problems.

In our exemplary simulation, we first prepared the initial state of the system to be the ground state of a Bose-Hubbard chain of size with average filling and , which we denote as . We then turn on the dissipation, with , and . For the simulations we have truncated the local Hilbert space to and used a bond dimension for the MPS, and the time step . To ensure that our choices of and are appropriate, we have also done simulations and computed the local density with , getting a difference of the order , and with , getting a difference of the order . In Fig.4(a), we plot the distribution of the quantum state in different number sectors, i.e.

(15) |

something readily done with dynamically symmetric MPS. The initial state has a total of bosons, then during the evolution, the state becomes a superposition of states from many number sectors due to the boundary dissipative driving. The distribution of the density operator in different number sectors becomes broader with time. In Fig.4(b), we plot the average occupation as a function of site at different times. Initially, since the evolution starts from the ground state, the distribution of is symmetric around the middle site. At later times the distribution becomes unbalanced due to the different drivings at the two edges.

## Iv conclusion

In this work we have proposed a Matrix Product States algorithm that can dynamically take into account of a global symmetry or of one of its subgroups. At the same time, this method can also deal with non-number conserving systems. We have shown how this method, based on dynamically symmetric Matrix Product States, can be applied to search for the ground state of a systems and also for time evolution.

We shall note that the conventional symmetric Matrix Product States method can be viewed as a special case of the one presented here. Moreover, for systems without symmetry, or, for instance, with only symmetry, the dynamically symmetric Matrix Product States method allows to readily acquire additional information about the distribution of the state in different number sectors.

This method could be very useful in some applications in which the presence of symmetries, and their type, changes within the evolution. We have studied two such examples both for unitary and dissipative systems. In both cases the use of dynamically symmetric Matrix Product States allows to readily follow the evolution of the system and benefit of the symmetries when present. The efficiency of the method depends on the system studied. Dynamically symmetric Matrix Product States could also be extended to systems with non-Abelian symmetries, however this is beyond the scope of the current work.

###### Acknowledgements.

C. G. acknowledges support from National Natural Science Foundation of China under Grants No. 11504430 and No. 11805279. D.P. acknowledges support from the Singapore Ministry of Education, Singapore Academic Research Fund Tier-II (project MOE2016-T2-1-065).## References

- (1) M. Fannes, B. Nachtergaele, and R. Werner, Commun. Math. Phys., 144, 443 (1992).
- (2) S. Ostlund and S. Rommer, Phys. Rev. Lett., 75, 3537 (1995).
- (3) D. Perez-Garcia, and F. Verstraete, M.M. Wolf, and J.I. Cirac, Quantum Inf. Comput., 7, 401 (2007).
- (4) U. Schollwoc̈k, Ann. Phys. 326, 96 (2011).
- (5) M. B. Hastings, J. Stat. Mech: Theory Exp. (2007) P08024.
- (6) G. Vidal, Phys. Rev. Lett., 101, 110501 (2008).
- (7) G. Vidal, Phys. Rev. Lett. 91, 147902 (2003).
- (8) G. Vidal, Phys. Rev. Lett. 93, 040502 (2004).
- (9) A. J. Daley, C. Kollath, U. Schollwöck, and G. Vidal, J. Stat. Mech. Theor. Exp., P04005 (2004).
- (10) S. R. White and A. E. Feiguin, Phys. Rev. Lett. 93, 076401 (2004).
- (11) F. Verstraete, J. J. García-Ripoll, and J. I. Cirac, Phys. Rev. Lett., 93, 207204 (2004).
- (12) A. J. Daley, Adv. Phys. 63, 77 (2014)
- (13) I. de Vega, and M.C. Banuls, Phys. Rev. A (2015)
- (14) C. Guo, I. de Vega, U. Schollwöck, and D. Poletti, Phys. Rev. A 97, 053610 (2018).
- (15) X. Xu, J. Thingna, C. Guo, and D. Poletti, Phys. Rev. A 99, 012106 (2019).
- (16) T. Prosen, and M. Znidaric, J. Stat. Mech. (2009) P02035.
- (17) M. Palmero, X. Xu, C. Guo, and D. Poletti, arXiv:1901.05145 (2019).
- (18) T. Prosen, and I. Pizorn, Phys. Rev. Lett. 101, 105701 (2008).
- (19) A.J. Daley, C. Kollath, U. Schollwoeck, G. Vidal, J. Stat. Mech.: Theor. Exp. (2004) P04005.
- (20) D. Perez-Garcia, M.M. Wolf, M. Sanz, F. Verstraete, and J.I. Cirac, Phys. Rev. Lett., 100, 167202 (2008).
- (21) S. Singh, R.N.C. Pfeifer, and G. Vidal, Phys. Rev. B 83, 115125 (2011).
- (22) I. P. McCulloch, J. Stat. Mech. (2007) P10014.
- (23) J. J. García-Ripoll, New J. Phys. 8, 305 (2006).
- (24) J.C. Halimeh, F. Kolley, and I.P. McCulloch, Phys. Rev. B 92, 115130 (2015).
- (25) M.P. Zaletel, R.S.K. Mong, C. Karrasch, J.E. Moore, and F. Pollmann, Phys. Rev. B 91, 165112 (2014).
- (26) C. Hubig, I. P. McCulloch, and U. Schollwoc̈k, Phys. Rev. B 95, 035129 (2017).
- (27) H.F. Trotter, Proc. Amer. Math. Soc. 10, 545 (1959).
- (28) M. Suzuki, Comm. Math. Phys. 51, 2 (1976).
- (29) C. Noh, and D. Angelakis, Rep. Prog. Phys. 80, 016401 (2016).
- (30) V. Gorini, A. Kossakowski, and E. C. G. Sudarshan, J. Math. Phys. 17, 821 (1976).
- (31) G. Lindblad, Commun. Math. Phys. 48, 119 (1976).
- (32) C. Guo, M. Mukherjee, and D. Poletti, Phys. Rev. A 92, 023637 (2015).
- (33) C. Guo, and D. Poletti, Phys. Rev. B 96, 165409 (2017).

## Appendix A Construction of d-MPS and d-MPO for XYZ chain

In this appendix we demonstrate the process of constructing of a and using the concrete example of a -site spin XYZ chain with the magnetization strength . We label the state with quantum number , and the state with . Assuming we have a initial state which is

(16) |

The correspond to state is

In fact on site there is only one possible quantum number flowing in , only one quantum number present on the site and hence only one possible output . Similarly, at site there are two possible values for which results in . At site , the local quantum number is only , but since takes two possible values, then .

In the non-symmetric case, the MPO corresponding to for can be straightforwardly written as

(17) |

Now we first rewrite in the following form

(18) |

which is a summation of terms (because there are two bonds). To write down the corresponding , a convenient way is to first write down the product MPO for each term and them use the MPO addition rule to sum them up ClaudiusSchollwock2017 (). The product MPO corresponding to the terms , , , , can be written as

We notice that for the terms , , , the auxiliary index on the left boundary , while we have for the term , and for the term . We can also construct the remaining terms of , , , , similarly. We note that if we encode the MPO with symmetry instead of symmetry, then for the terms and will also be since and are equivalent to modulo . The final is a summation of all these terms

(19) |