dimension 2 are skipped). The CGT dimensions are displayed in the third column
(shown including trailing singletons for non-abelian symmetries only; if a
QSpace has all-abelian symmetries, the CGT dimensions are skipped entirely).
Long listings are truncated, by default, as in line 9 by ': ...' for the
skipped records 7-8. For a full listing of all QSpace records, one needs to
explicitly call the underlying function display(A,'-f'), since this is the
only way to specify options (like '-f' for full display). The combined scale
factor out of the CGRs [cf. Eqs. (24) and (25)] is shown to the right. Since
it always represents the square root of some (product of) multiplet dimensions
$d\in\mathbb{N}$, it is printed for the sake of readability only for those
records where $d\neq 1$ as $\surd d$.

###  3.2 QSpace display

A summary of the content of a QSpace tensor is displayed by simply typing its
name at the Matlab command line prompt followed by enter (technically
speaking, this is realized by overloading the display function in @QSpace).

For example, Fig. 7 displays the $A$-tensor as generated in (LABEL:eq:A2),
with detailed additional explanations. The symmetry setting is a combination
of U(1) charge with SU(2) spin, with the order fixed by the call to
getLocalSpace in (LABEL:eq:gls:NC1). The display consists a header info
section, followed by a listing of the symmetry blocks including symmetry
labels and dimensional information of the underlying RMTs and CGRs. The pair
of header lines show global tensor information, including itags for all legs
together with their directions. In case of itags that are empty up to trailing
markers, the $q$-directions are shown in the +/- format as this enables color
coding. See also applications for more examples.

The square-root dimensional factors to the right in Fig. 7 are shown only if
the $w$-matrices provided with the CGRs have anomalous normalization, i.e.,
$w^{\dagger}w=d11$ with $d\neq 1$ (combined as a product across all
symmetries). In this case “{$\surd d$}” is included in the display as a
reminder of the altered normalization affecting the RMTs. From Eqs. (24) and
(25), it is generally expected that $d$ reflects some multiplet dimensions,
and therefore is an integer. This motivates the notation with the surd
($\surd$ ) for the sake of readability and compactness.

##  4 Collection of Examples

QSpace uses a compact syntax for tensor network commands with usability and
readability in mind. This compact notation, while it may be self-explanatory
to some extent, takes some time to digest and get used to, though. Hence it is
discussed in detail with the examples in this documentation. This section
provides a combined collection of examples referenced earlier in this
documentation. Essential syntax-highlighted Matlab code is shown in a gray
box. This is followed by a detailed discussion of output objects and sanity
checks on their structure. The source and output lines are given line numbers
to the left for ease of reference.

To provide QSpace examples, the first step is defining a set of exemplary
QSpace tensors. Hence to be specific, the examples below use the QSpace
operators that act within the local state space of a single spinful fermionic
site with U(1) charge and SU(2) spin symmetry,
⬇ >> [F,Z,S,IS]=getLocalSpace(’FermionS’,’Acharge,SU2spin’);  x The local
state space
$s\in\\{|0\rangle,|{\uparrow}\rangle,|{\downarrow}\rangle,|{\uparrow\downarrow}\rangle\\}$
represents empty, singly up/down, and doubly occupied states, respectively (in
contrast to earlier, e.g., Fig. 1, in what follows
$\sigma\in\\{\uparrow,\downarrow\\}\equiv\\{1,2\\}$ specifies the spin flavor
index that underlies the local many-body state space $s$). By not specifying
the number of flavors/channels above, by default, this assumes a single flavor
(for examples with NC=3 flavors, see further below). The above call to
getLocalSpace thus returns the fermionic annihilation operators
$\hat{F}\equiv(\hat{f}_{\uparrow},\hat{f}_{\downarrow})^{T}$ in irop index
order convention,

⬇

F = ?\mlbl{F1:start}?

Q: 3x { 2 x 2 } having ’A,SU2’ operator, { +– } ?\mlbl{F1:hdr1}?

data: 3-D double (224) 2 x 2 x 1 => 3 x 3 x 2 @ norm = 2 ?\mlbl{F1:hdr2}?

\par1. 1x1 | 1x2x2 [ -1 0 ; 0 1 ; -1 1 ] -1.41421 ?\mlbl{F1:rec1}?

2. 1x1 | 2x1x2 [ 0 1 ; 1 0 ; -1 1 ] -1. ?\mlbl{F1:rec2}? {?$\surd$?2}

the fermionic parity operator $\hat{Z}=(-1)^{\hat{n}}$, with $\hat{n}$ the
particle number operator [cf. Eq. (LABEL:eq:nloc)],

⬇

Z = ?\mlbl{Z:start}?

Q: 2x { 3 x 2 } having ’A,SU2’ { +- }

data: 2-D double (336) 3 x 3 => 4 x 4 @ norm = 2

\par1. 1x1 | 1x1 [ -1 0 ; -1 0 ] 1\.

2. 1x1 | 2x2 [ 0 1 ; 0 1 ] -1. {?$\surd$?2}

3. 1x1 | 1x1 [ 1 0 ; 1 0 ] 1\.

the spin irop as in Eq. (10) with
$\hat{S}_{a}=\sum_{\sigma\sigma^{\prime}}\hat{f}_{\sigma}^{\dagger}S_{a}^{\sigma\sigma^{\prime}}\hat{f}_{\sigma}^{\,}$
and $S_{\pm}\equiv S_{x}\pm\mathrm{i}\mkern 1.0muS_{y}$,

⬇

S = ?\mlbl{idx:Sop:SU2}?

Q: 3x { 1 x 2 } having ’A,SU2’ operator, { +– }

data: 3-D double (112) 1 x 1 x 1 => 2 x 2 x 3 @ norm = 1.225
?\mlbl{idx:Sop:SU2:hdr}?

\par1. 1x1 | 2x2x3 [ 0 1 ; 0 1 ; 0 2 ] -0.866025 {?$\surd$?2} ?\mlbl{Sop:SU2:1}?

and the structure IS with additional information,

x

⬇

IS = ?\mlbl{idx:NC1:IS}?

NC: 1

sym: ’Acharge,SU2spin’

E: [1x1 QSpace]

…

The latter contains the identity operator (IS.E) as a reference for the
complete local state space, a reminder of the number of flavors (I.NC), as
well as more verbose specification of the symmetries used (IS.sym) for later
idenfication. The remaining fields in IS are mostly for debugging purposes,
partly also duplicate data such as the parity Z, and hence can be safely
ignored. The above QSpace displays for F, Z, and S are simply obtained by
typing the respective name followed by enter, e.g., >> F, where the leading
'>>' indicates the Matlab prompt. The generated output conveniently summarizes
the content of the specified QSpace tensor (technically, by overloading the
display routine). The printed output is discussed in detail with Fig. 7.

In the case of the fermionic annihilation operator F, line LABEL:F1:hdr1
identifies the symmetries used in their specified order [A for abelian U(1),
and SU2 for SU(2) spin]. It also shows that F is an operator (see info.otype)
with standard index order +--. The symmetry labels in lines
LABEL:F1:rec1–LABEL:F1:rec2 consist of two numbers $q=[q_{1}\ q_{2}]$
specified for each leg, and separated by a semicolon. The first label $q_{1}$
encodes the particle number (filling) relative to half-filling. Here
$q_{1}\in[-1,0,1]$ on the first two legs which corresponds to $[0,1,2]$
particles occupying the fermionic level. The second label specifies the SU(2)
symmetry labels, given by the integers $q_{2}=2S$. The irop F transforms in
the labels of the third leg, namely $q^{\rm op}{=}$ [-1 1]: that is, the
operator F reduces particle number by one (hence $q_{1}=$ -1), and represents
a spinor of two components (with $q_{2}=2S=$ 1 it transforms like a spin
$S=1/2$ representing up and down spin). All of this is as expected for the set
of fermionic annihilation operators
$\hat{F}\equiv(\hat{f}_{\uparrow},\hat{f}_{\downarrow})^{T}$. Since the irop
index represents a single multiplet (line LABEL:F1:hdr2), the operator
$q$-labels are the same for all records. Hence the operator is irreducible,
i.e., it is an irop, indeed. For comparison, a pair-annihilation operator
$\hat{f}_{\uparrow}\hat{f}_{\downarrow}$ would have $q^{\rm op}_{1}=-2$.
Therefore by simply looking at the action of an operator based on its irop
symmetry labels allows one to (roughly) identify an operator.

By inspection, the first record for F (line LABEL:F1:rec1) describes the
matrix elements $\langle 0|\hat{f}_{\sigma}|\sigma\rangle$ from singly
occupied to empty, with the initial state on leg 2 and the final state on leg
1 (ket and bra state for matrix element, respectively). The corresponding CGT
has the printed symmetry labels $(0|1,1)$, bearing in mind $q=2S=1$ which
happens to correspond to a $1j$ symbol for the first record. It is of size
1x2x2. Conversely, the second record in line LABEL:F1:rec2 describes the
matrix elements from doubly to singly occupied state,
$\langle-\sigma|\hat{f}_{\sigma}|{\uparrow\downarrow}\rangle$ where $-\sigma$
denotes the spin opposite to $\sigma\in\\{\uparrow,\downarrow\\}$. Its CGT
hast the $q$-labels $(1|0,1)$ and dimensions 2x1x2. Therefore the CGT for the
second record corresponds to an identity matrix, when ignoring the singleton
on the second index. The respective reduced matrix elements (RMTs) for the
operator F turn out to be plain numbers, thus of size
1x1$\,[{\color[rgb]{0.45,0.45,0.45}\texttt{x1}}]$ (left column), skipping the
last trailing singleton since beyond rank 2). Their value is directly
displayed to the right. The signs (similar to the sign in line LABEL:Sop:SU2:1
for S) are determined by the sign convention of the underlying CGT for the
irop. They are of no further concern, otherwise. What is important, though, is
that the operator F is defined once and for all, so that the same object is
used for all sites that are identical.

The way getLocalSpace proceeds to obtain the QSpace representations above, is
to start from the Fock space description. For a single spinless fermionic
level, in the basis (0,1) for empty and occupied, the annihilation operator
$c$, fermionic parity operator $z$, and identity operator $e$ are given by

$\displaystyle c$ $\displaystyle=$ $\displaystyle\begin{pmatrix}0&1\\\
0&0\end{pmatrix}\ ,\qquad z=[c,c^{\dagger}]=\begin{pmatrix}1&0\\\
0&-1\end{pmatrix}\ ,\qquad e=\\{c,c^{\dagger}\\}=\begin{pmatrix}1&0\\\
0&1\end{pmatrix}$ (46)

with the sole non-zero matrix element of the annihilation operator being
$\langle 0|\hat{c}|1\rangle=1$. Hence when combining two fermionic levels with
different spin, for example, the matrix representation for the fermionic
annihilation operators can be written by the $4\times 4$ matrices (assuming
the first factor in the tensor product to be the fast index, consistent with
column major),

$\displaystyle f_{\uparrow}=c\otimes e=\begin{pmatrix}0&1&0&0\\\ 0&0&0&0\\\
0&0&0&1\\\ 0&0&0&0\end{pmatrix}\ ,\qquad f_{\downarrow}=z\otimes
c=\begin{pmatrix}0&0&1&0\\\ 0&0&0&-1\\\ 0&0&0&0\\\ 0&0&0&0\end{pmatrix}$ (47)

having assumed the fermionic order
$|{\uparrow\downarrow}\rangle\equiv\hat{c}^{\dagger}_{\uparrow}\hat{c}^{\dagger}_{\downarrow}|\rangle$
thus already leading to the fermionic sign via $z$ for $f_{\downarrow}$ above.
It is important for the fermionic anticommutator relations which must be
already built-in in the above matrix elements, e.g.,

$\displaystyle\\{f_{\uparrow},f_{\downarrow}\\}$ $\displaystyle=$
$\displaystyle(c\otimes e)(z\otimes c)+(z\otimes c)(c\otimes
e)=\\{c,z\\}\otimes c=0.$

since $cz=-zc$ by the very property of $c$ that it changes the particle number
by one. Combining the annihilation operators into the set
$F=(f_{\uparrow},f_{\downarrow})^{T}$, then for example

$\displaystyle N=F^{\dagger}\cdot F\ =\
f_{\uparrow}^{\dagger}f_{\uparrow}^{\,}+f_{\downarrow}^{\dagger}f_{\downarrow}^{\,}\
=\ (c^{\dagger}c\otimes e)+(c^{\dagger}c\otimes z^{\dagger}z)\ =\ (n\otimes
e)+(e\otimes n)\ \equiv\ N_{\uparrow}+N_{\downarrow}$

as expected, with $n=c^{\dagger}c$ the occupation operator. QSpace takes the
combined set F=(f↑,f↓)T above, organizes the state space based on symmetry
operators (Class/@SymOp), and then ‘compresses’ it via compactQS. This
diligently factorizes the CGTs as present in the RC_STORE, with the remainder
being in the RMTs.

##### Explicit data structure

By typing F at the Matlab prompt, this invokes @QSpace/display to conveniently
summarize the specified QSpace data, e.g., as shown for F in lines
LABEL:F1:start–LABEL:F1:rec2. The bare QSpace data structure as in App. A.3 is
obtained by stripping the tensor F of its QSpace class assignment via
⬇ >> struct(F)  x

⬇

struct with fields:

\parQ: {[2x2 double] [2x2 double] [2x2 double]}

data: {2x1 cell}

info: [1x1 struct]

which is considerably less readable and informative as compared to lines
LABEL:F1:start–LABEL:F1:rec2 above, despite referring to the same bare data.
The two RMTs $\texttt{i}=1,2$ in lines LABEL:F1:rec1–LABEL:F1:rec2 are stored
in F.data{i}

⬇

>> F.data’

\par1x2 cell array

\par{[-1.4142]} {[-1.0000]}

with the respective symmetry labels stored in F.Q{l}{i,:} for leg $l$,

⬇

>> F.Q{1}

-1 0

0 1

\par>> F.Q{2}

0 1

1 0

\par>> F.Q{3}

-1 1

-1 1

These are displayed as records i with lines LABEL:F1:rec1–LABEL:F1:rec2
above555As a Matlab technicality, note that F.Q{:} would not work here to
display all entries; rather one would have to use an intermediate assignment,
as in Q=F.Q; Q{:}.. The info structure

⬇

>> F.info ?\mlbl{F.info}?

qtype: ’A,SU2’

otype: ’operator’

itags: {” ’*’ ’*’}

ctime: 0

cgr: [2x2 struct]

is described in more detail with the data structure in App. A.3. It contains
symmetry-related content, such as the symmetry type F.info.qtype or the CGRs
in cgr. The field itags specifies the itag for each leg. Hence the number of
entries has to match numel(F.Q), or equivalently the rank of the tensor.

The field info.cgr is another structure array with one column per symmetry,
and as many rows as there are records in the QSpace. In case of the operator
F, for example, the entry F.info.cgr(1,2) contains the CGR for the QSpace
record $\texttt{i}=1$ and symmetry $2$, i.e., the SU(2) spin symmetry as per
the present setup in (LABEL:eq:gls:NC1),

⬇

>> F.info.cgr(1,2)

\partype: ’SU2’ qset: [0 1 1] ?\mlbl{F:cgr:qset}? qdir: ’+–’ cid: [1.50e+09
1.50e+09 690057 0 8] size: [1 2 2] ?\mlbl{F:cgr:size}? nnz: 2 cgw: 1 cgt: [1x0
double] \par

The symmetry labels qset in line LABEL:F:cgr:qset are those from the
referenced CGT and reflect the ones in X.Q{l}(1,2) for leg $l$ as listed with
the QSpace display of this record
$\texttt{i}={\color[rgb]{0.45,0.45,0.45}\texttt{1.}}$ in line LABEL:F1:rec1
earlier. The dimensions in line LABEL:F:cgr:size are also displayed there.

By contrast, the CGRs X.info.cgr(:,1) for the leading abelian U(1) symmetry
are all trivial. The underlying CGTs are simply $C=1.$, bearing in mind that
only non-zero CGTs are present in a QSpace tensor, i.e., only ones that are
permissible from a symmetry point of view. For example, the trivial CGR from
the first record reads

⬇

>> F.info.cgr(1,1)

\partype: ’A’

qset: [-1 0 -1]

qdir: ’+–’

cid: [0 0 0 0 1]

size: [1 1 1]

nnz: 0

cgw: 1

cgt: [1x0 double]

It reflects the plain number $1.$ by having size = 1, cgw = 1, nnz is ignored,
no internal IDs cid, etc. The symmetry labels qset = [-1 0 -1] are again
consistent with those in X.Q{l}(1,1) for leg $l$ as displayed in line
LABEL:F1:rec1 earlier.

##### State space from operators

The annihilation operator F is non-hermitian and contains only off-diagonal
symmetry blocks. Its full dimensions are shown as 3 x 3 x 2 in line
LABEL:F1:hdr2 which only references 3 out of the total of 4 states with either
of the first two indices. The reason being, that an annihilation operator
destroys an empty ket state $|0\rangle$, yet also has no matrix elements that
include the fully occupied bra state $\langle{\uparrow\downarrow}|$. However
their union of symmetry sectors as obtained by getIdentity when called on an
operator like F here,
⬇ >> E=getIdentity(F) ?\mlbl{E1:F1}?  x

⬇

Q: 2x { 3 x 2 } having ’A,SU2’ { +- }

data: 2-D double (336) 3 x 3 => 4 x 4 @ norm = 2 ?\mlbl{E1:hdr2}?

\par1. 1x1 | 1x1 [ -1 0 ; -1 0 ] 1\.

2. 1x1 | 2x2 [ 0 1 ; 0 1 ] 1. {?$\surd$?2}

3. 1x1 | 1x1 [ 1 0 ; 1 0 ] 1\.

\par

recovers the full 4 x 4 dimensional local state space (line LABEL:E1:hdr2) of
the single spinful fermionic level under consideration. By not having
specified an index with the operator F in line LABEL:E1:F1, the default
behavior of getIdentity assumes an operator as input with standard index
order, and hence takes the union of the state spaces of the first two indices
of F. The fact that the first two indices of the annihilation operator must
span the full state space of a site is also apparent from the realization that
its fermionic anticommutator relations must yield the identity [for an
explicit sanity check in this regard, see (LABEL:eq:acomm)].

Figure 8: Contraction $n=F^{\dagger}\cdot F$ to obtain the occupation operator
in some (local) state space $s$, having spin flavor
$\sigma\in\\{\uparrow,\downarrow\\}$. (b) The contraction in (LABEL:eq:nloc)
takes the conjugate tensor of F reverts all arrows and also reflects the
tensor, here vertically as indicated with the gray line and rotating arrow.
Hence the index order in $1$ and $2$ is already flipped up/down with
$F^{\ast}$. By contracting indices '13*' on the first, i.e., lower tensor, the
1 takes the transpose, thus completing the Hermitian conjugate. The remaining
open indices are collected (c), thus giving rise to the new index order in (d)
for the occupation operator. Uncontracted, and thus open indices are collected
in the order of the input QSpace tensors to contract. Therefore the chosen
order of input arguments in the command Eq. (LABEL:eq:nloc) already returns
the standard index order for operators, as indicated by the orange labels in
(d).

##### Occupation number operator

The operators returned by getLocalSpace in (LABEL:eq:gls:NC1) are considered
elementary for the description of the specified site. Any other local operator
can be constructed from these, usually by simple one-liners. By construction,
these operators are then already consistent with the chosen symmetries. For
example, the local occupation number operator is obtained by
⬇ >> nloc = contract(F,’13*’,F,’13’)  x

