text:
\section{Introduction}

Quantum impurity models, describing a quantum system with a small
number of discrete states, the impurity, coupled to a continuous
bath of fermionic or bosonic excitations, arise in a variety of
contexts. A particularly important example is the Anderson impurity
model,~\cite{Anderson1961} relevant for describing magnetic moments in
metals, transport through quantum dots, and for the treatment of
correlated lattice models using dynamical mean field theory.
While the standard version of this model has \SU 2 spin symmetry in the
absence of a magnetic field,  generalizations to settings with
higher symmetry have also been studied.
 The \SU N generalization of the Anderson model emerged
first in the context of heavy Fermion systems,\cite{Bickers} where
large $N$ expansions proved to be an efficient way to model and
describe magnetic atoms with orbital degeneracy.  Studying these
models in detail is not only useful in the context of heavy fermion
systems, but it also represents the first step to understand the
behavior of correlated cold atomic gases with \SU N symmetrical
interactions. \cite{Honerkamp:2004hq, Rapp:2007gl}

The \SU N Anderson model can also be realized in a controlled way. In
particular, the \SU 4 model has been realized in various mesoscopic
structures including carbon nanotubes,\cite{Kouwenhoven,Finkelstein}
vertical quantum dots,~\cite{Tarucha} and more recently in the
originally proposed double dot
structures.~\cite{Davidunpublished} Similarly, the \SU 3 Anderson
model could also be realized with quantum dot structures, though the
proposed set-up is maybe somewhat more complicated. \cite{Carmi}

The \SU N Anderson model is defined in terms of $N$ local orbitals
embedded in a conduction electron sea. Its Hamiltonian can be written
in terms of the corresponding creation operators,
$d^{\dagger}_{\alpha}$ ($\alpha=1,\dots,N$) and the number operator
$\hat n= \sum_{\alpha} d^{\dagger}_\alpha d_\alpha$ as \begin{eqnarray} H &=&
\varepsilon_d\;\hat n + \frac U 2 \hat n (\hat n-1) + V \sum_\alpha
(d^{\dagger}_{\alpha} \psi_\alpha(0)
+ {\rm h.c.}) \nonumber \\
&& \quad+\; H_{\rm chan}[\psi,\psi^\dagger]\;.
\label{eq:H_SUN}
\end{eqnarray} Here $\varepsilon_d$ and $U$ denote the position of the local orbital and
the strength of interaction on it, respectively, and the level
$d^\dagger$ hybridizes locally with the fermions at its position,
destroyed by $\psi_\alpha(0)$. The last term of the Hamiltonian
describes the kinetic energy of the conduction electrons. It generates
the dynamics of the field $\psi_\alpha(0)$, and amounts in a
broadening of the "atomic" level, $\varepsilon_d$.

In the present paper we show on the prototypical example of the
  \SU 3 symmetrical Anderson model how the numerical renormalization
  group\cite{Wilson,Bulla} (NRG) of Ken Wilson, one of the most
  versatile and reliable tools for treating quantum impurity models,
  can be adapted to fully take advantage of non-Abelian symmetries to
  reduce computational costs.  Within Wilson's procedure, one
rewrites \eqref{eq:H_SUN} as Hamiltonian of a semi-infinite chain,
and diagonalizes it iteratively.~\cite{Bulla} Using symmetries in the
course of this diagonalization procedure is crucial: it allows
computer memory to be used efficiently, and enables one to reach
the required numerical accuracy on relatively standard computers with
reasonable runtimes. Eq.~\eqref{eq:H_SUN} obviously possess a
$\SU N \times \text{U}(1)$ symmetry corresponding to rotations in spin space
and overall charge conservation. Here, we shall focus on the $N=3$
case, classify states and observables while exploiting these
symmetries, and determine the spectral functions of several local
observables.

In an earlier work, a general framework has been set up and
implemented to handle an arbitrary number of non-Abelian symmetries
dynamically.\cite{Toth} This formulation allowed us to build an
open-access flexible density-matrix NRG (DM-NRG)
code.\cite{OpenAccess} However, in Ref.~\onlinecite{Toth} we considered only
combinations of certain rather simple symmetries such as charge
and spin \SU 2 symmetries, $\text{Z}_2$ or $\text{U}(1)$ symmetries.  In a group
theoretical sense, these are simpler than \SU{N>2} and some other
discrete or Lie groups. For \SU 2, irreducible representations (irreps)
are labeled by the size of the spin.
When "adding" two \SU 2 spins, say $S_1$ and $S_2$, each possible
total spin $S_3$ that satisfies the angular momentum addition rule
$|S_1 - S_2| \le S_3 \le S_1 + S_2$ can be obtained in precisely one
way. More technically, in the decomposition of the direct product
of two \SU 2 irreps labeled by spins $S_1$ and $S_2$ into a
direct sum of irreps, the number of times $n_{S_1 S_2; S_3}$ that
any irrep labeled by spin $S_3$ occurs in the direct sum, the so-called
``outer multiplicity'', is either 0 or 1.  This is, however, not
true for \SU{N>2}: the decomposition of the direct product
$\Gamma_1 \otimes \Gamma_2$ of two \SU N irreps into a direct sum
can contain irreps with outer multiplicity $n_{{\Gamma}_1 {\Gamma}_2 ;{\Gamma}_3}$
larger than 1, in other words, there may be $n_{{\Gamma}_1 {\Gamma}_2 ; {\Gamma}_3}$
inequivalent ways to construct the irrep $\Gamma_3$ from ${\Gamma}_1$ and
${\Gamma}_2$.  Correspondingly, the Clebsch-Gordan coefficients have a
more complicated structure in this case, and the Wigner-Eckart
theorem, extensively used in the DM-NRG calculations, becomes also
somewhat more complicated. Here we show how a general framework can be
constructed to deal with this case,\cite{footnote0} and demonstrate it
on the specific example of the \SU 3 Anderson model.

We note that another general approach towards exploiting non-Abelian
symmetries such as \SU N or Sp($N$) within NRG, and more generally for
tensor network methods, has recently been published by
A. Weichselbaum.~\cite{Andreas} It is formulated in the language of
matrix product states, and may be regarded as complementary to our
own, which is phrased within the more traditional formulation of
NRG. We emphasize, though, that both approaches are fully
  equivalent, in that precisely the same NRG assumptions and
  approximations are made in both; they differ only in the data
  structures used for internal book-keeping in the numerical
  codes. Their relation is briefly scetched in App.~\ref{append:MPS}.

Our DM-NRG calculations require explicit knowledge of Clebsch-Gordan
coefficients. Whereas these are known in closed form for \SU 2, this
is not the case for \SU{N>2}. However, an efficient numerical
algorithm for their evaluation has recently been
developed,~\cite{ArneJan} which we use here.

This paper is structured as follows: In \Sec{sec:nonAbelianSymmetries}
we outline our approach for exploiting non-Abelian symmetries in
DM-NRG calculations. In \Sec{sec:SU3-Anderson} we apply it to the \SU 3
Anderson model; in particular, we present results for the
conductance through a quantum dot described by this model, and for
various local spectral functions. Our conclusions in
\Sec{sec:conclusions} are followed by four appendices, that summarize
some basic facts of \SU N representation theory, and some recursion
formulas involving Clebsch-Gordan coefficients, respectively.

\section{DM-NRG with non-Abelian symmetries}
\label{sec:nonAbelianSymmetries}

As stated in the introduction, the formalism presented in
Ref.~\onlinecite{Toth} applies only for a special (though relatively
large) class of symmetries. Therefore, let us review here the most
important formulas and the structure of the NRG calculations for the
more general case, where so-called outer multiplicities are also
considered. \cite{typos}

\subsection{Local symmetries on the Wilson chain}

