The Hierarchical Graphene model

Ian Jauslin

Department of Mathematics, Rutgers University

<EMAIL_ADDRESS>

These are the lecture notes for the summer school “Quantum Mechanics from
Condensed Matter to Computing”, organized by Niels Benedikter, Marcin
Napiórkowski, Jan Philip Solovej and Albert Werner in Copenhagen from June 13
to 17, 2022.

Abstract

The hierarchical graphene model is a simple toy model which is useful to
understand the mechanics of renormalization group flows in super-
renormalizable systems. It is based on a model of interacting electrons in
graphene, for which the renormalization group analysis was carried out by
Giuliani and Mastropietro. The analysis of the hierarchical graphene model is
significantly simpler than graphene, but one should not expect it to produce
good quantitative results about real-world graphene. Rather, the hierarchical
model is useful as a teaching tool to understand the core concepts of
renormalization group techniques. In this paper, we will first introduce a
model for electrons in graphene and set it up for a renormalization group
treatment by introducing its Grassmann representation and scale decomposition.
We then define the hierarchical graphene model and study it’s renormalization
group flow. From a renormalization group point of view, graphene is quite
simple: it is super-renormalizable. As an illustration of a more complicated
system, we repeat the analysis for the Kondo model, which is a strongly
coupled model with a non-trivial fixed point.

###### Contents

  1. 1 Introduction
  2. 2 Graphene
    1. 2.1 Lattice
    2. 2.2 Hamiltonian
    3. 2.3 Grassmann integral representation
    4. 2.4 Singularities of the propagator and scale decomposition
  3. 3 Hierarchical graphene
    1. 3.1 Definition of the hierarchical model
    2. 3.2 Beta function of the hierarchical model
    3. 3.3 Results for non-hierarchical graphene
  4. 4 Hierarchical Kondo model
    1. 4.1 Hamiltonian
    2. 4.2 Scale decomposition and hierarchical model
    3. 4.3 Beta function
  5. A Free Fermions
    1. A.1 Two-point correlation function
    2. A.2 Wick rule
  6. B Properties of Gaussian Grassman integrals
    1. B.1 Gaussian Grassmann integrals
    2. B.2 Wick rule for Gaussian Grassmann integrals
    3. B.3 Addition property

##  1 Introduction

The renormalization group is a powerful technique used to study a wide variety
of systems: from field theories to statistical mechanical systems. The
technique comes in many different flavors, some more appropriate to numerical
computations, others tailored to theoretical analyses. Some can be justified
mathematically, whereas for others, understanding why they work is still an
open problem. In this paper, our focus will be purposefully narrow, and we
will focus on the Wilsonian renormalization group [Wi65], which, to put it
simply, consists in separately considering different energy scales, and
studying how systems behave in each one. We will only discuss Fermionic
statistical field theories, that is, systems of many Fermions, for which
mathematically complete analyses abound. Our point of view will be that
developed by Benfatto and Gallavotti [BG90], which has been applied
successfully to a variety of systems: to name but a few, [BGe94, Ma11, GMP12,
GGM12, GMP17] among many others.

Of particular interest in this paper will be a model for interacting electrons
in graphene. Graphene is a two-dimensional crystal of carbon atoms in a
honeycomb structure, whose discovery in 2007 set off a flurry of interest
[NGe04] due to its unusual and potentially useful properties. To start out
with, two dimensional crystals (without a substrate) are quite rare. This two-
dimensional structure makes graphene crystals atom-thin and very flexible, but
since the carbon atoms in graphene are bound covalently, it can sustain very
high stresses. It is also an extremely good conductor. These properties give
graphene great potential in many technological applications, from flexible
displays to lightweight, conducting and tear-resistant plastics.

Graphene is often studied in the approximation that its electrons do not
interact with each other [Mc57, SW58]. In this case, the electronic properties
of graphene can be computed exactly. However, taking into account interactions
between electrons is more complicated. We will use the renormalization group
to accomplish just this, and show that the interactions, if they are weak
enough, do not change the physical picture much. In other words, the
renormalization group allows us to set up a perturbation theory for the
electrons in graphene, where the perturbation is the interaction. Perturbation
theory in many-Fermion systems is relatively straightforward when the
unperturbed Hamiltonian is “gapped”. This is not the case with graphene, so
standard perturbation theory does not work. However, the gap closes only at
two points in momentum space, called Fermi points, and near these Fermi
points, the bands are approximately conical. Both of these properties make
graphene, from a renormalization group point of view, super-renormalizable,
which one can understand as “easier to study using the renormalization group
than generic models”. This is the prime motivation for studying this model in
this paper: it will serve as an example in which to understand the core
concepts of Wilsonian renormalization.

The renormalization group analysis of interacting graphene was carried out in
[GM10, Gi10], and was used to prove the universality of the conductivity
[GMP12]. As is apparent from the length of these papers, despite the relative
simplicity of graphene, it is still not a trivial task to study it using the
renormalization group. In order to keep the discussion as simple as possible
and nevertheless get to the core ideas of the renormalization group, we will
simplify the graphene model, and introduce the hierarchical graphene model.
Hierarchical models have long been used: [Dy69, BCe78, GK81] as toy models in
which to understand renormalization group analyses in a simpler setting. The
renormalization group is notorious for requiring a large amount of
mathematical details to be worked out. Hierarchical models can be studied
without so many difficulties, and can help to grasp the conceptual core of the
renormalization group. Hierarchical models for Fermionic systems, first
studied by Dorlas [Do91], are even simpler than Bosonic ones: as was shown in
[BGJ15, GJ15], they are integrable, in that their renormalization group
analysis can be carried out explicitly (in Bosonic and non-hierarchical cases,
the renormalization group produces infinite power series). It is important to
emphasize that hierarchical models are toy models: they are not approximations
of their non-hierarchical counterparts, nor do they make good quantitative
predictions about them. They are only really useful to understand
renormalization group flows on a conceptual level.

For this paper, we will define the hierarchical graphene model, following the
ideas of [BGJ15]. (This model has not been introduced before, though this is
more likely because it is more useful as a teaching tool than as a physical
model.) The main idea is to eliminate everything from the graphene model
except for its scaling properties. This will leave us with a simple model, for
which we can compute the renormalization group flow exactly and explicitly,
which will help us understand what a renormlization group flow is, and how it
can be used. We will briefly discuss how to adapt the analysis to non-
hierarchical graphene, but few details will be given. Interested readers are
referred to the detailed presentation in [Gi10].

The hierarchical graphene model will allow us to get a handle on the
renormalization group analysis for graphene, which, as was mentioned above, is
a perturbative analysis (by which we mean that the interaction is a
perturbation; the term “perturbative” is sometimes used to mean “formal”, but
this is not what is meant here). Conceptually, the renormalization group has
much to say about non-perturbative systems. For instance, consider a system of
electrons in a superconducting phase. On short length scales (high energies),
the electrons are essentially independent from each other, but on large length
scales (low energies), they form Cooper pairs, which allows the system to
conduct electricity without resistance. From a renormalization group point of
view, we should see that at small distances, interactions between electrons
are not important, but at large distances, they change the behavior of
electrons qualitatively. In other words, the effective model at large length
scales should be very different from the non-interacting one (in fact, one
should see BCS theory emerging at these scales). However, the renormalization
group is so difficult to study away from the perturbative regime that this
has, so far, never been accomplished (at least not mathematically): not for
BCS theory, nor for any other strongly coupled system.

However, Fermionic hierarchical models can be studied exactly using the
renormalization group, even in strongly coupled situations. We give an example
of such a system at the end of this paper: the hierarchical Kondo model. The
Kondo model is a one-dimensional system of electrons on a lattice that
interact with a localizaed magnetic impurity. It was introduced [Ko64] as a
toy model to study conductance in disordered systems. It was studied rather
extensively by Anderson [An70], and later found to be exactly solvable by
Andrei [An80]. In developing his version of the renormalization group, Wilson
studied the Kondo model [Wi75], though his analysis has not yet been made
mathematically rigorous. One way to see that the Kondo model is strongly
coupled is through the Kondo effect, in which the magnetic impurity can be
shown to have a finite susceptibility at zero temperature. This means that, in
the lowest energy state, if one applies a magnetic field to the impurity, it
will not align perfectly with the field. In the absence of interactions with
the electrons, this would obviously not be true (the susceptibility would be
infinite). But, even if the interaction is arbitrarily small, as long as it is
ferromagnetic, the susceptibility comes out finite.

The Kondo effect can be seen in the renormalization group flow, which, at
small energy scales, goes to a non-trivial fixed point, that is, to an
effective theory that is qualitatively different from the non-interacting one.
This is similar, in essence, to the BCS question described above. Whereas the
Kondo effect can be proved by using the exact solvability of the Kondo model
[An80], it has never been shown using the renormalization group. It is,
however, provable in the hierarchical Kondo model [BGJ15, GJ15], which will be
briefly discussed at the end of this paper.

The rest of this paper is structured as follows. In section 2, we introduce
the (non-hierarchical) graphene model and set up its renormalization group
analysis by expressing observables in terms of Grassmann variables and
decomposing the system into scales. In section 3, we define the hierarchical
graphene model and study its renormalization group flow. In section 4, we
define the hierarchical Kondo model and discuss its renormalization group
flow. There are two appendices in which useful lemmas are proved: in appendix
A we compute properties of free Fermion systems, and in appendix B, we prove
some properties of Gaussian Grassmann integrals.

##  2 Graphene

In this section, we describe a model for the electrons in graphene. We will
use the tight-binding approximation, in which the electrons are assumed to be
bound to carbon atoms on a hexagonal lattice (see figure 2). For more details
on this model, see [Gi10, GM10].

Figure 1:  The hexagonal lattice. Each dot represents a carbon atom, and the
lines connect neighboring atoms, and representing the possible hoppings paths
for the electrons.

###  2.1 Lattice

The hexagonal lattice can be constructed by copying an elementary cell at
every integer combination of

$l_{1}:=\left(\frac{3}{2},\frac{\sqrt{3}}{2}\right),\quad
l_{2}:=\left(\frac{3}{2},-\frac{\sqrt{3}}{2}\right)$ (1)

where we have chosen the unit length to be equal to the distance between two
nearest neighbors. The elementary cell consists of two atoms at $(0,0)$ and at
$(0,1)$ (relative to the position of the cell).

Figure 2:  The cell decomposition of the hexagonal lattice. Each cell (blue
rhombus) contains two atoms. One of type $a$ (circle) and one of type $b$
(square).

We define the lattice

$\Lambda:=\left\\{n_{1}l_{1}+n_{2}l_{2},\
(n_{1},n_{2})\in\\{0,\cdots,L-1\\}^{2}\right\\}$ (2)

where $L$ is a positive integer that determines the size of the crystal, that
we will eventually send to infinity, with periodic boundary conditions. We
introduce the nearest neighbor vectors:

$\delta_{1}:=(1,0),\quad\delta_{2}:=\left(-\frac{1}{2},\frac{\sqrt{3}}{2}\right),\quad\delta_{3}:=\left(-\frac{1}{2},-\frac{\sqrt{3}}{2}\right).$
(3)

The dual of $\Lambda$ is

$\hat{\Lambda}:=\left\\{\frac{m_{1}}{L}G_{1}+\frac{m_{2}}{L}G_{2},\
(m_{1},m_{2})\in\\{0,\cdots,L-1\\}^{2}\right\\}$ (4)

with periodic boundary conditions, where

$G_{1}=\left(\frac{2\pi}{3},\frac{2\pi}{\sqrt{3}}\right),\quad
G_{2}=\left(\frac{2\pi}{3},-\frac{2\pi}{\sqrt{3}}\right).$ (5)

It is defined in such a way that $\forall x\in\Lambda$, $\forall
k\in\hat{\Lambda}$,

$e^{ikxL}=1.$ (6)

In the limit $L\to\infty$, the set $\hat{\Lambda}$ tends to the torus
$\hat{\Lambda}_{\infty}=\mathbb{R}^{2}/(\mathbb{Z}G_{1}+\mathbb{Z}G_{2})$,
also called the Brillouin zone, see figure 3.

Figure 3:  The Brillouin zone. The two dots are the Fermi points
$p_{F}^{(\pm)}$, see below.

###  2.2 Hamiltonian

We will consider a model of spin-$\frac{1}{2}$ electrons on the hexagonal
lattice. Given $x\in\Lambda$, we denote the Fermionic annihilation operators
with spin $\sigma$ at atoms of type $a$ and $b$ within the elementary cell
centered at $x$ respectively by $a_{x,\sigma}$ and $b_{x+\delta_{1},\sigma}$.

The Hamiltonian is split into two terms:

$\mathcal{H}=\mathcal{H}_{0}+\mathcal{H}_{I}$ (7)

where $\mathcal{H}_{0}$ is the free Hamiltonian, which describes the motion of
electrons from one atom to a neighbor, and $\mathcal{H}_{I}$ is the
interaction Hamiltonian, which describes the interaction between electrons.

1 -  Free Hamiltonian. The free Hamiltonian describes the hopping of electrons
from one atom to another:

$\mathcal{H}_{0}:=-\sum_{\sigma\in\\{\uparrow,\downarrow\\}}\sum_{\displaystyle\mathop{\scriptstyle
x\in\Lambda}_{j=1,2,3}}(a_{x,\sigma}^{\dagger}b_{x+\delta_{j},\sigma}+b_{x+\delta_{j},\sigma}^{\dagger}a_{x,\sigma})$
(8)

Equation (8) can be rewritten in Fourier space as follows. We define the
Fourier transform of the annihilation operators as

$\hat{a}_{k,\sigma}:=\frac{1}{\sqrt{|\Lambda|}}\sum_{x\in\Lambda}e^{ikx}a_{x,\sigma},\quad\hat{b}_{k,\sigma}:=\frac{1}{\sqrt{|\Lambda|}}\sum_{x\in\Lambda}e^{ikx}b_{x+\delta_{1},\sigma}$
(9)

where $|\Lambda|=L^{2}$. Note that, with this choice of normalization,
$\hat{a}_{k,\sigma}$ and $\hat{b}_{k,\sigma}$ satisfy the canonical
anticommutation relations:

$\\{a_{k,\sigma},a_{k^{\prime},\sigma^{\prime}}^{\dagger}\\}=\frac{1}{|\Lambda|}\sum_{x,x^{\prime}\in\Lambda}e^{ikx-
ik^{\prime}x^{\prime}}\\{a_{x,\sigma},a_{x^{\prime},\sigma^{\prime}}^{\dagger}\\}=\delta_{\sigma,\sigma^{\prime}}\frac{1}{|\Lambda|}\sum_{x^{\prime}\in\Lambda}e^{i(k-k^{\prime})x}=\delta_{\sigma,\sigma^{\prime}}\delta_{k,k^{\prime}}$
(10)

and similarly for $b$. We express $\mathcal{H}_{0}$ in terms of $\hat{a}$ and
$\hat{b}$:

$\mathcal{H}_{0}=-\sum_{\sigma\in\\{\uparrow,\downarrow\\}}\sum_{k\in\hat{\Lambda}}\hat{A}_{k,\sigma}^{\dagger}H_{0}(k)\hat{A}_{k,\sigma}$
(11)

$\hat{A}_{k,\sigma}$ is a column vector whose transpose is
$\hat{A}_{k,\sigma}^{T}=(\hat{a}_{k,\sigma},\hat{b}_{k,\sigma})$,

$H_{0}(k):=\left(\begin{array}[]{*{2}{c}}0&\Omega^{*}(k)\\\
\Omega(k)&0\end{array}\right)$ (12)

and

$\Omega(k):=\sum_{j=1}^{3}e^{ik(\delta_{j}-\delta_{1})}=1+2e^{-i\frac{3}{2}k_{x}}\cos({\textstyle\frac{\sqrt{3}}{2}k_{y}}).$
(13)

The eigenvalues of $H_{0}(k)$ are called the bands of non-interacting
graphene, and are

$\pm|\Omega(k)|=\pm\left(1+4\cos({\textstyle\frac{3}{2}k_{x}})\cos({\textstyle\frac{\sqrt{3}}{2}k_{y}})+4\cos^{2}({\textstyle\frac{\sqrt{3}}{2}k_{y}})\right)^{\frac{1}{2}}.$
(14)

These bands meet at $0$ at exactly two values of $k$: for $\omega=\pm$,

$p_{F}^{(\omega)}:=({\textstyle\frac{2\pi}{3},\omega\frac{2\pi}{3\sqrt{3}}})$
(15)