⬇

Q: 2x { 2 x 2 } having ’A,SU2’ { +- }

data: 2-D double (224) 2 x 2 => 3 x 3 @ norm = ?$\surd$?6

\par1. 1x1 | 2x2 [ 0 1 ; 0 1 ] 1. ?\mlbl{nloc:rec1}? {?$\surd$?2}

2. 1x1 | 1x1 [ 1 0 ; 1 0 ] 2. ?\mlbl{nloc:rec2}?

with a pictorial representation presented in Fig. 8. In addition to the
contraction of index 1 from the matrix product, the sum over spin corresponds
to the simultaneous contraction of the irop index 3 as specified by the
strings '13*' or '13' in compact notation for the first and second operator,
respectively. The above display for nloc shows two block entries corresponding
to RMTs with 1. or 2. particles as printed to the right in lines
LABEL:nloc:rec1 and LABEL:nloc:rec2, respectively. These values reflect
nloc.data{i}. There is no block record for 0. particles, as this would
represent a zero-block (more precisely, since block-diagonal zero blocks are
usually kept in rank-2 operators, in the present case this block simply does
not arise out of the contraction in the first place). Note that for the value
1. to appear in line LABEL:nloc:rec1, the normalization of the $w$-matrix had
been tweaked, as indicated by the $\\{\surd 2\\}$ factor to the very right
[cf. discussion with Eq. (24)].

The definition of nloc permits one to relate its trace ${\rm tr}(n_{\rm
loc})={\rm tr}(F^{\dagger}\cdot F)=\|F\|^{2}$ with the norm of
$\|{\color[rgb]{0.45,0.45,0.45}\texttt{F}}\|=2$ in line LABEL:F1:hdr2. This
can be quickly double-checked as follows,

⬇

>> [ trace(nloc), norm(F)^2 ]

\par4.0000 4.0000

\par

The trace above includes a factor of $2$ for the sector q=[0 1] in line
LABEL:nloc:rec1, since by representing a spin half ($2S=1$), this contains two
states. Hence ${\rm tr}(n_{\rm loc})=0+2\times 1+2=4$. The fermionic anti-
commutator relations can be checked by
⬇ >> E=acomm(F,F’)/2  x

⬇

Q: 2x { 3 x 2 } having ’A,SU2’ { +- }

data: 2-D double (336) 3 x 3 => 4 x 4 @ norm = 2

\par1. 1x1 | 1x1 [ -1 0 ; -1 0 ] 1\.

2. 1x1 | 2x2 [ 0 1 ; 0 1 ] 1. {?$\surd$?2}

3. 1x1 | 1x1 [ 1 0 ; 1 0 ] 1\.

\par>> isIdentityQS(E) ?\mlbl{E:isId}?

\par1

This needs to sum over $\sigma\in\\{\uparrow,\downarrow\\}$ since F exists as
an irop, where one cannot single out a particular component $\sigma$, as this
would break the symmetry. Line LABEL:E:isId checks whether E represents an
identity operator (solely within the symmetry space present in the QSpace E,
i.e., without knowledge of any potentially ‘absent’ zero blocks; this is no
issue here, though). The call to the QSpace routine acomm() above for the
anti-commutator is equivalent to performing the contraction in (LABEL:eq:nloc)
plus its respective counterpart to complete the anti-commutator [see comment
with (LABEL:eq:acomm)],

⬇

\par>> E_ = ( contract(F,’13*’,F,’13’ ) + …

contract(F,’23’, F,’23*’) ) / 2;

>> norm(E-E_) ?\mlbl{E1:diff}?

\par0

\par>> E==E_ ?\mlbl{E1:same}?

\par1

Lines LABEL:E1:diff and LABEL:E1:same verify that E and E_ are the same
tensors, indeed (line LABEL:E1:diff is somewhat more flexible as compared to
the overload in @QSpace/eq.m, as one may use the former to accept differences
on the order of double precision noise in a general context). Similar to the
fermionic occupation operator in Eq. (LABEL:eq:nloc), one can construct the
local Casimir operator from the spin operator,
⬇ >> S2=contract(S,’13*’,S,’13’)  x

⬇

Q: 2x { 1 x 2 } having ’A,SU2’ { +- }

data: 2-D double (112) 1 x 1 => 2 x 2 @ norm = 1.061 ?\mlbl{S2:hdr2}?

\par1. 1x1 | 2x2 [ 0 1 ; 0 1 ] 0.75 {?$\surd$?2}

Since the spin operator (10) is non-hermitian, the dagger is relevant. The
above display shows a single entry for the half-filled symmetry sector [0 1]
only, which being spin-half ($2S=1$), is the only one with a non-zero Casimir,
indeed, having the value 0.75 shown to the right. All other symmetry sectors,
namely empty and doubly occupied, have $S^{2}=0$. Thus being zero-blocks,
these are absent. Matter of fact, these are already absent in the spin
operator S that entered the contraction. The Casimir is thus already also
reflected in the RMTs of the spin operator, e.g., having the RMT
$\sqrt{3/4}=0.866025$ in line LABEL:Sop:SU2:1 up to a sign. The norm
$\|S^{2}\|$ is given by $\frac{3}{4}\|\ 11^{(2)}\|=\sqrt{9/8}=1.061$ as seen
in line LABEL:S2:hdr2. This is in contrast to $\|S\|=\sqrt{3/4}\
\|11^{(2)}\|=\sqrt{3/2}=1.225$ in line LABEL:idx:Sop:SU2:hdr above.
Based on the set of QSpace operators introduced above, the following proceeds
with explicit examples for the main text earlier. This starts with examples
for Sec. 2.1.

###  4.1 QSpace tensors of rank $r<2$

##### Rank-0 tensor

A tensor of rank $r=0$ has no open legs, and hence can only represent a scalar
number. It naturally occurs when all indices have been contracted, e.g., when
computing expectation values or overlaps. For example, computing the Frobenius
norm of a tensor translates to contracting a tensor fully with itself,

$\displaystyle\|\hat{f}\|^{2}$ $\displaystyle\equiv$ $\displaystyle{\rm
tr}(\hat{f^{\dagger}}\cdot\hat{f})=\sum_{\sigma}{\rm
tr}(\hat{f}_{\sigma}^{\dagger}\cdot\hat{f}_{\sigma})$

the result has no more open legs. In QSpace, this becomes

⬇

>> X=contract(F,’123*’,F,’123’)

x

⬇

Q: [] having ’A,SU2’ double scalar ?\mlbl{rank0:1}?

data: { 4 } ?\mlbl{rank0:2}?

\par>> getscalar(X) ?\mlbl{X:scalar}?

\par4.0000

Therefore it has empty $q$-labels X.Q, and only a single RMT X.data{1} [cf.
App. A.3] with value X.data{1} $=\|F\|^{2}=4$ (cf. line LABEL:F1:hdr2). For
such rank-0 tensors, the QSpace display switches to the simplified format in
lines LABEL:rank0:1–LABEL:rank0:2 above. The scalar value can be extracted by
getscalar (line LABEL:X:scalar) which also performs consistency checks (e.g.,
if X does not represent a scalar, it throws an error).

##### Rank-1 tensor

The only permissible rank-1 tensor is a tensor that carries the symmetry
sector of the vacuum state, namely $q=0$. For a general non-empty QSpace X of
rank $r\geq 2$, the rank-1 tensor representing the vacuum state can be
obtained by
⬇ >> V=getvac(F,’-1d’)  x

⬇

Q: 1x { 1 x 2 } having ’A,SU2’ { + }

data: 1-D double (112) 1 x 1 @ norm = 1

\par1. 1 | 1 [ 0 0 ] 1. ?\mlbl{V:rec1}?

With the option '-1d', this returns a vector of length $d=1$ (line
LABEL:V:rec1). Rank-1 tensors in $q=0$ with dimension $d>1$ 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$. Without the '-1d'
option, getvac in Eq. (LABEL:eq:rank1) returns a rank-2 tensor that represents
the identity operator in the vacuum state, and thus is of dimension $1\times
1$.

###  4.2 Permutations and index order

QSpace tensors support an arbitrary index order that can be permuted as
desired. While some tensors such as operators or $A$-tensors assume certain
index order conventions for coding purposes, these are relevant for
applications only, but not for QSpace tensors per se. For the sake of the
argument, consider an operator that describes fermionic hopping from site s02
and s01, i.e., $\hat{f}_{1}^{\dagger}\cdot\hat{f}_{2}\ =\
\sum_{\sigma}\hat{f}_{1\sigma}^{\dagger}\cdot\hat{f}_{2\sigma}$ with
$\sigma\in\\{\uparrow,\downarrow\\}$,
⬇ >> F2=contract(F,’-op:s01’,’*’, Z*F,’-op:s02’, [2 3 1 4]) ?\mlbl{F2:op}?  x

⬇

Q: 4x { 4 x 2 } having ’A,SU2’ { s01, s02, s01*, s02* } ?\mlbl{F2:hdr1}?

data: 4-D double (448) 2 x 2 x 2 x 2 => 3 x 3 x 3 x 3 @ norm = ?$\surd$?8

\par1. 1x1 | 2x1x1x2 [ 0 1 ; -1 0 ; -1 0 ; 0 1 ] 1.41421

2. 1x1 | 2x2x1x1 [ 0 1 ; 0 1 ; -1 0 ; 1 0 ] -1.41421

3. 1x1 | 1x1x2x2 [ 1 0 ; -1 0 ; 0 1 ; 0 1 ] -1.41421

4. 1x1 | 1x2x2x1 [ 1 0 ; 0 1 ; 0 1 ; 1 0 ] -1.41421

written as rank-4 tensor here by only contracting the operator index via the
dot product in the operators, thus having $\langle
s_{1}s_{2}|\hat{f}_{1}^{\dagger}\cdot\hat{f}_{2}|s^{\prime}_{1}s^{\prime}_{2}\rangle$.
By specifying a conjugate flag '*' with the first operator F, this takes its
Hermitian conjugate. The option '-op:s01' assigns operator itags for site
's01' on the fly. Specifically, this assigns the temporary itags
{'s01','s01*','op*'} to the preceding input operator [cf. Fig. 1(b)].
Similarly, the second operator Z*F gets assigned the operators itags for site
's02' (because in a simulation typically more than 10 sites are included,
QSpace frequently uses a 2-digit format for site or bond indices in itags,
like 's01' and 's02' above for the sake of aligned displays).

Fermionic signs are applied by contracting the parity operator
Z$\equiv(-1)^{\hat{n}}$ from the left onto the operator F acting on site $2$.
This assumes a fermionic order where site $1$ comes before site $2$ when
building a Fock state space, e.g.,
$\hat{f}^{\dagger}_{2}\hat{f}_{1}^{\dagger}|\rangle$ similar to Eq. (47).
Hence the matrix elements for $\hat{f}_{1}$ acquire fermionic signs depending
on the state of site $2$ which thus gets the parity operator applied (see also
appendix in [50] for more discussion on fermionic signs). Since both, Z and F,
are operators, one can simply use the operator * for matrix multiplication in
Z*F.666The matrix multiplication operator * is overloaded via @QSpace/mtimes.m
to perform a QSpace contraction; this also deals with potential rank-3
operators as input, like F here. As a safety measure, this requires the
explicit marking of the QSpace as an ‘operator’ via the flag
F.info.otype='operator'. Here this is automatically set via getLocalSpace,
e.g., see line LABEL:F1:hdr1 above. The default operator itag assigned above
by the option '-op:..' is 'op'. In the present example, this is the only
matching itag in the auto-contraction. Therefore the contraction results in a
rank-4 tensor that explicitly acts on the two separate sites s01 and s02.

Command (LABEL:eq:FdotF) also includes a trailing permutation that is applied
to the overall result. Therefore it is equivalent to
⬇ >> F2=contract(F,’-op:s01’,’*’,Z*F,’-op:s02’); ?\mlbl{F2:Z}?  >>
F2=permute(F2,’2314’)  x
The contraction in line LABEL:F2:Z collects the uncontracted legs of the first
operator conj(F) [note the '*' flag in the contraction], followed by the
uncontracted ones from the second operator, maintaining their original order.
The bra and ket indices of the rank-4 operator F2 out of the contraction can
therefore be properly grouped by the permutation [2 3 1 4] (line
LABEL:F2:perm, or trailing option in line LABEL:F2:op) with the effect
$(s^{\ast}_{1},s_{1},s_{2},s^{\ast}_{2})\to(s_{1},s_{2},s^{\ast}_{1}s^{\ast}_{2})$.
This finalizes the transpose for the first operator, which was conjugated by
the flag '*' in the contraction. With this, the itags of the final output for
F2 have $q$-directions ++-- as seen from the trailing conjugate flags '*' in
the itags in line LABEL:F2:hdr1.

The permutation in line LABEL:F2:perm accepts compact string notation aside
from the plain numerical representation. By contrast, line LABEL:F2:op insists
on a numerical representation of the permutation, in order to differentiate it
from other string options meant for contractions. Line LABEL:F2:perm accepts
any valid permutation to be applied to F2, including an optional trailing
conjugation flag '*'. For example, the following swaps bra with ket indices
and also applies overall tensor conjugation, i.e., obtains the Hermitian
conjugate $\langle
s_{1}s_{2}|\hat{f}_{2}^{\dagger}\cdot\hat{f}_{1}|s^{\prime}_{1}s^{\prime}_{2}\rangle$
for given rank-4 tensor,

⬇

>> F2_ = permute(F2,’3412*’)

\parQ: 4x { 4 x 2 } having ’A,SU2’ { s01, s02, s01*, s02* }

data: 4-D double (448) 2 x 2 x 2 x 2 => 3 x 3 x 3 x 3 @ norm = ?$\surd$?8

\par1. 1x1 | 1x2x2x1 [ -1 0 ; 0 1 ; 0 1 ; -1 0 ] 1.41421

2. 1x1 | 1x1x2x2 [ -1 0 ; 1 0 ; 0 1 ; 0 1 ] -1.41421

3. 1x1 | 2x2x1x1 [ 0 1 ; 0 1 ; 1 0 ; -1 0 ] -1.41421

4. 1x1 | 2x1x1x2 [ 0 1 ; 1 0 ; 1 0 ; 0 1 ] -1.41421

\par

While the itags of the resulting F2_ are identical to the ones from F2, their
content is clearly different, as can be simply checked by norm(F2-F2_) which
results in 4.

###  4.3 $1j$ tensor and trailing marker characters

The identity operator in Eq. (LABEL:eq:Id:F) is a regular operator with one
in- and one out-going index. Instead of the identity, one may request the $1j$
tensor by the same call to getIdentity, yet adding the option '-0' (‘fusing to
$q=0$’),
⬇ >> U=getIdentity(F,’-0’) ?\mlbl{U:1j}?  x

⬇

\parQ: 2x { 3 x 2 } having ’A,SU2’ { +?\textcolor{mcol}{+}? }
?\mlbl{U:1j:hdr1}?

data: 2-D double (336) 3 x 3 => 4 x 4 @ norm = 2 ?\mlbl{U:1j:hdr2}?

\par1. 1x1 | 1x1 [ -1 0 ; 1 0 ] 1\.

2. 1x1 | 2x2 [ 0 1 ; 0 1 ] 1. {?$\surd$?2}

3. 1x1 | 1x1 [ 1 0 ; -1 0 ] 1. ?\mlbl{U:1j:end}?

This shows a colored itag in line LABEL:U:1j:hdr1 when using the Matlab
terminal mode, for readability. The actual 'U.info.itags' are {'','’'}, with
the first itag empty, and the second itag just a prime (’). Since no string
delimiters are shown with itags in the QSpace display, empty strings would not
display at all. Hence, in case of empty strings, the $q$-directions +/- are
displayed, instead. In the presence of a trailing prime marker character, the
itag is colored as above. The color coding is understood by terminals only,
but not the Matlab desktop environment. For the latter by default, line
LABEL:U:1j:hdr1 above would display as

⬇

Q: 2x { 3 x 2 } having ’A,SU2’ { ++’ }