Let us start by first discussing the general structure of the symmetries of the Wilson
chain. The first step in Wilson's procedure of solving a quantum impurity problem is to
perform a Gram-Schmidt orthogonalization  and rewrite  the Hamiltonian in a ``tridiagonal''
form,\cite{Wilson}
\begin{equation}
{H} = {\cal H}_{0} +\sum_{n=0}^{\infty}\bigl( \tau_{n, n+1} +{\cal H}_{n+1}\bigr) \;.
\label{eq:H_Wilson_chain}
\end{equation}
Here ${\cal H}_{0}$ contains the local, interacting part of the Hamiltonian, while the rest of the chain represents the conduction electron (bath) degrees of freedom, coupled to it.
The on-site terms ${\cal H}_{n+1}$ are many times missing; they typically appear for more sophisticated electronic densities of states and can also account for superconducting correlations.
In our case, one could take, e.g.,
\begin{eqnarray}
{\cal H}_{0}  &=& \varepsilon_d\;\hat n + \frac U  2 \hat n (\hat n-1) + \tilde V
\sum_\alpha (d^{\dagger}_{\alpha} f^{[0]}_\alpha
+ {\rm h.c.}) \;,
\label{Wilson_H0}
\end{eqnarray}
as ${\cal H}_{0}$, with $f^{[0]}_\alpha\sim\psi_\alpha(0)$
a properly normalized on-site fermion,  and the hopping terms $\tau_{n, n+1}$ would read,
\begin{equation}
\tau_{n, n+1}= \sum_{\alpha} h^{[n]} \left(
 f^{[n]\dagger}_{\alpha} f^{[n+1]}_{\alpha} +{\rm h.c.} \right)\;.
\label{Wilson_Tau}
\end{equation}
Here $ f^{[n+1]}_{\alpha}$ annihilates  a fermion of \SU 3 spin $\alpha$ at site $[n]$, and
 the hopping amplitudes $h^{[n]}$ decay exponentially along chain,
thereby leading to energy scale separation.  Eq.~\eqref{eq:H_Wilson_chain} is
then diagonalized iteratively using the recursive relation
\begin{equation}
H_n= H_{n-1}+\tau_{n-1,n} + {\cal H}_n\;.
\label{eq:H_Wilson_recursive}
\end{equation}

In the following, we shall assume that $H$ (and $H_n$) are invariant
under the direct product of $n_S$ symmetry groups, \begin{equation} G = {\cal
  G}_{1}\times{\cal G}_{2}\times\dots\times{\cal G}_{n_S}\,.  \end{equation} This
means that for any group element $g_{\lambda}\in {\cal G}_{\lambda}$,
with $\lambda = 1, \dots, n_S$, there exists a unitary operator
${\cal U}_\lambda (g_{\lambda})$ on the Fock space, leaving $H$
invariant.  Here we do not need to make much restriction on the
groups: our considerations hold for any group which acts on the chain
locally at each lattice site $n$
\begin{equation} {\cal U}(g) =
\prod_{\lambda=1}^{n_S} {\cal U}_\lambda (g_{\lambda}) =
\prod_{\lambda=1}^{n_S}\prod_n {\cal U}_{\lambda,n} (g_{\lambda}) \,,
\end{equation}
and for which the Wigner-Eckart theorem holds.\cite{footnote1} Given
the above group structure, we can group all eigenstates and also
operators into {\it multiplets}, each of which transforms according to a certain
representation of $\cal G$,
\begin{equation} \Gamma \equiv \left\{{\Gamma}^1,
  {\Gamma}^2,...,{\Gamma}^{n_S}\right\}\leftrightarrow {\Gamma}^1\otimes {\Gamma}^2\otimes
\dots \otimes{\Gamma}^{n_S}\;.  \end{equation}
The "quantum numbers" $\Gamma^\lambda$, which label the
various irreducible representations (irreps) occurring in $\Gamma$, can
be spin labels, charges, or can label different
irreps of some point group.  States within a multiplet are
then labeled by a set of internal indices, $\gamma =\left\{ {\gamma}^1,
  {\gamma}^2,...,{\gamma}^{n_S} \right\}$, with the internal labels running from
$1\le {\gamma}^\lambda\le{\rm dim }({\Gamma}^\lambda)$. A given multiplet $i$
that transforms according to the representation $\Gamma_i$ thus
consists of \begin{equation} \dim (i)\equiv\dim
(\Gamma_i)=\displaystyle{\prod_{\lambda=1}^{n_S}}\dim ({\Gamma}_i^{\lambda})
\label{eq:dim_i}
\end{equation}
degenerate states.

States belonging to a product of two representations,
${\Gamma}_1^\lambda\otimes {\Gamma}_2^\lambda$, often appear in the
calculations. As outlined in the introduction, similar to spins, these
can be decomposed into irreps.
However, one irrep
may appear several times in this
decomposition,
\begin{equation}
{\Gamma}_s^\lambda\otimes {\Gamma}_p^\lambda \to ... \oplus
{\underbrace{{\Gamma}_q^\lambda \oplus ...\oplus {\Gamma}_q^\lambda}_{
n_{s,p;q}^{\lambda}
\mbox{ times}}}
\oplus ...\phantom{n}.
\label{eq:product_decomposition}
\end{equation}
Accordingly, in the most general case, the resulting states
$ {\Gamma}_s^\lambda\otimes {\Gamma}_p^\lambda \to \Gamma_q^\lambda$ must be labeled by a so-called "outer multiplicity" label, $\alpha_\lambda = 1,.., n_{s,p;\;q}^{\lambda}$. Correspondingly,
properly transforming multiplets  may be constructed from
the product of two multiplets as
\begin{equation}
| \Gamma, {\gamma}\rangle_\alpha =
\sum_{\lambda_1,\lambda_2}\left (   {\Gamma}_1,\gamma_1; {\Gamma}_2,\gamma_2| {\Gamma},\gamma   \right )_{\alpha }\;
|{\Gamma}_1,\gamma_1\rangle\otimes |{\Gamma}_2,\gamma_2 \rangle\;,
\label{eq:decomposition1}
\end{equation}
where $\alpha\equiv\{\alpha_\lambda\}$ denotes the composite multiplicity label, and the generalized Clebsch-Gordan coefficients are defined as
\begin{equation}
\left ( \Gamma_1, \gamma_1; \Gamma_2 ,\gamma_2 \mid  \Gamma,\gamma   \right )_{\alpha }
\equiv \prod_{\lambda=1}^{n_S}\left ( \Gamma_1^{\lambda},
\gamma_1^{\lambda}; \Gamma_2^{\lambda} ,\gamma_2^{\lambda} \mid  \Gamma^{\lambda},\gamma^{\lambda} \right )_{\alpha_{\lambda} }\;.
\label{eq:general_Clebsch}
\end{equation}

The outer multiplicity also appears in the Wigner-Eckart theorem. The
latter states that the matrix elements of an operator multiplet,
i.e. a set of operators $\{A_{\Gamma_A,\gamma_A}\}$ transforming under
transformations ${\cal U}(g)$ as a multiplet $\Gamma_A$, are determined
almost entirely by representation theory, and can be expressed in
terms of the Clebsch-Gordan coefficients as
\begin{eqnarray}
&&\langle i, \Gamma_i, \gamma_i  |  A_{\Gamma_A,\gamma_A}    | j, \Gamma_j,
\gamma_j  \rangle =
\nonumber\\
&&\quad\quad\sum_{\alpha}  \left ( \Gamma_A,
\gamma_A; \Gamma_j ,\gamma_j \mid  \Gamma_i,\gamma_i   \right )^{\ast}_{\alpha }\langle i\parallel A
\parallel j\rangle _{\alpha}\;.
\label{eq:WE}
\end{eqnarray}
Here multiplets $i$ and $j$ transform according to the
irreps, $\Gamma_i$ and $\Gamma_j$.
Thus, according to the Wigner-Eckart theorem, all matrix elements are
determined by only a few reduced matrix elements, $\langle i\parallel
A \parallel j\rangle _{\alpha}$, labeled just by the outer multiplicity
labels, $\alpha$ characterizing how many times the representation
$\Gamma_i$ appears in the product of $\Gamma_j$ and $\Gamma_A$. For many
commonly used symmetries as \SU 2, e.g., the outer multiplicity is
just always one and the label $\alpha$ can be dropped. However, it is needed
for, e.g., \SU{N\ge 3} or even for cubic point groups.