see figure 4. For $|k-p_{F}^{(\omega)}|\ll 1$,

$\pm|\Omega(k)|\sim\pm v_{F}|k-p_{F}^{(\omega)}|,\quad v_{F}=\frac{3}{2}.$
(16)

Figure 4:  The bands of graphene. The two bands meet at two points in a
conical intersection.

2 -  Interaction. We now define the interaction Hamiltonian which we take to
be of Hubbard form:

$\mathcal{H}_{I}:=U\sum_{x\in\Lambda}\sum_{\alpha\in\\{a,b\\}}\left(\alpha_{x+d_{\alpha},\uparrow}^{\dagger}\alpha_{x+d_{\alpha},\uparrow}-\frac{1}{2}\right)\left(\alpha_{x+d_{\alpha},\downarrow}^{\dagger}\alpha_{x+d_{\alpha},\downarrow}-\frac{1}{2}\right)$
(17)

where the $d_{\alpha}$ are the vectors that give the position of each atom
type with respect to the centers of the lattice $\Lambda$: $d_{a}:=0$,
$d_{b}:=\delta_{1}$.

###  2.3 Grassmann integral representation

The specific free energy on the lattice $\Lambda$ is defined by

$f_{\Lambda}:=-\frac{1}{\beta|\Lambda|}\log\left(\mathrm{Tr}\left(e^{-\beta\mathcal{H}}\right)\right)$
(18)

where $\beta$ is the inverse temperature. We define these at finite $\beta$
and $L$, but will take $\beta,L\to\infty$. A straightforward application of
the Trotter product formula implies that (see [Gi10, (4.1)])

$\mathrm{Tr}(e^{-\beta\mathcal{H}})=\mathrm{Tr}(e^{-\beta\mathcal{H}_{0}})+\sum_{N=1}^{\infty}\frac{(-\beta)^{N}}{N!}\int_{\beta\geqslant
t_{1}\geqslant\cdots\geqslant t_{N}\geqslant
0}\mathrm{Tr}\left(e^{-\beta\mathcal{H}_{0}}\mathcal{H}_{I}(t_{1})\cdots\mathcal{H}_{I}(t_{N})\right)$
(19)

where

$\mathcal{H}_{I}(t):=e^{t\mathcal{H}_{0}}\mathcal{H}_{I}e^{-t\mathcal{H}_{0}}.$
(20)

To compute this trace, we will use the Wick rule, which we will now descibe.
First, we define the free average:

$\left<A\right>:=\frac{\mathrm{Tr}(e^{-\beta\mathcal{H}_{0}}A)}{\mathrm{Tr}(e^{-\beta\mathcal{H}_{0}})}.$
(21)

Next, we define the imaginary time creation and annihilation operators: for
$\alpha\in\\{a,b\\}$ and $t\in[0,\beta)$,

$\alpha_{x,\sigma}^{+}(t):=e^{t\mathcal{H}_{0}}\alpha_{x,\sigma}^{\dagger}e^{-t\mathcal{H}_{0}},\quad\alpha_{x,\sigma}^{-}(t):=e^{t\mathcal{H}_{0}}\alpha_{x,\sigma}e^{-t\mathcal{H}_{0}}$
(22)

(note that $\alpha^{+}$ is not the adjoint of $\alpha^{-}$). The Wick rule can
be used to compute the free average of any polynomial of the creation and
annihilation operators: it is linear, and, for any $n\in\mathbb{N}$,
$\alpha^{(1)},\bar{\alpha}^{(1)},\cdots,\alpha^{(n)},\bar{\alpha}^{(n)}\in\\{a,b\\}$
$x_{1},\bar{x}_{1},\cdots,x_{n},\bar{x}_{n}\in\Lambda$,
$\sigma_{1},\bar{\sigma}_{1},\cdots,\sigma_{n},\bar{\sigma}_{n}\in\\{\uparrow,\downarrow\\}$,
$\beta\geqslant
t_{1}\geqslant\bar{t}_{1}>\cdots>t_{n}\geqslant\bar{t}_{n}\geqslant 0$,

$\left<\prod_{i=1}^{n}\alpha^{(i)-}_{x_{i},\sigma_{i}}(t_{i})\bar{\alpha}^{(i)+}_{\bar{x}_{i},\bar{\sigma}_{i}}(\bar{t}_{i})\right>=\sum_{\tau\in\mathcal{S}_{n}}(-1)^{\tau}\prod_{i=1}^{n}\left<\mathbf{T}\left(\alpha^{(i)-}_{x_{i},\sigma_{i}}(t_{i})\bar{\alpha}^{(\tau(i))+}_{\bar{x}_{\tau(i)},\bar{\sigma}_{\tau(i)}}(\bar{t}_{\tau(i)})\right)\right>$
(23)

(to alleviate the notation, we have replaced $\alpha_{x+d_{\alpha}}$ by
$\alpha_{x}$ here) where $\mathcal{S}_{n}$ is the set of permutations of
$\\{1,\cdots,n\\}$, $(-1)^{\tau}$ is the signature of $\tau$ and $\mathbf{T}$
is the time-ordering operator:

