Automated construction of $U(1)$-invariant matrix-product operators from graph representations

TL;DR

This paper introduces an automated method to construct $U(1)$-invariant matrix-product operators from finite-states machine representations, simplifying the process and reducing computational costs in tensor network algorithms.

Contribution

The authors develop an algorithmic scheme that automates the construction of $U(1)$-invariant MPOs from FSM representations, handling quantum number shifts and phase factors efficiently.

Findings

Automated MPO construction reduces computational complexity.

Exact MPO representations for Hamiltonian variances are generated.

Method applicable to various $U(1)$-invariant operators.

Abstract

We present an algorithmic construction scheme for matrix-product-operator (MPO) representations of arbitrary $U(1)$-invariant operators whenever there is an expression of the local structure in terms of a finite-states machine (FSM). Given a set of local operators as building blocks, the method automatizes two major steps when constructing a $U(1)$-invariant MPO representation: (i) the bookkeeping of auxiliary bond-index shifts arising from the application of operators changing the local quantum numbers and (ii) the appearance of phase factors due to particular commutation rules. The automatization is achieved by post-processing the operator strings generated by the FSM. Consequently, MPO representations of various types of $U(1)$-invariant operators can be constructed generically in MPS algorithms reducing the necessity of expensive MPO arithmetics. This is demonstrated by generating arbitrary products of operators in terms of FSM, from which we obtain exact MPO representations for the variance of the Hamiltonian of a $S=1$ Heisenberg chain.

Figures

Equations

$$\mathcal{H}$$

$$\hat{O}$$

$$\left[\hat{N},\hat{O}\right]$$

$$|\psi ⟩$$

$$T^{n_i}$$

$${\hat{A}}^{(i)}{\hat{B}}^{(j)}|\psi ⟩\equiv \sum _{n_1\cdots n_L}\sum _{n_i^{\prime },n_j^{\prime }}$$

$$\times \delta (n_i+n_j-(n_i^{\prime }+n_j^{\prime }))|n_1\cdots n_L⟩,$$

