Matrix Product States with adaptive global symmetries

TL;DR

This paper introduces an adaptive Matrix Product States algorithm that dynamically incorporates or analyzes symmetries like U(1) and Z2, enhancing quantum many-body simulations without prior symmetry knowledge.

Contribution

The work presents a novel MPS algorithm with adaptive global symmetries, allowing automatic symmetry detection and handling during simulations, improving flexibility and insight in quantum systems.

Findings

Efficiently captures symmetries in quantum systems.

Provides insights into symmetry-breaking dynamics.

Applicable to diverse models including XYZ and bosonic systems.

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 an adaptive $U(1)$ symmetry. This algorithm can take into account of, or benefit from, $U(1)$ or $Z_2$ 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.

Equations

$$\displaystyle\mathcal{O}=\{\mathbf{O}^{(\vec{n}_{1},I_{n_{1}}),(\vec{n}_{2},I_{n_{2}})}_{({\mathchoice{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\displaystyle\vec{}\mkern 4.0mu$}\cr\kern-3.01389pt\cr$\displaystyle n$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\textstyle\vec{}\mkern 4.0mu$}\cr\kern-3.01389pt\cr$\textstyle n$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\scriptstyle\vec{}\mkern 4.0mu$}\cr\kern-2.10971pt\cr$\scriptstyle n$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\scriptscriptstyle\vec{}\mkern 4.0mu$}\cr\kern-1.50694pt\cr$\scriptscriptstyle n$\cr}}}}_{3},I_{n_{3}}),({\mathchoice{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\displaystyle\vec{}\mkern 4.0mu$}\cr\kern-3.01389pt\cr$\displaystyle n$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\textstyle\vec{}\mkern 4.0mu$}\cr\kern-3.01389pt\cr$\textstyle n$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\scriptstyle\vec{}\mkern 4.0mu$}\cr\kern-2.10971pt\cr$\scriptstyle n$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\scriptscriptstyle\vec{}\mkern 4.0mu$}\cr\kern-1.50694pt\cr$\scriptscriptstyle n$\cr}}}}_{4},I_{n_{4}})}\}.$$

$$n_1+n_2-n_3-n_4=0,$$

$$\displaystyle\mathcal{M}_{l}=\{\mathbf{M}^{{\mathchoice{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\displaystyle\vec{}\mkern 4.0mu$}\cr\kern-3.01389pt\cr$\displaystyle\sigma$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\textstyle\vec{}\mkern 4.0mu$}\cr\kern-3.01389pt\cr$\textstyle\sigma$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\scriptstyle\vec{}\mkern 4.0mu$}\cr\kern-2.10971pt\cr$\scriptstyle\sigma$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\scriptscriptstyle\vec{}\mkern 4.0mu$}\cr\kern-1.50694pt\cr$\scriptscriptstyle\sigma$\cr}}}}_{l}}_{\vec{a}_{l},{\mathchoice{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\displaystyle\vec{}\mkern 4.0mu$}\cr\kern-3.01389pt\cr$\displaystyle a$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\textstyle\vec{}\mkern 4.0mu$}\cr\kern-3.01389pt\cr$\textstyle a$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\scriptstyle\vec{}\mkern 4.0mu$}\cr\kern-2.10971pt\cr$\scriptstyle a$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\scriptscriptstyle\vec{}\mkern 4.0mu$}\cr\kern-1.50694pt\cr$\scriptscriptstyle a$\cr}}}}_{l+1}}\},$$

$${\sigma }_l+a_{l+1}-a_l=0.$$

$$\mathrm{as}\text{-}\mathrm{MPS}=\sum _N\mathrm{s}\text{-}{\mathrm{MPS}}_N.$$

$$\displaystyle\mathcal{W}_{l}=\{\mathbf{W}^{\vec{\sigma}_{l},{\mathchoice{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\displaystyle\vec{}\mkern 4.0mu$}\cr\kern-3.01389pt\cr$\displaystyle\tau$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\textstyle\vec{}\mkern 4.0mu$}\cr\kern-3.01389pt\cr$\textstyle\tau$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\scriptstyle\vec{}\mkern 4.0mu$}\cr\kern-2.10971pt\cr$\scriptstyle\tau$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\scriptscriptstyle\vec{}\mkern 4.0mu$}\cr\kern-1.50694pt\cr$\scriptscriptstyle\tau$\cr}}}}_{l}}_{\vec{b}_{l},{\mathchoice{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\displaystyle\vec{}\mkern 4.0mu$}\cr\kern-3.01389pt\cr$\displaystyle b$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\textstyle\vec{}\mkern 4.0mu$}\cr\kern-3.01389pt\cr$\textstyle b$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\scriptstyle\vec{}\mkern 4.0mu$}\cr\kern-2.10971pt\cr$\scriptstyle b$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\scriptscriptstyle\vec{}\mkern 4.0mu$}\cr\kern-1.50694pt\cr$\scriptscriptstyle b$\cr}}}}_{l+1}}\},$$

$${\tau }_l+b_{l+1}-{\sigma }_l-b_l=0.$$

$$\mathrm{as}\text{-}\mathrm{MPO}=\sum _{N,N^{\prime }}\mathrm{s}\text{-}{\mathrm{MPO}}_N^{N^{\prime }},$$

$$\sum _{a_1,a_2,\dots ,a_{L+1}}{\mathcal{A}}_1\dots {\mathcal{A}}_{l-1}{\mathcal{M}}_l{\mathcal{M}}_{l+1}{\mathcal{B}}_{l+2}\dots {\mathcal{B}}_L,$$

$$\displaystyle\sum_{\sigma_{l},a_{l}}\mathcal{A}^{{\mathchoice{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\displaystyle\vec{}\mkern 4.0mu$}\cr\kern-3.01389pt\cr$\displaystyle\sigma$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\textstyle\vec{}\mkern 4.0mu$}\cr\kern-3.01389pt\cr$\textstyle\sigma$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\scriptstyle\vec{}\mkern 4.0mu$}\cr\kern-2.10971pt\cr$\scriptstyle\sigma$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\scriptscriptstyle\vec{}\mkern 4.0mu$}\cr\kern-1.50694pt\cr$\scriptscriptstyle\sigma$\cr}}}}_{l}}_{\vec{a}_{l},{\mathchoice{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\displaystyle\vec{}\mkern 4.0mu$}\cr\kern-3.01389pt\cr$\displaystyle a$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\textstyle\vec{}\mkern 4.0mu$}\cr\kern-3.01389pt\cr$\textstyle a$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\scriptstyle\vec{}\mkern 4.0mu$}\cr\kern-2.10971pt\cr$\scriptstyle a$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\scriptscriptstyle\vec{}\mkern 4.0mu$}\cr\kern-1.50694pt\cr$\scriptscriptstyle a$\cr}}}}_{l+1}}{\rm conjugate}(\mathcal{A})^{\vec{\sigma}_{l}}_{{\mathchoice{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\displaystyle\vec{}\mkern 4.0mu$}\cr\kern-3.01389pt\cr$\displaystyle a$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\textstyle\vec{}\mkern 4.0mu$}\cr\kern-3.01389pt\cr$\textstyle a$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\scriptstyle\vec{}\mkern 4.0mu$}\cr\kern-2.10971pt\cr$\scriptstyle a$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\scriptscriptstyle\vec{}\mkern 4.0mu$}\cr\kern-1.50694pt\cr$\scriptscriptstyle a$\cr}}}}_{l},\vec{a}^{\prime}_{l+1}}=\mathcal{I}_{{\mathchoice{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\displaystyle\vec{}\mkern 4.0mu$}\cr\kern-3.01389pt\cr$\displaystyle a$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\textstyle\vec{}\mkern 4.0mu$}\cr\kern-3.01389pt\cr$\textstyle a$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\scriptstyle\vec{}\mkern 4.0mu$}\cr\kern-2.10971pt\cr$\scriptstyle a$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\scriptscriptstyle\vec{}\mkern 4.0mu$}\cr\kern-1.50694pt\cr$\scriptscriptstyle a$\cr}}}}_{l+1},\vec{a}_{l+1}^{\prime}},$$

$$\displaystyle\sum_{\sigma_{l},a_{l+1}}\mathcal{B}^{{\mathchoice{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\displaystyle\vec{}\mkern 4.0mu$}\cr\kern-3.01389pt\cr$\displaystyle\sigma$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\textstyle\vec{}\mkern 4.0mu$}\cr\kern-3.01389pt\cr$\textstyle\sigma$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\scriptstyle\vec{}\mkern 4.0mu$}\cr\kern-2.10971pt\cr$\scriptstyle\sigma$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\scriptscriptstyle\vec{}\mkern 4.0mu$}\cr\kern-1.50694pt\cr$\scriptscriptstyle\sigma$\cr}}}}_{l}}_{\vec{a}_{l},{\mathchoice{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\displaystyle\vec{}\mkern 4.0mu$}\cr\kern-3.01389pt\cr$\displaystyle a$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\textstyle\vec{}\mkern 4.0mu$}\cr\kern-3.01389pt\cr$\textstyle a$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\scriptstyle\vec{}\mkern 4.0mu$}\cr\kern-2.10971pt\cr$\scriptstyle a$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\scriptscriptstyle\vec{}\mkern 4.0mu$}\cr\kern-1.50694pt\cr$\scriptscriptstyle a$\cr}}}}_{l+1}}{\rm conjugate}(\mathcal{B})^{\vec{\sigma}_{l}}_{{\mathchoice{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\displaystyle\vec{}\mkern 4.0mu$}\cr\kern-3.01389pt\cr$\displaystyle a$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\textstyle\vec{}\mkern 4.0mu$}\cr\kern-3.01389pt\cr$\textstyle a$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\scriptstyle\vec{}\mkern 4.0mu$}\cr\kern-2.10971pt\cr$\scriptstyle a$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\scriptscriptstyle\vec{}\mkern 4.0mu$}\cr\kern-1.50694pt\cr$\scriptscriptstyle a$\cr}}}}_{l}^{\prime},\vec{a}_{l+1}}=\mathcal{I}_{\vec{a}_{l},{\mathchoice{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\displaystyle\vec{}\mkern 4.0mu$}\cr\kern-3.01389pt\cr$\displaystyle a$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\textstyle\vec{}\mkern 4.0mu$}\cr\kern-3.01389pt\cr$\textstyle a$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\scriptstyle\vec{}\mkern 4.0mu$}\cr\kern-2.10971pt\cr$\scriptstyle a$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\scriptscriptstyle\vec{}\mkern 4.0mu$}\cr\kern-1.50694pt\cr$\scriptscriptstyle a$\cr}}}}_{l}^{\prime}}.$$