\subsection{Wilson's NRG with symmetries}

In course of the NRG procedure, one diagonalizes
Eq.~\eqref{eq:H_Wilson_chain} iteratively. The eigenstates of the
Hamiltonian $H_n$ of a chain of length $n$ can be grouped into
multiplets, with each multiplet $i$ transforming according to a
certain representation $\Gamma_i=\{\Gamma_i^1,..,\Gamma_i^{n_S} \}$.  Having
computed the approximate eigenstates (\emph{block states}) $|i, {\Gamma}_i,
{\gamma}_i, \rangle^{[n-1]}$ of $H_{n-1}$, one proceeds to construct
eigenstates of $H_{n}$.  To do that, one first appends to the chain
the multiplets $\{|\mu, {\Gamma}_\mu^{\rm loc}, {\gamma}_\mu^{\rm loc} \rangle \}
$, spanning \emph{local} Hilbert space at site $n$, and then
constructs properly transforming multiplets $\{| u, \Gamma_{u},
\gamma_{u}\rangle^{[n]}_{\alpha_{u}}\}$ by making use of the Clebsch-Gordan
coefficients, Eq.~\eqref{eq:decomposition1}:

 \begin{equation} | \mu, \Gamma_\mu^{\rm loc}, \gamma_\mu^{\rm loc}\rangle
\otimes | i, \Gamma_i,\gamma_i \rangle^{[n-1]}\to | u, \Gamma_{u},
\gamma_{u}\rangle^{[n]}_{\alpha_{u}}\;.  \end{equation}
Notice that a new multiplet  $u$ now also carries an outer multiplicity label, $\alpha_u$: This specifies the representation according to which $\Gamma_u$ has been produced from $\Gamma_i$ and $\Gamma_\mu$.
The advantage of using these states is that $H_n$ is diagonal both in $\Gamma_u$ and in the internal
labels, $\gamma_u$.  Therefore, it is sufficient to compute only the corresponding
irreducible matrix elements  $\langle u \parallel H_n\parallel v\rangle\big.^{[n]}$ in each
symmetry sector (block) separately, and diagonalize $H_n$ sectorwise by a unitary transformation
to obtain the  corresponding new eigenstates,
\begin{equation}
| u, \Gamma_{u}, \gamma_{u}\rangle ^{[n]}_{\alpha_{u}}  \to | {\tilde i}, {\Gamma}_{{\tilde i}}, {\gamma}_{{\tilde i}}  \rangle^{[n]}\;.
\label{eq:transformation}
\end{equation}
 As explained in App.~\ref{append:MPS}, this iterative procedure leads to a matrix product state (MPS) with a peculiar structure, reflecting the symmetry of the Hamiltonian.

The most difficult part in the procedure above is to determine  the matrix elements
$\langle u \parallel H_n\parallel v\rangle\big.^{[n]}$. These can be constructed by noticing that
each state $u$ and $v$ has been constructed from the eigenstates of $H_{n-1}$ and ${\cal H}_n$,
$i,\mu\to u$  and $j,\nu\to v$, and therefore
 \begin{equation}
 \langle u \parallel H_n\parallel v\rangle^{[n]} = \delta_{u,v}(E^{n-1}_i + \epsilon^{n}_\mu)
 + \langle u \parallel \tau_{n-1,n}\parallel v\rangle^{[n]} ,
 \nonumber
 \end{equation}
with $\epsilon^{n}_\mu$ being the  eigenenergy of ${\cal H}_n$. The matrix elements
of $\tau_{n-1,n}$ can be worked out by assuming  that the hopping part
consists of some fermionic or bosonic \emph{creation} operators $ C_{a,{\Gamma}_a,{\gamma}_a}^{[n]} $,
transforming again according to some
irreps ${\Gamma}_a$,
\begin{equation}
\tau_{n-1, n}= \sum_{a}
\left [
h_{a}^{[n-1]}\sum_{\gamma_a}\; C_{a,\Gamma_a,\gamma_a}^{[n-1]}(C_{a,\Gamma_a, \gamma_a}^{[n]})^\dagger +{\rm h.c.}
\right ].
\label{eq:hopping}
\end{equation} Here $a$ labels the different "hopping operators", and
$h_{a}^{[n-1]}$ the corresponding hopping amplitudes. Notice the
somewhat unusual way this hopping term is written:
$C_{a,\Gamma_a,\gamma_a}\leftrightarrow f^\dagger$ is a "creation operator",
which transforms according to the representation $\Gamma_a$, while
$(C_{a,\Gamma_a,\gamma_a})^\dagger\leftrightarrow f$ is an "annihilation
operator", transforming according to the conjugate representation,
${\Gamma_a}^*$.  We remark that for
charge \SU 2 symmetry, e.g., the "creation" operator multiplet
$C_{a,\Gamma_a,\gamma_a}$ is a Nambu spinor, and contains both $f$ and
$f^\dagger$ operators.\cite{Bulla1998} The number of hopping
operators may depend on the symmetry used: for a chain of spin 1/2
fermions treated in terms of $\SU 2 \times \text{U}(1)$ symmetry, e.g., one
has a single hopping operator of spin 1/2 and charge 1, while if only
the charge symmetry is used then one has two hopping operators of
charge 1, corresponding to the spin up and spin down directions.  In
our example of the \SU 3 Anderson model we have a single hopping
operator, and $C_{\Gamma_a,\gamma_a}^{[n]} \leftrightarrow \{
f^{[n]\dagger}_{\alpha}\}$.
 Assuming then that the reduced matrix elements of
 the creation operators acting on site
 $n-1$ of the chain,
$ \langle u \parallel C^{[n-1]}_a \parallel v\rangle_{\alpha}^{[n-1]}$,
   and those of the local creation operators at the added site,
   $\langle\nu\parallel C^{[n]}_a\parallel \mu\rangle_{\beta} $
  are known, one can use the Wigner-Eckart theorem to express  $\langle u \parallel \tau_{n-1,n}\parallel v\rangle^{[n]}$  as

 \begin{widetext}
\begin{multline}
\langle u\parallel\tau_{n-1,n}\parallel v\rangle\big.^{[n]}
=\delta_{\Gamma_{u},\Gamma_{v}} \displaystyle {\sum_{a,\alpha, \beta}}h_a^{[n-1]}\;
\langle i\parallel C^{[n-1]}_a \parallel j\rangle_{\alpha}^{[n-1]}\;
\langle\nu\parallel C^{[n]}_a\parallel
\mu\rangle_{\beta}^{\ast}\; D\left(a,  \alpha, \beta; u, v \right)
  + {\rm h.c.}\;, \phantom{nn}\\
  (u\leftarrow\mu,i; \; v\leftarrow \nu,j)
\;.
\end{multline}
\end{widetext}
 Here the outer multiplicity labels $\alpha$ ($\beta$) label inequivalent ways in which ${\Gamma}_i$ (${\Gamma}_\nu$)
appear in the product  ${\Gamma}_a\otimes {\Gamma}_j$ (${\Gamma}_a\otimes {\Gamma}_\mu$). The coefficients
$D\left(a,  \alpha, \beta; u, v \right)$ can be expressed in terms of Clebsch-Gordan coefficients, and are
given by Eq.~\eqref{eq:D} in Appendix~\ref{recursion_formulas}.
Similar expressions  hold for the matrix elements of "block operators" ($A$), acting somewhere  on
the first $n-1$ sites of the chain, and those of "local operators"   $A^{[n]}$,
acting only on the last site of the  chain.\cite{footnote2} For a block operator we
have, e.g.,
  \begin{equation}
\langle u \parallel A \parallel  v \rangle ^{[n]}_{\beta} =
\sum_{\alpha}
\langle i \parallel A \parallel   j \rangle^{[n-1]}_{\alpha}\;
 F\left (\alpha, \beta; u, v \right )\;\delta_{\mu, \nu}\; ,
 \label{eq:block_operator}
\end{equation}
while for the local operators the following equation holds,
\begin{equation}
\langle u \parallel A^{[n]} \parallel  v\rangle^{[n]}
_{\beta} =
\sum_{\alpha}
\langle \mu \parallel A^{[n]} \parallel  \nu \rangle _{\alpha}
\;
 K\left (\alpha, \beta; u,v \right )\;\delta_{i, j}\;,
\label{eq:local_operator}
\end{equation}
with the coefficients $F\left (\alpha, \beta; u, v \right )$ and $K\left (\alpha, \beta; u,v \right )$
given in Appendix~\ref{recursion_formulas}.  Here again, the outer multiplicity labels $\beta$  label inequivalent ways in which ${\Gamma}_u$  appears in ${\Gamma}_A\otimes {\Gamma}_v$, while $\alpha$
labels similarly inequivalent  ways how ${\Gamma}_i $ (${\Gamma}_\mu$) can be constructed
from   ${\Gamma}_A$ and ${\Gamma}_j$ (${\Gamma}_A$ and ${\Gamma}_\nu$).
Similar to $D$,
the coefficients
$F$ and $K$ are again determined only by symmetry, and can be expressed in terms of Clebsch-Gordan coefficients. As a last step of the iteration, the reduced matrix elements
$\langle u \parallel A^{[n]} \parallel  v\rangle^{[n]}_{\beta}$ and
$\langle u \parallel A \parallel  v\rangle^{[n]}_{\beta}$ need be transformed to the new basis,
$|{\tilde i}, {\Gamma}_{\tilde i}, {\gamma}_{\tilde i} \rangle $. This is performed by using precisely the same
unitary block transformations as the ones used to  diagonalize the Hamiltonian $H_n$,
\eqref{eq:transformation}, without affecting the outer multiplicity labels $\beta$.

Wilson's diagonalization procedure can be carried out then based upon the
equations above: In a given iteration, one
takes the lowest lying states of iteration $[n-1]$ and their matrix elements  $\langle i \parallel C_a^{[n-1]}\parallel j\rangle_\beta^{[n-1]}$, computes from these and from the matrix elements $\langle \mu \parallel C_a^{[n]}\parallel \nu\rangle_\beta$ the Hamiltonian $\langle u \parallel H_n\parallel v\rangle^{[n]}$. Then diagonalizing $\langle u \parallel H_n\parallel v\rangle^{[n]}$,  one obtains low-lying eigenstates of $H_n$ and determines their matrix elements  $\langle {\tilde i} \parallel C_a^{[n]}\parallel {\tilde j}\rangle_\beta^{[n]}$.
Continuing this procedure for larger and larger values $n$, one obtains better and better
approximations for the ground state of $H=H_{n\to\infty}$ and the lowest lying eigenstates.

\subsection{FDM-NRG approach}

So far, we discussed essentially Wilson's original NRG approach in
case of general symmetries. In practice, however, one often needs to
go beyond Wilson's RG and use the so called DM-NRG method,\cite{Hofstetter2000}
whereby a forward NRG run is first performed to obtain the density
matrix (DM) of the system, and then a backward NRG run is made to
compute physical observables. Moreover, to satisfy spectral sum rules,
a complete basis set\cite{AndersSchiller2005,Anders}  has to be used,
as first implemented in the context of DM-NRG
in Refs.~\onlinecite{Peters2006,Weichselbaum2007}.
In the full density-matrix NRG approach
(FDM-NRG) of Ref.~\onlinecite{Weichselbaum2007}, the
full density matrix of the entire chain is expressed
in the complete basis, which yields an improved treatment of
finite-temperature properties. Let us now
briefly discuss how symmetries can be implemented in the FDM-NRG
approach. (For a complimentary formulation of the same strategy
using matrix product states, see Ref.~\cite{Andreas}.)

First, to satisfy the necessary completeness relations, we consider a
chain of $N$ sites and introduce "environment" states $e$ for a each
state discarded in iteration $[n]$ ($i\in
D$),\cite{Peters2006,Weichselbaum2007} \begin{equation} |i, \Gamma_{i}\; \gamma_{i}
\rangle^{[n]}\rightarrow |i, \Gamma_{i}\; \gamma_{i}; e_n \rangle^{[n]}\;.
\label{complete}
\end{equation}
Here the states $e_n$ form an orthonormal basis for the remaining $N-n$ sites of the chain, and
their internal structure is irrelevant for the remaining discussion.

The states \eqref{complete} form a complete basis on the
Wilson chain,\cite{AndersSchiller2005,Anders} and can be used
to construct the density operator as follows,\cite{Weichselbaum2007}
\begin{eqnarray}
\varrho &=&\displaystyle \sum_{n=0}^{N}\; \varrho^{[n]}
\;,
\\
\varrho^{[n]} &\equiv& {\sum_{i\in D, e_n}} \sum_{\gamma_i}
\frac{e^{-\beta E_i^n}}{\cal Z}
|i, \Gamma_i,\gamma_{i}; e_n\rangle  ^{[n]}\, ^{[n]}
\langle i,\Gamma_i,\gamma_{i}; e_n |\;.\nonumber
\end{eqnarray}
Here $\beta =1/k_B T$ and the partition function is expressed as
\begin{eqnarray}
{\cal Z} =\displaystyle{\sum_{n=0}^{N}}\, \displaystyle{\sum_{i\in D}}\, \dim (i)\; e^{-\beta E_i^n}
d^{N-n}\;,
\end{eqnarray}
with $d$ the dimension of the Hilbert space at each added site of the Wilson chain, and $d^{N-n}$
the dimension of the space of the "environment" states, $e_n$. We remark that in the last iteration all states are
considered to be discarded, while in the first few iterations there are typically no discarded states yet.

To compute local observables and spectral functions of observables at the impurity site,
 one traces out step by step the
environment states, and introduces  the following set of truncated reduced density matrices,
\begin{equation}
R^{[n]}\equiv\underset{\left\{ e_n \right\}}{{\rm Tr\;}} \left\{\sum_{m \ge n} \varrho^{[m]}   \right\}\; .
\end{equation}
By symmetry, the reduced density matrices are invariant under the symmetries used, and have
therefore a blockdiagonal structure in the representation indices.\cite{Toth}
The matrices $R^{[n-1]}$ can be constructed iteratively from $R^{[n]}$ by tracing out site $n$ and then
adding the contribution of states  discarded  in iteration $(n-1)\to n$. The contribution of the kept states ($K$) reads:
 \begin{multline}
\langle i\parallel R^{[n-1]} \parallel  j \rangle^{[n-1]}_{i,j\in K} \\
=\displaystyle{\widetilde{ \sum_{{u, v},\mu}}}\;
\frac{\dim(u)}{\dim (i)}\;  \langle u \parallel {R}^{[n]} \parallel  v \rangle^{[n]}\;.
\end{multline}
Here the tilde indicates that the summation runs over states $u$ and $v$
having the same  symmetry, and
constructed  from states $i$ and $j$ by adding the \emph{same} local state,
$i\otimes \mu \to u$ and $j\otimes \mu \to v$.  The subscript indicates that $i$ and $j$ are both \emph{kept states }.
The discarded piece of $R^{[n-1]}$  is then simply
\begin{equation}
\langle i\parallel R^{[n-1]} \parallel  j \rangle^{[n-1]}_{i,j\in D}
= \delta_{i,j} \frac{d^{N+1-n}\; }{{\cal Z}}\;e^{-\beta E^{n-1}_i}.
\end{equation}

To gain insight to the dynamical properties of a quantum impurity, one usually computes the    retarded Green's functions for some operator multiplets  $A_{{\gamma}_A}$ and $B_{{\gamma}_B}$,

\begin{eqnarray}
G^{\rm ret}_{{{\gamma}_A},{{\gamma}_B}}
(\omega)
 \equiv  -\frac{i}{\hbar}\displaystyle{\int_{0}^{\infty}}\textrm{Tr}\Bigl\{{\varrho\;}\bigl[({ A}_{{\gamma}_A})^\dagger(t),{B}_{{\gamma}_B}(0)\bigr]_\xi\Bigr\}\;e^{i\omega t}dt\;,
 \nonumber
 \\
 \label{G}

\end{eqnarray}
with $\xi=-$ ($\xi=+$) referring to commutators (anticommutators) appearing for  bosonic (fermionic) operators. By symmetry, to have a non-vanishing value,
$A$ and $B$ must both transform according the same representation, ${\Gamma}_A\cong {\Gamma}_B$, and
${\gamma}_A\equiv {\gamma}_B$ must also be satisfied, $G^{\rm ret}_{{\gamma}_A,{\gamma}_B} = \delta_{{{\gamma}_A},{{\gamma}_B}}
G^{\rm ret}_{A^\dagger,B}$. Notice that Eq.~(\ref{G}) is defined in terms of
$({\hat A}_{{\gamma}_A})^\dagger$, transforming according to the conjugate
representation, ${\Gamma}_A^{\;\ast}$.\cite{footnote3}

The expression above can be evaluated in terms of the truncated  density matrices, $R^{[n]}$, and the reduced matrix elements
$\langle i \parallel A \parallel  j \rangle _{\alpha}^{[n]}$ and
$\langle i \parallel B \parallel  j \rangle _{\alpha}^{[n]}$
of the operators $A$ and $B$, to   obtain ~\cite{Anders,Toth}

\begin{widetext}
\begin{multline}
G^{\rm ret}_{A^\dagger,B}(z)=
\displaystyle{\sum^{N}_{n=0}}\;\displaystyle \sum_{i\in D,K}
{\sum_{(j,k)\notin(\rm K,\rm K)}}
\langle i \parallel R^{[n]}\parallel  j\rangle ^{[n]} \times
\left[
\sum_{\alpha}
\frac{\langle k \parallel A\parallel  j \rangle _{\alpha}^{[n]\ast}
\langle k  \parallel B\parallel  i\rangle_{\alpha}^{[n]}}
{z+\frac{1}{2}(E^n_i +E^n_j)-E^n_k}\frac{\dim(k)}{\dim(A)}\right .
\\
\left . -\xi\sum_{\alpha}   \frac{\langle j  \parallel B \parallel  k \rangle_{\alpha}^{[n]}
\langle i   \parallel A \parallel  k \rangle_{\alpha}^{[n]\ast}}
{ z-\frac{1}{2}(E^n_i +E^n_j)+E^n_k } \frac{\dim (i)}{\dim (A)} \right]\;.
\label{eq:green_function_explicite}
\end{multline}
\end{widetext}
This expression provides an efficient way to compute spectral
functions. Notice that it contains only the reduced matrix elements
and the dimensions of the multiplets involved.

\section{Study of the \SU 3 Anderson model}
\label{sec:SU3-Anderson}

To demonstrate how the scheme presented above works, we apply
it to study the repulsive
\SU 3-symmetrical Anderson model, defined already in the
Introduction. We perform our calculations for a conduction
band with a uniform local density of states between energies
$W>\epsilon>-W$ with the bandwidth set to $W\equiv1$, and use the
corresponding hopping amplitudes
$h^{[n]}\simeq(1/2) (1+\Lambda^{-1})\Lambda^{-n/2} $
in \eqref{Wilson_Tau}.
In this case, the width of the (noninteracting) level is approximately
given by
 ${\Delta} = \pi \varrho_c \tilde V ^2  $ with
$\varrho_c= 1/2W$ the local density of states at site $0$ of the Wilson
chain.
\begin{table}[h]
\begin{tabular}{|c|c|c|c|}
\hline
Operator & ($Q, F$) &  Dim & Components \\
\hline\hline
& & & \\
  &     &    &  ${\young (1) } \rightarrow f^{\dagger}_1 $\\
  & & & \\
 $f^{[n]\dagger}_\alpha$, $d^{\dagger}_\alpha$  & $\left( 1, \yng(1)\, \right ) $ & $1\times 3$ & ${\young ( 2)}\rightarrow f^{\dagger}_2$\\
 & & & \\
& & & ${\young ( 3)}\rightarrow f^{\dagger}_3$ \\
& & & \\
\hline
& & & \\
& & & ${\young(1,2)} \rightarrow d^{\dagger}_1 d^{\dagger}_2 $\\
& & & \\
$d_\alpha^{\dagger} d_\beta^{\dagger} $   &   $\left(2, \yng(1,1)\, \right)$  & $ 1\times 3$ &  ${\young
( 1,3)} \rightarrow d^{\dagger}_1 d^{\dagger}_3 $\\
& & & \\
& & & ${\young(2,3)} \rightarrow d^{\dagger}_2 d^{\dagger}_3 $\\
& & & \\
\hline
& & & \\
& & & ${\young(11,2)} \rightarrow d^{\dagger}_1 d_3 $\\
& & & \\
& & & ${\young (12,2)} \rightarrow d^{\dagger}_2 d_3 $\\
& & & \\
& & & ${\young (11,3)} \rightarrow - d^{\dagger}_1 d_2 $\\
& & & \\
& & & ${\young (12,3)} \rightarrow \frac{1}{\sqrt{2}}\left(d^{\dagger}_1 d_1-d^{\dagger}_2d_2\right) $\\
& & & \\
$d_\alpha^{\dagger} d_\beta$    &   $\left(0, \yng(2,1)\, \right)$  & $1\times 8$ &  ${\young (22,3)} \rightarrow d^{\dagger}_2 d_1 $\\
& & & \\
& & & ${\young (13,2)} \rightarrow \frac{1}{\sqrt{6}}\left(-d^{\dagger}_1 d_1-d^{\dagger}_2d_2 \right. $\\
& & & \\
& & & $\left . +2d^{\dagger}_3d_3\right)$\\
& & & \\
& & & ${\young (13,3)} \rightarrow -d^{\dagger}_3 d_2 $\\
& & & \\
& & & ${\young (23,3) } \rightarrow d^{\dagger}_3 d_1 $\\
& & & \\
\hline
\rule[-.5em]{0pt}{1.7em}
 $d_1^{\dagger}d_2^{\dagger}d_3^{\dagger}$  & $\left( 3, \bullet \right )$ &$ 1\times 1 $& ${ \bullet} \rightarrow d^{\dagger}_1d^{\dagger}_2d^{\dagger}_3$\\
\hline
\end{tabular}
\caption{Irreducible tensor operators for the \SU 3 symmetry.
\label{table:tensors}}
\end{table}

As mentioned before, the Hamiltonians \eqref{eq:H_SUN},
\eqref{Wilson_H0} and \eqref{Wilson_Tau} have a $\text{U}(1)\times \SU 3$
symmetry in the charge and flavor sectors, respectively.
Correspondingly, multiplets of the Hamiltonian
are characterized by a charge and a flavor quantum
number.  The charge quantum numbers $Q_i$ are simply
identical to the total charge,~\cite{footnote_charge}
\begin{equation} Q \equiv
\sum_{\alpha}\left \{\sum_{n=0}^{N}\left( f^{[n]\dagger}_{\alpha}
  f^{[n]}_{\alpha}-{3 \over 2}\right) +
  \left (d_\alpha^{\dagger}\,d_\alpha -{3\over 2}\right ) \right \},
\end{equation}
  conserved by
\eqref{eq:H_SUN}, \eqref{Wilson_H0}, and \eqref{Wilson_Tau}.

Labeling the \SU 3 representations and the states within an
\SU 3 multiplet is somewhat more complicated. The flavor quantum
numbers $F_i$ can be (and are usually) represented by Young tableaux,
characterized by two non-negative integers in case of \SU 3 (see
Appendix~\ref {appendix:SU3}). Young tableaux provide a nice
pictorial way to multiply and decompose representations, or calculate
their dimensions.  For our numerical calculations, however, we had to
construct explicitly the basis states of \SU 3 representations and to
compute the corresponding Clebsch-Gordan coefficients.\cite{footnote4}
This we carried out using the so-called Gelfand-Tsetlin patterns,
briefly discussed in
Appendix~\ref{appendix:SU3}. The states obtained this way are analogues of
the canonical \SU 2 basis states, created by the raising and lowering
spin operators, $S^\pm$.  Gelfand-Tsetlin patterns are in one-to-one
correspondence with the Young tableaux, but they allow for a simpler
explicit construction of the basis states.  For more details, we refer
the reader to Ref.~\onlinecite{ArneJan}.

\begin{figure}
  \begin{center}

      \includegraphics[width=0.9\columnwidth, clip]{Fig1.eps}

   \end{center}
\caption{\label{fig:occupation_conductance_T_0}
    (Color online) The occupation number (upper panel)
    and the conductance (lower panel) as function of $\varepsilon_d$ and for
	different broadening parameters ${\Delta}$. The temperature is $T=0$ and the
	Coulomb interaction is fixed to $U = 0.2\, W.$
	\label{fig:occupation_conductance}
}
\end{figure}

Similar to the  eigenstates and multiplets, we also need to group operators into \SU 3 multiplets and characterize them by appropriate \SU 3 quantum numbers.
This classification of the most important operators is summarized in Table~\ref{table:tensors}.
In  terms of \SU 3, there is only a single hopping operator,
$C^{[n]}_\gamma \leftrightarrow  {f^{[n]\dagger}_\alpha}$, which transforms according to the defining \SU 3 representation and has charge $Q=1$, similar to the creation operator of the localized level, $d^\dagger_\alpha$.  From $d^\dagger_\alpha$ we can also construct various local operators of interest. The  spin operators, $\sim d^\dagger_\alpha d_\beta$ form, e.g., an 8-dimensional charge $Q=0$
operator multiplet  in terms of $\text{U}(1)\times \SU 3$, while the charge $Q=2$ pairing
operators $\sim d^\dagger_\alpha d^\dagger_\beta$
transform according to a 3-dimensional \SU 3 representation. Finally, the "trion"
operator $d^\dagger_1d^\dagger_2d^\dagger_2$ has charge $Q=3$ and is an \SU 3 singlet.

\subsection{Numerical results}

In our runs we have kept about 250 multiplets at each
iteration, corresponding to approximately 1500 states on average, while
the Wilson parameter was fixed to $\Lambda =2$.

\begin{figure}[t]
  \includegraphics[width=0.98\columnwidth]{Fig2.eps}
  \caption{\label{fig:flow_diagram_half_filling}
 NRG finite size spectrum at half-filling and T=0.
 Upper (lower) panel represents the even
 (odd) part of the spectrum. In both plots we represent the lowest 50 energy levels. The parameters are: $\Lambda = 2$, $U = 0.2 \;W$, $\varepsilon_d = -U$ and
 $U/\pi{\Delta} = 5$. The convergence was reached after approximately 27 iterations.
 \label{fig:flow_diagram} }
\end{figure}

In Fig. \ref{fig:occupation_conductance}(a)  we
present the occupation of the localized level as a function of the
energy $\varepsilon_d$.  The localized level can accommodate up to 3 fermions.
For large enough $U/\Delta$'s, these fermions enter the local level one
by one as $\varepsilon_d$ is decreased, and a
Coulomb staircase is clearly visible. In the range $-U<\varepsilon_d<0$ there
is approximately one electron on the level. The \SU 3 spin of this
electron is screened by the conduction electrons, and an \SU 3 Kondo
state is formed. A similar, hole-like $\SU 3^*$ state emerges for
$-2U<\varepsilon_d<-U$.

\begin{table}[t]
Even iterations\vspace{0.3cm}
\begin{tabular}{|c|c|c|c|c|}
\hline
$E/({2\pi/L})$ & $E_{NRG}/ E_L$ &  $Q$ (charge) & \SU 3 & dimension\\
\hline\hline
0    &   0   &   $Q_0$ &  \textbullet &  1\\
\hline
& & & &  \\
$1/2$ & $0.5$ &  $Q_0+1$ &$ \yng (1)$ & 3\\
& & & &  \\
$1/2$ & $0.5$ & $Q_0-1$ &$ \yng (1,1)$ & 3\\
& & & &  \\
\hline
& & & &  \\
$1 $ & $1$ &  $Q_0+2$ &$ \yng(1,1)$ & 3\\
& & & &  \\
$1$ &  $1$ & $Q_0$  &$ \yng(2,1)$ & 8\\
& & & &  \\
$1$ & $1$ & $Q_0 $ & \textbullet & 1\\
& & & &  \\
$1$ & $1$ & $Q_0-2$ &$ \yng(1)$ & 3\\
& & & &  \\
\hline

\hline
\end{tabular}\\
\vspace{0.3cm}
Odd iterations\vspace{0.3cm}
\begin{tabular}{|c|c|c|c|c|}
\hline
$E/({2\pi/L})$ & $E_{NRG}/E_L$ &  $Q$ (charge) & \SU 3 & dimension\\
\hline\hline
0    &   0   &   $Q_0$ &  \textbullet &  1\\
0    &   0   &   $Q_0$ - 1&  $\yng(1,1) $ &  3\\
& & & &  \\
0    &   0   &   $Q_0$ - 2&  $\yng(1) $ &  3\\
0    &   0   &   $Q_0$ - 3&  \textbullet &  3\\
\hline
& & & &  \\
$1 $ & 0.98 & $Q_0$ + 1&$ \yng(1,1)$ & 3\\
& & & &  \\
$1$ & 0.98 & $Q_0$  &$ \yng(2,1) $ & 8\\
$1$ & 0.98 & $Q_0$  & \textbullet & 1\\
$1$ & 0.98 & $Q_0$ - 1 &$ \yng(1,1)$ & 3\\
& & & &  \\
$1$ & 0.98 & $Q_0$ - 1 &$ \yng(2)$ & 6\\
& & & &  \\
$1$ & 0.98 & $Q_0$ - 1 &$ \yng(1,1)$ & 3\\
& & & &  \\
$1$ & 0.98 & $Q_0$ + 1 &$ \yng(2,2)$ & 6\\
& & & &  \\
$1$ & 0.98 & $Q_0$ + 2 &$ \yng(1)$ & 3\\
& & & &  \\
\hline
\end{tabular}
\caption{ Energy spectrum for the electron-hole symmetric, mixed valence
point. $E_L$ stands for the finite size energy scale, $E_L = 2\pi/L = 1.311$.
$Q_0 =0\, (3/2) $ for even (odd) iterations.
\label{table:spectrum_mixed_valence}}
\end{table}

The most intriguing region is though the mixed valence
region, $\varepsilon_d\approx -U$. For $\varepsilon_d=-U$ the ground state of the
isolated impurity ($\Gamma=0$) would be sixfold degenerate due to
electron-hole symmetry,  connecting the two \SU 3 triplets,
$\{ d^\dagger_\alpha|0\rangle\}$ and $\{ d^\dagger_\alpha d^\dagger_\beta |0\rangle\}$.
 Valence fluctuations produce a state with
all these states strongly mixed by quantum fluctuations.
Fig.~\ref{fig:flow_diagram} and Table \ref{table:spectrum_mixed_valence} show the
flow diagram of the NRG levels, and the \SU 3 classification of the asymptotic
finite size spectrum, respectively. A detailed analysis reveals that this finite size spectrum
can simply be understood as the finite size spectrum of three
chiral fermions with a phase shift $\delta=\pi/2$ at the Fermi energy.
This phase shift is indeed in  full agreement with the Friedel
sum rule, $3\delta/\pi =\langle n\rangle $, and the occupation $\langle n\rangle=3/2$ assured by
electron-hole symmetry.

Similarly, the \SU 3 Kondo spectrum, displayed in Table \ref{table:spectrum_kondo}
can be understood as the finite size spectrum of three
chiral fermions with a phase shift $\delta\approx\pi/3$, implied by the Friedel sum rule
and the occupation $\langle n\rangle\approx 1$.

\begin{table}[t]
Even iterations\vspace{0.3cm}
\begin{tabular}{|c|c|c|c|c|}
\hline
$E/({2\pi/L})$ & $E_{NRG}/E_L$ &  $Q$ (charge) & \SU 3 & dimension\\
\hline\hline
0    &   0   &   $Q_0$ &  \textbullet &  1\\
\hline
& & & &  \\
$2/3 $ & 0.33 & $Q_0$ - 1 &$ \yng(1,1)$ & 3\\
& & & &  \\
\hline
& & & &  \\
$4/3$ & 0.66 & $Q_0$ + 1 &$ \yng(1)$ & 3\\
& & & &  \\
$4/3$ & 0.66 & $Q_0$ - 2 &$ \yng(1)$ & 3\\
& & & &  \\
\hline
& & & &  \\
$2 $ & 1 & $Q_0$  &$ \yng(2,1)$ & 8\\
$2 $ & 1 & $Q_0$  & \textbullet & 1\\
$2 $ & 1 & $Q_0$ - 3 & \textbullet & 1\\
& & & &  \\
\hline
\end{tabular}\\
\vspace{0.3cm}
Odd iterations\vspace{0.3cm}
\begin{tabular}{|c|c|c|c|c|}
\hline
$E/({2\pi/L})$ & $E_{NRG}/E_L$ &  $Q$ (charge) & \SU 3 & dimension\\
\hline
0    &   0   &   $Q_0$ &  \textbullet &  1\\
\hline
& & & &  \\
$ 1/3$ & $0.15$ &  $Q_0+1$  & $ \yng(1) $ & 3\\
& & & &  \\
\hline
& & & &  \\
$2/3$ & 0.32 & $Q_0 + 2$ & $ \yng(1,1)$ & 3\\
& & & &  \\
\hline
\end{tabular}
\caption{ Energy spectrum in the Kondo regime. $Q_0 = 0\,(-3/2)$ for
even (odd) iterations.
\label{table:spectrum_kondo}}
\end{table}

The crossover between the two \SU 3 Kondo regimes through the mixed valence regime is
maybe best captured by the local level's spectral function, shown in Fig.~\ref{fig:spectral_function_f_dagger}. In the \SU 3 and $\SU 3^*$ Kondo regimes we find
a  narrow Kondo resonance of an exponentially small width pinned somewhat
asymmetrically to the Fermi energy, $\omega=0$ (see inset), and two Hubbard peaks.
At the mixed valence point, $\varepsilon_d=-U$
on the other hand, a relatively broad and symmetrical resonance of width $\sim \Gamma$ appears
at the Fermi energy, and the charging peaks are absent.

\begin{figure}[t]
\includegraphics[width=0.93\columnwidth, clip]{Fig3}
\caption{\label{fig:spectral_function_f_dagger}
    (Color online) a) The normalized spectral function for the
    on-site creation operator $d^\dagger_\alpha$ for two different fillings $n = 1.5$
    (half filling, red dashed) and $n=1$ (1/3 filling, black solid).  In the singly occupied case, $n=1$,
    the dot is in the Kondo regime, and the resonance is shifted away from the Fermi
    energy, $\omega=0$. The value of the spectral function at the Fermi energy is determined by the Friedel sum rule, $\pi {\Delta} A(0)\to 3/4$.
    b) Evolution of the normalized spectral function
    as function of $U/\pi\Delta$
    in the Kondo regime, for    $\varepsilon_d /U = -0.5$. The deviation from the Friedel sum rule is
    less than 1\%.}