$${\hat{A}}^{(i)}{\hat{B}}^{(j)}|\psi ⟩=\sum _{n_i^{\prime },n_j^{\prime }}\{\cdots T^{n_{i-1}}[\sum _{\tau }$$

$$\times$$

$${\left({\hat{S}}^{\pm }\right)}^{n_i^{},n_i^{\prime }}\ne 0\iff n_i^{}-n_i^{\prime }=\pm 1\equiv \Delta ,$$

$${\hat{A}}^{(i)}$$

$${\hat{B}}^{(j)}$$

$$\displaystyle\left(A_{1},\cdots,A_{n}\right)\left(\begin{array}[]{ccc}D&\cdots&0\\ \vdots&\ddots&\vdots\\ 0&\cdots&D\end{array}\right)\left(\begin{array}[]{c}B_{1}\\ \vdots\\ B_{n}\end{array}\right)$$

$${\hat{R}}_{\Delta }^{(i)}$$

$${\hat{A}}^{(i)}{\hat{B}}^{(j)}|\psi ⟩$$

$$V^{n_i}$$

$$P^{n_k}$$

$$W^{n_j}$$

$${\hat{T}}^{(k)}$$

$${\hat{T}}_{a+{\Delta }_1+\cdots +{\Delta }_l}^{(k)}$$

$$\hat{H}=\sum _{\nu ,r}{\hat{H}}_{\nu ,r}=\sum _{\nu ,r}\sum _i{\hat{h}}_{\nu ,r}^{(i)},$$

$${\hat{H}}_{XX}$$

$$\hat{H}=\sum _{\nu ,r}{\hat{H}}_{\nu ,r}.$$

$$\oplus :\Lambda ({\hat{H}}_1+{\hat{H}}_2)=\Lambda ({\hat{H}}_1)\oplus \Lambda ({\hat{H}}_2).$$

$${\Lambda }_{\max }(\hat{H})=\underset{\nu ,r}{⨁}\Lambda ({\hat{H}}_{\nu ,r}).$$

$${\hat{o}}_{{\nu }_1}\dots {\hat{o}}_{{\nu }_r}⟶\left[{\hat{u}}_{{\nu }_1}{\hat{o}}_{{\nu }_1}\right]\dots \left[{\hat{u}}_{{\nu }_r}{\hat{o}}_{{\nu }_r}\right].$$

$${\hat{S}}_{(i)}^{+}\otimes {\widehat{\text{Id}}}^{(i+1)}\otimes {\hat{S}}_{(i+2)}^{-}⟶{\hat{S}}_{(i)}^{+}\otimes \left[{\hat{R}}_{+1}^{(i+1)}{\widehat{\text{Id}}}^{(i+1)}\right]\otimes \left[{\hat{R}}_{+1}^{(i+2)}{\hat{S}}_{(i+2)}^{-}\right]$$

$${\hat{a}}_{(i)}\to {\hat{a}}_{(i)}{\mathrm{e}}^{i\pi {\sum }_{j<i}{\hat{n}}^{(j)}}=\prod _{j<i}{\hat{\Pi }}^{(j)}{\hat{a}}_{(i)},$$

$${\hat{a}}_{(i)}^{†}\otimes {\widehat{\text{Id}}}^{(i+1)}\otimes {\hat{a}}_{(i+2)}^{}⟶\left[{\hat{\Pi }}^{(i)}{\hat{a}}_{(i)}^{†}\right]\otimes \left[{\hat{\Pi }}^{(i+1)}{\hat{R}}_{+1}^{(i+1)}{\widehat{\text{Id}}}^{(i+1)}\right]\otimes \left[{\hat{R}}_{+1}^{(i+2)}{\hat{a}}_{(i+2)}^{}\right].$$

$${\hat{H}}_1\cdot {\hat{H}}_2$$

$${\hat{H}}_{{\nu }^1,r^1}\cdot {\hat{H}}_{{\nu }^2,r^2}$$

$$\Lambda ({\hat{H}}_{{\nu }^1,r^1})\wedge \Lambda ({\hat{H}}_{{\nu }^2,r^2})$$

$$\Lambda ({\hat{H}}_{r^1,{\nu }^1}\cdot {\hat{H}}_{r^2,{\nu }^2})$$

Full text

Automated construction of $U(1)$-invariant matrix-product operators from graph representations

S. Paeckel*, T. Köhler, and S.R. Manmana

Institut für Theoretische Physik, Universität Göttingen, 37077 Göttingen, Germany

• sebastian.paeckel@theorie.physik.uni-goettingen.de

Abstract

**We present an algorithmic construction scheme for matrix-product-operator (MPO) representations of arbitrary $U(1)$-invariant operators whenever there is an expression of the local structure in terms of a finite-states machine (FSM). Given a set of local operators as building blocks, the method automatizes two major steps when constructing a $U(1)$-invariant MPO representation: (i) the bookkeeping of auxiliary bond-index shifts arising from the application of operators changing the local quantum numbers and (ii) the appearance of phase factors due to particular commutation rules. The automatization is achieved by post-processing the operator strings generated by the FSM. Consequently, MPO representations of various types of $U(1)$-invariant operators can be constructed generically in MPS algorithms reducing the necessity of expensive MPO arithmetics. This is demonstrated by generating arbitrary products of operators in terms of FSM, from which we obtain exact MPO representations for the variance of the Hamiltonian of a $S=1$ Heisenberg chain. **

Contents

1. 1 Introduction

2. 2 $U(1)$-invariant tensor networks and application of two-site gates

3. 3 $U(1)$-invariant MPO representation from FSMs

4. 3.1 MPO construction from FSMs

5. 3.2 Maximally branched representation and local transformations on graphs

6. 4 Graph arithmetics and MPO representation of the variance of operators

7. 5 Numerical behavior of the variance for antiferromagnetic $S=1$ Heisenberg chains

8. 6 Conclusions

9. A Construction of a graph representation for lattice-ordered string operator products

1 Introduction

The density-matrix renormalization group (DMRG) [1, 2, 3, 4, 5] has proven to be one of the most powerful tools to treat low-dimensional problems in strongly correlated quantum systems. Inspired by quantum information theory, a formulation in terms of matrix-product states (MPS) [6] and matrix-product operators (MPO) [5] boosted the development of related algorithms (e.g., [7, 8, 9, 10]). Numerous different techniques have been developed in the framework of MPS, most of which require the representation of the Hamiltonian of the system and of further observables in terms of MPOs [11, 12, 13]. Even though finding such a representation is a well-known and generally solved problem [14, 15, 16], it can become arbitrarily complicated as systems may incorporate complex lattice geometries, site-dependent interaction strengths, or transformations of local operators due to their commutation relations.

It is possible to build generic construction schemes for MPO representations of operators by means of finite-state machines (FSM) [15]. However, the numerical realization of these FSMs can be quite involved, especially when exploiting global symmetries [14, 17, 18].

In this paper, we use the fact that FSMs have an underlying graph structure to obtain a generic algorithmic construction scheme for MPO representations of operators conserving $U(1)$ symmetries. Consequently, most of the complexity can be unwrapped into tensor-network manipulations of strings of local operators, which can be mapped one-to-one to graph representations. The obtained construction scheme can be automatized so that an implementation is capable to efficiently create MPO representations. In particular, this also allows us to include commutation rules as well as conservation laws. We demonstrate how to map the operator arithmetics to operations on graph representations, leading to an algorithm for computing the sum as well as the product of two arbitrary operators. We use this algorithm to demonstrate how to compute the variance of a Hamiltonian $\mathop{var}\left(\hat{H}\right)=\left\langle\left(\hat{H}-\left\langle\hat{H}\right\rangle\right)^{2}\right\rangle$ with very high accuracy. This is achieved by reducing the effects of catastrophic cancellation [19] and the total amount of required floating-point operations by using exact operator arithmetics in terms of FSM graphs. In addition, this also minimizes the overall computational costs.

2 $U(1)$-invariant tensor networks and application of two-site gates

Consider an operator $\hat{O}$ acting on a Hilbert space $\mathcal{H}$, which decomposes into a tensor product of $L$-identical $d$-dimensional local Hilbert spaces $\mathcal{H}_{d}$ ($d,\,L\in\mathbb{N}$),

[TABLE]

and assume that this operator is invariant under a global $U(1)$ symmetry generated by local observables $\hat{n}_{i},i\in\left\{1,\ldots,L\right\}$111Note that by $\hat{n}_{i}$ we denote the local generators of the symmetry and not the local density operators, which play this role only in the case of conservation of the total particle number.,

[TABLE]

A state $\left|\psi\right\rangle\in\mathcal{H}$ can be expanded into an MPS with local basis states $\left\{\left|n_{i}\right\rangle\right\}$ spanned by the generators $\left\{\hat{n}_{i}\right\}$ of the global $U(1)$ symmetry

[TABLE]

with $T^{n_{i}}\equiv T^{n_{i}}_{\alpha_{i-1}\alpha_{i}}$ being rank-$3$ tensors acting on site $i$, which carry a physical bond index $n_{i}$ and two virtual bond indices $\alpha_{i-1},\alpha_{i}$. The coefficients $\left\langle n_{1}\cdots n_{L}\right.\left|\psi\right\rangle$ can be reobtained by contracting the corresponding tensors $T^{n_{1}}\cdots T^{n_{L}}$ over all their virtual indices. In a graphical notation, a tensor $T$ is represented by drawing a shape (circles, squares, …) with as many legs attached to it as there are tensor indices (see fig. 1). Contractions over indices are denoted by connecting the corresponding edges. This yields a tensor network.

Due to conservation of the global quantum number $\hat{N}$, these rank-$3$ site tensors $T^{n_{i}}$ are irreducible representations of the global $U(1)$ symmetry. The physical indices $n_{i}$ label the basis states of the local generators. As a consequence, the site tensors can be decomposed into blocks that are subject to an on-site conservation law [17, 18]

[TABLE]

with the block indices $(n_{i},\alpha_{i-1},\alpha_{i})$ labeling the irreducible representations of the global quantum number $\hat{N}$ on the particular site $i$. These tensor blocks in general are complex matrices $T^{n_{i}}_{\alpha_{i-1},\alpha_{i}}\in\mathbb{C}^{d_{\alpha_{i-1}}\times d_{\alpha_{i}}}$ of dimensions $d_{\alpha_{i-1}}\times d_{\alpha_{i}}$.222For brevity, in the following, we suppress bond indices whenever it is clear, which matrix contractions have to be performed, or if the contractions themselves are irrelevant for the discussion. Hence, thinking in terms of block-diagonal MPS, each site matrix $M^{n_{i}}$ is decomposed into the block diagonal form $M^{n_{i}}=\bigoplus_{a}T^{n_{i}}_{a}$ with the non-vanishing blocks $T^{n_{i}}_{a}=T^{n_{i}}_{\alpha_{i-1},\alpha_{i}}$ having block dimensions $d_{\alpha_{i-1}}\times d_{\alpha_{i}}$ and $a$ labeling the irreducible representations at site $i$.

There is a simple graphical representation for these local conservation laws acting on the tensor blocks, which is derived from the tensor-network framework [17]. Whenever there are indices with a (hidden) plus sign in the $\delta$ function, the corresponding bonds in the network are labeled by an ingoing arrow, whereas, for indices with a minus sign, an outgoing arrow is placed on the respective bond (see fig. 1).

Next, we take a more elaborate look at the action of a two-site gate on the symmetry conserving state, ignoring the common framework of irreducible representations and block diagonal structures for a moment (later we generalize the considerations to arbitrary operator expressions). Introducing operators $\hat{A}^{(i)},\hat{B}^{(j)}:\mathcal{H}^{(i,j)}_{d}\longrightarrow\mathcal{H}^{(i,j)}_{d}$ acting only on the sites $i,j$, we start from a generic $U(1)$-invariant expression of the form

[TABLE]

where the $\delta$ function ensures that the total quantum number $\hat{N}$ is conserved when applying the operator. A dummy index $\tau$ with $-(d-1)\leq\tau\leq(d-1)$ is introduced to factorize the $\delta$ function. Thus, suppressing the sums over the physical indices $n_{i}$ for a moment, we obtain

[TABLE]

which nearly restores the factorized shape of generic MPS, even though the physical sites $i,j$ are still connected by the sum over the dummy index $\tau$. It is desirable to restore a tensor-network form so that the action of an operator pair can be written as contraction of tensors acting on the local Hilbert spaces. Therefore, we need to take a closer look at the matrix elements $A^{n_{i}n^{\prime}_{i}}_{\gamma_{i-1}\gamma_{i}},B^{n_{j}n^{\prime}_{j}}_{\gamma_{j-1}\gamma_{j}}$ of the local operators $\hat{A}^{(i)},\hat{B}^{(j)}$. If the total operator expression is $U(1)$-invariant, then the local operators themselves have to be one-dimensional representations in the physical indices (acting on the local Hilbert space) so that each operator carries a unique mapping between states $\left|n_{i}\right\rangle\rightarrow\left|n^{\prime}_{i}\right\rangle$. In other words, for $U(1)$-invariant operator pairs each local operator is non-vanishing only for a certain value of the dummy index $\tau^{\prime}=\Delta$. For instance, in case of spin-ladder operators $\left(\hat{S}^{\pm}\right)^{(i)}$ the total value of $S^{z}$ is locally changed by $\pm 1$. Hence,

[TABLE]

where we introduced the change of local quantum numbers $\Delta$. Having this in mind, block-index conservation laws for the local operators can be realized via

[TABLE]

The case of two-site gates is obtained by setting $\tau_{i-1}=\tau_{j}\equiv 0$ as well as $\tau_{i}=\tau_{j-1}\equiv\tau$; the non-vanishing operator blocks $A^{n_{i}n^{\prime}_{i}}_{\tau_{i-1}\tau_{i}},B^{n_{j}n^{\prime}_{j}}_{\tau_{j-1}\tau_{j}}\in\mathbb{C}^{d_{\tau_{i-1,j-1}}\times d_{\tau_{i,j}}}$ now contain the reduced operator matrix elements acting on the local basis states $\left|n^{\prime}_{i,j}\right\rangle$. The contraction of the block index $\tau$ over the intermediate site tensors $T^{n_{k}},i<k<j$ can be recast into a matrix contraction using the matrix identity

[TABLE]

Therefore, the application of the above two-site gate is factorized completely if we define intermediate shift tensors that act on the sites $k$ with $i<k<j$ as identity, but are contracted over auxiliary bond indices $\tau_{k}$,

[TABLE]

Hence, the action of the $U(1)$-invariant operator pair can conveniently be written as tensor network,

[TABLE]

with the definitions

[TABLE]

When further local operators to the left and right of $\hat{A}^{(i)},\hat{B}^{(i)}$ are absent, the auxiliary indices are shrunk to dummy indices $\tau_{i-1}=\tau_{j}\equiv 0$. Thus, the sum over all remaining auxiliary block indices $\tau_{i\leq l\leq j}$ restores the global conservation law. Note that we have implicitly fixed a gauge freedom carried by the auxiliary indices $\tau_{i}$ to $\tau_{0}=\tau_{L}\equiv 0$ (we can even go further and permit any global change $\delta$ of the overall quantum number $\hat{N}$ by setting $\tau_{L}=\delta$).

Even though these expressions look a bit tedious, they can be represented compactly in form of a tensor network, which also reveals how useful this decoupling turns out to be (see fig. 2(a)). For example, it is possible to identify the well-known transformation law for $U(1)$-invariant MPO site tensors in the expressions above [17], yet the block indices $\tau_{k}$ are related to the change of the quantum number, captured by the shift tensors $R_{\Delta}^{(k)}$.

Generally, contractions over physical bond indices $n_{k},n^{\prime}_{k}$ of such shift tensors correspond to a mapping from one block $a$ of an irreducible representation $U^{(k)}(N)=\bigoplus_{a}U_{a}^{(k)}(N)$ on a site $k$ to another block $a+\Delta$. Hence, given a decomposition into the different quantum-number sectors of an $U(1)$-invariant MPS or MPO tensor

[TABLE]

the contraction along a chain of $l$ shift tensors over physical indices (with $\Delta_{1}\ldots\Delta_{l}$ the changes in the local quantum numbers) is given via

[TABLE]

From this point on it is clear how to generalize the considerations to arbitrary strings of local operators. In a given expression of local operators, we need to identify pairs of operators that conserve the global $U(1)$ quantum number but locally change the on-site quantum numbers $n^{\prime}_{i_{1}}\rightarrow n_{i_{1}},n^{\prime}_{i_{2}}\rightarrow n_{i_{2}}$. For each pair, shift tensors then have to be inserted acting on the intermediate sites $i_{1}<k<i_{2}$. Note that in a code there is no need to explicitly implement the shift tensors. It suffices to implement only their action on the virtual bonds of either the MPS or of another MPO, i.e., to collect all vertical strings of shift tensors $\left\{\hat{R}^{(k)}_{\Delta_{r}}\right\}_{1\leq r\leq l}$ in the network and to apply the shifting in the auxiliary block indices $\alpha_{k-1,k}\rightarrow\tilde{\alpha}_{k-1,k}=\alpha_{k-1,k}+\Delta_{1}+\cdots+\Delta_{l}$ (see fig. 2(b)).

3 $U(1)$-invariant MPO representation from FSMs

As discussed in [15, 16], the MPO representation of an operator on a tensor-product space $\mathcal{H}=\mathcal{H}_{d}^{\otimes L}$ can be obtained from the transition amplitudes of FSMs. In the following, the underlying graph structure of FSMs is used extensively. Thus, we first give a brief review on how to identify FSMs with MPO representations, following [15].

3.1 MPO construction from FSMs

Let $\mathcal{K}=\left\{\hat{o}^{(i)}_{1},\cdots,\hat{o}^{(i)}_{m}\right\}$ be a set of $m\in\mathbb{N}$ local operators $\hat{o}_{i}:\mathcal{H}_{d}\rightarrow\mathcal{H}_{d}$. Any global operator is of the general form

[TABLE]

with $\hat{h}^{(i)}_{\nu,r}=f^{r}_{\nu_{1}\ldots\nu_{n}}\hat{o}^{(i)}_{\nu_{1}}\cdots\hat{o}^{(i+r)}_{\nu_{n}}$ being strings of $n$ local operators $\hat{o}_{\nu}\in\mathcal{K}$ coupling lattice sites $i$ with range $r+1\geq n$, amplitude $f^{r}_{\nu_{1}\ldots\nu_{n}}\in\mathbb{C}$, and abbreviating the index set $\nu=(\nu_{1}\ldots\nu_{n})$. For later convenience we will call $\hat{h}^{(i)}_{\nu,r}$ lattice ordered $n$-point $r+1$-ranged operator strings.

The set of all lattice-ordered $n$-point $r+1$-ranged operator strings $\Sigma=\left\{\hat{h}^{(i)}_{\nu,r}\right\}_{\nu,r}$ defines a language of a FSM, i.e., there is a set of states $\mathcal{L}$ and a transition function $\delta:\mathcal{L}\times\mathcal{K}\longrightarrow\mathcal{L}$ so that with a proper initial and final state $I,F\in\mathcal{L}$, the FSM $M(\mathcal{K},\mathcal{L},\delta,I,F)$ is obtained with a generated language $\Sigma$. The corresponding graph is then a representation of $\hat{H}$ and is denoted by $\Lambda(\hat{H})$.

An example for the graph representation of the XX model is given in fig. 3(a) and in the following, we shortly explain how to obtain this graph.

At first, we have to rewrite the global operator in terms of lattice-ordered $n$-point $r+1$-ranged operator strings333In this example we change our convention and write site indices as lower indices.

[TABLE]

i.e., we have two distinct $2$-point $2$-ranged operator strings. Then, we define a default set of states $\mathcal{L}_{0}=\left\{I,F\right\}$ with $I,F$ being the initial and final state of the FSM, respectively. For convenience, we will highlight them in the corresponding graphs with green ($I$) and red ($F$) background. In general, FSMs are capable to generate sequences of symbols $\hat{o}\in\mathcal{K}$ by transitioning between states $a\in\mathcal{L}$ via permitted transitions, i.e., those with non-vanishing transition function $\delta(\hat{o},a)$. Thus, we define the initial/final state by demanding that all the previous/subsequent transitions have to be identities.

Next, we build the graph fig. 3(a) by constructing each operator string separately, so we may start with $\hat{h}^{+-}_{(i)}$. We therefore add $r=1$ intermediate states (here only one state $A_{1}\equiv A$) to the set of states: $\mathcal{L}=\mathcal{L}_{0}\cup\left\{A\right\}$ and define permitted transitions $\delta(I,\hat{S}^{+})=A$ and $\delta(A,\hat{S}^{-})=F$. Thus, the FSM can transition from the initial state (after placing an arbitrary number of identities) into state $A$ by appending an operator $\hat{S}^{+}$ onto the so far constructed operator string. Being in state $A$, the FSM has no other choice than to transition into state $F$ by appending a subsequent operator $\hat{S}^{-}$. A generalization of the construction of such local operator strings is shown in fig. 3(b) for the case of a general $2$-point $r$-ranged operator string.

From this point on, we will call graphs that generate only one distinct type of $n$-point $r+1$-ranged operator string single-branched graphs and identify them with the corresponding contribution in the global operator $\hat{H}_{\nu,r}$.

Returning to the XX model, we can finish the construction of the FSM by adding another single-branched graph transitioning from the initial state into an additional state $B$ via a local operator $\hat{S}^{-}$ and a transition into the final state via a local operator $\hat{S}^{+}$.

Finally, we can construct the MPO representation, i.e., the operator-valued matrices $\hat{W}^{n_{i}n^{\prime}_{i}}$, by assigning matrices $W^{n_{i}n^{\prime}_{i}}_{\tau_{i-1}\tau_{i}}$ of dimension $\left|\mathcal{L}\right|\times\left|\mathcal{L}\right|$ for all non-vanishing symmetry blocks $(n_{i},n^{\prime}_{i},\tau_{i-1},\tau_{i})$ (i.e., those with at least one local operator block $\hat{o}^{n_{i}n^{\prime}_{i}}_{\tau_{i-1}\tau_{i}}\neq 0$). Figure 4(a) sketches a general FSM and fig. 4(b) shows the generated MPO bulk tensor, in which the rows and columns are labeled by the states of the FSM. The corresponding boundary tensors are obtained by projecting out (a) the transition from the initial state into the bulk for $i=1$ and (b) the transitions from the bulk into the final state for $i=L$. We emphasize that the site-dependent coefficients $c^{(i)}_{ab}\in\mathbb{C}$ are free parameters and therefore can be chosen independently for every site.

3.2 Maximally branched representation and local transformations on graphs

As already mentioned, every global operator on $\mathcal{H}$ can be formulated as a sum over lattice-ordered $n$-point $r+1$-ranged operator strings

[TABLE]

Note that the graph representation via FSM is not unique; for every operator $\hat{H}$, there is a set of corresponding FSMs $\left\{\Lambda(\hat{H})\right\}_{\Lambda}$. Therefore, we are free to choose one representation $\Lambda(\hat{H})$, which makes it easier to perform operator arithmetics and then switch to another representation $\tilde{\Lambda}(\hat{H})$ to find the most compact MPO.

Referring to eq. 27, a natural translation to the graph representations of $\hat{H}$ can be obtained by introducing a commutative map $\oplus$ between graph representations of sums of operators $\hat{H}_{1},\hat{H}_{2}$ via

[TABLE]

The realization of $\oplus$ in terms of graphs is obtained by taking the graph representations of the operators $\hat{H}_{1},\hat{H}_{2}$ and by merging the initial and final states as depicted in fig. 5.

Next, we can define the notion of a maximally branched graph representation, which is given by the graph $\Lambda_{\rm max}(\hat{H})$ satisfying the conditions: a) the initial state $I$ is the only state with more than one child state, b) the final state $F$ is the only state with more than one parent state. $\Lambda_{\rm max}(\hat{H})$ satisfies the equation