$$\displaystyle\mathcal{L}_{\vec{a}_{l}^{\prime},{\mathchoice{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\displaystyle\vec{}\mkern 4.0mu$}\cr\kern-3.01389pt\cr$\displaystyle b$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\textstyle\vec{}\mkern 4.0mu$}\cr\kern-3.01389pt\cr$\textstyle b$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\scriptstyle\vec{}\mkern 4.0mu$}\cr\kern-2.10971pt\cr$\scriptstyle b$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\scriptscriptstyle\vec{}\mkern 4.0mu$}\cr\kern-1.50694pt\cr$\scriptscriptstyle b$\cr}}}}_{l},{\mathchoice{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\displaystyle\vec{}\mkern 4.0mu$}\cr\kern-3.01389pt\cr$\displaystyle a$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\textstyle\vec{}\mkern 4.0mu$}\cr\kern-3.01389pt\cr$\textstyle a$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\scriptstyle\vec{}\mkern 4.0mu$}\cr\kern-2.10971pt\cr$\scriptstyle a$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\scriptscriptstyle\vec{}\mkern 4.0mu$}\cr\kern-1.50694pt\cr$\scriptscriptstyle a$\cr}}}}_{l}}=\sum_{\begin{subarray}{c}{a_{l-1}^{\prime},b_{l-1},a_{l-1}}\\ {\sigma_{l-1},\sigma_{l-1}^{\prime}}\end{subarray}}\mathcal{L}_{\vec{a}_{l-1}^{\prime},{\mathchoice{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\displaystyle\vec{}\mkern 4.0mu$}\cr\kern-3.01389pt\cr$\displaystyle b$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\textstyle\vec{}\mkern 4.0mu$}\cr\kern-3.01389pt\cr$\textstyle b$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\scriptstyle\vec{}\mkern 4.0mu$}\cr\kern-2.10971pt\cr$\scriptstyle b$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\scriptscriptstyle\vec{}\mkern 4.0mu$}\cr\kern-1.50694pt\cr$\scriptscriptstyle b$\cr}}}}_{l-1},{\mathchoice{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\displaystyle\vec{}\mkern 4.0mu$}\cr\kern-3.01389pt\cr$\displaystyle a$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\textstyle\vec{}\mkern 4.0mu$}\cr\kern-3.01389pt\cr$\textstyle a$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\scriptstyle\vec{}\mkern 4.0mu$}\cr\kern-2.10971pt\cr$\scriptstyle a$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\scriptscriptstyle\vec{}\mkern 4.0mu$}\cr\kern-1.50694pt\cr$\scriptscriptstyle a$\cr}}}}_{l-1}}\;\;\times$$

$$\displaystyle\mathcal{A}_{\vec{a}_{l-1},{\mathchoice{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\displaystyle\vec{}\mkern 4.0mu$}\cr\kern-3.01389pt\cr$\displaystyle a$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\textstyle\vec{}\mkern 4.0mu$}\cr\kern-3.01389pt\cr$\textstyle a$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\scriptstyle\vec{}\mkern 4.0mu$}\cr\kern-2.10971pt\cr$\scriptstyle a$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\scriptscriptstyle\vec{}\mkern 4.0mu$}\cr\kern-1.50694pt\cr$\scriptscriptstyle a$\cr}}}}_{l}}^{{\mathchoice{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\displaystyle\vec{}\mkern 4.0mu$}\cr\kern-3.01389pt\cr$\displaystyle\sigma$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\textstyle\vec{}\mkern 4.0mu$}\cr\kern-3.01389pt\cr$\textstyle\sigma$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\scriptstyle\vec{}\mkern 4.0mu$}\cr\kern-2.10971pt\cr$\scriptstyle\sigma$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\scriptscriptstyle\vec{}\mkern 4.0mu$}\cr\kern-1.50694pt\cr$\scriptscriptstyle\sigma$\cr}}}}_{l-1}}\;\;{\rm conjugate}(\mathcal{A})_{{\mathchoice{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\displaystyle\vec{}\mkern 4.0mu$}\cr\kern-3.01389pt\cr$\displaystyle a$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\textstyle\vec{}\mkern 4.0mu$}\cr\kern-3.01389pt\cr$\textstyle a$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\scriptstyle\vec{}\mkern 4.0mu$}\cr\kern-2.10971pt\cr$\scriptstyle a$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\scriptscriptstyle\vec{}\mkern 4.0mu$}\cr\kern-1.50694pt\cr$\scriptscriptstyle a$\cr}}}}_{l-1}^{\prime},\vec{a}_{l}^{\prime}}^{\vec{\sigma}_{l-1}^{\prime}}\mathcal{W}_{\vec{b}_{l-1},{\mathchoice{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\displaystyle\vec{}\mkern 4.0mu$}\cr\kern-3.01389pt\cr$\displaystyle b$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\textstyle\vec{}\mkern 4.0mu$}\cr\kern-3.01389pt\cr$\textstyle b$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\scriptstyle\vec{}\mkern 4.0mu$}\cr\kern-2.10971pt\cr$\scriptstyle b$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\scriptscriptstyle\vec{}\mkern 4.0mu$}\cr\kern-1.50694pt\cr$\scriptscriptstyle b$\cr}}}}_{l}}^{\vec{\sigma}_{l-1},{\mathchoice{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\displaystyle\vec{}\mkern 4.0mu$}\cr\kern-3.01389pt\cr$\displaystyle\sigma$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\textstyle\vec{}\mkern 4.0mu$}\cr\kern-3.01389pt\cr$\textstyle\sigma$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\scriptstyle\vec{}\mkern 4.0mu$}\cr\kern-2.10971pt\cr$\scriptstyle\sigma$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\scriptscriptstyle\vec{}\mkern 4.0mu$}\cr\kern-1.50694pt\cr$\scriptscriptstyle\sigma$\cr}}}}_{l-1}^{\prime}},$$

$$\displaystyle\mathcal{R}_{{\mathchoice{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\displaystyle\vec{}\mkern 4.0mu$}\cr\kern-3.01389pt\cr$\displaystyle a$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\textstyle\vec{}\mkern 4.0mu$}\cr\kern-3.01389pt\cr$\textstyle a$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\scriptstyle\vec{}\mkern 4.0mu$}\cr\kern-2.10971pt\cr$\scriptstyle a$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\scriptscriptstyle\vec{}\mkern 4.0mu$}\cr\kern-1.50694pt\cr$\scriptscriptstyle a$\cr}}}}_{l+2}^{\prime},\vec{b}_{l+2},\vec{a}_{l+2}}=\sum_{\begin{subarray}{c}{a_{l+3}^{\prime},b_{l+3},a_{l+3}}\\ {\sigma_{l+2},\sigma_{l+2}^{\prime}}\end{subarray}}\mathcal{R}_{{\mathchoice{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\displaystyle\vec{}\mkern 4.0mu$}\cr\kern-3.01389pt\cr$\displaystyle a$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\textstyle\vec{}\mkern 4.0mu$}\cr\kern-3.01389pt\cr$\textstyle a$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\scriptstyle\vec{}\mkern 4.0mu$}\cr\kern-2.10971pt\cr$\scriptstyle a$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\scriptscriptstyle\vec{}\mkern 4.0mu$}\cr\kern-1.50694pt\cr$\scriptscriptstyle a$\cr}}}}_{l+3}^{\prime},\vec{b}_{l+3},\vec{a}_{l+3}}\;\;\times$$

