QSpace – An open-source tensor library for Abelian and non-Abelian symmetries

Andreas Weichselbaum$\star$

Department of Condensed Matter Physics and Materials Science,

Brookhaven National Laboratory, Upton, NY 11973-5000, USA

⋆<EMAIL_ADDRESS>
## Abstract

This is the documentation for the tensor library QSpace (v4.0), a toolbox to
exploit ‘quantum symmetry spaces’ in tensor network states in the quantum
many-body context. QSpace permits arbitrary combinations of symmetries
including the abelian symmetries $\mathbb{Z}_{n}$ and $U(1)$, as well as all
non-abelian symmetries based on the semisimple classical Lie algebras:
$A_{n}$, $B_{n}$, $C_{n}$, and $D_{n}$, or respectively, the special unitary
group SU($n$), the odd orthogonal group SO($2n+1$), the symplectic group
Sp($2n$), and the even orthogonal group SO($2n$). The code (C++ embedded via
the MEX interface into Matlab) is available open source as of QSpace v4.0 on
bitbucket [1] under the Apache 2.0 license.

QSpace is designed as a bottom-up approach for non-abelian symmetries. It
starts from the defining representation and the respective Lie algebra. By
explicitly computing and tabulating generalized Clebsch-Gordan coefficient
tensors, QSpace is versatile in the type of operations that it can perform
across all symmetries. At the level of an application, much of the symmetry-
related details are hidden within the QSpace C++ core libraries. Hence when
developing tensor network algorithms with QSpace, these can be coded (nearly)
as if there are no symmetries at all, despite being able to fully exploit
general non-abelian symmetries.

###### Contents

  1. 1 Introduction
    1. 1.1 Version history
    2. 1.2 Target audience and additional literature
    3. 1.3 Format conventions
    4. 1.4 Outline
  2. 2 General QSpace Approach and Conventions
    1. 2.1 QSpace stores tensors (not state spaces)
    2. 2.2 Composite index for state spaces
    3. 2.3 Direction of legs for tensors
      1. 2.3.1 Convention: Incoming leg to tensor is written as superscript index
      2. 2.3.2 Leg directions in QSpace
    4. 2.4 Conventions in pictorial representations
    5. 2.5 Tensor conjugation (denoted by '*')
    6. 2.6 $1j$ symbols (individual ireps)
    7. 2.7 Identity operator vs. $1j$ tensor and marker for dual state space
    8. 2.8 Generating scalar (interacting) Hamiltonian terms from irops
    9. 2.9 Compact symmetry labels for non-abelian symmetries
    10. 2.10 Outer multiplicity
    11. 2.11 General tensor decomposition
    12. 2.12 CGT orthonormalization convention
    13. 2.13 CGT data storage and precision
    14. 2.14 Order of legs (index order conventions)
    15. 2.15 Fully contracted CGT networks and $M$-symbols
    16. 2.16 Pairwise tensor contractions and X-symbols
    17. 2.17 Index tags (itags)
  3. 3 QSpace Essentials
    1. 3.1 Starting point: getLocalSpace()
      1. 3.1.1 Example: Single fermionic site
      2. 3.1.2 Example: Spin models
      3. 3.1.3 Beyond models implemented in getLocalSpace
    2. 3.2 QSpace display
  4. 4 Collection of Examples
    1. 4.1 QSpace tensors of rank $r<2$
    2. 4.2 Permutations and index order
    3. 4.3 $1j$ tensor and trailing marker characters
    4. 4.4 Spin-half spin operators
  5. 5 Simple Tutorials and Applications
    1. 5.1 Iterative build-up of many-body state spaces: Two fermionic sites
      1. 5.1.1 Changing local state space or symmetries
    2. 5.2 Swap operator for two spins $S=1$
    3. 5.3 Tweaking local state space
      1. 5.3.1 Project fermionic setup to spin model and complete operator basis
      2. 5.3.2 Appending additional symmetries
    4. 5.4 Local operators beyond getLocalSpace
    5. 5.5 Building rank-5 PEPS tensor
    6. 5.6 fdm-NRG: Iterative diagonalization of Wilson chain and spectral functions
  6. 6 Conclusion
    1. 6.1 Acknowledgements
    2. 6.2 Feedback and support
  7. A General Definitions
    1. A.1 Acronyms
    2. A.2 Terminology
    3. A.3 QSpace data structure
  8. B QSpace Installation and Environment
    1. B.1 System requirements
      1. B.1.1 Matlab and Unix-like environment
      2. B.1.2 C++ libraries required for MEX compilation
    2. B.2 Download and install
      1. B.2.1 System environment setup
      2. B.2.2 MEX compilation
      3. B.2.3 Matlab startup
    3. B.3 QSpace Environment within Matlab
    4. B.4 MEX environment
    5. B.5 HPC environment and parellization
    6. B.6 Environmental variables
    7. B.7 RC_STORE database
    8. B.8 Logging
  9. C QSpace applications in the repository
    1. C.1 DMRG simulations
    2. C.2 NRG impurity solver
  10. D Additional information and documentation
    1. D.1 Additional usage and help info
    2. D.2 Docu/ folder
  11. E Listing of MEX Routines
    1. E.1 QSpace core MEX routines
    2. E.2 NRG related hard-coded MEX applications
    3. E.3 QSpace MEX utility routines
  12. F QSpace Related Matlab Environment
    1. F.1 Class/@QSpace
    2. F.2 Class/@SymOp
    3. F.3 Class/@SymStore
    4. F.4 Class/@Hamilton1D

##  1 Introduction

The treatment of tensor network states in numerical simulations requires a set
of standardized elementary routines: ways to initialize tensors, Kronecker
products to build tensor product spaces, contractions to compute matrix
elements or expectation values, etc. Together these routines form a tensor
library that is designed to provide a user-friendly environment. While the
individual tensor operations are straightforward linear algebra operations per
se, the complexity of a tensor library increases significantly if one intends
to exploit symmetries. Yet symmetries are important, both for the sake of
numerical efficiency, but also to gain detailed symmetry-resolved physical
insights into the fabric of particular tensor network state realizations. By
now several mature open-source tensor libraries are available. Yet the
overwhelming majority of these only permit one to exploit abelian symmetries.
This includes, for example, iTensor [2, 3], Uni10 [4], ALPS [5], or TeNPy [6].
Implementations of non-abelian symmetries in tensor network libraries are
scarce, with TensorKit [7] a notable recent exception.

QSpace also started from plain abelian symmetries originally in version 1
(v1). Yet with v2, QSpace implemented non-abelian symmetries [8]. It was
designed as a bottom-up approach, based on the realization [9] that, as a
generalization of the Wigner-Eckart theorem [10, 11, 12, 13, 14], Clebsch-
Gordan coefficients of an arbitrary-rank symmetric tensor can be factorized in
the form of a tensor product with reduced matrix elements. QSpace explicitly
computes and tabulates Clebsch-Gordan coefficients for standard tensor product
decomposition (rank-3), as well as higher-rank Clebsch-Gordan tensors (CGTs)
obtained via contractions. It is in this sense that QSpace is considered a
bottom-up approach: by explicitly computing and utilizing CGTs, the tensor
representation is known in full numerical detail which then enables general
tensor operations in a self-contained fashion. Specifically, QSpace does not
rely on or use $6j$-symbols [9, 15], bearing in mind that these are known
analytically for SU(2) only, in which case these are simple numbers. For
general non-abelian symmetry, however, the $6j$ symbols become rank-4 tensors
in outer multiplicity (‘$4$M-symbols’) which introduces significant
complications.

Tensor network simulations typically deal with lattice models. These are
comprised of sites whose state space is considered small (microscopic). By
iteratively adding site after site, ever larger many-body Hilbert spaces can
be built. Eventually, these are truncated in tensor network simulations in a
controlled fashion based on the underlying entanglement structure. This
iterative build-up of a Hilbert space precisely also underlies the design
principle of QSpace for non-abelian symmetries: one starts from small local
units where all symmetry-related aspects can be generated, and thus are known
in a simple, transparent way. With these initial building blocks then one can
proceed to ‘play lego’: one builds ever larger more-complex tensor network
structures. From the QSpace perspective, this iterative nature automatically
generates the respective multiplet spaces as they occur on demand at each
step, based on well-defined elementary steps such as tensor product
decomposition, contractions, etc. Newly generated symmetry related data is
tabulated and stored much like in a database (the RC_STORE in QSpace). QSpace
thus automatically adapts to what is needed in the simulations that are run,
bearing in mind that it is impossible to store all symmetry-related data for
continuous non-abelian symmetries, since that data set is infinite.

###  1.1 Version history

The QSpace tensor library has a long-standing history that started around
2006. Since its inception, it has already been thoroughly scrutinized,
debugged, and optimized. A distinctive feature of QSpace has been to exploit
arbitrary symmetries from the very beginning for the sake of numerical
efficiency. The initial motivation of QSpace was at the interface between the
density matrix renormalization group (DMRG [16, 17]) and the numerical
renormalization group (NRG [18, 13]). Yet since it was designed as a general
tensor library from the very beginning, QSpace is also equally applicable in
other tensor network algorithms. QSpace went open-source as of version 4 [1].

  * •

QSpace v1.* (2006-2011) started out as a tensor library for arbitrary sets of
U(1) abelian symmetries, with applications in the realm of NRG (fdm-NRG [19,
20, 21]), DMRG [22, 23, 24], and their crossover [25, 26].

  * •

QSpace v2.* (2012-2015) newly introduced general continuous non-abelian
symmetries to its capability [8, 27, 28], then including the special unitary
group SU($N$) and the symplectic group Sp($2N$) where $N\geq 2$ is dealt with
as a parameter.

  * •

QSpace v3.* (2015-2022) completed the set of semi-simple symmetries by also
implementing the special orthogonal groups SO($N$), both for even and odd $N$,
noting that these represent different Lie algebras $D_{n}$ and $B_{n}$ with
$N=2n$ or $N=2n+1$, respectively. QSpace v3 also included significant
performance upgrades, with full tabulation of generalized Clebsch-Gordan
tensors (CGTs) as well as their contractions via the introduction of X-symbols
[15].

  * •

QSpace v4 (since 2022) finally was geared towards open source in response to
frequent requests from the community. QSpace v4 is mostly the state of QSpace
v3.2 (last version in v3.*), yet cleaned up in the sense of having completed
and merged still open development branches. It went open-source in 2022 as a
git-repository on bitbucket [1] under the Apache 2.0 license.

###  1.2 Target audience and additional literature

This documentation of QSpace assumes basic familiarity with tensor network
states for correlated quantum many-body states. This includes general tensor
network operations and strategies, as well as their graphical description in
the form of tensor network diagrams. Many excellent introductory papers,
tutorials, and reviews already exist in the literature in this regard, e.g.,
see [29, 30, 6, 31] or the lecture notes in Ref. [32] amongst many others.
This documentation does not intend to replicate these. Instead, it focuses on
the implementation of general symmetries in tensor network states, with a
strong emphasis on the distinguishing feature of QSpace, namely non-abelian
symmetries. This documentation builds on the earlier publications on QSpace in
Refs. [8, 15] which already also put an emphasis on a more self-contained
pedagogical presentation (e.g., see the extensive appendices in [8] including
App. A on Non-abelian symmetries 101). As such, these are highly recommended
additional literature that complement this documentation.

###  1.3 Format conventions

The format for the coding syntax (in this font and color coding) used
throughout this documentation is in intuitive compact Matlab semantics. The
main reason is that this also reflects the QSpace environment that the C++
QSpace core routines are embedded into (the many cross-references within this
documentation via hyperlinks are formatted in blue as shown; non-hyperlinked
terms or highlighted text are differentiated in a somewhat lighter color). A
detailed discussion of the structure of the vast low-level C++ code of the
core routines, on the other hand, is beyond the scope of this documentation.
The Matlab syntax used here is basic and concerns, for example, access to
fields in structures, like X.field, cell arrays Q{i}, the representation of
'strings', or matrix notation such as M(i,:) or M(:,j) for row $i$ or column
$j$, etc. The color format is changed when required for the sake of
differentiation or emphasis of Matlab-specific QSpace topics, e.g., like the
general (C++) QSpace implementation vs. its Matlab QSpace, i.e., @QSpace
counterpart. To visually indicate when the latter terms are hyperlinked, the
color changes to a darker shade. Optional arguments to functions are usually
bracketed as in [$\ldots$]. At the level of the operating system,
(environmental) variables like v are also referenced as $v or $(v) as in bash
or Makefile semantic.

In this documentation, snippets of Matlab code using QSpace will be displayed
based on the minted latex package that also enables syntax highlighting,

⬇

>> fprintf(1,’\nHello world from Matlab prompt!’); >>
[a,b]=helloworld(’some’,’input’,pi);

Here >> indicates the Matlab command line prompt. The second line, for
example, is a simple MEX test routine, that assigns a='some', b='input', etc.
As with the above example, QSpace always assumes Matlab strings in single
quotes ('...'). Textual output from Matlab commands such as the above is
displayed without any further syntax highlighting in the format

⬇

Hello world from Matlab prompt!

Hello world from MEX! (.mexmaci64: nargout=2, nargin=3)

The line numbers to the left are included for ease of reference in subsequent
discussions. These output displays may include further commands shown with a
leading Matlab prompt, together with their respective output. The output is
frequently adapted for the purpose of this documentation, e.g., by skipping
empty or trivial lines such as 'ans = ' indicating the answer to, i.e., result
of a command.

###  1.4 Outline

This documentation makes use of a broad range of semantics and definitions.
Essential terminology and concepts are introduced in detail in Sec. 2. This is
complemented by Apps. A.1 and A.2 which provide glossaries on frequently used
acronyms and general terminology for tensor networks and symmetries,
respectively. These short glossaries are grouped by their meaning, rather than
alphabetically, such that quickly glossing through them from top to bottom
should be meaningful. The bare-bones data structure of QSpace tensors is
discussed in App. A.3.

The general QSpace approach to implementing physical model systems is
summarized in Sec. 3. Because QSpace examples need QSpace tensors that need to
be generated first, more detailed examples for Secs. 2 and 3 are collected
into Sec. 4 for clarity, with cross-references as relevant. Section 5 provides
a compilation of more involved examples, intended as simple tutorials already
geared toward applications. The general part of this documentation concludes
in Sec. 6 which also comments on feedback and support.

System requirements are detailed in App. B.1. Download, compilation, and
QSpace environment are discussed in Apps. B.2 and B.3. Appendix C comments on
state-of-the-art NRG and DMRG applications that are already also present in
QSpace, but whose detailed documentation is beyond the scope of this
documentation. App. D comments on detailed usage information (like ‘man’
pages) and further documentation in the QSpace repository. Appendix E provides
a listing of the most relevant compiled binary QSpace MEX routines. Appendix F
finally gives an overview of the Matlab-intrinsic environment that QSpace is
embedded into.

##  2 General QSpace Approach and Conventions

The simple main idea underlying QSpace is that for any tensor $X$ with an
arbitrary but fixed number $r$ of legs or indices, referred to as (tensor)
rank $r$, the presence of a global symmetry leads to a block decomposition of
the tensor. Then for each block with well-defined symmetry sectors $q$ on all
its legs, the symmetry-related aspects factorize [9, 8],

$\displaystyle X=\bigoplus_{q}\ \|X\|_{q}\otimes C_{q}$ (1)