$\mathbf{T}\left(\alpha^{-}_{x,\sigma}(t)\bar{\alpha}^{+}_{\bar{x},\bar{\sigma}}(\bar{t})\right)=\left\\{\begin{array}[]{>\displaystyle
ll}\alpha^{-}_{x,\sigma}(t)\bar{\alpha}^{+}_{\bar{x},\bar{\sigma}}(\bar{t})&\mathrm{if\
}t\geqslant\bar{t}\\\
-\bar{\alpha}^{+}_{\bar{x},\bar{\sigma}}(\bar{t})\alpha^{-}_{x,\sigma}(t)&\mathrm{if\
}t<\bar{t}.\end{array}\right.$ (24)

The Wick rule can be proved by a direct computation, and follows from the fact
that $\mathcal{H}_{0}$ is quadratic in the annihilation operators, see lemma
A.2.

Fermionic creation and annihilation operators do not anticommute:
$\\{a_{i},a_{i}\dagger\\}=1$. However, the time-ordering operator effectively
makes them anticommute. We can make this more precise be re-expressing the
problem in terms of Grassmann variables.

3 -  Definition of the Grassmann algebra. We first define a Grassmann algebra
and an integration procedure on it. We will work in Fourier space, in both
space and time. Let us first define the Fourier transform in time: for every
$\alpha\in\\{a,b\\}$, $\sigma\in\\{\uparrow,\downarrow\\}$ for $k_{0}\in
2\pi\beta^{-1}(\mathbb{Z}+1/2)$, $k\in\hat{\Lambda}$, and
$\mathbf{k}\equiv(k_{0},k)$

$\hat{\alpha}_{\mathbf{k},\sigma}^{\pm}:=\int_{0}^{\beta}dt\ e^{\mp
itk_{0}}e^{\mathcal{H}_{0}t}\hat{\alpha}_{k,\sigma}^{\pm}e^{-\mathcal{H}_{0}t}\equiv\frac{1}{\sqrt{|\Lambda|}}\sum_{x\in\Lambda}\int_{0}^{\beta}dt\
e^{\mp(itk_{0}+ikx)}\alpha_{x,\sigma}^{\pm}(t)$ (25)

in which we use the shorthand
$\hat{\alpha}_{k,\sigma}^{-}\equiv\hat{\alpha}_{k,\sigma}$,
$\hat{\alpha}_{k,\sigma}^{+}\equiv\hat{\alpha}_{k,\sigma}^{\dagger}$. We
notice that
$\mathbf{k}\in\mathcal{B}_{\beta,L}:=(2\pi\beta^{-1}(\mathbb{Z}+1/2))\times\hat{\Lambda}$
varies in an infinite set. Since this will cause trouble when defining
Grassmann integrals, we shall impose a cutoff $M\in\mathbb{N}$: let
$\chi_{0}(\rho)$ be a smooth compact support function that returns $1$ if
$\rho\leqslant 1/3$ and $0$ if $\rho\geqslant 2/3$, and let

$\mathcal{B}_{\beta,L}^{*}:=\mathcal{B}_{\beta,L}\cap\\{(k_{0},k),\
\chi_{0}(2^{-M}|k_{0}|\neq 0)\\}.$ (26)

To every $\hat{\alpha}_{\mathbf{k},\sigma}^{\pm}$ for $\alpha\in\\{a,b\\}$ and
$\mathbf{k}\in\mathcal{B}_{\beta,L}^{*}$, we associate a Grassmann variable
$\hat{\psi}_{\mathbf{k},\alpha,\sigma}^{\pm}$, and we consider the finite
Grassmann algebra (an algebra in which the $\hat{\psi}$ anti-commute with each
other) generated by the collection
$\\{\hat{\psi}_{\mathbf{k},\alpha,\sigma}^{\pm}\\}_{\mathbf{k}\in\mathcal{B}_{\beta,L}^{*}}^{\alpha\in\\{a,b\\},\sigma\in\\{\uparrow,\downarrow\\}}$.
We define the Grassmann integral

$\int\prod_{\sigma\in\\{\uparrow,\downarrow\\}}\prod_{\alpha\in\\{a,b\\}}\prod_{\mathbf{k}\in{\mathcal{B}}^{*}_{\beta,L}}d\hat{\psi}_{\mathbf{k},\alpha,\sigma}^{+}d\hat{\psi}_{\mathbf{k},\alpha,\sigma}^{-}$
(27)

as the linear operator on the Grassmann algebra whose action on a monomial in
the variables $\hat{\psi}^{\pm}_{\mathbf{k},\alpha,\sigma}$ is $0$ except if
said monomial is
$\prod_{\sigma\in\\{\uparrow,\downarrow\\}}\prod_{\alpha\in\\{a,b\\}}\prod_{\mathbf{k}\in{\mathcal{B}}_{\beta,L}^{*}}\hat{\psi}^{-}_{\mathbf{k},\alpha,\sigma}\hat{\psi}^{+}_{\mathbf{k},\alpha,\sigma}$
up to a permutation of the variables, in which case the value of the integral
is determined using

$\int\prod_{\sigma\in\\{\uparrow,\downarrow\\}}\prod_{\alpha\in\\{a,b\\}}\prod_{\mathbf{k}\in{\mathcal{B}}_{\beta,L}^{*}}d\hat{\psi}_{\mathbf{k},\alpha,\sigma}^{+}d\hat{\psi}_{\mathbf{k},\alpha,\sigma}^{-}\left(\prod_{\sigma\in\\{\uparrow,\downarrow\\}}\prod_{\alpha\in\\{a,b\\}}\prod_{\mathbf{k}\in{\mathcal{B}}_{\beta,L}^{*}}\hat{\psi}^{-}_{\mathbf{k},\alpha,\sigma}\hat{\psi}^{+}_{\mathbf{k},\alpha,\sigma}\right)=1$
(28)

along with the anti-commutation of the $\hat{\psi}$.

Let us now define Gaussian Grassmann integrals. The Gaussian Grassmann measure
is specified by a propagator, which is a $2\times 2$ complex matrix
$\hat{g}(\mathbf{k})$:

$\begin{array}[]{@{}>{\displaystyle}l@{}}\hphantom{\hskip
423.62pt}\\\\[-14.22636pt]
P_{\hat{g}}(d\psi):=\left(\prod_{\mathbf{k}\in\mathcal{B}_{\beta,L}^{*}}(\beta^{2}\det\hat{g}(\mathbf{k}))^{2}\left(\prod_{\sigma\in\\{\uparrow,\downarrow\\}}\prod_{\alpha\in\\{a,b\\}}d\hat{\psi}_{\mathbf{k},\alpha}^{+}d\hat{\psi}_{\mathbf{k},\alpha}^{-}\right)\right)\cdot\\\\[14.22636pt]
\hfill\cdot\exp\left(-\frac{1}{\beta}\sum_{\sigma\in\\{\uparrow,\downarrow\\}}\sum_{\mathbf{k}\in\mathcal{B}_{\beta,L}^{*}}\hat{\psi}^{+}_{\mathbf{k},\cdot,\sigma}\cdot\hat{g}^{-1}(\mathbf{k})\hat{\psi}^{-}_{\mathbf{k},\cdot,\sigma}\right).\end{array}$
(29)

By a direct computation, one can prove that (see lemma B.1)

$\int P_{\hat{g}}(d\psi)\ 1=1,\quad\int P_{\hat{g}}(d\psi)\
\hat{\psi}_{\mathbf{k},\alpha,\sigma}^{-}\hat{\psi}_{\bar{\mathbf{k}},\bar{\alpha},\bar{\sigma}}^{+}=\beta\delta_{\mathbf{k},\bar{\mathbf{k}}}\delta_{\sigma,\bar{\sigma}}\hat{g}_{\alpha,\bar{\alpha}}(\mathbf{k}).$
(30)

In addition, the integral with respect to $P_{\hat{g}}(d\psi)$ satisfies the
Wick rule (see lemma B.2):

$\int
P_{\hat{g}}(d\psi)\prod_{i=1}^{n}\hat{\psi}_{\mathbf{k}_{i},\alpha_{i},\sigma_{i}}^{-}\hat{\psi}_{\bar{\mathbf{k}}_{i},\bar{\alpha}_{i},\bar{\sigma}_{i}}^{+}=\sum_{\tau\in\mathcal{S}_{n}}(-1)^{\tau}\prod_{i=1}^{n}\int
P_{\hat{g}}(d\psi)\
\hat{\psi}_{\mathbf{k}_{i},\alpha_{i},\sigma_{i}}^{-}\hat{\psi}_{\bar{\mathbf{k}}_{\tau(i)},\bar{\alpha}_{\tau(i)},\bar{\sigma}_{\tau(i)}}^{+}$
(31)

Finally, given two propagators $\hat{g}_{1}$ and $\hat{g}_{2}$, and any
polynomial $\mathfrak{P}(\psi)$ in the Grassmann variables, we have (see lemma
B.3),

$\int P_{\hat{g}_{1}+\hat{g}_{2}}(d\psi)\ \mathfrak{P}(\psi)=\int
P_{\hat{g}_{1}}(d\psi_{1})\int P_{\hat{g}_{2}}(d\psi_{2})\
\mathfrak{P}(\psi_{1}+\psi_{2}).$ (32)

4 -  Grassmann integrals and the free energy. Let us now make the connection
between the computation of the free energy and Grassmann integrals. As we have
seen, free averages of polynomials in the creation and annihilation operators
can be computed using the Wick rule, see (23). Gaussian Grassmann integrals of
Grassmann polynomials also satisfy the Wick rule, see (31), so Gaussian
Grassmann integrals and the free average share the same algebraic structure.
Thus, if we set $\hat{g}$ in such a way that

$\int P_{\hat{g}}(d\psi)\
\hat{\psi}_{\mathbf{k},\alpha,\sigma}^{-}\hat{\psi}_{\bar{\mathbf{k}},\bar{\alpha},\bar{\sigma}}^{+}=\frac{1}{|\Lambda|}\sum_{x,\bar{x}\in\Lambda}\int_{0}^{\beta}dt\int_{0}^{\beta}d\bar{t}\
e^{itk_{0}-i\bar{t}\bar{k}_{0}+ikx-i\bar{k}\bar{x}}\left<\mathbf{T}\left(\alpha_{x,\sigma}^{-}(t)\bar{\alpha}^{+}_{\bar{x},\bar{\sigma}}(\bar{t})\right)\right>$
(33)

then we can compute free averages using Gaussian Grassmann integrals.
Furthermore, by a direct computation (see lemma A.1),

$\frac{1}{|\Lambda|}\sum_{x,\bar{x}\in\Lambda}\int_{0}^{\beta}dt\int_{0}^{\beta}d\bar{t}\
e^{itk_{0}-i\bar{t}\bar{k}_{0}+ikx-i\bar{k}\bar{x}}\left<\mathbf{T}\left(\alpha_{x,\sigma}^{-}(t)\bar{\alpha}^{+}_{\bar{x},\bar{\sigma}}(\bar{t})\right)\right>=\beta\delta_{\mathbf{k},\bar{\mathbf{k}}}\delta_{\sigma,\bar{\sigma}}(-ik_{0}\mathds{1}-H_{0}(k))^{-1}$
(34)

so, by (30),

$\hat{g}(\mathbf{k})=(-ik_{0}\mathds{1}-H_{0}(k))^{-1}.$ (35)

Actually, since we cut off the momenta by $M$, in order to avoid introducing
Gibbs phenomena when inverting Fourier transforms, we define the propagator:

$\hat{g}_{\leqslant
M}(\mathbf{k}):=\chi_{0}(2^{-M}|k_{0}|)(-ik_{0}\mathds{1}-H_{0}(k))^{-1}$ (36)

and will take the limit $M\to\infty$ at the end of the computation. Thus, we
define the Gaussian Grassmann integration measure $P_{\leqslant
M}(d\psi)\equiv P_{\hat{g}_{\leqslant M}}(d\psi)$, which allows us to compute
the trace in (18):

$f_{\Lambda}=f_{0,\Lambda}-\lim_{M\to\infty}\frac{1}{\beta|\Lambda|}\log\int
P_{\leqslant M}(d\psi)\ e^{-\mathcal{V}(\psi)}$ (37)

where $f_{0,\Lambda}$ is the free energy in the $U=0$ case and

$\mathcal{V}(\psi)=U\sum_{\alpha\in\\{a,b\\}}\int_{0}^{\beta}dt\sum_{x\in\Lambda}\psi^{+}_{\mathbf{x},\alpha,\uparrow}\psi^{-}_{\mathbf{x},\alpha,\uparrow}\psi^{+}_{\mathbf{x},\alpha,\downarrow}\psi^{-}_{\mathbf{x},\alpha,\downarrow}$
(38)

in which $\mathbf{x}\equiv(t,x)$ and

$\psi^{\pm}_{\mathbf{x},\alpha,\sigma}:=\frac{1}{\beta\sqrt{|\Lambda|}}\sum_{\mathbf{k}\in{\mathcal{B}}^{*}_{\beta,L}}\hat{\psi}^{\pm}_{\mathbf{k},\alpha,\sigma}e^{\pm
i\mathbf{k}\mathbf{x}}.$ (39)

Note that we have dropped the $-\frac{1}{2}$ from the interaction Hamiltonian
when passing to the Grassmann variables. The reason for this is that, in
configuration space, the point $t=\bar{t}$ is important (this situation did
not arise above as this was of measure $0$). Recall that in lemma A.1, the two
point function $\left<\mathbf{T}(\alpha^{-}\alpha^{+})\right>$ at $t=\bar{t}$
has an extra $\frac{1}{2}$, and one can show that this exactly cancels with
the $-\frac{1}{2}$ in the interaction Hamiltonian.

###  2.4 Singularities of the propagator and scale decomposition

We thus have a clear strategy to compute the free energy of the electrons in
graphene: (37) reduces the computation to a Gaussian Grassmann integral, which
we can compute as a power series by expanding $e^{-\mathcal{V}}$. However,
when implementing this strategy, one almost immediately runs into a problem:
$\hat{g}$ is singular, at least in the limit $\beta,L\to\infty$. Indeed,
recall that the eigenvalues of $H_{0}$ are singular at $k=p_{F}^{(\omega)}$,
see (15), around which they behave like $v_{F}|k-p_{F}^{(\omega)}|$, see (16).
Thus, $\hat{g}_{\leqslant M}$ is singular at
$\mathbf{k}=\mathbf{p}_{F}^{(\omega)}:=(0,p_{F}^{(\omega)})$, near which the
eigenvalues of $\hat{g}_{\leqslant M}$ behave like

$\pm|v_{F}|k-p_{F}^{(\omega)}|+ik_{0}|^{-1}.$ (40)

(Note that $|k_{0}|\geqslant\frac{\pi}{\beta}$, so $\hat{g}$ is not singular
for finite $\beta$, but it becomes so in the limit.) Therefore, it is far from
obvious that our strategy to compute the free energy will work uniformly in
$\beta$, and so we may not be able to take the limit $\beta\to\infty$ (the
risk is that the series expansion does not converge in the limit).

We will need to proceed in a more subtle fashion, by performing a multiscale
decomposition (the multiscale decomposition is the heart of the
renormalization group method). The idea is to approach the singularities
$p_{F}^{(\omega)}$ slowly, by defining scale-by-scale propagators: for
$h\in-\mathbb{N}$, we define

$\Phi_{h}(\mathbf{k}-\mathbf{p}_{F}^{(\omega)}):=(\chi_{0}(2^{-h}|\mathbf{k}-\mathbf{p}_{F}^{(\omega)}|)-\chi_{0}(2^{-h+1}|\mathbf{k}-\mathbf{p}_{F}^{(\omega)}|))$
(41)

which is a smooth function that is supported in
$|\mathbf{k}-\mathbf{p}_{F}^{(\omega)}|\in[2^{h}\frac{1}{6},2^{h}\frac{2}{3}]$,
in other words, it localizes $\mathbf{k}$ to be at a distance from
$\mathbf{p}_{F}^{(\omega)}$ that is of order $2^{h}$, see figure 5. Since
$|k_{0}|\geqslant\frac{\pi}{\beta}$, we only need to consider

$h\geqslant-N_{\beta}:=\log_{2}\frac{\pi}{\beta}.$ (42)

We then define

$\hat{g}^{(h,\omega)}(\mathbf{k}):=\Phi_{h}(\mathbf{k}-\mathbf{p}_{F}^{(\omega)})(-ik_{0}\mathds{1}-H_{0}(k))^{-1}=\Phi_{h}(\mathbf{k}-\mathbf{p}_{F}^{(\omega)})\frac{1}{k_{0}^{2}+|\Omega(k)|^{2}}\left(\begin{array}[]{cc}ik_{0}&-\Omega^{*}(k)\\\
-\Omega(k)&ik_{0}\end{array}\right)$ (43)

(where we used (12) to compute the inverse). For $h$ sufficiently small
(sufficiently negative), we thus have

$\hat{g}^{(h,\omega)}=O(2^{-h}).$ (44)

To carry out the Gaussian Grassmann integral, we split

$\hat{g}_{\leqslant
M}=\sum_{h=0}^{-N_{\beta}}\sum_{\omega=\pm}\hat{g}^{(h,\omega)}(\mathbf{k})+\hat{g}_{\geqslant
0},\quad\hat{g}_{\geqslant 0}:=\hat{g}_{\leqslant
M}-\sum_{h=0}^{-N_{\beta}}\sum_{\omega=\pm}\hat{g}^{(h,\omega)}(\mathbf{k})$
(45)

(note that $\hat{g}_{\geqslant 0}$ is bounded uniformly in $\beta$) and use
(32) to compute the integrals with the propagators $g^{(h,\omega)}$ one at a
time.

Figure 5:  A sketch of the scale decomposition, the darker the color, the
smaller the scale.

This has been done in detail in [GM10, Gi10]. Carrying out this strategy in
practice is rather involved though, so here, we will simplify the problem by
considering a hierarchical model based on the graphene model.

##  3 Hierarchical graphene

The hierarchical graphene model is a simplification of the graphene model, in
which everything but the multiscale structure is culled. The philosophy of
using a hierarchical model is that renormalization group treatments are
typically quite involved, whereas hierarchical models behave in essentially
(from a renormalization group point of view) the same way, without having too
many details muddy the main ideas of the strategy. They have been used in many
settings [Dy69, BCe78, GK81], though most of those are Bosonic. Fermionic
hierarchical models were first studied in [Do91], and, in a more systematic
way, in [BGJ15, GJ15].

###  3.1 Definition of the hierarchical model

5 -  Boxes. The hierarchical model will be defined in configuration space (as
opposed to Fourier space). In the previous section, we introduced a multiscale
decomposition in Fourier space, in which we split Grassmann fields into
momentum shells, where $|\mathbf{k}-\mathbf{p}_{F}^{(\omega)}|\sim 2^{h}$. In
configuration space, this corresponds to considering Grassmann fields that are
constant on configuration-space boxes of size $\sim 2^{-h}$. These boxes are
obtained by doubling the size of the elementary cell of the hexagonal lattice
at every scale, see figure 6. We also have a time dimension, which we split
into boxes of size $2^{|h|}$. We thus define the set of boxes on scale
$h\in\\{-N_{\beta},\cdots,0\\}$ by

$\mathcal{Q}_{h}:=\left\\{[i2^{|h|},(i+1)2^{|h|})\times(\Lambda\cap\\{2^{|h|}(n_{1}+x_{1})l_{1}+2^{|h|}(n_{2}+x_{2})l_{2},\
x_{1},x_{2}\in[0,1)\\})\right\\}_{i,n_{1},n_{2}\in\mathbb{Z}}$ (46)

($\mathcal{Q}_{h}$ is a set of sets). For every
$(t,x)\in[0,\beta)\times\Lambda$ and $h\in\\{-N_{\beta},\cdots,0\\}$, there
exists a unique box $\Delta^{(h)}(t,x)\in\mathcal{Q}_{m}$ such that
$(t,x)\in\Delta^{(m)}(t,x)$. To simplify the computations further, we drop the
index $\omega$, which does not change the nature of the problem. In each box
$\Delta$, we have four Grassmann fields (and their conjugates), corresponding
to both choices of atom type ($a$ or $b$) and both choices of spin ($\uparrow$
and $\downarrow$): for $\alpha\in\\{a,b\\}$ and
$\sigma\in\\{\uparrow,\downarrow\\}$, $\psi_{\alpha,\sigma}^{[h]\pm}(\Delta)$.

Figure 6:  In this figure, we have quadrupled the size of the elementary cell.
Each large cell supports four Grassmann fields corresponding to both choices
of type of atom ($a$ and $b$) and of spin ($\uparrow$ and $\downarrow$).

6 -  Propagators. Let us now define the propagators associated to these
Grassmann fields. We will chose them to be similar, in a certain sense, to the
non-hierarchical propagators. We defined the single scale propagator in
Fourier space in (43), which, in configuration space, becomes

$g^{(h,\omega)}(t,x)=\frac{1}{\beta|\Lambda|}\sum_{\mathbf{k}\in\mathcal{B}_{\beta,L}}e^{-ik_{0}t-ikx}\Phi_{h}(\mathbf{k}-\mathbf{p}_{F}^{(\omega)})\frac{1}{k_{0}^{2}+|\Omega(k)|^{2}}\left(\begin{array}[]{cc}ik_{0}&-\Omega^{*}(k)\\\
-\Omega(k)&ik_{0}\end{array}\right).$ (47)

In the hierarchical approximation, we neglect propagators at different times,
so we can set $t=0$, at which, by symmetry, the diagonal terms vanish. Thus we
only have a non-trivial propagator in between atoms of type $a$ and $b$.
Furthermore,

$\frac{1}{\beta|\Lambda|}\sum_{\mathbf{k}\in\mathcal{B}_{\beta,L}}\Phi_{h}(\mathbf{k}-\mathbf{p}_{F}^{(\omega)})=O(2^{3h}),\quad\hat{g}^{(h,\omega)}=O(2^{-h})$
(48)

so

$g^{(h,\omega)}(t,x)=O(2^{2h}).$ (49)

In order to take this scaling factor into account, we rescale the Grassmann
fields $\psi$, and define, for $(t,x)\in[0,\beta)\times\Lambda$, (recall that
we have dropped the index $\omega$ in the hierarchical model)

$\psi_{\alpha,\sigma}^{\pm}(t,x)=\sum_{h=-N_{\beta}}^{0}2^{h}\psi_{\alpha,\sigma}^{[h]\pm}(\Delta^{(h)}(t,x)).$
(50)

We will take the propagators to be

$\int P^{[h]}(d\psi^{[h]})\
\psi_{a,\sigma}^{[h]-}(\Delta)\psi_{b,\sigma^{\prime}}^{[h]+}(\Delta^{\prime})=\delta_{\sigma,\sigma^{\prime}}\delta_{\Delta,\Delta^{\prime}}$
(51) $\int P^{[h]}(d\psi^{[h]})\
\psi_{b,\sigma}^{[h]-}(\Delta)\psi_{a,\sigma^{\prime}}^{[h]+}(\Delta^{\prime})=\delta_{\sigma,\sigma^{\prime}}\delta_{\Delta,\Delta^{\prime}}$
(52)

and all other propagators will be set to 0. We can now evaluate how well these
propagators approximate the non-hierarchical ones. Given $x,y$, let
$\eta(x,y)$ be the largest negative integer such that $x$ and $y$ are in the
same box. We have

$\int P(d\psi)\
\psi_{a,\sigma}^{-}(x,t)\psi_{b,\sigma^{\prime}}^{+}(y,t)=\delta_{\sigma,\sigma^{\prime}}\sum_{h=-N_{\beta}}^{\eta(x,y)}2^{2h}=\delta_{\sigma,\sigma^{\prime}}\frac{4}{3}2^{2\eta(x,y)}(1-4^{-N_{\beta}-\eta(x,y)})$
(53)

which has the same scaling as the non-hierarchical propagator (49) (and
similarly with $a$ and $b$ exchanged).

7 -  Effective potentials. In this hierarchical model, we compute

$\int P(d\psi)\ e^{-\mathcal{V}(\psi)}.$ (54)

We proceed inductively: for $\alpha\in\\{a,b\\}$,
$\sigma\in\\{\uparrow,\downarrow\\}$, $h\in\\{-N_{\beta},\cdots,0\\}$,
$\Delta\in\mathcal{Q}_{h}$, let $\bar{\Delta}\in\mathcal{Q}_{h-1}$ be such
that $\bar{\Delta}\supset\Delta$, we define

$\psi_{\alpha,\sigma}^{[\leqslant
h]\pm}(\Delta):=\frac{1}{2}\psi_{\alpha,\sigma}^{[\leqslant
h-1]\pm}(\bar{\Delta})+\psi_{\alpha,\sigma}^{[h]\pm}(\Delta).$ (55)

Note that (50) is then

$\psi_{\alpha,\sigma}^{\pm}(t,x)\equiv\psi_{\alpha,\sigma}^{[\leq
0]\pm}(\Delta^{(1)}(t,x)).$ (56)

We then define, for $h\in\\{-N_{\beta},\cdots 0\\}$,

$e^{\beta|\Lambda|c^{[h]}-\mathcal{V}^{[h-1]}(\psi^{\leqslant h-1]})}:=\int
P^{[h]}(d\psi^{[h]})\ e^{-\mathcal{V}^{[h]}(\psi^{[\leqslant
h]})},\quad\mathcal{V}^{[0]}(\psi^{[\leq 0]}):=\mathcal{V}(\psi^{[\leq 0]})$
(57)

in which $c^{[h]}\in\mathbb{R}$ is a constant and $\mathcal{V}^{[h-1]}$ has no
constant term. The function $\mathcal{V}^{[h]}$ is called the effective
potential on scale $h$, and it dictates the physical properties of the system
at distances $\sim 2^{-h}$. By a straightforward induction, we then find that

$\int P(d\psi)\
e^{-\mathcal{V}(\psi)}=\exp\left({\textstyle-\beta|\Lambda|\sum_{h=-N_{\beta}}^{0}c^{[h]}}\right).$
(58)

The specific free energy (37) is then

$f_{\Lambda}=f_{0,\Lambda}-\sum_{h=-N_{\beta}}^{0}c^{[h]}.$ (59)

We are then left with computing $\mathcal{V}^{[h]}$ and $c^{[h]}$ using (57).
By (38), $\mathcal{V}$ is local in $(t,x)$ and so, by induction, it takes the
form

$\mathcal{V}^{[h]}(\psi^{[\leqslant
h]})=-\sum_{\Delta\in\mathcal{Q}_{h}}v_{h}(\psi^{[\leqslant h]}(\Delta)).$
(60)

Because there are only four Grassmann fields and their conjugates per cell,
$v_{h}$ must be a polynomial in the Grassmann fields of order $\leqslant 8$.
In fact, by symmetry considerations, we find that $v_{h}$ must be of the form

$v_{h}(\psi)=\sum_{i=0}^{6}\ell_{i}^{(h)}O_{i}(\psi)$ (61)

with