$$\displaystyle\mathcal{B}_{\vec{a}_{l+2},{\mathchoice{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\displaystyle\vec{}\mkern 4.0mu$}\cr\kern-3.01389pt\cr$\displaystyle a$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\textstyle\vec{}\mkern 4.0mu$}\cr\kern-3.01389pt\cr$\textstyle a$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\scriptstyle\vec{}\mkern 4.0mu$}\cr\kern-2.10971pt\cr$\scriptstyle a$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\scriptscriptstyle\vec{}\mkern 4.0mu$}\cr\kern-1.50694pt\cr$\scriptscriptstyle a$\cr}}}}_{l+3}}^{{\mathchoice{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\displaystyle\vec{}\mkern 4.0mu$}\cr\kern-3.01389pt\cr$\displaystyle\sigma$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\textstyle\vec{}\mkern 4.0mu$}\cr\kern-3.01389pt\cr$\textstyle\sigma$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\scriptstyle\vec{}\mkern 4.0mu$}\cr\kern-2.10971pt\cr$\scriptstyle\sigma$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\scriptscriptstyle\vec{}\mkern 4.0mu$}\cr\kern-1.50694pt\cr$\scriptscriptstyle\sigma$\cr}}}}_{l+2}}{\rm conjugate}(\mathcal{B})_{{\mathchoice{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\displaystyle\vec{}\mkern 4.0mu$}\cr\kern-3.01389pt\cr$\displaystyle a$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\textstyle\vec{}\mkern 4.0mu$}\cr\kern-3.01389pt\cr$\textstyle a$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\scriptstyle\vec{}\mkern 4.0mu$}\cr\kern-2.10971pt\cr$\scriptstyle a$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\scriptscriptstyle\vec{}\mkern 4.0mu$}\cr\kern-1.50694pt\cr$\scriptscriptstyle a$\cr}}}}_{l+2}^{\prime},\vec{a}_{l+3}^{\prime}}^{\vec{\sigma}_{l+2}^{\prime}}\mathcal{W}_{\vec{b}_{l+2},{\mathchoice{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\displaystyle\vec{}\mkern 4.0mu$}\cr\kern-3.01389pt\cr$\displaystyle b$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\textstyle\vec{}\mkern 4.0mu$}\cr\kern-3.01389pt\cr$\textstyle b$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\scriptstyle\vec{}\mkern 4.0mu$}\cr\kern-2.10971pt\cr$\scriptstyle b$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\scriptscriptstyle\vec{}\mkern 4.0mu$}\cr\kern-1.50694pt\cr$\scriptscriptstyle b$\cr}}}}_{l+3}}^{\vec{\sigma}_{l+2},{\mathchoice{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\displaystyle\vec{}\mkern 4.0mu$}\cr\kern-3.01389pt\cr$\displaystyle\sigma$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\textstyle\vec{}\mkern 4.0mu$}\cr\kern-3.01389pt\cr$\textstyle\sigma$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\scriptstyle\vec{}\mkern 4.0mu$}\cr\kern-2.10971pt\cr$\scriptstyle\sigma$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\scriptscriptstyle\vec{}\mkern 4.0mu$}\cr\kern-1.50694pt\cr$\scriptscriptstyle\sigma$\cr}}}}_{l+2}^{\prime}}.$$

$$\displaystyle\hat{\rm H}_{{\rm eff}\;{\mathchoice{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\displaystyle\vec{}\mkern 4.0mu$}\cr\kern-3.01389pt\cr$\displaystyle a$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\textstyle\vec{}\mkern 4.0mu$}\cr\kern-3.01389pt\cr$\textstyle a$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\scriptstyle\vec{}\mkern 4.0mu$}\cr\kern-2.10971pt\cr$\scriptstyle a$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\scriptscriptstyle\vec{}\mkern 4.0mu$}\cr\kern-1.50694pt\cr$\scriptscriptstyle a$\cr}}}}_{l},\vec{\sigma}_{l},\vec{\sigma}_{l+1},\vec{a}_{l+2}}^{\;\;\;\;\;\vec{a}_{l}^{\prime},{\mathchoice{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\displaystyle\vec{}\mkern 4.0mu$}\cr\kern-3.01389pt\cr$\displaystyle\sigma$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\textstyle\vec{}\mkern 4.0mu$}\cr\kern-3.01389pt\cr$\textstyle\sigma$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\scriptstyle\vec{}\mkern 4.0mu$}\cr\kern-2.10971pt\cr$\scriptstyle\sigma$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\scriptscriptstyle\vec{}\mkern 4.0mu$}\cr\kern-1.50694pt\cr$\scriptscriptstyle\sigma$\cr}}}}_{l}^{\prime},{\mathchoice{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\displaystyle\vec{}\mkern 4.0mu$}\cr\kern-3.01389pt\cr$\displaystyle\sigma$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\textstyle\vec{}\mkern 4.0mu$}\cr\kern-3.01389pt\cr$\textstyle\sigma$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\scriptstyle\vec{}\mkern 4.0mu$}\cr\kern-2.10971pt\cr$\scriptstyle\sigma$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\scriptscriptstyle\vec{}\mkern 4.0mu$}\cr\kern-1.50694pt\cr$\scriptscriptstyle\sigma$\cr}}}}_{l+1}^{\prime},{\mathchoice{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\displaystyle\vec{}\mkern 4.0mu$}\cr\kern-3.01389pt\cr$\displaystyle a$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\textstyle\vec{}\mkern 4.0mu$}\cr\kern-3.01389pt\cr$\textstyle a$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\scriptstyle\vec{}\mkern 4.0mu$}\cr\kern-2.10971pt\cr$\scriptstyle a$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\scriptscriptstyle\vec{}\mkern 4.0mu$}\cr\kern-1.50694pt\cr$\scriptscriptstyle a$\cr}}}}_{l+2}^{\prime}}=$$

$$\displaystyle\mathcal{W}_{\vec{b}_{l},{\mathchoice{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\displaystyle\vec{}\mkern 4.0mu$}\cr\kern-3.01389pt\cr$\displaystyle b$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\textstyle\vec{}\mkern 4.0mu$}\cr\kern-3.01389pt\cr$\textstyle b$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\scriptstyle\vec{}\mkern 4.0mu$}\cr\kern-2.10971pt\cr$\scriptstyle b$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\scriptscriptstyle\vec{}\mkern 4.0mu$}\cr\kern-1.50694pt\cr$\scriptscriptstyle b$\cr}}}}_{l+1}}^{\vec{\sigma}_{l},{\mathchoice{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\displaystyle\vec{}\mkern 4.0mu$}\cr\kern-3.01389pt\cr$\displaystyle\sigma$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\textstyle\vec{}\mkern 4.0mu$}\cr\kern-3.01389pt\cr$\textstyle\sigma$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\scriptstyle\vec{}\mkern 4.0mu$}\cr\kern-2.10971pt\cr$\scriptstyle\sigma$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\scriptscriptstyle\vec{}\mkern 4.0mu$}\cr\kern-1.50694pt\cr$\scriptscriptstyle\sigma$\cr}}}}_{l}^{\prime}}\mathcal{W}_{\vec{b}_{l+1},{\mathchoice{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\displaystyle\vec{}\mkern 4.0mu$}\cr\kern-3.01389pt\cr$\displaystyle b$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\textstyle\vec{}\mkern 4.0mu$}\cr\kern-3.01389pt\cr$\textstyle b$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\scriptstyle\vec{}\mkern 4.0mu$}\cr\kern-2.10971pt\cr$\scriptstyle b$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\scriptscriptstyle\vec{}\mkern 4.0mu$}\cr\kern-1.50694pt\cr$\scriptscriptstyle b$\cr}}}}_{l+2}}^{\vec{\sigma}_{l+1},{\mathchoice{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\displaystyle\vec{}\mkern 4.0mu$}\cr\kern-3.01389pt\cr$\displaystyle\sigma$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\textstyle\vec{}\mkern 4.0mu$}\cr\kern-3.01389pt\cr$\textstyle\sigma$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\scriptstyle\vec{}\mkern 4.0mu$}\cr\kern-2.10971pt\cr$\scriptstyle\sigma$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\scriptscriptstyle\vec{}\mkern 4.0mu$}\cr\kern-1.50694pt\cr$\scriptscriptstyle\sigma$\cr}}}}_{l+1}^{\prime}}$$