\end{figure}

\begin{figure}[t]
    \begin{center}
       \includegraphics[width=0.6\columnwidth, clip]{Fig4.eps}
    \end{center}
\caption{
\label{fig:sketch}
(Color online)  Possible realization of \SU 3 Kondo states using four
    edge states and three quantum dots. The upper edge state is splitted and
    allows transport measurements between the two upper external leads.
    \cite{Carmi} }
\end{figure}

\subsection{Conductance}

Let us now consider the mesoscopic structure proposed in Ref.~\onlinecite{Carmi}, sketched in
Fig.~\ref{fig:sketch}.   As explained in Ref.~\onlinecite{Carmi}, this quantum dot structure would
possibly realize a Kondo state with an approximate \SU 3-symmetry.
The conductance through the device  can be directly related to the spectral functions of the
$d$-level,  and for a symmetrical device one finds \begin{equation}
G(T) = \pi \Delta \;{e^2 \over h } \int_{-\infty}^{\infty}{\rm d}\omega\;
A(\omega, T) \left (-\frac{\partial f(\omega)}{\partial \omega} \right)\;.
\label{eq:conductance}
\end{equation}

Here $A(\omega)$ is the spectral function of the $d_\alpha$ operator of the localized
level, $ A(\omega) = -\Im m\, G_{d_\alpha\,d_\beta^\dagger} (\omega)/\pi $, with
$G_{d_\alpha\, d_\beta^\dagger} (\omega)$ the Green's function defined in Eq.~\eqref{G}.
The corresponding $T=0$ temperature linear conductance is shown in
Fig.~\ref{fig:occupation_conductance}(b)
 as a function of the energy $\epsilon_d$.