$O_{0}(\psi):=\sum_{\sigma\in\\{\uparrow,\downarrow\\}}\left(\psi_{a,\sigma}^{+}\psi_{b,\sigma}^{-}+\psi_{b,\sigma}^{+}\psi_{a,\sigma}^{-}\right),\quad
O_{1}(\psi):=\sum_{\alpha\in\\{a,b\\}}\psi_{\alpha,\uparrow}^{+}\psi_{\alpha,\uparrow}^{-}\psi_{\alpha,\downarrow}^{+}\psi_{\alpha,\downarrow}^{-}$
(62)
$O_{2}(\psi):=\psi_{a,\uparrow}^{+}\psi_{a,\downarrow}^{-}\psi_{b,\downarrow}^{+}\psi_{b,\uparrow}^{-}+\psi_{b,\uparrow}^{+}\psi_{b,\downarrow}^{-}\psi_{a,\downarrow}^{+}\psi_{a,\uparrow}^{-}+\psi_{a,\downarrow}^{+}\psi_{a,\uparrow}^{-}\psi_{b,\uparrow}^{+}\psi_{b,\downarrow}^{-}+\psi_{b,\downarrow}^{+}\psi_{b,\uparrow}^{-}\psi_{a,\uparrow}^{+}\psi_{a,\downarrow}^{-}$
(63)
$O_{3}(\psi):=\sum_{\sigma\in\\{\uparrow,\downarrow\\}}\psi_{a,\sigma}^{+}\psi_{a,\sigma}^{-}\psi_{b,\sigma}^{+}\psi_{b,\sigma}^{-},\quad
O_{4}(\psi):=\psi_{a,\uparrow}^{+}\psi_{b,\uparrow}^{-}\psi_{a,\downarrow}^{+}\psi_{b,\downarrow}^{-}+\psi_{b,\uparrow}^{+}\psi_{a,\uparrow}^{-}\psi_{b,\downarrow}^{+}\psi_{a,\downarrow}^{-}$
(64) $\begin{array}[]{@{}>{\displaystyle}l@{}}\hphantom{\hskip
423.62pt}\\\\[-14.22636pt]
O_{5}(\psi):=\psi_{a,\uparrow}^{+}\psi_{a,\uparrow}^{-}\psi_{a,\downarrow}^{+}\psi_{b,\uparrow}^{-}\psi_{b,\uparrow}^{+}\psi_{b,\downarrow}^{-}+\psi_{a,\downarrow}^{+}\psi_{a,\downarrow}^{-}\psi_{a,\uparrow}^{+}\psi_{b,\downarrow}^{-}\psi_{b,\downarrow}^{+}\psi_{b,\uparrow}^{-}+\\\\[8.5359pt]
\hfill+\psi_{b,\uparrow}^{+}\psi_{b,\uparrow}^{-}\psi_{b,\downarrow}^{+}\psi_{a,\uparrow}^{-}\psi_{a,\uparrow}^{+}\psi_{a,\downarrow}^{-}+\psi_{b,\downarrow}^{+}\psi_{b,\downarrow}^{-}\psi_{b,\uparrow}^{+}\psi_{a,\downarrow}^{-}\psi_{a,\downarrow}^{+}\psi_{a,\uparrow}^{-}\end{array}$
(65)
$O_{6}(\psi):=\psi_{a,\uparrow}^{+}\psi_{a,\uparrow}^{-}\psi_{a,\downarrow}^{+}\psi_{a,\downarrow}^{-}\psi_{b,\uparrow}^{+}\psi_{b,\uparrow}^{-}\psi_{b,\downarrow}^{+}\psi_{b,\downarrow}^{-}.$
(66)

###  3.2 Beta function of the hierarchical model

8 -  Beta function. We have thus introduced a strategy to compute
$\mathcal{V}^{[h]}$ inductively: starting from $\mathcal{V}^{[h]}$, we have

$e^{\beta|\Lambda|c^{[h]}-\mathcal{V}^{[h-1]}(\psi^{[\leqslant h-1]})}=\int
P(d\psi^{[h]})\ e^{-\mathcal{V}^{[h]}(\psi^{[h]}+2^{-\gamma}\psi^{[\leqslant
h-1]})}$ (67)

where $\gamma\equiv 1$ is the scaling dimension of $\psi$ in (55). Now, by
(60), is

$e^{\beta|\Lambda|c^{[h]}+\sum_{\bar{\Delta}\in\mathcal{Q}_{h-1}}v_{h-1}(\psi^{[\leqslant
h-1]}(\bar{\Delta}))}=\prod_{\Delta\in\mathcal{Q}_{h}}\int
P(d\psi^{[h]}(\Delta))\
e^{v_{h}(\psi^{[h]}(\Delta)+2^{-\gamma}\psi^{[\leqslant
h-1]}(\bar{\Delta}))}.$ (68)

If we group the right side in boxes on scale $h-1$, we find

$e^{\beta|\Lambda|c^{[h]}+\sum_{\bar{\Delta}\in\mathcal{Q}_{h-1}}v_{h-1}(\psi^{[\leqslant
h-1]}(\bar{\Delta}))}=\prod_{\bar{\Delta}\in\mathcal{Q}_{h-1}}\prod_{\displaystyle\mathop{\scriptstyle\Delta\in\mathcal{Q}_{h}}_{\Delta\subset\bar{\Delta}}}\int
P(d\psi^{[h]}(\Delta))\
e^{v_{h}(\psi^{[h]}(\Delta)+2^{-\gamma}\psi^{[\leqslant
h-1]}(\bar{\Delta}))}.$ (69)

In addition, the integral over $\psi^{[h]}(\Delta)$ does not depend on
$\Delta$. Therefore, using the fact that $\bar{\Delta}$ contains
$2^{d+1}\equiv 8$ ($d\equiv 2$ is the dimension of the lattice) boxes
$\Delta$, we have, for all $\bar{\Delta}\in\mathcal{Q}_{h-1}$ and for any
choice of $\Delta\in\mathcal{Q}_{h}$ with $\Delta\subset\bar{\Delta}$,

$e^{\beta|\Lambda|c^{[h]}+v_{h-1}(\psi^{[\leqslant
h-1]}(\bar{\Delta}))}=\left(\int P(d\psi^{[h]}(\Delta))\
e^{v_{h}(\psi^{[h]}(\Delta)+2^{-\gamma}\psi^{[\leqslant
h-1]}(\bar{\Delta}))}\right)^{2^{d+1}}.$ (70)

We expand the exponential and use (61):

$\begin{array}[]{@{}>{\displaystyle}l@{}}\hphantom{\hskip
423.62pt}\\\\[-14.22636pt]
\beta|\Lambda|c^{[h]}+\sum_{i=0}^{6}\ell_{i}^{(h-1)}O_{i}(\psi^{[\leqslant
h-1]}(\bar{\Delta}))=\\\ \hfill=2^{d+1}\log\int
P(d\psi^{[h]}(\Delta))\sum_{n=0}^{\infty}\frac{1}{n!}\left(\sum_{i=0}^{6}\ell_{i}^{(h)}O_{i}\left(\psi^{[h]}(\Delta)+2^{-\gamma}\psi^{[\leqslant
h-1]}(\bar{\Delta})\right)\right)^{n}.\end{array}$ (71)

The computation is thus reduced to computing the map
$\ell^{(h)}\mapsto\ell^{(h-1)}$ using (71). The coefficients $\ell_{i}^{(h)}$
are called running coupling constants, and the map
$\ell^{(h)}\mapsto\ell^{(h-1)}$ is called the beta function of the model. The
running coupling constants play a very important role, as they specify the
effective potential on scale $h$, and thereby the physical properties of the
system at distances $\sim 2^{-h}$.

9 -  Feynman diagrams. Having defined the hierarchical model as we have, the
infinite sum in (71) is actually finite ($n\leqslant 4$), so to compute the
beta function, it suffices to compute Gaussian Grassmann integrals of a finite
number of Grassmann monomials. A convenient way to carry out this computation
is to represent each term graphically, using Feynman diagrams. First, let us
expand the power $n$ and graphically represent the terms that must be
integrated. For each $n$, we have $n$ possible choices of $\ell_{i}O_{i}$.
Now, $O_{i}$ can be quadratic in $\psi$ ($O_{0}$), quartic ($O_{1}$, $O_{2}$,
$O_{3}$, $O_{4}$), sextic ($O_{5}$) or octic ($O_{6}$). We will represent
$O_{i}$ by a vertex with the label $i$, from which two, four, six or eight
edges emanate, depending on the degree of $O_{i}$. Each edge corresponds to a
factor $\psi^{[h]}+2^{-\gamma}\psi^{[\leqslant h-1]}$. For each edge, we can
either choose the term $2^{-\gamma}\psi^{[\leqslant h-1]}$, which is not
integrated, and so the edge will be called external, or we can choose the term
$\psi^{[h]}$, in which case this edge will have to be integrated, and will be
called internal. Having made all of these choices, we have a set of vertices
as well as external and internal edges. We must now integrate $\psi^{[h]}$.
Recall that, by Wick’s rule (31), to integrate a monomial, we pair up
$\psi^{[h]}$’s (internal edges), and multiply the corresponding propagators.
Graphically, this is done by connecting internal edges. Each connection made
stands in for a propagator. The computation of the Gaussian Grassmann integral
in (71) can thus be reduced to a sum over graphs.

Let us make a few comments about this expansion, and how it is used for non-
hierarchical models.

  * •

The diagram expansion for the hierarchical model is finite, in that there are
only finitely many types of vertices, and finitely many of them, so the number
of possible graphs is finite. This diagrammatic representation can be done in
the non-hierarchical model as well, but there, the number of graphs is
infinite. In fact, even considering just one type of vertex, the number of
graphs grows as $(n!)^{2}$ as the number of vertices $n$ tends to $\infty$.
This makes the infinite sum in (71) absolutely divergent. However, in counting
these graphs, we are ignoring important signs: in the Wick rule (31), each
permutation comes with its signature. These signs lead to important
cancellations in the sum over $n$, which can be exploited using so-called Gram
bounds, or determinant expansions [BF78, BF84] (these are even important
numerically, as they can be used to compute the same integrals with many fewer
terms). Here, we will focus on the hierarchical model for which such
considerations are not needed. The Gram bounds for graphene are worked out in
detail in [Gi10, Appendix B].

  * •

The diagram expansion allows us to compute the integral in (71). One then has
to take the logarithm, which can also be done by an expansion. The sign in the
expansion of the logarithm leads to significant cancellations, which, it turns
out, cancels all of the graphs that are not connected. We will not dwell on
this fact here, as this is not so relevant for the hierarchical model, but
this is an important fact for non-hierarchical models.

  * •

The diagram expansion we have discussed performs the integral on a single
scale. In the grander scheme, the external edges will get contracted on lower
scales, and one can construct a larger diagram expansion that covers all
scales. Understanding the combinatorial properties of this full expansion
requires some extra tools to keep track of the scales of the edges. One way of
doing this is to use Gallavotti-Nicoló [GN85] trees, as detailed in [Gi10].
Another good reference for the tree expansion and its connection to Feynman
diagrams is [GJ16, Section 5] This is not needed for the hierarchical model.

10 -  Power counting. Returning to our hierarchical model, let us consider
some of the more important contributions to the Feynman diagram expansion:
bare vertices. These are the simplest possible graphs, in which we have one
vertex, and all edges are external. In other words, no integrating is taking
place. Let us denote the number of external edges by $2l$, which can either be
2, 4, 6 or 8. The contribution of this graph is (keeping track of the
$2^{d+1}$ factor in (71))

$2^{d+1-2l\gamma}\ell_{i}^{(h)}.$ (72)

Furthermore, this graph will contribute to the running coupling constant
$\ell_{i}$, and so, on scale $h-1$, we will have

$\ell_{i}^{(h-1)}=2^{d+1-2l\gamma}\ell_{i}^{(h)}+\cdots$ (73)

(the $\cdots$ stand in for the other terms contributing the this running
coupling constant, and will be discarded for the moment). Three things can
happen:

  * •

if $d+1-2l\gamma<0$, then this term will decrease exponentially as
$h\to-\infty$; such running coupling constants will be called irrelevant;

  * •

if $d+1-2l\gamma>0$, then this term will grow exponentially as $h\to-\infty$;
such running coupling constants will be called relevant;

  * •

if $d+1-2l\gamma=0$, then not much can be said without looking at the
neglected terms; such running coupling constants will be called marginal.

In the case of graphene, $d=2$ and $\gamma=1$, so if $2l=2$, then the running
coupling constant is relevant, and if $2l\geqslant 4$, it is irrelevant; there
are no marginal constants.

But what of the other terms we neglected? The notions of relevant, irrelevant
and marginal terms only apply to a perturbative setting, in which the running
coupling constants are small, and the question is whether they stay small. If
all running coupling constants are irrelevant, then if they start small at
scale $h=0$, they stay that way for all scales. When some are relevant or
marginal, there is a risk for some to grow and leave the perturbative regime.
Since the running coupling constants dictate the physical properties of the
system, leaving the perturbative regime means that, on the scales on which the
constants grow, the system will start behaving radically differently from the
non-interacting one (in which all running coupling constants are 0). In other
words, studying whether running coupling constants are relevant, irrelevant or
marginal, gives us information about whether the interacting system can be
approximated by the non-interacting one or not.

This power counting holds for the non-hierarchical model as well, as detailed
in [Gi10]. For a more general treatment of power counting in Fermionic models
with point-singularities, see [GJ16, Section 5.2].

In the case of graphene, we have one relevant coupling: $O_{0}$, which is
quadratic in the Grassmann fields. This is the only relevant coupling, and all
others stay small. However, since the relevant coupling is quadratic, it
merely shifts the non-interacting system (whose Hamiltonian is quadratic in
the Grassmann fields) to another system with a quadratic (that is, non-
interacting) Hamiltonian. Thus the relevant coupling does not imply that the
interactions are preponderant, but rather that the interaction terms shifts
the system from one non-interacting system to another. Since graphene only has
one relevant coupling, and that one is quadratic, graphene is called super-
renormalizable.

11 -  Hierarchical beta function. As was mentioned above, the beta function
can be computed explicitly for the hierarchical model, so the claims in the
previous paragraph can be verified rather easily. The exact computation
involves many terms, but it can be done easily using the meankondo software
package [mk]. The resulting beta function contains 888 terms, and will not be
written out here. A careful analysis of the beta function shows that there is
an equilibrium point at $\ell_{i}=0$ for $i=1,2,3,4,5,6$ and

$\ell_{0}\in\\{0,1\\}.$ (74)

The point with $\ell_{0}=0$ is unstable, whereas $\ell_{0}=1$ is stable.

Figure 7:  The projection of the directional vector field of the beta function
for hierarchical graphene onto the $(\ell_{0},\ell_{1})$ plane. (Each arrow
shows the direction of the vector field, the color corresponds to the
logarithm of the amplitude, with red being larger and blue smaller.) The
stable equilibrium point at $\ell_{0}=1$ and $\ell_{i}=0$ is clearly visible.

###  3.3 Results for non-hierarchical graphene

Using the Gram bounds and Gallavotti-Nicoló tree expansion hinted at above,
the renormalization group analysis can be carried out in full for (non-
hierarchical) graphene [Gi10, GM10]. The result is quite similar to what we
have found for the hierarchical model: the effective potential contains only
irrelevant terms, except for the quadratic terms, which one absorbs into the
quadratic non-interacting Hamiltonian. We thus find that the observables for
the interacting model are similar to those without interactions, provided the
free Hamiltonian is suitably changed (which is called renormalizing the non-
interacting Hamiltonian). In particular, the Fermi velocity $v_{F}$ in (16) is
renormalized.

This renormalization can also be quantified: the renormalization group
analysis provides a power series expansion in $U$, and, when using the Gram
bounds, this series is shown to be absolutely convergent. (In other words, the
observables are analytic in $U$.) Thus, we can estimate the effect of the
interaction by truncating the power series and bounding its remainder.

In the case of the non-hierarchical graphene model, the renormalization of the
free Hamiltonian is restricted by symmetries. This is a crucial ingredient of
the construction. Indeed, while iterating the renormalization group flow, the
free Hamiltonian is changed. However, the philosophy of the scale
decomposition is anchored in the singularities of the (inverse of the) non-
interacting Hamiltonian, so if these singularities were to change under the
flow, the entire theory might break down. The symmetries prevent that from
happening: in particular, they ensure that the singularities $p_{F}^{\omega}$
remain point singularities (they could a priori, turn them into curves). This
is crucial, as extended singularities would change the power counting, and
would not allow the renormalization group flow to extend to arbitrarily small
scales. (The symmetries also ensure that the singularities do not move, but if
they did, this could be dealt with by defining the scale decomposition
slightly differently, see [GJ16] for an analysis of bilayer graphene where
this happens).

##  4 Hierarchical Kondo model

The graphene model discussed above is perturbative, in the sense that the
interactions do not qualitatively change the behavior of the system. Let us
now discuss a different model where the interactions change the behavior of
the system drastically: the Kondo model. The Kondo model, first studied by
Kondo [Ko64], was one of the foundational models in the development stages of
the renormalization group [Wi75]. It is integrable [An80], but is nevertheless
an interesting model to study renormalization group flows in. The hierarchical
Kondo model was introduced and studied in [BGJ15, GJ15].