$$\displaystyle{\rm SVD}(\mathcal{M}_{\vec{a}_{l},{\mathchoice{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\displaystyle\vec{}\mkern 4.0mu$}\cr\kern-3.01389pt\cr$\displaystyle\sigma$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\textstyle\vec{}\mkern 4.0mu$}\cr\kern-3.01389pt\cr$\textstyle\sigma$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\scriptstyle\vec{}\mkern 4.0mu$}\cr\kern-2.10971pt\cr$\scriptstyle\sigma$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\scriptscriptstyle\vec{}\mkern 4.0mu$}\cr\kern-1.50694pt\cr$\scriptscriptstyle\sigma$\cr}}}}_{l},{\mathchoice{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\displaystyle\vec{}\mkern 4.0mu$}\cr\kern-3.01389pt\cr$\displaystyle\sigma$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\textstyle\vec{}\mkern 4.0mu$}\cr\kern-3.01389pt\cr$\textstyle\sigma$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\scriptstyle\vec{}\mkern 4.0mu$}\cr\kern-2.10971pt\cr$\scriptstyle\sigma$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\scriptscriptstyle\vec{}\mkern 4.0mu$}\cr\kern-1.50694pt\cr$\scriptscriptstyle\sigma$\cr}}}}_{l+1},{\mathchoice{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\displaystyle\vec{}\mkern 4.0mu$}\cr\kern-3.01389pt\cr$\displaystyle a$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\textstyle\vec{}\mkern 4.0mu$}\cr\kern-3.01389pt\cr$\textstyle a$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\scriptstyle\vec{}\mkern 4.0mu$}\cr\kern-2.10971pt\cr$\scriptstyle a$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\scriptscriptstyle\vec{}\mkern 4.0mu$}\cr\kern-1.50694pt\cr$\scriptscriptstyle a$\cr}}}}_{l+2}})=\sum_{s}\mathcal{A}_{\vec{a}_{l},{\mathchoice{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\displaystyle\vec{}\mkern 4.0mu$}\cr\kern-3.01389pt\cr$\displaystyle s$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\textstyle\vec{}\mkern 4.0mu$}\cr\kern-3.01389pt\cr$\textstyle s$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\scriptstyle\vec{}\mkern 4.0mu$}\cr\kern-2.10971pt\cr$\scriptstyle s$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\scriptscriptstyle\vec{}\mkern 4.0mu$}\cr\kern-1.50694pt\cr$\scriptscriptstyle s$\cr}}}}}^{{\mathchoice{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\displaystyle\vec{}\mkern 4.0mu$}\cr\kern-3.01389pt\cr$\displaystyle\sigma$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\textstyle\vec{}\mkern 4.0mu$}\cr\kern-3.01389pt\cr$\textstyle\sigma$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\scriptstyle\vec{}\mkern 4.0mu$}\cr\kern-2.10971pt\cr$\scriptstyle\sigma$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\scriptscriptstyle\vec{}\mkern 4.0mu$}\cr\kern-1.50694pt\cr$\scriptscriptstyle\sigma$\cr}}}}_{l}}\;\mathcal{S}_{\vec{s},{\mathchoice{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\displaystyle\vec{}\mkern 4.0mu$}\cr\kern-3.01389pt\cr$\displaystyle s$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\textstyle\vec{}\mkern 4.0mu$}\cr\kern-3.01389pt\cr$\textstyle s$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\scriptstyle\vec{}\mkern 4.0mu$}\cr\kern-2.10971pt\cr$\scriptstyle s$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\scriptscriptstyle\vec{}\mkern 4.0mu$}\cr\kern-1.50694pt\cr$\scriptscriptstyle s$\cr}}}}}\;\mathcal{B}_{\vec{s},{\mathchoice{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\displaystyle\vec{}\mkern 4.0mu$}\cr\kern-3.01389pt\cr$\displaystyle a$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\textstyle\vec{}\mkern 4.0mu$}\cr\kern-3.01389pt\cr$\textstyle a$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\scriptstyle\vec{}\mkern 4.0mu$}\cr\kern-2.10971pt\cr$\scriptstyle a$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\scriptscriptstyle\vec{}\mkern 4.0mu$}\cr\kern-1.50694pt\cr$\scriptscriptstyle a$\cr}}}}_{l+2}}^{{\mathchoice{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\displaystyle\vec{}\mkern 4.0mu$}\cr\kern-3.01389pt\cr$\displaystyle\sigma$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\textstyle\vec{}\mkern 4.0mu$}\cr\kern-3.01389pt\cr$\textstyle\sigma$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\scriptstyle\vec{}\mkern 4.0mu$}\cr\kern-2.10971pt\cr$\scriptstyle\sigma$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\scriptscriptstyle\vec{}\mkern 4.0mu$}\cr\kern-1.50694pt\cr$\scriptscriptstyle\sigma$\cr}}}}_{l+1}},$$

$$\displaystyle\mathcal{L}_{1}^{{\rm s}\text{-}{\rm MPS}}=\{\textbf{1}_{\vec{N},{\mathchoice{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\displaystyle\vec{}\mkern 4.0mu$}\cr\kern-3.01389pt\cr$\displaystyle 0$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\textstyle\vec{}\mkern 4.0mu$}\cr\kern-3.01389pt\cr$\textstyle 0$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\scriptstyle\vec{}\mkern 4.0mu$}\cr\kern-2.10971pt\cr$\scriptstyle 0$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\scriptscriptstyle\vec{}\mkern 4.0mu$}\cr\kern-1.50694pt\cr$\scriptscriptstyle 0$\cr}}}},{\mathchoice{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\displaystyle\vec{}\mkern 4.0mu$}\cr\kern-3.01389pt\cr$\displaystyle N$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\textstyle\vec{}\mkern 4.0mu$}\cr\kern-3.01389pt\cr$\textstyle N$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\scriptstyle\vec{}\mkern 4.0mu$}\cr\kern-2.10971pt\cr$\scriptstyle N$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\scriptscriptstyle\vec{}\mkern 4.0mu$}\cr\kern-1.50694pt\cr$\scriptscriptstyle N$\cr}}}}}\},$$

$$\displaystyle\mathcal{L}_{1}=\{\textbf{1}_{\vec{a}_{1}+\vec{b}_{1},{\mathchoice{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\displaystyle\vec{}\mkern 4.0mu$}\cr\kern-3.01389pt\cr$\displaystyle b$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\textstyle\vec{}\mkern 4.0mu$}\cr\kern-3.01389pt\cr$\textstyle b$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\scriptstyle\vec{}\mkern 4.0mu$}\cr\kern-2.10971pt\cr$\scriptstyle b$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\scriptscriptstyle\vec{}\mkern 4.0mu$}\cr\kern-1.50694pt\cr$\scriptscriptstyle b$\cr}}}}_{1},{\mathchoice{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\displaystyle\vec{}\mkern 4.0mu$}\cr\kern-3.01389pt\cr$\displaystyle a$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\textstyle\vec{}\mkern 4.0mu$}\cr\kern-3.01389pt\cr$\textstyle a$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\scriptstyle\vec{}\mkern 4.0mu$}\cr\kern-2.10971pt\cr$\scriptstyle a$\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{$\scriptscriptstyle\vec{}\mkern 4.0mu$}\cr\kern-1.50694pt\cr$\scriptscriptstyle a$\cr}}}}_{1}}\}.$$

$$\hat{\mathrm{H}}={\hat{\mathrm{H}}}_s+{\hat{\mathrm{H}}}_a.$$

$$\hat{\mathrm{U}}(\mathrm{dt})=\exp \left(\frac{-\mathrm{i}{\hat{\mathrm{H}}}_a\mathrm{dt}}{2}\right)\exp \left(-\mathrm{i}{\hat{\mathrm{H}}}_s\mathrm{dt}\right)\exp \left(-\mathrm{i}\frac{{\hat{\mathrm{H}}}_a\mathrm{dt}}{2}\right)$$

$${\hat{H}}_{\mathrm{XYZ}}(\gamma )=\sum _{l=1}^{L-1}$$

$$+\Delta {\hat{\sigma }}_l^z{\hat{\sigma }}_{l+1}^z]+h\sum _{l=1}^L{\hat{\sigma }}_l^z,$$

$$\mathcal{P}=\exp \left(\mathrm{i}\pi {\hat{S}}_T^z\right),$$

$$|\psi (t)⟩=\exp \left(-\mathrm{i}{\hat{H}}_{\mathrm{XYZ}}(\gamma =0.5)t\right)|G{S}_{\mathrm{XXZ}}⟩.$$

$$\frac{d\hat{\rho }}{dt}=\hat{\mathcal{L}}(\hat{\rho })=-\frac{\mathrm{i}}{\hslash }[{\hat{H}}_{\mathrm{BH}},\hat{\rho }]+\hat{\mathcal{D}}(\hat{\rho }),$$

$${\hat{H}}_{\mathrm{BH}}=$$

$$+\frac{U}{2}\sum _{l=1}^L{\hat{n}}_l({\hat{n}}_l-1).$$

$$\hat{\mathcal{D}}(\hat{\rho })=\sum _{l=1,L}$$

$$+{\bar{n}}_l\left(2{\hat{a}}_l^{†}\hat{\rho }{\hat{a}}_l-\{{\hat{a}}_l{\hat{a}}_l^{†},\hat{\rho }\}\right)].$$

$$P_N=\sum _{{\vec{n}}_N}⟨{\vec{n}}_N;N|\hat{\rho }|{\vec{n}}_N;N⟩,$$

$$|\psi ⟩=\frac{\sqrt{2}}{2}\left(|110⟩+|100⟩\right)$$

Full text

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 $U(1)$ symmetry. In other cases a system could have a discrete symmetry, like a $Z_{2}$ symmetry, which is a subgroup of $U(1)$. In such situations, it is often advantageous to explicitly encode the quantum symmetry into the tensors of the MPS DaleyVidal2004 ; Perez-GarciaVidal2008 ; SinghVidal2011 , which are also called symmetric Matrix Product States (${\rm s}\text{-}{\rm MPS}$). ${\rm s}\text{-}{\rm MPS}$ would usually result in more efficient numerical algorithms because it only explores the portion of the Hilbert space corresponding to the relevant symmetry sector. ${\rm s}\text{-}{\rm MPS}$ 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 ${\rm XXZ}$ chain, which has $U(1)$ symmetry, and then suddenly tune the parameter so that the system is modeled by a $Z_{2}-$symmetric ${\rm XYZ}$ chain. With conventional ${\rm s}\text{-}{\rm MPS}$, one would have to manually convert the ground state from $U(1)$ to $Z_{2}$ symmetry for the time evolution.

In this work we propose MPS with adaptive global $U(1)$ symmetry to overcome this difficulty, which we will refer to as adaptively symmetric Matrix Product States (${\rm as}\text{-}{\rm MPS}$). The central idea of our approach is to generalize $U(1)$ ${\rm s}\text{-}{\rm MPS}$ so that instead of a fixed total quantum number, it allows a superposition of different total quantum numbers only. Importantly, the amount of different quantum numbers required changes automatically during the calculations including only the quantum numbers required. As a result, a $U(1)$ ${\rm as}\text{-}{\rm MPS}$ is able to treat systems with $U(1)$ symmetry or systems whose symmetry group is a subgroup of $U(1)$, 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 $U(1)$ ${\rm as}\text{-}{\rm MPS}$, 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 ${\rm as}\text{-}{\rm MPS}$, one for a unitary system and the other for a dissipative system; In Sec.IV we draw our conclusions.

II method