The conductance reaches its maximal value in the mixed valence regime,
$\varepsilon_d\sim -U$, and displays \SU 3 Kondo effect related plateaus in the regions
$-2U<\varepsilon_d < -U$  and $-U<\varepsilon_d < 0$. The $T=0$ temperature conductance values observed
can be understood in terms of the Friedel sum rule, relating the total occupation of the $d$-level to the phase shift of the conduction electrons,
which yields $3 \delta/\pi = \langle n\rangle$ for the
\SU 3 Anderson model.~\cite{Carmi}
In a Fermi liquid state, --- at $T=0$ temperature, ---  the conductance can be computed using the Landauer-B\"uttiker formula, and in the geometry
of Fig.~\ref{fig:sketch}
is simply related to the phase shift as $G(T=0)=(e^2/h) \sin^2(\delta)$. This explains the value
$G_{\SU 3}\approx (3/4)\; e^2/h$ observed in the Kondo states; there the occupancies
are  $\langle n\rangle\approx 1$ and $\langle n\rangle\approx 2$, corresponding to phase shifts
$\delta=\pm \pi /3$, and the previously mentioned value of the conductance.
 At the mixed valence point,
$\varepsilon_d=-U$, on the other  hand, we have  $\langle n\rangle = 3/2$, implying a phase shift $\delta=\pi/2$, and a maximal conductance, $G_{\SU 3} = e^2/h$ .

