These are the lecture notes for the summer school {\rm``Quantum Mechanics from Condensed Matter to Computing''}, organized by {\rm Niels Benedikter, Marcin Napi\'orkowski, Jan Philip Solovej} and {\rm Albert Werner} in Copenhagen from June 13 to 17, 2022.
}
\vfill
\hfil{\bf Abstract}\par
{\small
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.
}
\vfill
\tableofcontents
\vfill
\eject
\setcounter{page}1
\pagestyle{plain}
\section{Introduction}
\indent
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 {\it Wilsonian} renormalization group\-~\cite{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\-~\cite{BG90}, which has been applied successfully to a variety of systems: to name but a few, \cite{BGe94,Ma11,GMP12,GGM12,GMP17} among many others.
\bigskip
\indent
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\-~\cite{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.
\bigskip
\indent
Graphene is often studied in the approximation that its electrons do not interact with each other\-~\cite{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, {\it 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.
\bigskip
\indent
The renormalization group analysis of interacting graphene was carried out in\-~\cite{GM10,Gi10}, and was used to prove the universality of the conductivity\-~\cite{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 {\it hierarchical graphene model}.
Hierarchical models have long been used:\-~\cite{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\-~\cite{Do91}, are even simpler than Bosonic ones: as was shown in\-~\cite{BGJ15,GJ15}, they are {\it integrable}, in that their renormalization group analysis can be carried out {\it explicitly} (in Bosonic and non-hierarchical cases, the renormalization group produces infinite power series).
It is important to emphasize that hierarchical models are {\it 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.
\bigskip
\indent
For this paper, we will define the hierarchical graphene model, following the ideas of\-~\cite{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 {\it exactly} and {\it 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\-~\cite{Gi10}.
\bigskip
\indent
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 {\it 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 {\it 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 {\it 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 {\it any} other strongly coupled system.
\bigskip
\indent
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\-~\cite{Ko64} as a toy model to study conductance in disordered systems.
It was studied rather extensively by Anderson\-~\cite{An70}, and later found to be exactly solvable by Andrei\-~\cite{An80}.
In developing his version of the renormalization group, Wilson studied the Kondo model\-~\cite{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 {\it 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.
\bigskip
\indent
The Kondo effect can be seen in the renormalization group flow, which, at small energy scales, goes to a {\it 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\-~\cite{An80}, it has never been shown using the renormalization group.
It is, however, provable in the hierarchical Kondo model\-~\cite{BGJ15,GJ15}, which will be briefly discussed at the end of this paper.
\bigskip
\indent
The rest of this paper is structured as follows.
In section\-~\ref{sec:graphene}, 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\-~\ref{sec:hierarchical_graphene}, we define the hierarchical graphene model and study its renormalization group flow.
In section\-~\ref{sec:hierarchical_kondo}, we define the hierarchical Kondo model and discuss its renormalization group flow.
There are two appendices in which useful lemmas are proved: in appendix\-~\ref{app:free_fermions} we compute properties of free Fermion systems, and in appendix\-~\ref{app:grassmann}, we prove some properties of Gaussian Grassmann integrals.
\section{Graphene}\label{sec:graphene}
\indent
In this section, we describe a model for the electrons in graphene.
We will use the {\it tight-binding} approximation, in which the electrons are assumed to be bound to carbon atoms on a hexagonal lattice (see figure\-~\ref{fig:graphene}).
For more details on this model, see\-~\cite{Gi10,GM10}.
\bigskip
\begin{figure}
\hfil\includegraphics[width=12cm]{graphene.pdf}
\caption{
The hexagonal lattice.
Each dot represents a carbon atom, and the lines connect neighboring atoms, and representing the possible hoppings paths for the electrons.
}
\label{fig:graphene}
\end{figure}
\subsection{Lattice}
\indent
The hexagonal lattice can be constructed by copying an elementary cell at every integer combination of
\begin{equation}
l_1:=\left(\frac{3}{2},\frac{\sqrt{3}}{2}\right)
,\quad
l_2:=\left(\frac{3}{2},-\frac{\sqrt{3}}{2}\right)
\label{laeo}
\end{equation}
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).
\bigskip
\begin{figure}
\hfil\includegraphics[width=12cm]{cell.pdf}
\caption{
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).
It is defined in such a way that $\forall x\in\Lambda$, $\forall k\in\hat\Lambda$,
\begin{equation}
e^{ikxL}=1.
\end{equation}
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 {\it Brillouin zone}, see figure\-~\ref{fig:brillouin}.
\bigskip
\begin{figure}
\hfil\includegraphics[width=4cm]{brillouin.pdf}
\caption{
The Brillouin zone.
The two dots are the Fermi points $p_F^{(\pm)}$, see below.
}
\label{fig:brillouin}
\end{figure}
\subsection{Hamiltonian}
\indent
We will consider a model of spin-$\frac12$ 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}$.
\bigskip
\indent
The Hamiltonian is split into two terms:
\begin{equation}
\mathcal H=\mathcal H_0+\mathcal H_I
\end{equation}
where $\mathcal H_0$ is the {\it free Hamiltonian}, which describes the motion of electrons from one atom to a neighbor, and $\mathcal H_I$ is the {\it interaction Hamiltonian}, which describes the interaction between electrons.
\bigskip
\point{\bf Free Hamiltonian.}
The free Hamiltonian describes the {\it hopping} of electrons from one atom to another:
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$.
The exponential of $\mathcal H$ is difficult to compute, due to the presence of both the kinetic term $\mathcal H_0$ and the interacting Hamiltonian $\mathcal H_I$.
We can split these from each other using the Trotter product formula:
\begin{equation}
e^{-\beta\mathcal H}=
\lim_{p\to\infty}
\left(
e^{-\frac\beta p\mathcal H_0}
e^{-\frac\beta p\mathcal H_I}
\right)^p
\end{equation}
(which follows from the Baker-Campbell-Hausdorff formula).
Taking the limit $p\to\infty$, we can expand the exponential of $\mathcal H_I$:
\begin{equation}
e^{-\beta\mathcal H}=
\lim_{p\to\infty}
\left(
e^{-\frac\beta p\mathcal H_0}
\left(1-\frac\beta p\mathcal H_I\right)
\right)^p
.
\end{equation}
We can expand the power by noting that, in each factor, either the term $1$ or $-\frac\beta p\mathcal H_I$ can be selected.
(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\geqslant0$,
(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:
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\-~\ref{lemma:wick}.
However, the time-ordering operator effectively makes them anticommute.
We can make this more precise be re-expressing the problem in terms of {\it Grassmann variables}.
\bigskip
\point{\bf 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\in2\pi\beta^{-1}(\mathbb Z+1/2)$, $k\in\hat\Lambda$, and $\mathbf k\equiv(k_0,k)$
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\leqslant1/3$ and $0$ if $\rho\geqslant2/3$, and let
To every $\hat\alpha_{\mathbf k,\sigma}^{\pm}$ for $\alpha\in\{a,b\}$ and $\mathbf k\in\mathcal B_{\beta,L}^*$, we associate a {\it 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\}}$.
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
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\-~\ref{lemma:grassmann_addition}),
\point{\bf 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\-~(\ref{wick}).
Gaussian Grassmann integrals of Grassmann polynomials also satisfy the Wick rule, see\-~(\ref{wick_grassmann}), so Gaussian Grassmann integrals and the free average share the same algebraic structure.
Actually, since we cut off the momenta by $M$, in order to avoid introducing Gibbs phenomena when inverting Fourier transforms, we define the {\it propagator}:
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\-~(\ref{freeen}):
Note that we have dropped the $-\frac12$ 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\-~\ref{lemma:schwinger}, the two point function $\left<\mathbf T(\alpha^-\alpha^+)\right>$ at $t=\bar t$ has an extra $\frac12$, and one can show that this exactly cancels with the $-\frac12$ in the interaction Hamiltonian.
\subsection{Singularities of the propagator and scale decomposition}
\indent
We thus have a clear strategy to compute the free energy of the electrons in graphene: (\ref{freeengrass}) 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\-~(\ref{fermipt}), around which they behave like $v_F|k-p_F^{(\omega)}|$, see\-~(\ref{dispersion}).
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
\begin{equation}
\pm|v_F|k-p_F^{(\omega)}|+ik_0|^{-1}
.
\end{equation}
(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).
\bigskip
\indent
We will need to proceed in a more subtle fashion, by performing a {\it 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
which is a smooth function that is supported in $|\mathbf k-\mathbf p_F^{(\omega)}|\in[2^h\frac16,2^h\frac23]$, 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\-~\ref{fig:scale}.
(note that $\hat g_{\geqslant0}$ is bounded uniformly in $\beta$) and use\-~(\ref{addprop}) to compute the integrals with the propagators $g^{(h,\omega)}$ one at a time.
\bigskip
\begin{figure}
\hfil\includegraphics[width=5cm]{scales.pdf}
\caption{
A sketch of the scale decomposition, the darker the color, the smaller the scale.
}
\label{fig:scale}
\end{figure}
\bigskip
\indent
This has been done in detail in\-~\cite{GM10,Gi10}.
Carrying out this strategy in practice is rather involved though, so here, we will simplify the problem by considering a {\it hierarchical model} based on the graphene model.
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\-~\cite{Dy69,BCe78,GK81}, though most of those are Bosonic.
Fermionic hierarchical models were first studied in\-~\cite{Do91}, and, in a more systematic way, in\-~\cite{BGJ15,GJ15}.
\bigskip
\subsection{Definition of the hierarchical model}
\point{\bf 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)}|\sim2^h$.
In configuration space, this corresponds to considering Grassmann fields that are constant on configuration-space boxes of size $\sim2^{-h}$.
These boxes are obtained by doubling the size of the elementary cell of the hexagonal lattice at every scale, see figure\-~\ref{fig:coarse}.
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
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)$.
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)
which has the same scaling as the non-hierarchical propagator\-~(\ref{propx}) (and similarly with $a$ and $b$ exchanged).
\bigskip
\point{\bf Effective potentials.}
In this hierarchical model, we compute
\begin{equation}
\int P(d\psi)\ e^{-\mathcal V(\psi)}
.
\end{equation}
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
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 {\it effective potential} on scale $h$, and it dictates the physical properties of the system at distances $\sim2^{-h}$.
Because there are only four Grassmann fields and their conjugates per cell, $v_h$ must be a {\it polynomial} in the Grassmann fields of order $\leqslant8$.
In fact, by symmetry considerations, we find that $v_h$ must be of the form
In addition, the integral over $\psi^{[h]}(\Delta)$ does not depend on $\Delta$.
Therefore, using the fact that $\bar\Delta$ contains $2^{d+1}\equiv8$ ($d\equiv2$ 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$,
The computation is thus reduced to computing the map $\ell^{(h)}\mapsto\ell^{(h-1)}$ using\-~(\ref{betadef}).
The coefficients $\ell_i^{(h)}$ are called {\it running coupling constants}, and the map $\ell^{(h)}\mapsto\ell^{(h-1)}$ is called the {\it 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 $\sim2^{-h}$.
\bigskip
\point{\bf Feynman diagrams.}
Having defined the hierarchical model as we have, the infinite sum in\-~(\ref{betadef}) is actually finite ($n\leqslant4$), 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 {\it Feynman diagrams}.
First, let us expand the power $n$ and graphically represent the terms that must be integrated.
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 {\it external}, or we can choose the term $\psi^{[h]}$, in which case this edge will have to be integrated, and will be called {\it 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\-~(\ref{wick_grassmann}), 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\-~(\ref{betadef}) can thus be reduced to a sum over graphs.
\bigskip
\indent
Let us make a few comments about this expansion, and how it is used for non-hierarchical models.
\begin{itemize}
\item
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\-~(\ref{betadef}) absolutely divergent.
However, in counting these graphs, we are ignoring important signs: in the Wick rule\-~(\ref{wick_grassmann}), each permutation comes with its signature.
These signs lead to important cancellations in the sum over $n$, which can be exploited using so-called {\it Gram bounds}, or {\it determinant expansions}\-~\cite{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\-~\cite[Appendix\-~B]{Gi10}.
\item
The diagram expansion allows us to compute the integral in\-~(\ref{betadef}).
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.
\item
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\'o\-~\cite{GN85} trees, as detailed in\-~\cite{Gi10}.
Another good reference for the tree expansion and its connection to Feynman diagrams is\-~\cite[Section\-~5]{GJ16}
This is not needed for the hierarchical model.
\end{itemize}
\bigskip
\point{\bf 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\-~(\ref{betadef}))
(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:
\begin{itemize}
\item
if $d+1-2l\gamma<0$, then this term will decrease exponentially as $h\to-\infty$; such running coupling constants will be called {\it irrelevant};
\item
if $d+1-2l\gamma>0$, then this term will grow exponentially as $h\to-\infty$; such running coupling constants will be called {\it relevant};
\item
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 {\it marginal}.
\end{itemize}
In the case of graphene, $d=2$ and $\gamma=1$, so if $2l=2$, then the running coupling constant is relevant, and if $2l\geqslant4$, it is irrelevant; there are no marginal constants.
\bigskip
\indent
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.
\bigskip
\indent
This {\it power counting} holds for the non-hierarchical model as well, as detailed in\-~\cite{Gi10}.
For a more general treatment of power counting in Fermionic models with point-singularities, see\-~\cite[Section\-~5.2]{GJ16}.
\bigskip
\indent
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 {\it 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 {\it super-renormalizable}.
\bigskip
\point{\bf Hierarchical beta function.}
As was mentioned above, the beta function can be computed {\it 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 {\tt meankondo} software package\-~\cite{mk}.
The resulting beta function contains 888 terms, and will not be written out here.
\subsection{Results for non-hierarchical graphene}
\indent
Using the Gram bounds and Gallavotti-Nicol\'o tree expansion hinted at above, the renormalization group analysis can be carried out in full for (non-hierarchical) graphene\-~\cite{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 {\it renormalizing} the non-interacting Hamiltonian).
In particular, the Fermi velocity $v_F$ in\-~(\ref{dispersion}) is renormalized.
\bigskip
\indent
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.
\bigskip
\indent
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\-~\cite{GJ16} for an analysis of bilayer graphene where this happens).
%\section{Hierarchical bilayer graphene}
%\indent
%Let us now discuss a model with less favorable scaling properties.
%We will consider {\it bilayer graphene} in AB stacking, which is a model that is defined similarly to graphene except that we add a second layer of graphene in the configuration depicted in figure\-~\ref{fig:bilayer_graphene}.
%\bigskip
%
%\begin{figure}
%\hfil\includegraphics[width=12cm]{bilayer.pdf}
%\caption{
% The bilayer graphene lattice.
% Each circle and square represents a carbon atom, and the lines connect neighboring atoms, and representing the intra-layer hoppings paths for the electrons.
% We also consider a single interlayer vertical hopping between full circles and empty squares.
%}
%\label{fig:bilayer_graphene}
%\end{figure}
%
%\subsection{Hamiltonian}
%\indent
%Here again, the Hamiltonian is split into two terms:
%The lattice can be constructed using the same elementary cell as single layer graphene, except that there are now four atoms in each cell.
%In one of the layers, the atoms are of type $a$ and $b$, as before, and in the other, they are of type $\tilde a$ and $\tilde b$, with the type $\tilde b$ atoms directly above type $a$ atoms.
%We will consider the same interlayer hoppings as for graphene, and add a single vertical hopping between type $a$ and type $\tilde b$ atoms:
%Two of these bands are bounded away from $0$: $\pm\frac{\gamma_1+\sqrt{\gamma_1^2+4|\Omega(k)|^2}}2$, whereas the other two meet at $0$ when $\Omega(k)=0$, that is, when $k=p_F^{(\omega)}$.
%By\-~(\ref{dispersion}), for $|k-p_F^{(\omega)}|\ll1$,
%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$.
%
%\subsection{Scale decomposition}
%\indent
%Since the free Hamiltonian vanishes at $p_F^{(\omega)}$, the propagator $\hat g$ is singular, and we will proceed in the same way as for graphene by performing a scale decomposition.
%There are two differences, first, only two of the bands induce singularities instead of all of them, and second, and most importantly, the bands vanish quadratically instead of linearly.
%\bigskip
%
%\point
%Dealing with the fact that only two of the bands are singular is simple: we only need the scale decomposition for the singular bands, so we will only perform it for those bands.
%In other words, we change variables to a basis in which $\hat{\bar H}_0$ is diagonal, split $\hat g$ into the terms corresponding to the non-singular bands and singular bands, and only define a scale decomposition for the singular bands.
%This is done in detail in\-~\cite{GJ16}.
%\bigskip
%
%\point
%We are left with the singular bands.
%By\-~(\ref{dispersion_bilayer}), near the singularities, the eigenvalues of $\hat g$ behave like
%however, when going back to configuration space, we have
%\begin{equation}
% \sum_{\mathbf k\in\mathcal B_{\beta,L}}
% \Phi_h(\mathbf k-\mathbf p_F^{(\omega)})
% =O(2^{4h})
%\end{equation}
%instead of $2^{3h}$ so we recover the estimate
%\begin{equation}
% g^{(h,\omega)}(t,x)=O(2^{2h})
%\end{equation}
%instead of $2^{2h}$.
%
%\subsection{Hierarchical beta function}
%\indent
%We define a hierarchical model for bilayer graphene, in which we only keep the bands that are singular (the non-singular bands have trivial scaling properties).
%(In this section, we use the indices $a,b$ for the Grassmann fields, even though they now represent the two singular bands rather than atoms of type $a$ and $b$.)
%We define the hierarchical model almost exactly as for graphene.
%We will only change the scaling the time-dimension of boxes to reflect the scaling in\-~(\ref{scaling_bilayer}): (\ref{boxes}) becomes
The graphene model discussed above is {\it 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\-~\cite{Ko64}, was one of the foundational models in the development stages of the renormalization group\-~\cite{Wi75}.
It is integrable\-~\cite{An80}, but is nevertheless an interesting model to study renormalization group flows in.
The hierarchical Kondo model was introduced and studied in\-~\cite{BGJ15,GJ15}.
\bigskip
\subsection{Hamiltonian}
\indent
The Kondo model is a one-dimensional model of spin-$\frac12$ electrons on a chain that interact with a localized magnetic impurity, which is represented as a spin-$\frac12$ 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:
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
\begin{equation}
\mathcal Q_h:=\left\{
[i2^{|h|},(i+1)2^{|h|})
\right\}_{i\in\mathbb Z}
.
\end{equation}
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
In accordance with the scaling in\-~(\ref{scaling_kondo}), we define, for $\sigma\in\{\uparrow,\downarrow\}$, $h\in\{-N_\beta,\cdots,-1\}$, $\Delta\in\mathcal Q_h$, $\eta\in\pm$
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.
}
\label{fig:sd_vector_field}
\end{figure}
\bigskip
\indent
Thus, for the hierarchical Kondo model, the flow can go to a {\it 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).
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)$,
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)}$.
\bigskip
\indent
Let us first prove a technical lemma.
\bigskip
\theo{Lemma}\label{lemma:fock_technical}
For $t\in\mathbb R$, and $\lambda_1,\cdots,\lambda_N\in\mathbb R$,
The proof of\-~(\ref{fock1}) proceeds as follows: we write $e^{-t\sum_i\lambda_ia_i^\dagger a_i}=\prod_ie^{-t\lambda_ia_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\geqslant1$.
Let us compute the Fourier transform of $s_{i,j}(t-t')$ with respect to $t-t'\in(-\beta,\beta)$: for $k_0\in\frac\pi\beta\mathbb Z$, by\-~(\ref{S1}) and\-~(\ref{S2}),
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, (\ref{wick_app}) is antisymmetric under exchanges of $a^{-}$'s and exchanges of $a^+$'s).
\bigskip
\point
We diagonalize $\mu$ by a unitary transform $U$, and define
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}\geqslant0$ be such that
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}^nb_{k_i}\prod_{i=n}^1b_{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
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}^nb_i$, but the $b_i^\dagger$ are ordered in the opposite order as in $\prod_{i=n}^1b_i^\dagger$.
As we will now discuss, reordering the $b$'s in this way yields
where $\mathds1_{\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}$.
First of all, $\prod_{i=n}^1b_{l_i}^\dagger$ is non-zero only if all $l_i$ are different from each other, and so is the right side of\-~(\ref{wick_notime}), 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\-~(\ref{wick_notime}), we can assume without loss of generality that $k_i=l_i$.
Let us now turn to a set of identities for Gaussian Grassmann integrals.
\bigskip
\theo{Lemma}\label{lemma:grassmann_id}
We have
\begin{equation}
\int P_\mu(d\psi)\ 1 = 1
\end{equation}
and
\nopagebreakaftereq
\begin{equation}
\int P_\mu(d\psi)\ \psi_i^-\psi_j^+=\mu_{i,j}
\end{equation}
where $\mu^{-1}$ is the inverse matrix of $\mu$.
\endtheo
\restorepagebreakaftereq
\bigskip
\indent\underline{Proof}:
This is a consequence of the change of variables lemma\-~\ref{lemma:change_vars_grassmann}.
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$:
in which we insert\-~(\ref{exp_trivial}) to find that $k$ must be equal to $l$ and
\begin{equation}
\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}
.
\end{equation}
\qed
\subsection{Wick rule for Gaussian Grassmann integrals}
\indent
Let us now turn to the Wick rule for Gaussian Grassmann integrals.
\bigskip
\theoname{Lemma}{Wick rule for Gaussian Grassmann integrals}\label{lemma:wick_grassmann}
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\}$,
\nopagebreakaftereq
\begin{equation}
\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)}}^+
.
\end{equation}
\endtheo
\restorepagebreakaftereq
\bigskip
\indent\underline{Proof}:
As before, we change variables using lemma\-~\ref{lemma:change_vars_grassmann}.
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$:
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.