Similarly to ${\rm s}\text{-}{\rm MPS}$, the building blocks of ${\rm as}\text{-}{\rm MPS}$ are symmetry protected tensors. The definition of a $U(1)$ 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 adaptively 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 $U(1)$ 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 $4$ dimensional symmetry protected tensor $\mathcal{O}$ could be written as

[TABLE]

Here \mathbf{O}^{(\vec{n}_{1},I_{n_{1}}),(\vec{n}_{2},I_{n_{2}})}_{({\mathchoice{\vbox{\offinterlineskip\halign{#\cr\reflectbox{\displaystyle\vec{}\mkern 4.0mu}\cr\kern-3.01389pt\cr\displaystyle n\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{\textstyle\vec{}\mkern 4.0mu}\cr\kern-3.01389pt\cr\textstyle n\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{\scriptstyle\vec{}\mkern 4.0mu}\cr\kern-2.10971pt\cr\scriptstyle n\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{\scriptscriptstyle\vec{}\mkern 4.0mu}\cr\kern-1.50694pt\cr\scriptscriptstyle n\cr}}}}_{3},I_{n_{3}}),({\mathchoice{\vbox{\offinterlineskip\halign{#\cr\reflectbox{\displaystyle\vec{}\mkern 4.0mu}\cr\kern-3.01389pt\cr\displaystyle n\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{\textstyle\vec{}\mkern 4.0mu}\cr\kern-3.01389pt\cr\textstyle n\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{\scriptstyle\vec{}\mkern 4.0mu}\cr\kern-2.10971pt\cr\scriptstyle n\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{\scriptscriptstyle\vec{}\mkern 4.0mu}\cr\kern-1.50694pt\cr\scriptscriptstyle n\cr}}}}_{4},I_{n_{4}})} is a $4$ dimensional dense tensor with dimensions $I_{n_{i}}$, and $n_{i}$ for $1\leq i\leq 4$ labels the four quantum number (note that here we don’t differentiate between superscripts and subscripts). The quantum numbers have two directions, $\vec{n}_{i}$ and {\mathchoice{\vbox{\offinterlineskip\halign{#\cr\reflectbox{\displaystyle\vec{}\mkern 4.0mu}\cr\kern-4.30554pt\cr\displaystyle n\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{\textstyle\vec{}\mkern 4.0mu}\cr\kern-4.30554pt\cr\textstyle n\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{\scriptstyle\vec{}\mkern 4.0mu}\cr\kern-3.01389pt\cr\scriptstyle n\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{\scriptscriptstyle\vec{}\mkern 4.0mu}\cr\kern-2.15277pt\cr\scriptscriptstyle n\cr}}}}_{i}, depending on whether, respectively, the quantum numbers flows in or out. The $U(1)$ symmetry can be explicitly enforced by the fusion rule

[TABLE]

where $n_{3}$ and $n_{4}$ are subtracted because their quantum numbers flow outside of the tensor $\mathcal{O}$, while $n_{1}$ and $n_{2}$ are added because their quantum numbers flow inside of the tensor $\mathcal{O}$. For the implementation of algorithms with symmetric Matrix Product States, it is convenient 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 $(\vec{n}_{i},I_{n_{i}})$ simply as $\vec{n}_{i}$ and ({\mathchoice{\vbox{\offinterlineskip\halign{#\cr\reflectbox{\displaystyle\vec{}\mkern 4.0mu}\cr\kern-4.30554pt\cr\displaystyle n\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{\textstyle\vec{}\mkern 4.0mu}\cr\kern-4.30554pt\cr\textstyle n\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{\scriptstyle\vec{}\mkern 4.0mu}\cr\kern-3.01389pt\cr\scriptstyle n\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{\scriptscriptstyle\vec{}\mkern 4.0mu}\cr\kern-2.15277pt\cr\scriptscriptstyle n\cr}}}}_{i},I_{n_{i}}) as {\mathchoice{\vbox{\offinterlineskip\halign{#\cr\reflectbox{\displaystyle\vec{}\mkern 4.0mu}\cr\kern-4.30554pt\cr\displaystyle n\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{\textstyle\vec{}\mkern 4.0mu}\cr\kern-4.30554pt\cr\textstyle n\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{\scriptstyle\vec{}\mkern 4.0mu}\cr\kern-3.01389pt\cr\scriptstyle n\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{\scriptscriptstyle\vec{}\mkern 4.0mu}\cr\kern-2.15277pt\cr\scriptscriptstyle n\cr}}}}_{i}, implicitly indicating the size $I_{n_{i}}$; the tensor in Eq.(17) would then be rewritten as \mathcal{O}=\{\mathbf{O}^{\vec{n}_{1},\vec{n}_{2}}_{{\mathchoice{\vbox{\offinterlineskip\halign{#\cr\reflectbox{\displaystyle\vec{}\mkern 4.0mu}\cr\kern-3.01389pt\cr\displaystyle n\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{\textstyle\vec{}\mkern 4.0mu}\cr\kern-3.01389pt\cr\textstyle n\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{\scriptstyle\vec{}\mkern 4.0mu}\cr\kern-2.10971pt\cr\scriptstyle n\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{\scriptscriptstyle\vec{}\mkern 4.0mu}\cr\kern-1.50694pt\cr\scriptscriptstyle n\cr}}}}_{3},{\mathchoice{\vbox{\offinterlineskip\halign{#\cr\reflectbox{\displaystyle\vec{}\mkern 4.0mu}\cr\kern-3.01389pt\cr\displaystyle n\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{\textstyle\vec{}\mkern 4.0mu}\cr\kern-3.01389pt\cr\textstyle n\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{\scriptstyle\vec{}\mkern 4.0mu}\cr\kern-2.10971pt\cr\scriptstyle n\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{\scriptscriptstyle\vec{}\mkern 4.0mu}\cr\kern-1.50694pt\cr\scriptscriptstyle n\cr}}}}_{4}}\}.

II.2 Adaptively symmetric Matrix Product States

An MPS representing a quantum state on $L$ sites is a chain with $3-$dimensional tensors. The tensor $\mathcal{M}_{l}$ on site $l$ of the MPS is labelled by a physical quantum number $\sigma_{l}$, an auxiliary quantum number $a_{l}$ which connects it to its left hand site, and $a_{l+1}$ which connects it to its right hand site. As a result, $\mathcal{M}_{l}$ can be written as

[TABLE]

and the fusion rule Eq.(18) becomes

[TABLE]

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

[TABLE]

As for typical ${\rm s}\text{-}{\rm MPS}$, the overall bond dimension (defined as the size of the auxiliary dimensions) of an ${\rm as}\text{-}{\rm MPS}$ 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 ${\rm as}\text{-}{\rm MPS}$.

II.3 Adaptively symmetric Matrix Product Operators

Adaptively symmetric Matrix Product Operators (${\rm as}\text{-}{\rm MPO}$) can be treated in a similar manner as they form a chain of $4-$dimensional symmetry protected tensors. The tensor $\mathcal{W}_{l}$ on the $l$-th site of an ${\rm as}\text{-}{\rm MPO}$ is labelled by two physical quantum numbers $\sigma_{l}$ and $\tau_{l}$, an auxiliary quantum number $b_{l}$ which connects it to its left hand site, and $b_{l+1}$ which connects it to its right hand site. Thus $\mathcal{W}_{l}$ can be written as

[TABLE]

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

[TABLE]

The flow of quantum numbers is represented in Fig.1 (top). For $U(1)$ symmetric systems, the MPO does not change the total quantum number of the quantum state. If the quantum state had total quantum number $N$ before being multiplied by the MPO, after the multiplication between MPO and MPS, the total quantum number for the quantum state is still $N$. As a result we can denote a symmetric MPO as ${\rm s}\text{-}{\rm MPO}=\sum_{N}{\rm s}\text{-}{\rm MPO}_{N}^{N}$, with the particular number preserving condition that ${b_{1}}={b_{L+1}}=0$. In order to generalize this concept, we allow the quantum number $b_{1}$ at the left boundary to take multiple values while still keeping $b_{L+1}=0$. Therefore the adaptively symmetric Matrix Product Operator can be written as

[TABLE]

where $N^{\prime}=N$ if and only if the system has $U(1)$ symmetry. In case the system has only a discrete symmetry, for example parity, a $Z_{2}$ symmetry, then $N$ and $N^{\prime}$ satisfy $N^{\prime}=N\pm E$ where $E$ is an even number. If the system has no symmetry at all, then $N$ and $N^{\prime}$ can be arbitrary. However, in this case, it may be more efficient to use non-symmetric ${\rm MPO}$ and ${\rm MPS}$ 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 ${\rm as}\text{-}{\rm MPS}$ and ${\rm as}\text{-}{\rm MPO}$.

II.4 Ground state search algorithm

Here we present the ground state search algorithm that we have implemented based on ${\rm as}\text{-}{\rm MPS}$. The core structure retraces that of a two-site 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

[TABLE]

where the $\mathcal{A}_{l}$ are called “left-canonical” because

[TABLE]

and the $\mathcal{B}_{l}$ are “right-canonical”, i.e.

[TABLE]

In both Eqs.(90,123) \mathcal{I}_{\vec{a}_{l},{\mathchoice{\vbox{\offinterlineskip\halign{#\cr\reflectbox{\displaystyle\vec{}\mkern 4.0mu}\cr\kern-3.01389pt\cr\displaystyle a\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{\textstyle\vec{}\mkern 4.0mu}\cr\kern-3.01389pt\cr\textstyle a\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{\scriptstyle\vec{}\mkern 4.0mu}\cr\kern-2.10971pt\cr\scriptstyle a\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{\scriptscriptstyle\vec{}\mkern 4.0mu}\cr\kern-1.50694pt\cr\scriptscriptstyle a\cr}}}}_{l}^{\prime}} represents the identity matrix for the two indices $\vec{a}_{l}$, {\mathchoice{\vbox{\offinterlineskip\halign{#\cr\reflectbox{\displaystyle\vec{}\mkern 4.0mu}\cr\kern-4.30554pt\cr\displaystyle a\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{\textstyle\vec{}\mkern 4.0mu}\cr\kern-4.30554pt\cr\textstyle a\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{\scriptstyle\vec{}\mkern 4.0mu}\cr\kern-3.01389pt\cr\scriptstyle a\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{\scriptscriptstyle\vec{}\mkern 4.0mu}\cr\kern-2.15277pt\cr\scriptscriptstyle a\cr}}}}_{l}^{\prime}. Here ${\rm conjugate}(\mathcal{O})$ means to take the element wise conjugate of the tensor $\mathcal{O}$, during which the direction of each quantum number also gets reversed. The two tenors with a commom quantum number $n$ which is summed over should have opposite directions. No restrictions is made on the tensors $\mathcal{M}$. A two-site ground state search variational MPS algorithm works by iteratively minimizing the energy of an effective Hamiltonian ${{\hat{\rm H}_{\rm eff}}}$ 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 ${{\hat{\rm H}_{\rm eff}}}$. To compute ${{\hat{\rm H}_{\rm eff}}}$ on sites $l$ and $l+1$ (also referred to as bond $l$), one needs to trace out the left environment consisting of the sites to the left of the $l$-th site to obtain the tensor $\mathcal{L}_{l}$

[TABLE]

and the right environment consisting of sites to the right of the $\left(l+1\right)$-th site to obtain the tensor $\mathcal{R}_{l+2}$

[TABLE]

With $\mathcal{L}_{l}$, $\mathcal{R}_{l+2}$ and $\mathcal{W}_{l}$, $\mathcal{W}_{l+1}$, one can construct ${{\hat{\rm H}_{\rm eff}}}$ as

[TABLE]

Grouping the tensor indexes ({\mathchoice{\vbox{\offinterlineskip\halign{#\cr\reflectbox{\displaystyle\vec{}\mkern 4.0mu}\cr\kern-4.30554pt\cr\displaystyle a\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{\textstyle\vec{}\mkern 4.0mu}\cr\kern-4.30554pt\cr\textstyle a\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{\scriptstyle\vec{}\mkern 4.0mu}\cr\kern-3.01389pt\cr\scriptstyle a\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{\scriptscriptstyle\vec{}\mkern 4.0mu}\cr\kern-2.15277pt\cr\scriptscriptstyle a\cr}}}}_{l},\vec{\sigma}_{l},\vec{\sigma}_{l+1},\vec{a}_{l+2}) and (\vec{a}_{l}^{\prime},{\mathchoice{\vbox{\offinterlineskip\halign{#\cr\reflectbox{\displaystyle\vec{}\mkern 4.0mu}\cr\kern-4.30554pt\cr\displaystyle\sigma\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{\textstyle\vec{}\mkern 4.0mu}\cr\kern-4.30554pt\cr\textstyle\sigma\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{\scriptstyle\vec{}\mkern 4.0mu}\cr\kern-3.01389pt\cr\scriptstyle\sigma\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{\scriptscriptstyle\vec{}\mkern 4.0mu}\cr\kern-2.15277pt\cr\scriptscriptstyle\sigma\cr}}}}_{l}^{\prime},{\mathchoice{\vbox{\offinterlineskip\halign{#\cr\reflectbox{\displaystyle\vec{}\mkern 4.0mu}\cr\kern-4.30554pt\cr\displaystyle\sigma\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{\textstyle\vec{}\mkern 4.0mu}\cr\kern-4.30554pt\cr\textstyle\sigma\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{\scriptstyle\vec{}\mkern 4.0mu}\cr\kern-3.01389pt\cr\scriptstyle\sigma\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{\scriptscriptstyle\vec{}\mkern 4.0mu}\cr\kern-2.15277pt\cr\scriptscriptstyle\sigma\cr}}}}_{l+1}^{\prime},{\mathchoice{\vbox{\offinterlineskip\halign{#\cr\reflectbox{\displaystyle\vec{}\mkern 4.0mu}\cr\kern-4.30554pt\cr\displaystyle a\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{\textstyle\vec{}\mkern 4.0mu}\cr\kern-4.30554pt\cr\textstyle a\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{\scriptstyle\vec{}\mkern 4.0mu}\cr\kern-3.01389pt\cr\scriptstyle a\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{\scriptscriptstyle\vec{}\mkern 4.0mu}\cr\kern-2.15277pt\cr\scriptscriptstyle a\cr}}}}_{l+2}^{\prime}), one could treat ${{\hat{\rm H}_{\rm eff}}}$ as a matrix and find its lowest eigenstate \mathcal{M}_{\vec{a}_{l},{\mathchoice{\vbox{\offinterlineskip\halign{#\cr\reflectbox{\displaystyle\vec{}\mkern 4.0mu}\cr\kern-3.01389pt\cr\displaystyle\sigma\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{\textstyle\vec{}\mkern 4.0mu}\cr\kern-3.01389pt\cr\textstyle\sigma\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{\scriptstyle\vec{}\mkern 4.0mu}\cr\kern-2.10971pt\cr\scriptstyle\sigma\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{\scriptscriptstyle\vec{}\mkern 4.0mu}\cr\kern-1.50694pt\cr\scriptscriptstyle\sigma\cr}}}}_{l},{\mathchoice{\vbox{\offinterlineskip\halign{#\cr\reflectbox{\displaystyle\vec{}\mkern 4.0mu}\cr\kern-3.01389pt\cr\displaystyle\sigma\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{\textstyle\vec{}\mkern 4.0mu}\cr\kern-3.01389pt\cr\textstyle\sigma\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{\scriptstyle\vec{}\mkern 4.0mu}\cr\kern-2.10971pt\cr\scriptstyle\sigma\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{\scriptscriptstyle\vec{}\mkern 4.0mu}\cr\kern-1.50694pt\cr\scriptscriptstyle\sigma\cr}}}}_{l+1},{\mathchoice{\vbox{\offinterlineskip\halign{#\cr\reflectbox{\displaystyle\vec{}\mkern 4.0mu}\cr\kern-3.01389pt\cr\displaystyle a\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{\textstyle\vec{}\mkern 4.0mu}\cr\kern-3.01389pt\cr\textstyle a\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{\scriptstyle\vec{}\mkern 4.0mu}\cr\kern-2.10971pt\cr\scriptstyle a\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{\scriptscriptstyle\vec{}\mkern 4.0mu}\cr\kern-1.50694pt\cr\scriptscriptstyle a\cr}}}}_{l+2}} and eigenvalue. \mathcal{M}_{{\mathchoice{\vbox{\offinterlineskip\halign{#\cr\reflectbox{\displaystyle\vec{}\mkern 4.0mu}\cr\kern-3.01389pt\cr\displaystyle a\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{\textstyle\vec{}\mkern 4.0mu}\cr\kern-3.01389pt\cr\textstyle a\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{\scriptstyle\vec{}\mkern 4.0mu}\cr\kern-2.10971pt\cr\scriptstyle a\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{\scriptscriptstyle\vec{}\mkern 4.0mu}\cr\kern-1.50694pt\cr\scriptscriptstyle a\cr}}}}_{l},\vec{\sigma}_{l},\vec{\sigma}_{l+1},\vec{a}_{l+2}} is then decomposed into two $3$ dimensional tensors by singular value decomposition (SVD)

[TABLE]

and the results are used to update $\mathcal{A}_{l}$ or $\mathcal{B}_{l+1}$ depending on whether the sweep is from left to right or from right to left.

To adapt the above algorithm to ${\rm as}\text{-}{\rm MPS}$, one notices that the only difference between ${\rm s}\text{-}{\rm MPS}$, ${\rm s}\text{-}{\rm MPO}$ and ${\rm as}\text{-}{\rm MPS}$, ${\rm as}\text{-}{\rm MPO}$ is that for the latter there is a much larger variety of possible values for the left-most index ($a_{1}$ or $b_{1}$) of the tensors at the left boundary. Conventionally, ${a_{1}}$ (auxiliary index on the ${\rm s}\text{-}{\rm MPS}$) can only take a single value $N$ and ${b_{1}}$ (auxiliary index on the ${\rm s}\text{-}{\rm MPO}$) has to be [math], therefore one can straightforwardly construct $\mathcal{L}_{1}$ as

[TABLE]

where 1 means a trivial dense tensor with a single element $1$ and where we have used the label ${\rm s}\text{-}{\rm MPS}$ to remind the reader that this concerns symmetric Matrix Product States. For ${\rm as}\text{-}{\rm MPS}$ and ${\rm as}\text{-}{\rm MPO}$, instead, $\mathcal{L}_{1}$ must be able to include all possible quantum numbers that could appear after the ${\rm as}\text{-}{\rm MPO}$ has been applied on the ${\rm as}\text{-}{\rm MPS}$. Thus a way to construct $\mathcal{L}_{1}$ is to look at the boundary indexes ${a_{1}}$ and ${b_{1}}$ and find all the relevant combinations $a_{1}+b_{1}$, namely

[TABLE]

With such a construction of $\mathcal{L}_{1}$, the variational ground state search algorithm would be able to automatically detect and preserve the symmetry in the system, whether it is $U(1)$ or, for example, $Z_{2}$.

II.5 Time evolution algorithm

Time evolution with ${\rm as}\text{-}{\rm MPS}$ 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 ${\rm as}\text{-}{\rm MPS}$. 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 ${\rm as}\text{-}{\rm MPS}$. In these algorithms, one splits the time evolution operator $\hat{{\rm U}}(dt)$ 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 $dt$ is divided in evolution on even and odd bonds (respectively with the Hamiltonians $\hat{H}_{e}$ and $\hat{H}_{o}$), and which can be parallelized (for more details on this $t-$MPS algorithm see, for instance, the review Schollwock2011 ). In the following, for clarity, we focus on a unitary evolution with a Hamiltonian $\hat{{\rm H}}$, however in the example discussed in Sec.III.2 we will show how this approach can also be readily applied to dissipative evolutions. For ${\rm as}\text{-}{\rm MPS}$, 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 $\hat{h}_{l}=\hat{\sigma}^{+}_{l}+\hat{\sigma}^{-}_{l}$ acting on the $l$-th site of an ${\rm as}\text{-}{\rm MPS}$ initially in a single number sector $N$, the resulting state will be a superposition of three number sectors $N-1,N,N+1$. This means that even if $\hat{h}_{l}$ only acts locally on site $l$, the left boundary tensor at site $l=1$ of the ${\rm as}\text{-}{\rm MPS}$ has to be updated and thus the effect on ${\rm as}\text{-}{\rm MPS}$ 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 $\hat{{\rm H}}$ into a symmetric, number-conserving part, $\hat{{\rm H}}_{s}$, and an asymmetric, non-number-conserving part, $\hat{{\rm H}}_{a}$:

[TABLE]

It is now possible to use a second order time evolution operator as (setting $\hbar=1$)

[TABLE]

as shown in Fig.2. For the symmetric portion, $\exp\left(-{\rm i}\hat{{\rm H}}_{s}dt\right)$, one can perform a Suzuki-Trotter based algorithm, while for $\exp\left(-{\rm i}\hat{{\rm H}}_{a}dt\right)$, 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 $D_{W}$ of $\exp\left(-{\rm i}\hat{{\rm H}}_{a}dt\right)$ is small which, for instance, is the case if the non-number-conserving term $\hat{{\rm H}}_{a}$ 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 ${\rm as}\text{-}{\rm 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 ${\rm XYZ}$ spin chain of size $L$ can be written as

[TABLE]

where $\hat{\sigma}^{x}_{l}$, $\hat{\sigma}^{y}_{l}$ and $\hat{\sigma}^{z}_{l}$ are the operators corresponding to Pauli matrices, $J_{{\rm XY}}$ is the tunneling, $\Delta$ is the interaction strength and $\gamma$ denotes the anisotropy ranging from [math] to $1$, $h$ is the strength of magnetization. For the special case $\gamma=0$, this model reduces to an ${\rm XXZ}$ chain with longitudinal field which conserves the total spin $\hat{S}^{z}_{T}=\sum_{l}\hat{\sigma}^{z}_{l}$ and has $U(1)$ symmetry, while for $\gamma\neq 0$, the Hamiltonian commutes with the parity operator $\mathcal{P}$ defined as

[TABLE]

thus having a discrete $Z_{2}$ symmetry. In the following we work in units for which $J_{{\rm XY}}=\hbar=1$. We first apply the ground state search algorithm to $\hat{H}_{{\rm XYZ}}$ with $L=50$. In the following we keep the interaction strength $\Delta=1.5$ and magnetization strength $h=0.5$, and we compute the ground state for different values of $\gamma$, starting the search of the ground state from a state with zero total magnetization $\hat{S}^{Z}_{T}=0$. For the case $\gamma=0$, 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 [math] total spin. However, for $\gamma\neq 0$, 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 $Z_{2}$ symmetric MPS for the $\gamma\neq 0$ case, one will obtain the same ground state but only in the number sector labelled by [math] (corresponding to $\mod\!(\hat{S}^{z}_{T},2)=0$ , i.e. total spin modulo $2$ is [math]). A major difference between using a $Z_{2}$ symmetric MPS and our ${\rm as}\text{-}{\rm 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 an adaptively symmetric MPS into many $U(1)$ symmetric MPSs with a fixed total spin as in Eq.(56), by just splitting the tensor on the left boundary. Therefore, with adaptively symmetric MPS, one can easily compute the distribution of the ground state in different number sectors $P_{S^{z}_{T}}$, which is shown in Fig.3(a). One can see clearly that for $\gamma=0$, the ground state has a single total spin (a kronecker delta at [math] represented by a dashed line), while the distribution becomes broader as $\gamma$ increases to $1$. we can also see that the non-zero values are in a period of $4$, which means the total number of spins pointing up (down) is even.

We benckmark our ${\rm as}\text{-}{\rm MPS}$ and ${\rm as}\text{-}{\rm MPO}$ based ground state search algorithm with the conventional ground state search algorithm based on a $Z_{2}$ ${\rm s}\text{-}{\rm MPS}$ and ${\rm s}\text{-}{\rm MPO}$. For the parameters which we have considered in Fig.3(a), the differences between the ground state energies computed from both approaches are of the order $10^{-9}$, with a bond dimension $D=2000$. These two approaches use around $1.7$ and $2.5$ Gb memory, $0.75$ and $2$ hours separately on a personal computer with $4$ cores of $3.7$ GHz frequency.

We then do a quench, starting from the ground state corresponding to $\gamma=0$, $|GS_{{\rm XXZ}}\rangle$, and then evolving it with the Hamiltonian (447) for $\gamma=0.5$, as if a sudden quench occurred. This evolution can be represented as

[TABLE]

We chose to perform the evolution using a fourth order Runge-Kutta method based on MPO$-$MPS multiplication. To ensure convergence of our results, we have done simulations with bond dimensions $D=1000,1500,2000$ and checked various observables, for instance, for the average spin occupation $\langle\Psi(t)|\hat{\sigma}^{z}_{l}|\Psi(t)\rangle$ we get a difference of the order of $10^{-9}$. We have also compared these observables with those computed with time-evolving matrix product states algorithm (t-MPS) Schollwock2011 using a $Z_{2}$ symmetric MPS and MPO, obtaining a difference of the order $10^{-4}$. We note that this difference is partially due to the fact that we have used a simple MPO based Runge-Kutta method for the time evolution, which is known to have poorer performance than other MPO based time evolution algorithms such as the approaches in ZaletelPollmann2014 . Instead, for the symmetric MPS with a $Z_{2}$ symmetry, we have used the standard t-MPS algorithm (by manually converting the initial state from $U(1)$ symmetry into $Z_{2}$ symmetry) which is known to be very stable and accurate. With the ${\rm as}\text{-}{\rm MPS}$ 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 $S^{z}_{T}$, 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 ${\rm as}\text{-}{\rm MPS}$ 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

[TABLE]

where we refer to $\hat{\mathcal{L}}$ as the Lindbladian, $\hat{\mathcal{D}}$ as dissipator, and where the Hamiltonian $\hat{H}_{{\rm BH}}$ can be written as

[TABLE]

The dissipator can be written as

[TABLE]

In Eqs.(451,452) we have used the notations $\hat{a}^{\dagger}_{l},\hat{a}_{l}$ for the bosonic creation and annihilation operators and $\hat{n}_{l}=\hat{a}^{\dagger}_{l}\hat{a}_{l}$. $J$ is the tunneling amplitude and $U$ is the local interaction strength. $\Lambda_{1,L}$ denotes the strength of the dissipation acting at the edges $l=1$ or $l=L$. If $J=0$, the dissipator will impose a local thermal distribution at the edges with average occupation $\bar{n}_{1,L}$ respectively. Similarly as for the case of the ${\rm XYZ}$ chain, in the following we work in units for which $J=\hbar=1$.

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 $d$. For unitary systems, for an interaction strength $U\geq 3$, it is usually sufficient to keep $d=4$. However, to faithfully represent a local bosonic thermal state, $d$ has to be much larger due to the long tail of the thermal distribution. In our simulations we have used $d=10$ and $d=11$ 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 $d^{2}$. If one uses non-symmetric Matrix Product States for the time evolution, for instance splitting the operator $\exp\left(\hat{\mathcal{L}}dt\right)$ into many two-sites operators, then each two-body operator would be a $d^{4}\times d^{4}$ matrix. For $d=10$, one such operator will contain $10^{8}$ (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 $Dd^{2}\times Dd^{2}$ matrix. For a reasonable bond dimension $D=100$ and local Hilbert space $d=10$, this tensor would have a size of $10^{4}\times 10^{4}$, 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 $d^{2}$ to $kd$ with $k<5$. This was possible for the dissipator used (same as the one in Eq.(452), and the system analyzed), as the local density matrix is mainly diagonal. Such approximation is however not sufficient to study larger systems. Using ${\rm as}\text{-}{\rm MPS}$ it is possible, as explained in the following, to study larger systems.

We note that the unitary part of $\hat{\mathcal{L}}$ is number conserving, while the dissipator in Eq.(452) is local, thus the MPO for the non-number-conserving part has a bond dimension equal to $1$, which is the smallest possible size. In addition, the dissipator in Eq.(452) has an ulterior property: if we write the density operator using the basis $|\vec{n}_{N};N\rangle\langle\vec{m}_{M};M|$ where $N$ is the total number of bosons in the ket ($M$ in the bra), and $\vec{n}_{N}$ is the vector detailing how the $N$ bosons are distributed over the $L$ sites (respectively $\vec{m}_{M}$ describes how the $M$ bosons are distributed between $L$ sites), then it appears clear that the dissipator couples the element corresponding to $|\vec{n}_{N};N\rangle\langle\vec{m}_{M};M|$ with $|\vec{n}_{N\pm 1};N\pm 1\rangle\langle\vec{m}_{M\pm 1};M\pm 1|$. This implies that if the initial condition belongs to a single number sector (e.g. at initial time $N=M$), then the density operator will have only non-zero terms for blocks in which $N=M$. 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 ${\rm as}\text{-}{\rm MPS}$ and ${\rm as}\text{-}{\rm MPO}$, 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 $L=20$ with average filling $0.5$ and $U=4$, which we denote as $|GS_{{\rm BH}}\rangle$. We then turn on the dissipation, with $\Lambda_{1}=\Lambda_{L}=1$, $\bar{n}_{1}=0.75$ and $\bar{n}_{L}=0.25$. For the simulations we have truncated the local Hilbert space to $d=10$ and used a bond dimension $D=1000$ for the MPS, and the time step $dt=0.01$. To ensure that our choices of $d$ and $D$ are appropriate, we have also done simulations and computed the local density $\langle\hat{n}_{j}\rangle={\rm tr}\left(\hat{n}_{l}\;\hat{\rho}\right)$ with $d=11$, getting a difference of the order of $10^{-5}$, and with $D=1500$, getting a difference of the order of $10^{-4}$. In Fig.4(a), we plot the distribution of the quantum state in different number sectors, i.e.

[TABLE]

something readily done with ${\rm as}\text{-}{\rm MPS}$. The initial state has a total of $10$ 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 $\langle\hat{n}_{l}\rangle$ as a function of site $l$ at different times. Initially, since the evolution starts from the ground state, the distribution of $\langle n_{l}\rangle$ 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 treat a global $U(1)$ symmetry or one of its subgroups on the same footing. At the same time, this method can also deal with non-number conserving systems. We have shown how this method, based on adaptively symmetric Matrix Product States, can be applied to search for the ground state of a system 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 $U(1)$ symmetry, or, for instance, with only $Z_{2}$ symmetry, the adaptively 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 adaptively 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. Adaptively 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).

Appendix A Construction of ${\rm as}\text{-}{\rm MPS}$ and ${\rm as}\text{-}{\rm MPO}$ for XYZ chain

In this appendix we demonstrate the process of constructing of an ${\rm as}\text{-}{\rm MPS}$ and an ${\rm as}\text{-}{\rm MPO}$ using the concrete example of a $3$-site spin XYZ chain with the magnetization strength $h=0$. We label the state $|0\rangle$ with quantum number [math], and the state $|1\rangle$ with $1$. Assuming we have an initial state which is

[TABLE]

The ${\rm as}\text{-}{\rm MPS}$ corresponding to state $|\psi\rangle$ is

[TABLE]

In fact on site $l=3$ there is only one possible quantum number flowing in {\mathchoice{\vbox{\offinterlineskip\halign{#\cr\reflectbox{\displaystyle\vec{}\mkern 4.0mu}\cr\kern-4.30554pt\cr\displaystyle a\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{\textstyle\vec{}\mkern 4.0mu}\cr\kern-4.30554pt\cr\textstyle a\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{\scriptstyle\vec{}\mkern 4.0mu}\cr\kern-3.01389pt\cr\scriptstyle a\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{\scriptscriptstyle\vec{}\mkern 4.0mu}\cr\kern-2.15277pt\cr\scriptscriptstyle a\cr}}}}_{4}=0, only one quantum number present on the site {\mathchoice{\vbox{\offinterlineskip\halign{#\cr\reflectbox{\displaystyle\vec{}\mkern 4.0mu}\cr\kern-4.30554pt\cr\displaystyle\sigma\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{\textstyle\vec{}\mkern 4.0mu}\cr\kern-4.30554pt\cr\textstyle\sigma\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{\scriptstyle\vec{}\mkern 4.0mu}\cr\kern-3.01389pt\cr\scriptstyle\sigma\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{\scriptscriptstyle\vec{}\mkern 4.0mu}\cr\kern-2.15277pt\cr\scriptscriptstyle\sigma\cr}}}}_{3}=0 and hence only one possible output $\vec{a}_{3}=0$. Similarly, at site $l=2$ there are two possible values for {\mathchoice{\vbox{\offinterlineskip\halign{#\cr\reflectbox{\displaystyle\vec{}\mkern 4.0mu}\cr\kern-4.30554pt\cr\displaystyle\sigma\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{\textstyle\vec{}\mkern 4.0mu}\cr\kern-4.30554pt\cr\textstyle\sigma\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{\scriptstyle\vec{}\mkern 4.0mu}\cr\kern-3.01389pt\cr\scriptstyle\sigma\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{\scriptscriptstyle\vec{}\mkern 4.0mu}\cr\kern-2.15277pt\cr\scriptscriptstyle\sigma\cr}}}}_{2}=0,\;1 which results in $\vec{a}_{2}=0,\;1$. At site $l=3$, the local quantum number is only {\mathchoice{\vbox{\offinterlineskip\halign{#\cr\reflectbox{\displaystyle\vec{}\mkern 4.0mu}\cr\kern-4.30554pt\cr\displaystyle\sigma\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{\textstyle\vec{}\mkern 4.0mu}\cr\kern-4.30554pt\cr\textstyle\sigma\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{\scriptstyle\vec{}\mkern 4.0mu}\cr\kern-3.01389pt\cr\scriptstyle\sigma\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{\scriptscriptstyle\vec{}\mkern 4.0mu}\cr\kern-2.15277pt\cr\scriptscriptstyle\sigma\cr}}}}_{1}=1, but since {\mathchoice{\vbox{\offinterlineskip\halign{#\cr\reflectbox{\displaystyle\vec{}\mkern 4.0mu}\cr\kern-4.30554pt\cr\displaystyle a\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{\textstyle\vec{}\mkern 4.0mu}\cr\kern-4.30554pt\cr\textstyle a\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{\scriptstyle\vec{}\mkern 4.0mu}\cr\kern-3.01389pt\cr\scriptstyle a\cr}}}{\vbox{\offinterlineskip\halign{#\cr\reflectbox{\scriptscriptstyle\vec{}\mkern 4.0mu}\cr\kern-2.15277pt\cr\scriptscriptstyle a\cr}}}}_{2} takes two possible values, then $\vec{a}_{1}=1,\;2$.

In the non-symmetric case, the MPO corresponding to $\hat{H}_{{\rm XYZ}}$ for $L=3$ can be straightforwardly written as

[TABLE]

Now we first rewrite $\hat{H}_{{\rm XYZ}}$ in the following form

[TABLE]

which is a summation of $10$ terms (because there are two bonds). To write down the corresponding ${\rm as}\text{-}{\rm MPO}_{\rm XYZ}$, it is convenient to first write the product MPO for each term and then use the MPO addition rule to sum them up ClaudiusSchollwock2017 . The product MPO corresponding to the terms $\hat{\sigma}^{+}_{1}\hat{\sigma}^{-}_{2}$, $\hat{\sigma}^{-}_{1}\hat{\sigma}^{+}_{2}$, $\hat{\sigma}^{+}_{1}\hat{\sigma}^{+}_{2}$, $\hat{\sigma}^{-}_{1}\hat{\sigma}^{-}_{2}$, $\hat{\sigma}^{z}_{1}\hat{\sigma}^{z}_{2}$ can be written as

[TABLE]

We notice that for the terms ${\rm MPO}_{\hat{\sigma}^{+}_{1}\hat{\sigma}^{-}_{2}}$, ${\rm MPO}_{\hat{\sigma}^{-}_{1}\hat{\sigma}^{+}_{2}}$, ${\rm MPO}_{\hat{\sigma}^{z}_{1}\hat{\sigma}^{z}_{2}}$, the auxiliary index on the left boundary ${b_{1}}=0$, while we have ${b_{1}}=2$ for the term ${\rm MPO}_{\hat{\sigma}^{+}_{1}\hat{\sigma}^{+}_{2}}$, and ${b_{1}}=-2$ for the term ${\rm MPO}_{\hat{\sigma}^{-}_{1}\hat{\sigma}^{-}_{2}}$. We can also construct the remaining terms of ${\rm MPO}_{\hat{\sigma}^{+}_{2}\hat{\sigma}^{-}_{3}}$, ${\rm MPO}_{\hat{\sigma}^{-}_{2}\hat{\sigma}^{+}_{3}}$, ${\rm MPO}_{\hat{\sigma}^{+}_{2}\hat{\sigma}^{+}_{3}}$, ${\rm MPO}_{\hat{\sigma}^{-}_{2}\hat{\sigma}^{-}_{3}}$, ${\rm MPO}_{\hat{\sigma}^{z}_{2}\hat{\sigma}^{z}_{3}}$ similarly. We note that if we encode the MPO with $Z_{2}$ symmetry instead of $U(1)$ symmetry, then ${b_{1}}$ for the terms ${\rm MPO}_{\hat{\sigma}^{+}_{1}\hat{\sigma}^{+}_{2}}$ and ${\rm MPO}_{\hat{\sigma}^{-}_{1}\hat{\sigma}^{-}_{2}}$ will also be [math] since $2$ and $-2$ are equivalent to [math] modulo $2$. The final ${\rm as}\text{-}{\rm MPO}$ is a summation of all these terms

[TABLE]