###  4.1 Hamiltonian

The Kondo model is a one-dimensional model of spin-$\frac{1}{2}$ electrons on
a chain that interact with a localized magnetic impurity, which is represented
as a spin-$\frac{1}{2}$ particle located at $x=0$. The Hilbert space for the
electrons is the usual Fock space $\mathcal{F}$, and the Hilbert space for the
impurity is $\mathbb{C}^{2}$. The Hamiltonian is split into the free
Hamiltonian and the interaction term:

$\bar{\mathcal{H}}_{0}=\bar{\mathcal{H}}_{0}+\bar{\mathcal{H}_{I}}$ (75)

with

$\bar{\mathcal{H}}_{0}=-\frac{1}{2}\sum_{\sigma\in\\{\uparrow,\downarrow\\}}\sum_{x\in\\{-\frac{L}{2}+1,\cdots,\frac{L}{2}\\}}a_{x,\sigma}^{\dagger}(a_{x+1,\sigma}+a_{x-1,\sigma})\otimes\mathds{1}$
(76)

and

$\bar{\mathcal{H}}_{I}=-U\sum_{\sigma_{1},\sigma_{2}\in\\{\uparrow,\downarrow\\}}\sum_{j=1,2,3}a_{0,\sigma_{1}}^{\dagger}S^{(j)}_{\sigma_{1},\sigma_{2}}a_{0,\sigma_{2}}\otimes
S^{(j)}$ (77)

where $a_{x,\sigma}$ is the annihilation operator at $x$ with spin $\sigma$
and $S^{(j)}$ are Pauli matrices:

$S^{(1)}=\left(\begin{array}[]{cc}0&1\\\ 1&0\end{array}\right),\quad
S^{(2)}=\left(\begin{array}[]{cc}0&-i\\\ i&0\end{array}\right),\quad
S^{(3)}=\left(\begin{array}[]{cc}1&0\\\ 0&-1\end{array}\right).$ (78)

In Fourier space, let

$\hat{\Lambda}_{L}:=\frac{2\pi}{L}\left\\{-\frac{L}{2}+1,\cdots,\frac{L}{2}\right\\}$
(79)

in term of which

$\bar{\mathcal{H}}_{0}=-\sum_{\sigma\in\\{\uparrow,\downarrow\\}}\frac{1}{L}\sum_{k\in\hat{\Lambda}_{L}}\cos(k)\hat{a}_{k,\sigma}^{\dagger}\hat{a}_{k,\sigma}\otimes\mathds{1}.$
(80)

The propagator is thus

$g(t,x)=\frac{1}{\beta
L}\sum_{k_{0}\in\frac{2\pi}{\beta}(\mathbb{Z}+\frac{1}{2})}\sum_{k\in\hat{\Lambda}_{L}}\
e^{-ik_{0}t}e^{-ikx}\frac{1}{-ik_{0}-\cos(k)}.$ (81)

###  4.2 Scale decomposition and hierarchical model

Since the interaction is localized at $x=0$, we only need propagators at
$x=0$:

$g(t,0)=\frac{1}{\beta
L}\sum_{k_{0}\in\frac{2\pi}{\beta}(\mathbb{Z}+\frac{1}{2})}\sum_{k\in\hat{\Lambda}_{L}}\
e^{-ik_{0}t}\frac{1}{-ik_{0}-\cos(k)}.$ (82)

We decompose the propagator into scales around the singularities
$(k_{0},k)=(0,\pm\frac{\pi}{2})$:

$g^{(h,\omega)}(t,0)=\frac{1}{\beta
L}\sum_{k_{0}\in\frac{2\pi}{\beta}(\mathbb{Z}+\frac{1}{2})}\sum_{k\in\hat{\Lambda}_{L}}\
e^{-ik_{0}t}\frac{1}{-ik_{0}-\cos(k)}\Phi_{h}(k_{0},k)$ (83)

where $\Phi_{h}$ is defined as in (41)

$\Phi_{h}(k_{0},k):=(\chi_{0}(2^{-h}|(k_{0},k-\omega{\textstyle\frac{\pi}{2}})|)-\chi_{0}(2^{-h+1}|(k_{0},k-\omega{\textstyle\frac{\pi}{2}})|)).$
(84)

On scale $h$ we have

$\frac{1}{|ik_{0}+\cos(k)|}=O(2^{-h}),\quad\sum_{k_{0},k}\Phi_{h}(k_{0},k)=O(2^{2h})$
(85)

so

$g^{(h,\omega)}(t,0)=O(2^{h})$ (86)

(as opposed to $2^{2h}$ for graphene).

We define the hierarchical model in a very similar way to graphene, except
that, since only the point $x=0$ matters, there is no spatial dependence. We
define the boxes as

$\mathcal{Q}_{h}:=\left\\{[i2^{|h|},(i+1)2^{|h|})\right\\}_{i\in\mathbb{Z}}.$
(87)

We further split each box into two halves, which is necessary otherwise the
model would come out trivial: for each box $\Delta=[i2^{|h|},(i+1)2^{|h|+1})$,
we define

$\Delta_{-}=[i2^{|h|},(i+{\textstyle\frac{1}{2}})2^{|h|}),\quad\Delta_{+}=[(i+{\textstyle\frac{1}{2}})2^{|h|},(i+1)2^{|h|})$
(88)

In accordance with the scaling in (86), we define, for
$\sigma\in\\{\uparrow,\downarrow\\}$, $h\in\\{-N_{\beta},\cdots,-1\\}$,
$\Delta\in\mathcal{Q}_{h}$, $\eta\in\pm$

$\psi_{\sigma}^{[\leqslant
h]\pm}(\Delta_{\eta}):=\frac{1}{\sqrt{2}}\psi^{[\leqslant
h-1]\pm}(\Delta)+\psi_{\sigma}^{[h]\pm}(\Delta_{\eta})$ (89)

(compare this to (55) in which the scaling factor is $\frac{1}{2}$ instead of
$\frac{1}{\sqrt{2}}$).

###  4.3 Beta function

For this model, it is more convenient to expand the exponential of the
potential: we denote

$e^{-\mathcal{V}^{[h]}(\psi)}=\mathcal{W}^{[h]}(\psi).$ (90)

One can show that the effective potentials $\mathcal{W}$ are of the form

$\mathcal{W}^{[h]}(\psi^{[\leqslant
h]})=\prod_{\Delta\in\mathcal{Q}_{h}}\prod_{\eta=\pm}(1+w_{h}(\psi^{[\leqslant
h]}(\Delta_{\eta}))).$ (91)

For this model, there are only two running coupling constants:

$w_{h}(\psi)=\sum_{i=0}^{1}\ell_{i}^{(h)}O_{i}(\psi)$ (92)

and

$O_{0}(\psi):=\frac{1}{2}\sum_{\sigma,\sigma^{\prime}\in\\{\uparrow,\downarrow\\}}\sum_{j=1,2,3}\psi_{\sigma}^{+}S^{(j)}_{\sigma,\sigma^{\prime}}\psi_{\sigma^{\prime}}^{-}\otimes
S^{(j)},\quad
O_{1}(\psi):=\frac{1}{2}\left(\sum_{\sigma,\sigma^{\prime}\in\\{\uparrow,\downarrow\\}}\sum_{j=1,2,3}\psi_{\sigma}^{+}S^{(j)}_{\sigma,\sigma^{\prime}}\psi_{\sigma^{\prime}}^{-}\right)^{2}\otimes\mathds{1}$
(93)

(the $O$ are a product of Grassmann variables and a matrix acting on the
impurity).

We compute the beta function in the same way as for graphene, see (71):