[TABLE]

This representation has several advantages; most importantly for our discussion, any local transformation of the operator can be mapped one-to-one to the transitions of the graph representation. Here local transformations are those transformations that map a lattice-ordered $n$-point $r+1$-ranged operator string into another lattice-ordered $m$-point $r+1$-ranged operator string $\hat{H}_{\nu,r}\rightarrow\hat{\tilde{H}}_{\tilde{\nu},r}$ without changing $r$. To clarify what is meant by this mapping we emphasize that each branch $\hat{H}_{\nu,r}$ generates exactly one type of local-operator string $\hat{o}_{\nu_{1}}\ldots\hat{o}_{\nu_{r}}$ and therefore the tensor representation of these strings factorizes on the local Hilbert spaces. Hence, if we can give a factorization of the local transformation $\hat{U}$ in terms of tensors acting on the local Hilbert spaces (e.g., $\hat{u}_{\nu_{1}}\ldots\hat{u}_{\nu_{r}}$), then we can represent the transformation directly by contracting the transformation tensors with the local operators over their physical indices

[TABLE]

Note that the transformations in eq. 19 forcing conservation of $U(1)$ quantum numbers for two-site gates are of exactly this kind and so is their generalization to arbitrary strings of local operators. Thus, conservation of $U(1)$ quantum numbers can be implemented by transforming the graph representation of an operator into its maximally branched representation. The local operator strings in each branch are transformed by applying shift tensors $\hat{R}_{\Delta}$. For example, let us consider the transformation for a next-to-nearest-neighbor spin-flip term