which is generalized and discussed in significantly more detail in Sec. 2.11
still. Here $\|X\|_{q}$ represents the reduced matrix element tensors (RMTs),
and $C_{q}$ the generalized Clebsch-Gordan coefficient tensors (CGTs). Both
mimic the same structure of the original tensor $X$, and thus, in particular,
are also of the same rank $r$. A scalar rank-2 tensor, like for example a
Hamiltonian, becomes block diagonal. The listing of such non-zero blocks is
denoted by the direct sum ($\oplus$) over symmetry block configurations that
are permitted by symmetry. These are labeled by $q$ here which is assumed
collective over all legs. Simply put, a QSpace tensor consists of a list of
RMTs $\|X\|_{q}$ with references to CGTs $C_{q}$. QSpace takes care of all the
symmetry-related aspects and the required book keeping. In this sense, the
tensor X is also referred to as QSpace X.

The decomposition in (1) is also applicable to abelian symmetries. This is
useful from the coding point of view for the sake of a shared consistent data
structure when having abelian and non-abelian symmetries present at the same
time. Fusing state spaces with abelian symmetries can also be cast into the
framework of CGTs, albeit trivially so by taking $C_{q}=1$ when the set of
labels $q$ collected over all legs is permissible from a fusion point of view,
and $C_{q}=0$ otherwise. These trivial factors can be created on the spot, of
course, with no need to tabulate.

In QSpace, the CGTs $C_{q}$ for non-abelian symmetries are explicitly
constructed on demand once and for all, and maintained in a database for later
reference. In the presence of multiple symmetries, the CGTs in the
decomposition in (1) further factorize, $C_{q}=\bigotimes_{s}C_{q_{s}}^{(s)}$,
with $C_{q_{s}}^{(s)}$ the CGT for symmetry $s$. The symmetry labels for any
symmetry sector on a particular leg then also become the collection of labels
across all $C_{q_{s}}^{(s)}$, i.e., $q\equiv\\{q_{s}\\}$ . All bookkeeping for
this additional layered structure is automatically taken care of by QSpace.

The remainder of this section summarizes frequently used terminology,
conventions, and acronyms that are extensively used with QSpace. For more
elementary tensor-network and symmetry-related definitions, it is advised to
briefly also glance over Apps. A.1 and A.2 before reading on.

###  2.1 QSpace stores tensors (not state spaces)

QSpace is a tensor library. Therefore it represents tensors (blobs in a
pictorial representation), but not state spaces per se (individual isolated
legs or lines in tensor networks, i.e., indices). QSpace operates on tensor
objects that have indices (legs) attached. From the point of view of
symmetries then, e.g., when incoming charge needs to be preserved into
outgoing charge, this necessitates at least two indices. In this sense, QSpace
necessarily represents tensors of rank $r\geq 2$. There are only two trivial
exceptions that permit rank $r<2$. These are:

#####  QSpace tensor with rank $r<2$

  * •

rank 0: when computing expectation values or overlaps, the resulting object is
fully contracted, i.e., has no more open legs. Hence this represents a tensor
of rank $r=0$. For an example, see (LABEL:eq:rank0) where a tensor is fully
contracted with itself.

  * •

rank 1: the only rank-1 tensor that is permitted from a symmetry point of view
is a tensor with a single leg in the vacuum symmetry sector i.e., with
symmetry labels $q=0$ (if there is no charge coming in, there is no need for
it to leave). For an example, see (LABEL:eq:rank1). A rank-1 tensor can occur
naturally out of contractions, e.g., when fully contracting a tensor of rank
$r$ on all its $r$ indices with a tensor of rank $r+1$. On general grounds,
the result of such a contraction can only have a non-zero contribution in the
$q=0$ symmetry sector on the single remaining open leg.

##### Referencing state spaces

To specify a particular state space in QSpace, this needs to be done by
referring to respective legs on an existing tensor (e.g., see getIdentity).
Since objects are referenced when handed over to routines (both Matlab and MEX
routines alike), this is efficient in the sense that tensors do not need to be
copied for this purpose [cf. Matlab’s just-in-time (JIT) concept where copies
are only generated once an object is actually altered]. Nevertheless, if
desired, the minimal setting that gets closest to storing a state space on leg
$l$ of QSpace X of some rank $r\geq l$, is by defining the identity tensor on
that leg in compact diagonal form: diag(getIdentity(X,l)). This also inherits
the dimensionality of each symmetry sector, as well as the index tag (itag)
for that leg.

###  2.2 Composite index for state spaces

In the presence of symmetries, any state space, e.g., as associated with legs
of tensors, needs to be organized into symmetry sectors. When exploiting non-
abelian symmetries, any state $s$ must be identifiable by answering the
following questions: (1) what symmetry sector does it belong to? Within that
symmetry sector, (2) which multiplet does it belong to? And in that particular
symmetry multiplet, (3) what state does it represent? Hence for the
description of any state space $s$ in the presence of non-abelian symmetries,
this naturally acquires a composite index structure that contains the answers
to the three questions above [8],

$\displaystyle|s\rangle\equiv|qn;q_{z}\rangle\ .$ (2)

Here $q$ denotes the combined set of symmetry labels that specify (1) the
symmetry sector for all symmetries included. The index $n$ then specifies (2)
the multiplet in that symmetry sector, and $q_{z}$ spans (3) the internal
structure of the individual multiplets in this symmetry sector $q$. To
emphasize that the symmetry labels can refer to any symmetry or combinations
thereof, QSpace generally uses the letter $q$ for symmetry labels, as in
QSpace for quantum symmetry spaces. For the particular case of a single SU(2)
symmetry, for example, $(q,q_{z})\equiv(2S,2S_{z})$ reflects the spin quantum
numbers.

The composite index structure in (2) can also be trivially applied to abelian
symmetries. Since each abelian ‘multiplet’ has only one state, this also
represents the ‘maximum weight’ state. Hence the $q_{z}$ label may be promoted
to the multiplet label, i.e., $q\equiv q_{z}$. For the case of all-abelian
symmetries, the redundant $q_{z}$ labels together with their respective
trivial CGTs can be skipped, altogether. In QSpace this explicitly happens for
the case of all-U(1) symmetries which results in an empty X.info.cgr field in
the QSpace data structure for QSpace X (this reflects the original setting as
of QSpace v1).

Within the QSpace data structure [see Eq. (LABEL:eq:F:struct) and subsequent
lines for an explicit example], the symmetry labels $q$ in Eq. (2) for a
QSpace tensor X are stored as rows in the matrices X.Q{l} where $l=1,\ldots,r$
indexes the leg of the tensor. There is a one-to-one correspondence of the
rows X.Q{l}(i,:) with the respective RMT X.data{i} which describes the i’th
non-zero block in QSpace X. The index $n$ in the decomposition (2) is not
explicitly stored per se, but is implicit as the index within the blocks of
the RMTs stored in the cell array X.data, e.g., when taking some matrix
element X.data{i}($\ldots,n_{l},\ldots$) with respect to leg $l$ [for an
example, see discussion of output to (LABEL:eq:SdotS)]. Similarly, $q_{z}$ in
Eq. (2) corresponds to an implicit index within the CGTs that also needs to
deal with inner multiplicity [8]. The CGTs are tabulated and stored in
RC_STORE/CStore since fixed by symmetry. Hence an actual QSpace tensor X only
stores the relevant reduced data based on the multiplet structure $(q,n)$.
Concerning the symmetry related aspects ($q_{z}$) it suffices to only
reference the underlying Clebsch-Gordan data, referred to as Clebsch-Gordan
references (CGRs) and stored in X.info.cgr. Overall, this switches the
description of a tensor from a state-based to a multiplet-based setting.

###  2.3 Direction of legs for tensors

All lines in a tensor network diagram are directed, i.e, carry an arrow (e.g.,
see Fig. 2; [15]). This binary scheme of whether an index is ‘incoming’ (+) or
‘outgoing’ (-) from a tensor derives from, and thus is equivalent to the
respective binary concepts of raised vs. lowered indices, or contra- vs. co-
variant indices, but also bra- vs. ket-states. This binary concept is
intrinsic to tensors. For this reason, the direction of legs is meticulously
enforced and tracked for all tensors in QSpace.

The semantical preference on the notion of direction is mainly due to
pictorial representations where directed lines are easier to visualize.
Besides, this also gives a simple pictorial representation of the Einstein
summation convention: a raised index that is summed over with a paired-up
lowered index simply corresponds to a well-defined direction of the line that
depicts that contracted index. This directed line necessarily leaves from one
tensor and and enters another. From the perspective of a tensor, the related
notation of using +/- for directions of its legs, as also used in the code
this way, is motivated from the point of view of symmetries: for example, for
U(1) charge conservation, the total combined incoming charge must equal the
total outgoing charge, i.e., $\sum_{\rm in}q_{\rm in}=\sum_{\rm out}q_{\rm
out}$ which may be trivially rewritten as $+\sum_{\rm in}q_{\rm in}-\sum_{\rm
out}q_{\rm out}=0$. This motivates that from the point of view of a tensor, an
incoming (outgoing) index adds (removes) charge, and hence is associated with
a + (-), respectively.

Figure 1:  (a) Fusion of the local state space $\sigma_{(i)}$ by a tensor
$A_{i}$ as in an MPS. Skipping the index $i$ for readability, the left state
space $l\in L$, and the local state space $\sigma\in s$ are fused into the
state space $r\in R$ to the right of site $i$. Irrespective of leg directions
on $l$ and $r$, $A$-tensors use an LRs index order convention as indicated by
the orange labels. (b) (Irreducible) operator $S$ like a spin operator as in
Eq. (10) acting on local state space $\sigma,\sigma^{\prime}\in s$. By
convention, operators use ($\sigma$,$\sigma^{\prime}$,op) index order
convention as indicated by the orange labels. The operator index (third leg)
may be skipped and thus absent for scalar operators where this represents a
trivial singleton dimension.

####  2.3.1 Convention: Incoming leg to tensor is written as superscript index

Consider an $A$-tensor that fuses a local state space $s_{(i)}$, e.g., with a
matrix product state (MPS) $|A\rangle$ and local state spaces $\sigma\in
s_{(i)}$ for some site $i$ in mind [Fig. 1(a)]. For simplicity, the index $i$
is skipped as well as any other indices present on $A$ for the sake of the
argument here. The local identity operator is given by

$\displaystyle\hat{11}$ $\displaystyle=$ $\displaystyle\sum_{\sigma\in s}\
\underbrace{|\sigma\rangle}_{\equiv|s_{\sigma}\rangle}\underbrace{\langle\sigma|}_{\equiv\langle
s^{\sigma}|}\ \equiv\ |s_{\sigma}\rangle\langle s^{\sigma}|\ \ \equiv\ \ \
\raisebox{-5.16663pt}{\includegraphics[width=108.405pt]{fig_Id}}$ (3)

with implicit summation over pairs of raised and lowered indices if not
explicitly specified (cf. Einstein summation convention). The above starts
with the convention that (i) ket-states have lowered indices, and hence are
considered covariant. Therefore bra-states are contra-variant and written with
a raised index. Applying the identity in Eq. (3) to $\hat{A}$, one obtains

$\displaystyle\hat{A}\ =\ \sum_{\sigma}|s_{\sigma}\rangle\underbrace{\langle
s^{\sigma}|\,\hat{A}}_{\equiv A^{\sigma}}\ \equiv
A^{\sigma}|\sigma\rangle\text{ .}$ (4)

In a numerical context, one always has to use ‘matrix elements’ $A^{\sigma}$
that are obtained in a particular basis. By convention, such tensor
coefficient spaces are written without hats. Now also with graphical
depictions in mind, raised or lowered indices will be assigned directions.
Intuitively from the perspective of a tensor, by convention, (ii) the index of
the local state space $\sigma$ is considered ‘incoming’ to the tensor $A$.
Hence in QSpace an incoming index to a tensor coefficient space [as in
$A^{\sigma}\equiv\langle s^{\sigma}|\,\hat{A}$ in Eq. (4)] corresponds to a
projection onto a bra state, and is written as a raised index. The identity in
Eq. (3) may thus be depicted graphically as shown at its r.h.s., which also
shows that a ket-state in itself is considered outgoing. Eventually, via
contractions, a lowered index needs to be paired with a raised index, or
equivalently, a ket- with a bra-state, or a co- with a contra-variant index,
or an outgoing with an incoming leg.

The incoming index in $A^{\sigma}$ in Eq. (4) is contracted with the outgoing
index of the ket state $|\sigma\rangle$ [see also right term in Eq. (3)]. As
such this represents the ‘operator’ $\hat{A}$ with no open index, thus written
with a hat. It is more of an abstract object that is independent of a
particular basis. In an actual numerical context, however, only matrix
elements occur, such as $A^{\sigma}=\langle s^{\sigma}|\hat{A}$. This way, the
index $\sigma$ as in in $A^{\sigma}$ becomes an open, i.e., uncontracted
index. Similarly, for example, the identity operator $\hat{11}$ as written in
Eq. (3), has no open indices. However, when expressed in a particular basis,
the matrix elements are $\langle
s^{\sigma}|\hat{11}|s_{\sigma^{\prime}}\rangle=\delta^{\sigma}_{\
\sigma^{\prime}}$. This ‘strips off’ the ket and bra from the r.h.s. in Eq.
(3), resulting in  with open indices $\sigma$ and $\sigma^{\prime}$.

Now when applying an operator $\hat{S}$ on the site that is fused by the
tensor $\hat{A}$, the matrix element corresponds to the projection onto the
local state $\sigma$ [cf. Fig. 1(b)],

$\displaystyle\langle\sigma|\,\hat{S}\hat{A}$ $\displaystyle=$
$\displaystyle\sum_{\sigma^{\prime}}\underbrace{\langle
s^{\sigma}|\hat{S}|s_{\sigma^{\prime}}\rangle}_{\equiv S_{\
\sigma^{\prime}}^{\sigma}}\underbrace{\langle
s^{\sigma^{\prime}}|\hat{A}}_{\equiv A^{\sigma^{\prime}}}\ \equiv\
S^{\sigma}_{\ \sigma^{\prime}}\,A^{\sigma^{\prime}}\text{ ,}$ (5)

having inserted the identity as in Eq. (3). Since by convention the index
$\sigma^{\prime}$ is incoming to $A$, it has to leave and thus be outgoing for
the tensor $S^{\sigma}_{\ \sigma^{\prime}}$ via the projection to the ket-
state $|s_{\sigma^{\prime}}\rangle$. Therefore $S^{\sigma}_{\
\sigma^{\prime}}\equiv\langle s^{\sigma}|\hat{S}|s_{\sigma^{\prime}}\rangle$
has bra-state $\sigma$ as incoming, and the ket-state $\sigma^{\prime}$ as
outgoing. The order of the bra indices relative to the ket indices is
irrelevant, in that it can be trivially changed. For physical operators,
nevertheless, by standard convention, bra indices are typically listed before
ket indices. Hence their index order convention is $\langle
1|\hat{S}|2\rangle\equiv S^{1}_{\ \ 2}\equiv S^{1}_{2}$, and therefore

$\displaystyle S^{\sigma}_{\sigma^{\prime}}\equiv S^{\sigma}_{\
\,\sigma^{\prime}}\equiv\langle
s^{\sigma}|\hat{S}|s_{\sigma^{\prime}}\rangle\text{ .}$ (6)

When following the direction of the indices in a pictorial representation, Eq.
(5) ‘flows’ from the left towards the right (upwards in Fig. 2). This is
reverse to how one would read such an expression, e.g., as part of an
expectation value from a ‘time-ordered’ perspective, namely right to left
(downwards in Fig. 2): one starts with the ket-state $|A\rangle$ at the right,
and then applies particular operators in a well-defined sequence, and in this
sense ‘time order’. This proceeds towards the left, until it meets a bra state
at the very left for the case of a matrix element or expectation value.

####  2.3.2 Leg directions in QSpace