instead, which shows the marker character with the second itag, ++' $\equiv$
{+,+'}. The same can also be enforced by setting setenv QS_LOG_COLOR 0 in the
environment, e.g., on the Matlab prompt itself. If getIdentity is called with
reference to a tensor leg that carries an itag, then that itag is inherited,

⬇

>> E=setitags(getIdentity(F),{’s01’,’s01’});

>> U=getIdentity(E,’-0’)

\parQ: 2x { 3 x 2 } having ’A,SU2’ { s01, ?\textcolor{mcol}{s01}? }
?\mlbl{U:1j:hdr1b}?

…

\par

again, with colored output, as in Matlab terminal mode. The remainder (...) is
the same as in lines LABEL:U:1j:hdr2–LABEL:U:1j:end. The actual itags
U.info.itags are {'s01', 's01’'}, as displayed in the Matlab desktop, where
line LABEL:U:1j:hdr1b would read

⬇

Q: 2x { 3 x 2 } having ’A,SU2’ { s01, s01’ } ?\mlbl{U:1j:hdr1c}?

If one were to take the conjugate tensor of $U$, this results in all outgoing
legs, with the second leg still also marked, which in terminal mode would
translate to the color coding for the sake of readability,

⬇

>> conj(U)

…

Q: 2x { 3 x 2 } having ’A,SU2’ { ?\textcolor{ccol}{s02}?,
?\textcolor{mccl}{s02}? }

…

where the an outgoing itag is printed in gray, and if marked, the original
brighter green in line LABEL:U:1j:hdr1b becomes a darker green. In the Matlab
desktop, the itags are fully spelled out

⬇

>> conj(U)

…

Q: 2x { 3 x 2 } having ’A,SU2’ { s02*, s02’* }

…

The trailing marker toggle inserted by getIdentity with the option '-0' may be
disabled by using the option '-z', instead (z as in zero in '-0' for
historical reasons; not recommended),

⬇

>> Uz=getIdentity(E,’-z’)

\parQ: 2x { 3 x 2 } having ’A,SU2’ { s01, s01 } ?\mlbl{Uz:hdr1}?

data: 2-D double (336) 3 x 3 => 4 x 4 @ norm = 2

\par1. 1x1 | 1x1 [ -1 0 ; 1 0 ] 1\.

2. 1x1 | 2x2 [ 0 1 ; 0 1 ] 1. {?$\surd$?2}

3. 1x1 | 1x1 [ 1 0 ; -1 0 ] 1\.

\par

In this case, the itags (line LABEL:Uz:hdr1) are the same for both legs and
become indistinguishable from the ones for the transpose of the tensor.
However, taking the transpose of Uz,

⬇

>> Uz_=permute(Uz,’21’)

\parQ: 2x { 3 x 2 } having ’A,SU2’ { s01, s01 }

data: 2-D double (336) 3 x 3 => 4 x 4 @ norm = 2

\par1. 1x1 | 1x1 [ 1 0 ; -1 0 ] 1\.

2. 1x1 | 2x2 [ 0 1 ; 0 1 ] 1. {-?$\surd$?2}

3. 1x1 | 1x1 [ -1 0 ; 1 0 ] 1\.

\par

give rise to a minus sign with the second entry (with the permutation absorbed
into the $w$-matrix, and hence showing up with the curly bracket to the very
right). Therefore Uz is clearly different from Uz_, having norm(Uz-Uz_) =
$\surd 8$, even though indistinguishable, say, pictorially when including
itags. Both are unitaries, though, and therefore can act as $1j$ symbols.
However, mixing them up when reverting directions of legs, e.g., by inserting
Uz * Uz_${}^{\dagger}\neq 11$ into a line in a tensor network, this is likely
to give wrong results because of sign errors, even though the contraction is
possible on principle grounds. To avoid such sources of errors that may be
difficult to spot, the trailing prime as a marker character was introduced
with QSpace v4. After all, the dual state space is different from its original
one, even if the considered state space may happen to be self-dual in the
sense that it maps onto itself as a whole, as in the present example. Trying
to contract a marked with a non-marked itag then, will result in an error
because such itags are considered different.

###  4.4 Spin-half spin operators

The spin operators for a single spin $S=1/2$ are obtained as in Eq. (39),
⬇ >> [S,IS] = getLocalSpace(’Spin’,1/2,’-A’);  x resulting in

⬇

>> S \par1. (U1) { +– } 3D double 0.7071 224 2x2x1

2. (U1) { +– } 3D double 0.7071 112 1x1x1 operator

3. (U1) { +– } 3D double -0.7071 112 1x1x1 operator

Since S is a QSpace array of more than two entries, the QSpace display
switches to a list view with one-liners for each entry in S that further
summarize the header lines in the detailed display for each S(i) with
$i=1,2,3$. To see the full QSpace display of all entries in the QSpace array
S, one needs to call the underlying display routine together with the option
'-v' that stands for more ‘verbose’,

⬇

>> display(S,’-v’)

\parS(1) =

Q: 3x { 2 x 1 } abelian U(1) { +– }

data: 3-D double (224) 2 x 2 x 1 @ norm = 0.7071

\par1. 1x1 [ 1 ; 1 ; 0 ] 0.5 ?\mlbl{Sop:1}?

2. 1x1 [ -1 ; -1 ; 0 ] -0.5

\parS(2) =

Q: 3x { 1 x 1 } abelian U(1) operator, { +– }

data: 3-D double (112) 1 x 1 x 1 @ norm = 0.7071

\par1. 1x1 [ -1 ; 1 ; -2 ] 0.707107 ?\mlbl{Sop:2}?

\parS(3) =

Q: 3x { 1 x 1 } abelian U(1) operator, { +– }

data: 3-D double (112) 1 x 1 x 1 @ norm = 0.7071

\par1. 1x1 [ 1 ; -1 ; 2 ] -0.707107 ?\mlbl{Sop:3}?

\par

By inspection, this shows that the returned spin operator S has three entries
corresponding to

$\displaystyle{\color[rgb]{0.45,0.45,0.45}\texttt{S}}\,=[\hat{S}_{z},\tfrac{1}{\sqrt{2}}\hat{S}_{-},-\tfrac{1}{\sqrt{2}}\hat{S}_{+}]\text{
.}$ (59)

For algorithmic convenience, this keeps $\hat{S}_{z}$ to the front by
permuting $\hat{S}_{+}$ to the end w.r.t. Eq. (10), such that S(1) is either
$\hat{S}_{z}$ or the full spin irop in case of SU(2) spin. The spin operator S
above is derived from the representation as $S=1$ irop in Eq. (10), here
corresponding to symmetry labels $S_{z}\in\\{0,-1,+1\\}$ or
${\hyperref@@ii[idx:SU2:qlabels]{q_{\rm op}=2S}}\in\\{0,-2,+2\\}$ for the irop
$q$-label with the third index in lines LABEL:Sop:1f, LABEL:Sop:2, and
LABEL:Sop:3, respectively [hats are used with the operators in (59) to
differentiate them from the symmetry labels $|S{=}1;\,S_{z}\in-1,0,1\rangle$
here for the $S=1$ multiplet underlying the irop]. This explains the
normalization and signs in the components of S up to the permutation of the
entries. Note that the inverse square root factors with the lowering and
raising operator lead to the same norm
$|{\color[rgb]{0.45,0.45,0.45}\texttt{S(i)}}|=\frac{1}{\sqrt{2}}=0.7071$
($i=2,3$) as compared to the $S_{z}$ operator ($i=1$). As such, they also
allow for a convenient $S^{\dagger}\cdot S$ contraction for isotropic spin
interactions.

Having U(1) spin, all three operators in S are distinguishable by their irop
symmetry labels. Therefore the three spin operators may be combined into a
single spin operator by simply ‘adding’ them,

⬇

>> Sop=sum(S) ?\mlbl{Sop:sum}?

\parQ: 3x { 4 x 1 } abelian U(1) { +– }

data: 3-D double (448) 2 x 2 x 3 @ norm = 1.225 ?\mlbl{Sop:hdr}?

\par1. 1x1 [ -1 ; -1 ; 0 ] -0.5

2. 1x1 [ -1 ; 1 ; -2 ] 0.707107

3. 1x1 [ 1 ; -1 ; 2 ] -0.707107

4. 1x1 [ 1 ; 1 ; 0 ] 0.5

\par

Because all irop symmetry sectors in S are different, addition translates into
a direct sum $\oplus$ here. This simply catenates all QSpace records into a
single list, since the entries are still differentiated by their U(1) irop
labels, thus preserving the RMT block structure. The resulting Sop is no
longer irreducible from a symmetry perspective, though. Therefore in contrast
to the individual entries S(i), the combined Sop no longer qualifies as an
‘irop’, since it has multiple operator symmetry labels listed with its third
index. It is still a perfectly well-defined operator, nevertheless. The norm
of the spin operator Sop,

$\displaystyle\|S\|^{2}$ $\displaystyle\equiv$ $\displaystyle{\rm
tr}(S^{\dagger}\cdot S)=\sum_{i=1,2,3}{\rm
tr}\,(S_{i}^{\dagger}S_{i})=\sum_{a=x,y,z}{\rm
tr}\,(S_{a}^{\dagger}S_{a})=3\times(2\,\tfrac{1}{2^{2}})\ =\ \tfrac{3}{2}$

agrees with line LABEL:Sop:hdr, having $|S|=\sqrt{3/2}=1.225$. The combined
spin operator Sop also permits one to easily obtain the Casimir
$\hat{S}^{2}\equiv\hat{S}^{\dagger}\cdot\hat{S}$,

⬇

>> S2=contract(Sop,’13*’,Sop,’13’) ?\mlbl{Sop:sum:S2}?

\parQ: 2x { 2 x 1 } abelian U(1) { +- }

data: 2-D double (224) 2 x 2 @ norm = 1.061

\par1. 1x1 [ -1 ; -1 ] 0.75

2. 1x1 [ 1 ; 1 ] 0.75

\par

resulting in $\frac{3}{4}11$ as expected for the underlying spin-half. The
contraction on the third index above implements the dot product in
$\hat{S}^{\dagger}\cdot\hat{S}$.

##### Switch to $Z_{2}$ symmetry

The option '-A' in the setup (LABEL:eq:gls:spinhalf) requested abelian U(1)
symmetry. By specifying '--Z2', instead, this switches to a $Z_{2}$ symmetry
representation [cf. Eq. (41)],

⬇

>> [S,IS] = getLocalSpace(’Spin’,1/2,’–Z2’) ?\mlbl{Sop:Z2}?

⬇

\parS(1) = ?\mlbl{Sop:Z2:1}?

Q: 3x { 2 x 1 } having ’Z2’ operator, { +– }

data: 3-D double (224) 2 x 2 x 1 @ norm = 0.7071

\par1. 1x1 [ 1 ; 1 ; 0 ] -0.5

2. 1x1 [ 0 ; 0 ; 0 ] 0.5

\parS(2) = ?\mlbl{Sop:Z2:2}?

Q: 3x { 1 x 1 } having ’Z2’ operator, { +– }

data: 3-D double (112) 1 x 1 x 1 @ norm = 0.7071

\par1. 1x1 [ 0 ; 1 ; 1 ] 0.707107

\parS(3) = ?\mlbl{Sop:Z2:3}?

Q: 3x { 1 x 1 } having ’Z2’ operator, { +– }

data: 3-D double (112) 1 x 1 x 1 @ norm = 0.7071

\par1. 1x1 [ 1 ; 0 ; 1 ] -0.707107

\par

By comparison to the earlier U(1) spin operator, the setup of the model
changed from $q=2S$ for U(1) spin to a $Z_{2}$ setting with
$q\in\\{0,1\\}\equiv\\{\uparrow,\downarrow\\}$, as seen from $S_{z}$ in lines
LABEL:Sop:Z2:1ff. The spin operator S(1) $=S_{z}$ carries $q_{\rm op}=0$,
whereas the operators $\hat{S}_{+}$ and $\hat{S}_{-}$ become indistinguishable
in terms of irop symmetry labels, both having $q_{\rm op}=1$ due to the
$Z_{2}$ symmetry. However, when including the symmetry labels for all legs,
all entries for the spin operator above are still distinguishable from a
symmetry sector point of view. Hence the earlier U(1) constructions in lines
LABEL:Sop:sum and LABEL:Sop:sum:S2 still also work.

##### Switch to no symmetry

By switching the symmetry option in (LABEL:eq:gls:spinhalf) from '-A' to '--
nosym' (‘no symmetries’), this effectively turns off symmetries [cf. Eq.
(43)],

⬇

>> [S,IS] = getLocalSpace(’Spin’,1/2,’–nosym’) ?\mlbl{Sop:nosym}?

⬇

\parS(1) =

Q: 2x { 1 x 1 } having ’A’ { +- }

data: 2-D double (136) 2 x 2 @ norm = 0.7071

\par1. 2x2 [ 0 ; 0 ] 32 b

\parS(2) =

Q: 2x { 1 x 1 } having ’A’ { +- }

data: 2-D double (136) 2 x 2 @ norm = 0.7071

\par1. 2x2 [ 0 ; 0 ] 32 b

\parS(3) =

Q: 2x { 1 x 1 } having ’A’ { +- }

data: 2-D double (136) 2 x 2 @ norm = 0.7071

\par1. 2x2 [ 0 ; 0 ] 32 b

\par

This is a tweak to abelian U(1), seen as 'A' in the header lines, that only
uses symmetry labels $q=0$. This trivial setting is equivalent to having no
symmetry at all. The RMTs are $2\times 2$ matrizes,

⬇

>> S(1).data{1}

\par0.5000 0

0 -0.5000

\par>> S(2).data{1}

\par0 0

0.7071 0

\par>> S(3).data{1}

\par0 -0.7071

0 0

now reflecting in detail the spin operators in terms of the Pauli matrices,
$[\frac{1}{2}\sigma_{z},\,\frac{1}{\sqrt{2}}\sigma_{-},\,-\frac{1}{\sqrt{2}}\sigma_{+}]$,
respectively. In the present case of no symmetry, there are no symmetry labels
that can describe the action of the operators. Hence all operators are rank-2
scalars. This way, all spin operators become indistinguishable in terms of
irop $q$-label. These spin operators can no longer be combined into a single
spin operator as in line LABEL:Sop:sum, since the RMTs are no longer
distinguishable by symmetry labels. Hence sum(S) would no longer perform the
direct sum ($\oplus$) by catenating the QSpace records, but actually add the
RMTs which is not meaningful.

##  5 Simple Tutorials and Applications

This section extends on the simple examples in Sec. 4 which supplemented the
earlier sections in the main text. Here the focus shifts more toward tensor
network applications. For this purpose, it should be noted that the public
QSpace repository already contains significantly more than the bare tensor
library documented here, e.g., see the listing of (MEX) applications in App.
C. The full documentation of this is left for the future. The standard help
explaining the purpose and usage of any function is available as part of the
repository, nevertheless. Yet their discussion as well as the description of
the related setup scripts is beyond the scope of this documentation which
rather focuses on QSpace as a tensor library. Still, the interested user may
find it rewarding to explore. Since fdm-NRG [19] has been one of the very
first applications of QSpace, it is discussed in more detailed in Sec. 5.6.

Sanity checks for tensors that have a simple interpretation and structure are
also frequently included. These are helpful in practice quite generally, as
they ensure that one has a good understanding of all the operators and tensors
that one is dealing with. Much of the same spirit also underlies the earlier
examples given, and QSpace as a whole. There are many consistency and
plausibility checks internally in the QSpace MEX core routines, as well as in
the Matlab environment. At negligible numerical overhead, these assertions are
safeguards to ensure overall consistency.

###  5.1 Iterative build-up of many-body state spaces: Two fermionic sites

An elementary step to building many-body state spaces is the iterative
addition of a new site. The example here shows how one can start this within
QSpace by adding a second site. To be specific, suppose one wants to describe
two interacting spinful fermionic sites with respective local state spaces
$\sigma_{i}$ and site index $i=1,2$ in a combined state space. The description
of a single site again starts with (LABEL:eq:gls:NC1) above. With this one can
go ahead and start taking iterative tensor products of that state space with
copies of itself to build a many-body basis.

Reminding oneself that QSpace is a tensor library that encodes tensors only,
state spaces are specified by pointing to particular legs of a given tensor.
Hence when expanding a Hilbert space via a tensor product, one needs to ensure
to pick operators that have complete state space, i.e., non-zero blocks in
every symmetry sector of the target state space. Natural choices are the
identity stored with IS.E [cf. Eq. (LABEL:eq:gls:NC1)], or the parity Z. When
building a tensor product of state spaces, the routine getIdentity takes two
QSpaces as input, both representing a state space of their own (here by having
identical sites, the operator Z is listed twice),
⬇ A=getIdentity(Z,Z,[1 3 2]); ?\mlbl{A2:Z}?
A=setitags(A,{’K01’,’K02’,’s02’}); ?\mlbl{A2:itags}?  x
In line LABEL:A2:Z, one may use Z $\to$ IS.E, instead. One even could have
used F instead of Z, since if no explicit leg index is provided with the
input, getIdentity assumes an operator, and hence combines the state spaces
from the first two legs [same argument as with Eq. (LABEL:eq:Id:F) above].
Either way, A encodes the tensor-product fusion of the state spaces of two
sites. The optional last argument [1 3 2] in line LABEL:A2:Z specifies a
permutation to be applied to the returned output $A$-tensor. By default, the
output of getIdentity in line LABEL:A2:Z has an index order that lists the two
input spaces first (in the order specified), and lists the combined state
space last (third index). By permuting this order as specified in line
LABEL:A2:Z, this returns the rank-3 $A$-tensor with index order that adheres
to the LRs index order convention [cf. Fig. 1]. As an implementational detail,
the permutation in line 1 requires numeric format, since a (compact) string
representation such as '132' would be interpreted as an itag for the fused
index, instead.

From an MPS perspective, the above tensor product of sites $\sigma_{1}$ and
$\sigma_{2}$ may represent the start of a physical chain. Virtual bond indices
specify the many-body state spaces that are kept from one iteration to the
next. Hence $\sigma_{1}\to K_{1}\equiv$ 'K01' may be considered the state
space that is kept from iteration 1 (cf. Fig. 2). QSpace, and hence also this
documentation, frequently uses a 2-digit format for site or bond indices for
the sake of aligning display output, bearing in mind that typically more than
10 sites are present in a simulation. Hence 'K01' instead of, e.g., 'K1'. By
adding site $\sigma_{2}\to$'s02', one arrives at the combined state space
'K02'. The itags in line LABEL:A2:itags are thus chosen as in Fig. 2 for
$A_{2}$ at the start of an MPS. With this, the $A$-tensor returned with line
LABEL:A2:itags reads (see Fig. 7 for detailed explanations of this very
example)

⬇

>> A ?\mlbl{A2:disp}?

Q: 3x {10 x 2 } having ’A,SU2’ A-matrix, { K01, K02*, s02 } ?\mlbl{A2:hdr1}?

data: 3-D double (1200) 3 x 10 x 3 => 4 x 16 x 4 @ norm = 4 ?\mlbl{A2:hdr2}?

\par1. 1x1 | 1x1x1 [ -1 0 ; -2 0 ; -1 0 ] 1. ?\mlbl{A2:rec1}?

2. 1x2 | 1x2x2 [ -1 0 ; -1 1 ; 0 1 ] 16 b {?$\surd$?2} ?\mlbl{A2:rec2}?

3. 1x2 | 2x2x1 [ 0 1 ; -1 1 ; -1 0 ] 16 b {?$\surd$?2}

4. 1x3 | 1x1x1 [ -1 0 ; 0 0 ; 1 0 ] 24 b

5. 1x3 | 2x1x2 [ 0 1 ; 0 0 ; 0 1 ] 24 b

6. 1x3 | 1x1x1 [ 1 0 ; 0 0 ; -1 0 ] 24 b

7. 1x1 | 2x3x2 [ 0 1 ; 0 2 ; 0 1 ] 1. {?$\surd$?3} ?\mlbl{A2:rec7}?

8. 1x2 | 2x2x1 [ 0 1 ; 1 1 ; 1 0 ] 16 b {?$\surd$?2}

9. 1x2 | 1x2x2 [ 1 0 ; 1 1 ; 0 1 ] 16 b {?$\surd$?2}

10. 1x1 | 1x1x1 [ 1 0 ; 2 0 ; 1 0 ] 1. ?\mlbl{A2:recl}?

The itags are shown in the header line LABEL:A2:hdr1. With K02 the outgoing
index, it has the trailing conjugate flag '*'. Since $q$-directions need to be
preserved, when setting itags via setitags there is no need to specify the
trailing asterisk with K02 in line LABEL:A2:itags (if trailing asterisks had
been specified, these would be ignored).

The norm in line LABEL:A2:hdr2 shows $\|A\|=4$, i.e., $\|A\|^{2}={\rm
tr}(11)=16$, which is consistent with the full Hilbert space dimension
reflected in the dimensions [4 16 4] of $A$ in the same line. The square root
factors to the very right, as always, are just reminders of the tweaked
normalization with the $w$-matrix [cf. Eq. (25)], so that the reduced matrix
elements read 1., as one may expect for a normalized mapping of state space,
represented as a sliced-up identity matrix in multiplet space. For example, in
line LABEL:A2:rec7, the factor $\sqrt{3}$ [which is due to having the fused
multiplet $q=2S=2$ of dimension $3$] is split off into the $w$-matrix, such
that the RMT can read 1. When the RMT has dimensions $>1$, this shows the size
of that RMT in bytes instead of its value. For example, for the 1x2[x1] RMT in
line LABEL:A2:rec2, 2 * (8 bytes for double) = 16 b. The trailing singleton
dimensions in the RMT on index 3 are not shown (second column).

With the 2-site Hilbert space described by $A$, one can now proceed to compute
matrix elements in the combined many-body state space. Consider, for example,
the spin interaction $S_{12}\equiv A^{\dagger}(S_{1}^{\dagger}\cdot
S_{2}^{\,})A$ when cast into the basis $A$ [see Fig. 10(d) or also Fig. 9 with
$(Z)F\to S$ for pictorial representations]
⬇ S12 = contract(A,’!2*’,{S,’-op:K01’,’*’,{A,S,’-op:^s’}}) ?\mlbl{S12:op}?  x

⬇

Q: 2x { 2 x 2 } having ’A,SU2’ { K02, K02* } ?\mlbl{S12:hdr1}?

data: 2-D double (288) 4 x 4 => 6 x 6 @ norm = 0.866 ?\mlbl{S12:hdr2}?

\par1. 3x3 | 1x1 [ 0 0 ; 0 0 ] 72 b ?\mlbl{S12:rec1}?

2. 1x1 | 3x3 [ 0 2 ; 0 2 ] 0.25 {?$\surd$?3} ?\mlbl{S12:rec2}?

\par>> S12.data{1}

\par0 0 0 ?\mlbl{S12:data1:1}?

0 -0.7500 0 ?\mlbl{S12:data1:2}?

0 0 0 ?\mlbl{S12:data1:3}?