[TABLE]

with $\hat{S}^{\pm}$ being the usual angular-momentum ladder operators with lower site indices. The transformed graph representation is given in fig. 6(a).

Conveniently, these rules can be extended to anticommuting operators by applying a Jordan-Wigner transformation444Note that this argument also holds in higher dimensions, because in MPS approaches a 1D path is used to sweep through the system, and the Jordan-Wigner transformation is also applicable in the presence of long-range interactions.. For $U(1)$-invariant operators, the Jordan-Wigner transformation is also local, as there has to be a corresponding annihilation operator for every fermionic creation operator appearing in a string operator – and vice versa – because of quantum-number conservation. The Jordan-Wigner transformation can be implemented as a product of parity operators via

[TABLE]

with $\hat{a}_{(i)}^{(\dagger)}$ annihilation (creation) operators for hard-core bosons at site $i$. Then we find that, for any $U(1)$-conserving product of fermionic creation and annihilation operators, the transformations act only within the operator strings. For instance, for a next-to-nearest-neighbor-hopping term we find

[TABLE]

Again, these transformations have a simple graph representation, see fig. 6(b).

To sum up, we have derived a construction scheme for MPO representations of generic $U(1)$-invariant operators that takes a FSM as input. Specifying the phase factor $e^{i\phi}$ for the commutation relations of the local operators (e.g., $\phi=2\pi$ for bosons, $\pi$ for fermions), the scheme automatizes the construction of the MPO site tensors so that we can identify the graph representation of the FSM with the MPO. This permits us to take advantage of the graph representation to improve operator arithmetics, which will be discussed in the following section.