The concept of raised or lowered indices applies to every level of a QSpace
tensor. It applies to the tensor as a whole, but also to all RMTs as well as
CGRs individually. The CGRs are references to sorted CGTs (these are tabulated
with sorted $q$-labels to avoid proliferation of entries). Bearing in mind
that the relative index order of raised relative to the lowered indices is
irrelevant, the raised and lowered $q$-labels in CGTs can be simply grouped
without paying attention to their relative order. For example, irops have CGTs
of the form $C^{q^{\prime}}_{qq_{\rm op}}\equiv\langle q^{\prime}|C_{q_{\rm
op}}|q\rangle$. For the case of a scalar operator with a singleton operator
dimension having $q_{\rm op}=0$, this can be reduced to $C_{q}^{q}$. Standard
CGCs have $C^{q_{1}q_{2}}_{q_{3}}$ as they fuse $(q_{1},q_{2})\to q_{3}$.
$1j$-symbols have CGTs $C^{q\bar{q}}$ or their conjugate $C_{q\bar{q}}$ which
fuse $(q,\bar{q})\to 0$ or vice versa, with $\bar{q}$ the dual irep to $q$.
Again since the relative index order of raised relative to the lowered indices
does not matter, this documentation also uses the alternative notation

$\displaystyle(q_{1}^{\,}\ldots q_{l}^{\,}q_{l+1}^{\ast}\ldots q_{r}^{\ast})$
$\displaystyle\equiv$ $\displaystyle(q_{1}^{\,}\ldots q_{l}^{\,}|q_{l+1}\ldots
q_{r})\ \equiv\ (q_{1}^{\,},\ldots,q_{l}^{\,};q_{l+1},\ldots,q_{r})$ (7)

for a rank-$r$ CGT $C^{q_{1}\ldots q_{l}}_{q_{l+1}\ldots q_{r}}$, where the
asterisk * on the l.h.s. explicitly denotes an outgoing index (see also tensor
conjugation). By listing all $l$ incoming legs first, followed by the
remaining $r-l$ outgoing legs with $l\in[0,r]$, this permits the notation on
the r.h.s. This index ordering is adopted in the RC_STORE/CStore to avoid a
proliferation of CGTs. On top of this simple grouping of indices, the sorted
CGTs in the RC_STORE also have the $q$-labels non-trivially sorted in a
lexicographical manner within the incoming as well as the outgoing indices.

QSpace permits an arbitrary permutation of legs of a tensor, with an example
shown in Eq. (LABEL:eq:F2:perm). For this the information on in- or outgoing
cannot be simply grouped as in Eq. (7), but needs to be specified explicitly
with each leg. It is stored with the QSpace tensor as a whole with the
tensor’s itags. By convention, a trailing '*' in the string X.info.itags{$l$}
indicates an outgoing index for leg $l$. If there is no trailing '*', the
index is considered ingoing, instead. QSpace insists that directions are
specified with the itags globally for a tensor, i.e., the field X.info.itags
must be set, at the very minimum containing the conjugate flags for every leg.
Having the directions specified with the itags is important for the purely
abelian setting, since in that case X.info.cgr is not required and can be set
empty. In the presence of non-abelian symmetries, however, the leg directions
specified with the itags is redundant since this information can be derived
from each individual CGR. The $q$-directions specified with the itags via
trailing asterisks for the tensor overall are thus fixed for a given tensor
and must be preserved when modifying its itags. To ensure this, QSpace
provides the function setitags(), with an example shown in Eq.
(LABEL:eq:itags:A).

The $q$-direction for CGTs as referenced by CGRs are collected and encoded as
plain strings qdir='$\eta_{1}\eta_{2}\ldots\eta_{r}$', with
$\eta_{i}\in\\{\mbox{{+}},\mbox{{-}}\\}$. For example, see
X.info.cgr(i,j).qdir for a QSpace X. The $q$-directions may also be derived in
string form from the itags of a QSpace tensor via qdir=getqdir(X,'-s'). For
example, this yields the strings

'+-'for scalar operators '+--'for irops (8) '+-+'for $A$-tensors, etc.

QSpace also uses this string notation to organize much of the folder structure
for each non-abelian symmetry in the RC_STORE. Having sorted CGTs there, this
groups all +’s to the front.

###  2.4 Conventions in pictorial representations

Pictorial representations are very useful in describing tensor network
approaches and algorithms [33, 17, 30]. For example, Fig. 2 shows a tensor
network diagram for an expectation value of some product of interactions
$\hat{S}_{i}^{\dagger}\cdot\hat{S}_{j}$. Based on this, general conventions in
graphical depictions of tensor networks with QSpace are as follows [8, 15]:

  * •

Tensors are denoted as some blob (box, circle, etc.). Each tensor has as many
legs attached as it has indices, the number of which defines its rank. The
legs of tensors represent starting or end points of lines.

  * •

Every line in a tensor network carries an arrow: from the point of view of a
tensor, arrows specify incoming/outgoing legs or, equivalently, raised/lowered
indices. A contracted line has no open ends: it starts from some tensor, and
terminates in another (possibly also the same) tensor. Such a contracted line
represents a summed-over index where, by construction, a lowered index leaving
some tensor needs to terminate as a raised index with some other tensor (cf.
Einstein summation convention).

  * •

The tensor network diagram as in Fig. 2 is read top to bottom: starting from
the ket-state $|A\rangle$ at the top, a set of operators, gates, or
interactions are applied (round objects $S_{i}^{\dagger}\cdot S_{j}$). The
order of the operators applied is relevant, up to trivial shifts along the
vertical line they act upon, as long as they do not cross with other operators
acting on the same line (state space). The diagram is terminated eventually by
a bra-state $\langle A|$ at the bottom. Overall, this results in the desired
matrix elements or expectation value.

  * •

By way of reading the diagram, this suggests an implicit ‘temporal’ order top
to bottom (big green arrow to the left), as defined by the particular order of
operators applied onto $|A\rangle$. With the intuitive convention that the
local state space enters a ket state, however, the ‘temporal’ direction in
Fig. 2 is opposite to the arrows shown with the lines that are associated with
the local state space $\sigma_{i}$.

  * •

When taking a ket-state $|A\rangle$, one implies its coefficients $\langle
s^{\sigma}|A\rangle$ with index $\sigma$. That index represents a projection
onto a bra-state [cf. Eq. (4)]. This reflects the general situation in
numerical simulations that tensors are necessarily always cast into a
particular basis. By attaching legs to a tensor in a tensor network diagram,
one implies that it is cast into a coefficient space in some arbitrary but
fixed basis that may be given a label (cf. itags). By convention the ket-state
like $|\sigma\rangle$ represents an outgoing index. Therefore the tensor
coefficients projected onto a bra-state are reverse and thus ingoing.

  * •

When a tensor is shown as a filled (empty) object such as a box, circle, etc.,
its conjugate tensor is shown as a mirrored empty (filled) shape,
respectively, for visual differentiation (e.g., see $A^{\ast}$, $S^{\ast}$ or
$S^{\dagger}$ in Fig. 2).

Figure 2:  Conventions in tensor network diagrams based on an MPS-based
example for the expectation value of some product of interactions
$\hat{S}_{i}^{\dagger}\cdot\hat{S}_{j}$ (pairs of connected circles) like for
Heisenberg spins. The corresponding mathematical expression is specified to
the left and vertically aligned with the various terms in the pictorial
representation. The definition of the operator $S^{\dagger}$ is depicted to
the right in terms of the conjugate tensor $S^{\ast}$. It is the mirror image
of $S$ reflected horizontally and leg directions reverted, such that the
operator index enters $S^{\ast}$ (and thus also $S^{\dagger}$) from the left,
but exits $S$ from the right (see also Fig. 9 for a similar context). The
state is given by $|A\rangle\equiv\prod_{i}^{L}A_{i}$ (ket-state described by
the MPS at the top) operating on the local state space $\sigma_{i}$ for sites
$i=1,\ldots,L$, having $L=5$. Bond states, e.g., up to and including site $i$
are denoted as $K_{i}$ (like ‘kept’ many-body states up to iteration $i$ in
the presence of truncation). Here $K_{0}$ is assumed the left vacuum state
with trivial bond dimension $D_{0}=1$. Similarly, $K_{5}$ describes the full
MPS $|A\rangle$. Being a single many-body state, hence also $D_{5}=1$. The
outer contractions (dashed lines) are thus trivial. The direction of ‘time’
(green arrow) indicates the order of operators applied to the initial ket-
state $|A\rangle$ at the top.

###  2.5 Tensor conjugation (denoted by '*')

The standard concept of Hermitian conjugate, or equivalently conjugate
transpose, is well-defined only for rank-2 tensors, i.e., matrices or
operators. For an irreducible operator $S$, for example, its Hermitian
conjugate $S^{\dagger}$ = permute($S$,'21*') 111 The potentially missing 3 in
the permutation '21*' $\to$ '213*' is implied in QSpace for trailing indices
that keep their position preserved, i.e., act like an identity permutation.
This way, S may be a rank-2 or rank-3 operator. consists of transposition on
the first two indices ('21'), and conjugation (the trailing '*'). For a
general tensor, however, transposition in itself is not well-defined, as it
generalizes to a non-unique permutation that is chosen depending on the
context. The situation is different, though, in pictorial representations
where the generalization of the transpose to arbitrary-rank tensors is well-
defined in that one simply takes the mirror image of the tensor (cf.
$A^{\ast}_{i}$ or $S^{\ast}$ in Fig. 2; [17, 30, 34]). Nevertheless, this
leaves the freedom to rotate the tensor as a whole in the pictorial setting.
But then this holds for any tensor, conjugate or not. More specifically,
tensor conjugation thus proceeds as follows [15]:

  1. (1)

draw the mirror image of the original tensor in pictorial representations

  2. (2)

complex conjugate the array entries (relevant for complex RMTs only)

  3. (3)

revert arrows on all legs (i.e., swap raised/lowered indices)

  4. (4)

index spaces remain the same, such that, e.g., also symmetry labels are left
intact.

Point (4) is added for emphasis only here concerning symmetry labels. For
better visual differentiation of conjugate tensors concerning point (1), one
may furthermore switch from filled tensor objects to empty (outlined) ones, or
vice versa (see conventions in Sec. 2.4, or $A^{\ast}_{i}$ or $S^{\ast}$ in
Fig. 2 as explicit examples).

The concept of point (1) only applies to pictorial representations, but not to
the practical numerical context. Because the permutation underlying point (1)
is not well-defined, tensor conjugation in QSpace adheres to points (2-4)
only. QSpace preserves the index order in tensor conjugations, unless
specified otherwise by the user (see Sec. 4.2 for examples). Point (1)
reflects the generalization of the transposition of matrices. Effectively,
this fully reverts the order of indices. But for tensors of arbitrary rank,
the starting point of what is considered the first index, is not well-defined
unless one adheres to some arbitrary but fixed convention. For example, there
are many different ways where the mirror plane can be drawn. In Fig. 2, for
example, the $A_{i}^{\ast}$ (bottom) are mirrored vs. their respective $A_{i}$
(top) by a horizontal plane, whereas for the interactions, the operator
$S^{\ast}$ is mirrored by a vertical plane w.r.t. their paired up $S$ operator
(see definition of $S^{\dagger}$ to the right of Fig. 2).

Due to (3), tensor conjugation (or ‘Hermitian conjugate’) where present must
always be explicitly specified and included in QSpace, even if the matrix
elements themselves in the RMTs are not complex (point 2). This is important
since reverting directions is equivalent to raising all lowered indices and
vice versa. This is an intrinsic concept of tensors that needs to be respected
for overall consistency. From the symmetry point of view, directions on legs
frequently also indicate types of orthonormalization. Hence by QSpace
insisting on well-defined directions of legs throughout the tensor network,
point (3) also helps to avoid errors in the coding of tensor network
algorithms.

When conjugating a tensor, the action in steps (3-4) will be also referred to
as ‘conjugating’ indices or legs. For a given leg and symmetry sector $q$,
this is denoted by $q\to q^{\ast}$. It needs to be differentiated from the
dual representation denoted by $\bar{q}$. The latter only comes into the
picture when reverting directions of individual legs by applying, i.e.,
contracting $1j$-symbols. This then not only reverts the direction of the
affected leg, but also maps symmetry labels to their duals, $q\to\bar{q}$.
Hypothetically, if one were to explicitly revert all arrows after tensor
conjugation back to their original direction by applying, i.e., contracting
the respective $1j$-symbols onto every leg, the resulting (QSpace) tensor is
generally different from the original one. Aside from having flipped symmetry
labels to their duals, for the case of complex RMTs the reduced matrix
elements would still also remain complex conjugated, bearing in mind that
$1j$-symbols only have non-trivial structure within the CGTs, but behave like
an identity within the RMTs.

###  2.6 $1j$ symbols (individual ireps)

In the literature of SU(2), $3j$ symbols refer to standard CGCs
$(j_{1}j_{2}|j_{3})$ that fuse angular momenta $(j_{1},j_{2})\to j_{3}$.
Within QSpace, symmetry labels are generally written as $j\to q$ to emphasize
that these may represent symmetry labels for a general (non-)abelian symmetry.
Hence a standard CGC is written as rank-3 CGT $C^{q_{1}q_{2}}_{q_{3}}$. Now
for the case that two multiplets $q_{1}$ and $q_{2}$ fuse into the scalar irep
$q_{3}=0$ [Fig. 3(c)], the respective CGT has a singleton dimension as third
index that may be skipped. This case is possible if and only if $q_{2}$ is the
dual representation to $q_{1}$, i.e., $q\equiv q_{1}=\bar{q}_{2}$. This
results in the CGT $C^{q\bar{q}}\equiv C^{q\bar{q}}_{0}$ which exists and is
unique for every irep $q$, since no outer multiplicity can occur for CGTs of
rank $\leq 2$, the CGTs $C^{q\bar{q}}$ are unique up to the standard CGC sign
convention. In QSpace, these particular rank-2 CGTs $C^{q\bar{q}}$ are
automatically computed and stored for every new irep $q$ that is generated in
a tensor product decomposition for non-abelian symmetries without necessarily
performing the full decomposition of $(q,\bar{q})$ (see Ref. [15] for
details). For abelian symmetries, the $1j$ symbols can be trivially generated
on the fly with no need to store.

Figure 3: Identity and $1j$ symbols – (a) Fusing a multiplet $q$ with a scalar
(singleton $0$) trivially leads to the same multiplet $q$. Hence the
underlying CGT must be the identity $11^{q}/\sqrt{|q|}$ up to the CGT
(ortho)normalization convention. Since the singleton index for the scalar irep
($0$) can be skipped, the same holds for rank-2 CGTs as in (b). The identity
$\sqrt{|q|}\,C^{q}_{q}=11^{q}$ can be trivially inserted into or removed from
any directed line. (c) Construction of the $1j$ symbol $U^{q}$ (‘one $q\to j$’
symbol [35, 36]) from a standard rank-3 CGC $C^{q\bar{q}}_{0}$ up to
normalization involving the dual representation $\bar{q}$. As apparent from
the construction in (d-f), $1j$ symbols represent orthogonal matrices
(unitaries).

Importantly, every $C^{q\bar{q}}$ exactly also represents an orthogonal
transformation up to a normalization factor [Fig. 3(c)],

$\displaystyle U^{q}\equiv\sqrt{|q|}\,C^{q\bar{q}}\qquad\text{ ($1j$
symbol)}\text{ ,}$ (9)