Before explaining the compact semantics of the contraction in Eq.
(LABEL:eq:SdotS), the output in QSpace S12 is examined. As expected for the
operator $S_{1}^{\dagger}\cdot S_{2}^{\,}$, there are only matrix elements at
half-filling in the above display for S12, having $q_{1}=0$ for the first
symmetry label. More restrictive still, non-zero matrix elements can only
derive from the half-filled space on either site. The respective tensor
product of two spin-halfs gives one singlet ($q_{2}=0$) and one triplet state
(${\hyperref@@ii[idx:SU2:qlabels]{q_{2}=2S}}=2$). The triplet state is
represented by the second record in line LABEL:S12:rec2 which has $\langle
S_{1}^{\dagger}\cdot S_{2}^{\,}\rangle=$ 0.25 as expected. The singlet across
both sites is buried within the first record in line LABEL:S12:rec1 where it
appears as the second entry (line LABEL:S12:data1:2) and carries the value
$\langle S_{1}^{\dagger}\cdot S_{2}^{\,}\rangle=$ -0.7500, as expected for the
singlet. In the global half-filled sector, there are two other states in the
scalar sector q=[0 0], namely the combination of completely empty at site 1
with completely filled at site 2, or vice versa. Their expectation value
$\langle S_{1}^{\dagger}\cdot S_{2}^{\,}\rangle=0$, thus corresponding to the
other two zero entries along the diagonal in lines LABEL:S12:data1:1 and
LABEL:S12:data1:3. Because they are part of the same global symmetry sector,
they also need to appear here within the same RMT indexed by $n$ as in Eq.
(2). Here this results in the RMT of dimension $3\times 3$ in lines
LABEL:S12:data1:1–LABEL:S12:data1:3. Since the Hamiltonian preserves spin
symmetry, the matrix elements of S12 must already be in a diagonal
representation, as seen above, indeed.

##### Contraction semantics

The nested contraction in Eq. (LABEL:eq:SdotS) consists of three levels, based
on the pairwise contraction pattern fixed by adding brackets to
$S_{12}=A^{\dagger}\,(S_{1}^{\dagger}\cdot(S_{2}^{\,}\,A))$,

$\displaystyle S_{12}$ $\displaystyle=$ $\displaystyle{\rm
contract}(\underbrace{A,{\color[rgb]{0.45,0.45,0.45}\texttt{\textquotesingle{!2*}\textquotesingle}},\underbrace{\\{S,{\color[rgb]{0.45,0.45,0.45}\texttt{\textquotesingle{-op:K01}\textquotesingle}},{\color[rgb]{0.45,0.45,0.45}\texttt{\textquotesingle{*}\textquotesingle}},\underbrace{\\{A,S,{\color[rgb]{0.45,0.45,0.45}\texttt{\textquotesingle{-op:\textasciicircum
s}\textquotesingle}}\\}}_{\equiv C_{2}}\\}}_{\equiv C_{1}}}_{\equiv C_{0}})$
(62)

[see Fig. 10(d) or also Fig. 9 with $(Z)F\to S$ for pictorial
representations]. The underlying MEX routine contractQS supports nested cell
structures based on pairwise contractions, also referred to as contraction
pattern. With this, one may contract any arbitrary number of tensors in a
single call to contract (the Matlab wrapper routine contractQS). For any
included pairwise contraction, the two involved tensors can be either an
existing QSpace specified as argument, or another cell $\\{\ldots\\}$ thus
adding another recursive ‘level’ to the contraction. The latter indicates a
nested contraction that needs to be carried out first, so that its result can
be used. If further cells are encountered, contractQS proceeds recursively
through the nested calls as they are encountered. Eventually, at the deepest
level of any set of nested calls must be two QSpace objects to be contracted
[like with $C_{2}$ in Eq. (62)] which then is carried out first. Each QSpace
input with empty itags may have an option that assigns itags on the fly,
typically via '-op:..' as above. Each QSpace or cell also accepts an option
that specifies what indices (not) to contract, and whether to take the
conjugate. This option is specified via a compact string with the format
[[!]cidx][*] where terms in square brackets are optional. The term cidx
represents a compact string that explicitly specifies which indices to
contract (or not to contract in the presence of the leading '!'). If cidx is
empty, an implicit auto-contract searches for all matching indices based on
their itags. The trailing * takes the conjugate of the preceding QSpace or
cell. For an explanation of all options available for contractions, see help
contractQS, or equivalently, contractQS -h.

The nested cell structure as provided by the input to contractQS fully
determines the order of contractions. Specifically, QSpace does not optimize
or restructure the order of contractions. The optimal order of contractions is
thus left for the user to determine. In any contraction like $T_{12}\equiv
T_{1}\ast T_{2}$, the non-contracted, i.e., kept indices of the first tensor
$T_{1}$ are collected before the ones in the paired up $T_{2}$, while
maintaining their original order. This holds for any pairwise construction in
the nested structure provided to contractQS. This way, the structure of the
input to contractQS also fully determines the index order of the final result.

Figure 9:  Computing the matrix elements for the fermionic hopping
$\hat{f}_{1}^{\dagger}\cdot\hat{f}_{2}$ in the combined basis $K_{2}$
corresponding to the nested cell contraction in (LABEL:eq:FdotF:T12). The
contraction includes four tensors as shown in (b). By taking the conjugate of
the tensor $F_{1}$, i.e., $F$ acting on $K_{1}\equiv\sigma_{1}$, similar to
Fig. 8 again the up/down reflected mirror image is drawn for $F^{\ast}$. This
is emphasized by the gray line with a rotating arrow and the reverted index
order 1 and 2. Auto-contraction based on itags automatically links the indices
correctly to take the transpose for the Hermitian conjugate $F_{1}^{\dagger}$
overall. Since $F_{2}$, i.e., $F$ applied to $\sigma_{2}$, acts on $|A\rangle$
first [panel (b)], the operator line op leads to a crossing of the line
$\sigma_{2}$ below $F_{2}$, indicated by the blue symbol. Because the operator
line op carries an odd charge parity, this crossing gives rise to a fermionic
swap gate [33]. Here this translates into the parity operator $Z$ (blue)
applied to $\sigma_{2}$ at the location of the blue symbol. This is equivalent
to taking F $\to$ ZF for site 2. If spin operators $S_{1}^{\dagger}\cdot
S_{2}$ had been considered instead of the fermionic operators
$\hat{f}_{1}^{\dagger}\cdot\hat{f}_{2}$ here, the precise location of the line
crossing (above or below $F_{2}\to S_{2}$) would be of no further importance
because the spin operators carry an even charge parity on op (hence they also
commute, $S_{1}^{\dagger}\cdot S_{2}=S_{2}\cdot S_{1}^{\dagger}$). The nested
cell contraction in (LABEL:eq:FdotF:T12) gives the sequence of pairwise
contraction indicated in panels (b-d). The contractions $C_{2}$, $C_{1}$, and
$C_{0}$ highlighted by the gray shaded area are with reference to Eq. (62).
The last contraction $C_{0}$ in panel (d) gives the resulting scalar operator
in the state space $K_{2}$.

Returning to the example in Eq. (62), the inner most cell $C_{2}$ occurs at
nested contraction level 2\. It is contracted first. Since S does not have any
itags assigned yet (after all, it can be applied to any site), they are
assigned on the fly by the provided ‘hint’ whose semantics is along the lines
of a simplified regular expression: by specifying '-op:..', the tensor S is
assumed to be an operator, i.e., of rank 2 or 3, that acts on a particular leg
on the other paired up tensor in the contraction (here A in $C_{2}$). By
looking for an itag that starts with s (denoted by ^s; see C++ std::regex for
more on this syntax), in the present case, this is sufficient to uniquely
identify leg s02 on A. In this sense, the syntax in the contraction internally
assigns on the fly the itags S for this particular contraction. Based on these
itags then, QSpace auto-contracts this spin operator onto s02 of A. The
temporary itags in S that are not contracted, are inherited by and kept with
the temporary tensor $C_{2}$. Given $C_{2}$ as a temporary object in memory,
the next cell to be contracted is $C_{1}$, at nested contraction level 1.
Analogous to $C_{2}$, this now contracts S onto K01. More precisely, because
of the presence of the '*' option, this applies the conjugate, and hence
effectively the Hermitian conjugate and thus the dagger of S onto K01. The
subsequent auto-contraction in $C_{1}$ also contracts the assigned default
operator itag 'op' for both instances of 'S'. Overall, therefore this
contracts $S_{1}^{\dagger}\cdot S_{2}$. The last contraction $C_{0}$ at base
level 0 uses the temporary QSpace tensor $C_{1}$. The result out of $C_{0}$ is
returned and assigned to S12 (technically, $C_{0}$ is already generated in the
memory space of the output S12 to avoid copying objects). For this last
contraction, however, all indices can be paired up, and hence auto-contraction
would contract them all. To avoid this, QSpace introduces the ‘not’ semantics
'!...'. For example, here by specifying '!2..' after the first operator, this
indicates that index 2 of that operator (here referring to the combined state
space K02) shall not be contracted and thus left open, even though the itags
could be paired up, in principle. The trailing *, finally, indicates to also
take the conjugate tensor of A. This corresponds to bra states, and therefore
concludes the calculation of the matrix elements in the combined state space
K02 of sites 1 and 2. Completely analogous to the spin interaction above, one
can also compute the matrix elements for the fermionic hopping across sites 1
and 2, $T_{12}\equiv f_{1}^{\dagger}\cdot f_{2}^{\,}\ \to\
{\color[rgb]{0.45,0.45,0.45}\texttt{F}}_{1}^{\ast}\cdot({\color[rgb]{0.45,0.45,0.45}\texttt{ZF}})_{2}$,
⬇ >> T12 = contract(A,’!2*’,{F,’-op:K01’,’*’,{A,(Z*F),’-op:^s’}});  x
This contraction is the same as in (LABEL:eq:SdotS), except for the operators
acting on the two sites, substituting (S,S) $\to$ (F, Z*F). Here by having
fermionic operators, one needs to incorporate fermionic signs [50, 33]. As a
short remnant from a Jordan Wigner string for the nearest-neighbor hopping
here, the fermionic parity operator Z is applied onto the second F operator.
In terms of fermionic order, this adopts the convention that the first site
acts on the vacuum state first, followed by the second site, etc. Hence
similar to (LABEL:eq:FdotF), as $f_{1}$ needs to be pulled past any
$f_{2}^{\dagger}$ when computing matrix elements in Fock space, it leaves a
trail of parity operators (here to be applied onto $f_{2}$ from the left). In
the pictorial representation of Fig. 9, the crossing [50, 33] of the operator
line op with $\sigma_{2}$ is precisely the location where the parity gate Z
needs to be applied if the local operator has a fermionic character. This
results in S $\to$ ZF for the operator acting on site $\sigma_{2}$. Adding the
Hermitian conjugate to $T_{12}$, one obtains a typical Hamiltonian term

⬇

>> H12 = T12 + T12’

\parQ: 2x { 3 x 2 } having ’A,SU2’ { K02, K02* } ?\mlbl{H12:hdr1}?

data: 2-D double (448) 7 x 7 => 11 x 11 @ norm = 4 ?\mlbl{H12:hdr2}?

\par1. 2x2 | 2x2 [ -1 1 ; -1 1 ] 32 b {?$\surd$?2} ?\mlbl{H12:rec1}?

2. 3x3 | 1x1 [ 0 0 ; 0 0 ] 72 b ?\mlbl{H12:rec2}?

3. 2x2 | 2x2 [ 1 1 ; 1 1 ] 32 b {?$\surd$?2} ?\mlbl{H12:rec3}?

\par

It is hermitian, as can be explicitly checked by computing norm(H12-H12’)
which yields 0. When looking at the overall dimensions of H12, one realizes a
subtlety: the combined Hilbert space for two spinful fermionic levels has a
dimension of 16, yet H12 shows a dimension of 11 x 11 in line LABEL:H12:hdr2.
Apparently, 5 out of 16 states went ‘missing’. The only way this can occur is
that these belong to symmetry sectors of H12 that represent zero blocks, i.e.,
have Frobenius norm below or comparable to numerical double precision noise
(it is not possible for individual states to be missing within otherwise
present symmetry sectors). For the efficiency of a tensor library, such zero-
blocks are not stored, by default. They are kept nevertheless, however, in
QSpace for scalar operators (cf. skipZerosQS). With diagonal zero blocks
absent here nevertheless, these sectors simply did not arise out of the
contraction itself. Indeed, this is a consequence of the symmetry sectors that
were present in the input tensors F and the contracted performed. From the
one-particle picture, there are two spin-degenerate single-particle levels for
the Hamiltonian H12 above at energies $\varepsilon_{\pm}=\pm 1$. When building
the many-body state space for this, the ground state represents the doubly-
filled level $\varepsilon_{-}$ at energy $E_{0}=-2$. It corresponds to a half-
filled state and hence belongs to the symmetry sector
$q={\color[rgb]{0.45,0.45,0.45}\texttt{[0 0]}}$ (line LABEL:H12:rec2). Adding
or removing particles then gives the many-body excitation spectrum
$E\in\\{-2,-1,0,1,2\\}$ (e.g., see first column in lines
LABEL:I12.ee:1–LABEL:I12.ee:end below). Here the completely empty or filled
sectors ($[\pm 2\ 0]$, both 1-dimensional) are missing in H12, since indeed,
these have hopping matrix elements equal to zero. By comparing to the symmetry
sectors on A in lines LABEL:A2:rec1–LABEL:A2:recl [or equivalently, but more
readable, by comparing to the output of E2 in line LABEL:E2:A:2 below], one
finds that the remaining symmetry sector that is missing is ([0 2]). One also
quickly realizes why: this is the triplet state at half-filling (having
$S=2/2=1$). Indeed, say with an up-spin on each site, there cannot be spin-
preserving hopping. This adds up to the total of $1+1+3=5$ ‘missing’ states,
or correspondingly, 3 symmetry sectors. They may be added to H12 by an
infinitesimal (with the value zero also permitted) to also include all-zero
eigenvalue blocks. Since H12 has non-diagonal RMTs, eigenvalue decomposition
yields,

⬇

>> E2=getIdentity(A,2); ?\mlbl{E2:A:2}?

[U12,E12,I12] = eig( H12 + 0*E2 ); ?\mlbl{H12:eig}?

Here eig is a wrapper to the MEX routine eigQS. Line LABEL:H12:eig adds the
identity operator of the full state space on leg 2 in A, i.e., K02, with
weight zero. Being a scalar operator, the resulting diagonal zero-blocks are
kept. The reference to the $A$-tensor ensures a complete state space, by
construction. The output of the above can be inspected, as usual, by simply
typing the object’s name,

⬇

>> U12

Q: 2x { 6 x 2 } having ’A,SU2’ { K02, K02* }

data: 2-D double (784) 10 x 10 => 16 x 16 @ norm = 4

\par1. 1x1 | 1x1 [ -2 0 ; -2 0 ] 1\.

2. 2x2 | 2x2 [ -1 1 ; -1 1 ] 32 b {?$\surd$?2}

3. 3x3 | 1x1 [ 0 0 ; 0 0 ] 72 b

4. 1x1 | 3x3 [ 0 2 ; 0 2 ] 1. {?$\surd$?3}

5. 2x2 | 2x2 [ 1 1 ; 1 1 ] 32 b {?$\surd$?2}

6. 1x1 | 1x1 [ 2 0 ; 2 0 ] 1\.

\par>> E12

Q: 2x { 6 x 2 } having ’A,SU2’ { K02, K02* }

data: 2-D double (704) 6 x 10 => 10 x 16 @ norm = 4

\par1. 1x1 | 1x1 [ -2 0 ; -2 0 ] 0\.

2. 1x2 | 2x2 [ -1 1 ; -1 1 ] 16 b {?$\surd$?2}

3. 1x3 | 1x1 [ 0 0 ; 0 0 ] 24 b

4. 1x1 | 3x3 [ 0 2 ; 0 2 ] 0. {?$\surd$?3}

5. 1x2 | 2x2 [ 1 1 ; 1 1 ] 16 b {?$\surd$?2}

6. 1x1 | 1x1 [ 2 0 ; 2 0 ] 0\.

\par>> I12.ee ?\mlbl{I12.ee}?

\par-2.0000 1.0000 ?\mlbl{I12.ee:1}?

-1.0000 2.0000

-1.0000 2.0000

0.0000 1.0000

0.0000 1.0000

0.0000 3.0000

0.0000 1.0000

1.0000 2.0000

1.0000 2.0000

2.0000 1.0000 ?\mlbl{I12.ee:end}?

\par

The columns of U12 encode the eigenstates that bring H12 into diagonal form
E12. The latter is stored in compact diagonal format (storing the diagonal
only for each RMT as a row vector, as seen by the RMT dimensions in the second
column). The compact diagonal format can be expanded to a regular diagonal
matrix via diag(E12). As seen from E12, there are three symmetry sectors that
have zero blocks (all of them containing a single multiplet which thus prints
the 0. to the right). The full eigenspectrum is provided in plain numeric
format with the info structure IS returned as third argument. The first column
in I12.ee in lines LABEL:I12.ee:1–LABEL:I12.ee:end shows the many-body
spectrum $E\in\\{-2,-1,0,1,2\\}$ as discussed above. Each row represents a
particular multiplet. The combined multiplet dimension over all symmetries,
i.e., the eigenenergy’s degeneracy is specified in the second column. In the
absence of non-abelian symmetries, the second column is omitted since in that
case it would contain only trivial 1’s. In the example here, the second column
adds up to 16. This confirms the full Hilbert space dimension.

####  5.1.1 Changing local state space or symmetries

In the above discussion, the overall symmetry was specified exactly once,
namely in Eq. (LABEL:eq:gls:NC1) which defines the type of ‘site’ to use.
Since symmetries are global, all symmetries to be included need to be defined
there once and for all. If one were to switch to a different symmetry
combination, line (LABEL:eq:gls:NC1) is the only one that needs to be changed
for the entire subsequent discussion. For example, one may use the same
fermionic model, yet with NC=1 $\to$ 3 fermionic channels (flavors),

⬇

[F,Z,S,IS]=getLocalSpace(’FermionS’,’Acharge,SU2spin,SUNchannel’,’NC’,3);
?\mlbl{idx:gls:NC3}?

[F,Z,S,IS]=getLocalSpace(’FermionS’,’Acharge,SU2spin,SU3channel’);

The local state space of a site now consists of $d_{\rm loc}=4^{3}=64$ states
which can be reduced to an effective $d_{\rm loc}^{\ast}=10$ multiplets, e.g.,
as seen by inspecting the identity operator [cf. (LABEL:eq:IS:struct)],

⬇

>> IS.E ?\mlbl{IS.E:NC3}?

\parQ: 2x {10 x 4 } having ’A,SU2,SU3’ { +- } ?\mlbl{IS.E:hdr1}?

data: 2-D double (1120) 10 x 10 => 64 x 64 @ norm = 8 ?\mlbl{IS.E:hdr2}?

\par1. 1x1 | 1x1 1x1 [ -3 0 00 ; -3 0 00 ] 1\.

2. 1x1 | 2x2 3x3 [ -2 1 10 ; -2 1 10 ] 1. {?$\surd$?6}

3. 1x1 | 1x1 6x6 [ -1 0 20 ; -1 0 20 ] 1. {?$\surd$?6}

4. 1x1 | 3x3 3x3 [ -1 2 01 ; -1 2 01 ] 1. {?$\surd$?9}

5. 1x1 | 2x2 8x8 [ 0 1 11 ; 0 1 11 ] 1. {?$\surd$?16} ?\mlbl{IS.E:rec5}?