4 Graph arithmetics and MPO representation of the variance of operators

Having derived a construction scheme that permits us to automatize the generation of MPO representations for $U(1)$-invariant operators, we now make use of the established connection between graphs and operators by replacing operator arithmetics with graph manipulations. We then demonstrate the power of this approach by employing the graph algebra to derive various expressions for the variance of a Hamiltonian $\hat{H}$. As we will show, this not only allows for more efficient calculations, but also addresses the problem of catastrophic cancellation that comes up in a naïve evaluation of $\mathop{var}(\hat{H})=\left\langle\hat{H}^{2}\right\rangle-\left\langle\hat{H}\right\rangle^{2}$ due to the need to subtract large numbers, which scale roughly as $\mathcal{O}(L^{2})$ in the system size $L$ [19].

Let us consider the product of two global operators $\hat{H}_{1},\hat{H}_{2}$ in terms of their maximally branched representations

[TABLE]

and in particular a single summand that is the product of two lattice-ordered string operators555We use the notation introduced in eq. 25.

[TABLE]

Note that we have introduced superscripts $\nu^{1,2},r^{1,2}$ to distinguish the index sets of the global operator666Expanding the indices one would have $\hat{H}_{i}=\sum_{\nu^{i}_{1}\ldots\nu^{i}_{n^{i}},r^{i}}\hat{H}_{\nu^{i}_{1}\ldots\nu^{i}_{n^{i}}r^{i}}$. A representation of this product in terms of a FSM and therefore its graph representation requires a reformulation in terms of lattice-ordered string operators. Although the product of two lattice-ordered operators $\hat{H}_{r^{1},\nu^{1}}\cdot\hat{H}_{\nu^{2},r^{2}}$ is no longer lattice-ordered, a careful inspection of the terms violating the lattice order reveals how to build a graph representation generating the product $\hat{H}_{r^{1},\nu^{1}}\cdot\hat{H}_{r^{2},\nu^{2}}$. It turns out to be useful to define a non-commutative $\wedge$ product, which maps two single-branched graphs to a single-branched graph via

[TABLE]

A graph realization of $\wedge$ is obtained by identifying the final state of $\Lambda(\hat{H}_{\nu^{1},r^{1}})$ with the initial state of $\Lambda(\hat{H}_{\nu^{2},r^{2}})$. For instance, see fig. 7(a) for an exemplary evaluation of $\Lambda(\sum_{i}\hat{o}^{(i)}_{11}\hat{o}^{(i+1)}_{12})\wedge\Lambda(\sum_{j}\hat{o}^{(j)}_{21}\hat{o}^{(j+1)}_{22})$.

As carried out in appendix A, an algorithm can be constructed that yields the following graph representation

[TABLE]

with the symmetrized wedge product

[TABLE]

and $\operatorname{sgn}(\hat{h}_{\nu^{1},r^{1}},\hat{h}_{\nu^{2},r^{2}})$ the sign of the commutation relation between local operators $\hat{h}_{\nu^{1},r^{1}}$ and $\hat{h}_{\nu^{2},r^{2}}$ acting on different sites. Employing linearity of the graph representation for addition, the product of two general operators can then be formulated via

[TABLE]

Despite the compact form, we emphasize that eq. 37 and eq. 39 describe a graph representation that is maximally expanded, so that there is no branching below the initial node, and the bond dimension of the generated MPO is very large. However, the size of the graph can be reduced very efficiently by shrinking it into its most compact form. The idea behind the shrinking is best demonstrated with a concrete example. Consider an expanded graph generated from merging two branches

[TABLE]

The number of nodes can be reduced by fusing edges sharing one node and carrying the same transitions. For example, in fig. 7(b) these are the edges connecting the states $I\rightarrow A_{1}$ and $I\rightarrow A_{2}$. In the same way, entire branches that are completely equal with respect to their transitions can be fused together by employing the linearity of $\oplus$

[TABLE]

But we can go even further by defining which local operator pairs acting on the same lattice site vanish identically. For example, this is the case for $\hat{S}^{+}\hat{S}^{+}\equiv 0$ for $S=1/2$ models. Generated branches that contain this type of transitions can be discarded completely. For more complex graphs, these shrinking procedures can be applied iteratively, so that only the reduced graph is stored and used in the actual calculations.

An illustrative example demonstrating the advantages of the approach above is to construct two different expressions for the variance of the Hamiltonian $\hat{H}$ of a system, where the goal is to avoid catastrophic cancellation as best as possible while keeping calculations as cheap as possible. For instance, the variance can be used as a control parameter in numerical simulations to test whether a state $\left|\psi\right\rangle$ is close to an eigenstate of $\hat{H}$. A naïve evaluation is obtained by directly calculating the expectation values in

[TABLE]

which eventually vanishes, if $\left|\psi\right\rangle$ is an exact eigenstate. However, as the expectation value of $\hat{H}^{2}$ scales as $L^{2}$, explicit evaluation of eq. 42 has major drawbacks when it comes to numerical calculations with finite-precision arithmetics such as catastrophic cancellation. Every (exact) number $z$ is numerically represented up to a certain precision [20], which is usually measured in orders $p$ of magnitudes so that we can mimic limited numerical precision by replacing $z\rightarrow z(1+\epsilon\times 10^{-p})$ with a random variable $\epsilon\in\left(-1,1\right)$. Let $z_{1,2}$ be numbers represented with the same numerical precision $p$ and $0<z_{1}-z_{2}=10^{-\delta}$ their exact difference. In finite-precision arithmetics, we then obtain

[TABLE]

For $\delta>0$, the second term cannot be represented due to finite precision. In our case, the values for $z_{1,2}$ are obtained from expectation values of operators acting on the whole system. Hence, if we estimate them by their leading-order contribution $z_{1,2}\sim L^{q}$ with magnitude $q$, with $\gamma\equiv\log_{10}(L)$, and $\epsilon\equiv\epsilon_{2}-\epsilon_{1}$, we obtain

[TABLE]

where $\epsilon$ is a random variable of order $\pm 1$. It follows that we need $\delta<p-\gamma\cdot q$ in order to have a reasonable numerical outcome. In case of the variance, i.e., $q=2$, with double-precision arithmetics, $p=p_{\text{num}}=16$, the naïve evaluation, $\delta<16-2\gamma$, yields an upper bound of a maximally possible precision of $10^{-12}$ for a lattice with $L=100$ sites. However, the numerical precision in general is not only bound by the exact numerical precision $p_{\text{num}}$, but it is also subject to round-off errors of preceding calculations. Thus, in actual calculations, we have an effective precision that depends on simulation parameters $p(L,\chi,...)\leq p_{\text{num}}$. Returning to eq. 42, we realize that the calculated variance is not only bound, but may even become negative (as $\epsilon$ in eq. 44 can be negative).