where $|q|$ represents the state space dimension of irep $q$, such that
$U^{q}U^{q\dagger}=U^{q\dagger}U^{q}=11$ (bearing in mind that with all CGTs
real, the dagger is equivalent to a simple transpose). These $U^{q}$ can be
used to revert arrows on lines in a tensor network state, as sketched in Fig.
3(d-f) [15] from the symmetry perspective of CGTs [in practice, reverting
arrows typically also affects the RMTs, e.g., when orthonormalizing state
spaces in the process of shifting the orthogonality center (OC), etc.]. By
inserting the identity $UU^{\dagger}$ into a fully contracted line for irep
$q$ in between, say, tensors $A$ and $B$, splitting this product, and
contracting $U$ and $U^{\dagger}$ onto $A$ and $B$, respectively [indicated by
the shaded gray arrows in the background of Fig. 3(f)], this results in a (i)
reverted arrow on that line, and also on the respective legs of the tensors
$A$ and $B$. Furthermore, by construction, this simultaneously (ii) switches
the symmetry sector on that line to the dual representation, $q\to\bar{q}$.
Because of its usefulness, $U^{q}$ is given a name. Since it only depends on a
single symmetry label $q$ [or ‘one $j$’ symbol for SU(2)] with $\bar{q}$
simply inferred, it is referred to as $1j$-symbol with reference to the
literature of SU(2) [35, 36, 15].

##### Caveat self-dual ireps

In the construction of the $1j$-symbol as rank-2 tensors, the singleton
dimension in the scalar representation $0$ [black dashed line in Fig. 3(e,f)]
was skipped. This results in a caveat, in that the truncated ‘stem’ of this
singleton dimension [remnant of the black dashed line in Fig. 3(f)] needs to
emerge on the same side for both $U^{q}$ and $U^{q\dagger}$, here emerging at
the top. Alternatively, one could have closed the black dashed line in Fig.
3(e) at the bottom.

More specifically, the caveat concerns self-dual ireps here, $q=\bar{q}$, like
all ireps in SU(2). In that case the $1j$-symbol $U^{q}$ becomes
indistinguishable from $U^{\bar{q}}$ in terms of symmetry labels and
$q$-directions, despite that they may differ by a sign [15], since
$U^{q}\equiv C^{q\bar{q}}=C^{qq}$, yet $U^{\bar{q}}\equiv
C^{\bar{q}q}\equiv(C^{q\bar{q}})^{T}=(C^{qq})^{T}=\pm U^{q}$. The last sign is
negative, for example, for all half-integer spins for SU(2). Hence when
inserting $U^{q}U^{q\dagger}=11$ into a contracted line of a tensor network
for the sake of reverting its arrow, proper care must be taken that the same
object is used twice, once with, and once without the dagger. Otherwise, this
can give rise to erroneous signs. As a safeguard in that regard, trailing
primes are introduced with itags to mark dual state spaces, as in
$C^{q\bar{q}}$ even if $q=\bar{q}$.

###  2.7 Identity operator vs. $1j$ tensor and marker for dual state space

The $1j$-symbol was introduced above for a particular multiplet $q$ in some
symmetry. Yet it can be easily generalized to cover full state spaces as
implemented with getIdentity. Based on its origin, this will be also referred
to as $1j$ symbol or $1j$ tensor, nevertheless. The $1j$ tensor is closely
related to the identity operator [cf. Fig. 3]. Hence both of them refer to the
state space associated with a particular leg of a tensor and are also obtained
by the same QSpace routine getIdentity. See Sec. 4.3 for detailed examples. To
start with, both objects are of rank $2$. However, while the identity operator
is a plain scalar operator, the $1j$ tensor is not an operator per se, as it
has both of its indices incoming (or both of them outgoing for its conjugate).
Therefore it is also not necessarily block-diagonal, but rather anti-block-
diagonal with respect to its symmetry sectors.

The identity operator $E$ as obtained by getIdentity adheres to the index
order convention of operators as in Fig. 1(b). By having a scalar rank-2
operator here, $E$ has one in- and one out-going index. The $1j$ tensor is
obtained by simply adding the option '-0' (‘fuse into $q=0$’) to the same call
to getIdentity otherwise. By contrast to the identity operator, this returns a
unitary $1j$-tensor $U$ on the same state space, yet with both indices
ingoing. It is real and therefore satisfies $UU^{T}=UU^{\dagger}=11$. The
QSpace $U$ returned by getIdentity has exactly the same RMTs as the identity
operator. However, instead of identity CGTs, it references the respective $1j$
CGTs required for a given state space, properly normalized via the respective
CGRs. In this sense, the $1j$-tensor $U$ defined over an entire state space of
a leg of a tensor is also frequently referred to as a $1j$ symbol itself.

##### Marking itag of dual state space

When calling U=getIdentity(A,l,'-0') this returns the full state space as
present in QSpace A on leg $l$ (first index in $U$) with its dual on the
second index in $U$. To safeguard against the above caveat, the itag of the
dual state space on the second index is marked with a trailing prime (’).
After all, the dual state space is different from the original state space,
and hence its itag also should be differentiated. The trailing marker
character serves as a toggle: reverting the arrow twice returns to the
original itag. Therefore adding a prime to an existing trailing prime,
annihilates both of them, i.e., foo’’ $\to$ foo. So if by the history of an
algorithm, the input space already has a mark on its itag, then that mark is
removed for the dual space on the second index. Detailed examples are given
and discussed in Sec. 4.3.

###  2.8 Generating scalar (interacting) Hamiltonian terms from irops

A Hamiltonian represents a scalar operator. Hence from the point of view of
symmetries, all terms constituting a Hamiltonian term must be scalars
themselves. With this in mind, general two-body interaction terms in a
Hamiltonian can always be written as or built from scalar bilinears
$\hat{S}_{i}^{\dagger}\cdot\hat{T}_{j}$ acting on sites $i$ and $j$, also
permitting $i=j$. Typically, $\hat{S}=\hat{T}$, with pictorial examples shown
in Fig. 2, and explicit QSpace contractions to describe such pairwise
interactions provided in Eq. (LABEL:eq:FdotF) or Sec. 5.1, e.g., Eqs.
(LABEL:eq:SdotS)–(LABEL:eq:FdotF:T12) together with Fig. 9.

From a symmetry point of view, the bilinear
$\hat{S}_{i}^{\dagger}\cdot\hat{T}_{j}$ requires (i) that the operators
$\hat{S}$ and $\hat{T}$ are irops, or more generally operators. In the case of
a pair of irops, (ii) these must transform according to the same operator irep
$q_{\rm op}$, as this is summed over (contracted) in the dot-product
$\hat{S}_{i}^{\dagger}\cdot\hat{T}_{j}$. In addition, for the contraction to
be valid, the index must be out-going from one irop, and incoming to the other
irop. Therefore the general natural structure includes the conjugate tensor,
i.e., always makes explicitly use of the dagger, as in
$\hat{S}_{i}^{\dagger}\cdot\hat{T}_{j}$. The dagger may equally well be
written with the second operator, instead. For the example, consider the SU(2)
spin operator written as the irop that transforms like an $S=1$ multiplet

$\displaystyle S$ $\displaystyle=$
$\displaystyle\begin{pmatrix}-\tfrac{1}{\sqrt{2}}S_{+}\\\ S_{z}\\\
\tfrac{1}{\sqrt{2}}S_{-}\end{pmatrix}$ (10)

where the minus sign with the first component originates from the application
of raising or lowering operators [8]. The inverse square root factors ensure
the simple dot-product structure for isotropic spin interactions, having
$\hat{S}_{i}^{\dagger}\cdot\hat{S}_{j}=\frac{1}{2}(\hat{S}_{i+}\hat{S}_{j-}+{\rm
H.c.})+\hat{S}_{iz}\hat{S}_{jz}=\hat{S}_{ix}\hat{S}_{jx}+\hat{S}_{iy}\hat{S}_{jy}+\hat{S}_{iz}\hat{S}_{jz}$.
The components of the spin operator in Eq. (10) are evidently non-hermitian,
and thus do require the dagger with $\hat{S}_{i}^{\dagger}\cdot\hat{S}_{j}$.
For general operators, if the scalar result out of
$\hat{S}_{i}^{\dagger}\cdot\hat{T}_{j}$ is non-hermitian [e.g., as is the case
for fermionic hopping described by $S\to F$ with $F$ some set of fermionic
annihilation operators; see Eq. (LABEL:eq:FdotF:T12) and subsequent
discussion], then the Hermitian conjugate
$(\hat{S}_{i}^{\dagger}\cdot\hat{T}_{j})^{\dagger}=\hat{T}_{j}^{\dagger}\cdot\hat{S}_{i}$
needs to be added if part of a Hamiltonian. Hence a typical Hamiltonian can be
written as a sum of scalar operator contributions,

$\displaystyle\hat{H}=\sum_{i\alpha}\varepsilon^{\alpha}_{i}\hat{n}_{\alpha
i}+\sum_{ij,\alpha}V_{ij}^{\alpha}\ \hat{S}_{\alpha
i}^{\dagger}\cdot\hat{T}_{\alpha j}\ \ [+\ {\rm H.c.}]$ (11)

which may be further extended to also include products of terms such as the
above. Here
$\hat{n}_{\alpha}\equiv\hat{c}_{\alpha}^{\dagger}\cdot\hat{c}_{\alpha}$
describes some set of local scalar operators, and $\hat{S}_{\alpha}$ and
$\hat{T}_{\alpha}$ some set of local irops, typically having $\hat{S}=\hat{T}$
for a particular interaction. The Hermitian conjugate is required for all
interactions whose operator irep $q$ is not self-dual. For example, it is not
required for SU(2) spin interactions $\hat{S}_{i}^{\dagger}\cdot\hat{S}_{j}$.
It is required for fermionic hopping, though, when using U(1) charge symmetry,
since the annihilation operators $F$ transform like charge $q=-1$ (they reduce
charge by $1$), whose dual are creation operators with $\bar{q}=+1$.

Taking the irop $S$, its daggered version $S^{\dagger}$ could be explicitly
constructed, in principle: the tensor $S$ would have to be conjugated which
reverts all legs, including the irop index. By subsequently applying a $1j$
symbol $U$ onto the irop index 3 and permuting indices 1 and 2, this brings
the irop $\tilde{S}\equiv US^{\dagger}$ back to canonical QSpace form for
irops, with the crucial difference, that the irop now transforms according to
$\bar{q}$. In practice, such a transformation to obtain a canconical form of
$S^{\dagger}$ by using $1j$ symbols etc. is never required when computing
$\hat{S}_{i}^{\dagger}\cdot\hat{T}_{j}$ [e.g., see Fig. 9 and corresponding
text]. Much to the contrary, it would be rather impractical and prone to
errors. Having both, $S$ and $\tilde{S}\equiv S^{\dagger}$ in canonical irop
form, the irop index is outgoing in both tensors. This would prevent a simple
contraction as in $\tilde{S}\cdot S$ on the irop index, because leg directions
are incompatible. Considering $\tilde{S}^{\dagger}\cdot S$, instead, this
would be contractible in principle for self-dual $q_{\rm op}$, such as for the
SU(2) spin operator, but wrong in most cases, in the sense that it gives a
different operator $\tilde{S}^{\dagger}\cdot S=(SU^{\dagger})\cdot S\neq
S^{\dagger}\cdot S$.

Hence there shall never be a need to explicitly construct $\tilde{S}\equiv
S^{\dagger}$ as canonical irop based on $1j$ symbols, etc., when computing
interaction terms in a Hamiltonian. What is solely required is to contract the
conjugate tensor of $S_{i}$. This can be indicated on the fly as an option
when calling contract itself, while ensuring that the ‘transposed’ indices are
correctly contracted. The latter is also dictated by the directions of the
legs [e.g., see Fig. 9 and corresponding text]. Hence much of this is taken
care of automatically when using auto-contraction based on itags.

###  2.9 Compact symmetry labels for non-abelian symmetries

The motivation of a ‘compact’ notation for symmetry labels ($q$-labels) of
non-abelian symmetries is mainly for readability purposes, e.g., to have well-
aligned displayed tables, but also a means for a simple standardized notation
in text. The compact notation also reflects the situation in practice in the
code, where an $n$-tupel of labels is just a vector of numbers without any
white space or separators. The contact notation permits a lean concise
specification of symmetry labels. It is specific to non-abelian symmetries
since ireps for any abelian symmetry are always described by a single symmetry
label only (see symmetry rank). Therefore plain signed integer notation is
used for abelian symmetries. This is also the only setting where negative
symmetry labels $q<0$ can occur.

#####  QSpace label convention for SU(2): $q=2S$

The concept of compact symmetry labels for non-abelian symmetries relies on
symmetry labels that are (i) non-negative and (ii) integer to start with. This
is the case for all simple Lie algebras, except for the half-integer spins in
SU(2) for historical reasons. The latter ‘exception’, however, is only by
convention. Hence QSpace adopts the alternative labeling convention for SU(2)
which is fully aligned with the labeling scheme for SU($N>2$) ireps more
generally, namely by choosing the SU(2) multiplet label
$q=2S\in\mathbb{N}_{0}$. This is consistent with the symmetry labels for
SU($N$) in general, in that they specify Young tableaus [e.g., see Eq.
(LABEL:eq:YTableau)]: an SU(2) spin $S$ has a Young tableau of a single row
with $q=2S$ boxes. This convention for SU(2) also permits the usage of plain
integer labels for symmetry labels, and thus avoids, e.g., the need to print
half-integers. Half-integer spins $S$ then simply map to odd integers for
$q=2S$.

As a corollary, if an SU(2) symmetry gets broken down to U(1), the symmetry
labels become $q=2S_{z}\in\\{-2S,-2S+1,\ldots,2S\\}$, thus also ensuring
integer labels in the U(1) context for the case of half-integer spins. For
example, a spin $S=1/2$ then acquires the U(1) spin labels $q_{S}\equiv
2S_{z}\in\\{-1,1\\}$.

#####  QSpace label convention for U(1) charge

From a practical point of view, it is desirable in tensor network simulations
to keep symmetry labels around $q\sim 0$. This way the low-energy symmetry
sector remains the same with increasing system size also towards the
thermodynamic limit. This avoids a ‘running’ symmetry label that keeps growing
proportional to block size. For counting particle number (charge), it is
therefore desirable to count particles relative to the average filling, i.e.,
$q\equiv n-n_{0}$ with $n_{0}$ the average filling per site for the entire
system, typically half-filling. For example, for a single spinful fermionic
level, the total occupation relative to half-filling $n_{0}=1$ is given by the
sectors $q_{C}\equiv n-n_{0}\in\\{-1,0,1\\}$, corresponding to empty, half-
filled, and fully occupied, respectively. This reflects the charge symmetry
labels as returned by getLocalSpace.

For the description of an odd number of flavors, however, such as a single
spinless fermionic level, the prescription above with a focus on ‘relative to
half-filling’ would result in half-integer labels for charge. Hence in the
case of an odd number of flavors per site, where $n_{0}$ for half-filling
becomes a half-integer, QSpace adopts the convention $q_{C}\to 2(n-n_{0})$,
which includes a factor of $2$ for the sake of having integer symmetry labels
for charge. For a single spinless fermionic level, therefore getLocalSpace
returns $q_{C}\in\\{-1,1\\}$. As for the case of SU(2) labels above, this
convention of applying a factor $2$ is mainly for convenience and readability,
yet also for consistency across symmetries. For example, SU(2) particle/hole
symmetry always takes the viewpoint of particle number relative to half-
filling. In this sense, when breaking down an SU(2) particle/hole symmetry to
U(1) charge, $q_{C}\to 2(n-n_{0})$ is consistent with the QSpace convention
$q_{S}\to 2S_{z}$ for spin.