6. 1x1 | 1x1 6x6 [ 1 0 02 ; 1 0 02 ] 1. {?$\surd$?6}

7. 1x1 | 4x4 1x1 [ 0 3 00 ; 0 3 00 ] 1. {?$\surd$?4} ?\mlbl{IS.E:rec7}?

8. 1x1 | 3x3 3x3 [ 1 2 10 ; 1 2 10 ] 1. {?$\surd$?9}

9. 1x1 | 2x2 3x3 [ 2 1 01 ; 2 1 01 ] 1. {?$\surd$?6}

10. 1x1 | 1x1 1x1 [ 3 0 00 ; 3 0 00 ] 1. ?\mlbl{IS.E:NC3:end}?

\par

Here every symmetry sector on each of the 2 legs in the total of 10 records
carries 4 symmetry labels $q\equiv(q_{1},q_{2},q_{3},q_{4})\equiv(q_{1}\
q_{2}\ q_{3}q_{4})$, as also indicated with the dimensions 2x {10 x 4} in line
LABEL:IS.E:hdr1. Based on the order of symmetries requested with
getLocalSpace, the first symmetry label describes the filling relative to
half-filling. Given NC=3 flavors, this has the range
$q_{1}\in\\{-3,-2,\ldots,3\\}$ from completely empty to completely filled,
respectively. The second symmetry label describes SU(2) spin, $q_{1}=2S$, with
the largest spin multiplet given by $q_{2}=3$ (line LABEL:IS.E:rec7), i.e.,
spin $S=3/2$. It derives from the spin-half for each half-filled level, and
hence has combined symmetry labels $[0\ 3\ 00]$. The last two symmetry labels
$(q_{3}q_{4})$ describe the SU(3) channel or flavor symmetry. The largest
SU(3) multiplet above (11) $\equiv{\bf 8}$ in line LABEL:IS.E:rec5 contains 8
states (hence referred to as octet; for SU(3) this is also the adjoint
representation). Having combined symmetry labels $[0\ 1\ 11]$, again this
describes states at half-filling. With NC=3 being odd, half-filling
necessarily has a non-zero spin, here $2S=1$. Hence the combined multiplet
dimension in $[0\ 1\ 11]$ is $2\times 8=16$ states (see CGT dimensions in line
LABEL:IS.E:rec5). The QSpace displays as the ones above are just for
information. Their understanding requires at least a rudimentary understanding
of the symmetries used, including their labeling structure. The display is
generated automatically by QSpace. Looking more closely at the individual
entries permits simple consistency checks like the ones above. This ensures
that one has a good understanding of the tensors under consideration. Once a
single site is defined together with the symmetries employed, one can proceed
identically to the earlier case of a single flavor to build a many-body state
space. For example, two sites can be combined identically as in (LABEL:eq:A2)
⬇ >> A=getIdentity(Z,Z,[1 3 2]);  >> A=setitags(A,{’K01’,’K02’,’s02’})  x
In the present case, however, this results in

⬇

A = ?\mlbl{A3:disp}?

\parQ: 3x {258 x 4 } having ’A,SU2,SU3’ A-matrix, { K01, K02*, s02 }
?\mlbl{A3:hdr1}?

data: 3-D double (38.08k) 10 x 260 x 10 => 64 x 4,096 x 64 @ norm = 64
?\mlbl{A3:hdr2}?

\par1. 1x1 | 1x1x1 1x1x1 [ -3 0 00 ; -6 0 00 ; -3 0 00 ] 1\.

2. 1x2 | 1x2x2 1x3x3 [ -3 0 00 ; -5 1 10 ; -2 1 10 ] 16 b {?$\surd$?6}

3. 1x2 | 2x2x1 3x3x1 [ -2 1 10 ; -5 1 10 ; -3 0 00 ] 16 b {?$\surd$?6}

4. 1x1 | 2x1x2 3x3x3 [ -2 1 10 ; -4 0 01 ; -2 1 10 ] 1. {?$\surd$?3}

5. 1x3 | 1x1x1 1x6x6 [ -3 0 00 ; -4 0 20 ; -1 0 20 ] 24 b {?$\surd$?6}

: … ?\mlbl{A3:rec:i}?

117. 1x8x1 @2 | 2x1x2x1 8x8x8x2 [ 0 1 11 ; 0 0 11 ; 0 1 11 ] 128 b {?$\surd$?8}?\mlbl{A3:rec117}?

138. 1x12x1 @2 | 2x3x2x1 8x8x8x2 [ 0 1 11 ; 0 2 11 ; 0 1 11 ] 192 b {?$\surd$?24}?\mlbl{A3:rec138}?

: … ?\mlbl{A3:rec:j}?

257. 1x2 | 1x2x2 1x3x3 [ 3 0 00 ; 5 1 01 ; 2 1 01 ] 16 b {?$\surd$?6}

258. 1x1 | 1x1x1 1x1x1 [ 3 0 00 ; 6 0 00 ; 3 0 00 ] 1\.

\par

As compared to the $A$-tensor for NC=1 in (LABEL:eq:A2) and its output
following line LABEL:A2:disp, the fusion of two sites here leads to
considerably more entries (258 vs. 10 earlier). By default, the listing is
truncated, strongly so in the present case, skipping records 6–256 in lines
LABEL:A3:rec:i–LABEL:A3:rec:j, except for two largest entries, either in size
or outer multiplicity. Lines LABEL:A3:rec117–LABEL:A3:rec138 show records 117
and 138, both of which are entries with outer multiplicity $M=2$, denoted via
“@2” with the dimensions of the RMTs. This outer multiplicity is also
reflected in the trailing dimensions of the respective CGTs [e.g., having a
trailing “x2” for the CGTs for SU(3), as in 8x8x8x2, while having no OM with
the CGTs for SU(2), hence showing a trailing “x1” with their dimensions, as in
2x1x2x1]. To display all records in A, one can type display(A,'-f'). The
overall tensor dimensions in the header (line LABEL:A3:hdr2) are as expected,
fusing $64\times 64\to 4,096$ states. In terms of multiplets, this fuses
$10\times 10\to 260$ multiplets, total which is about an order of magnitude
smaller. From the header in line LABEL:A3:hdr2 still, the norm is given by
$\|A\|=\sqrt{4,096}=64$ also reflecting the fused state space dimension. While
from the display of the identity for a single site in lines
LABEL:IS.E:NC3–LABEL:IS.E:NC3:end all input multiplets have unique symmetry
labels showing an RMT size of 1x ..., record 138 (line LABEL:A3:rec138) here
reports 12 multiplets in the fused sector $[0\ 2\ 11]$. This combines a spin
$S=1$ with an SU(3) multiplet (11), such that these multiplets contain
$3\times 8=24$ states each (cf. also $w$-factor to the very right). The other
commands earlier with NC=1, such as the spin interaction or fermionic hopping,
can be used here also for NC=3 identically without any change. Hence with the
symmetries only specified at the very beginning of a setup with getLocalSpace,
the actual tensor network commands can be implemented the same way
irrespective of the symmetry setting, whether abelian or non-abelian
symmetries were used, whether multiple symmetries are used in parallel, or
just individual ones, or no symmetry at all. Hence once the local state space
with the desired symmetries is defined, the subsequent code can proceed nearly
as if there had been no symmetries at all. This allows one to focus on tensor
network algorithms without having to worry much about symmetries, while fully
exploiting complex symmetry settings, nevertheless. Further increasing the
symmetry for NC=3 in line LABEL:idx:gls:NC3 above to symplectic (e.g., by
having particle/hole symmetry on a bipartite lattice as well as SU(NC) channel
symmetry [8]),
⬇ [F,Z,S,IS]=getLocalSpace(’FermionS’,’SU2spin,SpNchannel’,’NC’,3);
?\mlbl{idx:gls:Sp6}?
[F,Z,S,IS]=getLocalSpace(’FermionS’,’SU2spin,Sp6channel’);  x
As seen from the fermionic parity (or equally from the identity in IS.E),

⬇

Z =

Q: 2x { 4 x 4 } having ’SU2,Sp6’ { +- }

data: 2-D double (448) 4 x 4 => 64 x 64 @ norm = 8

\par1. 1x1 | 4x4 1x1 [ 3 000 ; 3 000 ] -1. {?$\surd$?4} ?\mlbl{Sp6:Z:1}?

2. 1x1 | 3x3 6x6 [ 2 100 ; 2 100 ] 1. {?$\surd$?18} ?\mlbl{Sp6:Z:2}?

3. 1x1 | 2x2 14x14 [ 1 010 ; 1 010 ] -1. {?$\surd$?28} ?\mlbl{Sp6:Z:3}?

4. 1x1 | 1x1 14x14 [ 0 001 ; 0 001 ] 1. {?$\surd$?14} ?\mlbl{Sp6:Z:4}?

in the symplectic case the local state space of $d=64$ states has been reduced
to an effective $d^{\ast}=4$ multiplets. There are still a total of four
symmetry labels here to identify a multiplet, $q=(q_{1}\ q_{2}q_{3}q_{4})$,
but the split up changed. The first label $q_{1}=2S$ now specifies the SU(2)
spin multiplet. Since Sp(6) is a rank ${\mathfrak{r}}=3$ symmetry, it carries
the remaining three labels $(q_{2}q_{3}q_{4})$. The scalar 000 $\equiv{\bf 1}$
is shown with line LABEL:Sp6:Z:1, the defining 100 $\equiv{\bf 6}$ in line
LABEL:Sp6:Z:2. The largest Sp(6) muliplet occurs with line LABEL:Sp6:Z:4,
namely 001 $\equiv{\bf 14}$. Note that Sp(6) like the leading spin SU(2) is
self-dual. The annihilation operator

⬇

F = ?\mlbl{NC3:F}?

Q: 3x { 6 x 4 } having ’SU2,Sp6’ operator, { +– }

data: 3-D double (672) 4 x 4 x 1 => 64 x 64 x 12 @ norm = ?$\surd$?384
?\mlbl{NC3:F:hdr2}?

\par1. 1x1 | 3x4x2 6x1x6 [ 2 100 ; 3 000 ; 1 100 ] 1.15470 {?$\surd$?18}

2. 1x1 | 4x3x2 1x6x6 [ 3 000 ; 2 100 ; 1 100 ] 2.44949 {?$\surd$?4}

3. 1x1 | 2x3x2 14x6x6 [ 1 010 ; 2 100 ; 1 100 ] -1.73205 {?$\surd$?28}

4. 1x1 | 3x2x2 6x14x6 [ 2 100 ; 1 010 ; 1 100 ] -2.16025 {?$\surd$?18}

5. 1x1 | 1x2x2 14x14x6 [ 0 001 ; 1 010 ; 1 100 ] -2.44949 {?$\surd$?14}

6. 1x1 | 2x1x2 14x14x6 [ 1 010 ; 0 001 ; 1 100 ] 1.73205 {?$\surd$?28}

now has become an irop with
$({\color[rgb]{0.45,0.45,0.45}\texttt{NC=3}})\times 4=12$ components (line
LABEL:NC3:F:hdr2) that transforms like a single multiplet in the symmetry
sector [1 100], consistently also of dimension $2\times 6=12$. Each of these
components contributes equally to the norm $\|F\|^{2}={\rm
tr}(F^{\dagger}F)=12*2^{6-1}=384$ (line LABEL:NC3:F:hdr2). The SU(2) spin
operator acts trivially in the symplectic sector having 000,

⬇

S = ?\mlbl{Sp6:S}?

Q: 3x { 3 x 4 } having ’SU2,Sp6’ operator, { +– }

data: 3-D double (336) 3 x 3 x 1 => 50 x 50 x 3 @ norm = ?$\surd$?72

\par1. 1x1 | 4x4x3 1x1x1 [ 3 000 ; 3 000 ; 2 000 ] 1.936491 {?$\surd$?4} ?\mlbl{Sp6:S:rec1}?

2. 1x1 | 3x3x3 6x6x1 [ 2 100 ; 2 100 ; 2 000 ] -1.414214 {?$\surd$?18} ?\mlbl{Sp6:S:rec2}?

3. 1x1 | 2x2x3 14x14x1 [ 1 010 ; 1 010 ; 2 000 ] -0.866025 {?$\surd$?28} ?\mlbl{Sp6:S:rec3}?

The spin sectors $S=0$ are absent since these correspond to zero blocks. The
norm of the individual RMTs in lines LABEL:Sp6:S:rec1–LABEL:Sp6:S:rec3 already
reflect the Casimir for the respective local spin, having values
$(\sqrt{\frac{15}{4}},-\sqrt{2},\sqrt{\frac{3}{4}}$), and thus
$\pm\sqrt{S(S+1)}$ for $q_{1}=2S\in\\{3,2,1\\}$. This emphasizes the
importance of the relative weight across the RMTs. The relative signs are also
equally important. The tensor product space of two sites can proceed the same
way as in the initial setting in (LABEL:eq:A2). The $A$-tensor now becomes

⬇

A = ?\mlbl{A:Sp6}?

Q: 3x {61 x 4 } having ’SU2,Sp6’ A-matrix, { K01, K02*, s02 }

data: 3-D double (7.828k) 4 x 61 x 4 => 64 x 4,096 x 64 @ norm = 64

\par1. 1x4 | 1x1x1 14x1x14 [ 0 001 ; 0 000 ; 0 001 ] 32 b

2. 1x4 | 2x1x2 14x1x14 [ 1 010 ; 0 000 ; 1 010 ] 32 b

3. 1x4 | 3x1x3 6x1x6 [ 2 100 ; 0 000 ; 2 100 ] 32 b

4. 1x4 | 4x1x4 1x1x1 [ 3 000 ; 0 000 ; 3 000 ] 32 b

5. 1x1 | 1x1x1 14x84x14 [ 0 001 ; 0 002 ; 0 001 ] 1. {?$\surd$?84}

6. 1x2 | 2x1x2 14x14x14 [ 1 010 ; 0 010 ; 1 010 ] 16 b {?$\surd$?14}

: …

35. 1x6 | 3x3x3 6x14x6 [ 2 100 ; 2 010 ; 2 100 ] 48 b {?$\surd$?42}?\mlbl{A:Sp6:intermediate:recs}?

36. 1x6 | 4x3x2 1x14x14 [ 3 000 ; 2 010 ; 1 010 ] 48 b {?$\surd$?42}

: …

60. 1x2 | 4x6x3 1x6x6 [ 3 000 ; 5 100 ; 2 100 ] 16 b {?$\surd$?36}

61. 1x1 | 4x7x4 1x1x1 [ 3 000 ; 6 000 ; 3 000 ] 1. {?$\surd$?7}

This reduces the full Hilbert space for the two sites from $D=64^{2}=4096$
states to an effective $D^{\ast}=61$ multiplets across 23 symmetry sectors
[the latter can be seen from getIdentity(A,2)]. The matrix elements in
(LABEL:eq:SdotS) or (LABEL:eq:FdotF:T12) can be computed as previously without
any change.

Figure 10:  (a) The swap operator $X_{12}$ acting on two sites $i\in\\{1,2\\}$
with respective local state spaces $\sigma_{i}$. (b) The swap operator is
invariant under arbitrary unitary transformations on the local state spaces
since $UU^{\dagger}=1$ annihilate each other along the diagonal lines. (c) The
matrix elements of the swap operator expressed in the combined state space as
described by an $A$-tensor. (d-e) The matrix elements of
$(S_{1}^{\dagger}\cdot S_{2}{\,})^{n}$ for $n=1$ and $2$, respectively. The
contraction in (d) is the same as in Fig. 9, except that spin operators are
applied here. Similarly, also the conjugate operator for the first site is
drawn as the vertical mirror image as indicated by the gray circular arrow and
line. The intermediate object in the contraction indicated by the blue dashed
box in (d), referred to as SAS in the text (line LABEL:SAS), can be reused in
(e).

###  5.2 Swap operator for two spins $S=1$

The following proceeds similar to the previous two-site examples [e.g., see
(LABEL:eq:A2)], yet switches from a fermionic state space for a site to a spin
$S=1$ in SU(2),

⬇

[S,IS]=getLocalSpace(’Spin’,’SU2’,2); ?\mlbl{S1:op}? A=getIdentity(S,S,[1 3
2]); ?\mlbl{A2:S}? A=setitags(A,{’K01’,’K02’,’s02’});

⬇

>> A

\parQ: 3x { 3 x 1 } having ’SU2’ A-matrix, { K01, K02*, s02 }

data: 3-D double (336) 1 x 3 x 1 => 3 x 9 x 3 @ norm = 3

\par1. 1x1 | 3x1x3 [ 2 ; 0 ; 2 ] 1\.

2. 1x1 | 3x3x3 [ 2 ; 2 ; 2 ] 1. {?$\surd$?3}

3. 1x1 | 3x5x3 [ 2 ; 4 ; 2 ] 1. {?$\surd$?5}

The state space of two sites as described by A results in a total spin
$S\in\\{0,2,4\\}/2=\\{0,1,2\\}$, as expected. The following then computes the
‘interaction’ $X$ that describes a plain swap of the state space of two sites,
also known as the permutation operator. This operator $X$ is Hermitian with
eigenvalues $\pm 1$ since $X^{2}=1$. Its matrix elements can be obtained as
follows,

⬇

A_= untag(A); ?\mlbl{Xop:A_}? X = contract(A_,’13*’,A_,’31’); ?\mlbl{Xop:ctr}?

⬇

>> X ?\mlbl{X:swp}?

Q: 2x { 3 x 1 } having ’SU2’ { +- }

data: 2-D double (336) 3 x 3 => 9 x 9 @ norm = 3

\par1. 1x1 | 1x1 [ 0 ; 0 ] 1. ?\mlbl{X:swp:rec1}?

2. 1x1 | 3x3 [ 2 ; 2 ] -1. {?$\surd$?3} ?\mlbl{X:swp:rec2}?

3. 1x1 | 5x5 [ 4 ; 4 ] 1. {?$\surd$?5} ?\mlbl{X:swp:recl}?

This swap operator $X$ is already diagonal with eigenvalues $\pm 1$ as shown
with the RMTs. The itags in the contracted local indices were temporarily
removed (untag in line LABEL:Xop:A_), since the swap operator contracts legs
with otherwise different itags [see Fig. 10(a)]. Leaving the itags in place
would result in an error due to itag mismatch [cf. Fig. 10(c)]. The
contraction in line LABEL:Xop:ctr uses compact index notation. It contracts
legs $1$ and $3$, i.e., $K_{1}\equiv\sigma_{1}$ and $\sigma_{2}$ of the
conjugate of the first tensor with legs $3$ and $1$, i.e., $\sigma_{2}$ and
$\sigma_{1}$ and thus reverse order, on the second tensor (tensor conjugation
leaves the order of legs intact). Now an interaction of two spin-$S$ sites can
be decomposed into powers $(S_{1}^{\dagger}\cdot S_{2}^{\,})^{n}$ of the spin-
spin interaction which forms a (non-orthonormalized) operator basis. Clearly,
by only having $2S+1$ RMTs, each of which is given by a single number, one
reaches a complete basis for $n\leq 2S$. Hence powers $n\geq 2S+1$ can be
expressed in terms of smaller powers. For example for spin-half sites the
identity and the spin operator already exhaust the complete set of local
operators, since $2^{2}=1+3$. For the swap operator here for $S=1$ sites, this
implies

$\displaystyle X=\sum_{n=0}^{2}a_{n}\ (S_{1}^{\dagger}\cdot
S_{2}^{\,})^{n}\text{ .}$ (66)

