Initial commitv1.0

34 files changed, 4340 insertions, 0 deletions

diff --git a/Jauslin_2022.tex b/Jauslin_2022.tex
new file mode 100644
index 0000000..f200478
--- /dev/null
+++ b/Jauslin_2022.tex
@@ -0,0 +1,2395 @@
+\documentclass{ian}
+
+\usepackage{graphicx}
+\usepackage{array}
+\usepackage{largearray}
+\usepackage{dsfont}
+\usepackage{iantheo}
+\usepackage[reset_at_theo]{point}
+
+\begin{document}
+
+\pagestyle{empty}
+
+\hbox{}
+\hfil{\bf\LARGE
+The Hierarchical Graphene model\par
+}
+\vfill
+
+\hfil{\bf\large Ian Jauslin}\par
+\hfil{\it Department of Mathematics, Rutgers University}\par
+\hfil{\tt\color{blue}\href{mailto:ian.jauslin@rutgers.edu}{ian.jauslin@rutgers.edu}}\par
+\vskip20pt
+
+\vfill
+
+{\it\small
+  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).
+}
+\label{fig:graphene}
+\end{figure}
+
+\indent
+We define the lattice 
+\begin{equation}
+  \Lambda:=\left\{n_1 l_1+n_2 l_2,\ (n_1,n_2)\in\{0,\cdots,L-1\}^2\right\}
+  \label{laet}
+\end{equation}
+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:
+\begin{equation}
+  \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).
+  \label{lea}
+\end{equation}
+\bigskip
+
+\indent
+The {\it dual} of $\Lambda$ is
+\begin{equation}
+  \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\}
+  \label{lae}
+\end{equation}
+with periodic boundary conditions, where
+\begin{equation}
+  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).
+  \label{laeG}
+\end{equation}
+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:
+\begin{equation}
+  \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})
+  \label{hamx}
+\end{equation}
+Equation\-~(\ref{hamx}) can be rewritten in Fourier space as follows.
+We define the Fourier transform of the annihilation operators as
+\begin{equation}
+  \hat a_{k,\sigma}:=\frac1{\sqrt{|\Lambda|}}\sum_{x\in\Lambda}e^{ikx}a_{x,\sigma}
+  ,\quad 
+  \hat b_{k,\sigma}:=\frac1{\sqrt{|\Lambda|}}\sum_{x\in\Lambda}e^{ikx}b_{x+\delta_1,\sigma}
+  .
+\end{equation} 
+Note that, with this choice of normalization, $\hat a_{k,\sigma}$ and $\hat b_{k,\sigma}$ satisfy the canonical anticommutation relations:
+\begin{equation}
+  \{a_{k,\sigma},a_{k',\sigma'}^\dagger\}
+  =
+  \frac1{|\Lambda|}\sum_{x,x'\in\Lambda}e^{ikx-ik'x'}
+  \{a_{x,\sigma},a_{x',\sigma'}^\dagger\}
+  =
+  \delta_{\sigma,\sigma'}
+  \frac1{|\Lambda|}\sum_{x'\in\Lambda}e^{i(k-k')x}
+  =
+  \delta_{\sigma,\sigma'}
+  \delta_{k,k'}
+\end{equation}
+and similarly for $b$.
+We express $\mathcal H_0$ in terms of $\hat a$ and $\hat b$:
+\begin{equation}
+  \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}
+  \label{hamk}
+\end{equation}
+where $|\Lambda|=L^2$, $\hat A_{k,\sigma}$ is a column vector whose transpose is $\hat A_{k,\sigma}^T=(\hat a_{k,\sigma},\hat{b}_{k,\sigma})$, 
+\begin{equation}
+  H_0( k):=
+  \left(\begin{array}{*{2}{c}}
+    0&\Omega^*(k)\\
+    \Omega(k)&0
+  \end{array}\right)
+  \label{hmat}
+\end{equation}
+and
+\begin{equation}
+  \Omega(k):=\sum_{j=1}^3e^{ik(\delta_j-\delta_1)}
+  =1+2e^{-i\frac32k_x}\cos({\textstyle\frac{\sqrt{3}}2k_y})
+  .
+  \label{Omega}
+\end{equation}
+The eigenvalues of $H_0(k)$ are called the {\it bands} of non-interacting graphene, and are
+\begin{equation}
+  \pm|\Omega(k)|
+  =\pm\left(
+    1+4\cos({\textstyle\frac32k_x})\cos({\textstyle\frac{\sqrt{3}}2k_y})
+    +4\cos^2({\textstyle\frac{\sqrt{3}}2k_y})
+  \right)^{\frac12}
+  .
+  \label{bands}
+\end{equation}
+These bands meet at $0$ at exactly two values of $k$: for $\omega=\pm$,
+\begin{equation}
+  p_F^{(\omega)}:=({\textstyle\frac{2\pi}3,\omega\frac{2\pi}{3\sqrt3}})
+  \label{fermipt}
+\end{equation}
+see figure\-~\ref{fig:bands}.
+For $|k-p_F^{(\omega)}|\ll1$,
+\begin{equation}
+  \pm|\Omega(k)|\sim\pm v_F|k-p_F^{(\omega)}|
+  ,\quad
+  v_F=\frac32
+  .
+  \label{dispersion}
+\end{equation}
+\bigskip
+
+\begin{figure}
+\hfil\includegraphics[width=8cm, trim={1cm 1cm 1cm 1cm}, clip]{bands.pdf}
+\caption{
+  The bands of graphene.
+  The two bands meet at two points in a conical intersection.
+}
+\label{fig:bands}
+\end{figure}
+\bigskip
+
+\point{\bf Interaction.}
+We now define the interaction Hamiltonian which we take to be of {\it Hubbard} form:
+\begin{equation}
+  \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)
+\label{hamintx}\end{equation}
+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{Grassmann integral representation}
+\indent
+The {\it specific free energy} on the lattice $\Lambda$ is defined by
+\begin{equation}
+  f_\Lambda:=-\frac{1}{\beta|\Lambda|}\log\left(\mathrm{Tr}\left( e^{-\beta\mathcal H}\right)\right)
+  \label{freeen}
+\end{equation}
+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\-~\cite[(4.1)]{Gi10})
+\begin{equation}
+  \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\geqslant0}\mathrm{Tr}\left(e^{-\beta\mathcal H_0}\mathcal H_I(t_1)\cdots\mathcal H_I(t_N)\right)
+\end{equation}
+where
+\begin{equation}
+  \mathcal H_I(t):=e^{t\mathcal H_0}\mathcal H_Ie^{-t\mathcal H_0}
+  .
+\end{equation}
+To compute this trace, we will use the {\it Wick rule}, which we will now descibe.
+First, we define the {\it free average}:
+\begin{equation}
+  \left<A\right>:=\frac{\mathrm{Tr}(e^{-\beta\mathcal H_0}A)}{\mathrm{Tr}(e^{-\beta\mathcal H_0})}
+  .
+\end{equation}
+Next, we define the imaginary time creation and annihilation operators: for $\alpha\in\{a,b\}$ and $t\in[0,\beta)$,
+\begin{equation}
+  \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}
+\end{equation}
+(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$,
+\begin{equation}
+  \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>
+  \label{wick}
+\end{equation}
+(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:
+\begin{equation}
+  \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.
+\end{equation}
+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}.
+\bigskip
+
+\indent
+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 {\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)$
+\begin{equation}
+  \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
+  \frac1{\sqrt{|\Lambda|}}\sum_{x\in\Lambda}\int_0^\beta dt\ e^{\mp(itk_0+ikx)}\alpha_{x,\sigma}^\pm(t)
+\end{equation}
+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
+\begin{equation}
+  \mathcal B_{\beta,L}^*:=\mathcal B_{\beta,L}\cap\{(k_0,k),\ \chi_0(2^{-M}|k_0|\neq0)\}
+  .
+\end{equation}
+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\}}$.
+We define the Grassmann integral
+\begin{equation}
+  \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}^-
+\end{equation}
+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
+\begin{equation}
+  \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
+\end{equation}
+along with the anti-commutation of the $\hat\psi$.
+\bigskip
+
+\indent
+Let us now define {\it Gaussian Grassmann integrals}.
+The Gaussian Grassmann measure is specified by a {\it propagator}, which is a $2\times2$ complex matrix $\hat g(\mathbf k)$:
+\begin{equation}
+  \begin{largearray}
+    P_{\hat g}(d\psi) := \left(
+      \prod_{\mathbf k\in\mathcal B_{\beta,L}^*}
+      (\beta\det\hat g(\mathbf k))^4
+      \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\\[0.5cm]\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)
+    .
+    \label{grassgauss}
+  \end{largearray}
+\end{equation}
+By a direct computation, one can prove that (see lemma\-~\ref{lemma:grassmann_id})
+\begin{equation}
+  \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)
+  .
+  \label{twopt_grassmann}
+\end{equation}
+In addition, the integral with respect to $P_{\hat g}(d\psi)$ satisfies the Wick rule (see lemma\-~\ref{lemma:wick_grassmann}):
+\begin{equation}
+  \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)}}^+
+  \label{wick_grassmann}
+\end{equation}
+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}),
+\begin{equation}
+  \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)
+  .
+  \label{addprop}
+\end{equation}
+\bigskip
+
+\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.
+Thus, if we set $\hat g$ in such a way that
+\begin{equation}
+  \int P_{\hat g}(d\psi)\ \hat\psi_{\mathbf k,\alpha,\sigma}^-\hat\psi_{\bar\mathbf k,\bar\alpha,\bar\sigma}^+
+  =
+  \frac1{|\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>
+\end{equation}
+then we can compute free averages using Gaussian Grassmann integrals.
+Furthermore, by a direct computation (see lemma\-~\ref{lemma:schwinger}),
+\begin{equation}
+  \frac1{|\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}
+\end{equation}
+so, by\-~(\ref{twopt_grassmann}),
+\begin{equation}
+  \hat g(\mathbf k)=(-ik_0\mathds 1-H_0(k))^{-1}
+  .
+\end{equation}
+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}:
+\begin{equation}
+  \hat g_{\leqslant M}(\mathbf k):=\chi_0(2^{-M}|k_0|) (-ik_0\mathds 1-H_0(k))^{-1}
+  \label{freeprop}
+\end{equation}
+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}):
+\begin{equation}
+  f_\Lambda=f_{0,\Lambda}-\lim_{M\to\infty}\frac{1}{\beta|\Lambda|}\log\int P_{\leqslant M}(d\psi)\ e^{-\mathcal V(\psi)}
+  \label{freeengrass}
+\end{equation}
+where $f_{0,\Lambda}$ is the free energy in the $U=0$ case and
+\begin{equation}
+  \mathcal V(\psi)=U\sum_{\alpha\in\{a,b\}}\frac1{|\Lambda|}\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}
+  \label{V_grassmann}
+\end{equation}
+in which $\mathbf x\equiv(t,x)$ and
+\begin{equation}
+  \psi^\pm_{\mathbf x,\alpha,\sigma}:=\frac1{\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}
+  .
+\end{equation}
+\bigskip
+
+\indent
+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
+\begin{equation}
+  \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)}|))
+  \label{fh}
+\end{equation}
+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, if localizes $\mathbf k$ to be at a distance from $\mathbf p_F^{(\omega)}$ that is of order $2^h$, see figure\-~\ref{fig:scale}.
+Since $|k_0|\geqslant\frac\pi\beta$, we only need to consider
+\begin{equation}
+  h\geqslant -N_\beta:=\log_2\frac\pi\beta
+  .
+\end{equation}
+We then define
+\begin{equation}
+  \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)})
+  \frac1{k_0^2+|\Omega(k)|^2}
+  \left(\begin{array}{cc}
+    ik_0&-\Omega^*(k)\\
+    -\Omega(k)&ik_0
+  \end{array}\right)
+  \label{prop}
+\end{equation}
+(where we used\-~(\ref{hmat}) to compute the inverse).
+For $h$ sufficiently small (sufficiently negative), we thus have
+\begin{equation}
+  \hat g^{(h,\omega)}=O(2^{-h}).
+  \label{propscale}
+\end{equation}
+To carry out the Gaussian Grassmann integral, we split
+\begin{equation}
+  \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)
+\end{equation}
+(note that $\hat g_{\geqslant 0}$ 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.
+
+\section{Hierarchical graphene}\label{sec:hierarchical_graphene}
+\indent
+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 $\sim 2^{-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
+\begin{equation}
+  \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}
+  \label{boxes}
+\end{equation}
+($\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)$.
+\bigskip
+
+\begin{figure}
+\hfil\includegraphics[width=12cm]{cell-coarse.pdf}
+\caption{
+  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$).
+}
+\label{fig:coarse}
+\end{figure}
+\bigskip
+
+\point{\bf Propagators.}
+Let us now define the propagators associated to these Grassmann fields.
+We will chose them to be {\it similar}, in a certain sense, to the non-hierarchical propagators.
+We defined the single scale propagator in Fourier space in\-~(\ref{prop}), which, in configuration space, becomes
+\begin{equation}
+  g^{(h,\omega)}(t,x)=
+  \frac1{\beta|\Lambda|}\sum_{\mathbf k\in\mathcal B_{\beta,L}}e^{-ik_0t-ikx}
+  \Phi_h(\mathbf k-\mathbf p_F^{(\omega)})
+  \frac1{k_0^2+|\Omega(k)|^2}
+  \left(\begin{array}{cc}
+    ik_0&-\Omega^*(k)\\
+    -\Omega(k)&ik_0
+  \end{array}\right)
+  .
+\end{equation}
+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,
+\begin{equation}
+  \frac1{\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})
+\end{equation}
+so
+\begin{equation}
+  g^{(h,\omega)}(t,x)=O(2^{2h})
+  .
+  \label{propx}
+\end{equation}
+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)
+\begin{equation}
+  \psi_{\alpha,\sigma}^\pm(t,x)=\sum_{h=-N_\beta}^02^h\psi_{\alpha,\sigma}^{[h]\pm}(\Delta^{(h)}(t,x))
+  .
+  \label{psi_hierarchical}
+\end{equation}
+We will take the propagators to be
+\begin{equation}
+  \int P^{[h]}(d\psi^{[h]})\ \psi_{a,\sigma}^{[h]-}(\Delta)\psi_{b,\sigma'}^{[h]+}(\Delta')=\delta_{\sigma,\sigma'}\delta_{\Delta,\Delta'}
+\end{equation}
+\begin{equation}
+  \int P^{[h]}(d\psi^{[h]})\ \psi_{b,\sigma}^{[h]-}(\Delta)\psi_{a,\sigma'}^{[h]+}(\Delta')=\delta_{\sigma,\sigma'}\delta_{\Delta,\Delta'}
+  .
+\end{equation}
+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
+\begin{equation}
+  \int P(d\psi)\ \psi_{a,\sigma}^-(x,t)\psi_{b,\sigma'}^+(y,t)
+  =
+  \delta_{\sigma,\sigma'}
+  \sum_{h=-N_\beta}^{\eta(x,y)}
+  2^{2h}
+  =
+  \delta_{\sigma,\sigma'}
+  \frac432^{2\eta(x,y)}(1-4^{-N_\beta-\eta(x,y)})
+\end{equation}
+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
+\begin{equation}
+  \psi_{\alpha,\sigma}^{[\leqslant h]\pm}(\Delta)
+  :=
+  \frac12\psi_{\alpha,\sigma}^{[\leqslant h-1]\pm}(\bar\Delta)
+  +
+  \psi_{\alpha,\sigma}^{[h]\pm}(\Delta)
+  .
+  \label{scaling_psi}
+\end{equation}
+Note that\-~(\ref{psi_hierarchical}) is then
+\begin{equation}
+  \psi_{\alpha,\sigma}^\pm(t,x)\equiv\psi_{\alpha,\sigma}^{[\le0]\pm}(\Delta^{(1)}(t,x))
+  .
+\end{equation}
+We then define, for $h\in\{-N_\beta,\cdots0\}$,
+\begin{equation}
+  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^{[\le0]}):=\mathcal V(\psi^{[\le0]})
+  \label{indV}
+\end{equation}
+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 $\sim 2^{-h}$.
+By a straightforward induction, we then find that
+\begin{equation}
+  \int P(d\psi)\ e^{-\mathcal V(\psi)}
+  =\exp\left({\textstyle-\beta|\Lambda|\sum_{h=-N_\beta}^0c^{[h]}}\right)
+  .
+\end{equation}
+The specific free energy\-~(\ref{freeengrass}) is then
+\begin{equation}
+  f_\Lambda=f_{0,\Lambda}-\sum_{h=-N_\beta}^0c^{[h]}
+  .
+\end{equation}
+We are then left with computing $\mathcal V^{[h]}$ and $c^{[h]}$ using\-~(\ref{indV}).
+By\-~(\ref{V_grassmann}), $\mathcal V$ is local in $(t,x)$ and so, by induction, it takes the form
+\begin{equation}
+  \mathcal V^{[h]}(\psi^{[\leqslant h]})=-\sum_{\Delta\in\mathcal Q_h}v_h(\psi^{[\leqslant h]}(\Delta))
+  .
+  \label{box_dcmp}
+\end{equation}
+\bigskip
+
+\indent
+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 $\leqslant 8$.
+In fact, by symmetry considerations, we find that $v_h$ must be of the form
+\begin{equation}
+  v_h(\psi)=
+  \sum_{i=0}^6\alpha_i^{(h)}O_i(\psi)
+  \label{vh_rcc}
+\end{equation}
+with
+\begin{equation}
+  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}^{-}
+\end{equation}
+\begin{equation}
+  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}^{-}
+\end{equation}
+\begin{equation}
+  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}^-
+\end{equation}
+\begin{equation}
+  \begin{largearray}
+    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}^-
+    +\\[0.3cm]\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{largearray}
+\end{equation}
+\begin{equation}
+  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}^-
+  .
+\end{equation}
+
+\subsection{Beta function of the hierarchical model}
+\point{\bf Beta function.}
+We have thus introduced a strategy to compute $\mathcal V^{[h]}$ inductively: starting from $\mathcal V^{[h]}$, we have
+\begin{equation}
+  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]})}
+\end{equation}
+where $\gamma\equiv1$ is the scaling dimension of $\psi$ in\-~(\ref{scaling_psi}).
+Now, by\-~(\ref{box_dcmp}), is
+\begin{equation}
+  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))}
+  .
+\end{equation}
+If we group the right side in boxes on scale $h-1$, we find
+\begin{equation}
+  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))}
+  .
+\end{equation}
+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\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$,
+\begin{equation}
+  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}}
+  .
+\end{equation}
+We expand the exponential and use\-~(\ref{vh_rcc}):
+\begin{equation}
+  \begin{largearray}
+    \beta|\Lambda|c^{[h]}
+    +
+    \sum_{i=0}^6\alpha_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
+    \frac1{n!}
+    \left(\sum_{i=0}^6\alpha_i^{(h)}O_i\left(\psi^{[h]}(\Delta)+2^{-\gamma}\psi^{[\leqslant h-1]}(\bar\Delta)\right)\right)^n
+    .
+  \end{largearray}
+  \label{betadef}
+\end{equation}
+The computation is thus reduced to computing the map $\alpha^{(h)}\mapsto\alpha^{(h-1)}$ using\-~(\ref{betadef}).
+The coefficients $\alpha_i^{(h)}$ are called {\it running coupling constants}, and the map $\alpha^{(h)}\mapsto\alpha^{(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\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 {\it 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 $\alpha_iO_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 {\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}))
+\begin{equation}
+  2^{d+1-2l\gamma}\alpha_i^{(h)}
+  .
+\end{equation}
+Furthermore, this graph will contribute to the running coupling constant $\alpha_i$, and so, on scale $h-1$, we will have
+\begin{equation}
+  \alpha_i^{(h-1)}=
+  2^{d+1-2l\gamma}\alpha_i^{(h)}
+  +
+  \cdots
+\end{equation}
+(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\geqslant 4$, 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.
+A careful analysis of the beta function shows that there is an equilibrium point at $\alpha_i=0$ for $i=1,2,3,4,5,6$ and
+\begin{equation}
+  \alpha_0\in\{0,1\}
+  .
+\end{equation}
+The point with $\alpha_0=0$ is unstable, whereas $\alpha_0=1$ is stable.
+\bigskip
+
+\begin{figure}
+\hfil\includegraphics[width=12cm]{graphene_vector_field.pdf}
+\caption{
+  The projection of the directional vector field of the beta function for hierarchical graphene onto the $(\alpha_0,\alpha_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 $\alpha_0=1$ and $\alpha_i=0$ is clearly visible.
+}
+\label{fig:vector_field}
+\end{figure}
+
+\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:
+%\begin{equation}
+%  \bar{\mathcal H}=\bar{\mathcal H}_0+\bar{\mathcal H}_I
+%  .
+%\end{equation}
+%\bigskip
+%
+%\point{\bf Free Hamiltonian.}
+%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:
+%\begin{equation}
+%  \begin{largearray}
+%    \bar{\mathcal H}_0:=
+%    -\sum_{\sigma\in\{\uparrow,\downarrow\}}\sum_{\displaystyle\mathop{\scriptstyle x\in\Lambda}_{j=1,2,3}}
+%    \left( a_{x,\sigma}^\dagger b_{x +\delta_j,\sigma}+b_{x +\delta_j,\sigma}^\dagger a_{x,\sigma}+
+%    \tilde b_{x,\sigma}^\dagger \tilde a_{x-\delta_j,\sigma} +\tilde a_{x-\delta_j,\sigma}^\dagger\tilde b_{x,\sigma}\right)
+%    -\\\hfill
+%    -\gamma_1\sum_{\sigma\in\{\uparrow,\downarrow\}}\sum_{x\in\Lambda}\left( a_{x,\sigma}^\dagger \tilde b_{x,\sigma}+\tilde b_{x,\sigma}^\dagger a_{x,\sigma}\right)
+%    \label{hamx_bilayer}
+%  \end{largearray}
+%\end{equation}
+%which we rewrite in Fourier space, defining
+%\begin{equation}
+%  \hat a_{k,\sigma}:=\sum_{x\in\Lambda}e^{ikx}a_{x,\sigma}
+%  ,\quad 
+%  \hat b_{k,\sigma}:=\sum_{x\in\Lambda}e^{ikx}b_{x+\bar\delta_1,\sigma}
+%  ,\quad
+%  \hat{\tilde b}_{k,\sigma}:=\sum_{x\in\Lambda}e^{ikx}\tilde b_{x+\delta_1,\sigma}\;,\quad
+%  ,\quad
+%  \hat{\tilde a}_{k,\sigma}:=\sum_{x\in\Lambda}e^{ikx}\tilde a_{x-\delta_1,\sigma}\;,\quad
+%\end{equation} 
+%we have
+%\begin{equation}
+%  \bar{\mathcal H}_0=-\frac{1}{|\Lambda|}\sum_{\sigma\in\{\uparrow,\downarrow\}}\sum_{k\in\hat\Lambda}\hat{\bar A}_{k,\sigma}^\dagger\bar H_0(k)\hat{\bar A}_{k,\sigma}
+%  \label{hamk_bilayer}
+%\end{equation}
+%where $|\Lambda|=L^2$, $\hat{\bar A}_{k,\sigma}$ is a column vector, whose transpose is $\hat A_k^T=(\hat a_{k,\sigma},\hat{b}_{k,\sigma},\hat{\tilde a}_{k,\sigma},\hat{\tilde b}_{k,\sigma})$, 
+%\begin{equation}H_0( k):=
+%  \left(\begin{array}{*{4}{c}}
+%    0&\Omega^*(k)&0&\gamma_1\\
+%    \Omega(k)&0&0&0\\
+%    0&0&0&\Omega^*(k)\\
+%    \gamma_1&0&\Omega(k)&0\\
+%  \end{array}\right)
+%  \label{hmat_bilayer}
+%\end{equation}
+%in which $\Omega(k)$ is the same as\-~(\ref{Omega}):
+%\begin{equation}
+%  \Omega(k):=\sum_{j=1}^3e^{ik(\delta_j-\delta_1)}
+%  =1+2e^{-i\frac32k_x}\cos({\textstyle\frac{\sqrt{3}}2k_y})
+%  .
+%\end{equation}
+%The bands of this model are
+%\begin{equation}
+%  \pm\frac{\gamma_1\pm\sqrt{\gamma_1^2+4|\Omega(k)|^2}}2
+%  .
+%\end{equation}
+%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$,
+%\begin{equation}
+%  \pm\frac{\gamma_1-\sqrt{\gamma_1^2+4|\Omega(k)|^2}}2
+%  \sim
+%  \mp\frac{v_F^2}{\gamma_1^2}|k-p_F^{(\omega)}|^2
+%  ,\quad
+%  v_F=\frac32
+%  \label{dispersion_bilayer}
+%\end{equation}
+%so the bands meet {\it quadratically} instead of linearly.
+%\bigskip
+%
+%\point{\bf Interaction.}
+%The interaction will be the same as for graphene:
+%\begin{equation}
+%  \bar{\mathcal H}_I:=
+%  U\sum_{x\in\Lambda}\sum_{\alpha\in\{a,b,\tilde a,\tilde 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)
+%\label{hamintx}\end{equation}
+%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
+%\begin{equation}
+%  \pm\left|\frac{v_F^2}{\gamma_1^2}|k-p_F^{(\omega)}|^2+ik_0\right|^{-1}
+%\end{equation}
+%so we need to scale $k_0$ and $k-p_F^{(\omega)}$ differently: let
+%\begin{equation}
+%  \bar \Phi_{h}(\mathbf k-\mathbf p_F^{(\omega)}):=
+%  (\chi_0(2^{-h}(|k-p_F^{(\omega)}|)+2^{-2h}|k_0|)
+%  -\chi_0(2^{-h+1}|k-p_F^{(\omega)}|+2^{-2(h+1)}|k_0|))
+%  \label{scaling_bilayer}
+%\end{equation}
+%and define
+%\begin{equation}
+%  \hat g^{(h,\omega)}(\mathbf k)
+%  :=
+%  \bar \Phi_h(\mathbf k-\mathbf p_{F^{(\omega)}})(-ik_0\mathds 1-\bar H_0(k))^{-1}
+%  .
+%\end{equation}
+%We now have
+%\begin{equation}
+%  \hat g^{(h,\omega)}(\mathbf k)=O(2^{-2h})
+%\end{equation}
+%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
+%\begin{equation}
+%  \mathcal Q_h:=\left\{
+%    [i2^{2|h|},(i+1)2^{2|h|+1})\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}
+%  .
+%\end{equation}
+%The rest is defined in the same way, and we find the renormalization group flow to be given by the analogue of\-~(\ref{betadef}):
+%\begin{equation}
+%  \begin{largearray}
+%    \sum_{i=0}^6\alpha_i^{(h-1)}O_i(\psi^{[\leqslant h-1]}(\bar\Delta))
+%    =\\\hfill=
+%    2^{d+2}\log
+%    \int P(d\psi^{[h]}(\Delta))
+%    \sum_{n=0}^\infty
+%    \frac1{n!}
+%    \left(\sum_{i=0}^6\alpha_i^{(h)}O_i\left(\psi^{[h]}(\Delta)+2^{-\gamma}\psi^{[\leqslant h-1]}(\bar\Delta)\right)\right)^n
+%  \end{largearray}
+%  \label{betadef_bilayer}
+%\end{equation}
+%(the only difference with\-~(\ref{betadef}) is that $2^{d+1}$ is replaced with $2^{d+2}$) with $d=2$ and $\gamma=1$.
+%The power counting is thus determined by
+%\begin{equation}
+%  d+2-2l\gamma
+%  \equiv
+%  4-2l
+%\end{equation}
+%and thus, we find that quadratic terms are {\it relevant}, terms of order $\geqslant 6$ are {\it irrelevant}, and quartic terms are {\it marginal}.
+%This is a significant difference from single layer graphene, in which there were no marginal terms.
+%Thus, $\alpha_0$ is relevant, $\alpha_1$, $\alpha_2$, $\alpha_3$ and $\alpha_4$ are marginal, and $\alpha_5$ and $\alpha_6$ are irrelevant.
+
+\section{Hierarchical Kondo model}\label{sec:hierarchical_kondo}
+\indent
+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:
+\begin{equation}
+  \bar{\mathcal H}_0=\bar{\mathcal H}_0+\bar{\mathcal H_I}
+\end{equation}
+with
+\begin{equation}
+  \bar{\mathcal H}_0=-\frac12\sum_{\sigma\in\{\uparrow,\downarrow\}}\sum_{x\in\{-\frac L2+1,\cdots,\frac L2\}}a_{x,\sigma}^\dagger(a_{x+1,\sigma}+a_{x-1,\sigma})
+  \otimes\mathds 1
+\end{equation}
+and
+\begin{equation}
+  \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)}
+\end{equation}
+where $a_{x,\sigma}$ is the annihilation operator at $x$ with spin $\sigma$ and $S^{(j)}$ are Pauli matrices:
+\begin{equation}
+  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)
+  .
+\end{equation}
+In Fourier space, let
+\begin{equation}
+  \hat\Lambda_L:=\frac{2\pi}L\left\{-\frac L2+1,\cdots,\frac L2\right\}
+\end{equation}
+in term of which
+\begin{equation}
+  \bar{\mathcal H}_0=-\sum_{\sigma\in\{\uparrow,\downarrow\}}\frac1L\sum_{k\in\hat\Lambda_L}\cos(k)\hat a_{k,\sigma}^\dagger\hat a_{k,\sigma}
+  \otimes\mathds 1
+  .
+\end{equation}
+The propagator is thus
+\begin{equation}
+  g(t,x)=\frac1{\beta L}\sum_{k_0\in\frac{2\pi}\beta(\mathbb Z+\frac12)}\sum_{k\in\hat\Lambda_L}\ e^{-ik_0t}e^{-ikx}\frac1{-ik_0-\cos(k)}
+  .
+\end{equation}
+
+\subsection{Scale decomposition and hierarchical model}
+\indent
+Since the interaction is localized at $x=0$, we only need propagators at $x=0$:
+\begin{equation}
+  g(t,0)=\frac1{\beta L}\sum_{k_0\in\frac{2\pi}\beta(\mathbb Z+\frac12)}\sum_{k\in\hat\Lambda_L}\ e^{-ik_0t}\frac1{-ik_0-\cos(k)}
+  .
+\end{equation}
+We decompose the propagator into scales around the singularities $(k_0,k)=(0,\pm\frac\pi2)$:
+\begin{equation}
+  g^{(h,\omega)}(t,0)=\frac1{\beta L}\sum_{k_0\in\frac{2\pi}\beta(\mathbb Z+\frac12)}\sum_{k\in\hat\Lambda_L}\ e^{-ik_0t}\frac1{-ik_0-\cos(k)}\Phi_h(k_0,k)
+\end{equation}
+where $\Phi_h$ is defined as in\-~(\ref{fh})
+\begin{equation}
+  \Phi_{h}(k_0,k):=(\chi_0(2^{-h}|(k_0,k-\omega{\textstyle\frac\pi2})|)-\chi_0(2^{-h+1}|(k_0,k-\omega{\textstyle\frac\pi2})|))
+  .
+\end{equation}
+On scale $h$ we have
+\begin{equation}
+  \frac1{|ik_0+\cos(k)|}=O(2^{-h})
+  ,\quad
+  \sum_{k_0,k}\Phi_h(k_0,k)
+  =O(2^{2h})
+\end{equation}
+so
+\begin{equation}
+  g^{(h,\omega)}(t,0)=O(2^h)
+  \label{scaling_kondo}
+\end{equation}
+(as opposed to $2^{2h}$ for graphene).
+\bigskip
+
+\indent
+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
+\begin{equation}
+  \Delta_-=[i2^{|h|},(i+{\textstyle\frac12})2^{|h|})
+  ,\quad
+  \Delta_+=[(i+{\textstyle\frac12})2^{|h|},(i+1)2^{|h|})
+\end{equation}
+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$
+\begin{equation}
+  \psi_{\sigma}^{[\leqslant h]\pm}(\Delta_\eta)
+  :=
+  \frac1{\sqrt2}\psi^{[\leqslant h-1]\pm}(\Delta)+\psi_{\sigma}^{[h]\pm}(\Delta_\eta)
+\end{equation}
+(compare this to\-~(\ref{scaling_psi}) in which the scaling factor is $\frac12$ instead of $\frac1{\sqrt2}$).
+
+
+\subsection{Beta function}
+\indent
+For this model, it is more convenient to expand the exponential of the potential: we denote
+\begin{equation}
+  e^{-\mathcal V^{[h]}(\psi)}=\mathcal W^{[h]}(\psi)
+  .
+\end{equation}
+One can show that the effective potentials $\mathcal W$ are of the form
+\begin{equation}
+  \mathcal W^{[h]}(\psi^{[\leqslant h]})
+  =
+  \prod_{\Delta\in\mathcal Q_h}\prod_{\eta=\pm}(1+w_h(\psi^{[\leqslant h]}(\Delta_\eta)))
+  .
+\end{equation}
+For this model, there are only two running coupling constants:
+\begin{equation}
+  w_h(\psi)=\sum_{i=0}^1\ell_i^{(h)}O_i(\psi)
+\end{equation}
+and
+\begin{equation}
+  O_0(\psi):=\frac12\sum_{\sigma,\sigma'\in\{\uparrow,\downarrow\}}\sum_{j=1,2,3}\psi_\sigma^+S^{(j)}_{\sigma,\sigma'}\psi_{\sigma'}^-
+  \otimes S^{(j)}
+  ,\quad
+  O_1(\psi):=\frac12\left(\sum_{\sigma,\sigma'\in\{\uparrow,\downarrow\}}\sum_{j=1,2,3}\psi_\sigma^+S^{(j)}_{\sigma,\sigma'}\psi_{\sigma'}^-\right)^2
+  \otimes\mathds 1
+\end{equation}
+(the $O$ are a product of Grassmann variables and a matrix acting on the impurity).
+\bigskip
+
+\indent
+We compute the beta function in the same way as for graphene, see\-~(\ref{betadef}):
+\begin{equation}
+  \begin{largearray}
+    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{largearray}
+\end{equation}
+with $\gamma\equiv\frac12$.
+In this expression, the square plays the same role as the prefactor $2^{d+1}$ in\-~(\ref{betadef}).
+The power counting is determined by the sign of
+\begin{equation}
+  1+2l\gamma\equiv1+l
+\end{equation}
+so quadratic terms ($\ell_0$) are {\it marginal} and quartic terms ($\ell_1$) are {\it irrelevant}.
+\bigskip
+
+\indent
+The beta function can be computed exactly for the hierarchical model:
+\begin{equation}
+  \begin{array}{r@{\ }>{\displaystyle}l}
+    C^{(h)}=&1+\frac32(\ell_0^{(h)})^2+9(\ell_1^{(h)})^2\\[0.3cm]
+    \ell_0^{(h-1)}=&\frac1{C^{(h)}}\Big(\ell_0^{(h)}+3 \ell_0^{(h)}\ell_1 -(\ell_0^{(h)})^2\Big)\\[0.3cm]
+    \ell_1^{(h-1)}=&\frac1{C^{(h)}}\Big( \frac12\ell_1^{(h)} +\frac18(\ell_0^{(h)})^2\Big)
+    .
+  \end{array}
+\end{equation}
+One can show that there are two equilibrium points: $\ell_0=\ell_1=0$, as well as
+\begin{equation}
+\ell^*_0=-x_0\frac{1+5x_0}{1-4x_0},\quad \ell^*_1=\frac{x_0}3
+\label{kondoeqfixedpt}\end{equation}
+where $x_0\approx0.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\-~\ref{fig:sd_vector_field}
+\bigskip
+
+\begin{figure}
+\hfil\includegraphics[width=12cm]{sd_vector_field.pdf}
+\caption{
+  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).
+
+
+\vfill
+\eject
+\appendix
+
+\section{Free Fermions}\label{app:free_fermions}
+Let
+\begin{equation}
+  \mathcal H_0=\sum_{i,j=1}^N \mu_{i,j}a^\dagger_ia_j
+\end{equation}
+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
+\begin{equation}
+  \left<A\right>:=\frac{\mathrm{Tr}(e^{-\beta\mathcal H_0}A)}{\mathrm{Tr}(e^{-\beta\mathcal H_0})}
+\end{equation}
+and, for $t\in[0,\beta)$, let
+\begin{equation}
+  a_i^-(t):=e^{t\mathcal H_0}a_ie^{-t\mathcal H_0}
+  ,\quad
+  a_i^+(t):=e^{t\mathcal H_0}a_i^\dagger e^{-t\mathcal H_0}
+  .
+\end{equation}
+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)$,
+\begin{equation}
+  \mathbf T\left(\prod_{i=1}^na^{\omega_i}_{j_i}(t_i)\right)
+  =
+  (-1)^\tau\prod_{i=1}^na^{\omega_{\tau(i)}}_{j_{\tau(i)}}(t_{\tau(i)})
+\end{equation}
+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$,
+  \begin{equation}
+    e^{-t\sum_{i=1}^N\lambda_ia_i^\dagger a_i}
+    =
+    \prod_{i=1}^N(1+(e^{-t\lambda_i}-1)a_i^\dagger a_i)
+    \label{fock1}
+  \end{equation}
+  and
+  \nopagebreakaftereq
+  \begin{equation}
+    e^{t\sum_{i=1}^N\lambda_ia_i^\dagger a_i}a_j^\dagger e^{-t\sum_{i=1}^N\lambda_ia_i^\dagger a_i}
+    =e^{t\lambda_j}a_j^\dagger
+    ,\quad
+    e^{t\sum_{i=1}^N\lambda_ia_i^\dagger a_i}a_j e^{-t\sum_{i=1}^N\lambda_ia_i^\dagger a_i}
+    =e^{-t\lambda_j}a_j
+    .
+    \label{fock4}
+  \end{equation}
+\endtheo
+\restorepagebreakaftereq
+\bigskip
+
+\indent\underline{Proof}:
+  To prove\-~(\ref{fock1}), 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\geqslant 1$.
+  By\-~(\ref{fock1}),
+  \begin{equation}
+    e^{-t\sum_{i=1}^N\lambda_ia_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
+  \end{equation}
+  and, commuting $a_j^\dagger$ through using $\{a_i,a_j^\dagger\}=\delta_{i,j}$ and $(a_j^\dagger)^2=0$, we find
+  \begin{equation}
+    e^{-t\sum_{i=1}^N\lambda_ia_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)
+    .
+    \label{fock2}
+  \end{equation}
+  Similarly,
+  \begin{equation}
+    e^{-t\sum_{i=1}^N\lambda_ia_i^\dagger a_i}a_j
+    =
+    \left(\prod_{i=1}^n(1+(e^{-t\lambda_i}-1)a_i^\dagger a_i)\right)a_j
+  \end{equation}
+  and so, using $a_i^2=0$, we find
+  \begin{equation}
+    e^{-t\sum_{i=1}^N\lambda_ia_i^\dagger a_i}a_j
+    =
+    a_j
+    \prod_{i\neq j}
+    (1+(e^{-t\lambda_i}-1)a_i^\dagger a_i)
+    .
+    \label{fock3}
+  \end{equation}
+  Furthermore, taking the $\dagger$ of\-~(\ref{fock3}), we find
+  \begin{equation}
+    a_j^\dagger e^{-t\sum_{i=1}^N\lambda_ia_i^\dagger a_i}
+    =
+    \left(
+      \prod_{i\neq j}
+      (1+(e^{-t\lambda_i}-1)a_i^\dagger a_i)
+    \right)
+    a_j^\dagger
+  \end{equation}
+  and since
+  \begin{equation}
+    (1+(e^{-t\lambda_j}-1)a_j^\dagger a_j)a_j^\dagger
+    =
+    e^{-t\lambda_j}a_j^\dagger
+  \end{equation}
+  we have
+  \begin{equation}
+    a_j^\dagger e^{-t\sum_{i=1}^N\lambda_ia_i^\dagger a_i}
+    =
+    e^{t\lambda_j}e^{-t\sum_{i=1}^N\lambda_ia_i^\dagger a_i}a_j^\dagger 
+    .
+    \label{fock2'}
+  \end{equation}
+  This implies the first of\-~(\ref{fock4}).
+  Taking the $\dagger$ of\-~(\ref{fock2}) yields
+  \begin{equation}
+    a_je^{-t\sum_{i=1}^N\lambda_ia_i^\dagger a_i}
+    =e^{-t\lambda_j}\prod_{i\neq j}
+    (1+(e^{-t\lambda_i}-1)a_i^\dagger a_i)a_j
+  \end{equation}
+  and since
+  \begin{equation}
+    (1+(e^{-t\lambda_i}-1)a_j^\dagger a_k)a_j
+    =a_j
+  \end{equation}
+  we have
+  \begin{equation}
+    a_je^{-t\sum_{i=1}^N\lambda_ia_i^\dagger a_i}
+    =e^{-t\lambda_j}e^{-t\sum_{i=1}^N\lambda_ia_i^\dagger a_i}a_j
+    .
+    \label{fock3'}
+  \end{equation}
+  This implies the second of\-~(\ref{fock4}).
+\qed
+\bigskip
+
+\subsection{Two-point correlation function}
+\indent
+We now compute the two-point correlation function of $\left<\cdot\right>$.
+\bigskip
+
+\theo{Lemma}\label{lemma:schwinger}
+  For $i,j\in\{1,\cdots,N\}$ and $t,t'\in[0,\beta)$, if $t\neq t'$, then
+  \begin{equation}
+    s_{i,j}(t-t'):=\left<\mathbf T(a_i^-(t)a_j^+(t'))\right>
+    =
+    \frac1\beta\sum_{k_0\in\frac{2\pi}\beta(\mathbb Z+\frac12)}
+    e^{-ik_0(t-t')}(-ik_0\mathds 1+\mu)^{-1}_{i,j}
+    \label{schwinger}
+  \end{equation}
+  and if $t=t'$,
+  \nopagebreakaftereq
+  \begin{equation}
+    s_{i,j}(0):=\left<\mathbf T(a_i^-(t)a_j^+(t))\right>
+    =
+    \frac12\delta_{i,j}+
+    \frac1\beta\sum_{k_0\in\frac{2\pi}\beta(\mathbb Z+\frac12)}
+    (-ik_0\mathds 1+\mu)^{-1}_{i,j}
+    .
+    \label{schwinger}
+  \end{equation}
+\endtheo
+\restorepagebreakaftereq
+\bigskip
+
+\indent\underline{Proof}:
+  We diagonalize $\mu$ by a unitary transform $U$, and define
+  \begin{equation}
+    b_i:=\sum_{j=1}^NU_{j,i}^*a_j
+    ,\quad
+    b_i^\dagger:=\sum_{j=1}^NU_{j,i}a_j^\dagger
+  \end{equation}
+  in terms of which
+  \begin{equation}
+    e^{-t\mathcal H_0}
+    =
+    e^{-t\sum_{i=1}^N\lambda_ib_i^\dagger b_i}
+    =
+    \prod_{i=1}^N(1+(e^{-t\lambda_i}-1)b_i^\dagger b_i)
+  \end{equation}
+  where $\lambda_i$ are the eigenvalues of $\mu$.
+  Note that
+  \begin{equation}
+    \{b_i,b_j^\dagger\}
+    =
+    \sum_{k,l=1}^NU_{k,i}^*U_{l,j}\{a_k,a_l^\dagger\}
+    =
+    \sum_{k=1}^NU_{k,i}^*U_{k,j}
+    =\delta_{i,j}
+    .
+  \end{equation}
+  Thus
+  \begin{equation}
+    \mathrm{Tr}(e^{-\beta\mathcal H_0})
+    =\prod_{i=1}^N(1+e^{-\beta\lambda_i})
+    .
+  \end{equation}
+  Let
+  \begin{equation}
+    b_i^{-}(s):=e^{t\mathcal H_0}b_ie^{-t\mathcal H_0}
+    =\sum_{j=1}^NU_{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}^NU_{j,i}a_j^{+}(s)
+    .
+  \end{equation}
+  By\-~(\ref{fock4}),
+  \begin{equation}
+    b_i^-(s)=e^{-s\lambda_i}b_i
+    ,\quad
+    b_i^+(s)=e^{s\lambda_i}b_i^\dagger
+    .
+  \end{equation}
+  \bigskip
+
+  \point
+  Now, if $t\geqslant t'$,
+  \begin{equation}
+    \mathbf T(a_i^-(t)a_j^+(t'))
+    =a_i^-(t)a_j^+(t')
+    =\sum_{k,l=1}^N
+    U_{i,k}U_{j,l}^*
+    e^{-t\lambda_k+t'\lambda_l}b_kb_l^\dagger
+    .
+  \end{equation}
+  In order for the trace to be non-zero, we must take $k=l$, so
+  \begin{equation}
+    \mathrm{Tr}(e^{-\beta\mathcal H_0}\mathbf T(a_i^-(t)a_j^+(t')))
+    =\sum_{k=1}^N
+    U_{i,k}U_{j,k}^*
+    e^{-(t-t')\lambda_k}
+    \mathrm{Tr}(e^{-\beta\mathcal H_0}b_kb_k^\dagger)
+  \end{equation}
+  and since $b_kb_k^\dagger=1-b_k^\dagger b_k$,
+  \begin{equation}
+    \mathrm{Tr}(e^{-\beta\mathcal H_0}\mathbf T(a_i^-(t)a_j^+(t')))
+    =\sum_{k=1}^N
+    U_{i,k}U_{j,k}^*
+    e^{-(t-t')\lambda_k}
+    \prod_{l\neq k}(1+e^{-\beta\lambda_l})
+    .
+  \end{equation}
+  Therefore,
+  \begin{equation}
+    s_{i,j}(t-t')
+    =
+    \sum_{k=1}^N
+    U_{i,k}U_{j,k}^*
+    \frac{e^{-(t-t')\lambda_k}}{1+e^{-\beta\lambda_k}}
+    .
+    \label{S1}
+  \end{equation}
+  \bigskip
+
+  \point If $t<t'$, then
+  \begin{equation}
+    \mathbf T(a_i^-(t)a_j^+(t'))
+    =-a_j^+(t')a_i^-(t)
+    =-e^{t'\mathcal H_0}a_j^\dagger e^{-(t'-t)\mathcal H_0}a_i e^{-t\mathcal H_0}
+  \end{equation}
+  so, using the cyclicity of the trace,
+  \begin{equation}
+    \mathrm{Tr}(e^{-\beta\mathcal H_0}\mathbf T(a_i^-(t)a_j^+(t')))
+    =-\mathrm{Tr}(e^{t\mathcal H_0}a_i e^{-(t+\beta-t')\mathcal H_0}a_j^\dagger e^{-t'\mathcal H_0})
+  \end{equation}
+  which is identical to the trace computed above, but with $t$ replaced by $t+\beta$, so
+  \begin{equation}
+    s_{i,j}(t-t')
+    =
+    -\sum_{k=1}^NU_{i,k}U_{j,k}^*\frac{e^{-(t-t'+\beta)\lambda_k}}{1+e^{-\beta\lambda_k}}
+    .
+    \label{S2}
+  \end{equation}
+  \bigskip
+
+  \point
+  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}),
+  \begin{equation}
+    \frac12\int_{-\beta}^\beta d\tau\ e^{ik_0\tau}s_{i,j}(t-t')
+    =
+    \sum_{k=1}^NU_{i,k}U_{j,k}^*
+    \frac1{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)
+  \end{equation}
+  and so
+  \begin{equation}
+    \frac12\int_{-\beta}^\beta d\tau\ e^{ik_0\tau}s_{i,j}(t-t')
+    =
+    \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)
+    .
+  \end{equation}
+  Thus, if $k_0\in\frac\pi\beta2\mathbb Z$, then
+  \begin{equation}
+    \frac12\int_{-\beta}^\beta d\tau\ e^{ik_0\tau}s_{i,j}(t-t')
+    =
+    0
+  \end{equation}
+  and if $k_0\in\frac{2\pi}\beta(\mathbb Z+\frac12)$, then
+  \begin{equation}
+    \frac12\int_{-\beta}^\beta d\tau\ e^{ik_0\tau}s_{i,j}(t-t')
+    =
+    \sum_{k=1}^NU_{i,k}U_{j,k}^*
+    \frac1{-ik_0+\lambda_k}
+    .
+  \end{equation}
+  Thus, if the Fourier transform can be inverted, then we have
+  \begin{equation}
+    s_{i,j}(t-t')
+    =
+    \frac1\beta\sum_{k_0\in\frac{2\pi}\beta(\mathbb Z+\frac12)}
+    e^{-ik_0(t-t')}(-ik_0\mathds 1+\mu)^{-1}_{i,j}
+  \end{equation}
+  and this is the case where $s_{i,j}(t-t')$ is continuous.
+  \bigskip
+
+  \point
+  Obviously, $s_{i,j}$ is continuous at $t-t'\neq0$.
+  At $t-t'=0$, we have
+  \begin{equation}
+    \lim_{t-t'\to 0_+}s_{i,j}(t-t')=\lim_{t-t'\to 0_-}s_{i,j}(t-t')+\delta_{i,j}
+  \end{equation}
+  so $s_{i,j}$ is continuous if $i\neq j$.
+  Now, if $i=j$ and $t=t'$, then by\-~(\ref{S1}),
+  \begin{equation}
+    s_{i,j}(0)=\sum_{k=1}^NU_{i,k}U_{i,k}^*\frac1{1+e^{-\beta\lambda_k}}
+    .
+  \end{equation}
+  Now,
+  \begin{equation}
+    s_{i,j}(t-t')-\frac12\mathrm{sgn}(t-t')
+  \end{equation}
+  is continuous at $t-t'=0$, and its Fourier transform is, by a similar computation to that above,
+  \begin{equation}
+    \frac12\int_{-\beta}^\beta d\tau\ \left(s_{i,j}(\tau)-\frac12\mathrm{sgn}(\tau)\right)
+    =
+    \sum_{k=1}^NU_{i,k}U_{j,k}^*
+    \left(\frac1{-ik_0+\lambda_k}+\frac1{ik_0}\right)
+  \end{equation}
+  which is absolutely summable in $k_0$.
+  Thus
+  \begin{equation}
+    \lim_{\tau\to0_+}\frac1\beta\sum_{k_0\in\frac{2\pi}\beta(\mathbb Z+\frac12)}
+    e^{-ik_0\tau}
+    \left((-ik_0\mathds 1+\mu)^{-1}+\frac1{ik_0}\mathds 1\right)
+    =
+    \frac1\beta\sum_{k_0\in\frac{2\pi}\beta(\mathbb Z+\frac12)}
+    \left((-ik_0\mathds 1+\mu)^{-1}+\frac1{ik_0}\mathds 1\right)
+  \end{equation}
+  and so
+  \begin{equation}
+    \lim_{\tau\to0_+}\left(s_{i,j}(\tau)-\frac12\mathrm{sgn}(\tau)\right)
+    =
+    \frac1\beta\sum_{k_0\in\frac{2\pi}\beta(\mathbb Z+\frac12)}
+    (-ik_0\mathds 1+\mu)^{-1}
+  \end{equation}
+  (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
+  \begin{equation}
+    s_{i,j}(0)=\frac12
+    +\frac1\beta\sum_{k_0\in\frac{2\pi}\beta(\mathbb Z+\frac12)}
+    (-ik_0\mathds 1+\mu)^{-1}
+    .
+  \end{equation}
+\qed
+
+\subsection{Wick rule}
+\indent
+This system of {\it free Fermions} satisfies the Wick rule, stated below.
+\bigskip
+
+\theoname{Lemma}{Wick rule}\label{lemma:wick}
+  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)$,
+  \nopagebreakaftereq
+  \begin{equation}
+    \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>
+    .
+    \label{wick_app}
+  \end{equation}
+\endtheo
+\restorepagebreakaftereq
+\bigskip
+
+\indent\underline{Proof}:
+  First, note that
+  \begin{equation}
+    \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}^1a_{\bar j_i}^+(\bar t_i)\right)\right)\right>
+    .
+    \label{fock_split_wick}
+  \end{equation}
+  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
+  \begin{equation}
+    b_i:=\sum_{j=1}^NU_{j,i}^*a_j
+    ,\quad
+    b_i^\dagger:=\sum_{j=1}^NU_{j,i}a_j^\dagger
+  \end{equation}
+  in terms of which
+  \begin{equation}
+    e^{-t\mathcal H_0}
+    =
+    e^{-t\sum_{i=1}^N\lambda_ib_i^\dagger b_i}
+    =
+    \prod_{i=1}^N(1+(e^{-t\lambda_i}-1)b_i^\dagger b_i)
+  \end{equation}
+  where $\lambda_i$ are the eigenvalues of $\mu$.
+  Note that
+  \begin{equation}
+    \{b_i,b_j^\dagger\}
+    =
+    \sum_{k,l=1}^NU_{k,i}^*U_{l,j}\{a_k,a_l^\dagger\}
+    =
+    \sum_{k=1}^NU_{k,i}^*U_{k,j}
+    =\delta_{i,j}
+  \end{equation}
+  Let
+  \begin{equation}
+    b_i^{-}(s):=e^{t\mathcal H_0}b_ie^{-t\mathcal H_0}
+    =\sum_{j=1}^NU_{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}^NU_{j,i}a_j^{+}(s)
+    .
+  \end{equation}
+  By\-~(\ref{fock4}),
+  \begin{equation}
+    b_i^-(s)=e^{-s\lambda_i}b_i
+    ,\quad
+    b_i^+(s)=e^{s\lambda_i}b_i^\dagger
+    .
+  \end{equation}
+  We have
+  \begin{equation}
+    \begin{largearray}
+      \left<\mathbf T\left(\left(\prod_{i=1}^{n}a_{j_i}^-(t_i)\right)\left(\prod_{i=n}^1a_{\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}^nU_{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}^1b_{l_i}^+(\bar t_i)\right)\right)\right>
+      .
+    \end{largearray}
+    \label{fock_diagonalization}
+  \end{equation}
+  \bigskip
+
+  \point
+  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
+  \begin{equation}
+    \mathbf T\left(\left(\prod_{i=1}^{n}b_{k_i}^-(t_i)\right)\left(\prod_{i=n}^1b_{l_i}^+(\bar t_j)\right)\right)
+    =
+    \sigma\prod_{i=1}^{2n}b_{\gamma_i}^{\omega_i}(s_i)
+  \end{equation}
+  and, by\-~(\ref{fock4}),
+  \begin{equation}
+    \mathbf T\left(\left(\prod_{i=1}^{n}b_{k_i}^-(t_i)\right)\left(\prod_{i=n}^1b_{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}
+    \label{fock_timeorder}
+  \end{equation}
+  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
+  \begin{equation}
+    \left|b^\dagger_{l_n}\cdots b^\dagger_{l_1}\right\}
+    :=
+    \left(\prod_{i=1}^nb_{k_i}\right)\left(\prod_{i=n}^1b^\dagger_{l_i}\right)
+  \end{equation}
+  and make the symbol $\left|\cdot\right\}$ multilinear: for instance,
+  \begin{equation}
+    \left|(b^\dagger_{l_n}-1)\cdots b^\dagger_{l_1}\right\}
+    \equiv
+    \left(\prod_{i=1}^nb_{k_i}\right)\left(\prod_{i=n}^1b^\dagger_{l_i}\right)
+    -
+    \left(\prod_{i=1}^{n-1}b_{k_i}\right)\left(\prod_{i=n-1}^1b^\dagger_{l_i}\right)
+    .
+  \end{equation}
+  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
+  \begin{equation}
+    \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\}
+    \label{fock_reorder}
+  \end{equation}
+  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\-~(\ref{fock_reorder}).
+  \bigskip
+
+  \point
+  We show that
+  \begin{equation}
+    \left<\left|\prod_{i=n}^1b_{l_i}^\dagger\right\}\right>
+    \equiv
+    \left<\left(\prod_{i=1}^nb_{k_i}\right)\left(\prod_{i=n}^1b_{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>
+    .
+    \label{wick_notime}
+  \end{equation}
+  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$.
+  In that case,
+  \begin{equation}
+    \mathrm{Tr}\left(e^{-\beta\mathcal H_0}\left(\prod_{i=1}^nb_{k_i}\right)\left(\prod_{i=n}^1b_{l_i}^\dagger\right)\right)
+    =
+    \mathrm{Tr}\left(e^{-\beta\mathcal H_0}\left(\prod_{i=1}^nb_{k_i}\right)\left(\prod_{i=n}^1b_{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)
+  \end{equation}
+  and so
+  \begin{equation}
+    \left<e^{-\beta\mathcal H_0}\left(\prod_{i=1}^nb_{k_i}\right)\left(\prod_{i=n}^1b_{l_i}^\dagger\right)\right>
+    =
+    \prod_{i=1}^n\left<b_{k_i}b_{k_i}^\dagger\right>
+  \end{equation}
+  which is equal to the right side of\-~(\ref{wick_notime}).
+  \bigskip
+
+  \point
+  By\-~(\ref{wick_notime}),
+  \begin{equation}
+    \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)
+    .
+    \label{fock_det}
+  \end{equation}
+  Now,
+  \begin{equation}
+    \left<b_{k_i}b_{l_{\tau(i)}}^\dagger\right>
+    =
+    \delta_{k_i,l_{\tau(i)}}
+    \frac1{1+e^{-\beta\lambda_{k_i}}}
+    \label{2pt_notime_+}
+  \end{equation}
+  and
+  \begin{equation}
+    \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}}}
+    .
+    \label{2pt_notime_-}
+  \end{equation}
+  Thus, inserting\-~(\ref{2pt_notime_+}), (\ref{2pt_notime_-}) into\-~(\ref{fock_det}), (\ref{fock_reorder}), (\ref{fock_timeorder}), (\ref{fock_diagonalization}) and\-~(\ref{fock_split_wick}), we find
+  \begin{equation}
+    \begin{largearray}
+      \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}^nU_{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{largearray}
+  \end{equation}
+  and, relabeling the $l_i$,
+  \begin{equation}
+    \begin{largearray}
+      \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{largearray}
+  \end{equation}
+  which, by\-~(\ref{S1}) and\-~(\ref{S2}), is
+  \begin{equation}
+    \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>
+    .
+  \end{equation}
+\qed
+
+\section{Properties of Gaussian Grassman integrals}\label{app:grassmann}
+\indent
+Consider a family $\psi_1^\pm,\cdots,\psi_N^\pm$ of Grassmann variables and the Gaussian Grassmann measure
+\begin{equation}
+  P_\mu(d\psi)=\det(\mu)\prod_{i=1}^Nd\psi_i^+d\psi_i^-
+  e^{-\sum_{i,j=1}^N\mu_{i,j}\psi_i^+\mu^{-1}_{i,j}\psi_i^-}
+\end{equation}
+where $\mu$ is a Hermitian, invertible matrix.
+\bigskip
+
+\indent
+We first prove that one can perform changes of variables in Grassmann integrals.
+
+\theoname{Lemma}{Change of variables in Grassmann integrals}\label{lemma:change_vars_grassmann}
+  Given a unitary $N\times N$ matrix $U$ and a polynomial $f$ in the Grassmann variables $\psi$, if
+  \begin{equation}
+    \varphi_i^-:=\sum_{j=1}^NU_{i,j}\psi_j^-
+    ,\quad
+    \varphi_i^+:=\sum_{j=1}^NU_{i,j}^*\psi_j^+
+  \end{equation}
+  then
+  \nopagebreakaftereq
+  \begin{equation}
+    \int\prod_{i=1}^Nd\varphi_i^+d\varphi_i^-\ f(U^\dagger\varphi^-,(U^\dagger)^*\varphi^+)
+    =
+    \int \prod_{i=1}^Nd\psi_i^+d\psi_i^-\ f(\psi^-,\psi^+)
+    .
+  \end{equation}
+\endtheo
+\restorepagebreakaftereq
+\bigskip
+
+\indent\underline{Proof}:
+  Without loss of generality, we can assume that the highest order term of $f$ is
+  \begin{equation}
+    \alpha\prod_{i=1}^N\psi_i^-\psi_i^+
+    .
+  \end{equation}
+  We have
+  \begin{equation}
+    \int \prod_{i=1}^Nd\psi_i^+d\psi_i^-\ f(\psi^-,\psi^+)=\alpha
+    .
+  \end{equation}
+  Furthermore,
+  \begin{equation}
+    \int\prod_{i=1}^Nd\varphi_i^+d\varphi_i^-\ f(U^\dagger\varphi^-,U\varphi^+)
+    =
+    \alpha\int\prod_{i=1}^Nd\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)
+  \end{equation}
+  and, avoiding repetitions in the products, we find that
+  \begin{equation}
+    \prod_{i=1}^N\left(
+      \sum_{j,k=1}^N
+      U_{i,j}^*\varphi_j^-
+      U_{i,k}\varphi_k^+
+    \right)
+    =
+    \sum_{\tau,\tau'\in\mathcal S_N}
+    \prod_{i=1}^NU^*_{i,\tau(i)}U_{i,\tau'(i)}\varphi_{\tau(i)}^-\varphi_{\tau'(i)}^+
+  \end{equation}
+  and so, reordering the $\varphi$,
+  \begin{equation}
+    \prod_{i=1}^N\left(
+      \sum_{j,k=1}^N
+      U_{i,j}^*\varphi_j^-
+      U_{i,k}\varphi_k^+
+    \right)
+    =
+    \sum_{\tau,\tau'\in\mathcal S_N}
+    (-1)^\tau(-1)^{\tau'}
+    \prod_{i=1}^NU^*_{i,\tau(i)},U_{i,\tau'(i)}\varphi_{i}^-\varphi_{i}^+
+    =|\det(U)|^2
+    \prod_{i=1}^N\varphi_i^-\varphi_i^+
+  \end{equation}
+  and, since $\det(U)=1$,
+  \begin{equation}
+    \int\prod_{i=1}^Nd\varphi_i^+d\varphi_i^-\ f(U^\dagger\varphi^-,U\varphi^+)
+    =
+    \alpha
+    =
+    \int \prod_{i=1}^Nd\psi_i^+d\psi_i^-\ f(\psi^-,\psi^+)
+    .
+  \end{equation}
+\qed
+\bigskip
+
+\subsection{Gaussian Grassmann integrals}
+\indent
+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$:
+  \begin{equation}
+    \int P_\mu(d\psi)\ 1
+    =
+    \det(\mu)\int \prod_{i=1}^Nd\varphi_i^+d\varphi_i^-\ e^{-\sum_{i=1}^N\lambda_i\varphi_i^+\varphi_i^-}
+  \end{equation}
+  where $\lambda_i$ are the eigenvalues of $\mu$.
+  Furthermore,
+  \begin{equation}
+    e^{-\sum_{i=1}^N\lambda_i\varphi_i^+\varphi_i^-}
+    =
+    \prod_{i=1}^Ne^{-\lambda_i\varphi_i^+\varphi_i^-}
+    =
+    \prod_{i=1}^N(1-\lambda_i\varphi_i^+\varphi_i^-)
+    \label{exp_trivial}
+  \end{equation}
+  so
+  \begin{equation}
+    \int P_\mu(d\psi)\ 1
+    =
+    \det(\mu)\int \prod_{i=1}^Nd\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
+    .
+  \end{equation}
+  \bigskip
+
+  \indent
+  Similarly,
+  \begin{equation}
+    \int P_\mu(d\psi)\ \psi_i^-\psi_j^+
+    =
+    \sum_{k,l=1}^N
+    \det(\mu)\int \prod_{i=1}^Nd\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^+
+  \end{equation}
+  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$:
+  \begin{equation}
+    \begin{largearray}
+      \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}^Nd\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{largearray}
+  \end{equation}
+  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{equation}
+    \begin{largearray}
+      \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}^Nd\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{largearray}
+  \end{equation}
+  and so
+  \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}\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)}}
+    .
+  \end{equation}
+\qed
+
+\subsection{Addition property}
+\indent
+Finally, let us prove the addition property of Gaussiann Grassmann integrals.
+\bigskip
+
+\theoname{Lemma}{Addition property}\label{lemma:grassmann_addition}
+  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$,
+  \nopagebreakaftereq
+  \begin{equation}
+    \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)
+    .
+  \end{equation}
+\endtheo
+\restorepagebreakaftereq
+\bigskip
+
+\indent\underline{Proof}:
+  It is sufficient to prove the lemma when $f$ is a monomial of the form
+  \begin{equation}
+    f(\psi)=
+    \prod_{i=1}^{n}
+    \psi_{j_i}^-\psi_{\bar j_i}^+
+    .
+  \end{equation}
+  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=\frac1{\sqrt2}(\psi_1+\psi_2)$ and $\varphi_2=\frac1{\sqrt2}(\psi_1-\psi_2)$ in terms of which
+  \begin{equation}
+    \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)
+  \end{equation}
+  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\-~\ref{lemma:wick_grassmann}, and thus, changing variables back to $\psi$,
+  \begin{equation}
+    \begin{largearray}
+      \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{largearray}
+  \end{equation}
+  Now,
+  \begin{equation}
+    \begin{largearray}
+      \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{largearray}
+  \end{equation}
+  so, by lemma\-~\ref{lemma:grassmann_id},
+  \begin{equation}
+    \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)}}
+  \end{equation}
+  and
+  \begin{equation}
+    \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)}}
+    .
+  \end{equation}
+  We conclude the proof of the lemma using the Wick rule for $\psi$:
+  \begin{equation}
+    \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)}}^+
+    .
+  \end{equation}
+\qed
+
+\vfill
+\eject
+
+
+\begin{thebibliography}{WWW99}
+\small
+\IfFileExists{bibliography/bibliography.tex}{\input bibliography/bibliography.tex}{}
+\end{thebibliography}
+
+
+\end{document}
diff --git a/Makefile b/Makefile
new file mode 100644
index 0000000..08cfaf2
--- /dev/null
+++ b/Makefile
@@ -0,0 +1,59 @@
+PROJECTNAME=$(basename $(wildcard *.tex))
+LIBS=$(notdir $(wildcard libs/*))
+FIGS=$(notdir $(wildcard figs/*.fig))
+
+PDFS=$(addsuffix .pdf, $(PROJECTNAME))
+SYNCTEXS=$(addsuffix .synctex.gz, $(PROJECTNAME))
+
+all: $(PROJECTNAME)
+
+$(PROJECTNAME): $(LIBS) $(FIGS)
+	pdflatex -file-line-error $@.tex
+	pdflatex -file-line-error $@.tex
+	pdflatex -synctex=1 $@.tex
+
+$(PROJECTNAME).aux: $(LIBS) $(FIGS)
+	pdflatex -file-line-error -draftmode $(PROJECTNAME).tex
+
+
+$(SYNCTEXS): $(LIBS) $(FIGS)
+	pdflatex -synctex=1 $(patsubst %.synctex.gz, %.tex, $@)
+
+
+libs: $(LIBS)
+
+$(LIBS):
+	ln -fs libs/$@ ./
+
+bibliography/bibliography.tex: $(PROJECTNAME).aux
+	BBlog -c bibliography/conf.BBlog -d $(BIBLIOGRAPHY) -b bibliography/bibliography.tex
+
+figs: $(FIGS)
+
+$(FIGS):
+	make -C figs/$@
+	for pdf in $$(find figs/$@/ -name '*.pdf'); do ln -fs "$$pdf" ./ ; done
+
+clean-aux: clean-figs-aux
+	rm -f $(addsuffix .aux, $(PROJECTNAME))
+	rm -f $(addsuffix .log, $(PROJECTNAME))
+	rm -f $(addsuffix .out, $(PROJECTNAME))
+	rm -f $(addsuffix .toc, $(PROJECTNAME))
+
+clean-libs:
+	rm -f $(LIBS)
+
+clean-figs:
+	$(foreach fig,$(addprefix figs/, $(FIGS)), make -C $(fig) clean clean-dat; )
+	rm -f $(notdir $(wildcard figs/*.fig/*.pdf))
+
+clean-figs-aux:
+	$(foreach fig,$(addprefix figs/, $(FIGS)), make -C $(fig) clean-aux; )
+
+clean-tex:
+	rm -f $(PDFS) $(SYNCTEXS)
+
+clean-bibliography:
+	rm -f bibliography/bibliography.tex
+
+clean: clean-aux clean-tex clean-libs clean-figs
diff --git a/README b/README
new file mode 100644
index 0000000..45215fa
--- /dev/null
+++ b/README
@@ -0,0 +1,42 @@
+This directory contains the source files to typeset the article, and generate
+the figures. This can be accomplished by running
+  make
+
+This document uses a custom class file, located in the 'libs' directory, which
+defines a number of commands. Most of these are drop-in replacements for those
+defined in the 'article' class.
+
+Some extra functionality is provided in custom style files, located in the
+'libs' directory.
+
+
+* Dependencies:
+
+  pdflatex
+  TeXlive packages:
+    amsfonts
+    array
+    color
+    graphics
+    hyperref
+    latex
+    marginnote
+    mathds
+    pgf
+    standalone
+  GNU make
+  gnuplot
+  meankondo v1.5
+
+
+* Files:
+
+  Jauslin_2022.tex:
+    main LaTeX file
+
+  libs:
+    custom LaTeX class file
+
+  figs:
+    source code for the figures
+
diff --git a/bibliography/bibliography.tex b/bibliography/bibliography.tex
new file mode 100644
index 0000000..8611cc5
--- /dev/null
+++ b/bibliography/bibliography.tex
@@ -0,0 +1,80 @@
+\bibitem[An70]{An70}P.W. Anderson - {\it A poor man's derivation of scaling laws for the Kondo problem}, Journal of Physics C: Solid State Physics, volume\-~3, page\-~2436, 1970,\par\penalty10000
+doi:{\tt\color{blue}\href{http://dx.doi.org/10.1088/0022-3719/3/12/008}{10.1088/0022-3719/3/12/008}}.\par\medskip
+ 
+\bibitem[An80]{An80}N. Andrei - {\it Diagonalization of the Kondo Hamiltonian}, Physical Review Letters, volume\-~45, issue\-~5, 1980,\par\penalty10000
+doi:{\tt\color{blue}\href{http://dx.doi.org/10.1103/PhysRevLett.45.379}{10.1103/PhysRevLett.45.379}}.\par\medskip
+ 
+\bibitem[BF84]{BF84}G.A. Battle, P. Federbush - {\it A note on cluster expansions, tree graph identities, extra $1/N!$ factors!!!}, Letters in Mathematical Physics, volume\-~8, pages\-~55-57, 1984,\par\penalty10000
+doi:{\tt\color{blue}\href{http://dx.doi.org/10.1007/BF00420041}{10.1007/BF00420041}}.\par\medskip
+ 
+\bibitem[BCe78]{BCe78}G. Benfatto, M. Cassandro, G. Gallavotti, F. Nicol\`o, E. Olivieri, E. Presutti, E. Scacciatelli - {\it Some probabilistic techniques in field theory}, Communications in Mathematical Physics, volume\-~59, issue\-~2, pages\-~143-166, 1978,\par\penalty10000
+doi:{\tt\color{blue}\href{http://dx.doi.org/10.1007/BF01614247}{10.1007/BF01614247}}.\par\medskip
+ 
+\bibitem[BG90]{BG90}G. Benfatto, G. Gallavotti - {\it 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,\par\penalty10000
+doi:{\tt\color{blue}\href{http://dx.doi.org/10.1007/BF01025844}{10.1007/BF01025844}}.\par\medskip
+ 
+\bibitem[BGe94]{BGe94}G. Benfatto, G. Gallavotti,  A.Procacci, B. Scoppola - {\it 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,\par\penalty10000
+doi:{\tt\color{blue}\href{http://dx.doi.org/10.1007/BF02099791}{10.1007/BF02099791}}.\par\medskip
+ 
+\bibitem[BGJ15]{BGJ15}G. Benfatto, G. Gallavotti, I. Jauslin - {\it Kondo Effect in a Fermionic Hierarchical Model}, Journal of Statistical Physics, volume\-~161, issue\-~5, pages\-~1203-1230, 2015,\par\penalty10000
+doi:{\tt\color{blue}\href{http://dx.doi.org/10.1007/s10955-015-1378-7}{10.1007/s10955-015-1378-7}}, arxiv:{\tt\color{blue}\href{http://arxiv.org/abs/1506.04381}{1506.04381}}.\par\medskip
+ 
+\bibitem[BF78]{BF78}D. Brydges, P. Federbush - {\it A new form of the Mayer expansion in classical statistical mechanics}, Journal of Mathematical Physics, volume\-~19, page\-~2064, 1978,\par\penalty10000
+doi:{\tt\color{blue}\href{http://dx.doi.org/10.1063/1.523586}{10.1063/1.523586}}.\par\medskip
+ 
+\bibitem[Do91]{Do91}T.C. Dorlas - {\it Renormalization group analysis of a simple hierarchical fermion model}, Communications in Mathematical Physics, volume\-~136, pages\-~169-194, 1991,\par\penalty10000
+doi:{\tt\color{blue}\href{http://dx.doi.org/10.1007/BF02096796}{10.1007/BF02096796}}.\par\medskip
+ 
+\bibitem[Dy69]{Dy69}F.J. Dyson - {\it Existence of a phase-transition in a one-dimensional Ising ferromagnet}, Communications in Mathematical Physics, volume\-~12, pages\-~91-107, 1969,\par\penalty10000
+doi:{\tt\color{blue}\href{http://dx.doi.org/10.1007/BF01645907}{10.1007/BF01645907}}.\par\medskip
+ 
+\bibitem[GJ15]{GJ15}G. Gallavotti, I. Jauslin - {\it Kondo Effect in the Hierarchical s-d Model}, Journal of Statistical Physics, volume\-~161, issue\-~5, pages\-~1231-1235, 2015,\par\penalty10000
+doi:{\tt\color{blue}\href{http://dx.doi.org/10.1007/s10955-015-1370-2}{10.1007/s10955-015-1370-2}}, arxiv:{\tt\color{blue}\href{http://arxiv.org/abs/1507.05678}{1507.05678}}.\par\medskip
+ 
+\bibitem[GN85]{GN85}G. Gallavotti, F. Nicol\`o - {\it Renormalization theory for four dimensional scalar fields I}, Communications in Mathematical Physics, volume\-~100, pages\-~545-590, 1985,\par\penalty10000
+doi:{\tt\color{blue}\href{http://dx.doi.org/10.1007/BF01217729}{10.1007/BF01217729}}.\par\medskip
+ 
+\bibitem[GK81]{GK81}K. Gawedzki, A. Kupiainen - {\it 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,\par\penalty10000
+doi:{\tt\color{blue}\href{http://dx.doi.org/10.1007/BF01237048}{10.1007/BF01237048}}.\par\medskip
+ 
+\bibitem[Gi10]{Gi10}A. Giuliani - {\it 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:{\tt\color{blue}\href{http://arxiv.org/abs/1102.3881}{1102.3881}}.\par\medskip
+ 
+\bibitem[GGM12]{GGM12}A. Giuliani, R.L. Greenblatt, V. Mastropietro - {\it The scaling limit of the energy correlations in non-integrable Ising models}, Journal of Mathematical Physics, volume\-~53, page\-~095214, 2012,\par\penalty10000
+doi:{\tt\color{blue}\href{http://dx.doi.org/10.1063/1.4745910}{10.1063/1.4745910}}, arxiv:{\tt\color{blue}\href{http://arxiv.org/abs/1204.4040}{1204.4040}}.\par\medskip
+ 
+\bibitem[GJ16]{GJ16}A. Giuliani, I. Jauslin - {\it The ground state construction of bilayer graphene}, Reviews in Mathematical Physics, volume\-~28, issue\-~8, page\-~1650018, 2016,\par\penalty10000
+doi:{\tt\color{blue}\href{http://dx.doi.org/10.1142/S0129055X16500185}{10.1142/S0129055X16500185}}, arxiv:{\tt\color{blue}\href{http://arxiv.org/abs/1507.06024}{1507.06024}}.\par\medskip
+ 
+\bibitem[GM10]{GM10}A. Giuliani, V. Mastropietro - {\it The two-dimensional Hubbard model on the honeycomb lattice}, Communications in Mathematical Physics, volume\-~293, pages\-~301-364, 2010,\par\penalty10000
+doi:{\tt\color{blue}\href{http://dx.doi.org/10.1007/s00220-009-0910-5}{10.1007/s00220-009-0910-5}}, arxiv:{\tt\color{blue}\href{http://arxiv.org/abs/0811.1881}{0811.1881}}.\par\medskip
+ 
+\bibitem[GMP12]{GMP12}A. Giuliani, V. Mastropietro, M. Porta - {\it Universality of conductivity in interacting graphene}, Communications in Mathematical Physics, volume\-~311, issue\-~2, pages\-~317-355, 2012,\par\penalty10000
+doi:{\tt\color{blue}\href{http://dx.doi.org/10.1007/s00220-012-1444-9}{10.1007/s00220-012-1444-9}}, arxiv:{\tt\color{blue}\href{http://arxiv.org/abs/1101.2169}{1101.2169}}.\par\medskip
+ 
+\bibitem[GMP17]{GMP17}A. Giuliani, V. Mastropietro, M. Porta - {\it Universality of the Hall Conductivity in Interacting Electron Systems}, Communications in Mathematical Physics, volume\-~349, issue\-~3, pages\-~1107-1161, 2017,\par\penalty10000
+doi:{\tt\color{blue}\href{http://dx.doi.org/10.1007/s00220-016-2714-8}{10.1007/s00220-016-2714-8}}, arxiv:{\tt\color{blue}\href{http://arxiv.org/abs/1511.04047}{1511.04047}}.\par\medskip
+ 
+\bibitem[Ko64]{Ko64}J. Kondo - {\it Resistance minimum in dilute magnetic alloys}, Progress of Theoretical Physics, volume\-~32, issue\-~1, 1964,\par\penalty10000
+doi:{\tt\color{blue}\href{http://dx.doi.org/10.1143/PTP.32.37}{10.1143/PTP.32.37}}.\par\medskip
+ 
+\bibitem[Ma11]{Ma11}V. Mastropietro - {\it Conductivity between Luttinger liquids: coupled chains and bilayer graphene}, Physical Review B, volume\-~84, page\-~035109, 2011,\par\penalty10000
+doi:{\tt\color{blue}\href{http://dx.doi.org/10.1103/PhysRevB.84.035109}{10.1103/PhysRevB.84.035109}}, arxiv:{\tt\color{blue}\href{http://arxiv.org/abs/1012.5736}{1012.5736}}.\par\medskip
+ 
+\bibitem[Mc57]{Mc57}J.W. McClure - {\it Band structure of graphite and de Haas-van Alphen effect}, Physical review, volume\-~108, pages\-~612-618, 1957,\par\penalty10000
+doi:{\tt\color{blue}\href{http://dx.doi.org/10.1103/PhysRev.108.612}{10.1103/PhysRev.108.612}}.\par\medskip
+ 
+\bibitem[NGe04]{NGe04}K.S. Novoselov, A.K. Geim, S. V.Morozov, D. Jiang, Y. Zhang, S.V. Dubonos, I.V. Grigorieva, A.A. Firsov - {\it Electric field effect in atomically thin carbon films}, Science, vol. 306, pages\-~666-669, 2004,\par\penalty10000
+doi:{\tt\color{blue}\href{http://dx.doi.org/10.1126/science.1102896}{10.1126/science.1102896}}.\par\medskip
+ 
+\bibitem[SW58]{SW58}J.C. Slonczewski, P.R. Weiss - {\it Band structure of graphite}, Physical Review, volume\-~109, pages\-~272-279, 1958,\par\penalty10000
+doi:{\tt\color{blue}\href{http://dx.doi.org/10.1103/PhysRev.109.272}{10.1103/PhysRev.109.272}}.\par\medskip
+ 
+\bibitem[Wi65]{Wi65}K.G. Wilson - {\it Model Hamiltonians for Local Quantum Field Theory}, Physical Review, volume\-~140, issue\-~2B, page B445-B457, 1965,\par\penalty10000
+doi:{\tt\color{blue}\href{http://dx.doi.org/10.1103/PhysRev.140.B445}{10.1103/PhysRev.140.B445}}.\par\medskip
+ 
+\bibitem[Wi75]{Wi75}K.G. Wilson - {\it The renormalization group: Critical phenomena and the Kondo problem}, Reviews of Modern Physics, volume\-~47, issue\-~4, 1975,\par\penalty10000
+doi:{\tt\color{blue}\href{http://dx.doi.org/10.1103/RevModPhys.47.773}{10.1103/RevModPhys.47.773}}.\par\medskip
+ 
+\bibitem[mk]{mk}{\tt meankondo} software package v1.5, I.\-~Jauslin,\par\penalty10000
+{\tt\color{blue}\href{http://ian.jauslin.org/software/meankondo}{http://ian.jauslin.org/software/meankondo}}.\par\medskip
+ 
diff --git a/figs/bands.fig/Makefile b/figs/bands.fig/Makefile
new file mode 100644
index 0000000..cb5ada5
--- /dev/null
+++ b/figs/bands.fig/Makefile
@@ -0,0 +1,25 @@
+PROJECTNAME=$(basename $(wildcard *.gnuplot))
+
+PDFS=$(addsuffix .pdf, $(PROJECTNAME))
+
+all: $(PDFS)
+
+$(PDFS):
+	gnuplot $(patsubst %.pdf, %.gnuplot, $@) > $(patsubst %.pdf, %.tikz.tex, $@)
+	pdflatex -jobname $(basename $@) -file-line-error $(patsubst %.pdf, %.tikz.tex, $@)
+
+install: $(PDFS)
+	cp $^ $(INSTALLDIR)/
+
+clean-aux:
+	rm -f $(addsuffix .tikz.tex, $(PROJECTNAME))
+	rm -f $(addsuffix .aux, $(PROJECTNAME))
+	rm -f $(addsuffix .log, $(PROJECTNAME))
+
+clean-dat:
+	rm -f $(DATS)
+
+clean-tex:
+	rm -f $(PDFS)
+
+clean: clean-aux clean-tex
diff --git a/figs/bands.fig/bands.gnuplot b/figs/bands.fig/bands.gnuplot
new file mode 100644
index 0000000..8922c00
--- /dev/null
+++ b/figs/bands.fig/bands.gnuplot
@@ -0,0 +1,33 @@
+# default output canvas size: 12.5cm x 8.75cm
+set term lua tikz size 8,6 standalone
+
+unset key
+
+unset xtics
+unset ytics
+unset ztics
+unset border
+
+unset colorbox
+
+set pointsize 1
+
+#set xrange [0:4*pi/3]
+#set yrange [-2*pi/sqrt(3):2*pi/sqrt(3)]
+
+set pm3d depthorder
+
+# make this odd and divisible by 7 to get the Fermi points
+set isosamples 49
+
+# rescale x and y components to render the Fermi points better
+# units of 2*pi/3
+set xrange [0:2]
+# units of 2*pi/(3*sqrt(3))
+set yrange [-3:3]
+
+set view 80,35
+
+splot sqrt(1+4*cos(pi*x)*cos(pi/3*y)+4*cos(pi/3*y)**2) w pm3d, -sqrt(1+4*cos(pi*x)*cos(pi/3*y)+4*cos(pi/3*y)**2) w pm3d
+
+#pause mouse close
diff --git a/figs/bilayer.fig/Makefile b/figs/bilayer.fig/Makefile
new file mode 100644
index 0000000..147037c
--- /dev/null
+++ b/figs/bilayer.fig/Makefile
@@ -0,0 +1,23 @@
+PROJECTNAME=$(basename $(basename $(wildcard *.tikz.tex)))
+
+PDFS=$(addsuffix .pdf, $(PROJECTNAME))
+
+all: $(PDFS)
+
+$(PDFS):
+	pdflatex -jobname $(basename $@) $(patsubst %.pdf, %.tikz.tex, $@)
+
+
+install: $(PROJECTNAME)
+	cp $^.pdf $(INSTALLDIR)/
+
+
+clean-aux:
+	rm -f $(addsuffix .aux, $(PROJECTNAME))
+	rm -f $(addsuffix .log, $(PROJECTNAME))
+	rm -f $(addsuffix .out, $(PROJECTNAME))
+
+clean-tex:
+	rm -f $(PDFS) $(SYNCTEXS)
+
+clean: clean-aux clean-tex
diff --git a/figs/bilayer.fig/bilayer.tikz.tex b/figs/bilayer.fig/bilayer.tikz.tex
new file mode 100644
index 0000000..a35b039
--- /dev/null
+++ b/figs/bilayer.fig/bilayer.tikz.tex
@@ -0,0 +1,40 @@
+\documentclass{standalone}
+
+\usepackage{tikz}
+\usepackage{graphene}
+
+% reflected graphene grid at #1 of width #2 and height #3
+\def\graphenereflected#1#2#3{
+  \foreach \i in {0,...,#2}{
+    \foreach \j in {0,...,#3}{
+      \cellreflected{#1++(\i*3,-2*\j*\sqrtThOT)}
+      \cellreflected{#1++(\i*3+1.5,\sqrtThOT-2*\j*\sqrtThOT)}
+    }
+  }
+}
+\def\cellreflected#1{
+  \draw[dotted]#1--++(0:-1);
+  \draw[dotted]#1--++(120:-1);
+  \draw[dotted]#1--++(240:-1);
+  \asite{#1}
+  \draw[dotted]#1++(-1,0)--++(0:1);
+  \draw[dotted]#1++(-1,0)--++(120:1);
+  \draw[dotted]#1++(-1,0)--++(240:1);
+  \bsite{#1++(-1,0)}
+}
+\def\square#1#2{\draw#1++(-#2,-#2)--++(#2,0)--++(#2,0)--++(0,#2)--++(0,#2)--++(-#2,0)--++(-#2,0)--++(0,-#2)--++(0,-#2);}
+
+
+\begin{document}
+\begin{tikzpicture}
+  % first layer
+  \graphene{(0,0)}{3}{3}
+
+  % redefine sites for second layer
+  \def\bsite#1{\draw#1circle(.1);}
+  \def\asite#1{\square{#1}{.17}}
+
+  % second layer
+  \graphenereflected{(0,0)}{3}{3}
+\end{tikzpicture}
+\end{document}
diff --git a/figs/bilayer.fig/graphene.sty b/figs/bilayer.fig/graphene.sty
new file mode 120000
index 0000000..477bbc1
--- /dev/null
+++ b/figs/bilayer.fig/graphene.sty
@@ -0,0 +1 @@
+../lib/graphene.sty
\ No newline at end of file
diff --git a/figs/brillouin.fig/Makefile b/figs/brillouin.fig/Makefile
new file mode 100644
index 0000000..147037c
--- /dev/null
+++ b/figs/brillouin.fig/Makefile
@@ -0,0 +1,23 @@
+PROJECTNAME=$(basename $(basename $(wildcard *.tikz.tex)))
+
+PDFS=$(addsuffix .pdf, $(PROJECTNAME))
+
+all: $(PDFS)
+
+$(PDFS):
+	pdflatex -jobname $(basename $@) $(patsubst %.pdf, %.tikz.tex, $@)
+
+
+install: $(PROJECTNAME)
+	cp $^.pdf $(INSTALLDIR)/
+
+
+clean-aux:
+	rm -f $(addsuffix .aux, $(PROJECTNAME))
+	rm -f $(addsuffix .log, $(PROJECTNAME))
+	rm -f $(addsuffix .out, $(PROJECTNAME))
+
+clean-tex:
+	rm -f $(PDFS) $(SYNCTEXS)
+
+clean: clean-aux clean-tex
diff --git a/figs/brillouin.fig/brillouin.tikz.tex b/figs/brillouin.fig/brillouin.tikz.tex
new file mode 100644
index 0000000..6b392af
--- /dev/null
+++ b/figs/brillouin.fig/brillouin.tikz.tex
@@ -0,0 +1,17 @@
+\documentclass{standalone}
+
+\usepackage{tikz}
+
+% sqrt(3)
+\def\sqrtTh{1.732}
+% pi
+\def\pival{3.14159}
+
+\begin{document}
+\begin{tikzpicture}
+  \draw(0,0)--(2*\pival/3,2*\pival/\sqrtTh)--(4*\pival/3,0)--(2*\pival/3,-2*\pival/\sqrtTh)--cycle;
+
+  \fill(2*\pival/3,2*\pival/3/\sqrtTh)circle(0.05);
+  \fill(2*\pival/3,-2*\pival/3/\sqrtTh)circle(0.05);
+\end{tikzpicture}
+\end{document}
diff --git a/figs/brillouin.fig/scales.tikz.tex b/figs/brillouin.fig/scales.tikz.tex
new file mode 100644
index 0000000..ebafd14
--- /dev/null
+++ b/figs/brillouin.fig/scales.tikz.tex
@@ -0,0 +1,47 @@
+\documentclass{standalone}
+
+\usepackage{tikz}
+\usepackage{xcolor}
+
+\definecolor{gray0}{HTML}{EEEEEE}
+\definecolor{gray1}{HTML}{CCCCCC}
+\definecolor{gray2}{HTML}{AAAAAA}
+\definecolor{gray3}{HTML}{888888}
+\definecolor{gray4}{HTML}{666666}
+\definecolor{gray5}{HTML}{444444}
+\definecolor{gray6}{HTML}{222222}
+\definecolor{gray7}{HTML}{000000}
+
+% sqrt(3)
+\def\sqrtTh{1.732}
+% pi
+\def\pival{3.14159}
+
+\begin{document}
+\begin{tikzpicture}
+  \draw(0,0)--(2*\pival/3,2*\pival/\sqrtTh)--(4*\pival/3,0)--(2*\pival/3,-2*\pival/\sqrtTh)--cycle;
+
+  \fill[color=gray0](0,0)--(2*\pival/3,2*\pival/\sqrtTh)--(4*\pival/3,0)--(2*\pival/3,-2*\pival/\sqrtTh)--cycle;
+
+  \fill[color=gray1](2*\pival/3, 2*\pival/3/\sqrtTh)circle(1);
+  \fill[color=gray1](2*\pival/3,-2*\pival/3/\sqrtTh)circle(1);
+
+  \fill[color=gray2](2*\pival/3, 2*\pival/3/\sqrtTh)circle(0.5);
+  \fill[color=gray2](2*\pival/3,-2*\pival/3/\sqrtTh)circle(0.5);
+
+  \fill[color=gray3](2*\pival/3, 2*\pival/3/\sqrtTh)circle(0.25);
+  \fill[color=gray3](2*\pival/3,-2*\pival/3/\sqrtTh)circle(0.25);
+
+  \fill[color=gray4](2*\pival/3, 2*\pival/3/\sqrtTh)circle(0.125);
+  \fill[color=gray4](2*\pival/3,-2*\pival/3/\sqrtTh)circle(0.125);
+
+  \fill[color=gray5](2*\pival/3, 2*\pival/3/\sqrtTh)circle(0.0626);
+  \fill[color=gray5](2*\pival/3,-2*\pival/3/\sqrtTh)circle(0.0625);
+
+  \fill[color=gray6](2*\pival/3, 2*\pival/3/\sqrtTh)circle(0.03125);
+  \fill[color=gray6](2*\pival/3,-2*\pival/3/\sqrtTh)circle(0.03125);
+
+  \fill[color=gray7](2*\pival/3, 2*\pival/3/\sqrtTh)circle(0.015625);
+  \fill[color=gray7](2*\pival/3,-2*\pival/3/\sqrtTh)circle(0.015625);
+\end{tikzpicture}
+\end{document}
diff --git a/figs/cell-coarse.fig/Makefile b/figs/cell-coarse.fig/Makefile
new file mode 100644
index 0000000..147037c
--- /dev/null
+++ b/figs/cell-coarse.fig/Makefile
@@ -0,0 +1,23 @@
+PROJECTNAME=$(basename $(basename $(wildcard *.tikz.tex)))
+
+PDFS=$(addsuffix .pdf, $(PROJECTNAME))
+
+all: $(PDFS)
+
+$(PDFS):
+	pdflatex -jobname $(basename $@) $(patsubst %.pdf, %.tikz.tex, $@)
+
+
+install: $(PROJECTNAME)
+	cp $^.pdf $(INSTALLDIR)/
+
+
+clean-aux:
+	rm -f $(addsuffix .aux, $(PROJECTNAME))
+	rm -f $(addsuffix .log, $(PROJECTNAME))
+	rm -f $(addsuffix .out, $(PROJECTNAME))
+
+clean-tex:
+	rm -f $(PDFS) $(SYNCTEXS)
+
+clean: clean-aux clean-tex
diff --git a/figs/cell-coarse.fig/cell-coarse.tikz.tex b/figs/cell-coarse.fig/cell-coarse.tikz.tex
new file mode 100644
index 0000000..f296cb5
--- /dev/null
+++ b/figs/cell-coarse.fig/cell-coarse.tikz.tex
@@ -0,0 +1,28 @@
+\documentclass{standalone}
+
+\usepackage{tikz}
+\usepackage{graphene}
+
+% rescale cells by
+\def\rr{4}
+
+\begin{document}
+\begin{tikzpicture}
+  %\graphene{(0,0)}{12}{12}
+  \graphene{(-3,3*\sqrtTh)}{17}{16}
+
+  \foreach \i in {0,...,3}{
+    \draw[line width=1pt, color=blue](-1*\rr,-\i*\sqrtTh*\rr)--++(1.5*2*\i*\rr+3*\rr,\sqrtThOT*2*\i*\rr+\sqrtTh*\rr);
+  }
+  \foreach \i in {0,...,3}{
+    \draw[line width=1pt, color=blue](0.5*\rr+3*\i*\rr,-3.5*\sqrtTh*\rr)--++(-1.5*2*\i*\rr+12*\rr,-\sqrtThOT*2*\i*\rr+4*\sqrtTh*\rr);
+  }
+
+  \foreach \i in {0,...,3}{
+    \draw[line width=1pt, color=blue](-1*\rr,-\i*\sqrtTh*\rr)--++(-1.5*2*\i*\rr+10.5*\rr,+\sqrtThOT*2*\i*\rr-3.5*\sqrtTh*\rr);
+  }
+  \foreach \i in {0,...,3}{
+    \draw[line width=1pt, color=blue](-1*\rr+3*\i*\rr,\sqrtTh*\rr)--++(-1.5*2*\i*\rr+13.5*\rr,\sqrtThOT*2*\i*\rr-4.5*\sqrtTh*\rr);
+  }
+\end{tikzpicture}
+\end{document}
diff --git a/figs/cell-coarse.fig/graphene.sty b/figs/cell-coarse.fig/graphene.sty
new file mode 120000
index 0000000..477bbc1
--- /dev/null
+++ b/figs/cell-coarse.fig/graphene.sty
@@ -0,0 +1 @@
+../lib/graphene.sty
\ No newline at end of file
diff --git a/figs/cell.fig/Makefile b/figs/cell.fig/Makefile
new file mode 100644
index 0000000..147037c
--- /dev/null
+++ b/figs/cell.fig/Makefile
@@ -0,0 +1,23 @@
+PROJECTNAME=$(basename $(basename $(wildcard *.tikz.tex)))
+
+PDFS=$(addsuffix .pdf, $(PROJECTNAME))
+
+all: $(PDFS)
+
+$(PDFS):
+	pdflatex -jobname $(basename $@) $(patsubst %.pdf, %.tikz.tex, $@)
+
+
+install: $(PROJECTNAME)
+	cp $^.pdf $(INSTALLDIR)/
+
+
+clean-aux:
+	rm -f $(addsuffix .aux, $(PROJECTNAME))
+	rm -f $(addsuffix .log, $(PROJECTNAME))
+	rm -f $(addsuffix .out, $(PROJECTNAME))
+
+clean-tex:
+	rm -f $(PDFS) $(SYNCTEXS)
+
+clean: clean-aux clean-tex
diff --git a/figs/cell.fig/cell.tikz.tex b/figs/cell.fig/cell.tikz.tex
new file mode 100644
index 0000000..ec66a7c
--- /dev/null
+++ b/figs/cell.fig/cell.tikz.tex
@@ -0,0 +1,30 @@
+\documentclass{standalone}
+
+\usepackage{tikz}
+\usepackage{graphene}
+
+\begin{document}
+\begin{tikzpicture}
+  \graphene{(0,0)}{3}{3}
+
+  \foreach \i in {0,...,3}{
+    \draw[color=blue](-1,-\i*\sqrtTh)--++(1.5*2*\i+3,\sqrtThOT*2*\i+\sqrtTh);
+  }
+  \foreach \i in {0,...,3}{
+    \draw[color=blue](0.5+3*\i,-3.5*\sqrtTh)--++(-1.5*2*\i+12,-\sqrtThOT*2*\i+4*\sqrtTh);
+  }
+
+  \foreach \i in {0,...,3}{
+    \draw[color=blue](-1,-\i*\sqrtTh)--++(-1.5*2*\i+10.5,+\sqrtThOT*2*\i-3.5*\sqrtTh);
+  }
+  \foreach \i in {0,...,3}{
+    \draw[color=blue](-1+3*\i,\sqrtTh)--++(-1.5*2*\i+13.5,\sqrtThOT*2*\i-4.5*\sqrtTh);
+  }
+
+  \draw[color=red] (6,-\sqrtTh)++(0.5,.6)node[rotate=30]{$l_1$};
+  \draw[color=red](6,-\sqrtTh)++(0.5,-.6)node[rotate=-30]{$l_2$};
+  \pgfsetarrowsend{latex};
+  \draw[color=red,line width=1.5pt](6,-\sqrtTh)--++\Lo;
+  \draw[color=red,line width=1.5pt](6,-\sqrtTh)--++\Lt;
+\end{tikzpicture}
+\end{document}
diff --git a/figs/cell.fig/graphene.sty b/figs/cell.fig/graphene.sty
new file mode 120000
index 0000000..477bbc1
--- /dev/null
+++ b/figs/cell.fig/graphene.sty
@@ -0,0 +1 @@
+../lib/graphene.sty
\ No newline at end of file
diff --git a/figs/graphene.fig/Makefile b/figs/graphene.fig/Makefile
new file mode 100644
index 0000000..147037c
--- /dev/null
+++ b/figs/graphene.fig/Makefile
@@ -0,0 +1,23 @@
+PROJECTNAME=$(basename $(basename $(wildcard *.tikz.tex)))
+
+PDFS=$(addsuffix .pdf, $(PROJECTNAME))
+
+all: $(PDFS)
+
+$(PDFS):
+	pdflatex -jobname $(basename $@) $(patsubst %.pdf, %.tikz.tex, $@)
+
+
+install: $(PROJECTNAME)
+	cp $^.pdf $(INSTALLDIR)/
+
+
+clean-aux:
+	rm -f $(addsuffix .aux, $(PROJECTNAME))
+	rm -f $(addsuffix .log, $(PROJECTNAME))
+	rm -f $(addsuffix .out, $(PROJECTNAME))
+
+clean-tex:
+	rm -f $(PDFS) $(SYNCTEXS)
+
+clean: clean-aux clean-tex
diff --git a/figs/graphene.fig/graphene.sty b/figs/graphene.fig/graphene.sty
new file mode 120000
index 0000000..477bbc1
--- /dev/null
+++ b/figs/graphene.fig/graphene.sty
@@ -0,0 +1 @@
+../lib/graphene.sty
\ No newline at end of file
diff --git a/figs/graphene.fig/graphene.tikz.tex b/figs/graphene.fig/graphene.tikz.tex
new file mode 100644
index 0000000..c78e3e9
--- /dev/null
+++ b/figs/graphene.fig/graphene.tikz.tex
@@ -0,0 +1,12 @@
+\documentclass{standalone}
+
+\usepackage{tikz}
+\usepackage{graphene}
+
+\def\bsite#1{\fill#1circle(.1);}
+
+\begin{document}
+\begin{tikzpicture}
+  \graphene{(0,0)}{3}{3}
+\end{tikzpicture}
+\end{document}
diff --git a/figs/hierarchical_graphene.fig/Makefile b/figs/hierarchical_graphene.fig/Makefile
new file mode 100644
index 0000000..e5e1f5d
--- /dev/null
+++ b/figs/hierarchical_graphene.fig/Makefile
@@ -0,0 +1,32 @@
+PROJECTNAME=graphene_vector_field
+
+DATS=graphene_vector_field.dat
+PDFS=$(addsuffix .pdf, $(PROJECTNAME))
+TEXS=$(addsuffix .tikz.tex, $(PROJECTNAME))
+
+all: $(PDFS)
+
+$(PDFS): $(DATS)
+	gnuplot $(patsubst %.pdf, %.gnuplot, $@) > $(patsubst %.pdf, %.tikz.tex, $@)
+	pdflatex -jobname $(basename $@) -file-line-error $(patsubst %.pdf, %.tikz.tex, $@)
+
+graphene_vector_field.dat:
+	meankondo -C graphene.mk > graphene-numkondo.mk
+	./gendat > $@
+
+install: $(PDFS)
+	cp $^ $(INSTALLDIR)/
+
+clean-aux:
+	rm -f $(addsuffix .tikz.tex, $(PROJECTNAME))
+	rm -f $(addsuffix .aux, $(PROJECTNAME))
+	rm -f $(addsuffix .log, $(PROJECTNAME))
+
+clean-dat:
+	rm -f $(DATS)
+	rm -f graphene-numkondo.mk
+
+clean-tex:
+	rm -f $(PDFS)
+
+clean: clean-aux clean-tex
diff --git a/figs/hierarchical_graphene.fig/gendat b/figs/hierarchical_graphene.fig/gendat
new file mode 100755
index 0000000..d905203
--- /dev/null
+++ b/figs/hierarchical_graphene.fig/gendat
@@ -0,0 +1,30 @@
+#!/bin/bash
+
+# location of numkondo config
+config=$(dirname $0)/graphene-numkondo.mk
+
+# range and number of points in each direction
+minx=-2
+maxx=2
+nx=20
+
+miny=-1
+maxy=1
+ny=20
+
+# compute each point
+for i in $(seq $nx); do
+  # alpha_0
+  alpha0=$(python -c "print($minx+($maxx-($minx))/($nx-1)*($i-1))")
+  for j in $(seq $ny); do
+    # alpha_1
+    alpha1=$(python -c "print($miny+($maxy-($miny))/($ny-1)*($j-1))")
+
+    # print alphas
+    echo -n "$alpha0 $alpha1 "
+
+    # compute step and select the relevant lines
+    dat=$(numkondo -N 1 -F -I "0:$alpha0,1:$alpha1" "$config" | sed -r '/^[^01]/d;s/^[0-9]*://;s/,//' | tr '\n' ' ')
+    echo "$dat"
+  done
+done
diff --git a/figs/hierarchical_graphene.fig/graphene.mk b/figs/hierarchical_graphene.fig/graphene.mk
new file mode 100644
index 0000000..7bd1cca
--- /dev/null
+++ b/figs/hierarchical_graphene.fig/graphene.mk
@@ -0,0 +1,148 @@
+## This computes the logarithm of the flow equation
+
+#!fields
+# extrenal fields: xy0 where x=alpha, y=sigma
+x:110,120,210,220
+# internal fields: xy1 where x=alpha, y=sigma
+i:111,121,211,221
+# all are Fermions
+f:110,120,210,220,111,121,211,221
+
+&
+#!groups
+# different spins are independent
+(111,211) (121,221)
+
+&
+#!preprocessor_variables
+
+# psi's with internal fields
+psiAU=(1/2)[f110]+[f111],
+psiAD=(1/2)[f120]+[f121],
+psiBU=(1/2)[f210]+[f211],
+psiBD=(1/2)[f220]+[f221],
+
+psiAU-=(1/2)[f-110]+[f-111],
+psiAD-=(1/2)[f-120]+[f-121],
+psiBU-=(1/2)[f-210]+[f-211],
+psiBD-=(1/2)[f-220]+[f-221],
+
+# operators with internal fields
+O0=
+<<$psiAU>*<$psiBU->>+
+<<$psiBU>*<$psiAU->>+
+<<$psiAD>*<$psiBD->>+
+<<$psiBD>*<$psiAD->>,
+
+O1=
+<<$psiAU>*<$psiAU->*<$psiAD>*<$psiAD->>+
+<<$psiBU>*<$psiBU->*<$psiBD>*<$psiBD->>,
+
+O2=
+<<$psiAU>*<$psiAD->*<$psiBD>*<$psiBU->>+
+<<$psiBU>*<$psiBD->*<$psiAD>*<$psiAU->>+
+<<$psiAD>*<$psiAU->*<$psiBU>*<$psiBD->>+
+<<$psiBD>*<$psiBU->*<$psiAU>*<$psiAD->>,
+
+O3=
+<<$psiAU>*<$psiAU->*<$psiBU>*<$psiBU->>+
+<<$psiAD>*<$psiAD->*<$psiBD>*<$psiBD->>,
+
+O4=
+<<$psiAU>*<$psiBU->*<$psiAD>*<$psiBD->>+
+<<$psiBU>*<$psiAU->*<$psiBD>*<$psiAD->>,
+
+O5=
+<<$psiAU>*<$psiAU->*<$psiAD>*<$psiBU->*<$psiBU>*<$psiBD->>+
+<<$psiAD>*<$psiAD->*<$psiAU>*<$psiBD->*<$psiBD>*<$psiBU->>+
+<<$psiBU>*<$psiBU->*<$psiBD>*<$psiAU->*<$psiAU>*<$psiAD->>+
+<<$psiBD>*<$psiBD->*<$psiBU>*<$psiAD->*<$psiAD>*<$psiAU->>,
+
+O6=<<$psiAU>*<$psiAU->*<$psiAD>*<$psiAD->*<$psiBU>*<$psiBU->*<$psiBD>*<$psiBD->>,
+
+
+# psi's without internal fields
+phiAU=[f110],
+phiAD=[f120],
+phiBU=[f210],
+phiBD=[f220],
+
+phiAU-=[f-110],
+phiAD-=[f-120],
+phiBU-=[f-210],
+phiBD-=[f-220],
+
+# operators without internal fields
+E0=
+<<$phiAU>*<$phiBU->>+
+<<$phiBU>*<$phiAU->>+
+<<$phiAD>*<$phiBD->>+
+<<$phiBD>*<$phiAD->>,
+
+E1=
+<<$phiAU>*<$phiAU->*<$phiAD>*<$phiAD->>+
+<<$phiBU>*<$phiBU->*<$phiBD>*<$phiBD->>,
+
+E2=
+<<$phiAU>*<$phiAD->*<$phiBD>*<$phiBU->>+
+<<$phiBU>*<$phiBD->*<$phiAD>*<$phiAU->>+
+<<$phiAD>*<$phiAU->*<$phiBU>*<$phiBD->>+
+<<$phiBD>*<$phiBU->*<$phiAU>*<$phiAD->>,
+
+E3=
+<<$phiAU>*<$phiAU->*<$phiBU>*<$phiBU->>+
+<<$phiAD>*<$phiAD->*<$phiBD>*<$phiBD->>,
+
+E4=
+<<$phiAU>*<$phiBU->*<$phiAD>*<$phiBD->>+
+<<$phiBU>*<$phiAU->*<$phiBD>*<$phiAD->>,
+
+E5=
+<<$phiAU>*<$phiAU->*<$phiAD>*<$phiBU->*<$phiBU>*<$phiBD->>+
+<<$phiAD>*<$phiAD->*<$phiAU>*<$phiBD->*<$phiBD>*<$phiBU->>+
+<<$phiBU>*<$phiBU->*<$phiBD>*<$phiAU->*<$phiAU>*<$phiAD->>+
+<<$phiBD>*<$phiBD->*<$phiBU>*<$phiAD->*<$phiAD>*<$phiAU->>,
+
+E6=<<$phiAU>*<$phiAU->*<$phiAD>*<$phiAD->*<$phiBU>*<$phiBU->*<$phiBD>*<$phiBD->>
+
+&
+#!propagator
+111;211: 1 ,121;221: 1 ,
+211;111: 1 ,221;121: 1
+
+&
+#!input_polynomial
+<%exp<
+  <<[l0]>*<$O0>>+
+  <<[l1]>*<$O1>>+
+  <<[l2]>*<$O2>>+
+  <<[l3]>*<$O3>>+
+  <<[l4]>*<$O4>>+
+  <<[l5]>*<$O5>>+
+  <<[l6]>*<$O6>>
+>>
+
+&
+#!postprocess_flow_equation
+<<8>*<%log_1<$FLOW>>>
+
+&
+#!id_table
+0: <$E0>,
+1: <$E1>,
+2: <$E2>,
+3: <$E3>,
+4: <$E4>,
+5: <$E5>,
+6: <$E6>
+
+&
+#!labels
+
+0:"O0" ,
+1:"O1" ,
+2:"O2" ,
+3:"O3" ,
+4:"O4" ,
+5:"O5" ,
+6:"O6"
diff --git a/figs/hierarchical_graphene.fig/graphene_vector_field.gnuplot b/figs/hierarchical_graphene.fig/graphene_vector_field.gnuplot
new file mode 100644
index 0000000..9eef976
--- /dev/null
+++ b/figs/hierarchical_graphene.fig/graphene_vector_field.gnuplot
@@ -0,0 +1,25 @@
+set ylabel "$\\alpha_1$" norotate
+set xlabel "$\\alpha_0$"
+
+# default output canvas size: 12.5cm x 8.75cm
+set term lua tikz size 8,6 standalone
+
+unset key
+
+set pointsize 1
+
+# color functions
+
+rgb(r,g,b) = 65536 * int(r) + 256 * int(g) + int(b)
+color(x) = rgb(x*255, 0 , (1-x)*255)
+
+# get min/max
+set yrange [-1000:1000]
+stats "graphene_vector_field.dat" u (log10(sqrt(($3-$1)**2+($4-$2)**2)))
+
+set xrange [-2:2]
+set yrange [-1:1]
+
+rescale=15
+
+plot "graphene_vector_field.dat" u 1:2:(($3-$1)/sqrt(($3-$1)**2+($4-$2)**2)/rescale):(($4-$2)/sqrt(($3-$1)**2+($4-$2)**2)/rescale):(color((log10(sqrt(($3-$1)**2+($4-$2)**2))-STATS_min)/(STATS_max-STATS_min))) with vec filled head linecolor rgb variable
diff --git a/figs/hierarchical_sd.fig/Makefile b/figs/hierarchical_sd.fig/Makefile
new file mode 100644
index 0000000..175f61b
--- /dev/null
+++ b/figs/hierarchical_sd.fig/Makefile
@@ -0,0 +1,32 @@
+PROJECTNAME=sd_vector_field
+
+DATS=sd_vector_field.dat
+PDFS=$(addsuffix .pdf, $(PROJECTNAME))
+TEXS=$(addsuffix .tikz.tex, $(PROJECTNAME))
+
+all: $(PDFS)
+
+$(PDFS): $(DATS)
+	gnuplot $(patsubst %.pdf, %.gnuplot, $@) > $(patsubst %.pdf, %.tikz.tex, $@)
+	pdflatex -jobname $(basename $@) -file-line-error $(patsubst %.pdf, %.tikz.tex, $@)
+
+sd_vector_field.dat:
+	meankondo -C sd.mk > sd-numkondo.mk
+	./gendat > $@
+
+install: $(PDFS)
+	cp $^ $(INSTALLDIR)/
+
+clean-aux:
+	rm -f $(addsuffix .tikz.tex, $(PROJECTNAME))
+	rm -f $(addsuffix .aux, $(PROJECTNAME))
+	rm -f $(addsuffix .log, $(PROJECTNAME))
+
+clean-dat:
+	rm -f $(DATS)
+	rm -f sd-numkondo.mk
+
+clean-tex:
+	rm -f $(PDFS)
+
+clean: clean-aux clean-tex
diff --git a/figs/hierarchical_sd.fig/gendat b/figs/hierarchical_sd.fig/gendat
new file mode 100755
index 0000000..08bbdd9
--- /dev/null
+++ b/figs/hierarchical_sd.fig/gendat
@@ -0,0 +1,30 @@
+#!/bin/bash
+
+# location of numkondo config
+config=$(dirname $0)/sd-numkondo.mk
+
+# range and number of points in each direction
+minx=-1
+maxx=0.5
+nx=20
+
+miny=-0.4
+maxy=0.4
+ny=20
+
+# compute each point
+for i in $(seq $nx); do
+  # alpha_0
+  alpha0=$(python -c "print($minx+($maxx-($minx))/($nx-1)*($i-1))")
+  for j in $(seq $ny); do
+    # alpha_1
+    alpha1=$(python -c "print($miny+($maxy-($miny))/($ny-1)*($j-1))")
+
+    # print alphas
+    echo -n "$alpha0 $alpha1 "
+
+    # compute step and select the relevant lines
+    dat=$(numkondo -N 1 -F -I "0:$alpha0,1:$alpha1" "$config" | sed -r '/^[^01]/d;s/^[0-9]*://;s/,//' | tr '\n' ' ')
+    echo "$dat"
+  done
+done
diff --git a/figs/hierarchical_sd.fig/sd.mk b/figs/hierarchical_sd.fig/sd.mk
new file mode 100644
index 0000000..f0bd584
--- /dev/null
+++ b/figs/hierarchical_sd.fig/sd.mk
@@ -0,0 +1,81 @@
+#!fields
+x:100,110
+i:101,102,111,112
+h:400,410,420
+f:100,110,101,102,111,112
+a:400,410,420
+
+&
+#!preprocessor_variables
+psiAu1=[f100]+[f101],
+psiAu2=[f100]+[f102],
+psiAd1=[f110]+[f111],
+psiAd2=[f110]+[f112],
+
+psiAu1-=[f-100]+[f-101],
+psiAu2-=[f-100]+[f-102],
+psiAd1-=[f-110]+[f-111],
+psiAd2-=[f-110]+[f-112],
+
+A11=<<$psiAu1>*<$psiAd1->>+<<$psiAd1>*<$psiAu1->>,
+A12=<<$psiAu2>*<$psiAd2->>+<<$psiAd2>*<$psiAu2->>,
+
+A21=<<(-1)s{-1}>*<$psiAu1>*<$psiAd1->>+<<s{-1}>*<$psiAd1>*<$psiAu1->>,
+A22=<<(-1)s{-1}>*<$psiAu2>*<$psiAd2->>+<<s{-1}>*<$psiAd2>*<$psiAu2->>,
+
+A31=<<$psiAu1>*<$psiAu1->>+<<-1>*<$psiAd1>*<$psiAd1->>,
+A32=<<$psiAu2>*<$psiAu2->>+<<-1>*<$psiAd2>*<$psiAd2->>,
+
+AB1=<<$A11>*<[f400]>>+<<$A21>*<[f410]>>+<<$A31>*<[f420]>>,
+AA1=<<$A11>*<$A11>>+<<$A21>*<$A21>>+<<$A31>*<$A31>>,
+
+AB2=<<$A12>*<[f400]>>+<<$A22>*<[f410]>>+<<$A32>*<[f420]>>,
+AA2=<<$A12>*<$A12>>+<<$A22>*<$A22>>+<<$A32>*<$A32>>,
+
+
+A1=[f100][f-110] + [f110][f-100],
+A2=((-1)s{-1})[f100][f-110] + (s{-1})[f110][f-100],
+A3=[f100][f-100] + (-1)[f110][f-110],
+
+AB=<<$A1>*<[f400]>>+<<$A2>*<[f410]>>+<<$A3>*<[f420]>>,
+AA=<<$A1>*<$A1>>+<<$A2>*<$A2>>+<<$A3>*<$A3>>
+
+&
+#!identities
+[f400][f400]=(1),
+[f400][f410]=(s{-1})[f420],
+[f420][f410]=((-1)s{-1})[f400],
+[f410][f410]=(1),
+[f410][f420]=(s{-1})[f400],
+[f400][f420]=((-1)s{-1})[f410],
+[f420][f420]=(1),
+[f420][f400]=(s{-1})[f410],
+[f410][f400]=((-1)s{-1})[f420]
+
+&
+#!propagator
+101;102: 1 ,111;112: 1 ,
+102;101: -1 ,112;111: -1
+&
+#!input_polynomial
+<
+  <1>
+  + <<(1/2)[l0]>*<$AB1>>
+  + <<(1/2)[l1]>*<$AA1>>
+>*<
+  <1>
+  + <<(1/2)[l0]>*<$AB2>>
+  + <<(1/2)[l1]>*<$AA2>>
+>
+
+&
+
+#!id_table
+0: <$AB>,
+1: <2>*<$AA>
+
+&
+
+#!labels
+0:"at" ,
+1:"aa" ,
diff --git a/figs/hierarchical_sd.fig/sd_vector_field.gnuplot b/figs/hierarchical_sd.fig/sd_vector_field.gnuplot
new file mode 100644
index 0000000..9952d15
--- /dev/null
+++ b/figs/hierarchical_sd.fig/sd_vector_field.gnuplot
@@ -0,0 +1,25 @@
+set ylabel "$\\ell_1$" norotate
+set xlabel "$\\ell_0$"
+
+# default output canvas size: 12.5cm x 8.75cm
+set term lua tikz size 8,6 standalone
+
+unset key
+
+set pointsize 1
+
+# color functions
+
+rgb(r,g,b) = 65536 * int(r) + 256 * int(g) + int(b)
+color(x) = rgb(x*255, 0 , (1-x)*255)
+
+# get min/max
+set yrange [-1000:1000]
+stats "sd_vector_field.dat" u (log10(sqrt(($3-$1)**2+($4-$2)**2)))
+
+set xrange [-1:0.5]
+set yrange [-0.4:0.4]
+
+rescale=40
+
+plot "sd_vector_field.dat" u 1:2:(($3-$1)/sqrt(($3-$1)**2+($4-$2)**2)/rescale):(($4-$2)/sqrt(($3-$1)**2+($4-$2)**2)/rescale):(color((log10(sqrt(($3-$1)**2+($4-$2)**2))-STATS_min)/(STATS_max-STATS_min))) with vec filled head linecolor rgb variable
diff --git a/figs/lib/graphene.sty b/figs/lib/graphene.sty
new file mode 100644
index 0000000..4f044ae
--- /dev/null
+++ b/figs/lib/graphene.sty
@@ -0,0 +1,48 @@
+% sqrt(3)/2
+\def\sqrtThOT{0.866}
+% sqrt(3)
+\def\sqrtTh{1.732}
+% pi
+\def\pival{3.14159}
+
+% base vectorsj
+\def\Lo{(1.5,\sqrtThOT)}
+\def\Lt{(1.5,-\sqrtThOT)}
+\def\mLo{(-1.5,-\sqrtThOT)}
+\def\mLt{(-1.5,\sqrtThOT)}
+
+% filled square at #1 of side-length #2
+\def\fullsquare#1#2{\fill#1++(-#2,-#2)--++(#2,0)--++(#2,0)--++(0,#2)--++(0,#2)--++(-#2,0)--++(-#2,0)--++(0,-#2)--++(0,-#2);}
+% in color
+\def\fullsquarecolor#1#2#3{\fill[color=#3]#1++(-#2,-#2)--++(#2,0)--++(#2,0)--++(0,#2)--++(0,#2)--++(-#2,0)--++(-#2,0)--++(0,-#2)--++(0,-#2);}
+
+% type-a site at #1
+\def\asite#1{\fill#1circle(.1);}
+% in color
+\def\asitec#1#2{\fill[color=#2]#1circle(.1);}
+% type-b site at #1
+\def\bsite#1{\fullsquare{#1}{.1}}
+% in color
+\def\bsitec#1#2{\fullsquarecolor{#1}{.1}{#2}}
+
+% base cell with a site at #1
+\def\cell#1{
+  \draw#1--++(0:1);
+  \draw#1--++(120:1);
+  \draw#1--++(240:1);
+  \asite{#1}
+  \draw#1++(1,0)--++(0:-1);
+  \draw#1++(1,0)--++(120:-1);
+  \draw#1++(1,0)--++(240:-1);
+  \bsite{#1++(1,0)}
+}
+
+% grpahene grid at #1 of width #2 and height #3
+\def\graphene#1#2#3{
+  \foreach \i in {0,...,#2}{
+    \foreach \j in {0,...,#3}{
+      \cell{#1++(\i*3,-2*\j*\sqrtThOT)}
+      \cell{#1++(\i*3+1.5,\sqrtThOT-2*\j*\sqrtThOT)}
+    }
+  }
+}
diff --git a/libs/ian.cls b/libs/ian.cls
new file mode 100644
index 0000000..e8d3add
--- /dev/null
+++ b/libs/ian.cls
@@ -0,0 +1,667 @@
+%%
+%% Ian's class file
+%%
+
+%% TeX format
+\NeedsTeXFormat{LaTeX2e}[1995/12/01]
+
+%% class name
+\ProvidesClass{ian}[2017/09/29]
+
+%% boolean to signal that this class is being used
+\newif\ifianclass
+\ianclasstrue
+
+%% options
+% no section numbering in equations
+\DeclareOption{section_in_eq}{\sectionsineqtrue}
+\DeclareOption{section_in_fig}{\sectionsinfigtrue}
+\DeclareOption{section_in_all}{\sectionsineqtrue\sectionsinfigtrue}
+\DeclareOption{subsection_in_eq}{\subsectionsineqtrue}
+\DeclareOption{subsection_in_fig}{\subsectionsinfigtrue}
+\DeclareOption{subsection_in_all}{\subsectionsineqtrue\subsectionsinfigtrue}
+\DeclareOption{no_section_in_eq}{\sectionsineqfalse}
+\DeclareOption{no_section_in_fig}{\sectionsinfigfalse}
+\DeclareOption{no_section_in_all}{\sectionsineqfalse\sectionsinfigfalse}
+\DeclareOption{no_subsection_in_eq}{\subsectionsineqfalse}
+\DeclareOption{no_subsection_in_fig}{\subsectionsinfigfalse}
+\DeclareOption{no_subsection_in_all}{\subsectionsineqfalse\subsectionsinfigfalse}
+
+\def\ian@defaultoptions{
+  \ExecuteOptions{section_in_all, no_subsection_in_all}
+  \ProcessOptions
+
+  %% required packages
+  \RequirePackage{color}
+  \RequirePackage{marginnote}
+  \RequirePackage{amssymb}
+  \PassOptionsToPackage{hidelinks}{hyperref}
+  \RequirePackage{hyperref}
+
+  \pagestyle{plain}
+}
+
+%% paper dimensions
+\setlength\paperheight{297mm}
+\setlength\paperwidth{210mm}
+
+%% fonts
+\input{size11.clo}
+\DeclareOldFontCommand{\rm}{\normalfont\rmfamily}{\mathrm}
+\DeclareOldFontCommand{\sf}{\normalfont\sffamily}{\mathsf}
+\DeclareOldFontCommand{\tt}{\normalfont\ttfamily}{\mathtt}
+\DeclareOldFontCommand{\bf}{\normalfont\bfseries}{\mathbf}
+\DeclareOldFontCommand{\it}{\normalfont\itshape}{\mathit}
+
+%% text dimensions
+\hoffset=-50pt
+\voffset=-72pt
+\textwidth=460pt
+\textheight=704pt
+
+
+%% remove default indentation
+\parindent=0pt
+%% indent command
+\def\indent{\hskip20pt}
+
+%% something is wrong with \thepage, redefine it
+\gdef\thepage{\the\c@page}
+
+%% array lines (to use the array environment)
+\setlength\arraycolsep{5\p@}
+\setlength\arrayrulewidth{.4\p@}
+
+
+%% correct vertical alignment at the end of a document
+\AtEndDocument{
+  \vfill
+  \eject
+}
+
+
+%% hyperlinks
+% hyperlinkcounter
+\newcounter{lncount}
+% hyperref anchor
+\def\hrefanchor{%
+\stepcounter{lncount}%
+\hypertarget{ln.\thelncount}{}%
+}
+
+%% define a command and write it to aux file
+\def\outdef#1#2{%
+  % define command%
+  \expandafter\xdef\csname #1\endcsname{#2}%
+  % hyperlink number%
+  \expandafter\xdef\csname #1@hl\endcsname{\thelncount}%
+  % write command to aux%
+  \immediate\write\@auxout{\noexpand\expandafter\noexpand\gdef\noexpand\csname #1\endcsname{\csname #1\endcsname}}%
+  \immediate\write\@auxout{\noexpand\expandafter\noexpand\gdef\noexpand\csname #1@hl\endcsname{\thelncount}}%
+}
+
+%% can call commands even when they are not defined
+\def\safe#1{%
+  \ifdefined#1%
+    #1%
+  \else%
+    {\color{red}\bf?}%
+  \fi%
+}
+
+%% define a label for the latest tag
+%% label defines a command containing the string stored in \tag
+\def\deflabel{
+  \def\label##1{\expandafter\outdef{label@##1}{\safe\tag}}
+
+  \def\ref##1{%
+    % check whether the label is defined (hyperlink runs into errors if this check is omitted)
+    \ifcsname label@##1@hl\endcsname%
+      \hyperlink{ln.\csname label@##1@hl\endcsname}{{\color{blue}\safe\csname label@##1\endcsname}}%
+    \else%
+      \ifcsname label@##1\endcsname%
+	{\color{blue}\csname ##1\endcsname}%
+	\else%
+	{\bf ??}%
+      \fi%
+    \fi%
+  }
+}
+
+
+%% make a custom link at any given location in the document
+\def\makelink#1#2{%
+  \hrefanchor%
+  \outdef{label@#1}{#2}%
+}
+
+
+%% section command
+% counter
+\newcounter{sectioncount}
+% space before section
+\newlength\secskip
+\setlength\secskip{40pt}
+% a prefix to put before the section number, e.g. A for appendices
+\def\sectionprefix{}
+% define some lengths
+\newlength\secnumwidth
+\newlength\sectitlewidth
+\def\section#1{
+  % reset counters
+  \stepcounter{sectioncount}
+  \setcounter{subsectioncount}{0}
+  \ifsectionsineq
+    \setcounter{seqcount}0
+  \fi
+  \ifsectionsinfig
+    \setcounter{figcount}0
+  \fi
+
+  % space before section (if not first)
+  \ifnum\thesectioncount>1
+    \vskip\secskip
+    \penalty-1000
+  \fi
+
+  % hyperref anchor
+  \hrefanchor
+  % define tag (for \label)
+  \xdef\tag{\sectionprefix\thesectioncount}
+
+  % get widths
+  \def\@secnum{{\bf\Large\sectionprefix\thesectioncount.\hskip10pt}}
+  \settowidth\secnumwidth{\@secnum}
+  \setlength\sectitlewidth\textwidth
+  \addtolength\sectitlewidth{-\secnumwidth}
+
+  % print name
+  \parbox{\textwidth}{
+  \@secnum
+  \parbox[t]{\sectitlewidth}{\Large\bf #1}}
+
+  % write to table of contents
+  \iftoc
+    % save lncount in aux variable which is written to toc
+    \immediate\write\tocoutput{\noexpand\expandafter\noexpand\edef\noexpand\csname toc@sec.\thesectioncount\endcsname{\thelncount}}
+    \write\tocoutput{\noexpand\tocsection{#1}{\thepage}}
+  \fi
+
+  %space
+  \par\penalty10000
+  \bigskip\penalty10000
+}
+
+%% subsection
+% counter
+\newcounter{subsectioncount}
+% space before subsection
+\newlength\subsecskip
+\setlength\subsecskip{30pt}
+\def\subsection#1{
+  % counters
+  \stepcounter{subsectioncount}
+  \setcounter{subsubsectioncount}{0}
+  \ifsubsectionsineq
+    \setcounter{seqcount}0
+  \fi
+  \ifsubsectionsinfig
+    \setcounter{figcount}0
+  \fi
+
+  % space before subsection (if not first)
+  \ifnum\thesubsectioncount>1
+    \vskip\subsecskip
+    \penalty-500
+  \fi
+
+  % hyperref anchor
+  \hrefanchor
+  % define tag (for \label)
+  \xdef\tag{\sectionprefix\thesectioncount.\thesubsectioncount}
+
+  % get widths
+  \def\@secnum{{\bf\large\hskip.5cm\sectionprefix\thesectioncount.\thesubsectioncount.\hskip5pt}}
+  \settowidth\secnumwidth{\@secnum}
+  \setlength\sectitlewidth\textwidth
+  \addtolength\sectitlewidth{-\secnumwidth}
+  % print name
+  \parbox{\textwidth}{
+  \@secnum
+  \parbox[t]{\sectitlewidth}{\large\bf #1}}
+
+  % write to table of contents
+  \iftoc
+    % save lncount in aux variable which is written to toc
+    \immediate\write\tocoutput{\noexpand\expandafter\noexpand\edef\noexpand\csname toc@subsec.\thesectioncount.\thesubsectioncount\endcsname{\thelncount}}
+    \write\tocoutput{\noexpand\tocsubsection{#1}{\thepage}}
+  \fi
+
+  % space
+  \par\penalty10000
+  \medskip\penalty10000
+}
+
+%% subsubsection
+% counter
+\newcounter{subsubsectioncount}
+% space before subsubsection
+\newlength\subsubsecskip
+\setlength\subsubsecskip{20pt}
+\def\subsubsection#1{
+  % counters
+  \stepcounter{subsubsectioncount}
+
+  % space before subsubsection (if not first)
+  \ifnum\thesubsubsectioncount>1
+    \vskip\subsubsecskip
+    \penalty-500
+  \fi
+
+  % hyperref anchor
+  \hrefanchor
+  % define tag (for \label)
+  \xdef\tag{\sectionprefix\thesectioncount.\thesubsectioncount.\thesubsubsectioncount}
+
+  % get widths
+  \def\@secnum{{\bf\hskip1.cm\sectionprefix\thesectioncount.\thesubsectioncount.\thesubsubsectioncount.\hskip5pt}}
+  \settowidth\secnumwidth{\@secnum}
+  \setlength\sectitlewidth\textwidth
+  \addtolength\sectitlewidth{-\secnumwidth}
+  % print name
+  \parbox{\textwidth}{
+  \@secnum
+  \parbox[t]{\sectitlewidth}{\large\bf #1}}
+
+  % write to table of contents
+  \iftoc
+    % save lncount in aux variable which is written to toc
+    \immediate\write\tocoutput{\noexpand\expandafter\noexpand\edef\noexpand\csname toc@subsubsec.\thesectioncount.\thesubsectioncount.\thesubsubsectioncount\endcsname{\thelncount}}
+    \write\tocoutput{\noexpand\tocsubsubsection{#1}{\thepage}}
+  \fi
+
+  % space
+  \par\penalty10000
+  \medskip\penalty10000
+}
+
+%% itemize
+\newlength\itemizeskip
+% left margin for items
+\setlength\itemizeskip{20pt}
+\newlength\itemizeseparator
+% space between the item symbol and the text
+\setlength\itemizeseparator{5pt}
+% penalty preceding an itemize
+\newcount\itemizepenalty
+\itemizepenalty=0
+% counter counting the itemize level
+\newcounter{itemizecount}
+
+% item symbol
+\def\itemizept#1{
+  \ifnum#1=1
+    \textbullet
+  \else
+    $\scriptstyle\blacktriangleright$
+  \fi
+}
+
+
+\newlength\current@itemizeskip
+\setlength\current@itemizeskip{0pt}
+\def\itemize{%
+  \par\expandafter\penalty\the\itemizepenalty\medskip\expandafter\penalty\the\itemizepenalty%
+  \addtocounter{itemizecount}{1}%
+  \addtolength\current@itemizeskip{\itemizeskip}%
+  \leftskip\current@itemizeskip%
+}
+\def\enditemize{%
+  \addtocounter{itemizecount}{-1}%
+  \addtolength\current@itemizeskip{-\itemizeskip}%
+  \par\expandafter\penalty\the\itemizepenalty\leftskip\current@itemizeskip%
+  \medskip\expandafter\penalty\the\itemizepenalty%
+}
+
+% item, with optional argument to specify the item point
+% @itemarg is set to true when there is an optional argument
+\newif\if@itemarg
+\def\item{%
+  % check whether there is an optional argument (if there is none, add on empty '[]')
+  \@ifnextchar [{\@itemargtrue\@itemx}{\@itemargfalse\@itemx[]}%
+}
+\newlength\itempt@total
+\def\@itemx[#1]{
+  \if@itemarg
+    \settowidth\itempt@total{#1}
+  \else
+    \settowidth\itempt@total{\itemizept\theitemizecount}
+  \fi
+  \addtolength\itempt@total{\itemizeseparator}
+  \par
+  \medskip
+  \if@itemarg
+    \hskip-\itempt@total#1\hskip\itemizeseparator
+  \else
+    \hskip-\itempt@total\itemizept\theitemizecount\hskip\itemizeseparator
+  \fi
+}
+
+%% prevent page breaks after itemize
+\newcount\previtemizepenalty
+\def\nopagebreakafteritemize{
+  \previtemizepenalty=\itemizepenalty
+  \itemizepenalty=10000
+}
+%% back to previous value
+\def\restorepagebreakafteritemize{
+  \itemizepenalty=\previtemizepenalty
+}
+
+%% enumerate
+\newcounter{enumerate@count}
+\def\enumerate{
+  \setcounter{enumerate@count}0
+  \let\olditem\item
+  \let\olditemizept\itemizept
+  \def\item{
+    % counter
+    \stepcounter{enumerate@count}
+    % set header
+    \def\itemizept{\theenumerate@count.}
+    % hyperref anchor
+    \hrefanchor
+    % define tag (for \label)
+    \xdef\tag{\theenumerate@count}
+    \olditem
+  }
+  \itemize
+}
+\def\endenumerate{
+  \enditemize
+  \let\item\olditem
+  \let\itemizept\olditemizept
+}
+
+
+%% equation numbering
+% counter
+\newcounter{seqcount}
+% define possible prefix to equation
+\def\eqprefix{}
+% booleans (write section or subsection in equation number)
+\newif\ifsectionsineq
+\newif\ifsubsectionsineq
+\def\seqcount{
+  \stepcounter{seqcount}
+  % the output
+  \edef\seqformat{\eqprefix\theseqcount}
+  % add subsection number
+  \ifsubsectionsineq
+    \let\tmp\seqformat
+    \edef\seqformat{\thesubsectioncount.\tmp}
+  \fi
+  % add section number
+  \ifsectionsineq
+    \let\tmp\seqformat
+    \edef\seqformat{\sectionprefix\thesectioncount.\tmp}
+  \fi
+  % define tag (for \label)
+  \xdef\tag{\seqformat}
+  % write number
+  \marginnote{\hfill(\seqformat)}
+}
+%% equation environment compatibility
+\def\equation{\hrefanchor$$\seqcount}
+\def\endequation{$$\@ignoretrue}
+
+
+%% figures
+% counter
+\newcounter{figcount}
+% booleans (write section or subsection in equation number)
+\newif\ifsectionsinfig
+\newif\ifsubsectionsinfig
+% width of figures
+\newlength\figwidth
+\setlength\figwidth\textwidth
+\addtolength\figwidth{-2.5cm}
+% caption
+\def\defcaption{
+  \long\def\caption##1{
+    \stepcounter{figcount}
+
+    % hyperref anchor
+    \hrefanchor
+
+    % the number of the figure
+    \edef\figformat{\thefigcount}
+    % add subsection number
+    \ifsubsectionsinfig
+      \let\tmp\figformat
+      \edef\figformat{\thesubsectioncount.\tmp}
+    \fi
+    % add section number
+    \ifsectionsinfig
+      \let\tmp\figformat
+      \edef\figformat{\sectionprefix\thesectioncount.\tmp}
+    \fi
+
+    % define tag (for \label)
+    \xdef\tag{\figformat}
+
+    % write
+    \hfil fig \figformat: \parbox[t]{\figwidth}{\leavevmode\small##1}
+
+    % space
+    \par\bigskip
+  }
+}
+%% short caption: centered
+\def\captionshort#1{
+  \stepcounter{figcount}
+
+  % hyperref anchor
+  \hrefanchor
+
+  % the number of the figure
+  \edef\figformat{\thefigcount}
+  % add section number
+  \ifsectionsinfig
+  \let\tmp\figformat
+  \edef\figformat{\sectionprefix\thesectioncount.\tmp}
+  \fi
+
+  % define tag (for \label)
+  \xdef\tag{\figformat}
+
+  % write
+  \hfil fig \figformat: {\small#1}
+
+  %space
+  \par\bigskip
+}
+
+%% environment
+\def\figure{
+  \par
+  \vfil\penalty100\vfilneg
+  \bigskip
+}
+\def\endfigure{
+  \par
+  \bigskip
+}
+
+
+%% start appendices
+\def\appendix{
+  \vfill
+  \pagebreak
+
+  % counter
+  \setcounter{sectioncount}0
+
+  % prefix
+  \def\sectionprefix{A}
+
+  % write
+  {\bf \LARGE Appendices}\par\penalty10000\bigskip\penalty10000
+
+  % add a mention in the table of contents
+  \iftoc
+    \immediate\write\tocoutput{\noexpand\tocappendices}\penalty10000
+  \fi
+
+  %% uncomment for new page for each appendix
+  %\def\seqskip{\vfill\pagebreak}
+}
+
+
+%% bibliography
+% size of header
+\newlength\bibheader
+\def\thebibliography#1{
+  \hrefanchor
+
+  % add a mention in the table of contents
+  \iftoc
+    % save lncount in aux variable which is written to toc
+    \immediate\write\tocoutput{\noexpand\expandafter\noexpand\edef\noexpand\csname toc@references\endcsname{\thelncount}}
+    \write\tocoutput{\noexpand\tocreferences{\thepage}}\penalty10000
+  \fi
+
+  % write
+  {\bf \LARGE References}\par\penalty10000\bigskip\penalty10000
+  % width of header
+  \settowidth\bibheader{[#1]}
+  \leftskip\bibheader
+}
+% end environment
+\def\endthebibliography{
+  \par\leftskip0pt
+} 
+
+%% bibitem command
+\def\bibitem[#1]#2{%
+  \hrefanchor%
+  \outdef{label@cite#2}{#1}%
+  \hskip-\bibheader%
+  \makebox[\bibheader]{\cite{#2}\hfill}%
+}
+
+%% cite command
+% @tempswa is set to true when there is an optional argument
+\newif\@tempswa
+\def\cite{%
+  % check whether there is an optional argument (if there is none, add on empty '[]')
+  \@ifnextchar [{\@tempswatrue\@citex}{\@tempswafalse\@citex[]}%
+}
+% command with optional argument
+\def\@citex[#1]#2{\leavevmode%
+  % initialize loop
+  \let\@cite@separator\@empty%
+  % format
+  \@cite{%
+    % loop over ',' separated list
+    \@for\@cite@:=#2\do{%
+      % text to add at each iteration of the loop (separator between citations)
+      \@cite@separator\def\@cite@separator{,\ }%
+      % add entry to citelist
+      \@writecitation{\@cite@}%
+      \ref{cite\@cite@}%
+    }%
+  }%
+  % add optional argument text (as an argument to '\@cite')
+  {#1}%
+}
+\def\@cite#1#2{%
+  [#1\if@tempswa , #2\fi]%
+}
+%% add entry to citelist after checking it has not already been added
+\def\@writecitation#1{%
+  \ifcsname if#1cited\endcsname%
+  \else%
+    \expandafter\newif\csname if#1cited\endcsname%
+    \immediate\write\@auxout{\string\citation{#1}}%
+  \fi%
+}
+
+%% table of contents
+% boolean
+\newif\iftoc
+\def\tableofcontents{
+  \hfil{\bf Table of contents}\par\penalty10000\smallskip\penalty10000
+
+  % copy content from file
+  \IfFileExists{\jobname.toc}{\input{\jobname.toc}}{{\tt error: table of contents missing}}
+
+  % open new toc
+  \newwrite\tocoutput
+  \immediate\openout\tocoutput=\jobname.toc
+
+  \toctrue
+}
+%% close file
+\AtEndDocument{
+  % close toc
+  \iftoc
+    \immediate\closeout\tocoutput
+  \fi
+}
+
+
+%% fill line with dots
+\def\leaderfill{\leaders\hbox to 1em {\hss. \hss}\hfill}
+
+%% same as sectionprefix
+\def\tocsectionprefix{}
+
+%% toc formats
+\newcounter{tocsectioncount}
+\def\tocsection #1#2{
+  \stepcounter{tocsectioncount}
+  \setcounter{tocsubsectioncount}{0}
+  \setcounter{tocsubsubsectioncount}{0}
+  % write
+  \smallskip\hyperlink{ln.\csname toc@sec.\thetocsectioncount\endcsname}{{\bf \tocsectionprefix\thetocsectioncount}.\hskip5pt {\color{blue}#1}\leaderfill#2}\par
+}
+\newcounter{tocsubsectioncount}
+\def\tocsubsection #1#2{
+  \stepcounter{tocsubsectioncount}
+  \setcounter{tocsubsubsectioncount}{0}
+  % write
+  {\hskip10pt\hyperlink{ln.\csname toc@subsec.\thetocsectioncount.\thetocsubsectioncount\endcsname}{{\bf \thetocsubsectioncount}.\hskip5pt {\color{blue}\small #1}\leaderfill#2}}\par
+}
+\newcounter{tocsubsubsectioncount}
+\def\tocsubsubsection #1#2{
+  \stepcounter{tocsubsubsectioncount}
+  % write
+  {\hskip20pt\hyperlink{ln.\csname toc@subsubsec.\thetocsectioncount.\thetocsubsectioncount.\thetocsubsubsectioncount\endcsname}{{\bf \thetocsubsubsectioncount}.\hskip5pt {\color{blue}\small #1}\leaderfill#2}}\par
+}
+\def\tocappendices{
+  \medskip
+  \setcounter{tocsectioncount}0
+  {\bf Appendices}\par
+  \smallskip
+  \def\tocsectionprefix{A}
+}
+\def\tocreferences#1{
+  \medskip
+  {\hyperlink{ln.\csname toc@references\endcsname}{{\color{blue}\bf References}\leaderfill#1}}\par
+  \smallskip
+}
+
+
+%% definitions that must be loaded at begin document
+\let\ian@olddocument\document
+\def\document{
+  \ian@olddocument
+
+  \deflabel
+  \defcaption
+}
+
+%% end
+\ian@defaultoptions
+\endinput
diff --git a/libs/iantheo.sty b/libs/iantheo.sty
new file mode 100644
index 0000000..d33a93d
--- /dev/null
+++ b/libs/iantheo.sty
@@ -0,0 +1,162 @@
+%%
+%% iantheorem package:
+%%   Ian's customized theorem command
+%%
+
+%% boolean to signal that this package was loaded
+\newif\ifiantheo
+
+%% TeX format
+\NeedsTeXFormat{LaTeX2e}[1995/12/01]
+
+%% package name
+\ProvidesPackage{iantheo}[2016/11/10]
+
+%% options
+\newif\ifsectionintheo
+\DeclareOption{section_in_theo}{\sectionintheotrue}
+\DeclareOption{no_section_in_theo}{\sectionintheofalse}
+\newif\ifsubsectionintheo
+\DeclareOption{subsection_in_theo}{\subsectionintheotrue}
+\DeclareOption{no_subsection_in_theo}{\subsectionintheofalse}
+
+\def\iantheo@defaultoptions{
+  \ExecuteOptions{section_in_theo, no_subsection_in_theo}
+  \ProcessOptions
+
+  %%% reset at every new section
+  \ifsectionintheo
+    \let\iantheo@oldsection\section
+    \gdef\section{\setcounter{theocount}{0}\iantheo@oldsection}
+  \fi
+
+  %% reset at every new subsection
+  \ifsubsectionintheo
+    \let\iantheo@oldsubsection\subsection
+    \gdef\subsection{\setcounter{theocount}{0}\iantheo@oldsubsection}
+  \fi
+}
+
+
+%% delimiters
+\def\delimtitle#1{
+  \par%
+  \leavevmode%
+  \raise.3em\hbox to\hsize{%
+    \lower0.3em\hbox{\vrule height0.3em}%
+    \hrulefill%
+    \ \lower.3em\hbox{#1}\ %
+    \hrulefill%
+    \lower0.3em\hbox{\vrule height0.3em}%
+  }%
+  \par\penalty10000%
+}
+
+%% callable by ref
+\def\delimtitleref#1{
+  \par%
+%
+  \ifdefined\ianclass%
+    % hyperref anchor%
+    \hrefanchor%
+    % define tag (for \label)%
+    \xdef\tag{#1}%
+  \fi%
+%
+  \leavevmode%
+  \raise.3em\hbox to\hsize{%
+    \lower0.3em\hbox{\vrule height0.3em}%
+    \hrulefill%
+    \ \lower.3em\hbox{\bf #1}\ %
+    \hrulefill%
+    \lower0.3em\hbox{\vrule height0.3em}%
+  }%
+  \par\penalty10000%
+}
+
+%% no title
+\def\delim{
+  \par%
+  \leavevmode\raise.3em\hbox to\hsize{%
+    \lower0.3em\hbox{\vrule height0.3em}%
+    \hrulefill%
+    \lower0.3em\hbox{\vrule height0.3em}%
+  }%
+  \par\penalty10000%
+}
+
+%% end delim
+\def\enddelim{
+  \par\penalty10000%
+  \leavevmode%
+  \raise.3em\hbox to\hsize{%
+    \vrule height0.3em\hrulefill\vrule height0.3em%
+  }%
+  \par%
+}
+
+
+%% theorem
+% counter
+\newcounter{theocount}
+% booleans (write section or subsection in equation number)
+\def\theo#1{
+  \stepcounter{theocount}
+  \ifdefined\ianclass
+    % hyperref anchor
+    \hrefanchor
+  \fi
+  % the number
+  \def\formattheo{\thetheocount}
+  % add subsection number
+  \ifsubsectionintheo
+    \let\tmp\formattheo
+    \edef\formattheo{\thesubsectioncount.\tmp}
+  \fi
+  % add section number
+  \ifsectionintheo
+    \let\tmp\formattheo
+    \edef\formattheo{\sectionprefix\thesectioncount.\tmp}
+  \fi
+  % define tag (for \label)
+  \xdef\tag{\formattheo}
+  % write
+  \delimtitle{\bf #1 \formattheo}
+}
+\let\endtheo\enddelim
+%% theorem headers with name
+\def\theoname#1#2{
+  \theo{#1}\hfil({\it #2})\par\penalty10000\medskip%
+}
+
+
+%% qed symbol
+\def\qedsymbol{$\square$}
+\def\qed{\penalty10000\hfill\penalty10000\qedsymbol}
+
+
+%% compatibility with article class
+\ifdefined\ianclasstrue
+  \relax
+\else
+  \def\thesectioncount{\thesection}
+  \def\thesubsectioncount{\thesubsection}
+  \def\sectionprefix{}
+\fi
+
+
+%% prevent page breaks after displayed equations
+\newcount\prevpostdisplaypenalty
+\def\nopagebreakaftereq{
+  \prevpostdisplaypenalty=\postdisplaypenalty
+  \postdisplaypenalty=10000
+}
+%% back to previous value
+\def\restorepagebreakaftereq{
+  \postdisplaypenalty=\prevpostdisplaypenalty
+}
+
+
+%% end
+\iantheo@defaultoptions
+\endinput
diff --git a/libs/largearray.sty b/libs/largearray.sty
new file mode 100644
index 0000000..ad5753b
--- /dev/null
+++ b/libs/largearray.sty
@@ -0,0 +1,19 @@
+%%
+%% largearray package:
+%%   Array spanning the entire line
+%%
+
+%% TeX format
+\NeedsTeXFormat{LaTeX2e}[1995/12/01]
+
+%% package name
+\ProvidesPackage{largearray}[2016/11/10]
+
+\RequirePackage{array}
+
+%% array spanning the entire line
+\newlength\largearray@width
+\setlength\largearray@width\textwidth
+\addtolength\largearray@width{-10pt}
+\def\largearray{\begin{array}{@{}>{\displaystyle}l@{}}\hphantom{\hspace{\largearray@width}}\\[-.5cm]}
+\def\endlargearray{\end{array}}
diff --git a/libs/point.sty b/libs/point.sty
new file mode 100644
index 0000000..4a761b7
--- /dev/null
+++ b/libs/point.sty
@@ -0,0 +1,114 @@
+%%
+%% Points package:
+%%   \point commands
+%%
+
+%% TeX format
+\NeedsTeXFormat{LaTeX2e}[1995/12/01]
+
+%% package name
+\ProvidesPackage{point}[2017/06/13]
+
+%% options
+\newif\ifresetatsection
+\DeclareOption{reset_at_section}{\resetatsectiontrue}
+\DeclareOption{no_reset_at_section}{\resetatsectionfalse}
+\newif\ifresetatsubsection
+\DeclareOption{reset_at_subsection}{\resetatsubsectiontrue}
+\DeclareOption{no_reset_at_subsection}{\resetatsubsectionfalse}
+\newif\ifresetatsubsubsection
+\DeclareOption{reset_at_subsubsection}{\resetatsubsubsectiontrue}
+\DeclareOption{no_reset_at_subsubsection}{\resetatsubsubsectionfalse}
+\newif\ifresetattheo
+\DeclareOption{reset_at_theo}{\resetattheotrue}
+\DeclareOption{no_reset_at_theo}{\resetattheofalse}
+
+\def\point@defaultoptions{
+  \ExecuteOptions{reset_at_section, reset_at_subsection, reset_at_subsubsection, no_reset_at_theo}
+  \ProcessOptions
+
+  %% reset at every new section
+  \ifresetatsection
+    \let\point@oldsection\section
+    \gdef\section{\resetpointcounter\point@oldsection}
+  \fi
+  %% reset at every new subsection
+  \ifresetatsubsection
+    \let\point@oldsubsection\subsection
+    \gdef\subsection{\resetpointcounter\point@oldsubsection}
+  \fi
+  %% reset at every new subsubsection
+  \ifresetatsubsubsection
+    \let\point@oldsubsubsection\subsubsection
+    \gdef\subsubsection{\resetpointcounter\point@oldsubsubsection}
+  \fi
+
+  %% reset at every new theorem
+  \ifresetattheo
+    \ifdefined\iantheotrue
+      \let\point@oldtheo\theo
+      \gdef\theo{\resetpointcounter\point@oldtheo}
+    \fi
+  \fi
+}
+
+
+%% point
+% counter
+\newcounter{pointcount}
+\def\point{
+  \stepcounter{pointcount}
+  \setcounter{subpointcount}{0}
+  % hyperref anchor (only if the class is 'ian')
+  \ifdefined\ifianclass
+    \hrefanchor
+    % define tag (for \label)
+    \xdef\tag{\thepointcount}
+  \fi
+  % header
+  \indent{\bf \thepointcount\ - }
+}
+
+%% subpoint
+% counter
+\newcounter{subpointcount}
+\def\subpoint{
+  \stepcounter{subpointcount}
+  \setcounter{subsubpointcount}0
+  % hyperref anchor (only if the class is 'ian')
+  \ifdefined\ifianclass
+    \hrefanchor
+    % define tag (for \label)
+    \xdef\tag{\thepointcount-\thesubpointcount}
+  \fi
+  % header
+  \indent\hskip.5cm{\bf \thepointcount-\thesubpointcount\ - }
+}
+
+%% subsubpoint
+% counter
+\newcounter{subsubpointcount}
+\def\subsubpoint{
+  \stepcounter{subsubpointcount}
+  % hyperref anchor (only if the class is 'ian')
+  \ifdefined\ifianclass
+    \hrefanchor
+    % define tag (for \label)
+    \xdef\tag{\thepointcount-\thesubpointcount-\thesubsubpointcount}
+  \fi
+  \indent\hskip1cm{\bf \thepointcount-\thesubpointcount-\thesubsubpointcount\ - }
+}
+
+
+%% reset point counters
+\def\resetpointcounter{
+  \setcounter{pointcount}{0}
+  \setcounter{subpointcount}{0}
+  \setcounter{subsubpointcount}{0}
+}
+
+
+
+%% end
+\point@defaultoptions
+\endinput