On general grounds, U(1) labels may undergo an arbitrary linear map $n\to\
q(n)=bn-a$ with $a,b\in\mathbb{N}$. Thus if the desired target value for the
average filling $n_{0}=a/b>0$ is rational, this ensures that the integer
charge labels $q\in\mathbb{Z}$ remain centered around $q\sim 0$ when
increasing block or system size. Since this is a rather specialized setting,
though, it is not implemented in getLocalSpace as is. Nevertheless, one may
simply tweak its output in this regard.

##### Compact alpha-numeric labeling scheme

By QSpace convention, all symmetry labels for non-abelian symmetries are
compact non-negative integers. Therefore any set of $q$-labels
$q=(q_{1},q_{2},\ldots,q_{\mathfrak{r}})$ with $q_{i}\in\mathbb{N}_{0}$ for a
particular symmetry of rank ${\mathfrak{r}}$ is a meaningful irep. As already
pointed out above, for the example of SU($N$), these $q$-labels directly
specify the Young tableau in the sense that $q_{i}\geq 0$ specifies the offset
(i.e., difference) of boxes from rows $i$ to $i+1$,

This represents the Dynkin symmetry labels for SU($N$) ireps. For a
rank-${\mathfrak{r}}$ symmetry, there are at most ${\mathfrak{r}}$ rows since
any $q_{r+1}$ of completely filled columns, i.e., with ${\mathfrak{r}}+1$
boxes, if present, can be split off as a scalar factor in a product-state like
fashion: a completely filled column represents the largest set of $N$
particles that can be brought into a completely antisymmetrized state which
represents a singlet state. Hence from an irep point of view, the label
$q_{r+1=N}$ is irrelevant and can be skipped, as indicated in Eq.
(LABEL:eq:YTableau) by fading out the leftmost $q_{r+1}$ completely filled
columns.

When exploiting non-abelian symmetries, one usually assumes that there is no
spontaneous symmetry breaking on physical grounds.222 If a non-abelian
symmetry of a Hamiltonian is broken spontaneously in its low-energy regime,
non-abelian multiplets become macroscopically large. In this case, however,
there is no gain in exploiting non-abelian symmetries. Rather, the symmetry
should be reduced, e.g., to an abelian U(1) symmetry in the simulation itself
for the sake of numerical efficiency based on the underlying physics. As a
corollary, the tendency of a system to generate ever larger multiplets during
initialization of a tensor network wave function, e.g., via iterative
diagonalization in an NRG-like fashion even for a uniform system, is a strong
indicator for a spontaneously symmetry broken ground state. Therefore assuming
a ‘non-magnetic’ state, the symmetry labels $q_{i}$ do not grow
macroscopically large with system size, but are rather spread around
$q_{i}\geq 0$. Then if the values for $q_{i}\in\mathbb{N}_{0}$ do not grow
macroscopically large, QSpace adopts the following extension to the
hexadecimal format for the sake of a compact notation: the alphabetic
characters A-Z are used do describe the range $q_{i}\in[10,35]$. If the
operating system supports a case-sensitive file system, the range is further
extended to $q_{i}\in[36,61]$ which maps to the lower-case alphabet a-z,
respectively (the reference to file system comes into play here since
$q$-labels are also inserted into file names within the file-based database
RC_STORE). For example then, for SU(2) the irep q='Y'$\equiv 34$ corresponds
to a multiplet with spin $S=34/2=17$. Therefore if a given a set of $q$-labels
satisfies $q_{i}\in\mathbb{N}_{0}\leq q_{1}^{\rm max}$ for all
$i=1,\ldots,{\mathfrak{r}}$ with