\begin{figure}[b]
\begin{center}
\includegraphics[width=1\columnwidth, clip]{Fig5.eps}
\end{center}
\caption{   \label{fig:conductance_finite_T}
 (Color online)  Finite temperature conductance as function of
    $\varepsilon_d$. Inset: The temperature dependence of the conductance
    on a logarithmic scale,	for two different fillings,
    $n=1.5$ (half filling, red dashed) and $n=1$ (1/3 filling, black solid).
    In both panels the Coulomb energy and the broadening to the
    leads were fixed to $U = 0.2\,W$ and $U/\pi{\Delta} = 5$.}
\end{figure}

\begin{figure}
      \includegraphics[width=0.8\columnwidth, clip]{Fig6a.eps}
      \vskip0.2cm
      \includegraphics[width= 0.8\columnwidth, clip]{Fig6b.eps}
       \caption{
      \label{fig:spin}
     (Color online) Upper panel: Spectral function of the spin operator $d^\dagger_\alpha d_\beta$ on a
    logarithmic scale for two different fillings, $\langle n\rangle= 1.5$
    (half filling, red dashed) and $\langle n\rangle \approx 1$ (1/3 filling, black solid) for
    $U = 0.2\,W$.
    In the singly occupied case, $\langle n\rangle \approx 1$,  the dot is in the Kondo regime. The inset indicates the linear decay
    in $A_{S}(\omega)$ in the small frequency limit.
    The blue, dashed-dotted line is a guideline for the eye.
    Lower panel: The spectral function of the spin operator
    in the Kondo regime, $\langle n\rangle \approx 1$, and for three different ratios $U/\pi{\Delta}$: $U/\pi{\Delta}=15$ (solid black line), $U/\pi{\Delta}=10$ (dashed red line) and
      $U/\pi{\Delta}=5$ (dashed-dotted green line). The inset indicates the
    universal scaling  collapse of the spin spectral function in
    the Kondo regime.
      }