$\begin{array}[]{@{}>{\displaystyle}l@{}}\hphantom{\hskip
423.62pt}\\\\[-14.22636pt]
C^{(h)}\left(1+\sum_{i=0}^{1}\ell_{i}^{(h-1)}O_{i}(\psi^{[\leqslant
h-1]}(\Delta))\right)=\\\ \hfill=\prod_{\eta=\pm}\int
P(d\psi^{[h]}(\Delta_{\eta})\left(\sum_{i=0}^{1}\ell_{i}^{(h)}O_{i}\left(\psi^{[h]}(\Delta)+2^{-\gamma}\psi^{[\leqslant
h-1]}(\Delta_{\eta})\right)\right)^{2}\end{array}$ (94)

with $\gamma\equiv\frac{1}{2}$. In this expression, the square plays the same
role as the prefactor $2^{d+1}$ in (71). The power counting is determined by
the sign of

$1+2l\gamma\equiv 1+l$ (95)

so quadratic terms ($\ell_{0}$) are marginal and quartic terms ($\ell_{1}$)
are irrelevant.

The beta function can be computed exactly for the hierarchical model:

$\begin{array}[]{r@{\
}>{\displaystyle}l}C^{(h)}=&1+\frac{3}{2}(\ell_{0}^{(h)})^{2}+9(\ell_{1}^{(h)})^{2}\\\\[8.5359pt]
\ell_{0}^{(h-1)}=&\frac{1}{C^{(h)}}\Big{(}\ell_{0}^{(h)}+3\ell_{0}^{(h)}\ell_{1}-(\ell_{0}^{(h)})^{2}\Big{)}\\\\[8.5359pt]
\ell_{1}^{(h-1)}=&\frac{1}{C^{(h)}}\Big{(}\frac{1}{2}\ell_{1}^{(h)}+\frac{1}{8}(\ell_{0}^{(h)})^{2}\Big{)}.\end{array}$
(96)

One can show that there are two equilibrium points: $\ell_{0}=\ell_{1}=0$, as
well as

$\ell^{*}_{0}=-x_{0}\frac{1+5x_{0}}{1-4x_{0}},\quad\ell^{*}_{1}=\frac{x_{0}}{3}$
(97)

where $x_{0}\approx 0.15878626704216...$ is the real root of
$4-19x-22x^{2}-107x^{3}=0$. The trivial equilibrium at $(0,0)$ is stable if
and only if $\ell_{0}>0$ (that is, if it is approached from the right).
Otherwise, the flow goes to the non-trivial equilibrium point, see figure 8

Figure 8:  The directional vector field of the beta function for the
hierarchical Kondo model. (Each arrow shows the direction of the vector field,
the color corresponds to the logarithm of the amplitude, with red being larger
and blue smaller.) The stable and unstable equilibrium points are clearly
visible.

Thus, for the hierarchical Kondo model, the flow can go to a non-trivial
equilibrium point, at which the interaction is preponderant. Physically, this
translates to the behavior of the interacting system being qualitatively
different from the non-interacting one (in this case, this manifests itself by
the magnetic susceptibility of the impurity being finite at $0$-temperature,
which means that the impurity resists an external magnetic field even in the
ground state, which is not the case for the non-interacting system).

##  Appendix A Free Fermions

Let

$\mathcal{H}_{0}=\sum_{i,j=1}^{N}\mu_{i,j}a^{\dagger}_{i}a_{j}$ (98)

where $(a_{i})_{i}$ is a family of Fermionic annihilation operators and
$\mu_{i,j}=\mu_{j,i}^{*}$. Given an operator $A$ on Fock space and $\beta>0$,
let

$\left<A\right>:=\frac{\mathrm{Tr}(e^{-\beta\mathcal{H}_{0}}A)}{\mathrm{Tr}(e^{-\beta\mathcal{H}_{0}})}$
(99)

and, for $t\in[0,\beta)$, let

$a_{i}^{-}(t):=e^{t\mathcal{H}_{0}}a_{i}e^{-t\mathcal{H}_{0}},\quad
a_{i}^{+}(t):=e^{t\mathcal{H}_{0}}a_{i}^{\dagger}e^{-t\mathcal{H}_{0}}.$ (100)

Furthermore, we define the Fermionic time ordering operator as a linear
operator on polynomials of $a_{i}(t)$ such that, for any
$n\in\\{1,\cdots,N\\}$, $j_{1},\cdots,j_{n}\in\\{1,\cdots,N\\}$,
$\omega_{1},\cdots,\omega_{n}\in\\{-,+\\}$, $t_{1},\cdots,t_{n}\in[0,\beta)$,

$\mathbf{T}\left(\prod_{i=1}^{n}a^{\omega_{i}}_{j_{i}}(t_{i})\right)=(-1)^{\tau}\prod_{i=1}^{n}a^{\omega_{\tau(i)}}_{j_{\tau(i)}}(t_{\tau(i)})$
(101)

where $\tau$ is the permutation of $\\{1,\cdots,n\\}$ that orders the
$t_{i}$’s in reverse order, and in the case of equality of $t$’s, $\tau$
places $\omega=-$ to the left of $\omega=+$ and orders the $j$’s. Formally,
$\tau$ is such that $t_{\tau(i)}\geqslant t_{\tau(i+1)}$ and if
$t_{\tau(i)}=t_{\tau(i+1)}$, then either $\omega_{\tau(i)}=-$ and
$\omega_{\tau(i+1)}=+$ or $j_{\tau(i)}<j_{\tau(i+1)}$.

Let us first prove a technical lemma.

Lemma A.1

For $t\in\mathbb{R}$, and $\lambda_{1},\cdots,\lambda_{N}\in\mathbb{R}$,

$e^{-t\sum_{i=1}^{N}\lambda_{i}a_{i}^{\dagger}a_{i}}=\prod_{i=1}^{N}(1+(e^{-t\lambda_{i}}-1)a_{i}^{\dagger}a_{i})$
(102)

and

$e^{t\sum_{i=1}^{N}\lambda_{i}a_{i}^{\dagger}a_{i}}a_{j}^{\dagger}e^{-t\sum_{i=1}^{N}\lambda_{i}a_{i}^{\dagger}a_{i}}=e^{t\lambda_{j}}a_{j}^{\dagger},\quad
e^{t\sum_{i=1}^{N}\lambda_{i}a_{i}^{\dagger}a_{i}}a_{j}e^{-t\sum_{i=1}^{N}\lambda_{i}a_{i}^{\dagger}a_{i}}=e^{-t\lambda_{j}}a_{j}.$
(103)

Proof: To prove (102), we write
$e^{-t\sum_{i}\lambda_{i}a_{i}^{\dagger}a_{i}}=\prod_{i}e^{-t\lambda_{i}a_{i}^{\dagger}a_{i}}$
and expand the exponential, using the fact that
$(a_{i}^{\dagger}a_{i})^{n}=a_{i}^{\dagger}a_{i}$ for any $n\geqslant 1$. By
(102),

$e^{-t\sum_{i=1}^{N}\lambda_{i}a_{i}^{\dagger}a_{i}}a_{j}^{\dagger}=\left(\prod_{i=1}^{n}(1+(e^{-t\lambda_{i}}-1)a_{i}^{\dagger}a_{i})\right)a_{j}^{\dagger}$
(104)

and, commuting $a_{j}^{\dagger}$ through using
$\\{a_{i},a_{j}^{\dagger}\\}=\delta_{i,j}$ and $(a_{j}^{\dagger})^{2}=0$, we
find

$e^{-t\sum_{i=1}^{N}\lambda_{i}a_{i}^{\dagger}a_{i}}a_{j}^{\dagger}=e^{-t\lambda_{j}}a_{j}^{\dagger}\prod_{i\neq
j}(1+(e^{-t\lambda_{i}}-1)a_{i}^{\dagger}a_{i})$ (105)

and since $(a_{j}^{\dagger})^{2}=0$,

$e^{-t\sum_{i=1}^{N}\lambda_{i}a_{i}^{\dagger}a_{i}}a_{j}^{\dagger}=e^{-t\lambda_{j}}a_{j}^{\dagger}\prod_{i}(1+(e^{-t\lambda_{i}}-1)a_{i}^{\dagger}a_{i})=e^{-t\lambda_{j}}a_{j}^{\dagger}e^{-t\sum_{i=1}^{N}\lambda_{i}a_{i}^{\dagger}a_{i}}.$
(106)

Taking the $\dagger$ of (106), we find

$a_{j}e^{-t\sum_{i=1}^{N}\lambda_{i}a_{i}^{\dagger}a_{i}}=e^{-t\lambda_{j}}e^{-t\sum_{i=1}^{N}\lambda_{i}a_{i}^{\dagger}a_{i}}a_{j}.$
(107)

Combining (106) and (107), we find (103). $\square$

###  A.1 Two-point correlation function

We now compute the two-point correlation function of $\left<\cdot\right>$.

Lemma A.2

For $i,j\in\\{1,\cdots,N\\}$ and $t,t^{\prime}\in[0,\beta)$, if $t\neq
t^{\prime}$, then

$s_{i,j}(t-t^{\prime}):=\left<\mathbf{T}(a_{i}^{-}(t)a_{j}^{+}(t^{\prime}))\right>=\frac{1}{\beta}\sum_{k_{0}\in\frac{2\pi}{\beta}(\mathbb{Z}+\frac{1}{2})}e^{-ik_{0}(t-t^{\prime})}(-ik_{0}\mathds{1}+\mu)^{-1}_{i,j}$
(108)

and if $t=t^{\prime}$,

$s_{i,j}(0):=\left<\mathbf{T}(a_{i}^{-}(t)a_{j}^{+}(t))\right>=\frac{1}{2}\delta_{i,j}+\frac{1}{\beta}\sum_{k_{0}\in\frac{2\pi}{\beta}(\mathbb{Z}+\frac{1}{2})}(-ik_{0}\mathds{1}+\mu)^{-1}_{i,j}.$
(109)

Proof: We diagonalize $\mu$ by a unitary transform $U$, and define

$b_{i}:=\sum_{j=1}^{N}U_{j,i}^{*}a_{j},\quad
b_{i}^{\dagger}:=\sum_{j=1}^{N}U_{j,i}a_{j}^{\dagger}$ (110)

in terms of which

$e^{-t\mathcal{H}_{0}}=e^{-t\sum_{i=1}^{N}\lambda_{i}b_{i}^{\dagger}b_{i}}=\prod_{i=1}^{N}(1+(e^{-t\lambda_{i}}-1)b_{i}^{\dagger}b_{i})$
(111)

where $\lambda_{i}$ are the eigenvalues of $\mu$. Note that

$\\{b_{i},b_{j}^{\dagger}\\}=\sum_{k,l=1}^{N}U_{k,i}^{*}U_{l,j}\\{a_{k},a_{l}^{\dagger}\\}=\sum_{k=1}^{N}U_{k,i}^{*}U_{k,j}=\delta_{i,j}.$
(112)

Thus

$\mathrm{Tr}(e^{-\beta\mathcal{H}_{0}})=\prod_{i=1}^{N}(1+e^{-\beta\lambda_{i}}).$
(113)

Let

$b_{i}^{-}(s):=e^{t\mathcal{H}_{0}}b_{i}e^{-t\mathcal{H}_{0}}=\sum_{j=1}^{N}U_{j,i}^{*}a_{j}^{-}(s),\quad
b_{i}^{+}(s):=e^{t\mathcal{H}_{0}}b_{i}^{\dagger}e^{-t\mathcal{H}_{0}}=\sum_{j=1}^{N}U_{j,i}a_{j}^{+}(s).$
(114)

By (103),

$b_{i}^{-}(s)=e^{-s\lambda_{i}}b_{i},\quad
b_{i}^{+}(s)=e^{s\lambda_{i}}b_{i}^{\dagger}.$ (115)

1 -  Now, if $t\geqslant t^{\prime}$,

$\mathbf{T}(a_{i}^{-}(t)a_{j}^{+}(t^{\prime}))=a_{i}^{-}(t)a_{j}^{+}(t^{\prime})=\sum_{k,l=1}^{N}U_{i,k}U_{j,l}^{*}e^{-t\lambda_{k}+t^{\prime}\lambda_{l}}b_{k}b_{l}^{\dagger}.$
(116)

In order for the trace to be non-zero, we must take $k=l$, so

$\mathrm{Tr}(e^{-\beta\mathcal{H}_{0}}\mathbf{T}(a_{i}^{-}(t)a_{j}^{+}(t^{\prime})))=\sum_{k=1}^{N}U_{i,k}U_{j,k}^{*}e^{-(t-t^{\prime})\lambda_{k}}\mathrm{Tr}(e^{-\beta\mathcal{H}_{0}}b_{k}b_{k}^{\dagger})$
(117)

and since $b_{k}b_{k}^{\dagger}=1-b_{k}^{\dagger}b_{k}$,

$\mathrm{Tr}(e^{-\beta\mathcal{H}_{0}}\mathbf{T}(a_{i}^{-}(t)a_{j}^{+}(t^{\prime})))=\sum_{k=1}^{N}U_{i,k}U_{j,k}^{*}e^{-(t-t^{\prime})\lambda_{k}}\prod_{l\neq
k}(1+e^{-\beta\lambda_{l}}).$ (118)

Therefore,

$s_{i,j}(t-t^{\prime})=\sum_{k=1}^{N}U_{i,k}U_{j,k}^{*}\frac{e^{-(t-t^{\prime})\lambda_{k}}}{1+e^{-\beta\lambda_{k}}}.$
(119)

2 -  If $t<t^{\prime}$, then

$\mathbf{T}(a_{i}^{-}(t)a_{j}^{+}(t^{\prime}))=-a_{j}^{+}(t^{\prime})a_{i}^{-}(t)=-e^{t^{\prime}\mathcal{H}_{0}}a_{j}^{\dagger}e^{-(t^{\prime}-t)\mathcal{H}_{0}}a_{i}e^{-t\mathcal{H}_{0}}$
(120)

so, using the cyclicity of the trace,

$\mathrm{Tr}(e^{-\beta\mathcal{H}_{0}}\mathbf{T}(a_{i}^{-}(t)a_{j}^{+}(t^{\prime})))=-\mathrm{Tr}(e^{t\mathcal{H}_{0}}a_{i}e^{-(t+\beta-t^{\prime})\mathcal{H}_{0}}a_{j}^{\dagger}e^{-t^{\prime}\mathcal{H}_{0}})$
(121)

which is identical to the trace computed above, but with $t$ replaced by
$t+\beta$, so

$s_{i,j}(t-t^{\prime})=-\sum_{k=1}^{N}U_{i,k}U_{j,k}^{*}\frac{e^{-(t-t^{\prime}+\beta)\lambda_{k}}}{1+e^{-\beta\lambda_{k}}}.$
(122)

3 -  Let us compute the Fourier transform of $s_{i,j}(t-t^{\prime})$ with
respect to $t-t^{\prime}\in(-\beta,\beta)$: for
$k_{0}\in\frac{\pi}{\beta}\mathbb{Z}$, by (119) and (122),

$\frac{1}{2}\int_{-\beta}^{\beta}d\tau\
e^{ik_{0}\tau}s_{i,j}(t-t^{\prime})=\sum_{k=1}^{N}U_{i,k}U_{j,k}^{*}\frac{1}{1+e^{-\beta\lambda_{k}}}\left(\int_{0}^{\beta}d\tau\
e^{ik_{0}\tau}e^{-\tau\lambda_{k}}-\int_{-\beta}^{0}d\tau\
e^{ik_{0}\tau}e^{-(\tau+\beta)\lambda_{k}}\right)$ (123)

and so

$\frac{1}{2}\int_{-\beta}^{\beta}d\tau\
e^{ik_{0}\tau}s_{i,j}(t-t^{\prime})=\sum_{k=1}^{N}\frac{U_{i,k}U_{j,k}^{*}}{1+e^{-\beta\lambda_{k}}}\left(\frac{e^{(ik_{0}-\lambda_{k})\beta}-1}{ik_{0}-\lambda_{k}}-\frac{e^{-\beta\lambda_{k}}-e^{-ik_{0}\beta}}{ik_{0}-\lambda_{k}}\right).$
(124)

Thus, if $k_{0}\in\frac{\pi}{\beta}2\mathbb{Z}$, then

$\frac{1}{2}\int_{-\beta}^{\beta}d\tau\ e^{ik_{0}\tau}s_{i,j}(t-t^{\prime})=0$
(125)

and if $k_{0}\in\frac{2\pi}{\beta}(\mathbb{Z}+\frac{1}{2})$, then

$\frac{1}{2}\int_{-\beta}^{\beta}d\tau\
e^{ik_{0}\tau}s_{i,j}(t-t^{\prime})=\sum_{k=1}^{N}U_{i,k}U_{j,k}^{*}\frac{1}{-ik_{0}+\lambda_{k}}.$
(126)

Thus, if the Fourier transform can be inverted, then we have

$s_{i,j}(t-t^{\prime})=\frac{1}{\beta}\sum_{k_{0}\in\frac{2\pi}{\beta}(\mathbb{Z}+\frac{1}{2})}e^{-ik_{0}(t-t^{\prime})}(-ik_{0}\mathds{1}+\mu)^{-1}_{i,j}$
(127)

and this is the case where $s_{i,j}(t-t^{\prime})$ is continuous.

4 -  Obviously, $s_{i,j}$ is continuous at $t-t^{\prime}\neq 0$. At
$t-t^{\prime}=0$, we have

$\lim_{t-t^{\prime}\to 0_{+}}s_{i,j}(t-t^{\prime})=\lim_{t-t^{\prime}\to
0_{-}}s_{i,j}(t-t^{\prime})+\delta_{i,j}$ (128)

so $s_{i,j}$ is continuous if $i\neq j$. Now, if $i=j$ and $t=t^{\prime}$,
then by (119),

$s_{i,j}(0)=\sum_{k=1}^{N}U_{i,k}U_{i,k}^{*}\frac{1}{1+e^{-\beta\lambda_{k}}}.$
(129)

Now,

$s_{i,j}(t-t^{\prime})-\frac{1}{2}\mathrm{sgn}(t-t^{\prime})$ (130)

is continuous at $t-t^{\prime}=0$, and its Fourier transform is, by a similar
computation to that above,

$\frac{1}{2}\int_{-\beta}^{\beta}d\tau\
\left(s_{i,j}(\tau)-\frac{1}{2}\mathrm{sgn}(\tau)\right)=\sum_{k=1}^{N}U_{i,k}U_{j,k}^{*}\left(\frac{1}{-ik_{0}+\lambda_{k}}+\frac{1}{ik_{0}}\right)$
(131)

which is absolutely summable in $k_{0}$. Thus

$\lim_{\tau\to
0_{+}}\frac{1}{\beta}\sum_{k_{0}\in\frac{2\pi}{\beta}(\mathbb{Z}+\frac{1}{2})}e^{-ik_{0}\tau}\left((-ik_{0}\mathds{1}+\mu)^{-1}+\frac{1}{ik_{0}}\mathds{1}\right)=\frac{1}{\beta}\sum_{k_{0}\in\frac{2\pi}{\beta}(\mathbb{Z}+\frac{1}{2})}\left((-ik_{0}\mathds{1}+\mu)^{-1}+\frac{1}{ik_{0}}\mathds{1}\right)$
(132)

and so

$\lim_{\tau\to
0_{+}}\left(s_{i,j}(\tau)-\frac{1}{2}\mathrm{sgn}(\tau)\right)=\frac{1}{\beta}\sum_{k_{0}\in\frac{2\pi}{\beta}(\mathbb{Z}+\frac{1}{2})}(-ik_{0}\mathds{1}+\mu)^{-1}$
(133)

(the sum is to be understood as a principal part:
$\sum_{k_{0}}(-ik_{0}\mathds{1}+\mu)^{-1}\equiv\sum_{k_{0}}\mu(k_{0}^{2}\mathds{1}+\mu^{2})^{-1}$)
and thus

$s_{i,j}(0)=\frac{1}{2}+\frac{1}{\beta}\sum_{k_{0}\in\frac{2\pi}{\beta}(\mathbb{Z}+\frac{1}{2})}(-ik_{0}\mathds{1}+\mu)^{-1}.$
(134)

$\square$

###  A.2 Wick rule

This system of free Fermions satisfies the Wick rule, stated below.

Lemma A.3

(Wick rule)

For any $n\in\\{1,\cdots,N\\}$, $j_{1},\cdots,j_{n}\in\\{1,\cdots,N\\}$,
$\bar{j}_{1},\cdots,\bar{j}_{n}\in\\{1,\cdots,N\\}$,
$t_{1},\cdots,t_{n}\in[0,\beta)$,
$\bar{t}_{1},\cdots,\bar{t}_{n}\in[0,\beta)$,

$\left<\mathbf{T}\left(\prod_{i=1}^{n}a_{j_{i}}^{-}(t_{i})a_{\bar{j}_{i}}^{+}(\bar{t}_{i})\right)\right>=\sum_{\tau\in\mathcal{S}_{n}}(-1)^{\tau}\prod_{i=1}^{n}\left<\mathbf{T}\left(a_{i}^{-}(t_{i})a_{\bar{j}_{\tau(i)}}^{+}(\bar{t}_{\tau(i)})\right)\right>.$
(135)

Proof: First, note that

$\left<\mathbf{T}\left(\prod_{i=1}^{n}a_{j_{i}}^{-}(t_{i})a_{\bar{j}_{i}}^{+}(\bar{t}_{i})\right)\right>=\left<\mathbf{T}\left(\left(\prod_{i=1}^{n}a_{j_{i}}^{-}(t_{i})\right)\left(\prod_{i=n}^{1}a_{\bar{j}_{i}}^{+}(\bar{t}_{i})\right)\right)\right>.$
(136)

Without loss of generality, we may assume that $t_{1}\geqslant\cdots\geqslant
t_{n}$ and $\bar{t}_{1}\geqslant\cdots\geqslant\bar{t}_{n}$ (indeed, (135) is
antisymmetric under exchanges of $a^{-}$’s and exchanges of $a^{+}$’s).

1 -  We diagonalize $\mu$ by a unitary transform $U$, and define

$b_{i}:=\sum_{j=1}^{N}U_{j,i}^{*}a_{j},\quad
b_{i}^{\dagger}:=\sum_{j=1}^{N}U_{j,i}a_{j}^{\dagger}$ (137)

in terms of which

$e^{-t\mathcal{H}_{0}}=e^{-t\sum_{i=1}^{N}\lambda_{i}b_{i}^{\dagger}b_{i}}=\prod_{i=1}^{N}(1+(e^{-t\lambda_{i}}-1)b_{i}^{\dagger}b_{i})$
(138)

where $\lambda_{i}$ are the eigenvalues of $\mu$. Note that

$\\{b_{i},b_{j}^{\dagger}\\}=\sum_{k,l=1}^{N}U_{k,i}^{*}U_{l,j}\\{a_{k},a_{l}^{\dagger}\\}=\sum_{k=1}^{N}U_{k,i}^{*}U_{k,j}=\delta_{i,j}$
(139)

Let

$b_{i}^{-}(s):=e^{t\mathcal{H}_{0}}b_{i}e^{-t\mathcal{H}_{0}}=\sum_{j=1}^{N}U_{j,i}^{*}a_{j}^{-}(s),\quad
b_{i}^{+}(s):=e^{t\mathcal{H}_{0}}b_{i}^{\dagger}e^{-t\mathcal{H}_{0}}=\sum_{j=1}^{N}U_{j,i}a_{j}^{+}(s).$
(140)

By (103),

$b_{i}^{-}(s)=e^{-s\lambda_{i}}b_{i},\quad
b_{i}^{+}(s)=e^{s\lambda_{i}}b_{i}^{\dagger}.$ (141)

We have

$\begin{array}[]{@{}>{\displaystyle}l@{}}\hphantom{\hskip
423.62pt}\\\\[-14.22636pt]
\left<\mathbf{T}\left(\left(\prod_{i=1}^{n}a_{j_{i}}^{-}(t_{i})\right)\left(\prod_{i=n}^{1}a_{\bar{j}_{i}}^{+}(\bar{t}_{i})\right)\right)\right>=\\\
\hfill=\sum_{k_{1},\cdots,k_{n}=1}^{N}\sum_{l_{1},\cdots,l_{n}=1}^{N}\left(\prod_{i=1}^{n}U_{j_{i},k_{i}}U_{\bar{j}_{i},l_{i}}^{*}\right)\left<\mathbf{T}\left(\left(\prod_{i=1}^{n}b_{k_{i}}^{-}(t_{i})\right)\left(\prod_{i=n}^{1}b_{l_{i}}^{+}(\bar{t}_{i})\right)\right)\right>.\end{array}$
(142)

2 -  Let $\sigma\in\\{1,-1\\}$ and
$\gamma_{1},\cdots,\gamma_{2n}\in\\{1,\cdots,N\\}$,
$\omega_{1},\cdots,\omega_{2n}\in\\{+,-\\}$, $\beta\geqslant
s_{1}\geqslant\cdots\geqslant s_{2n}\geqslant 0$ be such that

$\mathbf{T}\left(\left(\prod_{i=1}^{n}b_{k_{i}}^{-}(t_{i})\right)\left(\prod_{i=n}^{1}b_{l_{i}}^{+}(\bar{t}_{j})\right)\right)=\sigma\prod_{i=1}^{2n}b_{\gamma_{i}}^{\omega_{i}}(s_{i})$
(143)

and, by (103),

$\mathbf{T}\left(\left(\prod_{i=1}^{n}b_{k_{i}}^{-}(t_{i})\right)\left(\prod_{i=n}^{1}b_{l_{i}}^{+}(\bar{t}_{j})\right)\right)=\left(\prod_{i=1}^{n}e^{-\lambda_{k_{i}}t_{i}+\lambda_{l_{i}}\bar{t}_{i}}\right)\sigma\prod_{i=1}^{2n}b_{\gamma_{i}}^{\omega_{i}}$
(144)

where $b_{\gamma}^{+}\equiv b_{\gamma}^{\dagger}$ and $b_{\gamma}^{-}\equiv
b_{\gamma}$. We now commute the $b_{\gamma}$ back into the order
$\prod_{i=1}^{n}b_{k_{i}}\prod_{i=n}^{1}b_{l_{i}}^{\dagger}$, in such a way
that the $\sigma$ sign is canceled. This is not entirely straightforward:
whenever a $b^{\dagger}_{l_{i}}$ is commuted with a $b_{k_{j}}$, it either
moves through or yields $\delta_{l_{i},k_{j}}$. To keep track of these terms,
we introduce the following notation: given a product of $b$’s, we define

$\left|b^{\dagger}_{l_{n}}\cdots
b^{\dagger}_{l_{1}}\right\\}:=\left(\prod_{i=1}^{n}b_{k_{i}}\right)\left(\prod_{i=n}^{1}b^{\dagger}_{l_{i}}\right)$
(145)