with coefficients $a_{n}$ to be determined. The operator for $n=0$ refers to
the identity which has the same QSpace structure as displayed with X in line
LABEL:X:swp, except that it has all +1’s in the right column for the RMTs. The
case $n=1$ gives the Heisenberg interaction SS1 [cf. Fig. 10(d)],

⬇

>> SAS = contract(S,’-op:K01’,’*’,{A,S,’-op:^s’}); ?\mlbl{SAS}?

>> SS1 = contract(A,’!2*’,SAS) ?\mlbl{SS1}?

\parQ: 2x { 3 x 1 } having ’SU2’ { K02, K02* }

data: 2-D double (336) 3 x 3 => 9 x 9 @ norm = ?$\surd$?12

\par1. 1x1 | 1x1 [ 0 ; 0 ] -2.

2. 1x1 | 3x3 [ 2 ; 2 ] -1. {?$\surd$?3}

3. 1x1 | 5x5 [ 4 ; 4 ] 1. {?$\surd$?5}

and $n=2$ the quadrupolar term SS2 [cf. Fig. 10(e)],

⬇

>> SS2 = contract(A,’!2*’,{S,’-op:K01’,’*’,{SAS,S,’-op:^s’}})

\parQ: 2x { 3 x 1 } having ’SU2’ { K02, K02* }

data: 2-D double (336) 3 x 3 => 9 x 9 @ norm = ?$\surd$?24

\par1. 1x1 | 1x1 [ 0 ; 0 ] 4\.

2. 1x1 | 3x3 [ 2 ; 2 ] 1. {?$\surd$?3}

3. 1x1 | 5x5 [ 4 ; 4 ] 1. {?$\surd$?5}

Now by collecting the RMTs for $n=0,1,2$ , as well as for X in lines
LABEL:X:swp:rec1–LABEL:X:swp:recl into columns, one can determine the
coefficients $a_{n}$ in Eq. (66),

$\displaystyle\begin{pmatrix}1\\\ -1\\\
1\end{pmatrix}\overset{!}{=}\begin{pmatrix}1&-2&4\\\ 1&-1&1\\\
1&1&1\end{pmatrix}\begin{pmatrix}a_{0}\\\ a_{1}\\\
a_{2}\end{pmatrix}\quad\Rightarrow\quad\begin{pmatrix}a_{0}\\\ a_{1}\\\
a_{2}\end{pmatrix}=\begin{pmatrix}-1\\\ 1\\\ 1\end{pmatrix}$ (67)

and therefore

$\displaystyle X=-1+S_{1}^{\dagger}\cdot S_{2}^{\,}+(S_{1}^{\dagger}\cdot
S_{2}^{\,})^{2}\text{ .}$ (68)

A simple countercheck confirms the coefficients in Eq. (68), indeed,

⬇

>> E2 = getIdentity(A,2);

>> X_ = -E2 + SS1 + SS2

\parQ: 2x { 3 x 1 } having ’SU2’ { +- }

data: 2-D double (336) 3 x 3 => 9 x 9 @ norm = 3

\par1. 1x1 | 1x1 [ 0 ; 0 ] 1\.

2. 1x1 | 3x3 [ 2 ; 2 ] -1. {?$\surd$?3}

3. 1x1 | 5x5 [ 4 ; 4 ] 1. {?$\surd$?5}

with norm(X_-X) returning numerical double precision noise $<10^{-14}$. This
not being exactly zero is due to the square root factors shown to the right.
For this reason, X==X_ returns false. The QSpace routine sameas(X,X_) based on
norm(X_-X) returns true.

##### Physical interpretation

By using a simple swap operator on local state spaces of dimension $d_{\rm
loc}=3$, an arbitrary 3-dimensional unitary $U$ can be inserted as
$UU^{\dagger}$ into any of the two lines in Fig. 10(a). This takes the
interpretation of transforming the local state space for each site. However,
with $X_{12}\to(U\otimes U)\,X_{12}(U^{\dagger}\otimes
U^{\dagger})=(UU^{\dagger}\otimes 1)(1\otimes UU^{\dagger})X_{12}=X_{12}$
being identical to $X$ [most easily seen by considering the pictorial
representation in Fig. 10(a)], this implies that the swap operation in the
above Hamiltonian actually has an enlarged symmetry, namely SU($d_{\rm loc}$).
For the case of the spin $S=1$ here, this is the well-known SU(3) symmetric
point of the bilinear-biquadratic Heisenberg model for $S=1$ at the phase
boundary of the Haldane phase [51, 52].

###  5.3 Tweaking local state space

The local state space of a site is set up via a call to getLocalSpace which
has many standard symmetry setups already implemented. Nevertheless,
situations can arise that may go beyond to what is already provided by
getLocalSpace. The typical way to proceed then is to tweak its output starting
from a symmetry setup that may be considered closest to target setup. To be
specific, consider for example SU(2) spin-halfs with an additional two
symmetric orbital flavors, i.e., ${\rm SU}(2)_{\rm spin}\otimes{\rm
SU}(2)_{\rm orbital}$. This does not exist ready-made in getLocalSpace. Two
ways to deal with this situation are discussed in what follows.

####  5.3.1 Project fermionic setup to spin model and complete operator basis

Starting from a spinfull fermionic model with NC=2 orbital flavors
⬇ [F,Z,S,IS]=getLocalSpace(’FermionS’,’Acharge,SU2spin,SU2channel’);  x
the state space can be inspected by Z or IS.E,

⬇

>> E=IS.E

\parQ: 2x { 6 x 3 } having ’A,SU2,SU2’ { +- }

data: 2-D double (672) 6 x 6 => 16 x 16 @ norm = 4

\par1. 1x1 | 1x1 1x1 [ -2 0 0 ; -2 0 0 ] 1\.

2. 1x1 | 2x2 2x2 [ -1 1 1 ; -1 1 1 ] 1. {?$\surd$?4} ?\mlbl{EP:rec2}?

3. 1x1 | 1x1 3x3 [ 0 0 2 ; 0 0 2 ] 1. {?$\surd$?3}

4. 1x1 | 3x3 1x1 [ 0 2 0 ; 0 2 0 ] 1. {?$\surd$?3}

5. 1x1 | 2x2 2x2 [ 1 1 1 ; 1 1 1 ] 1. {?$\surd$?4}

6. 1x1 | 1x1 1x1 [ 2 0 0 ; 2 0 0 ] 1\.

Looking for the spin $S=1/2$ sector at a filling of $n=1$ particle on the
site, this corresponds to $q_{1}=-1$ relative to half-filling. This contains a
single multiplet in line LABEL:EP:rec2 which is thus the one of interest. It
can be projected by just picking this record from the identity using getsub,

⬇

>> q0=-1; i=find(Z.Q{1}(:,1)==q0); ?\mlbl{gls:tweak:q0}? >> E=getsub(E,i)
\parQ: 2x { 1 x 3 } having ’A,SU2,SU2’ { +- }

data: 2-D double (112) 1 x 1 => 4 x 4 @ norm = 2

\par1. 1x1 | 2x2 2x2 [ -1 1 1 ; -1 1 1 ] 1. {?$\surd$?4}

\par

Since now only a trivial single abelian symmetry sector is present in terms of
charge, it can be skipped (removed) as follows,

⬇

>> E=rmAbelian(E,1) ?\mlbl{E:rmAbelian}?

\parE =

Q: 2x { 1 x 2 } having ’SU2,SU2’ { +- }

data: 2-D double (112) 1 x 1 => 4 x 4 @ norm = 2

\par1. 1x1 | 2x2 2x2 [ 1 1 ; 1 1 ] 1. {?$\surd$?4}

\par>> clear F Z ?\mlbl{clearFZ}?

The syntax in line LABEL:E:rmAbelian removes the symmetry at position 1 from
E, i.e., the abelian U(1) charge. Projecting to a single charge sector renders
the fermionic operators irrelevant (line LABEL:clearFZ). The spin operator,
finally, can be projected the same way as the identity above,

⬇

>> S=getsub(find(S.Q{1}(:,1)==q0)); >> S=rmAbelian(S,1) ?\mlbl{Sop:tweak}?

\parS =

Q: 3x { 1 x 2 } having ’SU2,SU2’ operator, { +– }

data: 3-D double (112) 1 x 1 x 1 => 4 x 4 x 3 @ norm = ?$\surd$?3

\par1. 1x1 | 2x2x3 2x2x1 [ 1 1 ; 1 1 ; 2 0 ] -0.866025 {?$\surd$?4} ?\mlbl{Sop:tweak:rec1}?

which transforms like a $S=q/2$ $=1$ irop. This is the same spin operator as
for a single fermionic flavor [cf. (LABEL:eq:gls:NC1), and subsequent display
in lines LABEL:idx:Sop:SU2ff] that acts trivially, i.e., like an identity in
the orbital space. The pseudo-spin operator in the orbital space can be
obtained from a complete operator basis for the local state space. For its
construction, one needs to obtain the tensor product space in E with the dual
of its state space as obtained by the $1j$-symbol in line LABEL:Xb:U,

⬇

>> U=getIdentity(E,’-0’); ?\mlbl{Xb:U}?

>> X=getIdentity(U,1,U,2); ?\mlbl{Xb:1}?

>> X=contract(U,’!1*’,X,[2 1 3]) ?\mlbl{Xb:2}?

\parQ: 3x { 4 x 2 } having ’SU2,SU2’ { +– } ?\mlbl{Xb:hdr1}?

data: 3-D double (448) 1 x 1 x 4 => 4 x 4 x 16 @ norm = 4 ?\mlbl{Xb:hdr2}?

\par1. 1x1 | 2x2x1 2x2x1 [ 1 1 ; 1 1 ; 0 0 ] 0.5 {?$\surd$?4} ?\mlbl{Xb:rec1}?

2. 1x1 | 2x2x1 2x2x3 [ 1 1 ; 1 1 ; 0 2 ] 0.866025 {?$\surd$?4} ?\mlbl{Xb:rec2}?

3. 1x1 | 2x2x3 2x2x1 [ 1 1 ; 1 1 ; 2 0 ] 0.866025 {?$\surd$?4} ?\mlbl{Xb:rec3}?

4. 1x1 | 2x2x3 2x2x3 [ 1 1 ; 1 1 ; 2 2 ] 1.5 {?$\surd$?4} ?\mlbl{Xb:recl}?

The tensor product in line LABEL:Xb:1 derives from the fact that the full
operator space of a state space of dimension $d$ is described by $d^{2}\equiv
d\otimes d$ operators. However, the incoming dual state space, i.e., primed
index in X still needs to be reverted by contracting the same $1j$-symbol in
line LABEL:Xb:2, so that X represents a standard operator [cf. Fig. 1]. The
trailing permutation [2 1 3] restores the operator index order +--, as seen in
line LABEL:Xb:hdr1. Having a dimension of 1x1x4 in terms of multiplets (line
LABEL:Xb:hdr2), the object X is clearly not an irop, though. Rather it
contains a (direct) sum of four irops, which spans the full operator space of
a single site. This consists of the four records in above display, each of
which represents an irop with a unique operator $q$-label:

  * •

Record 1. (line LABEL:Xb:rec1) transforms like a scalar $q_{\rm op}=[0\ 0]$.
Since this is the only operator of this type and bearing in mind that X
represents a complete operator basis, it must be the identity operator up to
the scale factor given by RMT value 0.5 shown to the right. Via the tensor-
product in getIdentity in line LABEL:Xb:1 each operator has Frobenius norm 1
which explains this factor, since $\|\frac{1}{2}11^{(4)}\|=1$.

  * •

Record 3. (line LABEL:Xb:rec3) represents an irop that transforms according to
$q_{\rm op}=[2\ 0]$. This identifies it as the spin operator like in line
LABEL:Sop:tweak:rec1 above up to a sign (note that the sign is irrelevant in
bilinear products, as long as the same operator S is used consistently
throughout). Incidentally, in the present setting the value of its RMT
0.866025$=\frac{\sqrt{3}}{2}$ is consistent already with both, the
normalization of a spin-half operator ($\hat{S}^{2}=\frac{3}{4}\,11^{(4)}$),
as well as the normalization out of the tensor-product in getIdentity in line
LABEL:Xb:1, having $\|\hat{S}\|^{2}={\rm tr}(\frac{3}{4}11^{(4)})=3$
representing three normalized operators, indeed.

  * •

Record 2. (line LABEL:Xb:rec2) represents an irop that transforms according to
$q_{\rm op}=[0\ 2]$. Therefore this entry describes the pseudo-spin operator
in orbital space.

  * •

Record 4. (line LABEL:Xb:recl), finally, represents an irop that transforms
according to $q_{\rm op}=[2\ 2]$. As such this is a spin-orbit interaction ST
up to a factor $2$ when comparing to RMT value of 1.5 to the standard SU(2)
spin normalization
$\sqrt{(\hat{S}\hat{T})^{2}}=\sqrt{\hat{S}^{2}\hat{T}^{2}}=\frac{3}{4}\,11^{(4)}$.

By construction, the above full operator basis in X permits one to extract all
(other) irreducible operators within the present symmetry setting. Here this
includes the pseudo-spin operator T in orbital space (choosing the same sign
convention as in S),

⬇

>> T=-getsub(X,3) ?\mlbl{gls:X:T}?

\parQ: 3x { 1 x 2 } having ’SU2,SU2’ { +– }

data: 3-D double (112) 1 x 1 x 1 => 4 x 4 x 3 @ norm = ?$\surd$?3

\par1. 1x1 | 2x2x3 2x2x1 [ 1 1 ; 1 1 ; 2 0 ] -0.866025 {?$\surd$?4}

\par

The remaining last operator within the present symmetry setting is the spin-
orbit term

⬇

>> ST= getsub(X,4)/2; ?\mlbl{gls:X:ST}?

\parQ: 3x { 1 x 2 } having ’SU2,SU2’ { +– }

data: 3-D double (112) 1 x 1 x 1 => 4 x 4 x 9 @ norm = 1.5

\par1. 1x1 | 2x2x3 2x2x3 [ 1 1 ; 1 1 ; 2 2 ] 0.75 {?$\surd$?4}

The combination of these four irops that describe spinors of dimension
$1+3+3+3\times 3=16=4^{2}$ thus, indeed, exhausts the operator basis of given
local state space.

####  5.3.2 Appending additional symmetries

Alternative to the previous projective approach, one may tweak the output of
getLocalSpace by adding other symmetries still. Starting with a plain SU(2)
spin-half, one can manually add an additional SU(2) symmetry, here to the
identity operator returned by getLocalSpace. For the example above, this
starts here with the simpler spin model,

⬇

[S,IS]=getLocalSpace(’Spin’,’SU2’,1); E=addSymmetry(IS.E,’SU2’,’q’,1)
?\mlbl{E:addSym}?

⬇

E =

Q: 2x { 1 x 2 } having ’SU2,SU2’ { +- }

data: 2-D double (112) 1 x 1 => 4 x 4 @ norm = 2 ?\mlbl{gls:X:E:hdr2}?

\par1. 1x1 | 2x2 2x2 [ 1 1 ; 1 1 ] 1. {?$\surd$?4}

By default, addSymmetry appends the specified symmetry in the scalar irep
$q=0$. By explicitly specifying $q=1$ in line LABEL:E:addSym, this adds the
second SU(2) symmetry in the form of an $S=1/2$ multiplet. It is interpreted
here as SU(2) symmetric orbital. The local state space therefore becomes
$2\times 2=4$ dimensional as seen in line LABEL:gls:X:E:hdr2. Since
addSymmetry always applies a square-root factor of the multiplet dimension
added, the identity E already has the correct normalization (which matter of
fact, is the motivation as to why addSymmetry behaves this way). Next the spin
operators are tweaked. First, this generates the orbital spin operator T in
line LABEL:T:addSym based on S, before adapting the spin operator S itself in
line LABEL:S:addSym,

⬇

T=addSymmetry(S,’SU2’,’pos’,1,’q’,1); ?\mlbl{T:addSym}?
S=addSymmetry(S,’SU2’,’q’,1); ?\mlbl{S:addSym}?

Line LABEL:T:addSym uses the option 'pos',1 which inserts the new symmetry at
location 1, i.e., prepends the symmetry rather than appending it (see help
addSymmetry for more detailed usage information in this regard). The display
of the operators above reads

⬇

\parT =

Q: 3x { 1 x 2 } having ’SU2,SU2’ operator, { +– }

data: 3-D double (112) 1 x 1 x 1 => 4 x 4 x 3 @ norm = ?$\surd$?3
?\mlbl{T:hdr2}?

\par1. 1x1 | 2x2x1 2x2x3 [ 1 1 ; 1 1 ; 0 2 ] -0.866025 {?$\surd$?4}

\parS =

Q: 3x { 1 x 2 } having ’SU2,SU2’ operator, { +– }

data: 3-D double (112) 1 x 1 x 1 => 4 x 4 x 3 @ norm = ?$\surd$?3
?\mlbl{S:hdr2}?

\par1. 1x1 | 2x2x3 2x2x1 [ 1 1 ; 1 1 ; 2 0 ] -0.866025 {?$\surd$?4}

which is consistent with the alternative approach earlier. With respect to the
operator index (index 3) the above shows that $S$ ($T$) acts like a spin
operator in the first (second) symmetry, respectively, thus confirming ${\rm
SU}(2)_{\rm spin}\otimes{\rm SU}(2)_{\rm orbital}$ in this order. Similarly to
E above, the operators S and T are again also already correctly normalized
($\surd 3$ in header lines LABEL:T:hdr2 and LABEL:S:hdr2), since, e.g.,
$\|S\|^{2}={\rm tr}(\ 1^{(2)}\otimes\frac{3}{4}1^{(2)}\,)=3$.

###  5.4 Local operators beyond getLocalSpace

The routine getLocalSpace returns a set of elementary operators that are
representative for the chosen site. Many operators can be derived in simple
one-liners as already shown earlier. However, there are exceptions as already
partly discussed in the previous section, e.g., when generating the orbital
pseudy spin operator T in line LABEL:gls:X:T or the spin-orbit operator TS in
line LABEL:gls:X:ST from a complete operator basis. This approach is reviewed
once more in a more elaborate example for orbital spin operators below,
together with an alternative strategy for obtaining local spin operators based
on a Schrieffer-Wolff-like construction. Both represent instructive approaches
for generating additional local operators. The following discussion is based
on the fermionic state space for NC=2 orbitals in (LABEL:eq:gls:SU2x2) with a
local Hilbert space of $d_{\rm loc}=4^{2}=16$ states. The spin operator S for
the entire fermionic state space as returned by (LABEL:eq:gls:SU2x2) reads

⬇

>> S = ?\mlbl{gls:NC2:S}?

Q: 3x { 3 x 3 } having ’A,SU2,SU2’ operator, { +– }

data: 3-D double (336) 3 x 3 x 1 => 11 x 11 x 3 @ norm = ?$\surd$?12

\par1. 1x1 | 2x2x3 2x2x1 [ -1 1 1 ; -1 1 1 ; 0 2 0 ] -0.866025 {?$\surd$?4}?\mlbl{gls:NC2:S:rec1}?

2. 1x1 | 3x3x3 1x1x1 [ 0 2 0 ; 0 2 0 ; 0 2 0 ] -1.41421 {?$\surd$?3}?\mlbl{gls:NC2:S:rec2}?

3. 1x1 | 2x2x3 2x2x1 [ 1 1 1 ; 1 1 1 ; 0 2 0 ] -0.866025 {?$\surd$?4}?\mlbl{gls:NC2:S:rec3}?

\par