\end{figure}

Increasing the temperature, the conductance is quickly suppressed in
the Kondo regimes, and three Coulomb blockade conductance peaks emerge
at the points of charge degeneracy, as shown in
Fig.~\ref{fig:conductance_finite_T}.  The central peak corresponds to
the transition $n=1\leftrightarrow 2$ while the two side peaks
correspond to charge fluctuations $n=1\leftrightarrow 0$ and
$n=2\leftrightarrow 3$, respectively.  In
Fig.~\ref{fig:conductance_finite_T} we also display\scrap{ed} the
temperature dependence of the conductance at the mixed valence point
and in the \SU 3 Kondo regime. The Kondo temperature of the \SU 3
Kondo state is clearly much smaller than the mixed valence energy
scale even for these moderate interactions.  This Kondo temperature can be estimated
by first doing perturbation theory in $\tilde V$ and performing a Schrieffer-Wolff transformation and
then carrying out a renormalization group analysis. This analysis yields a Kondo temperature
of
\begin{equation}
T^{\SU 3}_K \approx D_0 \; e^{-1/3 \lambda} \;,
\label{eq:TK}
\end{equation}
with the dimensionless coupling $\lambda$ expressed as
\begin{equation}
\lambda = \frac{\Delta}{\pi E_+}  + \frac{\Delta}{\pi  E_-}\;
\end{equation}
in terms of the level width $\Delta$ and the  "ionization energies"
$E_+ = U+\epsilon_d$ and $ E_- =- \epsilon_d$, and
$D_0\approx \min (E_+, E_-)$ a high energy cut-off.

\begin{figure}[b]
     \includegraphics[width= 0.8\columnwidth, clip]{Fig7.eps}
\caption{
      (Color online)
    Spectral function of the superconducting operator for
    $U = 0.2\,D$ and $\langle n\rangle= 1.5$
    (half filling, red dashed) and $\langle n\rangle \approx 1$ (1/3 filling, black solid).
          \label{fig:SC}
    }
\end{figure}

\subsection{Spectral functions of local operators}

The behavior of strongly correlated
cold atomic and heavy fermion lattice systems can often be understood in terms of a self-consistent
quantum impurity model (dynamical mean field theory). Within this picture,  the local response functions of the quantum impurities  may drive superconducting or magnetic phase transitions, or lead to even more exotic quantum phases.
In this subsection, let us therefore analyze the spectral properties
of the \SU 3 Anderson model, and investigate
  the local spectral and response functions
  of  its spin, pairing, and trion operators.

\begin{figure}[b]

      \includegraphics[width= 0.9\columnwidth, clip]{Fig8.eps}
\caption{
     The spectral function for the trion operator
    for two different fillings $n = 1.5$
    (half filling, red dashed) and $n=1$ (1/3 filling, black solid).
    The inset represents the same data on a
    logarithmic scale. The trionic spectral function
    scales as $\omega^2$ in the small energy limit.
    \label{fig:trion}
    }