and make the symbol $\left|\cdot\right\\}$ multilinear: for instance,

$\left|(b^{\dagger}_{l_{n}}-1)\cdots
b^{\dagger}_{l_{1}}\right\\}\equiv\left(\prod_{i=1}^{n}b_{k_{i}}\right)\left(\prod_{i=n}^{1}b^{\dagger}_{l_{i}}\right)-\left(\prod_{i=1}^{n-1}b_{k_{i}}\right)\left(\prod_{i=n-1}^{1}b^{\dagger}_{l_{i}}\right).$
(146)

We move the $b^{\dagger}$ to the right of the $b$ and put them back in the
order $b_{l_{n}}^{\dagger}\cdots b_{l_{1}}^{\dagger}$. Note that, since
$\bar{t}_{1}\geqslant\cdots\bar{t}_{n}$ and $t_{1}\geqslant\cdots\geqslant
t_{n}$, the $b_{i}$ are already ordered as in $\prod_{i=1}^{n}b_{i}$, but the
$b_{i}^{\dagger}$ are ordered in the opposite order as in
$\prod_{i=n}^{1}b_{i}^{\dagger}$. As we will now discuss, reordering the $b$’s
in this way yields

$\sigma\prod_{i=1}^{2n}b_{\gamma_{i}}^{\omega_{i}}=\left|\prod_{i=n}^{1}\left(b_{l_{i}}^{\dagger}-{\textstyle\sum_{j=1}^{n}}(-1)^{j+i}\delta_{l_{i},k_{j}}\mathds{1}_{\bar{t}_{i}>t_{j}}\right)\right\\}$
(147)

where $\mathds{1}_{\bar{t}_{i}>t_{j}}\in\\{0,1\\}$ is equal to 1 if and only
if $\bar{t}_{i}>t_{j}$. When moving the $b^{\dagger}$’s, the operators either
anti-commute, or, when passing $b_{i}^{\dagger}$ over $b_{i}$, the two
operators may destroy each other (the $\delta_{i,j}$ term in
$\\{a_{i},a_{j}^{\dagger}\\}=\delta_{i,j}$). When the operators destroy each
other, they become numbers, and commute with all other operators. After this,
they no longer produce the signs necessary to cancel out $\sigma$. The sign in
$\sigma$ due to those terms is $-1$ to the power of the number of positions
between the location where the operators were destroyed and the location where
they end up. Because of the ordering of the $t$’s and $\bar{t}$’s, if
$b^{\dagger}_{l_{i}}$ and $b_{k_{j}}$ destroyed each other, then there are
$n-j$ $b$’s left to go over and $n-i$ $b^{\dagger}$’s. Thus, the contribution
of $\sigma$ for this term is $(-1)^{n-j+n-i}=(-1)^{i+j}$. This justifies
(147).

3 -  We show that

$\left<\left|\prod_{i=n}^{1}b_{l_{i}}^{\dagger}\right\\}\right>\equiv\left<\left(\prod_{i=1}^{n}b_{k_{i}}\right)\left(\prod_{i=n}^{1}b_{l_{i}}^{\dagger}\right)\right>=\sum_{\tau\in\mathcal{S}_{n}}(-1)^{\tau}\prod_{i=1}^{n}\left<b_{k_{i}}b_{l_{\tau(i)}}^{\dagger}\right>.$
(148)

First of all, $\prod_{i=n}^{1}b_{l_{i}}^{\dagger}$ is non-zero only if all
$l_{i}$ are different from each other, and so is the right side of (148), so
we can assume that all $l_{i}$ are different from each other. Next, since
$\mathcal{H}_{0}=\sum_{\lambda_{i}}b_{i}^{\dagger}b_{i}$, the
$b_{l_{i}}^{\dagger}$ and $b_{k_{i}}$ must be paired up. Using the
antisymmetry of both sides of (148), we can assume without loss of generality
that $k_{i}=l_{i}$. In that case,

$\mathrm{Tr}\left(e^{-\beta\mathcal{H}_{0}}\left(\prod_{i=1}^{n}b_{k_{i}}\right)\left(\prod_{i=n}^{1}b_{l_{i}}^{\dagger}\right)\right)=\mathrm{Tr}\left(e^{-\beta\mathcal{H}_{0}}\left(\prod_{i=1}^{n}b_{k_{i}}\right)\left(\prod_{i=n}^{1}b_{k_{i}}^{\dagger}\right)\right)=\mathrm{Tr}\left(e^{-\beta\mathcal{H}_{0}}\prod_{i=1}^{n}b_{k_{i}}b_{k_{i}}^{\dagger}\right)$
(149)

and so

$\left<e^{-\beta\mathcal{H}_{0}}\left(\prod_{i=1}^{n}b_{k_{i}}\right)\left(\prod_{i=n}^{1}b_{l_{i}}^{\dagger}\right)\right>=\prod_{i=1}^{n}\left<b_{k_{i}}b_{k_{i}}^{\dagger}\right>$
(150)

which is equal to the right side of (148).

4 -  By (148),

$\left<\left|\prod_{i=n}^{1}\left(b_{l_{i}}^{\dagger}-{\textstyle\sum_{j=1}^{n}}(-1)^{i+j}\delta_{l_{i},k_{j}}\mathds{1}_{\bar{t}_{i}>t_{j}}\right)\right\\}\right>=\sum_{\tau\in\mathcal{S}_{n}}(-1)^{\tau}\prod_{i=1}^{n}\left(\left<b_{k_{i}}b_{l_{\tau(i)}}^{\dagger}\right>-\delta_{k_{i},l_{\tau(i)}}\mathds{1}_{\bar{t}_{\tau(i)}>t_{i}}\right).$
(151)

Now,

$\left<b_{k_{i}}b_{l_{\tau(i)}}^{\dagger}\right>=\delta_{k_{i},l_{\tau(i)}}\frac{1}{1+e^{-\beta\lambda_{k_{i}}}}$
(152)

and

$\left<b_{k_{i}}b_{l_{\tau(i)}}^{\dagger}\right>-\delta_{k_{i},l_{\tau(i)}}=-\delta_{k_{i},l_{\tau(i)}}\frac{e^{-\beta\lambda_{k_{i}}}}{1+e^{-\beta\lambda_{k_{i}}}}.$
(153)

Thus, inserting (152), (153) into (151), (147), (144), (142) and (136), we
find

$\begin{array}[]{@{}>{\displaystyle}l@{}}\hphantom{\hskip
423.62pt}\\\\[-14.22636pt]
\left<\mathbf{T}\left(\prod_{i=1}^{n}a_{j_{i}}^{-}(t_{i})a_{\bar{j}_{i}}^{+}(\bar{t}_{i})\right)\right>=\\\
\hfill=\sum_{k_{1},\cdots,k_{n}=1}^{N}\sum_{l_{1},\cdots,l_{n}=1}^{N}\left(\prod_{i=1}^{n}U_{j_{i},k_{i}}U_{\bar{j}_{i},l_{i}}^{*}\right)\sum_{\tau\in\mathcal{S}_{n}}(-1)^{\tau}\prod_{i=1}^{n}\delta_{k_{i},l_{\tau(i)}}e^{-\lambda_{k_{i}}(t_{i}-\bar{t}_{\tau(i)})}\frac{-\mathds{1}_{\bar{t}_{\tau(i)}>t_{i}}e^{-\beta\lambda_{k_{i}}}+\mathds{1}_{\bar{t}_{\tau(i)}<t_{i}}}{1-e^{-\beta\lambda_{k_{i}}}}\end{array}$
(154)

and, relabeling the $l_{i}$,

$\begin{array}[]{@{}>{\displaystyle}l@{}}\hphantom{\hskip
423.62pt}\\\\[-14.22636pt]
\left<\mathbf{T}\left(\prod_{i=1}^{n}a_{j_{i}}^{-}(t_{i})a_{\bar{j}_{i}}^{+}(\bar{t}_{i})\right)\right>=\\\
\hfill=\sum_{\tau\in\mathcal{S}_{n}}(-1)^{\tau}\prod_{i=1}^{n}\left(\sum_{k=1}^{N}U_{j_{i},k}U_{\bar{j}_{\tau(i)},k}^{*}e^{-\lambda_{k}(t_{i}-\bar{t}_{\tau(i)})}\frac{-\mathds{1}_{\bar{t}_{\tau(i)}>t_{i}}e^{-\beta\lambda_{k}}+\mathds{1}_{\bar{t}_{\tau(i)}<t_{i}}}{1-e^{-\beta\lambda_{k}}}\right)\end{array}$
(155)

which, by (119) and (122), is

$\left<\mathbf{T}\left(\prod_{i=1}^{n}a_{j_{i}}^{-}(t_{i})a_{\bar{j}_{i}}^{+}(\bar{t}_{i})\right)\right>\sum_{\tau\in\mathcal{S}_{n}}(-1)^{\tau}\prod_{i=1}^{n}\left<\mathbf{T}(a_{j_{i}}^{-}(t_{i})a_{\bar{j}_{\tau(i)}}^{+}(\bar{t}_{\tau(i)}))\right>.$
(156)

$\square$

##  Appendix B Properties of Gaussian Grassman integrals

Consider a family $\psi_{1}^{\pm},\cdots,\psi_{N}^{\pm}$ of Grassmann
variables and the Gaussian Grassmann measure

$P_{\mu}(d\psi)=\det(\mu)\prod_{i=1}^{N}d\psi_{i}^{+}d\psi_{i}^{-}e^{-\sum_{i,j=1}^{N}\mu_{i,j}\psi_{i}^{+}\mu^{-1}_{i,j}\psi_{i}^{-}}$
(157)

where $\mu$ is a Hermitian, invertible matrix.

We first prove that one can perform changes of variables in Grassmann
integrals.

Lemma B.4

(Change of variables in Grassmann integrals)

Given a unitary $N\times N$ matrix $U$ and a polynomial $f$ in the Grassmann
variables $\psi$, if

$\varphi_{i}^{-}:=\sum_{j=1}^{N}U_{i,j}\psi_{j}^{-},\quad\varphi_{i}^{+}:=\sum_{j=1}^{N}U_{i,j}^{*}\psi_{j}^{+}$
(158)

then

$\int\prod_{i=1}^{N}d\varphi_{i}^{+}d\varphi_{i}^{-}\
f(U^{\dagger}\varphi^{-},(U^{\dagger})^{*}\varphi^{+})=\int\prod_{i=1}^{N}d\psi_{i}^{+}d\psi_{i}^{-}\
f(\psi^{-},\psi^{+}).$ (159)

Proof: Without loss of generality, we can assume that the highest order term
of $f$ is

$\alpha\prod_{i=1}^{N}\psi_{i}^{-}\psi_{i}^{+}.$ (160)

We have

$\int\prod_{i=1}^{N}d\psi_{i}^{+}d\psi_{i}^{-}\ f(\psi^{-},\psi^{+})=\alpha.$
(161)

Furthermore,

$\int\prod_{i=1}^{N}d\varphi_{i}^{+}d\varphi_{i}^{-}\
f(U^{\dagger}\varphi^{-},U\varphi^{+})=\alpha\int\prod_{i=1}^{N}d\varphi_{i}^{+}d\varphi_{i}^{-}\
\prod_{i=1}^{N}\left(\sum_{j,k=1}^{N}U_{i,j}^{*}\varphi_{j}^{-}U_{i,k}\varphi_{k}^{+}\right)$
(162)

and, avoiding repetitions in the products, we find that

$\prod_{i=1}^{N}\left(\sum_{j,k=1}^{N}U_{i,j}^{*}\varphi_{j}^{-}U_{i,k}\varphi_{k}^{+}\right)=\sum_{\tau,\tau^{\prime}\in\mathcal{S}_{N}}\prod_{i=1}^{N}U^{*}_{i,\tau(i)}U_{i,\tau^{\prime}(i)}\varphi_{\tau(i)}^{-}\varphi_{\tau^{\prime}(i)}^{+}$
(163)

and so, reordering the $\varphi$,

$\prod_{i=1}^{N}\left(\sum_{j,k=1}^{N}U_{i,j}^{*}\varphi_{j}^{-}U_{i,k}\varphi_{k}^{+}\right)=\sum_{\tau,\tau^{\prime}\in\mathcal{S}_{N}}(-1)^{\tau}(-1)^{\tau^{\prime}}\prod_{i=1}^{N}U^{*}_{i,\tau(i)},U_{i,\tau^{\prime}(i)}\varphi_{i}^{-}\varphi_{i}^{+}=|\det(U)|^{2}\prod_{i=1}^{N}\varphi_{i}^{-}\varphi_{i}^{+}$
(164)

and, since $\det(U)=1$,

$\int\prod_{i=1}^{N}d\varphi_{i}^{+}d\varphi_{i}^{-}\
f(U^{\dagger}\varphi^{-},U\varphi^{+})=\alpha=\int\prod_{i=1}^{N}d\psi_{i}^{+}d\psi_{i}^{-}\
f(\psi^{-},\psi^{+}).$ (165)

$\square$

###  B.1 Gaussian Grassmann integrals

Let us now turn to a set of identities for Gaussian Grassmann integrals.

Lemma B.5

We have

$\int P_{\mu}(d\psi)\ 1=1$ (166)

and

$\int P_{\mu}(d\psi)\ \psi_{i}^{-}\psi_{j}^{+}=\mu_{i,j}$ (167)

where $\mu^{-1}$ is the inverse matrix of $\mu$.

Proof: This is a consequence of the change of variables lemma B. We
diagonalize $\mu^{-1}$ by a unitary transform $U$, and change variables to
$\varphi^{-}_{i}=\sum_{j}U_{j,i}^{*}\psi^{-}_{j}$ and
$\varphi^{+}_{i}=\sum_{j}U_{j,i}\psi^{+}_{j}$:

$\int P_{\mu}(d\psi)\
1=\det(\mu)\int\prod_{i=1}^{N}d\varphi_{i}^{+}d\varphi_{i}^{-}\
e^{-\sum_{i=1}^{N}\lambda_{i}\varphi_{i}^{+}\varphi_{i}^{-}}$ (168)

where $\lambda_{i}$ are the eigenvalues of $\mu$. Furthermore,

$e^{-\sum_{i=1}^{N}\lambda_{i}\varphi_{i}^{+}\varphi_{i}^{-}}=\prod_{i=1}^{N}e^{-\lambda_{i}\varphi_{i}^{+}\varphi_{i}^{-}}=\prod_{i=1}^{N}(1-\lambda_{i}\varphi_{i}^{+}\varphi_{i}^{-})$
(169)

so

$\int P_{\mu}(d\psi)\
1=\det(\mu)\int\prod_{i=1}^{N}d\varphi_{i}^{+}d\varphi_{i}^{-}\
\prod_{i=1}^{N}(1+\lambda_{i}\varphi_{i}^{-}\varphi_{i}^{+})=\det(\mu)\prod_{i=1}^{N}\lambda_{i}=1.$
(170)