As with the earlier symplectic example in the symplectic case, the norm of the
individual RMTs already reflect the Casimir for the respective local spin,
having values $(\sqrt{\frac{3}{4}},\sqrt{2},\sqrt{\frac{3}{4}}$), and thus
$\sqrt{S(S+1)}$ for $q_{1}=2S\in\\{1,2,1\\}$. The normalization as well as the
(relative) signs of the records in the irop S are important. While the signs
here all the same, this is not always the case, e.g., as seen in the earlier
example. The global sign of any spin operator is fixed after defining the sign
of the spin operator in the defining representation. Furthermore, since the
spin operator is well-defined in the defining representation, so it must be in
any other many-body state space derived from it. Hence the spin operator,
representing ‘total spin’, is well-defined and thus unique in any state space.
The spin operator above has irop labels $q_{S}=$[0 2 0]. Based on
(LABEL:eq:gls:SU2x2) with $q=(q_{1},q_{2},q_{3})$, these symmetry labels are
$q_{1}$ for abelian U(1) charge, q2=2S for the spin SU(2), and q3=2S for the
orbital SU(2). As seen from the irop labels $q_{S}$, the spin operator acts
non-trivially only in the SU(2) spin sector. There it reflects the adjoint
representation ($S=1)$ that incorporates the full set of generators for that
symmetry.

##### Complete operator basis

Any local operator lives within the complete set of local operators. While
constructing the latter may appear overkill for specific local operators, this
approach is instructive, nevertheless. The complete operator basis Eop is
built starting from the tensor product of the local state space with its dual
via a $1j$-tensor (line LABEL:SU2x2:U), followed by a contraction of the same
$1j$-tensor in the line below [using the local operators out of getLocalSpace
as in (LABEL:eq:gls:SU2x2); same construction as in lines
LABEL:Xb:U–LABEL:Xb:2]

⬇

>> U = getIdentity(IS.E,’-0’); ?\mlbl{SU2x2:U}?

>> Eop = contract(getIdentity(U,1,U,2),2,U,’2*’,[1 3 2]) ?\mlbl{SU2x2:Eop}?

⬇

Q: 3x {60 x 3 } having ’A,SU2,SU2’ { +– }

data: 3-D double (7.734k) 6 x 6 x 60 => 16 x 16 x 256 @ norm = 16

\par1. 1x1x6 | 1x1x1 1x1x1 [ -2 0 0 ; -2 0 0 ; 0 0 0 ] 48 b

2. 1x1x6 | 2x1x2 2x1x2 [ -1 1 1 ; -2 0 0 ; 1 1 1 ] 48 b {?$\surd$?4}

3. 1x1x3 | 1x1x1 3x1x3 [ 0 0 2 ; -2 0 0 ; 2 0 2 ] 24 b {?$\surd$?3}

4. 1x1x3 | 3x1x3 1x1x1 [ 0 2 0 ; -2 0 0 ; 2 2 0 ] 24 b {?$\surd$?3}

5. 1x1x2 | 2x1x2 2x1x2 [ 1 1 1 ; -2 0 0 ; 3 1 1 ] 16 b {?$\surd$?4}

6. 1x1 | 1x1x1 1x1x1 [ 2 0 0 ; -2 0 0 ; 4 0 0 ] 1\.

: …

59. 1x1x6 | 2x1x2 2x1x2 [ 1 1 1 ; 2 0 0 ; -1 1 1 ] 48 b {?$\surd$?4}

60. 1x1x6 | 1x1x1 1x1x1 [ 2 0 0 ; 2 0 0 ; 0 0 0 ] 48 b

By construction, Eop encodes all $(d_{\rm loc}=16)^{2}=256$ local operators.
It represents 60 irops that transform in a total of 24 different ireps, as
seen from getIdentity(Eop,3). If one were to project the local state space to
a particular single multiplet, this would generate an effective spin model.
Only a global normalization of irops has to be fixed then. For the full local
state space, the spin operator as returned by getLocalSpace must also be
contained in Eop . A search for all records in Eop that have the symmetry
labels $q_{S}=$[0 2 0] on the operator index $3$ via getsub results in,

⬇

>> S_ = getsub(Eop,[0 2 0],3)

\parQ: 3x { 3 x 3 } having ’A,SU2,SU2’ { +– }

data: 3-D double (384) 3 x 3 x 3 => 11 x 11 x 9 @ norm = 3

\par1. 1x1x3 | 2x2x3 2x2x1 [ -1 1 1 ; -1 1 1 ; 0 2 0 ] 24 b {?$\surd$?4}

2. 1x1x3 | 3x3x3 1x1x1 [ 0 2 0 ; 0 2 0 ; 0 2 0 ] 24 b {?$\surd$?3}

3. 1x1x3 | 2x2x3 2x2x1 [ 1 1 1 ; 1 1 1 ; 0 2 0 ] 24 b {?$\surd$?4}

\par

It describes an operator subspace of three entries (having dimensions 1x1x3
for all RMTs) which is solely due to having three QSpace records here. The
spin operator can be built as some linear superposition within that subspace.
However, a priori, it is unclear how to obtain the correct coefficients for
this linear superposition as in lines
LABEL:gls:NC2:S:rec1–LABEL:gls:NC2:S:rec3 from the present setting. The spin
operator for the SU(2) spin symmetry is known, of course, because it is
explicitly returned by getLocalSpace. This has the benefit that one has a
reference to compare various approaches to. If one needs the pseudo-spin
operator in the full local state space above in the orbital space, or the
intertwined spin-orbit operator, these are no longer available from
getLocalSpace. For these, one needs to address the problem of how to make the
spin operator unique, at least up to a global scale factor. For this purpose,
one may target the spin operators more concretely as discussed next.

##### Spin operators from second order fermionic hopping

An alternative route to generalized ‘spin operators’ is physically motivated
by Schrieff-Wolff transformations, where a second-order perturbation in the
tunneling maps fermionic cotunneling to an effective ‘spin model’. There with
$\hat{f}_{i}$ the complete set of local fermionic operators with
$i=1,\ldots,N_{f}$ indexing all flavors or flavor combinations, one needs to
decompose the pool of operators
$\\{\hat{f}_{i}^{\dagger}\hat{f}_{i^{\prime}}^{\,}\\}$ into irops. Hence
rather than building the full operator space for the local state space, one
can constrain oneself to the operators obtained out of the bilinear
F${}^{\dagger}\otimes$F, here with a Kronecker rather than a dot-product on
the irop indices. Since F† as a conjugate operator has the irop index ingoing,
this procedure is equivalent to the construction of the complete operator
basis in lines LABEL:SU2x2:U–LABEL:SU2x2:Eop, except that one just takes the
irop space of F, instead. Replacing IS.E in line LABEL:SU2x2:U by F on leg 3,
one obtains,

⬇

>> U_ = getIdentity(F,3,’-0’); ?\mlbl{SU2x2:Eop:F}?

Eop_ = contract(getIdentity(U_,1,U_,2),2,U_,’2*’,[1 3 2]) ?\mlbl{SU2x2:Eop_}?

⬇

Q: 3x { 4 x 3 } having ’A,SU2,SU2’ { +– }

data: 3-D double (448) 1 x 1 x 4 => 4 x 4 x 16 @ norm = 4
?\mlbl{SU2x2:Eop_:hdr2}?

\par1. 1x1 | 2x2x1 2x2x1 [ -1 1 1 ; -1 1 1 ; 0 0 0 ] 0.5 {?$\surd$?4}?\mlbl{SU2x2:Eop_:rec1}?

2. 1x1 | 2x2x1 2x2x3 [ -1 1 1 ; -1 1 1 ; 0 0 2 ] 0.866025 {?$\surd$?4}?\mlbl{SU2x2:Eop_:rec2}?

3. 1x1 | 2x2x3 2x2x1 [ -1 1 1 ; -1 1 1 ; 0 2 0 ] 0.866025 {?$\surd$?4}?\mlbl{SU2x2:Eop_:rec3}?

4. 1x1 | 2x2x3 2x2x3 [ -1 1 1 ; -1 1 1 ; 0 2 2 ] 1.5 {?$\surd$?4}?\mlbl{SU2x2:Eop_:rec4}?

\par

The spinfull NC=2 model here based on (LABEL:eq:gls:SU2x2) has a combined
total of $N_{f}=2\times 2=4$ flavors. Hence the pool of operators
$\\{\hat{f}_{i}^{\dagger}\hat{f}_{i^{\prime}}^{\,}\\}$ has $4\times 4=16$
entries. These are cast into three irops here that span the full space of the
operator pool, as seen from the overall dimensions 4 x 4 x 16 in line
LABEL:SU2x2:Eop_:hdr2. By having RMTs of dimension 1x1 [x1] only, each record
now corresponds to an irop that is unique up to normalization. The tensor Eop_
still needs to be contracted onto the two irop indices in F†F. By choosing the
SU(2) spin irop in line LABEL:SU2x2:Eop_:rec2, one obtains the spin operator

⬇

>> x = getsub(Eop_,[0 2 0],3); Sop = contract({F,’1*’,F,1},’24’,x,’21’)

⬇

Q: 3x { 3 x 3 } having ’A,SU2,SU2’ { +– }

data: 3-D double (336) 3 x 3 x 1 => 11 x 11 x 3 @ norm = ?$\surd$?12
?\mlbl{SW:Sop}?

\par1. 1x1 | 2x2x3 2x2x1 [ -1 1 1 ; -1 1 1 ; 0 2 0 ] -0.866025 {?$\surd$?4}

2. 1x1 | 3x3x3 1x1x1 [ 0 2 0 ; 0 2 0 ; 0 2 0 ] -1.41421 {?$\surd$?3}

3. 1x1 | 2x2x3 2x2x1 [ 1 1 1 ; 1 1 1 ; 0 2 0 ] -0.866025 {?$\surd$?4}

\par

This is precisely the spin operator returned by getLocalSpace in lines
LABEL:gls:NC2:S–LABEL:gls:NC2:S:rec3, already also with the correct
normalization, as confirmed by norm(S-Sop)$\leq 10^{-14}$. The normalization
can be easily double checked, bearing in mind that the spin operator returned
by getLocalSpace describes the total spin of the local state space. For NC=2
fermionic channels with $S_{i}$ the spin operator for channel (orbital) $i$,
this leads to $\|S\|^{2}=\|S_{1}+S_{2}\|^{2}=2\,\|S_{1}\|^{2}\,{\rm
tr}(11^{(4)}_{2})=2\cdot(\frac{3}{4}2)\cdot 4=12$, in agreement with line
LABEL:SW:Sop. The present approach now also allows to extend to other
generalized spin operators. For example, one may choose the SU(2) pseudo-spin
irop in the orbital sector (line LABEL:SU2x2:Eop_:rec3),

⬇

>> x = getsub(Eop_,[0 0 2],3); Top = contract({F,’1*’,F,1},’24’,x,’21’)

⬇

Q: 3x { 3 x 3 } having ’A,SU2,SU2’ { +– }

data: 3-D double (336) 3 x 3 x 1 => 11 x 11 x 3 @ norm = ?$\surd$?12

\par1. 1x1 | 2x2x1 2x2x3 [ -1 1 1 ; -1 1 1 ; 0 0 2 ] -0.866025 {?$\surd$?4}

2. 1x1 | 1x1x1 3x3x3 [ 0 0 2 ; 0 0 2 ; 0 0 2 ] -1.41421 {?$\surd$?3}

3. 1x1 | 2x2x1 2x2x3 [ 1 1 1 ; 1 1 1 ; 0 0 2 ] -0.866025 {?$\surd$?4}

which shows the same RMT values as the previous spin operator 'Sop', or the
spin-orbital sector in line LABEL:SU2x2:Eop_:rec4,

⬇

>> x = getsub(Eop_,[0 2 2],3); STop = contract({F,’1*’,F,1},’24’,x,’21’) / 2

⬇

Q: 3x { 4 x 3 } having ’A,SU2,SU2’ { +– }

data: 3-D double (448) 4 x 4 x 1 => 14 x 14 x 9 @ norm = 3

\par1. 1x1 | 2x2x3 2x2x3 [ -1 1 1 ; -1 1 1 ; 0 2 2 ] 0.75 {?$\surd$?4}

2. 1x1 | 1x3x3 3x1x3 [ 0 0 2 ; 0 2 0 ; 0 2 2 ] 0.866025 {?$\surd$?3}

3. 1x1 | 3x1x3 1x3x3 [ 0 2 0 ; 0 0 2 ; 0 2 2 ] 0.866025 {?$\surd$?3}

4. 1x1 | 2x2x3 2x2x3 [ 1 1 1 ; 1 1 1 ; 0 2 2 ] -0.75 {?$\surd$?4}

In the last case the norm was adapted to yield the expected value
$\|{\color[rgb]{0.45,0.45,0.45}\texttt{STop}}\|^{2}=2\cdot
2\cdot(\frac{3}{4}2)^{2}=9$, having $S=S_{1}+S_{2}$, $T=T_{1}+T_{2}$, where
$S_{i}$ and $T_{i}$ are spin and orbital operators, respectively with $i=1,2$
indexing either orbital or spin. The same factor $1/2$ was already also
encountered earlier in a similar context with the operator ST on p.
LABEL:gls:X:ST.

###  5.5 Building rank-5 PEPS tensor

Projected entangled pair states (PEPS [53, 54, 30, 34]) can be used to tile a
tensor network state in one or higher dimension. If a set of tensors is used
to tile the lattice to infinity, this is referred to as iPEPS (or iTEBD in 1D
[55]). As with any lattice model, one needs to (i) define the local state
space of a site (in QSpace in terms of getLocalSpace), followed by (ii) an
initialization of the tensor network state as a whole. Assuming a 2D square
lattice, an $A$-tensor in iPEPS becomes a rank-5 tensor that combines the
virtual indices left (L), top (T), right (R), and bottom (B) with the local
state space of a site (s) [cf. Fig. 11].

Figure 11: Example for building rank-5 tensor for 2D PEPS, which ties the
virtual state space at the $L$(eft), $T$(op), $R$(ight), and $B$(ottom) with
the local state space of a site ($s$). Contracting the tensor decomposition in
(a) based on rank-3 tensors yields the rank-5 tensor in (b). Assuming periodic
tiling of a 2D-PEPS based on a unit cell comprised of a single tensor $A$, the
state spaces $L$ and $R$ must match, yet with reverse direction, and ditto for
$T$ and $B$.

For example, starting with a local state space described by an SU(3) spin in
the defining representation (10) $\equiv{\bf 3}$ [cf. Eq. (21)],

⬇

>> [S,IS]=getLocalSpace(’Spin’,’SU3’,’10’);

one may proceed to initialize the rank-5 tensor as sketched in Fig. 11(a).
Assuming, for simplicity, that the state spaces on all virtual indices are the
same, one nevertheless needs to make an initial choice for this state space.
To this end, one may build some arbitrary but fixed combination of symmetry
sectors derived from the local state space E,

⬇

>> E=IS.E; U=getIdentity(E,’-0’);

A2 = getIdentity(U,1,U,2);

E2 = getIdentity(U,1) + getIdentity(U,2) + getIdentity(A2,3);
?\mlbl{PEPS:E2:1}?

E2.data=repmat({eye(2)},size(E2.data)); ?\mlbl{PEPS:E2:2}?

⬇

E2 =

Q: 2x { 4 x 2 } having ’SU3’ { +- } ?\mlbl{E2:hdr1}?

data: 2-D double (544) 8 x 8 => 30 x 30 @ norm = ?$\surd$?30 ?\mlbl{E2:hdr2}?

\par1. 2x2 | 1x1 [ 00 ; 00 ] 32 b

2. 2x2 | 3x3 [ 01 ; 01 ] 32 b {?$\surd$?3}

3. 2x2 | 3x3 [ 10 ; 10 ] 32 b {?$\surd$?3}

4. 2x2 | 8x8 [ 11 ; 11 ] 32 b {?$\surd$?8}

Line LABEL:PEPS:E2:1 combines the local state space, its dual via the
$1j$-symbol U on index 2, and their fused state space via A2 into an identity
tensor E2. Since there is no overlap of symmetry sectors of the terms added
here, one can combine these spaces as a direct sum ($\oplus$) simply by adding
them up. Line LABEL:PEPS:E2:2 then sets the dimension to two multiplets per
symmetry sector by setting all RMTs in E2 to eye(2) $=11^{(2)}$. With the
state space E2 for all virtual indices assumed the same, one can combine two
of these, like $L$ and $T$, via the tensor $V$ as depicted in Fig. 11(a),

⬇

>> V=getIdentity(E2,E2);

X=getIdentity(V,3,E);

Since the fused state space in X as in Fig. 11(a) shall be contracted to the
conjugate of V (which then links to $B$ and $R$), their dimensions need to be
matched,

⬇

>> [Q3,dd]=getQDimQS(V,3); ?\mlbl{Q3:qdim}?

[I,J]=matchIndex(X.Q{3},Q3); ?\mlbl{Q3:match}?

for i=1:numel(I)

s=size(X.data{I(i)});

X.data{I(i)}=randn([s(1:2), dd(J(i))]); ?\mlbl{X:randn}?

end

Here getQDimQS in line LABEL:Q3:qdim returns the symmetry sectors (rows in Q3)
and their respective number of multiplets (integer vector dd) for leg $3$ in
V. Line LABEL:Q3:match finds all matching sectors in X on the third leg with
the entries in Q3, returned as I and J, respectively. The loop then matches
the dimension on the third leg of the respective RMT to the one in V via dd,
while randomizing their content using randn. With this, all ingredients in
Fig. 11(a) are defined and contractable. For transparency, one may add itags
now,

⬇

>> LT=setitags(V,{’L’,’T’,’LT’}); ?\mlbl{LT:def}?
BR=setitags(V,{’B’,’R’,’BR’}); ?\mlbl{BR:def}? X
=setitags(getsub(X,I),{’LT’,’s’,’BR’}); ?\mlbl{X:tags}?

The QSpace LT describes the upper tensor V in Fig. 11(a) that combines $L$ and
$T$. Similarly, BR describes the lower tensor V that combines $B$ and $R$ (up
to tensor conjugation which will reverse arrows in the contraction). Line
LABEL:X:tags selects the subset of records in X that matched BR (this step is
optional since non-matching symmetry sectors drop out in the contraction
anyway), and also assigns itags. Based on the assigned itags, the contraction
into the rank-5 tensor A5 in Fig. 11(b) can be carried out making use of auto-
contraction,

⬇

>> A5=contract({LT,X3},BR,’*’,[1 2 4 5 3]); ?\mlbl{A5:ctr}? d2=getDimQS(E2);
?\mlbl{A5:norm:d2}?

A5=A5*(d2(end)/norm(A5)); ?\mlbl{A5:norm}?

⬇

A5 = ?\mlbl{A5:disp}?

Q: 5x {85 x 2 } having ’SU3’ { L, T, B*, R*, s } ?\mlbl{A5:hdr1}?

data: 5-D double (39.26k) 8 x 8 x 8 x 8 x 1 => 30 x 30 x 30 x 30 x 3 @ norm =
30?\mlbl{A5:hdr2}?

\par1. 2x2x2x2 | 1x1x1x3x3 [ 00 ; 00 ; 00 ; 10 ; 10 ] 128 b

2. 2x2x2x2 | 1x1x3x3x3 [ 00 ; 00 ; 01 ; 01 ; 10 ] 128 b

3. 2x2x2x2 | 1x1x3x1x3 [ 00 ; 00 ; 10 ; 00 ; 10 ] 128 b

4. 2x2x2x2 | 1x1x3x8x3 [ 00 ; 00 ; 10 ; 11 ; 10 ] 128 b