$\displaystyle q_{1}^{\rm max}=\left\\{\begin{array}[c]{ll}61&(=9+2\times 26)\
\text{ for a fully case-sensitive file system, e.g., Linux}\\\ 35&(=9+1\times
26)\ \text{ otherwise, like MacOS}\end{array}\right.$ (15)

the above compact extended hexadecimal notation is adopted for this particular
set of symmetry labels in QSpace.

The range in Eq. (15) covers the vast majority of tensor network simulations.
For example for SU(2), typical tensor network ground state calculations in the
singlet symmetry sector $q=0$ stay within $q\lesssim 10$ i.e., $S\lesssim 5$.
This also generalizes to general non-abelian symmetries, typically having
$q_{i}\lesssim 10$, yet with significantly more possible combinations due to
the presence of $i=1,\ldots,{\mathfrak{r}}$ symmetry labels. Hence, the
compact symmetry labels overwhelmingly stay within the ‘numeric’ range
$0\ldots 9$ for all $q_{i}$, yet permit to exceed this range significantly
within the compact notation by including alphabetic characters. Overall this
motivates the compact alpha-numeric notation of symmetry labels

$\displaystyle q$ $\displaystyle\equiv$
$\displaystyle(q_{1},q_{2},\ldots,q_{\mathfrak{r}})\ \equiv\ (q_{1}q_{2}\ldots
q_{\mathfrak{r}})\ \equiv\ \texttt{\textquotesingle}q_{1}q_{2}\ldots
q_{\mathfrak{r}}\texttt{\textquotesingle}$ (16)

for a particular irep $q$, with all spacing (commas) in between the $q_{i}$
removed, and written with or without brackets, like simple strings. For the
case that $q_{i}>q_{1}^{\rm max}$, a more extended display with regular
integers is chosen with the $q_{i}$ separated by spaces.

For example, the scalar representation is given by $q_{\rm scalar}=$ $(0\ldots
0_{\mathfrak{r}})\equiv 0$, which is also conveniently denoted by a single $0$
in this documentation for any symmetry for readability. The defining
representation is given by $q_{\rm def}=(10\ldots 0_{\mathfrak{r}})$. For
example, for SU($N$),

$\displaystyle q_{\rm def}$ $\displaystyle=$
$\displaystyle\left\\{\begin{array}[c]{lll}(1)&\text{equivalent to
}S=q/2=1/2&\text{SU(2)}\\\ (10)&\equiv(1,0)&\text{SU(3) }\\\
(100)&\equiv(1,0,0)&\text{SU(4) }\\\ \ \ \ldots&\quad\ \
\ldots&\text{etc.}\end{array}\right.$ (21)

The dual representation for SU($N$) simply flips the order of labels in the
standardized scheme above, i.e.,
$\bar{q}=(q_{\mathfrak{r}}q_{{\mathfrak{r}}-1}\ldots q_{1})$.

Figure 4:  Outer multiplicity (OM; figure adapted from [15]) based on the
simple specific example of a rank-4 CGT of SU(2) in panel (a) where the CGT
contraction $(SS|S_{i})\ast_{i}(SS|S_{i})^{\ast}\equiv(SS|SS)_{i}$ on the
intermediate index for fixed intermediate multiplet $S_{i}\in 0,\ldots,2S$
gives rise to outer multiplicity for the rank-4 CGT
$\tilde{C}^{qq}_{qq,S_{i}\to\tilde{\mu}}$ with $q=2S$ and
$\tilde{\mu}=1,\ldots,M$ having $M=2S+1$. (b) More generally, contracting two
CGTs $C$ and $C^{\prime}$ on two shared indices with fixed intermediate
symmetry sectors $q_{i}$ and $q_{j}$. Each of the CGTs $C$ and $C^{\prime}$
can have OM on their own, for general non-abelian symmetries already so at the
level of rank-3 CGTs as in (a). The combined multiplicity is therefore
determined by obtaining the resulting CGTs $\tilde{C}$ for all combinations
$(\mu,\mu^{\prime},q_{i},q_{j})$. If multiple intermediate indices are
contracted simultaneously, like $q_{i}$ and $q_{j}$ above, then the plain
contractions $\tilde{C}(\mu,\mu^{\prime},q_{i},q_{j})$ are generally not
orthogonal and typically overcomplete. Hence these need to be orthonormalized
into an arbitrary but fixed combined OM basis indexed by $\tilde{\mu}$ [cf.
Eq. (27)].

###  2.10 Outer multiplicity

Outer multiplicity (OM) is ubiquitous when dealing with non-abelian symmetries
in tensor network algorithms. It adds significant overhead in the book-keeping
when dealing with symmetries, yet is fully taken care of by QSpace in its MEX
core routines. The discussion here therefore serves as background information
only. While it is important to be aware of on general grounds when
interpreting QSpace tensors, the actual book-keeping is hidden from, and in
this sense of no further concern for the typical QSpace user.

To be specific, consider a rank-4 tensor having SU(2) symmetry with two
incoming legs and two outgoing legs which, for simplicity, all carry the same
spin $S$. This is described by $(S,S|S,S)$ or, equivalently, the CGT
$C^{qq}_{qq}$ with $q=2S$ in QSpace [cf. Fig. 4]. Now since $(S,S)$ can be
fused into $M=2S+1$ intermediate spins $S_{i}=0,1,\ldots,2S$, each of these
can again be ‘unfused’ symmetrically into the outgoing indices. Hence one may
contract $(SS|S_{i})\ast_{i}(SS|S_{i})^{\ast}\equiv(SS|SS)_{i}$ on the
intermediate index for any fixed $S_{i}$ (this represents a projector onto the
intermediate multiplet $S_{i}$ of dimension $2S_{i}+1$). Then
$i\to\mu=1,\ldots,M$ takes the role of an outer multiplicity index: for given
CGT $C^{qq}_{qq}$ with fixed ireps on all legs, there exist multiple
orthogonal CGT’s permissible by symmetry. Hence the CGT $C_{qq}^{qq,\mu}$
naturally acquires an additional index, the OM index $\mu$ which spans the OM
space of dimension $M$. The CGT $C_{qq}^{qq}$ is unique then, up to an
arbitrary orthogonal rotation in OM space [this generalizes the sign
convention of standard rank-3 CGCs for SU(2)]. The choice of basis within this
OM space in QSpace is discussed with the CGT (ortho)normalization convention.
Overall then, one has $M>1$ in the presence of OM, $M\leq 1$ in the absence of
OM, and $M=0$ if a particular CGT does not exist at all, i.e., is not
permissible from a symmetry point of view.

For abelian symmetries, one always has $M\leq 1$, i.e., outer multiplicity
does not occur. It is clear from the above example, however, that OM is
already ubiquitous even for SU(2) for tensors of rank $r\geq 4$. For general
non-abelian symmetries of symmetry rank ${\mathfrak{r}}>1$, OM already occurs
at the fundamental level of standard CGCs, i.e, CGTs of rank $r=3$ [37, 38, 8]
(there can never be OM for $r\leq 2$, like scalar operators).

##### OM index is specific to tensor (not particular leg)

Standard CGC decomposition $C^{q_{1}q_{2}}_{q_{3}}$ may tempt one to assign
the OM index $\mu=1,\ldots,M$ with the fused multiplet $q_{3}$, as there are
$M$ of these occurring in the decomposition. However, by applying
$1j$-symbols, any other raised index in $C^{q_{1}q_{2}}_{q_{3}}$ can be
lowered, while raising $q_{3}$. While this changes symmetry labels to their
duals, this nevertheless also demonstrates that, by construction, precisely
the same multiplicity will occur in all these other instances, too. Therefore
the OM index may be associated with any of the legs, depending on the context.
In this sense, the OM index really belongs to a CGT as a whole, and not to any
particular leg of it. This disregard for the direction of the legs is also
reflected in the CGT norm convention in Eq. (26).

##### OM can become arbitrarily large

For non-abelian symmetries, OM typically grows exponentially with the number
of legs, i.e., the tensor rank [15]. Yet even at the level of standard rank-3
CGCs, OM can become arbitrarily large. For example, for SU(3) the CGT
$C^{qq}_{q}$ with $q=(nn)$ with $n\in\mathbb{N}_{0}$ has an outer multiplicity
$M=n+1$. That is, fusing two SU(3) multiplets $(nn)$ results in $n+1$ distinct
multiplets with the same fused multiplet label $(nn)$ amongst others, i.e.,
$(nn)\otimes(nn)=(00)\oplus(nn)^{n+1}\oplus\ldots$ , where the exponent
indicates OM. The multiplet $(nn)$ has dimension $|q|=(n+1)^{3}$, and is self-
dual, such that the scalar $(00)$ also occurs in the decomposition above. To
be specific, fusing two ireps in the adjoint representation (11), i.e., $n=1$
above, one obtains the full decomposition
$(11)\otimes(11)=(00)\oplus(03)\oplus(11)^{2}\oplus(22)\oplus(30)$, this gives
the adjoint (11) twice. In terms of multiplet dimensionality this reads ${\bf
8}\otimes{\bf 8}={\bf 1}+{\bf 10}+(2\times{\bf 8})+{\bf 27}+{\bf 10}$,
respectively, which adds up to the expected total of $8^{2}=64$ states. Based
on the unboundedness of OM even for rank-3 CGTs, there appears no simple way
to fix a canonical OM basis, e.g., w.r.t. to permutations of legs, etc., since
that symmetry is necessarily finite given the finite number of legs. To the
best of our knowledge, there is no superior systematic way to deal with
arbitrarily large OM for general non-abelian symmetries. Hence a bottom-up
constructive approach is adopted in QSpace that is based on demand, and thus
arbitrary but fixed, at the expense of a history dependent RC_STORE as
discussed with Sec. 2.12.

###  2.11 General tensor decomposition

The general decomposition of a QSpace X in the presence of symmetries was
already briefly touched upon in the introduction [cf. Eq. (1)] in the absence
of OM. In the presence of OM, this becomes [8, 39, 40, 15],

$\displaystyle
X=\bigoplus_{q}\Bigl{[}\|X\|_{q}^{{\color[rgb]{0.439,0.678,0.278}\bf\mu}}\otimes\bigl{(}w^{{\color[rgb]{0.439,0.678,0.278}\bf\mu^{\prime}}}_{{\color[rgb]{0.439,0.678,0.278}\bf\mu}}C_{q{\color[rgb]{0.439,0.678,0.278}\bf\mu^{\prime}}}\bigr{)}\Bigr{]}\text{
,}$ (22)

where the OM indices $\mu$ and $\mu^{\prime}$ are implicitly summed over
pairwise in a plain sum
($\sum_{{\color[rgb]{0.439,0.678,0.278}\mu\mu^{\prime}}}$). With this, the OM
indices are fully contracted in Eq. (22) as visualized in Fig. 5). In
practice, however, the CGTs
$C_{q{\color[rgb]{0.439,0.678,0.278}\bf\mu^{\prime}}}$ are referenced, and the
summation over multiplicity is postponed depending on the context. For
pairwise contractions, for example, this summation is translated to generating
X-symbols [15] from the CGTs and contracting these onto the RMTs
$\|X\|_{q}^{{\color[rgb]{0.439,0.678,0.278}\bf\mu}}$. That is, the tensor-
product in Eq. (22) across RMTs and CGTs is never explicitly expanded per se,
as this would defeat the purpose of exploiting symmetries.

The leading direct sum ($\oplus$) over $q$ in Eq. (22) implies block
structure, which translates to the record index i in the QSpace tensor
structure. Via the Wigner Eckart theorem one obtains the tensor product
($\otimes$) of the RMTs $\|X\|_{q\mu}$ with the CGTs $C_{q\mu}$ (assuming
$w=11$ for simplicity). The RMTs and CGTs can be dealt with separately, and
hence are stored separately. The RMTs can be chosen freely depending on the
physics or the algorithm, whereas the CGTs are fixed by symmetry, and hence
referenced via the CGRs. The CGRs thus store the matrices
$w^{\mu^{\prime}}_{\mu}$, while the RMTs $\|X\|_{q}^{\mu}$ carry trailing open
OM indices $\mu$ for $M>1$.

Figure 5: Pictorial representation of tensor decomposition in Eq. (22) for the
case of an $A$-tensor for fixed symmetry sectors and hence a particular
contribution $q$ for simplicity (figure adapted from [15]). The bare indices
$(l,r,\sigma)$ in the absence of symmetries become composite indices [cf. Eq.
(2)]. The RMTs $\|A\|_{q}$ are stored with a QSpace tensor, whereas the CGTs
$C_{q}$ are only referenced via CGRs which point into the RC_STORE database.
The CGRs also store the matrices $w_{q}$ with the tensor.

#####  $w$-matrix

An important element in the tensor decomposition (22) is the direct link
across the tensor product structure of the RMT and CGT via a plain regular sum
on the OM indices. It simply states that for every OM instance there can be a
different RMT (the additional space within the RMT due to OM only needs to be
explicitly allocated when fusing state spaces [8] which is taken care of by
the getIdentity routine). When summing over the OM index denoted in green in
Eq. (22), one may insert an additional matrix, referred to as the $w$-matrix,
which lives directly on the OM link that connects the RMT to the CGT [cf. Fig.
5]. The $w$-matrix can always be absorbed, i.e., contracted onto the RMT by
carrying out the summation over $\mu$.333 The $w$-matrix cannot be absorbed
onto the CGT, since CGTs are defined once and for all. Besides, the data
format is different: CGTs are stored in higher precision, whereas the RMTs,
same as the $w$-matrices, have plain double precision arithmetic. However,
there is practical value in keeping the $w$-matrix nevertheless: it permits
one to efficiently absorb rotations in OM space, e.g., as they occur out of
permutations of the legs of the tensor, but also to absorb scale factors. For
this reason the $w$-matrix is stored with each QSpace tensor via the CGR. It
is a real matrix with the default normalization convention

$\displaystyle w^{\dagger}w=11\text{ ,}$ (23)

typically taking the first non-zero matrix element in $w$ to be positive. When
$w$ is contracted onto the RMT, it can be simply replaced by $w\to 11$
afterwards. By the convention that the superscript index is the first index,
i.e., $w^{\mu^{\prime}}_{\mu}\equiv w^{\mu^{\prime}}_{\ \ \mu}$, the index
$\mu^{\prime}$ specifies its rows. Then by Eq. (22) the row index of $w$ links
to the CGT, whereas the column index links to the RMT. Typically, the OM
dimension on the RMT is the same as for the CGT, but it may also be smaller.
At any given time, the maximal OM dimension is clearly set by the CGT (OM may
grow dynamically during simulations, until the full OM is exhausted). Yet the
RMT of a given tensor may only relate to a particular subset of OM directions,
hence may have lower OM. For this reason, $w$ can have fewer columns than
rows, i.e., it is a matrix of dimensions $M\times M^{\prime}$ with
$M^{\prime}\leq M$. If $M^{\prime}<M$, $ww^{\dagger}\neq 11$, whereas Eq. (23)
holds in any case.

There are only two exceptions to the normalization convention in Eq. (23).
They are intended for the sake of readability on the Matlab interface only (by
contrast, the QSpace core routines always adhere to Eq. (23) internally). In
these cases the normalization in Eq. (23) is only proportional to the
identity, $w^{\dagger}w=|w|^{2}11$. For QSpace records $q$ where $|w_{q}|\neq
1$, these normalization factors are printed to the right when displaying a
QSpace tensor (see display in Fig. 7 for a specific example):

  * •

rank-2 tensors (e.g., scalar operators): for the sake of the argument consider
an identity operator $X$. In this case, for the sake of readability and also
to avoid confusion, one would like the RMTs themselves to be identities, i.e.,
$\|X_{q}\|=11$ an identity operator acting within symmetry sector $q$. But the
CGTs have fixed CGT normalization convention themselves. In particular, for a
scalar operator like the identity, they are trivially given by
$C_{q}^{q}=\frac{1}{\sqrt{|q|}}11^{|q|}$ with $|q|$ the multiplet dimension.
Therefore in order to reconcile the desired behavior of the RMTs with the CGT
normalization convention, the remaining scale factor can be absorbed at the
level of the CGR with the $w$ matrix by setting

$\displaystyle w_{q}=\sqrt{|q|}\text{ .}$ (24)

This cancels the normalization factor of the CGT, thus resulting in a combined
identity operator, $wC=11$. Since there is no OM for rank-2 tensors, the $w$
in Eq. (24) is always a plain number. Aside from identity operators, the
altered normalization convention on the $w$-matrix for the case of rank-2
tensors intuitively shows the ‘actual’ value of matrix elements in the RMTs
for scalar operators such as Hamiltonians or density matrices.

  * •

rank-3 tensors with only one in- or outgoing index: similar to the previous
case for rank-2 consider here for the sake of the argument the identity
$A$-tensor [8] that describes the fusion of two incoming state spaces without
truncation into their combined state space. Intuitively, here too, one would
like to see parts of sliced up identity matrix in the RMTs. That is, the
entries in the RMTs should be either 0 or 1, as this indicates a map of input
multiplets to fused multiplets. For this, however, the rank-3 CGTs would have
to be normalized like standard CGCs, thus desiring
$\sqrt{|q_{3}|}C^{q_{1}q_{2}}_{q_{3}}$. Again this can be simply achieved from
the CGR point of view by setting

$\displaystyle w_{q}=\sqrt{|q_{3}|}\ 11\text{ .}$ (25)

Now an $A$-tensor is a rank-3 tensor with one outgoing index only. Conversely,
irops are rank-3 tensors, with only one incoming index. In this sense, for
generality, the normalization in Eq. (25) is adopted for all rank-3 tensors
that have one incoming or one outgoing index only. The irep on this leg then
defines the dimensional factor in Eq. (25). In the absence of OM for these
rank-3 tensors [such as for SU(2) throughout], $w$ is again just a plain
number. Otherwise, the $11$ in Eq. (25) is the identity matrix in the OM space
of the underlying $C^{q_{1}q_{2}}_{q_{3}}$ or $C_{q_{1}q_{2}}^{q_{3}}$.

##### Avoiding high-rank tensors

The block decomposition via $q$ in Eq. (22) only includes permissible symmetry
combinations for the given number of legs. However, it is important to bear in
mind in the context of tensor networks that this number of possible
combinations $q$ (or equivalently, the number of records in a QSpace) can
become quickly very large (exponentially so) if tensors have many open legs,
i.e., are high-rank.

This proliferation of non-zero blocks for large-rank tensors is already
present for plain abelian symmetries. Consider, for example, a tensor of even
rank $r$, with $r/2$ incoming and $r/2$ outgoing indices and where for the
sake of the argument, all legs have the same state space, say, with $n_{s}$
symmetry sectors each. Then all $n_{s}^{r/2}$ symmetry combinations of the
$r/2$ incoming legs are permissible, since they can again be ‘unfused’ into
the outgoing legs. Therefore the number of blocks (RMTs) grows exponentially
with the rank. In practice, there will be typically many more combinations
still that satisfy $\sum_{\rm in}q_{\rm in}=\sum_{\rm out}q_{\rm out}$. The
proliferation of non-zero blocks with tensor rank eventually holds generally.
It is also irrespective of the direction of legs since these can be easily
flipped by contracting $1j$-symbols which does not affect the number of non-
zero blocks.

For non-abelian symmetries, the situation is further exacerbated: the fusion
of multiplets typically gives rise to several combined ireps, which thus
further increases the number of internal fusion channels and hence non-zero
blocks (e.g., hundreds of thousands for rank-6 tensors [34]). Furthermore,
when building larger-rank tensors by contraction, also the OM of the
contracted tensors themselves will grow exponentially with increasing tensor
rank [15]. If then one tries, for example, to contract two such large-rank
tensors, this quickly translates into extremely long lists of elementary
pairwise block contractions (billions when contracting rank-6 tensors [34]).
While not necessarily memory consuming, this may take far too long.

Hence reducing the tensor rank of large-rank tensors becomes mandatory. This
can be achieved in two ways: (i) simply fusing legs which in itself has a
significant positive impact on the overall performance [34]. Alternatively, if
the situation permits, one may also (ii) split up large-rank tensors into a
tensor network of several lower-rank tensors such as tensor trains [41]. The
latter has strong links to data compression which is of central importance for
tensor network algorithms such as DMRG/MPS [42, 17, 43] all along, even in the
absence of any symmetry.

###  2.12 CGT orthonormalization convention

All CGTs in QSpace adhere to the Frobenius-like normalization convention,

$\displaystyle\|C_{q}\|^{2}\equiv{\rm tr}(C_{q}^{\dagger}C_{q})\equiv{\rm
tr}(C^{q}C_{q})=1\text{ ,}$ (26)

where the tensor trace is to be interpreted such that it fully contracts
$C_{q}$ with the conjugate of itself on all legs, with only the OM index left
open if present. To start with, Eq. (26) assumes no OM, i.e., $M=1$. The
tensor $C_{q}$ denotes a general CGT with symmetry labels $q$ combined over
all legs with certain directions which are of no further interest for the
argument here. For simplicity, this collective $q$ is written as a subscript
here, where nevertheless its conjugate is symbolically written as $C^{q}$ with
raised $q$.

Importantly, the norm in Eq. (26) is invariant under an arbitrary orthogonal
transformation on any of the tensor’s legs. As this includes reverting
directions of legs based on $1j$-symbols, the normalization in Eq. (26) is
therefore independent of the directions of the legs (while also switching to
dual ireps), fully analogous to Wigner-$3j$ symbols. This bears significant
practical advantages.

In the presence of outer multiplicity, CGTs acquire an additional OM index
$\mu=1,\ldots,M$. By convention, this is taken as trailing and thus slowest
index at position $r+1$ with $r$ the rank of the tensor,444QSpace uses
generalized column-major storage for all tensors, throughout, RMTs as well as
CGTs. While RMTs are stored as full tensors, CGTs are stored in sparse format
generalized to arbitrary rank tensors [8]. In either case, listing the OM
index as trailing index implies that data is blocked in storage w.r.t. the OM
index since slowest. thus having $C_{q}\to C_{q\mu}$. In the presence of OM
then, i.e., $M>1$, the r.h.s. of Eq. (26) becomes an identity matrix in OM
space. That is, with $1\to 11^{(M)}$, the individual OM components can be
fully orthonormalized [15],

$\displaystyle\mathrm{Tr}\bigl{(}C_{q\mu}^{\dagger}C_{q\mu^{\prime}}\bigr{)}\equiv\mathrm{Tr}\bigl{(}C^{q\mu}C_{q\mu^{\prime}}\bigr{)}=\delta^{\mu}_{\mu^{\prime}}\text{
.}$ (27)

As apparent from the above, CGTs are only unique up to an arbitrary but fixed
orthogonal transformation in OM space. By contrast, the basis within the state
space of any leg is defined and thus fixed with the respective multiplet
decomposition which may have to deal with inner multiplicity [8]. In the
absence of OM [such as in standard CGCs for SU(2)], the normalization
convention in Eq. (26) reduces to a sign-convention. This sign convention is
still adopted in any case in the standard fashion of CGCs, namely that the
first non-zero entry in the CGT $C_{q\mu}$ for any $\mu$ is positive.

Aside from the sign-convention on the basis vectors in OM space, this still
leaves considerable arbitrariness in the presence of OM, i.e., $M>1$. This
arbitrariness is fixed within QSpace by its operational procedure: CGTs are
constructed once and for all with iteratively increasing $M^{\prime}$ as they
occur (first come first serve) and then stored in the RC_STORE/CStore [15],
until eventually at $M^{\prime}=M$ the full OM space is exhausted, and thus
complete. This procedure is motivated by the fact, that for larger-rank
tensors, depending on the type of contractions performed, the full OM
dimension $M$ space may never be explored and thus is not required. In this
sense, QSpace does not preemptively generate the full OM space for any CGT
$C_{q}$ the first time it occurs. This is different for standard rank-3 CGCs,
which are derived from an explicit numerical tensor product decomposition. By
construction, their OM space is always complete from the first time they are
encountered. For other rank-3 CGTs that derive out of contractions or higher-
rank tensors, the OM space is built as it occurs out of contractions. If a
rank-3 CGT was first generated via a contraction, and later on is completed
via a tensor product decomposition, the existing CGT always has preference to
ensure consistency within the existing RC_STORE. In this case, any additional
OM space derived from the tensor product decomposition is properly projected
and renormalized.

##### History dependent RC_STORE

By the operational procedure above, namely the deliberate decision to not
generate the full OM space in a unique well-defined manner for any CGT
encountered, the RC_STORE becomes history dependent. As such, the RC_STORE
cannot be simply switched across different tensor network simulations. This
makes access to a global centrally maintained RC_STORE essential. By
meticulously applying and checking IDs together with high-resolution time
stamps with CGTs [15], QSpace insists on a particular RC_STORE database within
a given tensor network simulation in order to ensure overall consistency.

##### Norm of a QSpace tensor

As a corollary, the norm of a QSpace X is defined by the full tensor trace of
$X$ with the conjugate of itself, i.e., the Frobenius norm

$\displaystyle\|X\|^{2}\equiv{\color[rgb]{0.837,0.414,0.135}\texttt{norm(X)}}\verb|^2|\equiv{\rm
tr}(X^{\dagger}X)\text{ ,}$ (28)

With the CGTs normalized as in Eq. (27), when contracted with themselves, they
result in simple identity matrix in OM space. Then if the $w$ matrices on the
OM links also follow their own orthonormalzation convention, $w^{\dagger}w=11$
[Eq. (23)], all that remains for obtaining the norm of a QSpace tensor is
simply to add up the square of the Frobenius norm for all RMTs (for the case
that $w^{\dagger}w\neq 11$ [e.g., see Eqs. (24) or (25)], this needs to be
contracted (or multiplied if a simple number $\propto 11$) onto the RMTs when
computing their norm squared). The normalization conventions therefore
considerably simplify obtaining the norm of a QSpace tensor. In any case, the
norm can be simply obtained using the class routine norm which automatically
takes care of all the above.

##### Representation of OM in QSpace tensors

The tensor decomposition in Eq. (22) has an RMT
$\|X\|_{q}^{{\color[rgb]{0.439,0.678,0.278}\bf\mu}}$ for every symmetry
configuration $q$, i.e., with symmetry labels fixed for every leg. The
combined set of labels for all legs as represented by $q$ thus is unique in
the direct sum ($\oplus$). Consequently, also for a QSpace X the combined,
i.e, catenated symmetry labels Q=X.Q; Q=[Q{:}]; (cf. QSpace data structure and
also subsref) are unique for any QSpace X. That is, QSpace records as of
QSpace v4 have unique symmetry labels for each record (this is in contrast to
v2 [8] or v3 which had OM split up as multiple records in a QSpace tensor with
‘degenerate’ $q$-labels). In the presence of outer multiplicity, i.e., $M>1$,
the OM index $\mu$ is carried as additional open trailing index in the RMTs
(single index fused over the combined OM for all non-abelian symmetries
present). Hence it occurs at index position $r+1$ with $r$ the rank of the
tensor. In the QSpace display, it is indicated as @M with the RMT dimensions
[see Eq. (LABEL:eq:itags:A) and subsequent display for an example]. This OM
index eventually is connected via the intermediary $w$-matrix to the
referenced CGTs which then also carry trailing OM indices.

###  2.13 CGT data storage and precision

All CGTs are real, highly sparse, and computed starting from better than quad-
precision by linking QSpace with the GMP/MPFR multi-precision library [44, 45]
during MEX compilation. The CGTs are also stored in a higher precision format
in RC_STORE/CStore. The major reason for higher precision is the iterative
nature of QSpace [15] that mimics the iterative nature of any tensor network
simulation: starting out from the defining representation, an initial set of
smaller ireps is generated which usually more than suffices for the
description of a physical lattice site. Starting from this local state space
(cf. getLocalSpace), ever more complex tensors can be built by adding sites to
the physical Hilbert space and thus to the tensor network. The process of
tensor product decomposition is fully numerical [8], and hence subject to
numerical precision error. The generalized raising/lowering operators (in the
literature also referred to as simple roots in the root space of a Lie
algebra, [46, 8]) are also computed fully numerically within QSpace [cf.
RC_STORE/RStore] subsequent to a tensor product decompositions via their
projection into the space of the newly encountered ireps [the generalized
raising/lowering operators could be generated directly, e.g., based on
Gelfand-Tsetlin pattern calculus for SU($N$) [38]; but with general non-
abelian symmetries also beyond SU($N$) in mind, this approach is not adopted
in QSpace]. The numerical approach in QSpace also automatically ensures
consistency w.r.t. inner multiplicity (which describes the situation
ubiquitous in rank ${\mathfrak{r}}>1$ symmetries, namely where within a given
multiplet the same $q_{z}$-labels can occur multiple times [37, 38, 8]).
Therefore, by the design choice of QSpace as a true bottom-up approach for
non-abelian symmetries, any inaccuracy at an earlier stage is inherited to any
later stage. In this sense, the accuracy of CGTs deteriorates as the RC_STORE
becomes ever larger. Empirically, however, by starting better than quad
precision (about 40 digits), all CGTs obtained in the history of QSpace so far
have been accurate at least up to 24 digits by monitoring numerical noise. In
this sense all CGTs are exact within double precision accuracy $10^{-16}$.
This is also important when exploiting the sparseness of CGTs since for that
purpose, one needs to be able to reliably distinguish small CGT coefficients
from numerical noise [for example, for SU(4) CGT coefficients as small as
$10^{-9}$ have been encountered].

The sparse CGT data is stored in GMP/MPFR format. The MPFR library does not
support binary storage per se. Hence the numerical data needs to be converted
to int8 strings in a specified base[44, 45]. In QSpace, for compactness, this
numerical base is chosen as 2*3*5 = 30 [8] (for comparison, hexadecimal has
base $16$). The reason for this choice is that it permits exact representation
of fractions based on the prime numbers $2$, $3$, and $5$, which thus includes
binary as well as decimal fractions. Within the mat-binary format, these
strings are automatically compressed, such that eventually the storage
overhead of the string representation is expected to be minor.

##### Sparse CGTs of arbitrary rank

All CGTs are stored in sparse format. For rank $r>2$, the tensors are reshaped
by grouping leading vs. trailing indices, such that the resulting ‘matrix’ is
closest to square. Thereafter the sparse data can be stored as a standard
sparse matrix. By also storing the actual dimensions for all legs of the
original tensor, this permits one to restore the entire sparse tensor of
arbitrary rank by returning to the original rank and recomputing sparse index
positions.

###  2.14 Order of legs (index order conventions)

QSpace supports arbitrary order of legs in a tensor. When permuting the index
order (legs) of a QSpace X, this applies the relevant permutation onto each of
the RMTs X.data{i}, and permutes the symmetry labels for the legs in X.Q{:},
as well as the itags. The general index order is also supported by the CGRs
stored with a QSpace.

##### Index order and CGRs

The CGRs store references to sorted CGTs in the RC_STORE. These CGTs have
sorted $q$-labels to avoid the proliferation of entries. The required
permutation (and possible conjugation) to match the current index order of a
QSpace X is auto-determined and stored with the CGRs internally in the MEX
core routines. When handed over to the user space in Matlab, however, this
permutation is applied to the CGRs in X.info.cgr and hence not directly
visible. An arbitrary permutation of the legs on a QSpace X is translated then
to the respective permutation relative to the underlying sorted CGT. The order
of raised indices relative to the lowered indices is irrelevant and can be
trivially permuted. Subtleties arise, however, for QSpace records with the
same $q$-labels and directions on multiple legs, referred to as degenerate
$q$-labels, as this may lead to rotations in OM space (cf. 2M-symbols). In any
case, all of this is taken care of in the QSpace core routine permute[QS] (MEX
routine permuteQS as called by the Matlab class wrapper routine permute).

##### Compact notation for permutations (with optional conjugation)

The order of legs can be changed using permutations. A permutation of length
$n\leq r$ permutes the first $n$ legs (or all legs for $n=r$). In this sense,
the specified permutation can be shorter than the rank of the tensor. At the
same time, one can also conjugate the tensor as a whole by specifying a
trailing '*' with the permutation. For permutations of length $n\leq 9$, which
is nearly always the case, QSpace accepts a compact string notation for
specifying permutations. For example, '21*' denotes a transposition of the
first two legs, i.e., the permutation p=[2 1] together with overall
conjugation '*'. In QSpace syntax, for example, the Hermitian conjugate of an
irop $S$ can be obtained then via permute(S,'21*') which is equivalent to the
non-compact semantic permute(S,[2 1],'conj').

##### Compact notation for index sets (with optional conjugation)

Specifying a subset of legs permits the same compact notation as the
permutations of legs above, including an optional trailing conjugation flag
'*'. This is relevant, for example, for contractions when several indices are
contracted simultaneously. For example, in Eq. (LABEL:eq:nloc) '13*' specifies
the simultaneous contraction of the preceding tensor on indices $1$ and $3$
while also taking its conjugate.

##### Index order for operators and leg directions

An irreducible operator (irop) consists of an irreducible set of operators
that transform like a multiplet, yet with symmetry operations applied via a
commutator. The set of operators thus acquires an additional index (referred
to as irop index) which, by convention, is always listed last, i.e., third.
The natural representation of an operator in the presence of symmetries
therefore is a rank-3 tensor [cf. Fig. 1(b)]. If the operator transforms like
a scalar, i.e., with symmetry labels $q=0$ on the irop index, this index
represents a singleton dimension which can be skipped. Thus the irop for
scalar operators, such as a Hamiltonian term or density matrix, can be reduced
to rank $2$.

The action of a particular operator is characterized by the symmetry labels on
its irop index. Consider, for example, the matrix elements of an annihilation
operator, e.g., $\langle s|\hat{c}|s^{\prime}\rangle$. Ignoring spin for
simplicity, this operator reduces charge $q$ by 1, i.e.,
$q_{s}=q_{s^{\prime}}-1$. Based on this intuitive notion, namely that the
action of the operator is to lower the charge, this implies the symmetry label
of the irop $q_{\rm op}=-1$. For this to hold, however, the irop index must be
grouped with the ket state, and therefore is also outgoing. Together with the
index order convention $S^{1}_{23}\equiv\langle 1|S_{[3]}|2\rangle$ for an
operator $S$, the leg directions of an irop are thus +--.

Not all rank-3 tensors with leg directions +-- are necessarily operators that
act within a particular state space. Therefore to explicitly declare an
arbitrary QSpace tensor X of rank $3$ with leg directions +-- an irop or
operator more generally, QSpace recommends (and in certain applications
insists on) the additional flag 'operator' with rank-3 operators (see field
X.info.otype in Sec. A.3).

##### Index order convention with $A$-tensors (MPS)

The $A$-tensors as with an MPS adhere to the LRs index order convention [cf.
Fig. 1(a)], having $\langle l|A^{[\sigma]}|r\rangle$ where $l\in L$, $r\in R$
refer to left and right block states, respectively, with $\sigma\in s$ the
local state space. When L$\to$R orthonormalized, $A^{l\sigma}_{r}$ has leg
directions +-+, when L$\leftarrow$R orthonormalized, $A^{r\sigma}_{l}$ has leg
directions -++, and when at the orthogonality center (OC), $A^{lr\sigma}$ has
leg directions +++ (all in).

##### Index order out of contractions

The contraction of a pair of tensors, symbolically written as $A*B$, always
preserves the index order of the non-contracted indices. The ones of the first
tensor ($A$) are collected first, followed by the non-contracted indices of
the second tensor ($B$). Hence the index order out of contractions is fixed by
the order of arguments handed over to contract[QS].

###  2.15 Fully contracted CGT networks and $M$-symbols

When projecting a CGT network onto a particular OM basis of the resulting CGT,
this corresponds to the evaluation of a fully contracted tensor network. It is
relevant, for example, for pairwise contractions of tensors [cf. X-symbols
further below]. With this in mind, consider then such a fully contracted
tensor network of $n$ CGTs. Each participating CGT may carry an OM index $\mu$
that needs to be kept open since it links to an RMT [cf. Eq. (22)] that is
outside the discussion here. In this sense, a fully contracted tensor network
of $n$ CGTs has $n$ open OM indices, and thus is referred to as an $n$M-symbol
with $n\in\mathbb{N}$. The $q$-labels for all contracted legs are considered
additional metadata.

##### 1M and 2M-symbols

The normalization convention in Eq. (27) represents a 2M-symbol
$\delta^{\mu}_{\mu^{\prime}}$. However, since the multiplicity indices always
refer to the OM space of the same CGT, in the spirit of $1j$-symbols, one may
also refer to this as a 1M-symbol. It is always simply a square identity
matrix in OM space. More non-trivially, in the presence of degenerate
$q$-labels on multiple legs with the same directions, a given CGT may be
contractable with a permuted version of itself. This can result in an
orthogonal matrix and thus an actual 2M-symbol.

##### 3M-symbols

When contracting a pair of CGTs and projecting the result onto the target CGT,
this corresponds to a 3M-symbol. In the specific context of tensor
contractions, this is equivalently referred to as X-symbols. Note that by
construction ($n=3$)M-symbols with $n$ an odd number require the presence of
at least some non-rank-3 CGTs in its network, i.e., the presence of non-CGCs.

##### 4M = $6j$-symbols, 6M = $9j$-symbols, etc

In the context of SU(2), a fully contracted CGC network consists of CGTs of
rank $3$ only. For this to be fully contractable, the total number $n$ of CGCs
must be necessarily even. It has a total of $3n/2\equiv 3\tilde{n}$ contracted
legs, each of which carries a multiplet label. Using $j$ with reference to
SU(2), therefore a $(3\tilde{n})$-$j$ symbol is a fully contracted network of
$n=2\tilde{n}$ CGCs. For example, a $6j$-symbol ($\tilde{n}=2$) contracts
$n=4$ CGCs, and thus represents a 4M-symbol. Similarly, a $9j$-symbol
($\tilde{n}=3$) represents a 6M-symbol, etc. In general, a $(3\tilde{n})$-$j$
symbol is a $(n=2\tilde{n})$-M symbol, with all $j$-labels on the contracted
legs considered additional metadata [15]. Therefore fully contracted networks
of rank-3 CGTs correspond to $n$M-symbols with $n$ an even number.

For the simple case of SU(2) where CGTs for rank $r\leq 3$ do not have OM,
$(3\tilde{n})$-$j$ symbols are plain numbers. However, for a more general non-
abelian symmetry of symmetry rank ${\mathfrak{r}}>1$ like SU($N>2$), for
example, the $6j$-symbols are rank-4 tensors with open OM indices $\mu$ on all
its legs, hence 4M-symbols. Since there is no simple canonical basis in OM
space, besides that the OM space even for CGCs (rank 3) can become arbitrarily
large already for SU(3), this makes it considerably more cumbersome to
generalize $6j$-symbols, etc., to arbitrary non-abelian symmetries. Therefore
analytical results for $6j$ symbols are only available for SU(2), but not for
rank ${\mathfrak{r}}>1$ symmetries. For this very reason, QSpace does not use
$6j$-symbols at all. It rather uses 3M-symbols (X-symbols), instead.

###  2.16 Pairwise tensor contractions and X-symbols

When contracting a set of tensors, in practice, this is always carried out as
a sequence of pairwise contractions. For a particular pairwise contraction
then, respective contractions occur in parallel at the level of the RMTs, but
also at the level of the CGTs for each non-abelian symmetry included in the
tensor representation [8].

The contraction of a pair of CGTs $C\ast D$ on a subset of indices results in
another CGT $E$ [cf. Fig. 6]. By projecting the outcome of $C\ast D$ onto $E$
(while completing the OM space of $E$ whenever relevant), the fully contracted
${\rm tr}(E^{\ast}(CD))$ corresponds to a 3M-symbol [Fig. 6(b)]. In the
context of pairwise contractions (‘X’), QSpace refers to these as X-symbols
[15].

In practice, the CGTs $C$ and $D$ in Fig. 6(a) as well as $E$ in Fig. 6(c) are
referenced in QSpace tensors and stored in the RC_STORE once and for all once
generated. Hence all that is required eventually for a pairwise contraction of
QSpace tensors are the X-symbols in Fig. 6(b,c). These are contracted onto the
OM indices of the input pair of RMTs [15], and in this sense fuse their OM
indices.

Figure 6: X-symbols (figure adapted from [15]) – When contracting a pair of
CGTs $C$ and $D$ on some subset of legs [panel (a)] with OM indices indicated
in green, this gives rise to a new CGT $(CD)\equiv C\ast D$. The $q$-labels
and $q$-directions are assumed arbitrary but fixed for all non-green lines,
but not shown for readability. The resulting contracted tensor is expected to
have a single OM index, hence the two green lines in (a) need to be fused, as
indicated in the remainder of the panels. If a CGT $E$ with the same rank and
$q$-labels as $(CD)$ already exists in the RC_STORE (up to sorting) the
contraction in (a) may be simply projected into the existing OM space of CGT
$E$ by inserting the identity $11=E^{\ast}\cdot E$, with the sum in the
product performed over the OM space of $E$ [thick green line in (b)]. If the
CGT $E$ does not yet exist or is incomplete, the additional OM components are
generated via Schmidt decomposition out of $(CD)$. The lower gray-shaded area
in (b) now encompasses a fully contracted network for three CGTs, i.e., a 3M
symbol. For the explicit purpose of pairwise contractions, these are referred
to as X-symbols.

###  2.17 Index tags (itags)

QSpace supports user-defined names on tensor legs with strings up to a length
of 8 characters. Since they represent tags that are assigned to indices, they
are referred to as itags. The usage of itags makes QSpace tensors
significantly more user-friendly. For one, they allow one to identify
individual QSpace tensors more easily in a simulation. Hence the itags are
also shown when displaying a QSpace tensor. Moreover, QSpace supports auto-
contraction based on itags, in that a matching pair of itags with opposite
directions (i.e., one in- and one outgoing) are auto-contracted unless
indicated otherwise. See help on contractQS for more on this.

Within the QSpace core routines, itags are converted to unsigned long
integers, which thus sets the size limit of 8 bytes. If an itag is shorter
than 8 characters, it is padded by zeros. For a QSpace tensor X of rank $r$,
its itags are stored in the mandatory field X.info.itags{$l$} for leg
$l=1,\ldots,r$.

Itags are expected to be printable, mostly alpha-numeric characters in ASCII
code range 33-126, excluding itag separators [,;|] (these characters are
sometimes used to separate multiple itags combined into a single string). It
is recommended to also avoid the trailing marker characters [’*] (see below)
within an itag string, in order to avoid confusion. The above ASCII range is
case-sensitive and excludes white space. The sign bit for each character is
exploited in QSpace core routines for other internally purposes, such as
flagging itags that shall not be contracted. With this then, for example, the
following itags are all different from each other: 's1', 's01', 'S01', 's:1',
's-1', etc.

##### Trailing marker characters in itags

The itags are also used to store the direction of legs for a given QSpace
tensor: a trailing ‘conjugate’ character '*' indicates outgoing legs, whereas
the absence of a trailing '*' indicates incoming. These directions must be
maintained and present even if the itags are empty otherwise. In case of
inconsistency of the leg directions specified with the itags of a QSpace
tensor as compared to the $q$-directions in the CGRs (if present), QSpace will
throw an error.

Moreover, QSpace uses trailing primes to mark dual state spaces. When
reverting the direction of a leg, this changes the interpretation of the
associated block state space. Therefore an automated modification of the itag
is natural. It is useful for better interpretability, but also as a safety
measure that guards against contracting wrong indices that would carry the
same itag otherwise (see also $1j$ symbols).

In summary, QSpace thus assigns special meaning to the following trailing
marker characters in itags:

  * •

A trailing asterisk (*) indicates a conjugate, i.e., outgoing index.
Internally within the C++ core routines, this is split off as $q$-directions
for an entire tensor. Therefore they do not count towards the 8-char limit of
an itag.

  * •

A trailing prime (’) indicates dual state space which occurs when a leg
direction is reverted. The prime is considered integral part of an itag, and
hence counts towards the 8-character length limit. For this reason, a trailing
prime also occurs before a trailing conjugate flag '*' if present.

Both of the trailing marker characters above act as a toggle. Since the dual
of the dual of a state space is the original state space, marking an already
marked itag removes the trailing prime. Similarly, taking the conjugate of an
already conjugate index reverts its direction twice, and therefore results in
the original direction of that leg. Hence in this case, a trailing conjugate
character '*' is removed.
Remark on notation – In this documentation, an itag 'foo' with a trailing
prime is frequently printed as 'foo’'. The outer single quotes '...' denote a
string, and differ in font from the prime (’), for the sake of the discussion.
They are considered the same, though, on the coding level. Specifically, with
Matlab syntax in mind, 'foo’' $\equiv$ 'foo''' $\equiv$ "foo'", yet noting
that as of QSpace v4, Matlab strings are always assumed to use single quotes.

##### Color coding

When using QSpace in interactive Matlab terminal mode, trailing marker
characters in itags are replaced by color coding for the sake of readability:
outgoing indices become gray-ish (while skipping the trailing * character),
and itags marked by the trailing prime (’), acquire a color other than the
default text color of the terminal (darker green as of QSpace v4, thus
replacing the trailing prime; see Sec. 4.3 for examples). This way, for
example, within an MPS, bond itags for $A$-tensors associated with the left
block (left relative to the OC) acquire a different color as compared to the
itags for $A$ tensors associated with the right block.

The color coding of itags, of course, is applicable to non-empty itags only
(i.e., itags that contain actual printable string after removing trailing
marker characters). If all itags are empty, then the +/- notation for
$q$-direction is used as strings when displaying itags (e.g., see Fig. 7) with
color coding applied. When only some itags are empty, no color substitution is
performed, with all trailing markers shown as is. By default, no color coding
is applied throughout for non-interactive batch or cluster jobs where log
output is expected to be redirected to log files. The behavior can be changed
via the environmental variable QS_LOG_COLOR.

##  3 QSpace Essentials

###  3.1 Starting point: getLocalSpace()

The starting point of any tensor network setup in QSpace is getLocalSpace.
This permits one to define a lattice site via a set of operators that act
within that state space. However, because these operators are expressed as
matrix elements in a particular symmetry eigenbasis of the site, this
implicitly also represents its full state space, e.g., via the identity
operator. The way getLocalSpace operates is that it builds the local state
space in Fock space without symmetries first, since everything can be simply
constructed there including the symmetry operators for all requested
symmetries. Based on this, getLocalSpace rotates and groups the basis into
symmetry multiplet spaces, and then ‘compresses’ the data via compactQS based
on the Wigner-Eckart theorem to obtain the RMTs while splitting off the CGCs
as tensor product factors [cf. Eq. (22)].

The overall argument structure of getLocalSpace() is as follows (for more
detailed usage information, see  help getLocalSpace):

[op1,...,opN,Isym] = getLocalSpace(type [, symmetries] [, options]) (29)

The input specifies the type of a site, together with its symmetries, and
further options. It returns a set of representative operators op1,...,opN for
the specified site (therefore N depends on the input), together with an info
structure Isym. The type of a site is a string that can be any of the
following:

The argument symmetries consists of a comma-separated catenated string that
specifies the type and also order of symmetries to use.

##### Spinfull fermions

For spinful fermionic models, i.e., type='FermionS', the following symmetry
entries are supported,

The order in which these symmetry labels are specified also determines the
order of symmetries used, and hence the order of the symmetry labels for any
state space. The last three symmetries only apply in the presence of multiple
channels (or flavors). The number of channels can be specified via the
optional argument 'NC', or directly with the symmetry itself. For example,
'SUNchannel' with 'NC',3 is equivalent to 'SU3channel'. For the case of
${\color[rgb]{0.45,0.45,0.45}\texttt{NC}}>1$ channels, abelian (non-abelian)
charge in all channels individually can be requested via 'Acharge(:)' or
'SU2charge(:)', respectively.

##### Spinless fermions

Spinless fermionic models, i.e., type='Fermion', support the symmetries
'Acharge' and 'SUNchannel' similar to the spinful case above. A subtlety
concerns the case of an odd number of channels NC. There the particle number
relative to half-filling becomes half-integer for the state space of a single
site or, more generally, an odd number of sites. In order to avoid this (also
bearing in mind that QSpace insists on integer symmetry labels), a factor $2$
is applied to the charge symmetry labels. See Sec. 2.9 for more on this.

The fermionic models require that fermionic signs are explicitly taken care of
in the applications via the fermionic parity operator Z. That is, graded
tensors [47, 48] are not implemented in QSpace yet which eventually may allow
one to take care of fermionic signs in an automated fashion. Bosonic degrees
of freedom with a finite local state space may be implemented as tweaks of the
U(1) fermionic models while stripping any fermionic signs based on the parity
operator Z. Bosonic degrees of freedom with infinite state space such as
harmonic oscillators need to be truncated. This can be done efficiently, e.g.,
via a shifted optimal boson basis [49].

##### Spin models

Plain spin models such as SU(2) spin-$S$ Heisenberg models can be setup with
type='Spin'. This typically only allows a single symmetry that is specified as
an additional argument, with the default being 'SU2' (see Sec. 3.1.2 below for
examples). To obtain other more complex spin models based on several
symmetries, the fermionic models may be tweaked e.g., via projections (see
Sec. 5.3.1 for an example).

##### Special case: No symmetries

QSpace also permits simulations in the absence of any symmetry. Yet since this
is a rather atypical case when using QSpace, this proceeds via a tweak. For
example, no symmetries can be implemented via an auxiliary trivial abelian
symmetry such as U(1), where only the scalar symmetry label $q=0$ of the
vacuum is present on any leg. With this, every tensor contains a single QSpace
record only (thus a single full block as RMT). An example is provided with the
spin setting in line (43) below, which returns objects of the above type. As
shown in Sec. 4.4, for spin $S=1/2$, S(i).data{1} contains the Pauli matrices
$[\sigma_{z},\,\frac{1}{\sqrt{2}}\sigma_{-},\,-\frac{1}{\sqrt{2}}\sigma_{+}]$
for $i=1,2,3$, respectively. The order and normalization of operators can
always be simply double-checked by inspection.

####  3.1.1 Example: Single fermionic site

The following call to getLocalSpace returns QSpace operators for a spinful
fermionic site with a single level (i.e., $d=4$ states)

(31)

The detailed structure of the returned QSpace tensors as displayed within
Matlab is discussed by the example of F in Sec. 3.2 (e.g., see Fig. 7) and
also with (LABEL:eq:gls:NC1). The operators returned by (31) are the fermionic
annihilation operators F (irop with two components, one for spin up, and one
for spin down), the fermionic parity operator Z, and the total spin operator S
(hence an irop with three components).

The last argument IS returned by (31) is an info structure that contains
additional detailed information concerning the chosen symmetry setup. As seen
in (LABEL:eq:IS:struct), it contains the identity operator IS.E, or IS.sym as
a more verbose description of the symmetries as specified with the input to
getLocalSpace. This permits one to look up the order of symmetries later,
e.g., to reconfirm whether the charge or spin symmetry is listed first, etc.,
since the symmetries specified with a QSpace tensor, such as A or SU2, no
longer contain that type of information. Other fields in IS also contain low-
level information relevant to the internal workings of getLocalSpace which may
be safely ignored.

The symmetries of a model hold throughout a tensor network simulation. In this
spirit, the call to getLocalSpace defines the symmetry once and for all during
the setup of a model. Hence symmetries can only be switched during the setup
for different simulations. For example, for the same local state space of a
single fermionic site, it may be changed then across the following settings
depending on the symmetries of the model Hamiltonian based on the underlying
model parameters:

All describe the same local state space of $d=4$ states. With increasing
symmetry from top to bottom, however, the local state space dimension is
effectively reduced from $d^{\ast}=4\to 2$ multiplets. For comparison, the
more involved example for NC=3 channels with symplectic symmetry,

[F,Z,S,IS]=getLocalSpace('FermionS','SU2spin,SpNchannel','NC',3);

contains 4*NC=12 components in the irop F (creation and annihilation operators
for each spin and channel). The total spin irop S still has three components.
The scalar parity operator Z shares much of the structure of an identity
matrix except that its diagonal entries are $\pm 1$. As such it also reflects
the full state space of a site. By contrast, the spin operator only acts
within partially occupied sectors since the fully occupied or empty states are
singlets with spin matrix elements all-zero. Yet zero blocks are not stored.
In the present example of NC=3 spinful levels, the largest non-abelian
symmetry is given by the symplectic Sp(2*NC=6) symmetry. Within this setting,
the local space dimension $d=4^{3}=64$ states is reduced to an effective
dimension of a mere $d^{\ast}=4$ multiplets!

####  3.1.2 Example: Spin models

For spin-models, i.e., type='Spin' in Eq. (29), the argument symmetries is
only allowed to contain a single symmetry, assuming SU(2) if not specified.
The options can specify the particular irep to use, assuming the defining
irep, by default. For example, the following two lines are equivalent for
SU(2),

(33) (35)

These return the spin operator for a spin $S$ multiplet in the local state
space with $q$-label $q=2S$. The usage in line (35) permits to switch to other
non-abelian symmetries. For example,

(37)

returns the ‘spin’ operator [synonymous for the irop in the adjoint
representation] for SU(3), with the local state space in the 6-dimensional
representation q='02'. Here the input supports both, compact string ('02') and
non-compact numeric notation ([0 2]).

##### Switch to abelian or no symmetry

In the case of SU(2), the usage in line (33) also accepts the following
switches to abelian symmetries,

(39) (41) (43)

where the option '-A' returns the spin operator in abelian U(1) symmetries
(equivalent to $Z_{\infty}$), '--Z2' returns it with $Z_{2}$ symmetry, and '--
nosym' with no symmetry at all (effectively implemented as trivial ‘U(1)’ with
$q=0$ only). In all cases, the returned spin operator S has three entries
corresponding to S
$=[\hat{S}_{z},\,\tfrac{1}{\sqrt{2}}\hat{S}_{-},\,-\tfrac{1}{\sqrt{2}}\hat{S}_{+}]$.
See Sec. 4.4 for a more detailed discussion.