\end{figure}

As discussed before, the spin operators, $d^\dagger_\alpha d_\beta$ transform according to an 8-dimensional \SU 3 representation. As shown in  Fig.~\ref{fig:spin},
their spectral function displays Fermi liquid properties (see inset), and behaves very similarly to the spin spectral function of a standard \SU 2
Anderson model.\cite{Hewson2006} In the mixed valence regime, for ${\Delta}\ll U$
charge fluctuations to the state $n=0$ and $n=3$ are frozen out, and at low energies
the only relevant energy scale is  ${\Delta}$; correspondingly, the spectral function exhibits a broad resonance at $\omega\sim {\Delta}$ (extending up to $\omega\sim U$), and decays linearly to zero for small frequencies, $A_S^{n\approx 3/2}(\omega)\sim \omega/{\Delta}^2$.
By Hilbert transform, this amounts in a spin
susceptibility, $\chi_S\sim 1/{\Delta}$. In the Kondo regime, $\langle n\rangle\approx 1$, on the other hand,
two separate scales can be distinguished. Below $\omega\sim \min(U,|\epsilon_d|)$ charge fluctuations are frozen and a clear Kondo anomaly can be observed as a logarithmic increase of the spectral
function, as indicated by the arrow in Fig.~\ref{fig:spin}.
The Fermi liquid behavior only emerges below  the Kondo scale, $\omega<T_K\ll {\Delta},U$.
In this Kondo regime $A_S^{\rm Kondo}(\omega)\sim \omega/T_K^2$, and correspondingly,
a susceptibility $\chi_S\sim 1/T_K$ is found.

The spin spectral function becomes universal in the Kondo limit, $T_K\ll \Delta,U$ in the sense that
the $T=0$ temperature dynamical susceptibility scales as
 $\chi_S(\omega)=  (1/T_K)\; f(\omega/T_K,\langle n\rangle)$, with $f(\omega/T_K,\langle n\rangle)$
a function, which only slightly depends on the occupation of the level,  $\langle n\rangle\approx 1$.
This is demonstrated in the lower panel of Fig.~\ref{fig:spin} for the imaginary part of the susceptibility,
$\chi^{\prime\prime}_S(\omega) =  - \pi A_S(\omega)$, computed for different values of $U/\pi{\Delta}$.
Numerically we define $T_K$ as the half width at half maximum of $A(\omega)$, the
spectral function of the $d$ operator of the localized level.

The correlations of the pairing operators,  $d^\dagger_\alpha d^\dagger_\beta$, behave somewhat similar to those of the spin in the sense that at small frequencies a linear frequency dependence is found, corresponding to a Fermi liquid state with a constant pairings susceptibility (see inset of Fig.~\ref{fig:SC}). However, the local pairing operators  lead out of the local low energy
charging states both in the mixed valence and in the Kondo regimes. Therefore, the
amplitudes of their spectral functions as well as the
corresponding pairing susceptibilities
obtained by a Hilbert transform,
are typically small, and they are not expected to drive any transition.

As a final example, we display the spectral function of the \SU 3
singlet trion operator, $T^\dagger =
d^\dagger_1d^\dagger_2d^\dagger_3$ in Fig.~\ref{fig:trion}. The trion
operator plays an important role in the attractive
case,\cite{Klingschat2010} however, in this repulsive model it
is a highly suppressed operator. In the Kondo regimes, it has a
non-zero spectral function only because the Kondo states have a small
 ($\sim {\Delta}/U$) admixture of the empty and the triply occupied
states, respectively. This explains why the amplitude of the signal is
in this case smaller in the mixed valence regime. It also explains the
strong electron-hole asymmetry in the Kondo regimes. In the $\langle
n\rangle \approx 1$ regime, shown in Fig.~\ref{fig:trion}, e.g., the
admixture of the $n=0$ state is relatively large, $\sim {\Delta}/U$,
while the $n=3$ state has a much smaller weight, $\sim
({\Delta}/U)^2$. As a consequence, most of the spectral weight appears
on the particle-like side of the spectral function, $\omega>0$. At
small frequencies the spectral function decays as $\sim
\omega^2$. This is in agreement with Fermi liquid theory, which would
predict a $\langle T(t) T^\dagger(0)\rangle \sim 1/t^3 $ decay of the
trionic correlation functions at very long times.

\section{Conclusions}
\label{sec:conclusions}

In this paper, we showed how to extend the DM-NRG scheme of
Ref.~\onlinecite{Toth} to symmetries with outer multiplicities. As an
application, we performed a detailed DM-NRG study of the \SU 3
symmetrical Anderson model by first incorporating \SU N
symmetries\cite{ArneJan} in the Open Access Budapest DM-NRG
code,\cite{OpenAccess} and then performing the numerical calculations
using the complete $\text U(1)\times \SU 3$ symmetry of the model.  A
similar extension has been carried out within the matrix product state
(MPS) approach parallel to this work.\cite{Andreas}

The properties of the \SU 3 Anderson model do not differ so much from those of the original Anderson model. As also discussed in  Refs.~\onlinecite{Carmi}, for $U>{\Delta}$
four distinct charging regions appear: the featureless empty and
fully occupied  regions ($\langle n\rangle \approx 3$ and $\langle n\rangle \approx 0$), and
two Kondo regions of occupancies $\langle n\rangle \approx 1$
and $\langle n\rangle \approx 2$, respectively.
The two \SU 3 Kondo regions  behave similarly: they are characterized by phase shifts
$\delta \sim \pm \pi/3$, as verified from the finite size spectrum, and correspondingly, a Kondo resonance shifted away from the Fermi energy.
In these Kondo regimes, the susceptibility has a universal
 form, $\chi(\omega)= f(\omega/T_K, \langle n\rangle) / T_K$, with a scaling function $f(x,\langle n\rangle)$ very similar to the one emerging in the \SU 2
Anderson and Kondo  models.
For completeness, we also studied the spectral properties of other local operators such as pairing or  the trion operators. Both of them turn out to have a  small amplitude for ${\Delta}\ll U$, and exhibit simple Fermi liquid properties below the mixed valence  and  Kondo scales, respectively.
Therefore, away from half filling,  a magnetic instability is expected to prevail on a lattice
in the \SU 3 Hubbard model, in general agreement with the results of Gutzwiller calculations
at low temperatures.\cite{RappRosch}

The \SU 3 Kondo regions are separated by a mixed valence state, which again has a Fermi liquid
character with a Fermi liquid scale  of the order of the level width, ${\Gamma}$. Here we find a phase shift
$\delta = \pi/2$, in agreement with the expectations based upon the Friedel sum rule, but apart from that, and the emerging electron-hole symmetry at this point, the properties of the mixed valence state appear to be very quite similar to those of the Kondo states.

We also investigated the conductance properties of the \SU 3
arrangement, proposed in Ref.~\onlinecite{Carmi}. At high temperatures
we observe in the side-conductance three charging peaks, corresponding
to the three charging steps. As the temperature is lowered, the
Coulomb blockade valleys are gradually filled up, and two conductance
shoulders emerge in the Kondo regime with a conductance $G\approx
(3/4)(e^2/h)$.

The methods and the computations presented here represent a first and
important step to perform DM-NRG and DMFT calculations for more
elaborate \SU N or $\text{Sp}(N)$ lattice models.  However, further work is
necessary to optimize
these DMFT calculations. While we definitely gain
enormous storage space by using \SU N symmetries, the evaluation of
the reduced matrix elements and multiple sums over internal
representation labels are currently not carried out with maximal efficiency.
 Since for large $N$'s the dimensions of irreducible \SU N
representations grow very fast, these summations quickly become the
bottleneck for DM-NRG calculations, and further work is needed to
increase the efficiency of the corresponding subroutines.

\emph{Acknowledgment:}
We acknowledge numerous fruitful discussions with A. Weichselbaum.
We thank the National Institute for Theoretical Physics in Stellenbosch,
South Africa, where this collaboration was started, for hospitality,
and the DFG (De730/7-1) for financially supporting our visit there.
This research has been supported by Hungarian research funds
OTKA and NKTH under Grant Nos.~K73361 and CNK80991,
and the EU-NKTH GEOMDISS project.
JvD and AA acknowledge financial support from SFB-TR12.
CPM has benefit from financial support from
UEFISCDI under French-Romanian Grant
DYMESYS (ANR 2011-IS04-001-01 and Contract
No. PN-II-ID-JRP-2011-1).