The problem can be addressed by a minimization of $\gamma\cdot q$, for instance by constructing an MPO representation for the operator $\hat{V}=(\hat{H}-\left\langle\psi\right|\hat{H}\left|\psi\right\rangle)^{2}$, (e.g., by distributing the expectation value $\left\langle\psi\right|\hat{H}\left|\psi\right\rangle$ over the lattice sites). Unfortunately, this comes at the cost of having constructed an operator that is dependent on its expectation value. Hence, its MPO representation has to be rebuilt for every new state $\left|\psi\right\rangle$, which requires an efficient way of obtaining MPO representations of powers of operators while keeping the numerical effort low.

To analyze both problems, we investigate two MPO representations $\Lambda_{1}(\mathop{var}(\hat{H}))$ (see fig. 8(a) for the corresponding graph) and $\Lambda_{2}(\mathop{var}(\hat{H}))$, which are obtained from the graph representations

[TABLE]

Let us briefly discuss the properties of the two graph representations and their generated MPO. We start with the numerical implementation and its costs. Both representations are obtained by loading pre-computed graph products at runtime. Then, we explicitly calculate $\langle\hat{H}\rangle$ with numerical costs scaling as $\mathcal{O}(LdD^{2}D^{2}_{W})$ with $D$ and $D_{W}$ being the maximal matrix dimensions of the MPS and MPO site tensors, respectively. Construction of the MPO representation for $\Lambda_{1,2}(\mathop{var}(\hat{H}))$ then only requires allocation of the site tensors with costs scaling as $\mathcal{O}(L[dD_{W_{1,2}}]^{2})$. Here, $D_{W_{1,2}}$ denotes the maximal bond dimensions of the MPO site-tensor matrices generated from the graphs $\Lambda_{1,2}(\mathop{var}(\hat{H}))$. Subsequently, the MPO tensors of the graphs are compressed by using an SVD with numerical costs scaling as $\mathcal{O}(L\cdot[d^{2}\cdot D_{W_{1,2}}]^{3})$.

To sum up, the numerical costs for obtaining the MPO representation of the graphs $\Lambda_{1,2}(\mathop{var}(\hat{H}))$ are to leading order governed by the expenses when calculating the expectation value $\langle\hat{H}\rangle$, as long as the MPS bond dimensions are the cost-determining factor. This demonstrates a major advantage of mapping the MPO arithmetics onto the graph representations: As the expectation value $\langle\hat{H}\rangle$, in general, is already available, the additional numerical costs for constructing the MPO representation are not significant.

Next, we take a look at the numerical stability. The graph $\Lambda_{2}(\mathop{var}(\hat{H}))$ generates the MPO representation of the operator $\langle\hat{H}-\langle\hat{H}\rangle\rangle$ multiplied by itself. Thus, the graph of $\hat{H}$ is expanded so that the expectation value $\langle\hat{H}\rangle$ is equally distributed over all lattice sites with an on-site value of $\langle\hat{H}\rangle/L$. As operator arithmetics are represented exactly by means of the constructed graph $\Lambda_{2}(\mathop{var}(\hat{H}))$, the only relevant source of catastrophic cancellation is the evaluation of $\langle\hat{H}-\langle\hat{H}\rangle\rangle$ along the lattice that compares terms of order $\mathcal{O}(L)$, hence $q=1$. Thus, we expect the variance to be bound from below by $10^{16-\gamma}$. Yet, the graph $\Lambda_{1}(\mathop{var}(\hat{H}))$ in general also suffers from catastrophic cancellation with $q=2$, as in the naïve evaluation of the variance. This is best seen by evaluating the structure of the generated matrix representation for non-vanishing tensor blocks $W^{n_{i}n^{\prime}_{i}}_{\tau_{i-1}\tau_{i}}$ (see fig. 8(b)). As the latter contains the complete matrix representation of $\hat{H}^{2}$, we end up at the final site by comparing numbers of order $\mathcal{O}(L^{2})$ when performing the tensor contractions to evaluate the variance. Therefore, $q=2$ yields the variance to be bound from below by $10^{16-2\gamma}$.

In addition, graph arithmetics are exact, whereas MPO arithmetics need a vast amount of numerical operations completed before expectation values can be calculated. Therefore, using MPO arithmetics is much more prone to collecting round-off errors, which, in drastic cases, reduces the numerical precision $p$ by orders of magnitudes [5].

5 Numerical behavior of the variance for antiferromagnetic $S=1$ Heisenberg chains

In this section, we test the behavior of the variance obtained from site-tensor representations of $\Lambda_{1,2}(\mathop{var}(\hat{H}_{\rm Heisenberg}))$ for $S=1$ antiferromagnetic Heisenberg chains [21, 22]

[TABLE]

with $\mathbf{\hat{S}}_{i}$ the vector of $S=1$ spin operators on lattice site $i$. For this purpose, we exploit the total magnetization of a system with $L$ lattice sites as the conserved $U(1)$ quantity. For the ground state search, we sweeps through the system (using open boundary conditions) and optimize the site tensors via a standard Lanczos algorithm, see [23, 24, 25]. All MPS contractions are formulated in two-site representation [5].

The largest observed bond dimensions $D_{W_{1,2}}$ of the site-tensor matrices $W^{n_{i}n^{\prime}_{i}}_{\tau_{i-1}\tau_{i}}$ for the MPO representations of $\Lambda_{1,2}(\mathop{var}(\hat{H}_{\rm Heisenberg}))$ are $D_{W_{1},max}=12$ and $D_{W_{2},max}=10$, respectively. Note that these bond dimensions are generally smaller than MPS bond dimensions used during the simulation. We conclude from $D_{W_{2},max}<D_{W_{1},max}$ that the distribution of the constant energy term $\langle\hat{H}_{\rm Heisenberg}\rangle$ over the lattice sites in $\Lambda_{2}(\mathop{var}(\hat{H}_{\rm Heisenberg}))$ ensures a more efficient compression of the resulting MPO representation. Consistently, for common MPS bond dimensions, we find that $\Lambda_{2}(\mathop{var}(\hat{H}_{\rm Heisenberg}))$ is evaluated faster than $\Lambda_{1}(\mathop{var}(\hat{H}_{\rm Heisenberg}))$. In addition, fig. 9 displays the build time of the MPO representation from the graphs $\Lambda_{1,2}(\mathop{var}(\hat{H}_{\rm Heisenberg}))$ for various system sizes. We find a perfect linear scaling $\sim\mathcal{O}(L)$ for the dependency of the build time, which moreover has no impact on the overall computation time.

We have performed simulations, in which we varied either the total bond dimension per MPS site tensor $\chi_{\text{max}}$ or the discarded weight $w$ while keeping the respective other fixed.777We define the total bond dimension $\chi_{\text{max}}$ as the number of singular values kept per site and the discarded weight as the sum over all squares of neglected singular values $w=\sum^{D}_{k=\chi_{\text{max}}+1}S^{2}_{k}$. Figures 10(a) and 10(b) show results for various values of $\chi_{\text{max}}$ and $w$. An important criterion for consistency of the calculations is the independence of the threshold at which catastrophic cancellation sets in. We find that varying both parameters $\chi_{\text{max}}$ and $w$ yields a constant value of