5. 2x2x2x2 | 1x1x8x3x3 [ 00 ; 00 ; 11 ; 10 ; 10 ] 128 b

6. 2x2x2x2 | 1x3x1x1x3 [ 00 ; 01 ; 00 ; 00 ; 10 ] 128 b

: … ?\mlbl{A5:dots}?

84. 2x2x2x2x1 @10 | 8x8x3x8x3x10 [ 11 ; 11 ; 10 ; 11 ; 10 ] 1.2 k ?\mlbl{A5:rec84}?

85. 2x2x2x2x1 @10 | 8x8x8x3x3x10 [ 11 ; 11 ; 11 ; 10 ; 10 ] 1.2 k ?\mlbl{A5:recl}?

The trailing permutation in line LABEL:A5:ctr sets the index order LTBRs, as
also seen in the display in line LABEL:A5:hdr1. Lines
LABEL:A5:norm:d2–LABEL:A5:norm apply some normalization to A5 given the
earlier randomization in line LABEL:X:randn [getDimQS in line LABEL:A5:norm:d2
returns the multiplet and state space dimensions for all legs in E2, with
d2(end)$=30$ the state space dimension of the last leg; this stands for the
total number of states on the bond indices [LTRB] as described by E2; cf. line
LABEL:E2:hdr2]. The tensor A5 has multiplet dimensions 8 x 8 x 8 x 8 x 1 (line
LABEL:A5:hdr2) that correspond to a respective 30 x 30 x 30 x 30 x 3 states.
The 1 multiplet q=(10) on the last leg describes the local state space
$\sigma\in s$ with 3 states. It is the same throughout all 85 records which
reflect permissible symmetry combinations across all five legs. By default,
the full listing is truncated for readability (dots in line LABEL:A5:dots).
The display also shows two entries with largest outer multiplicity [$M=10$,
indicated by the trailing @10 with the RMT dimensions in lines
LABEL:A5:rec84–LABEL:A5:recl, as well as the trailing 6-th dimension x10 with
the respective CGTs]. They also happen to be the largest and last two records,
otherwise these would be shown as intermediate records.

###  5.6 fdm-NRG: Iterative diagonalization of Wilson chain and spectral
functions

The following more elaborate example demonstrates how to use QSpace for the
iterative diagonalization in the context of the numerical renormalization
group (NRG [18, 13]). It is included with this documentation since the NRG was
one of the first applications of QSpace [19, 50, 8]. The example here is based
on the standard single impurity Anderson model (SIAM). The following code
represents the script NRG/xnrgQS.m in the repository, yet with additional
comments and with a focus on auto-contractions based on itags. It reciprocates
the hard-coded MEX routine NRGWilsonQS. For transparency and the sake of the
discussion here, however, the iterative diagonalization can also be simply
coded in Matlab itself, while nevertheless utilizing the QSpace toolbox via
the MEX routines for tensor product, contraction, diagonalization, etc. It
also uses QSpace Matlab utility routines in the lib/ directory of the
repository, such as setdef, isset, add2struct, cto, etc., whose purpose is
simple and intuitive as will be indicated. A full explanation can be found in
their respective help. Since the following source mimics the algorithm behind
NRGWilsonQS, the format of the output data is also chosen in a compatible way.
This then permits, for example, subsequent usage of the MEX routine fdmNRG_QS
to compute correlation functions (see fdm-NRG below).

##### NRG basics

The NRG is a long-standing powerful impurity solver that tackles $0+1$
dimensional physical problems [18, 13]. For this, the impurity Hamiltonian is
mapped onto the so-called Wilson chain. It is described by the general
Hamiltonian

$\displaystyle\hat{H}$ $\displaystyle=$
$\displaystyle\hat{H}_{0}(\hat{\mathbf{d}},\hat{\mathbf{f}}_{0})+\sum_{n=0}^{\infty}\bigl{(}t_{n}\hat{\mathbf{f}}_{n}^{\dagger}\cdot\hat{\mathbf{f}}_{n+1}+{\rm
H.c.}\bigr{)}$ (70)

where $\hat{d}_{\sigma}\equiv(\hat{\mathbf{d}})_{\sigma}$ are the fermionic
operators for the impurity (‘$d$-level’ for historical reasons), where
$\sigma$ denotes spin. The fermionic operators for the bath
$(\hat{\mathbf{f}}_{n})_{\sigma}\equiv\hat{f}_{n\sigma}$ are organized into a
semi-infinite chain where the impurity couples to the first Wilson site
$\hat{\mathbf{f}}_{0}$ only. From a symmetry point of view, any term in the
Hamiltonian needs to represent a scalar. Hence the hopping term on the right
naturally includes the scalar dot product
$\hat{\mathbf{f}}_{n}^{\dagger}\cdot\hat{\mathbf{f}}_{n+1}\equiv\sum_{\sigma}\hat{f}_{n\sigma}^{\dagger}\hat{f}_{n+1,\sigma}$
implemented via the contraction in line LABEL:Hcpl:ctr below. The hopping
amplitudes $t_{n}\sim\Lambda^{-n/2}$ in Eq. (70) derive from the NRG coarse-
graining with its dimensionless discretization parameter $\Lambda\gtrsim 2$.
The resulting exponential decay with increasing $n$ is essential to the NRG,
as this justifies the iterative diagonalization. In this sense, the chain
structure in NRG does not represent real space, but rather energy scales
which, nevertheless, may be interpreted as inverse length scale [56, 57]. Due
to the effective one-dimensional character of the Hamiltonian in Eq. (70), NRG
is naturally suited to an MPS description. As a characteristic of NRG, the MPS
is not obtained variationally, though, but rather via iterative
diagonalization [58, 50]. Nevertheless, it is arguably in close proximity to a
variational approach [25, 26, 24].

##### Setup of model (SIAM)

The starting point of any NRG is the Hamiltonian H0 [$=\hat{H}_{0}$ in (70)]
that defines the impurity ($\hat{\bf d}$) together with its interaction with
the bath site at the location of the impurity ($\hat{\bf f}_{0}$). The matrix
elements of H0 are written in the basis encoded in the $A$-tensor A0. With
index order LRs in mind, it fuses the impurity ($d\to$ L) with the first
Wilson site ($f_{0}\to$ s) where H0 is expressed in R. Together this sets the
stage for the iterative diagonalization, fully defined by the fermionic
operators F (with Z) for a Wilson site and the hopping amplitudes ff. For the
SIAM, the impurity is described by a single spinfull fermionic level that
already shares the same structure as a Wilson site $\hat{\bf f}_{n}$.
Moreover, for the SIAM all interaction is contained within the impurity itself
via the Hubbard $U$, whereas the ‘interaction’ with the bath becomes a plain
fermionic hopping via the impurity’s hybridization function. Therefore for the
case of the SIAM, the NRG setup can be simplified in that H0 describes the
impurity ($\hat{\bf d}$) only. The bath-site $\hat{\bf f}_{0}$ can be already
taken care of as an iterative NRG step. With this the sum in Eq. (70) can be
extended to start from $n=-1$, with the additional hopping term taking care of
the hybridization of the impurity via
$t_{-1}\equiv\sqrt{\frac{2D\Gamma}{\pi}}$, with the $D:=1$ the half-bandwidth
of the bath. The latter sets the unit of energy, unless specified otherwise.
In any case, the iterative NRG diagonalization only requires the intial A0 and
H0, the fermionic operators F (with Z), and the couplings ff in the same units
as H0. This does not care about the precise starting point of the iterative
diagonalization which is also in the spirit of the MEX routine NRGWilsonQS.
For the SIAM, therefore the ‘left’ state space of A0 can be taken as the
vacuum state ($|\rangle\to$ L), and the local state space as the impurity
itself ($d\to$ s). Then H0 describes the impurity Hamiltonian only. The
hopping amplitude $\Gamma$ is prepended to the couplings as ff(1)$=t_{-1}$.
Since getNRGcoupling takes $\Gamma$ as input parameter, the hoppings are
already returned this way (line LABEL:nrg:ff). The setup for the SIAM thus
chooses the following adapted version of the Hamiltonian in Eq. (70),

$\displaystyle\hat{H}$ $\displaystyle=$
$\displaystyle\underbrace{\varepsilon_{d}\hat{n}_{d}+B\hat{S}_{z}+\tfrac{U}{2}\hat{n}_{d}(\hat{n}_{d}-1)}_{\quad\
\ \equiv H_{\rm imp}\ \to\ {\color[rgb]{0.45,0.45,0.45}\texttt{H0}}\ \
\text{(lines~{}\ref{SIAM:H0:hdr}--\ref{SIAM:H0:2})}}\quad+\underset{\text{(lines~{}\ref{NRG:iter:hdr}--\ref{NRG:iter:end})}}{\sum_{k=1}^{N}}\underbrace{\bigl{(}{\color[rgb]{0.45,0.45,0.45}\texttt{ff(k)}}\,\hat{\mathbf{f}}_{k-2}^{\dagger}\cdot\hat{\mathbf{f}}_{k-1}+{\rm
H.c.}\bigr{)}}_{\text{(lines~{}\ref{Hcpl:hdr}--\ref{Hcpl:end})}}$ (71)

with ff(k)$\equiv t_{k-2}$ and $\mathbf{f}_{-1}\equiv\mathbf{d}$.
In the script NRG/xnrgQS.m, the initialization of the SIAM reads,

⬇

\parsetdef(’U’,0.12,’Gamma’,0.01,’B’,0); ?\mlbl{nrg:sdef:1}?
setdef(’epsd’,-U/2); ?\mlbl{nrg:sdef:2}?
\par\parsetdef(’Lambda’,2,’N’,55,’Nkeep’,1024,’Etrunc’,7,’z’,0);
?\mlbl{nrg:sdef}?

\parglobal param

param=add2struct(’-’,U,epsd,Gamma,B,N,Lambda,z); ?\mlbl{nrg:param}?

\parff=getNRGcoupling(Gamma,Lambda,N,’z’,z); ?\mlbl{nrg:ff}? \parif
~isset(’B’)

SYM=’Acharge,SU2spin’;

else SYM=’Acharge,Aspin’; ?\mlbl{SIAM:SYM}? end

[FF,Z,SS,IS]=getLocalSpace(’FermionS’,SYM,’-v’);

\parN0=QSpace(size(FF)); nF=numel(FF);

for i=1:nF N0(i)=contract(FF(i),’13*’,FF(i),’13’); end N0=sum(N0);

The command setdef in lines LABEL:nrg:sdef:1, LABEL:nrg:sdef:2, or
LABEL:nrg:sdef sets default values for the specified variables if not yet set,
i.e., if they are not defined or empty. The command add2struct in line
LABEL:nrg:param adds specified variables as fields to the structure param
where the leading input argument '-' indicates to start from an empty
structure. This collects physical model and discretization parameters for
general information purposes later. The call to getNRGcoupling in line
LABEL:nrg:ff returns the Wilson chain hopping parameters ff, starting with
ff(1)$=\sqrt{2\Gamma/\pi}$. The above also sets the stage for the definition
of A0 and H0,

⬇

A0 = getIdentity(getvac(IS.E),IS.E,[1 3 2]); ?\mlbl{nrg:getvac}?

\par?\mlbl{SIAM:H0:hdr}?H0 = epsd*N0 + (U/2)*N0*(N0-IS.E) + 0*IS.E;
?\mlbl{SIAM:H0:1}?

if isset(’B’), H0 = H0 - B*SS(end); end ?\mlbl{SIAM:H0:2}?

The command getvac in line LABEL:nrg:getvac returns the vacuum state space for
the symmetry setting of the input. If a finite value for the magnetic field is
set, the setup via SYM in line LABEL:SIAM:SYM switches to abelian U(1) spin
symmetry. In this case, as can be checked by inspection of SS, the impurity
$S_{z}$ operator is returned as SS(3) $\equiv$ SS(end) in line LABEL:nrg:S3.
The magnetic field is applied as $-BS_{z}^{d}$ at the impurity only [it is
ignored for the bath, assuming $B\ll(D=1)$], where the minus sign ensures that
for $B>0$ one also obtains a positive magnetization $\langle
S_{z}^{d}\rangle$. The above concludes the setup of the impurity. Based on A0,
H0, and ff, this could proceed with the MEX routine NRGWilsonQS from here
[hard-coded C++ code, e.g., called via the rnrg (‘run NRG’) script which for
the SIAM calls setup/setupSIAM_SU2x2.m as a more elaborate version of the
setup above]. For the sake of transparency, however, the iterative
diagonalization is explicitly transcribed into the Matlab script xnrgQS
discussed here.

##### Iterative diagonalization

To proceed with the iterative diagonalization, the setup above is further
adapted. The couplings ff from line LABEL:nrg:ff above are in physical units
of the half-bandwidth ($D=1$). Since finite-size spectra should be compared in
rescaled units of order 1, the exponentially decaying couplings
ff$\sim\Lambda^{-k/2}$ are also rescaled to order 1,

⬇

rL=sqrt(Lambda); ff=ff.* rL.^(1:N-1);

Since at the first iteration ff(1) incorporates the first Wilson site ${\bf
f}_{0}$, this leaves the Hamiltonian H0 at iteration ‘0’ in its original
units. By default, the data from the iterative diagonalization is stored to
files (one file for each iteration, starting from *_00.mat for iteration ‘0’
which just diagonalizes the input H0),

⬇

setdef(’sflag’,1); ?\mlbl{sflag}? if sflag

fout={’NRG’ ’/myNRG’};

cto lma; ?\mlbl{nrg:cto}?if ~isdir([’./’ fout{1}]), system([’mkdir ’
fout{1}]); end

fout=[fout{:}];

wblog(’I/O’,’using end

Storing the data to files allows one to use the data later, e.g., for
computing spectral functions via fdmNRG_QS further below. The code proceeds
with initializing variables and data spaces,

⬇

k=0; AK=setitags(setitags(A0,’-A’,k),1,’Lvac’); ?\mlbl{nrg:AK:itags}?
HK=setitags(H0,’-op:K’,k); ?\mlbl{nrg:HK:itags}? AD=QSpace;
\parFKK=QSpace(1,nF);

\parEX=nan(N,2); E0=zeros(1,N); N4=[];
\parInrg=struct(’N’,N,’Lambda’,Lambda,’EK’,[],’HK’,QSpace(1,N));
?\mlbl{Inrg:init}?

\parfprintf(1,’\n’); disp(param); ?\mlbl{param:disp}?

fprintf(1,’\n’);

By iteratively keeping all eigenstates below some energy (Etrunc) or count
(Nkeep) threshold (cf. line LABEL:nrg:sdef), this partitions states into kept
(AK$\equiv\,A_{K}$) and discarded (AD$\equiv\,A_{D}$) states for any
iteration. Ultimately, the collection of AK for all iterations together with
the lowest energy eigenstate for the last iteration $k=N$ constitutes the
ground state MPS for given Wilons chain of length $N$. Expressing operator
matrix elements in terms of kept and discarded (based on AK and AD,
respectively), leads to sectors KK, KD, DK, and DD (e.g. see operator F $\to$
FKK $\equiv\,F_{KK}$). For scalar operators, there are no off-diagonal
contributions, hence a single label K or D suffices, such as for HK
$\equiv\,H_{K}$. By using itags as in lines
LABEL:nrg:AK:itags–LABEL:nrg:HK:itags this will simplify the specification of
contractions, but also characterize QSpace tensors more generally, while
making the output more readable and identifiable. Basic NRG data will be
collected into the info structure Inrg initialized in line LABEL:Inrg:init.
Line LABEL:param:dispf summarizes the setup by displaying the parameters used.
The following then proceeds with the actual iterative diagonalization,

⬇

?\mlbl{NRG:iter:hdr}?for k=1:N if k<N || ~sflag

o={’Nkeep’,Nkeep,’Etrunc’,Etrunc}; ?\mlbl{iter:Nkeep}?

else o={’Nkeep’,0}; ?\mlbl{iter:Etrunc}? end

\parq=getDimQS(HK); N4(k,:,2)=q(:,2);
?\mlbl{nrg:Dtot}?\par[ee,I]=eigQS(HK,o{:}); ?\mlbl{nrg:eig}?

\paree=ee(:,1); EK=QSpace(I.EK); ED=QSpace(I.ED); \parif k>1, E0(k)=min(ee);
ee=ee-E0(k); end

\parif EK Dk=getDimQS(EK); N4(k,:,1)=Dk(:,2); ?\mlbl{nrg:DK}?

EK=EK-E0(k); q=EK.data; EX(k,1)=max([q{:}]); end

if ED ED=ED-E0(k); q=ED.data;

EX(k,2)=min([q{:}]); end

\parif sflag AD=contract(AK,I.AD,[1 3 2]); end

AK=contract(AK,I.AK,[1 3 2]); if sflag

q=struct(’AK’,AK,’AD’,AD,’HK’,EK,’HD’,ED,’E0’,E0(k));

save(sprintf(’end

\parm=min(Nkeep,length(ee));

EE(1:m,k)=ee(1:m); Inrg.HK(k)=EK; \parif k==N, fprintf(1,’\n\n’); break; end

\parfor i=1:nF

FKK(i)=contract(AK,’!2*’,{AK,FF(i),’-op:^s’});

end

\parAK=getIdentity(I.AK,2,IS.E,[1 3 2]); AK=setitags(AK,’-A’,k);
\parHK=diag(EK*rL); HK=contract(AK,’!2*’,{HK,AK});

\par?\mlbl{Hcpl:hdr}? for i=1:nF, q=ff(k)*Z*FF(i);
Q=contract(AK,’!2*’,{FKK(i),’*’,{AK,q,’-op:^s’}}); ?\mlbl{Hcpl:ctr}?

HK=HK+Q+Q’; end ?\mlbl{Hcpl:end}?

\parHK=HK+0*getIdentity(AK,2); ?\mlbl{nrg:iter:done}?

\parfprintf(1,’ \parend ?\mlbl{NRG:iter:end}?

When adding a new Wilson site in lines LABEL:nrg:add:site–LABEL:nrg:iter:done
for the next iteration, the extended state space is exactly diagonalized in
line LABEL:nrg:eig. This also carries out the truncation into kept (K) and
discarded (D) states [cf. lines LABEL:I.EK and LABEL:I.ED]. Two parameters
govern this truncation as provided with the options o{:}, namely Etrunc
$\equiv E_{\rm trunc}$ and Nkeep $\equiv N_{\rm keep}$. This keeps at most
$N_{\rm keep}$ states (multiplets in the presence of non-abelian symmetries),
or fewer if the energy threshold $E_{\rm trunc}$ is reached earlier. Near-
degenerate eigenstates are never cut across by the truncation but kept
together (or truncated) in full to avoid artificial symmetry breaking of
symmetries that are not explicitly enforced. As per lines
LABEL:iter:Nkeep–LABEL:iter:Etrunc, both truncation parameters, $N_{\rm keep}$
and $E_{\rm trunc}$, are used for all iterations except for the last one,
where all states are discarded. Hence with no states kept, this way the Wilson
chain is naturally terminated [59]. NRG converges exponentially with
increasing $E_{\rm trunc}$ [58]. On empirical grounds, typically a value of
Etrunc $\equiv E_{\rm trunc}\gtrsim 7$ (in rescaled units as in the source
above) leads to good convergence. By truncating based on an energy threshold,
this mimics truncation at about a comparable discarded weight and hence
comparable accuracy along the Wilson chain [58]. Therefore it is typically
recommended to truncate based on $E_{\rm trunc}$, while keeping $N_{\rm keep}$
sufficiently large, yet set nevertheless as a safeguard on what is affordable
numerically. With this in mind, when increasing the complexity of the impurity