Similarly,

$\int P_{\mu}(d\psi)\
\psi_{i}^{-}\psi_{j}^{+}=\sum_{k,l=1}^{N}\det(\mu)\int\prod_{i=1}^{N}d\varphi_{i}^{+}d\varphi_{i}^{-}\
e^{-\sum_{i=1}^{N}\lambda_{i}\varphi_{i}^{+}\varphi_{i}^{-}}U_{i,k}\varphi_{k}^{-}U_{j,l}^{*}\varphi_{l}^{+}$
(171)

in which we insert (169) to find that $k$ must be equal to $l$ and

$\int P_{\mu}(d\psi)\
1=\det(\mu)\sum_{k=1}^{N}U_{k,j}^{*}U_{i,k}\lambda_{k}^{-1}\prod_{i=1}^{N}\lambda_{i}=\sum_{k=1}^{N}U_{j,k}^{*}U_{i,k}\lambda_{k}^{-1}=\mu_{i,j}.$
(172)

$\square$

###  B.2 Wick rule for Gaussian Grassmann integrals

Let us now turn to the Wick rule for Gaussian Grassmann integrals.

Lemma B.6

(Wick rule for Gaussian Grassmann integrals)

For $n\in\\{1,\cdots,N\\}$, $j_{1},\cdots,j_{n}\in\\{1,\cdots,N\\}$ and
$\bar{j}_{1},\cdots,\bar{j}_{n}\in\\{1,\cdots,N\\}$,

$\int P_{\mu}(d\psi)\
\prod_{i=1}^{n}\psi_{j_{i}}^{-}\psi_{\bar{j}_{i}}^{+}=\sum_{\tau\in\mathcal{S}_{n}}(-1)^{\tau}\prod_{i=1}^{n}\int
P_{\mu}(d\psi)\ \psi_{j_{i}}^{-}\psi_{\bar{j}_{\tau(i)}}^{+}.$ (173)

Proof: As before, we change variables using lemma B. We diagonalize $\mu^{-1}$
by a unitary transform $U$, and change variables to
$\varphi^{-}_{i}=\sum_{j}U_{j,i}^{*}\psi^{-}_{j}$ and
$\varphi^{+}_{i}=\sum_{j}U_{j,i}\psi^{+}_{j}$:

$\begin{array}[]{@{}>{\displaystyle}l@{}}\hphantom{\hskip
423.62pt}\\\\[-14.22636pt] \int P_{\mu}(d\psi)\
\prod_{i=1}^{n}\psi_{j_{i}}^{-}\psi_{\bar{j}_{i}}^{+}=\\\
\hfill=\det(\mu)\sum_{k_{1},\cdots,k_{n}=1}^{N}\sum_{l_{1},\cdots,l_{n}=1}^{N}\int\prod_{i=1}^{N}d\varphi_{i}^{+}d\varphi_{i}^{-}\
\left(\prod_{i=1}^{n}(1+\lambda_{i}\varphi_{i}^{-}\varphi_{i}^{+})\right)\prod_{i=1}^{n}U_{j_{i},k_{i}}\varphi_{k_{i}}^{-}U_{\bar{j}_{i},l_{i}}^{*}\varphi_{l_{i}}^{+}.\end{array}$
(174)

Now, $\prod_{i}\varphi_{k_{i}}^{-}\varphi_{l_{i}}^{+}$ must be permutable into
products of pairs with identical indices, otherwise the resulting integral
would be 0. Summing over all such possibilities, we find

$\begin{array}[]{@{}>{\displaystyle}l@{}}\hphantom{\hskip
423.62pt}\\\\[-14.22636pt] \int P_{\mu}(d\psi)\
\prod_{i=1}^{n}\psi_{j_{i}}^{-}\psi_{\bar{j}_{i}}^{+}=\det(\mu)\cdot\\\
\hfill\cdot\sum_{k_{1},\cdots,k_{n}=1}^{N}\sum_{\tau\in\mathcal{S}_{n}}(-1)^{\tau}\left(\prod_{i=1}^{n}U_{j_{i},k_{i}}U_{\bar{j}_{\tau(i)},k_{i}}^{*}\right)\int\prod_{i=1}^{N}d\varphi_{i}^{+}d\varphi_{i}^{-}\
\left(\prod_{i=1}^{n}(1+\lambda_{i}\varphi_{i}^{-}\varphi_{i}^{+})\right)\prod_{i=1}^{n}\varphi_{k_{i}}^{-}\varphi_{k_{i}}^{+}\end{array}$
(175)

and so

$\int P_{\mu}(d\psi)\
\prod_{i=1}^{n}\psi_{j_{i}}^{-}\psi_{\bar{j}_{i}}^{+}=\sum_{\tau\in\mathcal{S}_{n}}(-1)^{\tau}\prod_{i=1}^{n}\left(\sum_{k=1}^{N}U_{j_{i},k}U_{\bar{j}_{\tau(i)},k}^{*}\lambda_{k}^{-1}\right)=\sum_{\tau\in\mathcal{S}_{n}}(-1)^{\tau}\prod_{i=1}^{n}\mu_{j_{i},j_{\tau(i)}}.$
(176)

$\square$

###  B.3 Addition property

Finally, let us prove the addition property of Gaussiann Grassmann integrals.

Lemma B.7

(Addition property)

Consider two invertible matrices $\mu_{1}$ and $\mu_{2}$ and its two
associated Gaussian Grassmann measure $P_{\mu_{1}}(d\psi_{1})$ and
$P_{\mu_{2}}(d\psi_{2})$. For any polynomial $f$,

$\int P_{\mu_{1}+\mu_{2}}(d\psi)\ f(\psi)=\int P_{\mu_{1}}(d\psi_{1})\int
P_{\mu_{2}}(d\psi_{2})\ f(\psi_{1}+\psi_{2}).$ (177)

Proof: It is sufficient to prove the lemma when $f$ is a monomial of the form

$f(\psi)=\prod_{i=1}^{n}\psi_{j_{i}}^{-}\psi_{\bar{j}_{i}}^{+}.$ (178)

First, note that $P_{\mu_{1}}(d\psi_{1})P_{\mu_{2}}(d\psi_{2})$ is
proportional to a Gaussian Grassmann measure. Next, change variables to
$\varphi_{1}=\frac{1}{\sqrt{2}}(\psi_{1}+\psi_{2})$ and
$\varphi_{2}=\frac{1}{\sqrt{2}}(\psi_{1}-\psi_{2})$ in terms of which

$\int P_{\mu_{1}}(d\psi_{1})\int P_{\mu_{2}}(d\psi_{2})\
f(\psi_{1}+\psi_{2})=\int P_{\nu_{1}}(d\varphi_{1})\int
P_{\nu_{2}}(d\varphi_{2})\ f(\varphi_{1})$ (179)

where $\nu_{1}$ and $\nu_{2}$ can be computed from the change of variables,
but this is not necessary here. Since $f$ is a monomial, it can be computed
using the Wick rule B.2, and thus, changing variables back to $\psi$,

$\begin{array}[]{@{}>{\displaystyle}l@{}}\hphantom{\hskip
423.62pt}\\\\[-14.22636pt] \int P_{\mu_{1}}(d\psi_{1})\int
P_{\mu_{2}}(d\psi_{2})\ f(\psi)=\\\
\hfill=\sum_{\tau\in\mathcal{S}_{n}}(-1)^{\tau}\prod_{i=1}^{n}\int
P_{\mu_{1}}(d\psi_{1})\int P_{\mu_{2}}(d\psi_{2})\
(\psi_{1,j_{i}}^{-}+\psi_{2,j_{i}}^{-})(\psi_{1,\bar{j}_{\tau(i)}}^{+}+\psi_{2,\bar{j}_{\tau(i)}}^{+}).\end{array}$
(180)

Now,

$\begin{array}[]{@{}>{\displaystyle}l@{}}\hphantom{\hskip
423.62pt}\\\\[-14.22636pt] \int P_{\mu_{1}}(d\psi_{1})\int
P_{\mu_{2}}(d\psi_{2})\
(\psi_{1,j_{i}}^{-}+\psi_{2,j_{i}}^{-})(\psi_{1,\bar{j}_{\tau(i)}}^{+}+\psi_{2,\bar{j}_{\tau(i)}}^{+})=\\\
\hfill=\int P_{\mu_{1}}(d\psi_{1})\
\psi_{1,j_{i}}^{-}\psi_{1,\bar{j}_{\tau(i)}}^{+}+\int P_{\mu_{2}}(d\psi_{2})\
\psi_{2,j_{i}}^{-}\psi_{2,\bar{j}_{\tau(i)}}^{+}\end{array}$ (181)

so, by lemma B.1,

$\int P_{\mu_{1}}(d\psi_{1})\int P_{\mu_{2}}(d\psi_{2})\
(\psi_{1,j_{i}}^{-}+\psi_{2,j_{i}}^{-})(\psi_{1,\bar{j}_{\tau(i)}}^{+}+\psi_{2,\bar{j}_{\tau(i)}}^{+})=(\mu_{1}+\mu_{2})_{j_{i},\bar{j}_{\tau(i)}}$
(182)

and

$\int P_{\mu_{1}+\mu_{2}}(d\psi)\
\psi_{j_{i}}^{-}\psi_{\bar{j}_{\tau(i)}}^{+}=(\mu_{1}+\mu_{2})_{j_{i},\bar{j}_{\tau(i)}}.$
(183)

We conclude the proof of the lemma using the Wick rule for $\psi$:

$\int P_{\mu_{1}+\mu_{2}}(d\psi)\
f(\psi)=\sum_{\tau\in\mathcal{S}_{n}}(-1)^{\tau}\prod_{i=1}^{n}\int
P_{\mu_{1}+\mu_{2}}(d\psi)\ \psi_{j_{i}}^{-}\psi_{\bar{j}_{\tau(i)}}^{+}.$
(184)

$\square$

## References

  * [An70] P.W. Anderson - A poor man’s derivation of scaling laws for the Kondo problem, Journal of Physics C: Solid State Physics, volume 3, page 2436, 1970,doi:10.1088/0022-3719/3/12/008.
  * [An80] N. Andrei - Diagonalization of the Kondo Hamiltonian, Physical Review Letters, volume 45, issue 5, 1980,doi:10.1103/PhysRevLett.45.379.
  * [BF84] G.A. Battle, P. Federbush - A note on cluster expansions, tree graph identities, extra $1/N!$ factors!!!, Letters in Mathematical Physics, volume 8, pages 55-57, 1984,doi:10.1007/BF00420041.
  * [BCe78] G. Benfatto, M. Cassandro, G. Gallavotti, F. Nicolò, E. Olivieri, E. Presutti, E. Scacciatelli - Some probabilistic techniques in field theory, Communications in Mathematical Physics, volume 59, issue 2, pages 143-166, 1978,doi:10.1007/BF01614247.
  * [BG90] G. Benfatto, G. Gallavotti - Perturbation theory of the Fermi surface in a quantum liquid - a general quasiparticle formalism and one-dimensional systems, Journal of Statistical Physics, volume 59, issue 3-4, pages 541-664, 1990,doi:10.1007/BF01025844.
  * [BGe94] G. Benfatto, G. Gallavotti, A.Procacci, B. Scoppola - Beta function and Schwinger functions for a many Fermions system in one dimension - Anomaly of the Fermi surface, Communications in Mathematical Physics, volume 160, pages 93-171, 1994,doi:10.1007/BF02099791.
  * [BGJ15] G. Benfatto, G. Gallavotti, I. Jauslin - Kondo Effect in a Fermionic Hierarchical Model, Journal of Statistical Physics, volume 161, issue 5, pages 1203-1230, 2015,doi:10.1007/s10955-015-1378-7, arxiv:1506.04381.
  * [BF78] D. Brydges, P. Federbush - A new form of the Mayer expansion in classical statistical mechanics, Journal of Mathematical Physics, volume 19, page 2064, 1978,doi:10.1063/1.523586.
  * [Do91] T.C. Dorlas - Renormalization group analysis of a simple hierarchical fermion model, Communications in Mathematical Physics, volume 136, pages 169-194, 1991,doi:10.1007/BF02096796.
  * [Dy69] F.J. Dyson - Existence of a phase-transition in a one-dimensional Ising ferromagnet, Communications in Mathematical Physics, volume 12, pages 91-107, 1969,doi:10.1007/BF01645907.
  * [GJ15] G. Gallavotti, I. Jauslin - Kondo Effect in the Hierarchical s-d Model, Journal of Statistical Physics, volume 161, issue 5, pages 1231-1235, 2015,doi:10.1007/s10955-015-1370-2, arxiv:1507.05678.
  * [GN85] G. Gallavotti, F. Nicolò - Renormalization theory for four dimensional scalar fields I, Communications in Mathematical Physics, volume 100, pages 545-590, 1985,doi:10.1007/BF01217729.
  * [GK81] K. Gawedzki, A. Kupiainen - Renormalization group study of a critical lattice model I: Convergence to the line of fixed points, Communications in Mathematical Physics, volume 82, issue 3, pages 407-433, 1981,doi:10.1007/BF01237048.
  * [Gi10] A. Giuliani - The Ground State Construction of the Two-dimensional Hubbard Model on the Honeycomb Lattice, Quantum Theory from Small to Large Scales, lecture notes of the Les Houches Summer School, volume 95, Oxford University Press, 2010, arxiv:1102.3881.
  * [GGM12] A. Giuliani, R.L. Greenblatt, V. Mastropietro - The scaling limit of the energy correlations in non-integrable Ising models, Journal of Mathematical Physics, volume 53, page 095214, 2012,doi:10.1063/1.4745910, arxiv:1204.4040.
  * [GJ16] A. Giuliani, I. Jauslin - The ground state construction of bilayer graphene, Reviews in Mathematical Physics, volume 28, issue 8, page 1650018, 2016,doi:10.1142/S0129055X16500185, arxiv:1507.06024.
  * [GM10] A. Giuliani, V. Mastropietro - The two-dimensional Hubbard model on the honeycomb lattice, Communications in Mathematical Physics, volume 293, pages 301-364, 2010,doi:10.1007/s00220-009-0910-5, arxiv:0811.1881.
  * [GMP12] A. Giuliani, V. Mastropietro, M. Porta - Universality of conductivity in interacting graphene, Communications in Mathematical Physics, volume 311, issue 2, pages 317-355, 2012,doi:10.1007/s00220-012-1444-9, arxiv:1101.2169.
  * [GMP17] A. Giuliani, V. Mastropietro, M. Porta - Universality of the Hall Conductivity in Interacting Electron Systems, Communications in Mathematical Physics, volume 349, issue 3, pages 1107-1161, 2017,doi:10.1007/s00220-016-2714-8, arxiv:1511.04047.
  * [Ko64] J. Kondo - Resistance minimum in dilute magnetic alloys, Progress of Theoretical Physics, volume 32, issue 1, 1964,doi:10.1143/PTP.32.37.
  * [Ma11] V. Mastropietro - Conductivity between Luttinger liquids: coupled chains and bilayer graphene, Physical Review B, volume 84, page 035109, 2011,doi:10.1103/PhysRevB.84.035109, arxiv:1012.5736.
  * [Mc57] J.W. McClure - Band structure of graphite and de Haas-van Alphen effect, Physical review, volume 108, pages 612-618, 1957,doi:10.1103/PhysRev.108.612.
  * [NGe04] K.S. Novoselov, A.K. Geim, S. V.Morozov, D. Jiang, Y. Zhang, S.V. Dubonos, I.V. Grigorieva, A.A. Firsov - Electric field effect in atomically thin carbon films, Science, vol. 306, pages 666-669, 2004,doi:10.1126/science.1102896.
  * [SW58] J.C. Slonczewski, P.R. Weiss - Band structure of graphite, Physical Review, volume 109, pages 272-279, 1958,doi:10.1103/PhysRev.109.272.
  * [Wi65] K.G. Wilson - Model Hamiltonians for Local Quantum Field Theory, Physical Review, volume 140, issue 2B, page B445-B457, 1965,doi:10.1103/PhysRev.140.B445.
  * [Wi75] K.G. Wilson - The renormalization group: Critical phenomena and the Kondo problem, Reviews of Modern Physics, volume 47, issue 4, 1975,doi:10.1103/RevModPhys.47.773.
  * [mk] meankondo software package v1.5, I. Jauslin,http://ian.jauslin.org/software/meankondo.