[TABLE]

which corresponds to the value at which the graph representation $\Lambda_{1}(\mathop{var}(\hat{H}_{\rm Heisenberg}))$ saturates. Employing eq. 44, these results suggest for the actual numerical precision an ansatz of the form $p=p_{\text{num}}-p_{r}(L)$ with a correction $p_{r}(L)$, which to first order depends only on the lattice size $L$. Repeating the calculations for $S=1$ Heisenberg chains with various system sizes, we can thus extract $p_{r}$ from eq. 44 by estimating $\delta=\log_{10}(\mathop{var}(\hat{H}_{\rm Heisenberg}))$ from the saturated value for the variances. For the two graph representations we consider the estimator for the actual numerical precision

[TABLE]

For both graphs, we perform a linear fit obtaining

[TABLE]

which is shown in Fig. 11. We emphasize that the observed behavior is perfectly consistent with the ansatz above for constant double-precision arithmetics $p_{\text{num}}=16$ and residual numeric precision $p_{r}(L)\approx\gamma(L)=\log_{10}(L)$. We hence find that, aside from catastrophic cancellation, the dominating contribution to the loss of numerical precision is proportional to the lattice size, which can be associated with inevitable rounding errors generating an error $\mathcal{O}(10^{-p_{\text{num}}})$ per lattice site.

To complete this section, we now turn to the method discussed by Hubig et al. in [14], which suggests a complementary access for introducing operator arithmetics for MPOs. In short, Hubig et al. present a scheme to build MPO representations numerically via in-code evaluation of direct sums and Kronecker products of local operator strings acting on individual lattice sites. To account for the growth in the MPO bond dimension, they discuss various numerical compression schemes and benchmark the obtained MPOs for a $100$-site $S=1$ Heisenberg chain. This scheme is somewhat more straightforward, as it only involves elementary matrix operations.

However, the FSM approach has additional advantages when it comes to numerical calculations. In particular, when comparing our results for the variance with those obtained by Hubig et al. for the same numerical simulation parameters888Note that Hubig et al. employ conservation of non-abelean quantum numbers. Therefore, the MPS bond dimensions used there do not directly compare to the ones used here. However, this does not affect the maximally achievable resolution of $\left\langle\mathop{var}(\hat{H})\right\rangle$ due to finite precision arithmetics., we find that the maximally obtained precision is two orders of magnitudes larger when using the MPOs obtained from graph arithmetics. The reason lies in the above observation that the graph representations only require a minimal number of floating-point operations. In fact, the MPO representation for $\mathop{var}(\hat{H})$ constructed from numerical MPO arithmetics as suggested by Hubig et al. already requires at least $L$ Kronecker products and sums to be evaluated for the advantageous MPO representation of the variance. In contrast, the FSM is generated from abstract graph operations and is therefore an exact representation.

Thus, before the actual variance calculation (namely the contraction of the tensor network representation of $\left\langle\psi\right|\mathop{var}(\hat{H})\left|\psi\right\rangle$) can begin, the MPO representation obtained from numerical MPO arithmetics already picked up numerical round-off errors of magnitude $\mathcal{O}\sim L=10^{\gamma}$. Therefore, in the example of the $100$-site $S=1$ Heisenberg chain, the decrease in precision can be estimated to be $p_{\text{loss}}\sim\gamma=2$ magnitudes, which is exactly what we found in our calculations.

Another nice side effect is that performing graph compression as described in section 4 already covers the deparallelisation compression method which was first introduced in [26].

Finally, note that implementing local transformations on an abstract level as discussed in section 3 maps the whole complexity of index shifting, as required by quantum-number conservation or in an implementation of fermionic anticommutation rules, from the code to an input level. In other words, there is no need for the programmer to hard-code these features when MPOs are generated from FSMs. Instead, they can be incorporated by designing the corresponding graphs, which in turn enormously increases the flexibility of the code.

6 Conclusions

We have formulated an optimized algorithmic construction scheme for efficient MPO representations of $U(1)$-invariant operators generated by FSMs. This scheme allows implementations to automatize the application of local transformations of operators in MPO representations, e.g., while exploiting $U(1)$ symmetries by propagating local changes in quantum numbers or by tackling the fermionic sign via a Jordan-Wigner transformation. As a consequence, graph representations for FSMs can be interpreted directly as representations of operators. Based on this, operator arithmetics are then mapped to transformations of the underlying graphs to generate exact MPO representations of operator sums and products. This permits us to exactly calculate operator arithmetics, which can be stored and quickly loaded in the course of simulations. We demonstrated the effectiveness of our approach by considering two graph representations $\Lambda_{1,2}(\mathop{var}(\hat{H}))$ of the variance of a system’s Hamiltonian $\hat{H}$. Investigating their numerical properties in a ground-state calculation for a $S=1$ Heisenberg chain with $L=100$ lattice sites, both representations behave numerically consistent and stable with a resolution for the variances $\mathop{var}_{1,2}(\hat{H}_{\rm Heisenberg})$ of the graph representations of $\mathop{var}_{1}(\hat{H}_{\rm Heisenberg})\stackrel{{\scriptstyle\sim}}{{=}}\mathcal{O}(10^{-10})$ and $\mathop{var}_{2}(\hat{H}_{\rm Heisenberg})\stackrel{{\scriptstyle\sim}}{{=}}\mathcal{O}(10^{-12})$, respectively. Investigating the dependence of the numerical breakdown on the lattice size shows that the graph representations achieve a numerical precision of $p\stackrel{{\scriptstyle\sim}}{{=}}\mathcal{O}(p_{\text{num}}-\log_{10}(L))$ in the calculations. We conclude that this high numerical precision is due to the exact graph representation of the operator arithmetics and comes without any significant additional computational costs during runtime.

Finally, we note that wrapping operator arithmetics into re-usable graph representations helps to obtain efficient and exact MPO representations, which are useful for various applications. For instance, consider variational problems of the form

[TABLE]

to find highly excited eigenstates and eigenvalues. Solving this minimization problem benefits significantly from the increased numerical precision and gains in efficiency of the calculation of the expectation value of $(\hat{H}-E)^{2}$.

Further applications of the introduced approach cover the efficient representation of long-ranged swap gates, and the concepts can be generalized to higher-dimensional tensor networks (e.g. PEPS) [27]. In this case, the transition amplitudes get additional degrees of freedom (’color’), corresponding to either transversal or longitudinal auxiliary indices. Similar benefits are to be expected as in 1D MPS, allowing for more precise and more flexible implementations of tensor network methods in higher dimensions.

Acknowledgements

We thank M. Marahrens and B. Lenz for helpful discussions, and I. Köhler for carefully proofreading the manuscript. We gratefully acknowledge financial support by the Deutsche Forschungsgemeinschaft (DFG) through Research Unit FOR 1807 (project P7) and SFB/CRC 1073 (project B03).

Appendix A Construction of a graph representation for lattice-ordered string operator products