##### Normalization of spin operator

The spin operator S is a specific irop that transforms in the adjoint
representation. It has as many components as there are generators in the Lie
algebra. The spin operators out of getLocalSpace [e.g., Eq. (37)] adhere to
standard normalization conventions in the literature. This can be easily
double-checked, e.g., via normQS(). For example, for an SU(2) spin in the
defining irep, $\|S\|^{2}\equiv{\rm tr}(S^{\dagger}\cdot S)=(3/4)\,{\rm
tr}(11^{(2)})=3/2$. Therefore it must hold for the respective QSpace spin
operator that norm(S)`^2` = 3/2, seen as $\|S\|=\sqrt{3/2}=1.225$ in line
LABEL:idx:Sop:SU2:hdr below.

The normalization and sign of the spin operator for all ireps is fully
determined by the definition of the spin operator in the defining
representation. After all, the spin operators in any other irep can be
constructed by combining defining ireps and then compute the matrix elements
of $S_{\rm tot}$ in the combined Hilbert space. This construction can be used
to reach any other irep. With this the normalization and sign of $S_{\rm tot}$
is fixed by specifying these for the defining irep. Hence one needs to be
careful when manually constructing spin operators directly for an irep other
than the defining representation, e.g., based on existing CGTs.

Based on the above constructive approach, the spin operator is also well-
defined in the presence of OM for larger ireps: while the spin operator is
unique in the defining irep (this particular case never has OM), OM can arise
for larger local multiplets. For example, consider the SU(3) spin operator in
the local irep q='11' also being in the adjoint representation. This has
$M=2$. Yet from the constructive approach above, there is also only one well-
defined spin operator in this setting, since it can be expressed as $S_{\rm
tot}$ in the tensor product space of three defining ireps (while there will be
two ireps in '11' in this construction, both replicate the same spin
operators, since for both the RMT links the same way to their shared CGR).