Consider two lattice-ordered string operators $\hat{H}_{\nu^{1,2},r^{1,2}}$ as introduced in the main text (eqn. 25) with single-branch graph representations $\Lambda(\hat{H}_{\nu^{1},r^{1}})$ and $\Lambda(\hat{H}_{\nu^{2},r^{2}})$, which are parts of maximally branched operators $\hat{H}_{1,2}$. In the following, we present an algorithm to construct the graph representation of the operator product $\hat{\tilde{H}}=\Lambda(\hat{H}_{\nu^{1},r^{1}})\cdot\Lambda(\hat{H}_{\nu^{2},r^{2}})$ in terms of generating a new graph. This procedure can then be applied to all branches to construct a new graph for $\hat{H}=\hat{H}_{1}\cdot\hat{H}_{2}$.

Let $\mathcal{T}(\mathcal{K},b=1,n)$ be the set of all single-branch graphs representing lattice-ordered $n$-point operators. Then, with a proper $\beta\geq b$, we look for a realization of the non-commutative map

[TABLE]

with $\mathcal{K}\times\mathcal{K}^{\prime}$ denoting the symmetrized on-site tensor-product set of $\mathcal{K}$ and $\mathcal{K}^{\prime}$. For this purpose, we apply the definition of $\oplus$ in eq. 28 and search for a graph representation of $\Lambda(\hat{H}_{\nu^{1},r^{1}}\cdot\hat{H}_{\nu^{2},r^{2}})$ by ordering the appearing types of terms in the resulting sum of the operator product according to the lattice treated. We construct single-branch graph representations for all different types of generated lattice-ordered operator strings, which we denote by $\hat{\gamma}$. Then, a graph representation is obtained by summing up all these strings

[TABLE]

From now on, we focus on the special case of $2$-point operators, i.e., operators of the form

[TABLE]

Nevertheless, the generalization to arbitrary $n$-point string operators is straightforward: Simply replace identities with additional local operators. Decomposing the operator product, we find

[TABLE]

Next, the commutation relation of local operators acting on different sites $i_{1}\neq i_{2}$ fulfills

[TABLE]

with $\epsilon_{\nu_{1}\nu_{2}}$ the sign according to the commutation relation. Then, we can decompose the product into lattice-ordered sums by commuting local operators acting on strictly unequal sites.

The first lattice-ordered contribution is given via $\hat{H}_{A}$. The corresponding diagram is a single-branch graph obtained by identifying the final state of the graph $\Lambda(\hat{H}_{\nu^{1},r^{1}})$ with the initial state of the graph $\Lambda(\hat{H}_{\nu^{2},r^{2}})$ by introducing an intermediate state $E$ (see fig. 12(a)). We now make use of the wedge product $\wedge$ for single-branched graphs as defined in eq. 36 to rewrite $\hat{H}_{A}$ in short as

[TABLE]

Swapping the operators, we obtain another lattice-ordered sum by commuting all local operator contributions, so that the corresponding graph picks up two factors, $\epsilon_{\nu_{1}^{1}\nu_{1}^{2}}\epsilon_{\nu_{2}^{1}\nu_{1}^{2}}$ for commuting $\hat{o}^{(i_{2})}_{\nu_{1}^{2}}$ and $\epsilon_{\nu_{1}^{1}\nu_{2}^{2}}\epsilon_{\nu_{2}^{1}\nu_{2}^{2}}$ for commuting $\hat{o}^{(i_{2}+r^{2})}_{\nu_{2}^{2}}$ with $\hat{o}^{(i_{1})}_{\nu_{1}^{1}}\hat{o}^{(i_{1}+r^{1})}_{\nu_{2}^{1}}$ (see fig. 12(b)). Again, the corresponding graph can be expressed via a wedge product,

[TABLE]

where we have introduced $\operatorname{sgn}(\hat{h}_{\nu^{1}_{1}\nu^{1}_{2},r^{1}},\hat{h}_{\nu^{2}_{1}\nu^{2}_{2},r^{2}})\equiv\epsilon_{\nu_{1}^{1}\nu_{1}^{2}}\epsilon_{\nu_{2}^{1}\nu_{1}^{2}}\epsilon_{\nu_{1}^{1}\nu_{2}^{2}}\epsilon_{\nu_{2}^{1}\nu_{2}^{2}}$ for brevity. Note that, for $U(1)$-invariant operators, the identity $\operatorname{sgn}(\hat{h}_{\nu^{1}_{1}\nu^{1}_{2},r^{1}},\hat{h}_{\nu^{2}_{1}\nu^{2}_{2},r^{2}})\equiv 1$ holds.

The remaining sums over $\hat{h}_{B}(i_{1},i_{i_{2}})$ correspond to overlapping interaction terms, i.e., all those lattice indices $i_{1},i_{2}$ that fulfill

[TABLE]

In order to generate all these terms using one algorithm, we write the $2$-point operator expressions graphically by representing a local operator $\hat{o}_{\nu}$ with a cross (“$\times$”) and the intermediate vacant operator sites with a circle (“$\circ$”), e.g.,

[TABLE]

Employing this condensed notation, all lattice-ordered combinations of local operator strings can be generated by placing the graphical representations of $\hat{h}_{\nu^{1}_{1}\nu^{1}_{2},r^{1}}^{(i_{1})}$ and $\hat{h}_{\nu^{2}_{1}\nu^{2}_{2},r^{2}}^{(i_{2})}$ next to each other by aligning the last operator of $\hat{h}_{\nu^{1}_{1}\nu^{1}_{2},r^{1}}^{(i_{1})}$ with the first operator of $\hat{h}_{\nu^{2}_{1}\nu^{2}_{2},r^{2}}^{(i_{2})}$ (see fig. 12(c)). Subsequently, the right string is shifted upward until the initial operator of $\hat{h}_{\nu^{1}_{1}\nu^{1}_{2},r^{1}}^{(i_{1})}$ is aligned with the final operator of $\hat{h}_{\nu^{2}_{1}\nu^{2}_{2},r^{2}}^{(i_{2})}$ (see fig. 12(f)). While the right string’s position is shifted by one step, the local operators $\hat{o}_{\nu_{j}^{1}}$ of the left operator string pick up a sign factor whenever they pass a local operator $\hat{o}_{\nu_{i}^{2}}$ in the right string, which is denoted by adding a factor $\epsilon_{\nu_{i}^{1}\nu_{j}^{2}}$ to the left condensed representation (see figs. 12(d) and 12(e)).

For each such step $\Delta$, with $0\leq\Delta\leq r^{1}+r^{2}$, this algorithm generates a string $\hat{\gamma}_{\Delta}$ of local operators by merging the shifted right string into the left and resubstituting the original operators: $\times\rightarrow\hat{o}_{\nu}$. Missing sites are replaced with identities, whereas two local operators per site are contracted into a new local site operator $\hat{u}_{\nu_{12}}=\hat{o}_{\nu_{1}^{1}}\hat{o}_{\nu_{1}^{2}}$. Note that, in the latter case, the local operators are ordered in a way that avoids evaluation of on-site commutators $\left[\hat{o}^{(k)}_{\nu_{i}^{1}},\hat{o}^{(k)}_{\nu_{j}^{2}}\right]$. Finally, each string $\hat{\gamma}_{\Delta}$ is converted into a single-branch graph by introducing a set of states $\left\{A^{\Delta}_{i}\right\}_{i}$ with the transitions between the $A^{\Delta}_{i}$’s properly chosen from the corresponding new local site operators $\left\{\hat{u}_{\nu_{12}}\right\}$ (see fig. 12(g) for an example).