####  3.1.3 Beyond models implemented in getLocalSpace

The implemented settings for fermionic or spin models cover a wide range of
applications. Yet, in practice, situations will arise where the desired
symmetry setting is not yet implemented in getLocalSpace. In these cases, it
is often sufficient to tweak the QSpace output of getLocalSpace that captures
the closest most relevant symmetry setting. For this, however, it is important
that one fully understands in detail the general QSpace data structure. See
Sec. 5.3 for detailed examples. As a last resort, one may contact QSpace
support for help. Typical strategies of tweaking the output of getLocalSpace
may include:

  * •

appending additional symmetries [see QSpace/addSymmetry.m].

  * •

projecting fermionic models to a fixed charge sector, thus effectively
obtaining a spin model. If the charge label (relative to half-filling) is non-
zero, the respective symmetry labels (columns in X.Q) may be simply set to
zero because they are all the same due to the projection.

  * •

converting abelian symmetries, like U(1) labels to $Z_{N}$ labels by changing
the symmetry name (cf. X.info.qtype) and converting the symmetry labels in
X.Q) accordingly.

  * •

the absence of any symmetries was already discussed above.

Figure 7:  Typical QSpace display, e.g., as printed when typing the name of a
QSpace tensor on the Matlab prompt (>>) in the user interface. This display
summarizes the tensor’s data structure into a more readable form. The example
here is based on the $A$-tensor as obtained in (LABEL:eq:A2). It fuses the
state spaces of two spinful fermionic sites with $D=4\times 4\to 16$ states [4
x 16 x 4 in line 2 after permutation to standard LRs index order; cf.
(LABEL:eq:A2)]. The symmetry setting of U(1) charge (abelian 'A') with SU(2)
spin is shown as 'A,SU2' in line 1. Hence all multiplets carry two symmetry
labels $q=(q_{1}\ q_{2})$, as reflected in the dimensions of Q in line 1, and
the columns of $q$-labels listed for each leg. The charge $q_{1}$ represents
particle number relative to half-filling and hence can become negative,
whereas ${\hyperref@@ii[idx:SU2:qlabels]{q_{2}=2S}}\geq 0$ reflects the SU(2)
spin multiplet. Due to the presence of SU(2) spin, switching to the multiplet-
based description reduced the Hilbert space to an effective dimension
$D^{\ast}=3\otimes 3\to 10$ multiplets (3 x 10 x 3 in line 2). Since each
fused multiplet carries unique symmetry labels here, this also translates to a
total of 10 QSpace records. The record index is shown to the left for
reference. The RMTs A.data{i} represent reduced matrix elements. Hence they
are tensors with effective dimensions as shown in the second column (this
adheres to Matlab conventions where trailing singleton dimensions beyond