15gij/Giuliani_Jauslin_2015.tex
2016-09-11 23:42:37 +00:00

3912 lines
326 KiB
TeX

\documentclass{kiss}
% load packages
\usepackage{header}
% bibliography commands
\usepackage{BBlog}
% miscellaneous commands
\usepackage{toolbox}
% main style file
\usepackage{iansecs}
\usepackage{local}
\begin{document}
\pagestyle{empty}
\hfil{\LARGE\bf The ground state construction of bilayer graphene}\par
\vskip20pt
{\bf\hfil Alessandro Giuliani, Ian Jauslin}\par
\vskip20pt
\hfil 2015
\hugeskip
\indent We consider a model of half-filled bilayer graphene, in which the three dominant Slonczewski-Weiss-McClure hopping parameters are retained, in the presence of short range interactions. Under a smallness assumption on the interaction strength $U$ as well as on the inter-layer hopping $\epsilon$, we construct the ground state in the thermodynamic limit, and prove that the pressure and two point Schwinger function, away from its singularities, are analytic in $U$, uniformly in $\epsilon$. The interacting Fermi surface is degenerate, and consists of eight Fermi points, two of which are protected by symmetries, while the locations of the other six are renormalized by the interaction, and the effective dispersion relation at the Fermi points is conical. The construction reveals the presence of different energy regimes, where the effective behavior of correlation functions changes qualitatively. The analysis of the crossover between regimes plays an important role in the proof of analyticity and in the uniform control of the radius of convergence. The proof is based on a rigorous implementation of fermionic renormalization group methods, including determinant estimates for the renormalized expansion.\par
\hugeskip
\tableofcontents
\vfill\eject
\pagestyle{plain}
\setcounter{page}1
\section{Introduction}
\label{introduction}
\indent Graphene, a one-atom thick layer of graphite, has captivated a large part of the scientific community for the past decade. With good reason: as was shown by A.~Geim's team, graphene is a stable two-dimensional crystal with very peculiar electronic properties~[\cite{NGe04}]. The mere fact that a two-dimensional crystal can be synthesized, and manipulated, at room temperature without working inside a vacuum~[\cite{Ge11}] is, in and of itself, quite surprising. But the most interesting features of graphene lay within its electronic properties. Indeed, electrons in graphene were found to have an extremely high mobility~[\cite{NGe04}], which could make it a good candidate to replace silicon in microelectronics; and they were later found to behave like massless Dirac Fermions~[\cite{NGe05}, \cite{ZTe05}], which is of great interest for the study of fundamental Quantum Electro-Dynamics. These are but a few of the intriguing features~[\cite{GN07}] that have prompted a lively response from the scientific community.\par
\indent These peculiar electronic properties stem from the particular energy structure of graphene. It consists of two energy bands, that meet at exactly two points, called the {\it Fermi points}~[\cite{Wa47}]. Graphene is thus classified as a {\it semi-metal}: it is not a {\it semi-conductor} because there is no gap between its energy bands, nor is it a {\it metal} either since the bands do not overlap, so that the density of charge carriers vanishes at the Fermi points. Furthermore, the bands around the Fermi points are approximately conical~[\cite{Wa47}], which explains the masslessness of the electrons in graphene, and in turn their high mobility.\par
\bigskip
\indent Graphene is also interesting for the mathematical physics community: its free energy and correlation functions, in particular its conductivity, can be computed non-perturbatively using constructive Renormalization Group (RG) techniques~[\cite{GM10}, \cite{GMP11}, \cite{GMP12}], at least if it is at {\it half-filling}, the interaction is {\it short-range} and its strength is {\it small enough}. This is made possible, again, by the special energy structure of graphene. Indeed, since the {\it propagator} (in the quantum field theory formalism) diverges at the Fermi points, the fact that there are only two such singularities in graphene instead of a whole line of them (which is what one usually finds in two-dimensional theories), greatly simplifies the RG analysis. Furthermore, the fact that the bands are approximately conical around the Fermi points, implies that a short-range interaction between electrons is {\it irrelevant} in the RG sense, which means that one need only worry about the renormalization of the propagator, which can be controlled.\par
\indent Using these facts, the formalism developed in~[\cite{BG90}] has been applied in~[\cite{GM10}, \cite{GMP12}] to express the free energy and correlation functions as convergent series.\par
\indent Let us mention that the
case of Coulomb interactions is more difficult, in that the effective interaction is marginal in an RG sense. In this case, the theory has been constructed at all orders in renormalized perturbation theory
[\cite{GMP10}, \cite{GMP11b}], but a non-perturbative construction is still lacking.
\par
\bigskip
\indent In the present work, we shall extend the results of~[\cite{GM10}] by performing an RG analysis of half-filled {\it bilayer} graphene with short-range interactions. Bilayer graphene consists of two layers of graphene in so-called {\it Bernal} or {\it AB}
stacking (see below). Before the works of A.~Geim et al.~[\cite{NGe04}], graphene was mostly studied in order to understand the properties of graphite, so it was natural to investigate the properties of multiple layers of
graphene, starting with the bilayer~[\cite{Wa47}, \cite{SW58}, \cite{Mc57}]. A common model for hopping electrons on graphene bilayers is the so-called {\it Slonczewski-Weiss-McClure} model,
which is usually studied by retaining only certain hopping terms, depending on the energy regime one is interested in: including more hopping terms corresponds to probing the system at lower energies.
The fine structure of the Fermi surface and the behavior of the
dispersion relation around it depends on which hoppings are considered or, equivalently, on the range of energies under inspection.
\indent In a first approximation, the energy structure of bilayer graphene is similar to that of the monolayer: there are only two Fermi points,
and the dispersion relation is approximately conical around them. This picture is valid for energy scales larger than the transverse hopping between the two layers, referred to in the following as the {\it first regime}.
At lower energies, the effective dispersion relation around the two Fermi points appears to be approximately {\it parabolic}, instead of conical. This implies that the effective mass of the electrons in bilayer graphene does not vanish, unlike those in the monolayer, which has been observed
experimentally~[\cite{NMe06}].\par
\indent From an RG point of view, the parabolicity implies that the electron interactions are {\it marginal} in bilayer graphene, thus making the RG analysis non-trivial. The flow of the effective couplings
has been studied by O.~Vafek~[\cite{Va10}, \cite{VY10}], who has found that it diverges logarithmically, and has identified the most divergent channels,
thus singling out which of the possible quantum instabilities are dominant (see also~[\cite{TV12}]).
However, as was mentioned earlier, the assumption of parabolic dispersion relation is only an approximation,
valid in a range of energies between the scale of the transverse hopping and a second threshold, proportional to the cube of the transverse hopping (asymptotically, as this hopping goes to zero).
This range will be called the {\it second regime}.\par
\indent By studying the smaller energies in more detail, one finds~[\cite{MF06}] that around each of the Fermi points, there are three extra Fermi points, forming a tiny equilateral triangle around the original ones. This is referred to in the literature as {\it trigonal warping}. Furthermore, around each of the now eight Fermi points, the energy bands are approximately conical. This means that, from an RG perspective, the
logarithmic divergence studied in [\cite{Va10}] is cut off at the energy scale where the conical nature of the eight Fermi points becomes observable (i.e. at the end of the second regime). At lower energies
the electron interaction is irrelevant in the RG sense, which implies that the flow of the effective interactions remains bounded at low energies. Therefore, the analysis of [\cite{Va10}] is meaningful
only if the flow of the effective constants has grown significantly in the second regime.
\par
\indent However, our analysis shows that the flow of the effective couplings in this regime does not grow at all, due to their smallness after integration over the first regime,
which we quantify in terms both of the bare couplings and of the transverse hopping. This puts into question the physical relevance of the ``instabilities'' coming from the logarithmic divergence in the second regime,
at least in the case we are treating, namely small interaction strength and small interlayer hopping.\par
\indent The transition from a normal phase to one with broken symmetry as the interaction strength is increased from small to intermediate values was studied in~[\cite{CTV12}] at second order in perturbation theory. Therein, it was found that while at small bare couplings the infrared flow is convergent, at larger couplings it tends to increase, indicating a transition towards an {\it electronic nematic state}.\par
\indent Let us also mention that the third regime is not believed to give an adequate description of the system at arbitrarily small energies: at energies smaller than a third threshold (proportional to the fourth power of the transverse hopping) one finds~[\cite{PP06}] that the six extra Fermi points around the two original ones, are actually microscopic ellipses. The analysis of the thermodynamic properties of the system in this regime (to be called the fourth regime) requires new ideas and techniques, due to the extended nature of the singularity, and goes beyond the scope of this paper. It may be possible to adapt the ideas of~[\cite{BGM06}] to this regime, which would allow us to carry out the analysis down to temperatures that are exponentially low in the coupling constant, and we hope to come back to this issue in a future publication.
\par
\bigskip
\indent To summarize, at weak coupling and small transverse hopping, we can distinguish four energy regimes: in the first, the system behaves like two uncoupled monolayers, in the second, the energy bands are approximately parabolic, in the third, the trigonal warping is taken into account and the bands are approximately conical, and in the fourth, six of the Fermi points become small curves. We shall treat the first, second and third regimes,
which corresponds to retaining only the three dominant Slonczewski-Weiss-McClure hopping parameters.
Informally, we will prove that {\it the interacting half-filled system is analytically close to the non-interacting one} in these regimes, and that the effect of the interaction is merely to renormalize the
hopping parameters. The proof depends on a sharp multiscale control of the crossover regions separating one regime from the next.
\par
\bigskip
\indent We will now give a quick description of the model, and a precise statement of the main result of the present work, followed by a brief outline of its proof.\par
\subsection{Definition of the model}
\label{intromodelsec}
\indent We shall consider a crystal of bilayer graphene, which is made of two honeycomb lattices in {\it Bernal} or {\it AB} stacking, as shown in figure~\ref{baregraph}. We can identify four inequivalent types of sites in the lattice, which we denote by $a$, $\tilde b$, $\tilde a$ and $b$ (we choose this peculiar order for practical reasons which will become apparent in the following).\par
\bigskip
\begin{figure}
\hfil\hskip-1.5cm\includegraphics{Figs/basegrid.pdf}\par
\caption{\ta\ and \tb\ represent atoms of type $a$ and $b$ on the lower layer and \tat\ and \tbt\ represent atoms of type $\tilde a$ and $\tilde b$ on the upper layer. Full lines join nearest neighbors within the lower layer and dashed lines join nearest neighbors within the upper layer.}
\label{baregraph}\end{figure}
\indent We consider a Hamiltonian of the form
\begin{equation}\mathcal H=\mathcal H_0+\mathcal H_I\label{fullham}\end{equation}
where the {\it free Hamiltonian} $\mathcal H_0$ plays the role of a kinetic energy for the electrons, and the {\it interaction Hamiltonian} $\mathcal H_I$ describes the interaction between electrons.\par
\bigskip
\indent $\mathcal H_0$ is given by a {\it tight-binding} approximation, which models the movement of electrons in terms of {\it hoppings} from one atom to the next. There are four inequivalent types of hoppings which we shall consider here, each of which will be associated a different {\it hopping strength} $\gamma_i$. Namely, the hoppings between neighbors of type $a$ and $b$, as well as $\tilde a$ and $\tilde b$ will be associated a hopping strength $\gamma_0$; $a$ and $\tilde b$ a strength $\gamma_1$; $\tilde a$ and $b$ a strength $\gamma_3$; $\tilde a$ and $a$, and $\tilde b$ and $b$ a strength $\gamma_4$ (see figure~\ref{intergraph}). We can thus express $H_0$ in {\it second quantized} form in {\it momentum space} at {\it zero chemical potential} as~[\cite{Wa47}, \cite{SW58}, \cite{Mc57}]
\begin{equation}
\mathcal H_0=\frac{1}{|\hat\Lambda|}\sum_{ k\in\hat \Lambda}\hat A_{ k}^\dagger H_0( k)\hat A_{ k}\label{hamkintro}\end{equation}
\begin{equation}\hat A_{ k}:=\left(\begin{array}{c}\hat a_{ k}\\\hat{\tilde b}_{ k}\\\hat{\tilde a}_{ k}\\\hat b_{ k}\end
{array}\right)\mathrm{\ and\ }H_0( k):=
-\left(\begin{array}{*{4}{c}}
\Delta&\gamma_1&\gamma_4\Omega(k)&\gamma_0\Omega^*(k)\\[0.2cm]
\gamma_1&\Delta&\gamma_0\Omega(k)&\gamma_4\Omega^*(k)\\[0.2cm]
\gamma_4\Omega^*(k)&\gamma_0\Omega^*(k)&0&\gamma_3\Omega(k)e^{3ik_x}\\[0.2cm]
\gamma_0\Omega(k)&\gamma_4\Omega(k)&\gamma_3\Omega^*(k)e^{-3ik_x}&0
\end{array}\right)\label{hmatintro}\end{equation}
in which $\hat a_k$, $\hat{\tilde b}_k$, $\hat{\tilde a}_k$ and $\hat b_k$ are {\it annihilation operators} associated to atoms of type $a$, $\tilde b$, $\tilde a$ and $b$, $k\equiv(k_x,k_y)$, $\hat\Lambda$ is the {\it first Brillouin zone}, and $\Omega(k):=1+2e^{-i\frac32k_x}\cos\left(\frac{\sqrt3}2k_y\right)$.
These objects will be properly defined in section~\ref{modelsec}. The $\Delta$ parameter in $H_0$ models a shift in the chemical potential around atoms of type $a$ and $\tilde b$~[\cite{SW58}, \cite{Mc57}]. We choose the energy unit in such a way that $\gamma_0=1$. The hopping strengths have been measured experimentally in graphite~[\cite{DD02}, \cite{TDD77}, \cite{MMD79}, \cite{DDe79}] and in bilayer graphene samples~[\cite{ZLe08}, \cite{MNe07}]; their values are given in the following table:
\begin{equation}
\begin{array}{|c|c|c|}
\cline{2-3}
\multicolumn1{c|}{}&\mathrm{bilayer\ graphene~[\cite{MNe07}]}&\mathrm{graphite~[\cite{DD02}]}\\
\hline
\gamma_1&0.10&0.12\\
\gamma_3&0.034&0.10\\
\gamma_4&0.041&0.014\\
\Delta&0.006\ \mathrm{[\cite{ZLe08}]}&-0.003\\
\hline
\end{array}\label{tabgamma}\end{equation}
We notice that the relative order of magnitude of $\gamma_3$ and $\gamma_4$ is quite different in graphite and in bilayer graphene. In the latter, $\gamma_1$ is somewhat small, and $\gamma_3$ and $\gamma_4$ are of the same order, whereas $\Delta$ is of the order of $\gamma_1^2$. We will take advantage of the smallness of the hopping strengths and treat $\gamma_1=:\epsilon$ as a small parameter: we fix
\begin{equation} \frac{\gamma_1}{\epsilon}=1,\ \frac{\gamma_3}{\epsilon}=0.33,\ \frac{\gamma_4}{\epsilon}=0.40,\ \frac{\Delta}{\epsilon^2}=0.58\label{relge}\end{equation}
and assume that $\epsilon$ is as small as needed.
\par
\bigskip
{\bf Remark:} The symbols used for the hopping parameters are standard. The reason why $\gamma_2$ was omitted is that it refers to next-to-nearest layer hopping in graphite.
In addition, for simplicity, we have neglected the intra-layer next-to-nearest neighbor hopping $\gamma_0'\approx 0.1\gamma_1$, which is known to play an analogous role to $\gamma_4$ and $\Delta$ [\cite{ZLe08}].
\par\bigskip
\begin{figure}
\hfil\includegraphics{Figs/hoppings.pdf}\par
\caption{The different types of hopping. From top-left to bottom-right: $a\leftrightarrow b$, $\tilde a\leftrightarrow\tilde b$, $a\leftrightarrow\tilde b$, $b\leftrightarrow\tilde a$, $a\leftrightarrow \tilde a$, $b\leftrightarrow\tilde b$. Atoms of type $a$ and $\tilde a$ are represented by spheres and those of type $b$ and $\tilde b$ by cubes; the interaction is represented by solid red (color online) cylinders; the interacting atoms are displayed either in purple or in blue.}\label{intergraph}\end{figure}
\indent The interactions between electrons will be taken to be of {\it extended Hubbard} form, i.e.
\begin{equation}
\mathcal H_I:=U\sum_{( x, y)}v( x- y)\left( n_{ x}-\frac{1}{2}\right)\left( n_{ y}-\frac{1}{2}\right)
\label{hamintxintro}\end{equation}
where $n_x:=\alpha_x^\dagger\alpha_x$ in which $\alpha_x$ is one of the annihilation operators $a_x$, $\tilde b_x$, $\tilde a_x$ or $b_x$; the sum over $(x,y)$ runs over all pairs of atoms in the lattice; $v$ is a short range interaction potential (exponentially decaying); $U$ is the {\it interaction strength} which we will assume to be small.\par
\bigskip
\indent We then define the {\it Gibbs average} as
$$\left<\cdot\right>:=\frac{1}{Z}\mathrm{Tr}\left( e^{-\beta\mathcal H}\cdot\right)$$
where
$$Z:=\mathrm{Tr}\left( e^{-\beta\mathcal H}\right)=:e^{-\beta|\Lambda| f}.$$
The physical quantities we will study here are the {\it free energy} $f$, and the {\it two-point Schwinger function} defined as the $4\times4$ matrix
\begin{equation}\check s_2(\mathbf x_1-\mathbf x_2):=\left(\left<\mathbf T(\alpha_{\mathbf x_1}'\alpha^\dagger_{\mathbf x_2})\right>\right)_{(\alpha',\alpha)\in\{a,\tilde b,\tilde a,b\}^2},
\label{schwindefintro}\end{equation}
where $\mathbf x_1:=(t_1,x_1)$ and $\mathbf x_2:=(t_2,x_2)$ includes an extra {\it imaginary time} component, $t_{1,2}\in[0,\beta)$, which is introduced in order to compute $Z$ and Gibbs averages,
$$\alpha_{t,x}:=e^{\mathcal H_0t}\alpha_{x}e^{-\mathcal H_0t}\quad\mathrm{for\ }\alpha\in\{a,\tilde b,\tilde a,b\}$$
and $\mathbf T$ is the {\it Fermionic time ordering operator}:
\begin{equation}
\mathbf T(\alpha'_{t_1,x_1}\alpha^\dagger_{t_2,x_2})=\left\{\begin{array}l\alpha_{t_1,x_1}'\alpha^\dagger_{t_2,x_2}\mathrm{\ if\ }t_1>t_2\\
-\alpha^\dagger_{t_2,x_2}\alpha'_{t_1,x_1}\mathrm{\ if\ }t_1\leq t_2\end{array}\right..
\label{timeorderdef}\end{equation}
We denote the Fourier transform of $\check s_2(\mathbf x)$ (or rather of its anti-periodic extension in imaginary time for $t$'s beyond $[0,\beta)$) by $s_2(\mathbf k)$ where $\mathbf k:=(k_0,k)$, and $k_0\in\frac{2\pi}{\beta}(\mathbb Z+\frac12)$.\par
\subsection{Non-interacting system}
\indent In order to state our main results on the interacting two-point Schwinger function,
it is useful to first review the scaling properties of the non-interacting one,
$$s_2^{(0)}(\mathbf k)=-(ik_0\mathds1+H_0(k))^{-1},$$
including a
discussion of the structure of its singularities in momentum space.
\bigskip
\point{\bf Non-interacting Fermi surface.} If $H_0(k)$ is not invertible, then $s_2^{(0)}(0,k)$ is divergent. The set of quasi-momenta $\mathcal F_0:=\{k,\ \det H_0(k)=0\}$ is called the non-interacting {\it Fermi surface}
at zero chemical potential, which has the following structure: it contains two isolated points located at
\begin{equation}p_{F,0}^{\omega}:=\left(\frac{2\pi}{3},\omega\frac{2\pi}{3\sqrt3}\right),\ \omega\in\{-1,+1\}\label{fermdefz}\end{equation}
around each of which there are three very small curves that are approximately elliptic (see figure~\ref{figferm}). The whole singular set $\mathcal F_0$ is contained within
two small circles (of radius $O(\epsilon^2)$), so that on scales larger than $\epsilon^2$, $\mathcal F_0$ can be approximated by just two points, $\{p_{F,0}^\pm\}$, see figure~\ref{figferm}. As we zoom in, looking at smaller
scales, we realize that each small circle contains four Fermi points: the central one, and three secondary points around it, called $\{p_{F,j}^\pm,\ j\in\{1,2,3\}\}$.
A finer zoom around the secondary points reveals that they are actually curves of size $\epsilon^3$. \par
\bigskip
\begin{figure}
\hfil\includegraphics{Figs/figferm.pdf}\par
\caption{Schematic representation of the Fermi points. Each dotted square represents a zoom into the finer structure of the Fermi points. The secondary Fermi points are labeled as indicated in the figure. In order not to clutter the drawing, only one of the zooms around the secondary Fermi point was represented.
}\label{figferm}\end{figure}\par
\point{\bf Non-interacting Schwinger function.} We will now make the statements about approximating the Fermi surface more precise, and discuss the behavior of $s_2^{(0)}$ around its singularities. We will identify four regimes in which the Schwinger function behaves differently.\par
\bigskip
\subpoint{\bf First regime.} One can show that, if $\mathbf p_{F,0}^\pm:=(0,p_{F,0}^\pm)$, and
$$\|(k_0,k_x',k_y') \|_{\mathrm I}:=\sqrt{k_0^2+(k_x')^2+(k_y')^2}$$
then
\begin{equation}
s_2^{(0)}(\mathbf p_{F,0}^\pm+\mathbf k' )=\left(\mathfrak L_{\mathrm{I}}\hat A(\mathbf p_{F,0}^\pm+\mathbf k' )\right)^{-1}\left(\mathds1+O(\|\mathbf k' \|_{\mathrm{I}},\epsilon\|\mathbf k' \|_{\mathrm{I}}^{-1})\right)
\label{freeschwino}\end{equation}
in which
$\mathfrak L_{\mathrm{I}}\hat A$ is a matrix, independent of $\gamma_1$, $\gamma_3$, $\gamma_4$ and $\Delta$, whose eigenvalues vanish {\it linearly} around $\mathbf p_{F,0}^\pm$ (see figure~\ref{figbands}b). We thus identify a {\it first regime}:
$$\epsilon\ll\|\mathbf k' \|_{\mathrm{I}}\ll1$$
in which the error term in~(\ref{freeschwino}) is {\it small}. In this first regime, $\gamma_1$, $\gamma_3$, $\gamma_4$ and $\Delta$ are negligible, and the Fermi surface is approximated by $\{p_{F,0}^\pm\}$, around which the Schwinger function diverges {\it linearly}.\par
\bigskip
\begin{figure}
\hfil\includegraphics{Figs/bands.pdf}\par
\caption{Eigenvalues of $H_0(k)$. The sub-figures b,c,d are finer and finer zooms around one of the Fermi points.}
\label{figbands}
\end{figure}
\subpoint{\bf Second regime.} Now, if
$$\|(k_0,k_x',k_y') \|_{\mathrm{II}}:=\sqrt{k_0^2+\frac{(k_x')^4}{\gamma_1^2}+\frac{(k_y')^4}{\gamma_1^2}}$$
then
\begin{equation}
s_2^{(0)}(\mathbf p_{F,0}^\pm+\mathbf k' )=\left(\mathfrak L_{\mathrm{II}}\hat A(\mathbf p_{F,0}^\pm+\mathbf k' )\right)^{-1}\left(\mathds1+O(\epsilon^{-1}\|\mathbf k' \|_{\mathrm{II}},\epsilon^{3/2}\|\mathbf k' \|_{\mathrm{II}}^{-1/2})\right)
\label{freeschwint}\end{equation}
in which $\mathfrak L_{\mathrm{II}}\hat A$ is a matrix, independent of $\gamma_3$, $\gamma_4$ and $\Delta$. Two of its eigenvalues vanish {\it quadratically} around $\mathbf p_{F,0}^\pm$ (see figure~\ref{figbands}c) and two are bounded away from 0. The latter correspond to {\it massive} modes, whereas the former to {\it massless} modes. We thus identify a {\it second regime}:
$$\epsilon^3\ll\|\mathbf k' \|_{\mathrm{II}}\ll\epsilon$$
in which $\gamma_3$, $\gamma_4$ and $\Delta$ are negligible, and the Fermi surface is approximated by $\{p_{F,0}^\pm\}$, around which the Schwinger function diverges {\it quadratically}.\par
\bigskip
\subpoint{\bf Third regime.} If $\mathbf p_{F,j}^\pm:=(0,p_{F,j}^\pm)$, $j=0,1,2,3$, and
$$\|(k_0,k'_{j,x},k'_{j,y})\|_{\mathrm{III}}:=\sqrt{k_0^2+\gamma_3^2(k'_{j,x})^2+\gamma_3^2(k'_{j,y})^2}$$
then
\begin{equation}
s_2^{(0)}(\mathbf p_{F,j}^\pm+\mathbf k'_{j})=\left(\mathfrak L_{\mathrm{III},j}\hat A(\mathbf p_{F,j}^\pm+\mathbf k'_{j})\right)^{-1}\left(\mathds1+O(\epsilon^{-3}\|\mathbf k'_{j}\|_{\mathrm{III}},\epsilon^{4}\|\mathbf k'_{j}\|_{\mathrm{III}}^{-1})\right)
\label{freeschwinth}\end{equation}
in which $\mathfrak L_{\mathrm{III},j}\hat A$ is a matrix, independent of $\gamma_4$ and $\Delta$, two of whose eigenvalues vanish {\it linearly} around $\mathbf p_{F,j}^\pm:=(0,p_{F,j}^\pm)$ (see figure~\ref{figbands}d) and two are bounded away from 0. We thus identify a {\it third regime}:
$$\epsilon^4\ll\|\mathbf k'_{j}\|_{\mathrm{III}}\ll\epsilon^3$$
in which $\gamma_4$ and $\Delta$ are negligible, and the Fermi surface is approximated by $\{p_{F,j}^\pm\}_{j\in\{0,1,2,3\}}$, around which the Schwinger function diverges {\it linearly}.
\par
\bigskip
{\bf Remark:} If $\gamma_4=\Delta=0$, then the error term $O(\epsilon^{4}\|\mathbf k'_{j}\|_{\mathrm{III}}^{-1})$ in (\ref{freeschwinth}) vanishes identically, which allows us to extend the third regime to all momenta satisfying
$$\|\mathbf k'_{j}\|_{\mathrm{III}}\ll\epsilon^3.$$
\subsection{Main Theorem}\label{maintheosec}
\indent We now state the Main Theorem, whose proof will occupy the rest of the paper. Roughly, our result is that as long as $|U|$ and $\epsilon$ are small enough and $\gamma_4=\Delta=0$ (see the remarks following the statement for an explanation of why this is assumed), the free energy and the two-point Schwinger function are well defined in the thermodynamic and zero-temperature limit $|\Lambda|,\beta\to\infty$, and that the two-point Schwinger function is analytically close to that with $U=0$. The effect of the interaction is shown to merely {\it renormalize} the constants of the non-interacting Schwinger function.\par
\bigskip
\indent We define
$$
\mathcal B_\infty:=\mathbb{R}\times\left(\mathbb{R}^2/(\mathbb ZG_1+\mathbb ZG_2)\right),\quad
G_1:=\left(\frac{2\pi}{3},\frac{2\pi}{\sqrt{3}}\right),\mathrm\ G_2:=\left(\frac{2\pi}{3},-\frac{2\pi}{\sqrt{3}}\right),
$$
where the physical meaning of $\mathbb{R}^2/(\mathbb ZG_1+\mathbb ZG_2)$ is that of the {\it first Brillouin zone}, and
$G_{1,2}$ are the generators of the dual lattice. \bigskip
\delimtitle{\bf Main Theorem}
If $\gamma_4=\Delta=0$, then there exists $U_0>0$ and $\epsilon_0>0$ such that for all $|U|<U_0$ and $\epsilon<\epsilon_0$, the specific ground state energy
$$e_0:=-\lim_{\beta\to\infty}\lim_{|\Lambda|\to\infty}\frac{1}{\beta|\Lambda|}\log(\mathrm{Tr}(e^{-\beta\mathcal H}))$$
exists and is analytic in $U$. In addition, there exist eight Fermi points $\{\tilde{\mathbf p}_{F,j}^\omega\}_{\omega=\pm,j=0,1,2,3}$ such that:
\begin{equation} \tilde{\mathbf p}_{F,0}^\omega={\mathbf p}_{F,0}^\omega, \qquad |\tilde{\mathbf p}_{F,j}^\omega-{\mathbf p}_{F,j}^\omega|\leqslant (\mathrm{const}.)\ |U|\epsilon^2, \ j=1,2,3,\end{equation}
and, $\forall\mathbf k\in\mathcal B_\infty\setminus\{\tilde{\mathbf p}_{F,j}^\omega\}_{\omega=\pm,j=0,1,2,3}$, the thermodynamic and zero-temperature limit of the two-point Schwinger function, $\lim_{\beta\to\infty}\lim_{|\Lambda|\to\infty}s_2(\mathbf k)$, exists and is analytic in $U$.
\endtheo
\bigskip
{\bf Remarks}:
\def\itemizepenalty{10000}
\begin{itemize}
\item The theorem requires $\gamma_4=\Delta=0$. As we saw above, those quantities play a negligible role in the non-interacting theory as long as we do not move beyond the third regime. This suggests that the theorem should hold with $\gamma_4,\Delta\neq0$ under the condition that $\beta$ is not too large,
i.e., smaller than $(\mathrm{const}.)\ \epsilon^{-4}$. However, that case presents a number of extra technical complications, which we will spare the reader.
\item The conditions that $|U|<U_0$ and $\epsilon<\epsilon_0$ are independent, in that we do not require any condition on the relative values of $|U|$ and $\epsilon$. Such a result calls for tight bounds on the integration over the first regime. If we were to assume that $|U|\ll\epsilon$, then the discussion would be greatly simplified, but such a condition would be artificial, and we will not require it be satisfied. L.~Lu~[\cite{Lu13}] sketched
the proof of a result similar to our Main Theorem, without discussing the first two regimes, which requires such an artificial condition on $U/\epsilon$.
The renormalization of the secondary Fermi points is also ignored in that reference.
\end{itemize}
\def\itemizepenalty{0}
\bigskip
\indent In addition to the Main Theorem, we will prove that the dominating part of the two point Schwinger function is qualitatively the same as the non-interacting one, with renormalized constants. This result is detailed in Theorems~\ref{theoo}, \ref{theot} and~\ref{theoth} below, each of which refers to one of the three regimes.\par
\bigskip
\point{\bf First regime.} Theorem~\ref{theoo} states that in the first regime, the two-point Schwinger function behaves at dominant order like the non-interacting one with renormalized factors.\par
\bigskip
\theo{Theorem}\label{theoo}
Under the assumptions of the Main Theorem, if $C\epsilon\leqslant \|\mathbf k-\mathbf p_{F,0}^\omega \|_{\mathrm I}\leqslant C^{-1}$ for a suitable $C>0$, then, in the thermodynamic and zero-temperature limit,
\begin{equation}
s_2(\mathbf k)=-\frac{1}{\tilde k_0\bar k_0+|\bar\xi|^2}\left(\begin{array}{*{4}{c}}
-i\bar k_0&0&0&\bar\xi^*\\
0&-i\bar k_0&\bar\xi&0\\
0&\bar\xi^*&-i\tilde k_0&0\\
\bar\xi&0&0&-i\tilde k_0\end{array}\right)
(\mathds1+r(\mathbf k))
\label{schwino}\end{equation}
where
\begin{equation} r(\mathbf p^\omega_{F,0}+\mathbf k')=O\big(\big(1+|U||\log\|\mathbf k' \|_{\mathrm{I}}|\big)\|\mathbf k' \|_{\mathrm{I}},\epsilon\|\mathbf k' \|_{\mathrm{I}}\big),\end{equation}
and, for $(k_0,k_{x}',k_{y}'):=\mathbf k-\mathbf p^\omega_{F,0}$,
\begin{equation}
\bar k_0:=z_1 k_0,\quad
\tilde k_0:=\tilde z_1 k_0,\quad
\bar \xi:=\frac{3}{2}v_1(ik'_x+\omega k'_y)
\label{rcco}\end{equation}
in which $(\tilde z_1,z_1,v_1)\in\mathbb{R}^3$ satisfy
\begin{equation}
|1-\tilde z_1|\leqslant C_1|U|,\quad
|1-z_1|\leqslant C_1|U|,\quad
|1-v_1|\leqslant C_1|U|
\label{ineqrcco}\end{equation}
for some constant $C_1>0$ (independent of $U$ and $\epsilon$).
\endtheo
\bigskip
{\bf Remarks}:
\def\itemizepenalty{10000}
\begin{itemize}
\item The singularities of $s_2$ are approached linearly in this regime.
\item By comparing (\ref{schwino}) with its non-interacting counterpart~(\ref{freepropzo}), we see that the effect of the interaction is to {\it renormalize} the constants in front of $k_0$ and $\xi$ in~(\ref{freepropzo}).
\item The {\it inter-layer correlations}, that is the $\{a,b\}\times\{\tilde a,\tilde b\}$ components of the dominating part of $s_2(\mathbf k)$ vanish. In this regime, the Schwinger function of bilayer graphene behave like that of two independent graphene layers.
\end{itemize}
\def\itemizepenalty{0}
\bigskip
\point{\bf Second regime.} Theorem~\ref{theot} states a similar result for the second regime. As was mentioned earlier, two of the components are {\it massive} in the second (and third) regime, and we first perform a change of variables to isolate them, and state the result on the massive and massless components, which are denoted below by $\bar s_2^{(M)}$ and $\bar s_2^{(m)}$ respectively.\par
\bigskip
\theo{Theorem}
\label{theot}
Under the assumptions of the Main Theorem, if $C\epsilon^3\leqslant \|\mathbf k-\mathbf p_{F,0}^\omega \|_{\mathrm{II}}\leqslant C^{-1}\epsilon$ for a suitable $C>0$, then, in the thermodynamic and zero-temperature limit,
\begin{equation}
s_2(\mathbf k)
=\left(\begin{array}{*{2}{c}}\mathds1&M(\mathbf k)^\dagger\\0&\mathds1\end{array}\right)\left(\begin{array}{*{2}{c}}\bar s_2^{(M)}&0\\
0&\bar s_2^{(m)}(\mathbf k)
\end{array}\right)\left(\begin{array}{*{2}{c}}\mathds1&0\\ M(\mathbf k)&\mathds1\end{array}\right)
(\mathds1+r(\mathbf k))
\label{schwint}\end{equation}
where:
\begin{equation}
r(\mathbf p_{F,0}^\omega+\mathbf k')=
O(\epsilon^{-1/2}\|\mathbf k' \|_{\mathrm{II}}^{1/2},
\epsilon^{3/2}\|\mathbf k' \|_{\mathrm{II}}^{-1/2},
|U|\epsilon\, |\log\epsilon|),
\label{errort}\end{equation}
\begin{equation}
\bar s_2^{(m)}(\mathbf k)
=\frac{1}{\bar\gamma_1^2\bar k_0^2+|\bar\xi|^4}
\left(\begin{array}{*{2}{c}}
i\bar\gamma_1^2\bar k_0&
\bar\gamma_1(\bar\xi^*)^2\\
\bar\gamma_1\bar\xi^2&
i\bar\gamma_1^2\bar k_0
\end{array}\right),\quad
\bar s_2^{(M)}=-\left(\begin{array}{*{2}{c}}0&\bar\gamma_1^{-1}\\\bar\gamma_1^{-1}&0\end{array}\right),
\label{schwintcomps}\end{equation}
\begin{equation}
M(\mathbf k):=-\frac1{\bar\gamma_1}\left(\begin{array}{*{2}{c}}\bar\xi^*&0\\0&\bar\xi\end{array}\right)
\label{theoMt}\end{equation}
and, for $(k_0,k_{x}',k_{y}'):=\mathbf k-\mathbf p^\omega_{F,0}$,
\begin{equation}
\bar\gamma_1:=\tilde m_2 \gamma_1,\quad
\bar k_0:=z_2 k_0,\quad
\bar \xi:=\frac{3}{2}v_2(ik'_x+\omega k'_y)
\label{rcct}\end{equation}
in which $(\tilde m_2,z_2,v_2)\in\mathbb{R}^3$ satisfy
\begin{equation}
|1-\tilde m_2|\leqslant C_2|U|,\quad
|1-z_2|\leqslant C_2|U|,\quad
|1-v_2|\leqslant C_2|U|
\label{ineqrcct}\end{equation}
for some constant $C_2>0$ (independent of $U$ and $\epsilon$).
\endtheo
\bigskip
{\bf Remarks}:
\def\itemizepenalty{10000}
\begin{itemize}
\item The {\it massless} components $\{\tilde a,b\}$ are left invariant under the change of basis that block-diagonalizes $s_2$. Furthermore, $M$ is {\it small} in the second regime, which implies that the {\it massive} components are {\it approximately} $\{a,\tilde b\}$.
\item As can be seen from~(\ref{schwintcomps}), the {\it massive} part $\bar s_2^{(M)}$ of $s_2$ is not singular in the neighborhood of the Fermi points, whereas the {\it massless} one, i.e. $\bar s_2^{(m)}$, is.
\item The massless components of $s_2$ approach the singularity {\it quadratically} in the spatial components of $\mathbf k$.
\item Similarly to the first regime, by comparing~(\ref{schwintcomps}) with~(\ref{freepropzt}), we find that the effect of the interaction is to {\it renormalize} constant factors.
\end{itemize}
\def\itemizepenalty{0}
\bigskip
\point{\bf Third regime.} Theorem~\ref{theoth} states a similar result as Theorem~\ref{theot} for the third regime, though the discussion is made more involved by the presence of the extra Fermi points.\par
\bigskip
\theo{Theorem}
\label{theoth}
For $j=0,1$, under the assumptions of the Main Theorem, if $\|\mathbf k-\tilde{\mathbf p}_{F,j}^\omega\|_{\mathrm{III}}\leqslant C^{-1}\epsilon^3$ for a suitable $C>0$, then
\begin{equation}
s_2(\mathbf k)
=\left(\begin{array}{*{2}{c}}\mathds1&M(\mathbf k)^\dagger\\0&\mathds1\end{array}\right)\left(\begin{array}{*{2}{c}}\bar s_2^{(M)}&0\\
0&\bar s_2^{(m)}(\mathbf k)
\end{array}\right)\left(\begin{array}{*{2}{c}}\mathds1&0\\ M(\mathbf k)&\mathds1\end{array}\right)
(\mathds1+r(\mathbf k))
\label{schwinth}\end{equation}
where
\begin{equation}
r(\tilde{\mathbf p}_{F,j}^\omega+\mathbf k'_{j})=O(
\epsilon^{-3}\|\mathbf k'_{j}\|_{\mathrm{III}}(1+\epsilon|\log\|\mathbf k'_{j}\|_{\mathrm{III}}||U|),
\epsilon(1+|\log\epsilon||U|)
)
\label{errorth}\end{equation}
\begin{equation}
\bar s_2^{(m)}(\mathbf k)
=\frac{1}{\bar k_{0,j}^2+\gamma_3^2|\bar x_j|^2}
\left(\begin{array}{*{2}{c}}
i\bar k_{0,j}&
\gamma_3\bar x_j^*\\
\gamma_3\bar x_j&
i\bar k_{0,j}
\end{array}\right),\quad
\bar s_2^{(M)}=-\left(\begin{array}{*{2}{c}}0&\bar\gamma_{1,j}^{-1}\\\bar\gamma_{1,j}^{-1}&0\end{array}\right),
\label{schwinthcomps}\end{equation}
\begin{equation}
M(\mathbf k):=-\frac1{\bar\gamma_{1,j}}\left(\begin{array}{*{2}{c}}\bar\Xi_j^*&0\\0&\bar\Xi_j\end{array}\right)
\label{theoMth}\end{equation}
and, for $(k_0,k_{x}',k_{y}'):=\mathbf k-\mathbf p^\omega_{F,j}$,
\begin{equation}\begin{array}c
\bar k_{0,j}:=z_{3,j}k_0,\quad
\bar\gamma_{1,j}:=\tilde m_{3,j}\gamma_1,\quad
\bar x_0:=\tilde v_{3,0}\frac{3}{2}(ik'_x-\omega k'_y)=:-\bar\Xi_0^*\\[0.5cm]
\bar x_{1}:=\frac{3}{2}\left(3\tilde v_{3,1}ik'_x+\tilde w_{3,1}\omega k'_x\right),\quad
\bar\Xi_{1}:=m_{3,1}\gamma_1\gamma_3+\bar v_{3,1}ik'_x+\bar w_{3,1}k'_y
\end{array}\label{rccthj}\end{equation}
in which $(\tilde m_{3,j},m_{3,j},z_{3,j},\bar v_{3,j},\tilde v_{3,j},\bar w_{3,j},\tilde w_{3,j})\in\mathbb{R}^7$ satisfy
\begin{equation}\begin{array}c
|m_{3,j}-1|+
|\tilde m_{3,j}-1|
\leqslant C_{3}|U|,\quad
|z_{3,j}-1|\leqslant C_{3}|U|,\\[0.5cm]
|\bar v_{3,j}-1|+
|\tilde v_{3,j}-1|\leqslant C_{3}|U|,\quad
|\bar w_{3,j}-1|+
|\tilde w_{3,j}-1|\leqslant C_{3}|U|
\end{array}\label{ineqrccthj}\end{equation}
for some constant $C_{3}>0$ (independent of $U$ and $\epsilon$).
\endtheo
\bigskip
Theorem~\ref{theoth} can be extended to the neighborhoods of $\tilde{\mathbf p}_{F,j}^\omega$ with $j=2,3$,
by taking advantage of the symmetry of the system under rotations of angle $2\pi/3$:\par
\medskip
\delimtitle{\bf Extension to $j=2,3$}
For $j=2,3$, under the assumptions of the Main Theorem, if $\|\mathbf k-\tilde{\mathbf p}_{F,j}^\omega\|_{\mathrm{III}}\leqslant C^{-1}\epsilon^3$ for a suitable $C>0$, then
\begin{equation}
s_2(\mathbf k'_{j}+\tilde{\mathbf p}_{F,j}^\omega)=
\left(\begin{array}{*{2}{c}}\mathds1&0\\0&\mathcal T_{T\mathbf k'_{j}+\tilde{\mathbf p}_{F,j-\omega}^\omega}\end{array}\right) s_2(T\mathbf k'_{j}+\tilde{\mathbf p}_{F,j-\omega}^{\omega})\left(\begin{array}{*{2}{c}}\mathds1&0\\0&\mathcal T^\dagger_{T\mathbf k'_{j}+\tilde{\mathbf p}_{F,j-\omega}^\omega}\end{array}\right)
\label{schwinthj}\end{equation}
where $T(k_0,k_x,k_y)$ denotes the rotation of the $k_x$ and $k_y$ components by an angle $2\pi/3$, $\mathcal T_{(k_0,k_x,k_y)}:=e^{-i(\frac32k_x-\frac{\sqrt3}2k_y)\sigma_3}$, and $\tilde{\mathbf p}_{F,4}^{-}\equiv \tilde{\mathbf p}_{F,1}^{-}$.
\endtheo
\bigskip
{\bf Remarks}:
\def\itemizepenalty{10000}
\begin{itemize}
\item The remarks below Theorem~\ref{theot} regarding the massive and massless fields hold here as well.
\item The massless components of $s_2$ approach the singularities {\it linearly}.
\item By comparing~(\ref{schwinth}) with (\ref{freepropzth}) and~(\ref{freepropoth}), we find that the effect of the interaction is to {\it renormalize} the constant factors.
\end{itemize}
\def\itemizepenalty{0}
\subsection{Sketch of the proof}
\indent In this section, we give a short account of the main ideas behind the proof of the Main Theorem.\par
\bigskip
\point{\bf Multiscale decomposition.} The proof relies on a {\it multiscale} analysis of the model, in which the free energy and Schwinger function are expressed as successive integrations over individual scales.
Each scale is defined as a set of $\mathbf k$'s contained inside an annulus at a distance of $2^h$ for $h\in\mathbb Z$ around the singularities located at $\mathbf p_{F,j}^\omega$.
The positive scales correspond to the ultraviolet regime, which we analyze in a multiscale fashion because of the (very mild) singularity of the free propagator at equal imaginary times.
It may be possible to avoid the decomposition by employing ideas in the spirit of [\cite{SS08}]. The negative scales are treated differently, depending on the regimes they belong to (see below),
and they contain the essential difficulties of the problem, whose nature is intrinsically infrared.
\par
\bigskip
\point{\bf First regime.} In the first regime, i.e. for $-1\gg h\gg h_\epsilon:=\log_2\epsilon$, the system behaves like two uncoupled graphene layers, so the analysis carried out in~[\cite{GM10}] holds. From a renormalization group perspective, this regime is {\it super-renormalizable}: the scaling dimension of diagrams with $2l$ external legs is $3-2l$, so that only the two-legged diagrams are relevant whereas all of the others are irrelevant (see section~\ref{powercountingsec} for precise definitions of scaling dimensions, relevance and irrelevance). This allows us to compute a strong bound on four-legged contributions:
$$|\hat W_4^{(h)}(\mathbf k)|\leqslant (\mathrm{const.})\ |U|2^{2h}$$
whereas a naive power counting argument would give $|U|2^{h}$ (recall that with our conventions $h$ is negative).\par
\bigskip
\indent The super-renormalizability in the first regime stems from the fact that the Fermi surface is 0-dimensional and that $H_0$ is linear around the Fermi points. While performing the multiscale integration, we deal with the two-legged terms by incorporating them into $H_0$, and one must therefore prove that by doing so, the Fermi surface remains 0-dimensional and that the singularity remains linear. This is guaranteed by a symmetry argument, which in particular shows the invariance of the Fermi surface.\par
\bigskip
\point{\bf Second regime.} In the second regime, i.e. for $3h_\epsilon\ll h\ll h_\epsilon$, the singularities of $H_0$ are quadratic around the Fermi points, which changes the {\it power counting} of the renormalization group analysis: the scaling dimension of $2l$-legged diagrams becomes $2-l$ so that the two-legged diagrams are still relevant, but the four-legged ones become marginal. One can then check~[\cite{Va10}] that they are actually marginally relevant, which means that their contribution increases proportionally to $|h|$. This turns out not to matter: since the second regime is only valid for $h\gg 3h_\epsilon$, $|\hat W_4^{(h)}|$ may only increase by $3|h_\epsilon|$, and since the theory is super-renormalizable in the first regime, there is an extra factor $2^{h_\epsilon}$ in $\hat W_4^{(h_\epsilon)}$, so that $\hat W_4^{(h)}$ actually increases from $2^{h_\epsilon}$ to $3|h_\epsilon|2^{h_\epsilon}$, that is to say it barely increases at all if $\epsilon$ is small enough.\par
\bigskip
\indent Once this essential fact has been taken into account, the renormalization group analysis can be carried out without major difficulties. As in the first regime, the invariance of the Fermi surface is guaranteed by a symmetry argument.\par
\bigskip
\point{\bf Third regime.} In the third regime, i.e. for $h\ll 3h_\epsilon$, the theory is again super-renormalizable (the scaling dimension is $3-2l$). There is however an extra difficulty with respect to the first regime, in that the Fermi surface is no longer invariant under the renormalization group flow, but one can show that it does remain 0-dimensional, and that the only effect of the multiscale integration is to move $p_{F,j}^\omega$ along the line between itself and $p_{F,0}^\omega$.\par
\subsection{Outline}
\indent The rest of this paper is devoted to the proof of the Main Theorem and of Theorems~\ref{theoo}, \ref{theot} and~\ref{theoth}. The sections are organized as follows.
\begin{itemize}
\item In section~\ref{themodelsec}, we define the model in a more explicit way than what has been done so far; then we show how to compute the free energy and Schwinger function using a Fermionic path integral formulation and a {\it determinant expansion}, due to Battle, Brydges and Federbush~[\cite{BF78}, \cite{BF84}], see also [\cite{BK87}, \cite{AR98}]; and finally we discuss the symmetries of the system.
\item In section~\ref{proppropsec}, we discuss the non-interacting system. In particular, we derive detailed formulae for the Fermi points and for the asymptotic behavior of the propagator around its singularities.
\item In section~\ref{schemesec}, we describe the multiscale decomposition used to compute the free energy and Schwinger function.
\item In section~\ref{treeexpsec}, we state and prove a {\it power counting} lemma, which will allow us to compute bounds for the effective potential in each regime. The lemma is based on the Gallavotti-Nicol\`o tree expansion~[\cite{GN85}], and follows~[\cite{BG90}, \cite{GM01}, \cite{Gi10}]. We conclude this section by showing how to compute the two-point Schwinger function from the effective potentials.
\item In section~\ref{uvsec}, we discuss the integration over the {\it ultraviolet regime}, i.e. scales $h>0$.
\item In sections~\ref{osec}, \ref{tsec} and~\ref{thsec}, we discuss the multiscale integration in the first, second and third regimes,
and complete the proofs of the Main Theorem, as well as of Theorems \ref{theoo}, \ref{theot}, \ref{theoth}.
\end{itemize}
\section{The model}
\label{themodelsec}
\hfil\framebox{\bf From this point on, we set $\gamma_4=\Delta=0$.}
\bigskip
\indent In this section, we define the model in precise terms, re-express the free energy and two-point Schwinger function in terms of Grassmann integrals and truncated expectations, which we will subsequently explain how to compute, and discuss the symmetries of the model and their representation in this formalism.\par
\subsection{Precise definition of the model}
\label{modelsec}
\indent In the following, some of the formulae are repetitions of earlier ones, which are recalled for ease of reference. This section complements section~\ref{intromodelsec},
where the same definitions were anticipated in a less verbose form. The main novelty lies in the momentum-real space correspondence, which is
made explicit.
\par
\bigskip
\point{\bf Lattice.} As mentioned in section~\ref{introduction}, the atomic structure of bilayer graphene consists in two honeycomb lattices in so-called {\it Bernal} or {\it AB} stacking, as was shown in figure~\ref{baregraph}. It 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},0\right),\mathrm{\ } l_2:=\left(\frac{3}{2},-\frac{\sqrt{3}}{2},0\right)
\label{laeo}\end{equation}
where we have chosen the unit length to be equal to the distance between two nearest neighbors in a layer (see figure~\ref{cellgraph}). The elementary cell consists of four atoms at the following coordinates
$$(0,0,0);\mathrm{\ }(0,0,c);\mathrm{\ }(-1,0,c);\mathrm{\ }(1,0,0)$$
given relatively to the center of the cell. $c$ is the spacing between layers; it can be measured experimentally, and has a value of approximately 2.4~[\cite{TMe92}].\par
\bigskip
\begin{figure}
\hfil\hskip-1.5cm\includegraphics{Figs/cellgraph.pdf}\par
\label{cellgraph}
\caption{decomposition of the crystal into elementary cells, represented by the blue (color online) rhombi. There are four atoms in each elementary cell: \ta of type $a$ at $(0,0,0)$, \tbt of type $\tilde b$ at $(0,0,c)$, \tat of type $\tilde a$ at $(-1,0,c)$ and \tb of type $b$ at $(0,0,c)$.}
\end{figure}\par\medskip
\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 intra-layer nearest neighbor vectors:
\begin{equation}
\delta_1:=(1,0,0),\quad\delta_2:=\left(-\frac{1}{2},\frac{\sqrt{3}}{2},0\right),\quad\delta_3:=\left(-\frac{1}{2},-\frac{\sqrt{3}}{2},0\right).
\label{lea}\end{equation}\par
\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}},0\right),\mathrm\ G_2=\left(\frac{2\pi}{3},-\frac{2\pi}{\sqrt{3}},0\right).
\label{laeG}\end{equation}
It is defined in such a way that $\forall x\in\Lambda$, $\forall k\in\hat\Lambda$,
$$e^{ikxL}=1.$$
Since the third component of vectors in $\hat\Lambda$ is always 0, we shall drop it and write vectors of $\hat\Lambda$ as elements of $\mathbb{R}^2$. 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}.
\par
\point{\bf Hamiltonian.} Given $ x\in\Lambda$, we denote the Fermionic annihilation operators at atoms of type $a$, $\tilde b$, $\tilde a$ and $b$
within the elementary cell centered at $x$ respectively by $a_{ x}$, $\tilde b_{ x}$, $\tilde a_{ x-\delta_1}$ and $ b_{ x+\delta_1}$. The corresponding creation operators are their adjoint operators.\par
\bigskip
\indent We recall the Hamiltonian~(\ref{fullham})
$$\mathcal H=\mathcal H_0+\mathcal H_I$$
where $\mathcal H_0$ is the {\it free Hamiltonian} and $\mathcal H_I$ is the {\it interaction Hamiltonian}.\par
\bigskip
\subpoint{\bf Free Hamiltonian.} As was mentioned in section~\ref{introduction}, the free Hamiltonian describes the {\it hopping} of electrons from one atom to another. Here, we only consider the hoppings $\gamma_0,\gamma_1,\gamma_3$, see figure~\ref{intergraph}, so that $\mathcal H_0$ has the following expression in $x$ space:
\begin{equation}\begin{array}{r@{\ }>{\displaystyle}l}
\mathcal H_0:=&-\gamma_0\sum_{\displaystyle\mathop{\scriptstyle x\in\Lambda}_{j=1,2,3}}\left( a_{ x}^\dagger b_{ x +\delta_j}+b_{ x +\delta_j}^\dagger a_{ x}+
\tilde b_x^\dagger \tilde a_{x-\delta_j} +\tilde a_{x-\delta_j}^\dagger\tilde b_x\right)
-\gamma_1\sum_{ x\in\Lambda}\left( a_{ x}^\dagger \tilde b_{ x}+\tilde b_{ x}^\dagger a_{ x}\right)
\\[0.2cm]
&-\gamma_3\sum_{\displaystyle\mathop{\scriptstyle x\in\Lambda}_{j=1,2,3}}\left( \tilde a_{ x-\delta_1}^\dagger b_{ x-\delta_1-\delta_j}+b_{ x-\delta_1-\delta_j}^\dagger\tilde a_{ x-\delta_1}\right)\\[0.2cm]
\end{array}\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}:=\sum_{x\in\Lambda}e^{ikx}a_{x}\;,\quad
\hat{\tilde b}_{k}:=\sum_{x\in\Lambda}e^{ikx}\tilde b_{x+\delta_1}\;,\quad
\hat{\tilde a}_{k}:=\sum_{x\in\Lambda}e^{ikx}\tilde a_{x-\delta_1}\;,\quad
\hat b_{k}:=\sum_{x\in\Lambda}e^{ikx}b_{x+\delta_1}\;\end{equation}
in terms of which
\begin{equation}
\mathcal H_0=-\frac{1}{|\Lambda|}\sum_{ k\in\hat\Lambda}\hat A_{ k}^\dagger H_0( k)A_{ k}\label{hamk}
\end{equation}
where $|\Lambda|=L^2$, $\hat A_k$ is a column vector, whose transpose is $\hat A_k^T=(\hat a_{ k},\hat{\tilde b}_{ k},\hat{\tilde a}_{ k},\hat{b}_{ k})$,
\begin{equation}H_0( k):=
\left(\begin{array}{*{4}{c}}
0&\gamma_1&0&\gamma_0\Omega^*(k)\\
\gamma_1&0&\gamma_0\Omega(k)&0\\
0&\gamma_0\Omega^*(k)&0&\gamma_3\Omega(k)e^{3ik_x}\\
\gamma_0\Omega(k)&0&\gamma_3\Omega^*(k)e^{-3ik_x}&0
\end{array}\right)\label{hmat}\end{equation}
and
$$
\Omega( k):=\sum_{j=1}^3e^{ik(\delta_j-\delta_1)}=1+2e^{-i\frac32k_x}\cos\left(\frac{\sqrt{3}}2k_y\right).
$$
We pick the energy unit in such a way that $\gamma_0=1$.
\bigskip
\subpoint{\bf Interaction.} We now define the interaction Hamiltonian. We first define the number operators $n_x^\alpha$ for
$\alpha\in\{a,\tilde b,\tilde a,b\}$ and $x\in\Lambda$ in the following way:
\begin{equation} n_x^a=a^\dagger_x a_x\;,\quad
n^{\tilde b}_x=\tilde b_{ x}^\dagger \tilde b_x\;,\quad
n^{\tilde a}_x=\tilde a_{ x-\delta_1}^\dagger \tilde a_{ x-\delta_1}\;,\quad
n^b_x=b_{ x+\delta_1}^\dagger b_{ x+\delta_1}\;\end{equation}
and postulate the form of the interaction to be of an extended {\it Hubbard} form:
\begin{equation}
\mathcal H_I:=U\sum_{( x, y)\in\Lambda^2}\sum_{(\alpha,\alpha')\in\{a,\tilde b,\tilde a,b\}^2}v( x+d_\alpha- y-d_{\alpha'})\left( n^\alpha_{ x}-\frac{1}{2}\right)\left( n^{\alpha'}_{ y}-\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_{\tilde b}:=(0,0,c),\ d_{\tilde a}:=(0,0,c)-\delta_1,\ d_b:=\delta_1$ and $v$ is a bounded, rotationally invariant function, which decays exponentially fast to zero at infinity. In our spin-less case, we can assume without loss of generality that $v(0)=0$.
\subsection{Schwinger function as Grassmann integrals and expectations}
\label{tracespathintssec}
\indent The aim of the present work is to compute the {\it specific free energy} and the {\it two-point Schwinger function}. These quantities are defined for finite lattices 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 inverse temperature and
\begin{equation}
\check s_{\alpha',\alpha}(\mathbf x_1-\mathbf x_2):=\left<\mathbf T(\alpha_{\mathbf x_1}'\alpha^\dagger_{\mathbf x_2})\right>:=\frac{\mathrm{Tr}(e^{-\beta\mathcal H}\mathbf T(\alpha'_{\mathbf x_1}\alpha_{\mathbf x_2}^\dagger))}{\mathrm{Tr}(e^{-\beta\mathcal H})}
\label{schwindef}\end{equation}
in which $(\alpha,\alpha')\in\mathcal A^2:=\{a,\tilde b,\tilde a,b\}^2$; $\mathbf x_{1,2}=(t_{1,2},x_{1,2})$ with $t_{1,2}\in[0,\beta)$; $\alpha_{\mathbf x}=e^{\mathcal Ht}\alpha_xe^{-\mathcal Ht}$; and $\mathbf T$ is the {\it Fermionic time ordering operator} defined in~(\ref{timeorderdef}). Our strategy essentially consists in deriving convergent expansions for $f_\Lambda$ and $\check s$, uniformly in $|\Lambda|$ and $\beta$, and then to take $\beta,|\Lambda|\to\infty$.\par
\bigskip
\indent However, the quantities on the right side of~(\ref{freeen}) and~(\ref{schwindef}) are somewhat difficult to manipulate. In this section, we will re-express $f_\Lambda$ and $\check s$ in terms of {\it Grassmann integrals} and {\it expectations}, and show how such quantities can be computed using a {\it determinant expansion}. This formalism will lay the groundwork for the procedure which will be used in the following to express $f_\Lambda$ and $\check s$ as series, and subsequently prove their convergence.\par
\bigskip
\point{\bf Grassmann integral formulation.} We first describe how to express~(\ref{freeen}) and~(\ref{schwindef}) as Grassmann integrals. The procedure is well known and details can be found in many references, see e.g. [\cite{GM10}, appendix~B] and [\cite{Gi10}] for a discussion adapted to the case of graphene, or [\cite{GM01}] for a discussion adapted to general low-dimensional Fermi systems,
or [\cite{BG95}] and [\cite{Sa13}] and references therein for an even more general picture. \par
\bigskip
\subpoint{\bf Definition.} We first define a Grassmann algebra and an integration procedure on it. We move to Fourier space: for every $\alpha\in\mathcal A:=\{a,\tilde b,\tilde a,b\}$, the operator $\alpha_{(t,x)}$ is associated
$$\hat\alpha_{\mathbf k=(k_0,k)}:=\frac{1}{\beta}\int_0^\beta dt\ e^{itk_0}e^{\mathcal H_0 t}\hat\alpha_ke^{-\mathcal H_0 t}$$
with $k_0\in2\pi\beta^{-1}(\mathbb Z+1/2)$ (notice that because of the $1/2$ term, $k_0\neq0$ for finite $\beta$). We notice that $\mathbf k\in\mathcal B_{\beta,L}:=(2\pi\beta^{-1}(\mathbb Z+1/2))\times\hat\Lambda$ varies in an infinite set. Since this will cause trouble when defining Grassmann integrals, we shall impose a cutoff $M\in\mathbb N$: let $\chi_0(\rho)$ be a smooth compact support function that returns $1$ if $\rho\leqslant 1/3$ and $0$ if $\rho\geqslant 2/3$, and let
$${\mathcal B}_{\beta,L}^*:=\mathcal B_{\beta,L}\cap\{(k_0,k),\ \chi_0(2^{-M}|k_0|\neq0)\}.$$
To every $(\hat\alpha_{\mathbf k},\hat\alpha_{\mathbf k}^\dagger)$ for $\alpha\in\mathcal A$ and $\mathbf k\in\mathcal B_{\beta,L}^*$, we associate a pair of {\it Grassmann variables} $(\hat\psi_{\mathbf k,\alpha}^-,\hat\psi_{\mathbf k,\alpha}^+)$, and we consider the finite Grassmann algebra (i.e. an algebra in which the $\hat\psi$ anti-commute with each other) generated by the collection $\{\hat\psi_{\mathbf k,\alpha}^\pm\}_{\mathbf k\in\mathcal B_{\beta,L}^*}^{\alpha\in\mathcal A}$. We define the Grassmann integral
$$\int
\prod^{\alpha\in{\mathcal A}}_{{\bf k}\in{\mathcal B}^*_{\beta,L}}d\hat\psi_{{\bf k},\alpha}^+ d\hat\psi_{{\bf k},\alpha}^-$$
as the linear operator on the Grassmann algebra whose action on a monomial in the variables $\hat\psi^\pm_{{\bf k},\alpha}$ is $0$ except if said monomial is $\prod^{\alpha\in{\mathcal A}}_{{\bf k}\in{\mathcal B}_{\beta,L}^*} \hat\psi^-_{{\bf k},\alpha} \hat\psi^+_{{\bf k},\alpha}$ up to a permutation of the variables, in which case the value of the integral is determined using
\begin{equation} \int
\prod^{\alpha\in{\mathcal A}}_{{\bf k}\in{\mathcal B}_{\beta,L}^*}
d\hat\psi_{{\bf k},\alpha}^+
d\hat\psi_{{\bf k},\alpha}^-
\left(\prod_{{\bf k}\in{\mathcal B}_{\beta,L}^*}^{\alpha\in{\mathcal A}}
\hat\psi^-_{{\bf k},\alpha}
\hat\psi^+_{{\bf k},\alpha}\right)=1
\label{2.1}\end{equation}
along with the anti-commutation of the $\hat\psi$.\par
\bigskip
\indent In the following, we will express the free energy and Schwinger function as {\it Grassmann integrals}, specified by a {\it propagator} and a {\it potential}.
The propagator is a $4\times4$ complex matrix $\hat g(\mathbf k)$, supported on some set $\mathcal B\subset\mathcal B_{\beta,L}^*$, and is associated with
the {\it Gaussian Grassmann integration measure}
\begin{equation}
P_{\hat g}(d\psi) := \left(\prod_{\mathbf k\in\mathcal B}
(\beta|\Lambda|)^{4}\det\hat g(\mathbf k)
\left(\prod_{\alpha\in\mathcal A}d\hat\psi_{\mathbf k,\alpha}^+d\hat\psi_{\mathbf k,\alpha}^-\right)\right)
\exp\left(-\frac{1}{\beta|\Lambda|}\sum_{\mathbf k\in\mathcal B}\hat\psi^{+}_{\mathbf k}\hat g^{-1}(\mathbf k)\hat\psi^{-}_{\mathbf k}\right).
\label{grassgauss}\end{equation}
Gaussian Grassmann integrals satisfy the following {\it addition principle}: given two propagators $\hat g_1$ and $\hat g_2$, and any polynomial $\mathfrak P(\psi)$ in the Grassmann variables,
\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
\subpoint{\bf Free energy.} We now express the free energy as a Grassmann integral. We define the {\it free propagator}
\begin{equation}
\hat g_{\leqslant M}({\bf k}):=\chi_0(2^{-M}|k_0|) (-ik_0\mathds 1-H_0(k))^{-1}
\label{freeprop}\end{equation}
and the Gaussian integration measure $P_{\leqslant M}(d\psi)\equiv P_{\hat g_{\leqslant M}}(d\psi)$.
One can prove (see e.g. [\cite{GM10}, appendix~B]) that if
\begin{equation} \frac{1}{\beta|\Lambda|}\log\int P_{\leqslant M}(d\psi)\ e^{-\mathcal V(\psi)}\label{log}\end{equation}
is analytic in $U$, uniformly as $M\to\infty$, a fact we will check a posteriori, then the finite volume free energy can be written as
\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, using the symbol $\int d\mathbf x$ as a shorthand for $\int_{0}^\beta dt \sum_{x\in \Lambda}$,
\begin{equation}
\mathcal V(\psi)=U\sum_{(\alpha,\alpha')\in{\mathcal A}^2}\int d\mathbf xd\mathbf y\ w_{\alpha,\alpha'}( \mathbf x-\mathbf y) \psi^+_{\mathbf x,\alpha}\psi^{-}_{\mathbf x,\alpha}
\psi^+_{\mathbf y,\alpha'}\psi^{-}_{\mathbf y,\alpha'}
\label{2.3a}\end{equation}
in which $w_{\alpha,\alpha'}(\mathbf x):=\delta(x_0)v(x+d_\alpha-d_{\alpha'})$, where $\delta(x_0)$ denotes the $\beta$-periodic Dirac delta function,
and
\begin{equation} \psi^\pm_{\mathbf x,\alpha}:=\frac1{\beta|\Lambda|}\sum_{{\bf k}\in{\mathcal B}^*_{\beta,L}}\hat \psi^{\pm}_{{\bf k},\alpha}e^{\pm i{\bf k}\mathbf x}\;.\end{equation}
Notice that the expression of $\mathcal V(\psi)$ in~(\ref{2.3a}) is very similar to that of $\mathcal H_I$, with an added imaginary time $(x_0,y_0)$ and the $\alpha_{\mathbf x}$ replaced by $\psi_{\mathbf x,\alpha}$, except that $(\alpha_\mathbf x^\dagger\alpha_\mathbf x-1/2)$ becomes $\psi_{\mathbf x,\alpha}^+\psi_{\mathbf x,\alpha}^-$. Roughly, the reason why we ``drop the 1/2'' is because of the difference between the anti-commutation rules of $\alpha_\mathbf x$ and $\psi_{\mathbf x,\alpha}$ (i.e., $\{\alpha_\mathbf x,\alpha_\mathbf x^\dagger\}=1$, vs. $\{\psi_{\mathbf x,\alpha}^+,\psi_{\mathbf x,\alpha}^-\}=0$).
More precisely, taking $\mathbf x=(x_0,x)$ with $x_0\in(-\beta,\beta)$,
it is easy to check that the limit as $M\to\infty$ of $g_{\leqslant M}(\mathbf x):=\int P_{\leqslant M}(d\psi)\psi^-_\mathbf x\psi^+_{\bf 0}$ is equal to $\check s(\mathbf x)$, if $\mathbf x\neq {\bf 0}$,
and equal to $\check s({\bf 0})+1/2$, otherwise. This extra $+1/2$ accounts for the ``dropping of the 1/2" mentioned above.\par
\bigskip
\subpoint{\bf Two-point Schwinger function.} The two-point Schwinger function can be expressed as a Grassmann integral as well: under the condition that
\begin{equation}
\frac{\int P_{\leqslant M}(d\psi)\ e^{-\mathcal V(\psi)}\hat \psi_{\mathbf k,\alpha_1}^-\hat \psi_{\mathbf k,\alpha_2}^+}{\int P_{\leqslant M}(d\psi)\ e^{-\mathcal V(\psi)}}\label{sc}
\end{equation}
is analytic in $U$ uniformly in $M$, a fact we will also check a posteriori, then one can prove (see e.g. [\cite{GM10}, appendix~B]) that the two-point Schwinger function can be written as
\begin{equation}
s_{\alpha_1,\alpha_2}(\mathbf k)=\lim_{M\to\infty}\frac{\int P_{\leqslant M}(d\psi)\ e^{-\mathcal V(\psi)}\hat \psi_{\mathbf k,\alpha_1}^-\hat \psi_{\mathbf k,\alpha_2}^+}{\int P_{\leqslant M}(d\psi)\ e^{-\mathcal V(\psi)}}.
\label{schwingrass}\end{equation}
In order to facilitate the computation of the right side of~(\ref{schwingrass}), we will first rewrite it as
\begin{equation}
s_{\alpha_1,\alpha_2}(\mathbf k)
=\lim_{M\to\infty}\int d\hat J_{\mathbf k,\alpha_1}^-d\hat J_{\mathbf k,\alpha_2}^+\log\int P_{\leqslant M}(d\psi)e^{-\mathcal V(\psi)+\hat J_{\mathbf k,\alpha_1}^+\hat\psi_{\mathbf k,\alpha_1}^-+\hat\psi_{\mathbf k,\alpha_2}^+\hat J_{\mathbf k,\alpha_2}^-}
\label{schwinform}\end{equation}
where $\hat J_{\mathbf k,\alpha}^-$ and $\hat J_{\mathbf k,\alpha'}^+$ are extra Grassmann variables introduced for the purpose of the computation (note here that the Grassmann integral over the variables $\hat J_{\mathbf k,\alpha_1}^-,\hat J_{\mathbf k,\alpha_2}^+$ acts as a functional derivative with respect to the same variables, due to the Grassmann integration/derivation rules).
We define the {\it generating functional}
\begin{equation}
\mathcal W(\psi,\hat J_{\mathbf k,\underline\alpha}):=\mathcal V(\psi)-\hat J_{\mathbf k,\alpha_1}^+\hat\psi_{\mathbf k,\alpha_1}^--\hat\psi_{\mathbf k,\alpha_2}^+\hat J_{\mathbf k,\alpha_2}^-.
\label{Weffpotdef}\end{equation}
\bigskip
\point{\bf Expectations.} We have seen that the free energy and Schwinger function can be computed as Grassmann integrals, it remains to see how one computes such integrals.
We can write (\ref{log}) as
\begin{equation}
\log\int P_{\leqslant M}(d\psi)e^{-\mathcal V(\psi)}=\sum_{N=1}^\infty\frac{(-1)^N}{N!}\mathcal E_{\leqslant M}^T(\underbrace{\mathcal V,\cdots,\mathcal V}_{N\mathrm{\ times}})=:\sum_{N=1}^\infty\frac{(-1)^N}{N!}\mathcal E_{\leqslant M}^T(\mathcal V;N).
\label{grasstruncexp}\end{equation}
where the {\it truncated expectation} is defined as
\begin{equation}\mathcal E_{\leqslant M}^T(\mathcal V_1,\cdots,\mathcal V_N):=\left.\frac{\partial^N}{\partial\lambda_1\cdots\partial\lambda_N}\log\int P_{\leqslant M}(d\psi)\ e^{\lambda_1\mathcal V_1+\cdots+\lambda_N\mathcal V_N}\right|_{\lambda_1=\cdots=\lambda_N=0}\label{truncexpdef}\;.\end{equation}
in which $(\mathcal V_1,\cdots,\mathcal V_N)$ is a collection of commuting polynomials and the index $_{\leqslant M}$ refers to the propagator of $P_{\leqslant M}(d\psi)$. A similar formula holds for~(\ref{sc}).\par
\bigskip
\indent The purpose of this rewriting is that we can compute truncated expectations in terms of a {\it determinant expansion}, also known as the Battle-Brydges-Federbush formula~[\cite{BF78}, \cite{BF84}], which expresses it as the determinant of a Gram matrix. The advantage of this writing is that, provided we first re-express the propagator $\hat g_{\leqslant M}(\mathbf k)$ in $\mathbf x$-space, the afore-mentioned Gram matrix can be bounded effectively (see section~\ref{powercountingsec}). We therefore first define an $\mathbf x$-space representation for $\hat g(\mathbf k)$:
\begin{equation}
g_{\leqslant M}(\mathbf x):=\frac{1}{\beta|\Lambda|}\sum_{\mathbf k\in\mathcal B^*_{\beta,L}}e^{i\mathbf k\cdot\mathbf x}\hat g_{\leqslant M}(\mathbf k).
\label{freepropx}\end{equation}
The determinant expansion is given in the following lemma, the proof of which can be found in [\cite{GM01}, appendix~A.3.2], [\cite{Gi10}, appendix~B].\par
\bigskip
\theo{Lemma}\label{detexplemma}
\indent Consider a family of sets $\mathbf P=(P_1,\cdots,P_s)$ where every $P_j$ is an ordered collection of Grassmann variables, we denote the product of the elements in $P_j$ by $\Psi_{P_j}:=\prod_{\psi\in P_j}\psi$.\par
\smallskip
\indent We call a pair $(\psi_{\mathbf x,\alpha}^-,\psi_{\mathbf x',\alpha'}^+)\in\mathbf P^2$ a {\it line}, and define the set of {\it spanning trees} $\mathbf T(\mathbf P)$ as the set of collections $T$ of lines that are such that upon drawing a vertex for each $P_i$ in $\mathbf P$ and a line between the vertices corresponding to $P_i$ and to $P_j$ for each line $(\psi_{\mathbf x,\alpha}^-,\psi_{\mathbf x',\alpha'}^+)\in T$ that is such that $\psi_{\mathbf x,\alpha}^-\in P_i$ and $\psi_{\mathbf x',\alpha'}^+\in P_j$, the resulting graph is a tree that connects all of the vertices.\par
\smallskip
\indent For every spanning tree $T\in\mathbf T(\mathbf P)$, to each line $l=(\psi_{\mathbf x,\alpha}^-,\psi_{\mathbf x',\alpha'}^+)\in T$ we assign a {\it propagator} $g_l:=g_{\alpha,\alpha'}(\mathbf x-\mathbf x')$.\par
\smallskip
\indent If $\mathbf P$ contains $2(n+s-1)$ Grassmann variables, with $n\in\mathbb N$, then there exists a probability measure $dP_T(\mathbf t)$ on the set of $n\times n$ matrices of the form $\mathbf t=M^TM$ with $M$ being a matrix whose columns are unit vectors of $\mathbb R^n$, such that
\begin{equation}\mathcal E_{\leqslant M}^T(\Psi_{P_1},\cdots,\Psi_{P_s})=\sum_{T\in\mathbf T(\mathbf P)}\sigma_T\prod_{l\in T}g_l\int dP_{T}(\mathbf t)\ \det G^{(T)}(\mathbf t)\label{detexp}\end{equation}
where $\sigma_T\in\{-1,1\}$ and $G^{(T)}(\mathbf t)$ is an $n\times n$ complex matrix each of whose components is indexed by a line $l\not\in T$ and is given by
$$G^{(T)}_{l}(\mathbf t)=\mathbf t_lg_l$$
(if $s=1$, then $\mathbf T(\mathbf P)$ is empty and both the sum over $T$ and the factor $\sigma_T\prod_{l\in T}g_l$ should be dropped from the right side of~(\ref{detexp})).
\endtheo
\bigskip
\indent Lemma~\ref{detexplemma}\ gives us a formal way of computing the right side of~(\ref{grasstruncexp}). However, proving that this formal expression is correct, in the sense that it is not divergent, will require a control over the quantities involved in the right side of~(\ref{detexp}), namely the propagator $g_{\leqslant M}$.
Since, as was discussed in the introduction, $g_{\leqslant M}$ is singular, controlling the right side of~(\ref{grasstruncexp}) is a non-trivial task that will require a multiscale analysis described in section~\ref{schemesec}.\par
\subsection{Symmetries of the system}
\label{symsec}
\indent In the following, we will rely heavily on the symmetries of the system, whose representation in terms of Grassmann variables is now discussed.\par
\bigskip
\indent A {\it symmetry} of the system is a map that leaves {\it both}
\begin{equation}
h_0:=\sum_{\mathbf x,\mathbf y}\psi^+_{\mathbf x}g^{-1}(\mathbf x-\mathbf y)\psi_{\mathbf y}^-
\label{hzdef}\end{equation}
and $\mathcal V(\psi)$ invariant ($\mathcal V(\psi)$ was defined in (\ref{2.3a})). We define
\begin{equation}
\hat\xi_{\mathbf k}^+:=\left(\begin{array}{*{2}{c}}\hat\psi_{\mathbf k,a}^+&\hat\psi_{\mathbf k,\tilde b}^+\end{array}\right),\quad
\hat\xi_{\mathbf k}^-:=\left(\begin{array}{c}\hat\psi_{\mathbf k,a}^-\\\hat\psi_{\mathbf k,\tilde b}^-\end{array}\right),\quad
\hat\phi_{\mathbf k}^+:=\left(\begin{array}{*{2}{c}}\hat\psi_{\mathbf k,\tilde a}^+&\hat\psi_{\mathbf k,b}^+\end{array}\right)
\hat\phi_{\mathbf k}^-:=\left(\begin{array}{c}\hat\psi_{\mathbf k,\tilde a}^-\\\hat\psi_{\mathbf k,b}^-\end{array}\right)
\label{axiphidef}\end{equation}
as well as the Pauli matrices
$$\sigma_1:=\left(\begin{array}{*{2}{c}}0&1\\1&0\end{array}\right),\quad\sigma_2:=\left(\begin{array}{*{2}{c}}0&-i\\i&0\end{array}\right),\quad\sigma_3:=\left(\begin{array}{*{2}{c}}1&0\\0&-1\end{array}\right).$$
We now enumerate the symmetries of the system, and postpone their proofs to appendix~\ref{symapp}.\par
\bigskip
\point{\bf Global $U(1)$.} For $\theta\in\mathbb{R}/(2\pi\mathbb Z)$, the map
\begin{equation}\left\{\begin{array}l
\hat\xi_{\mathbf k}^\pm\longmapsto e^{\pm i\theta}\hat\xi_{\mathbf k}^\pm\\[0.2cm]
\hat\phi_{\mathbf k}^\pm\longmapsto e^{\pm i\theta}\hat\phi_{\mathbf k}^\pm
\end{array}\right.\label{aglobaluo}\end{equation}
is a symmetry.\par
\bigskip
\point{\bf $2\pi/3$ rotation.} Let $T\mathbf k:=(k_0,e^{-i\frac{2\pi}3\sigma_2}k)$, $l_2:=(3/2,-\sqrt3/2)$ and $\mathcal T_{\mathbf k}:=e^{-i(l_2\cdot k)\sigma_3}$, the mapping
\begin{equation}\left\{\begin{array}l
\hat\xi_{\mathbf k}^\pm\longmapsto\hat\xi_{T\mathbf k}^\pm\\[0.2cm]
\hat\phi_{\mathbf k}^-\longmapsto\mathcal T_{T\mathbf k}\hat\phi_{T\mathbf k}^-,\ \hat\phi^+_{\mathbf k}\longmapsto\hat\phi_{T\mathbf k}^+\mathcal T_{T\mathbf k}^\dagger
\end{array}\right.\label{arotation}\end{equation}\par
is a symmetry.\par
\bigskip
\point{\bf Complex conjugation.} The map in which
\begin{equation}\left\{\begin{array}l
\hat\xi_{\mathbf k}^\pm\longmapsto\hat\xi_{-\mathbf k}^\pm\\[0.2cm]
\hat\phi_{\mathbf k}^\pm\longmapsto\hat\phi_{-\mathbf k}^\pm.
\end{array}\right.\label{acomplexconjugation}\end{equation}
and every complex coefficient of $h_0$ and $\mathcal V$ is mapped to its complex conjugate is a symmetry.
\bigskip
\point{\bf Vertical reflection.} Let $R_v\mathbf k=(k_0,k_x,-k_y)$,
\begin{equation}
\left\{\begin{array}l
\hat\xi_{\mathbf k}^\pm\longmapsto \hat\xi_{R_v\mathbf k}^\pm\\[0.2cm]
\hat\phi_{\mathbf k}^\pm\longmapsto \hat\phi_{R_v\mathbf k}^\pm
\end{array}\right.\label{averticalreflection}\end{equation}
is a symmetry.
\bigskip
\point{\bf Horizontal reflection.} Let $R_h\mathbf k=(k_0,-k_x,k_y)$,
\begin{equation}\left\{\begin{array}l
\hat\xi_{\mathbf k}^-\longmapsto \sigma_1\hat\xi_{R_h\mathbf k}^-,\ \hat\xi_{\mathbf k}^+\longmapsto\hat\xi_{R_h\mathbf k}^+\sigma_1\\[0.2cm]
\hat\phi_{\mathbf k}^-\longmapsto \sigma_1\hat\phi_{R_h\mathbf k}^-,\ \hat\phi_{\mathbf k}^+\longmapsto\hat\phi_{R_h\mathbf k}^+\sigma_1
\end{array}\right.\label{ahorizontalreflection}\end{equation}
is a symmetry.\par
\bigskip
\point{\bf Parity.} Let $P\mathbf k=(k_0,-k_x,-k_y)$,
\begin{equation}\left\{\begin{array}l
\hat\xi_{\mathbf k}^\pm\longmapsto i(\hat\xi_{P\mathbf k}^{\mp})^T\\[0.2cm]
\hat\phi_{\mathbf k}^\pm\longmapsto i(\hat\phi_{P\mathbf k}^{\mp})^T
\end{array}\right.\label{aparticlehole}\end{equation}
is a symmetry.
\bigskip
\point{\bf Time inversion.} Let $I\mathbf k=(-k_0,k_x,k_y)$, the mapping
\begin{equation}\left\{\begin{array}l
\hat\xi_{\mathbf k}^-\longmapsto -\sigma_3\hat\xi_{I\mathbf k}^-,\ \hat\xi_{\mathbf k}^+\longmapsto\hat\xi_{I\mathbf k}^+\sigma_3\\[0.2cm]
\hat\phi_{\mathbf k}^-\longmapsto-\sigma_3\hat\phi_{I\mathbf k}^-,\ \hat\phi_{\mathbf k}^+\longmapsto\hat\phi_{I\mathbf k}^+\sigma_3
\end{array}\right.\label{ainversiont}\end{equation}
is a symmetry.\par
\bigskip
\section{Free propagator}
\label{proppropsec}
\indent In section~\ref{tracespathintssec}, we showed how to express the free energy and the two-point Schwinger function as a formal series of truncated expectations~(\ref{grasstruncexp}). Controlling the convergence of this series is made difficult by the fact that the propagator $\hat g_{\leqslant M}$ is singular, and will require a finer analysis. In this section, we discuss which are the singularities of $\hat g_{\leqslant M}$ and how it behaves close to them, and identify three regimes in which the propagator behaves differently.\par
\subsection{Fermi points}
\indent The free propagator is singular if $k_0=0$ and $k$ is such that $H_0(k)$ is not invertible. The set of such $k$'s is called the {\it Fermi surface}.
In this subsection, we study the properties of this set. We recall the definition of $H_0$ in (\ref{hmat}),
$$ H_0( k):=
-\left(\begin{array}{*{4}{c}}
0&\gamma_1&0&\Omega^*(k)\\
\gamma_1&0&\Omega(k)&0\\
0&\Omega^*(k)&0&\gamma_3\Omega(k)e^{3ik_x}\\
\Omega(k)&0&\gamma_3\Omega^*(k)e^{-3ik_x}&0
\end{array}\right)$$
so that, using corollary~\ref{Anokzprop} (see appendix~\ref{inversapp}),
\begin{equation}
\det H_0(k)=\left|\Omega^2(k)-\gamma_1\gamma_3\Omega^*(k)e^{-3ik_x}\right|^2.
\label{detHz}\end{equation}
It is then straightforward to compute the solutions of $\det H_0(k)=0$ (see appendix~\ref{fermiapp} for details): we find that as long as $0<\gamma_1\gamma_3<2$, there are 8 Fermi points:
\begin{equation}\left\{\begin{array}l
p_{F,0}^\omega:=\left(\frac{2\pi}{3},\omega\frac{2\pi}{3\sqrt3}\right)\\[0.5cm]
p_{F,1}^\omega:=\left(\frac{2\pi}{3},\omega\frac{2}{\sqrt3}\arccos\left(\frac{1-\gamma_1\gamma_3}{2}\right)\right)\\[0.5cm]
p_{F,2}^\omega:=\left(\frac{2\pi}{3}+\frac{2}{3}\arccos\left(\frac{\sqrt{1+\gamma_1\gamma_3}(2-\gamma_1\gamma_3)}{2}\right),\omega\frac{2}{\sqrt3}\arccos\left(\frac{1+\gamma_1\gamma_3}{2}\right)\right)\\[0.5cm]
p_{F,3}^\omega:=\left(\frac{2\pi}{3}-\frac{2}{3}\arccos\left(\frac{\sqrt{1+\gamma_1\gamma_3}(2-\gamma_1\gamma_3)}{2}\right),\omega\frac{2}{\sqrt3}\arccos\left(\frac{1+\gamma_1\gamma_3}{2}\right)\right).
\end{array}\right.\label{fermdef}\end{equation}
for $\omega\in\{-,+\}$. Note that
\begin{equation}\begin{array}c
p_{F,1}^\omega=p_{F,0}^\omega+\left(0,\omega\frac{2}{3}\gamma_1\gamma_3\right)+O(\epsilon^4),\quad
p_{F,2}^\omega=p_{F,0}^\omega+\left(\frac{1}{\sqrt3}\gamma_1\gamma_3,-\omega\frac{1}{3}\gamma_1\gamma_3\right)+O(\epsilon^4),\\[0.5cm]
p_{F,3}^\omega=p_{F,0}^\omega+\left(-\frac{1}{\sqrt3}\gamma_1\gamma_3,-\omega\frac{1}{3}\gamma_1\gamma_3\right)+O(\epsilon^4).
\end{array}\label{eq:3.2bis}\end{equation}
The points $p_{F,j}^\omega$ for $j=1,2,3$ are labeled as per figure~\ref{figferm}.\par
\subsection{Behavior around the Fermi points}
\indent In this section, we compute the dominating behavior of $\hat g(\mathbf k)$ close to its singularities, that is close to $\mathbf p_{F,j}^\omega:=(0,p_{F,j}^\omega)$. We recall that $\hat A(\mathbf k):=(-ik_0\mathds1+H_0(k))$ and $\hat g(\mathbf k)=\chi_0(2^{-M}|k_0|)\hat A^{-1}(\mathbf k)$.
\bigskip
\point{\bf First regime.} We define $k':=k-p_{F,0}^\omega=(k'_x,k'_y)$, $\mathbf k':=(k_0,k')$. We have
\begin{equation}\Omega(p_{F,0}^\omega+k')=\frac{3}{2}(ik'_x+\omega k'_y)+O(|k'|^2)=:\xi+O(|k'|^2)\label{eq:omega}\end{equation}
so that, by using (\ref{detAexpr}) with $(\mathfrak a,\mathfrak b,\mathfrak c,\mathfrak x,\mathfrak z)=-(\gamma_1,\Omega(k),\gamma_3\Omega(k)e^{3ik_x}, k_0,k_0)$,
\begin{equation}
\det\hat A(\mathbf p_{F,0}^\omega+\mathbf k')=(k_0^2+|\xi|^2)^2+O(\|\mathbf k'\|_{\mathrm I}^5,\epsilon^2\|\mathbf k'\|_{\mathrm I}^2)
\label{detao}\end{equation}
where
\begin{equation}
\|\mathbf k'\|_{\mathrm I}:=\sqrt{k_0^2+|\xi|^2}
\label{normo}\end{equation}
in which the label $_{\mathrm I}$ stands for ``first regime''. If
\begin{equation}
\kappa_1\epsilon<\|\mathbf k'\|_{\mathrm I}<\bar\kappa_0
\label{odef}\end{equation}
for suitable constants $\kappa_1, \bar\kappa_0>0$, then the remainder term in (\ref{detao}) is smaller than the explicit term, so that~(\ref{detao}) is adequate in this regime, which we call the ``first regime''.
\indent We now compute the dominating part of $\hat A^{-1}$ in this regime. The computation is carried out in the following way: we neglect terms of order $\gamma_1$, $\gamma_3$ and $|k'|^2$ in $\hat A$, invert the resulting matrix using~(\ref{invAexpr}), prove that this inverse is bounded by $(\mathrm{const}.)\ \|\mathbf k'\|_{\mathrm{I}}^{-1}$, and deduce a bound on the error terms. We thus find
\begin{equation}
\hat A^{-1}(\mathbf p_{F,0}^\omega+\mathbf k')=
-\frac{1}{k_0^2+|\xi|^2}
\left(\begin{array}{*{4}{c}}
-ik_0&0&0&\xi^*\\
0&-ik_0&\xi&0\\
0&\xi^*&-ik_0&0\\
\xi&0&0&-ik_0
\end{array}\right)
\left(\mathds1+O(\|\mathbf k'\|_{\mathrm{I}},\epsilon\|\mathbf k'\|_{\mathrm{I}}^{-1})\right)
\label{freepropzo}\end{equation}
and
\begin{equation}
|\hat A^{-1}(\mathbf p_{F,0}^\omega+\mathbf k')|\leqslant(\mathrm{const}.)\ \|\mathbf k'\|_{\mathrm{I}}^{-1}.
\label{boundfreepropzo}\end{equation}
Note that, recalling that the basis in which we wrote $A^{-1}$ is $\{a,\tilde b,\tilde a,b\}$, each graphene layer is decoupled from the other in the dominating part of~(\ref{freepropzo}).\par
\bigskip
\point{\bf Ultraviolet regime.} The regime in which $\|\mathbf k'\|_{\mathrm{I}}\geqslant\bar\kappa_0$ for both $\omega=\pm$, and is called the {\it ultraviolet} regime.
For such $\mathbf k'=:\mathbf k-\mathbf p_{F,0}^\omega$, one easily checks that
\begin{equation}
|\hat A^{-1}(\mathbf k)|\leqslant(\mathrm{const}.)\ |\mathbf k|^{-1}.
\label{boundfreepropuv}\end{equation}
\bigskip
\point{\bf Second regime.} We now go beyond the first regime: we assume that $\|\mathbf k'\|_{\mathrm I}\leqslant \kappa_1\epsilon$ and, using again (\ref{eq:omega}) and (\ref{detAexpr}), we write
\begin{equation}
\det\hat A(\mathbf p_{F,0}^\omega+\mathbf k')=\gamma_1^2k_0^2+|\xi|^4+O(\epsilon^{7/2}\|\mathbf k'\|_{\mathrm{II}}^{3/2},\epsilon^5\|\mathbf k'\|_{\mathrm{II}},\epsilon\|\mathbf k'\|_{\mathrm{II}}^3)\label{detat}\end{equation}
where
\begin{equation}
\|\mathbf k'\|_{\mathrm{II}}:=\sqrt{k_0^2+\frac{|\xi|^4}{\gamma_1^2}}.
\label{normt}\end{equation}
If
\begin{equation}
\kappa_2\ \epsilon^3<\|\mathbf k'\|_{\mathrm{II}}<\bar\kappa_1\ \epsilon
\label{tdef}\end{equation}
for suitable constants $\kappa_2,\bar\kappa_1>0$, then the remainder in~(\ref{detat}) is smaller than the explicit term, and we thus define the ``second regime'', for which~(\ref{detat}) is appropriate.\par
\bigskip
\indent We now compute the dominating part of $\hat A^{-1}$ in this regime. To that end, we define the dominating part $\mathfrak L_{\mathrm{II}}\hat A$ of $\hat A$ by neglecting the terms of order $\gamma_3$ and $|k'|^2$ in $\hat A$ as well as the elements $\hat A_{aa}$ and $\hat A_{\tilde b\tilde b}$ (which are both equal to $-ik_0$), block-diagonalize it using proposition~\ref{blockdiagprop} (see appendix~\ref{diagapp}) and invert it:
\begin{equation}
\left(\mathfrak L_{\mathrm{II}}\hat A(\mathbf k)\right)^{-1}=\left(\begin{array}{*{2}{c}}\mathds1&M_{\mathrm{II}}(\mathbf k)^\dagger\\0&\mathds1\end{array}\right)\left(\begin{array}{*{2}{c}}a_{\mathrm{II}}^{(M)}&0\\0&a_{\mathrm{II}}^{(m)}(\mathbf k)\end{array}\right)\left(\begin{array}{*{2}{c}}\mathds1&0\\M_{\mathrm{II}}(\mathbf k)&\mathds1\end{array}\right)
\label{LfpAt}\end{equation}
where
\begin{equation}
a_{\mathrm{II}}^{(M)}:=-\left(\begin{array}{*{2}{c}}0&\gamma_1^{-1}\\\gamma_1^{-1}&0\end{array}\right),\quad
a_{\mathrm{II}}^{(m)}(\mathbf p_{F,0}^\omega+\mathbf k'):=\frac{\gamma_1}{\gamma_1^2k_0^2+|\xi|^4}\left(\begin{array}{*{2}{c}}
i\gamma_1k_0&(\xi^*)^2\\
\xi^2&i\gamma_1k_0\end{array}\right)
\label{compsLfpAt}\end{equation}
and
\begin{equation}
M_{\mathrm{II}}(\mathbf p_{F,0}^\omega+\mathbf k'):=-\frac1{\gamma_1}\left(\begin{array}{*{2}{c}}\xi^*&0\\0&\xi\end{array}\right).
\label{diagLfpAt}\end{equation}
We then bound the right side of~(\ref{LfpAt}), and find
\begin{equation}|\left(\mathfrak L_{\mathrm{II}}\hat A(\mathbf p_{F,0}^\omega+\mathbf k')\right)^{-1}|\leqslant (\mathrm{const}.)\ \left(\begin{array}{*{2}{c}}
\epsilon^{-1}&\epsilon^{-1/2}\|\mathbf k'\|_{\mathrm{II}}^{-1/2}\\
\epsilon^{-1/2}\|\mathbf k'\|_{\mathrm{II}}^{-1/2}&\|\mathbf k'\|_{\mathrm{II}}^{-1}\end{array}\right),\label{LA-1}\end{equation}
in which the bound should be understood as follows: the upper-left element in~(\ref{LA-1}) is the bound on the upper-left $2\times2$ block of $(\mathfrak L_{\mathrm{II}}\hat A)^{-1}$, and similarly for the upper-right, lower-left and lower-right. Using this bound in $$\hat A^{-1}(\mathbf k)=
\left(\mathfrak L_{\mathrm{II}}\hat A(\mathbf k)\right)^{-1}\Big[\mathds{1}+\big(\hat A(\mathbf k)-\mathfrak L_{\mathrm{II}}\hat A(\mathbf k)\big)\left(\mathfrak L_{\mathrm{II}}\hat A(\mathbf k)\right)^{-1}\Big]^{-1}$$
we deduce a bound on the error term in square brackets and find
\begin{equation}\begin{largearray}
\hat A^{-1}(\mathbf p_{F,0}^\omega+\mathbf k')=\left(\begin{array}{*{2}{c}}\mathds1&M_{\mathrm{II}}(\mathbf k)^\dagger\\0&\mathds1\end{array}\right)\left(\begin{array}{*{2}{c}}a_{\mathrm{II}}^{(M)}&0\\0&a_{\mathrm{II}}^{(m)}(\mathbf k)\end{array}\right)\left(\begin{array}{*{2}{c}}\mathds1&0\\M_{\mathrm{II}}(\mathbf k)&\mathds1\end{array}\right)\cdot\\[0.5cm]
\hfill\cdot(\mathds1+O(\epsilon^{-1/2}\|\mathbf k'\|_{\mathrm{II}}^{1/2},\epsilon^{3/2}\|\mathbf k'\|_{\mathrm{II}}^{-1/2}))
\end{largearray}\label{freepropzt}\end{equation}
which implies the analogue of (\ref{LA-1}) for $\hat A^{-1}$,
\begin{equation}
|\hat A^{-1}(\mathbf p_{F,0}^\omega+\mathbf k')|\leqslant(\mathrm{const}.)\ \left(\begin{array}{*{2}{c}}
\epsilon^{-1}&\epsilon^{-1/2}\|\mathbf k'\|_{\mathrm{II}}^{-1/2}\\
\epsilon^{-1/2}\|\mathbf k'\|_{\mathrm{II}}^{-1/2}&\|\mathbf k'\|_{\mathrm{II}}^{-1}\end{array}\right) .
\label{boundfreepropzt}\end{equation}
\bigskip
{\bf Remark:} Using the explicit expression for $\hat A^{-1}(\mathbf p_{F,0}^\omega+\mathbf k')$ obtained by applying proposition~\ref{matinvprop} (see appendix~\ref{inversapp}), one can show that the error term on the right side
of (\ref{freepropzt})
can be improved to $O(\epsilon^{-1}\|\mathbf k'\|_{\mathrm{II}},\epsilon^{3/2}\|\mathbf k'\|_{\mathrm{II}}^{-1/2}))$. Since we will not need this improved bound in the following,
we do not belabor further details.
\par
\bigskip
\point{\bf Intermediate regime.} In order to derive~(\ref{freepropzt}), we assumed that $\|\mathbf k'\|_{\mathrm{II}}<\bar\kappa_1\epsilon$ with $\bar\kappa_1$ small enough. In the intermediate regime defined by $\bar\kappa_1\epsilon<\|\mathbf k'\|_{\mathrm{II}}$ and $\|\mathbf k'\|_{\mathrm I}<\kappa_1\epsilon$, we have that $\|\mathbf k'\|_{\mathrm I}\sim \|\mathbf k'\|_{\mathrm{II}}\sim \epsilon$ (given two positive functions $a(\epsilon)$ and $b(\epsilon)$, the symbol $a\sim b$ stands for $c b\leqslant a\leqslant C b$ for some universal constants $C>c>0$). Moreover,
\begin{equation}
\det\hat A(\mathbf p_{F,0}^\omega+\mathbf k')=(k_0^2+|\xi|^2)^2+\gamma_1^2k_0^2+O(\epsilon^5)
\label{detaot}\end{equation}
therefore $|\det\hat A|>(\mathrm{const}.)\ \epsilon^4$ and
\begin{equation}
|\hat A^{-1}(\mathbf p_{F,0}^\omega,\mathbf k')|\leqslant(\mathrm{const}.)\ \epsilon^{-1}
\label{boundfreepropzot}\end{equation}
which is identical to the bound at the end of the first regime and at the beginning of the second.\par
\bigskip
\point{\bf Third regime.} We now probe deeper, beyond the second regime, and assume that $\|\mathbf k'\|_{\mathrm{II}}\leqslant \kappa_2\epsilon^3$.
Since we will now investigate the regime in which $|k'|<(\mathrm{const}.)\ \epsilon^2$, we will need to consider all the Fermi points $p_{F,j}^\omega$ with $j\in\{0,1,2,3\}$.\par
\bigskip
\subpoint{\bf Around $\mathbf p_{F,0}^\omega$} We start with the neighborhood of $\mathbf p_{F,0}^\omega$:
\begin{equation}
\det\hat A(\mathbf p_{F,0}^\omega+\mathbf k')=\gamma_1^2(k_0^2+\gamma_3^2|\xi|^2)+O(\epsilon^{-1}\|\mathbf k'\|_{\mathrm{III}}^3)
\label{detath}\end{equation}
where
\begin{equation}
\|\mathbf k'\|_{\mathrm{III}}:=\sqrt{k_0^2+\gamma_3^2|\xi|^2}.
\label{normth}\end{equation}
The third regime around $\mathbf p_{F,0}^\omega$ is defined by
\begin{equation}
\|\mathbf k'\|_{\mathrm{III}}<\bar\kappa_2\epsilon^3
\label{thdef}\end{equation}
for some $\bar\kappa_2<\kappa_2$.
The computation of the dominating part of $\hat A^{-1}$ in this regime around $\mathbf p_{F,0}^\omega$ is similar to that in the second regime, but for the fact that we only neglect the terms of order $|k'|^2$ in $\hat A$ as well as the elements $\hat A_{aa}$ and $\hat A_{\tilde b\tilde b}$. In addition, the terms that are of order $\epsilon^{-3}\|\mathbf k'\|_{\mathrm{III}}^2$ that come out of the computation of the dominating part of $\hat A$ in block-diagonal form are also put into the error term. We thus find
\begin{equation}
\hat A^{-1}(\mathbf k)=\left(\begin{array}{*{2}{c}}\mathds1&M_{\mathrm{III},0}(\mathbf k)^\dagger\\0&\mathds1\end{array}\right)\left(\begin{array}{*{2}{c}}a_{\mathrm{III},0}^{(M)}&0\\0&a_{\mathrm{III},0}^{(m)}(\mathbf k)\end{array}\right)\left(\begin{array}{*{2}{c}}\mathds1&0\\M_{\mathrm{III},0}(\mathbf k)&\mathds1\end{array}\right)
(\mathds1+O(\epsilon^{-3}\|\mathbf k'\|_{\mathrm{III}}))
\label{freepropzth}\end{equation}
where
\begin{equation}
a_{\mathrm{III},0}^{(M)}:=-\left(\begin{array}{*{2}{c}}0&\gamma_1^{-1}\\\gamma_1^{-1}&0\end{array}\right),\quad
a_{\mathrm{III},0}^{(m)}(\mathbf p_{F,0}^\omega+\mathbf k'):=-\frac{1}{k_0^2+\gamma_3^2|\xi|^2}\left(\begin{array}{*{2}{c}}
-ik_0&\gamma_3\xi\\
\gamma_3\xi^*&-ik_0\end{array}\right)
\label{compsLfpAth}\end{equation}
and
\begin{equation}
M_{\mathrm{III},0}(\mathbf p_{F,0}^\omega+\mathbf k'):=-\frac1{\gamma_1}\left(\begin{array}{*{2}{c}}\xi^*&0\\0&\xi\end{array}\right)
\label{diagLfpAth}\end{equation}
and
\begin{equation}
|\hat A^{-1}(\mathbf p_{F,0}^\omega+\mathbf k')|\leqslant(\mathrm{const}.)\ \left(\begin{array}{*{2}{c}}
\epsilon^{-1}&\epsilon^{-2}\\
\epsilon^{-2}&\|\mathbf k'\|_{\mathrm{III}}^{-1}\end{array}\right).
\label{boundfreepropzth}\end{equation}
\bigskip
\subpoint{\bf Around $\mathbf p_{F,1}^\omega$} We now discuss the neighborhood of $\mathbf p_{F,1}^\omega$. We define $k'_{1}:=k-p_{F,1}^\omega=(k'_{1,x},k'_{1,y})$ and $\mathbf k'_{1}:=(k_0,k'_{1})$. We have
\begin{equation}\Omega(p_{F,1}^\omega+k'_{1})=\gamma_1\gamma_3+\xi_1+O(\epsilon^2|k'_{1}|)\label{eq:O1}\end{equation}
where
$$\xi_1:=\frac{3}{2}(ik'_{1,x}+\omega k'_{1,y}).$$
Using (\ref{detAexpr}) and (\ref{detnokz}), we obtain
\begin{equation}
\det\hat A(\mathbf p_{F,1}^\omega+\mathbf k'_{1})=\gamma_1^2 k_0^2+|\Omega^2-\gamma_1\gamma_3\Omega^*e^{-3ik'_{1,x}}|^2+O(\epsilon^4 |k_0|^2)
\end{equation}
where $\Omega$ is evaluated at $p_{F,1}^\omega+k'_{1}$. Injecting (\ref{eq:O1}) into this equation, we find
\begin{equation}
\det\hat A(\mathbf p_{F,1}^\omega+\mathbf k'_{1})=\gamma_1^2 (k_0^2+\gamma_3^2|x_1|^2)
+O(\epsilon^4\|\mathbf k'_{1}\|_{\mathrm{III}}^2,\epsilon^{-1}\|\mathbf k'_{1}\|_{\mathrm{III}}^3)
\end{equation}
where
$$x_1:=\frac{3}{2}(3ik'_{1,x}+\omega k'_{1,y}).$$
The third regime around $p_{F,1}^\omega$ is therefore defined by
$$\|\mathbf k'_{1}\|_{\mathrm{III}}<\bar\kappa_{2}\ \epsilon^3$$
where $\bar\kappa_2$ can be assumed to be the same as in (\ref{thdef}) without loss of generality.
The dominating part of $\hat A^{-1}$ in this regime around $\mathbf p_{F,1}^\omega$ is similar to that around $\mathbf p_{F,0}^\omega$, except that we neglect the terms of order $\epsilon^2k'_{1}$ in $\hat A$ as well as the elements $\hat A_{aa}$ and $\hat A_{\tilde b\tilde b}$. As around $\mathbf p_{F,0}^\omega$, the terms of order $\epsilon^{-3}\|\mathbf k'_{1}\|_{\mathrm{III}}^2$ are put into the error term. We thus find
\begin{equation}\begin{largearray}
\hat A^{-1}(\mathbf k)=\left(\begin{array}{*{2}{c}}\mathds1&M_{\mathrm{III},1}(\mathbf k)^\dagger\\0&\mathds1\end{array}\right)\left(\begin{array}{*{2}{c}}a_{\mathrm{III},1}^{(M)}&0\\0&a_{\mathrm{III},1}^{(m)}(\mathbf k)\end{array}\right)\left(\begin{array}{*{2}{c}}\mathds1&0\\M_{\mathrm{III},1}(\mathbf k)&\mathds1\end{array}\right)\cdot\\[0.5cm]
\hfill\cdot(\mathds1+O(\epsilon,\epsilon^{-3}\|\mathbf k'\|_{\mathrm{III}}))
\end{largearray}\label{freepropoth}\end{equation}
where
\begin{equation}
a_{\mathrm{III},1}^{(M)}:=-\left(\begin{array}{*{2}{c}}0&\gamma_1^{-1}\\\gamma_1^{-1}&0\end{array}\right),\quad
a_{\mathrm{III},1}^{(m)}(\mathbf p_{F,1}^\omega+\mathbf k'):=\frac{1}{k_0^2+\gamma_3^2|x_1|^2}\left(\begin{array}{*{2}{c}}
ik_0&\gamma_3 x_1^*\\
\gamma_3x_1&ik_0\end{array}\right)
\label{compsLfpAthj}\end{equation}
and
\begin{equation}
M_{\mathrm{III},1}(\mathbf p_{F,1}^\omega+\mathbf k'_{1}):=-\gamma_3\mathds1-\frac1{\gamma_1}\left(\begin{array}{*{2}{c}}\xi_1^*&0\\0&\xi_1\end{array}\right)
\label{diagLfpAthj}\end{equation}
and
\begin{equation}
|\hat A^{-1}(\mathbf p_{F,1}^\omega+\mathbf k'_{1})|\leqslant(\mathrm{const}.)\ \left(\begin{array}{*{2}{c}}
\epsilon^{2}\|\mathbf k'_{1}\|_{\mathrm{III}}^{-1}&\epsilon\|\mathbf k'_{1}\|_{\mathrm{III}}^{-1}\\
\epsilon\|\mathbf k'_{1}\|_{\mathrm{III}}^{-1}&\|\mathbf k'_{1}\|_{\mathrm{III}}^{-1}\end{array}\right).
\label{boundfreepropzthj}\end{equation}
\bigskip
\subpoint{\bf Around $\mathbf p_{F,j}^\omega$} The behavior of $\hat g(\mathbf k)$ around $p_{F,j}^\omega$ for $j\in\{2,3\}$ can be deduced from~(\ref{freepropoth}) by using the symmetry (\ref{arotation})
under $2\pi/3$ rotations: if we define $k'_{j}:=k-p_{F,j}^\omega=(k'_{j,x},k'_{j,y})$, $\mathbf k'_{j}:=(k_0,k'_{j})$ then, for $j=2,3$ and $\omega\pm$,
\begin{equation}
\hat A^{-1}(\mathbf k'_{j}+\mathbf p_{F,j}^\omega)=
\left(\begin{array}{*{2}{c}}\mathds1&0\\0&\mathcal T_{T\mathbf k'_{j}+\mathbf p_{F,j-\omega}^\omega}\end{array}\right)\hat A^{-1}(T\mathbf k'_{j}+\mathbf p_{F,j-\omega}^{\omega})\left(\begin{array}{*{2}{c}}\mathds1&0\\0&\mathcal T^\dagger_{T\mathbf k'_{j}+\mathbf p_{F,j-\omega}^\omega}\end{array}\right)
\label{freeproptth}\end{equation}
where $T$ and $\mathcal T_{\mathbf k}$ were defined above~(\ref{arotation}), and $\mathbf p_{F,4}^{-}\equiv \mathbf p_{F,1}^{-}$.
In addition, if $\mathbf k'_{2}$ and $\mathbf k'_{3}$ are in the third regime, then
$\mathcal T_{T\mathbf k'_{j}+\mathbf p_{F,j}^\omega}=e^{-i\omega\frac{2\pi}3\sigma_3}+O(\epsilon^2)$.
\bigskip
\point{\bf Intermediate regime.} We are left with an intermediate regime between the second and third regimes, defined by
\begin{equation}
\bar\kappa_2\epsilon^3<\|\mathbf k'\|_{\mathrm{III}}\ ,\quad\|\mathbf k'\|_{\mathrm{II}}<\kappa_2\epsilon^3\quad\mathrm{and}\quad\bar\kappa_{2}\epsilon^3<\|\mathbf k'_{j}\|_{\mathrm{III}},\
\forall j\in\{1,2,3\},
\label{inttth}\end{equation}
which implies
$$\|\mathbf k'\|_{\mathrm{III}}\sim \|\mathbf k'\|_{\mathrm{II}}\sim \|\mathbf k'_{j}\|_{\mathrm{III}}\sim \epsilon^3$$
and
\begin{equation}
\det\hat A(\mathbf p_{F,0}^\omega+\mathbf k')=\gamma_1^2k_0^2+\left|\xi^2-\gamma_1\gamma_3\xi^*\right|^2+O(\epsilon^{10}).
\label{detatth}\end{equation}
One can prove (see appendix~\ref{boundproptthapp}) that injecting~(\ref{inttth}) into~(\ref{detatth}) implies that $|\det\hat A|\geqslant(\mathrm{const}.)\ \epsilon^8$, which in turn implies that
\begin{equation}
|\hat A^{-1}(\mathbf p_{F,0}^\omega+\mathbf k')|\leqslant (\mathrm{const}.)\ \left(\begin{array}{*{2}{c}}
\epsilon^{-1}&\epsilon^{-2}\\
\epsilon^{-2}&\epsilon^{-3}\end{array}\right)
\label{boundfreepropztth}\end{equation}
which is identical to the bound at the end of the second regime and at the beginning of the third.\par\bigskip
\point{\bf Summary.} Let us briefly summarize this sub-section: we defined the norms
\begin{equation}
\|\mathbf k'\|_{\mathrm I}:=\sqrt{k_0^2+|\xi|^2},\quad
\|\mathbf k'\|_{\mathrm{II}}:=\sqrt{k_0^2+\frac{|\xi^4|}{\gamma_1^2}},\quad \|\mathbf k'\|_{\mathrm{III}}:=\sqrt{k_0^2+\gamma_3^2|\xi|^2},
\label{normsdef}\end{equation}
and identified an {\it ultraviolet} regime and three {\it infrared} regimes in which the free propagator $\hat g(\mathbf k)$ behaves differently:
\begin{itemize}
\item for $\| \mathbf k'\|_{\mathrm{I}}>\bar\kappa_0$, (\ref{boundfreepropuv}) holds.
\item for $\kappa_1\epsilon<\|\mathbf k'\|_{\mathrm I}<\bar\kappa_0$, (\ref{freepropzo}) holds.
\item for $\kappa_2\epsilon^3<\|\mathbf k'\|_{\mathrm{II}}<\bar\kappa_1\epsilon$, (\ref{freepropzt}) holds.
\item for $\|\mathbf k'\|_{\mathrm{III}}<\bar\kappa_2\epsilon^3$, (\ref{freepropzth}) holds, for $\|\mathbf k'_{1}\|_{\mathrm{III}}<\bar\kappa_{2}\epsilon^3$, (\ref{freepropoth}) holds,
and similarly for the $j=2,3$ cases.
\end{itemize}
\section{Multiscale integration scheme}
\label{schemesec}
\indent In this section, we describe the scheme that will be followed in order to compute the right side of~(\ref{grasstruncexp}). We will first define a {\it multiscale decomposition} in each regime which will play an essential role in showing that the formal series in~(\ref{grasstruncexp}) converges. In doing so, we will define {\it effective} interactions and propagators, which will be defined in $\mathbf k$-space, but since we wish to use the determinant expansion in lemma~\ref{detexplemma}\ to compute and bound the effective truncated expectations, we will have to define the effective quantities in $\mathbf x$-space as well. Once this is done, we will write bounds for the propagator in terms of scales.\par
\subsection{Multiscale decomposition}
\label{multiscalesec}
\indent We will now discuss the scheme we will follow to compute the Gaussian Grassmann integrals in terms of which the free energy and two-point Schwinger function were expressed in~(\ref{freeengrass}) and~(\ref{schwinform}). The main idea is to decompose them into scales, and compute them one scale at a time. The result of the integration over one scale will then be considered as an {\it effective theory} for the remaining ones.\par
\bigskip
\indent Throughout this section, we will use a smooth cutoff function $\chi_0(\rho)$, which returns 1 for $\rho\leqslant1/3$ and 0 for $\rho\geqslant2/3$.\par
\point{\bf Ultraviolet regime.} Let $\bar{\mathfrak h}_0:=\lfloor\log_2(\bar\kappa_0)\rfloor$ (in which $\bar\kappa_0$ is the constant that appeared after~(\ref{normsdef}) which defines the inferior bound of the ultraviolet regime). For $h\in\{\bar{\mathfrak h}_0,\cdots,M\}$ and $h'\in\{\bar{\mathfrak h}_0+1,\cdots,M\}$, we define
\begin{equation}\begin{array}c
f_{\leqslant h'}(\mathbf k):=\chi_0(2^{-h'}|k_0|),\quad
f_{\leqslant\bar{\mathfrak h}_0}(\mathbf k):=\sum_{\omega=\pm}\chi_0(2^{-\bar{\mathfrak h}_0}\|\mathbf k-\mathbf p_{F,0}^\omega\|_{\mathrm{I}}),\\[0.2cm]
f_{h'}(\mathbf k):=f_{\leqslant h'}(\mathbf k)-f_{\leqslant h'-1}(\mathbf k)\\[0.2cm]
\mathcal B_{\beta,L}^{(\leqslant h)}:=\mathcal B_{\beta,L}\cap\mathrm{supp}f_{\leqslant h},\quad
\mathcal B_{\beta,L}^{(h')}:=\mathcal B_{\beta,L}\cap\mathrm{supp}f_{h'},\quad
\end{array}\label{indicuvh}\end{equation}
in which $\|\cdot\|_{\mathrm{I}}$ is the norm defined in~(\ref{normsdef}). In addition, we define
\begin{equation}
\hat g_{h'}(\mathbf k):=f_{h'}(\mathbf k)\hat A^{-1}(\mathbf k),\quad \hat g_{\leqslant h}(\mathbf k):=f_{\leqslant h}(\mathbf k)\hat A^{-1}(\mathbf k)
\label{inddresspropexpuv}\end{equation}
so that, in particular,
$$\hat g_{\leqslant M}(\mathbf k)=\hat g_{\leqslant M-1}(\mathbf k)+\hat g_M(\mathbf k).$$
Furthermore, it follows from the addition property~(\ref{addprop}) that for all $h\in\{\bar{\mathfrak h}_0,\cdots,M-1\}$,
\begin{equation}\left\{\begin{array}{>{\displaystyle}l}
\int P_{\leqslant M}(d\psi)\ e^{-\mathcal V(\psi)}=e^{-\beta|\Lambda|F_{h}}\int P_{\leqslant h}(d\psi^{(\leqslant h)})\ e^{-\mathcal V^{(h)}(\psi^{(\leqslant h)})}\\[0.5cm]
\int P_{\leqslant M}(d\psi)\ e^{-\mathcal W(\psi,\hat J_{\mathbf k,\underline\alpha})}=e^{-\beta|\Lambda|F_{h}}\int P_{\leqslant h}(d\psi^{(\leqslant h)})\ e^{-\mathcal W^{(h)}(\psi^{(\leqslant h)},\hat J_{\mathbf k,\underline\alpha})}
\end{array}\right.\label{freeengrassuv}\end{equation}
where $P_{\leqslant h}(d\psi^{(\leqslant h)})\equiv P_{\hat g_{\leqslant h}}(d\psi^{(\leqslant h)})$,
\begin{equation}\begin{array}{r@{\ }>{\displaystyle}l}
-\beta|\Lambda|F_{h}-\mathcal V^{(h)}(\psi^{(\leqslant h)}):=&-\beta|\Lambda|F_{h+1}+\log\int P_{h+1}(d\psi^{(h+1)})\ e^{-\mathcal V^{(h+1)}(\psi^{(h+1)}+\psi^{(\leqslant h)})}\\[0.5cm]
=&-\beta|\Lambda|F_{h+1}+\sum_{N=1}^\infty\frac{(-1)^N}{N!}\mathcal E_{h+1}^T(\mathcal V^{(h+1)}(\psi^{(h+1)}+\psi^{(\leqslant h)});N)
\end{array}\label{effpotuv}\end{equation}
and
\begin{equation}\begin{largearray}
-\beta|\Lambda|(F_{h}-F_{h+1})-\mathcal W^{(h)}(\psi^{(\leqslant h)},\hat J_{\mathbf k,\underline\alpha})\\[0.2cm]
\hfill:=\sum_{N=1}^\infty\frac{(-1)^N}{N!}\mathcal E_{(h+1)}^T(\mathcal W^{(h+1)}(\psi^{(h+1)}+\psi^{(\leqslant h)},\hat J_{\mathbf k,\underline\alpha});N)
\end{largearray}\label{effpotWuv}\end{equation}
in which the induction is initialized by
$$\mathcal V^{(M)}:=\mathcal V,\quad\mathcal W^{(M)}:=\mathcal W,\quad F_M:=0.$$
\bigskip
\point{\bf First regime.} We now decompose the first regime into scales. The main difference with the ultraviolet regime is that we incorporate the quadratic part of the effective potential into the propagator at each step of the multiscale integration. This is necessary to get satisfactory bounds later on. The propagator will therefore be changed, or {\it dressed}, inductively at every scale, as discussed below.\par
\bigskip
\indent Let $\mathfrak h_1:=\lceil\log_2(\kappa_1\epsilon)\rceil$ (in which $\kappa_1$ is the constant that appears after~(\ref{normsdef}) which defines the inferior bound of the first regime), and $\|\cdot\|_{\mathrm I}$ be the norm defined in~(\ref{normsdef}). We define for $h\in\{\mathfrak h_1,\cdots,\bar{\mathfrak h}_0\}$ and $h'\in\{\mathfrak h_1+1,\cdots,\bar{\mathfrak h}_0\}$,
\begin{equation}
\begin{array}c
f_{\leqslant h,\omega}(\mathbf k):=\chi_0(2^{-h}\|\mathbf k-\mathbf p_{F,0}^\omega\|_{\mathrm I}),\quad f_{h',\omega}(\mathbf k):=f_{\leqslant h',\omega}(\mathbf k)-f_{\leqslant h'-1,\omega}(\mathbf k)\\[0.2cm]
\mathcal B_{\beta,L}^{(\leqslant h,\omega)}:=\mathcal B_{\beta,L}\cap\mathrm{supp}f_{\leqslant h,\omega},\quad
\mathcal B_{\beta,L}^{(h',\omega)}:=\mathcal B_{\beta,L}\cap\mathrm{supp}f_{\leqslant h',\omega}
\end{array}\label{indicoh}\end{equation}
and
\begin{equation}
\hat g_{h',\omega}(\mathbf k):=f_{h',\omega}(\mathbf k)\hat A^{-1}(\mathbf k),\quad
\hat g_{\leqslant h,\omega}(\mathbf k):=f_{\leqslant h,\omega}(\mathbf k)\hat A^{-1}(\mathbf k).
\label{indundresspropexpo}\end{equation}
For $h\in\{\mathfrak h_1,\cdots,\bar{\mathfrak h}_0-1\}$, we define
\begin{equation}\begin{largearray}
-\beta|\Lambda|(F_h-F_{h+1})-\mathcal Q^{(h)}(\psi^{(\leqslant h)})-\bar{\mathcal V}^{(h)}(\psi^{(\leqslant h)})\\[0.2cm]
\hfill:=\sum_{N=1}^\infty\frac{(-1)^N}{N!}\bar{\mathcal E}_{h+1}^T(\bar{\mathcal V}^{(h+1)}(\psi^{(h+1)}+\psi^{(\leqslant h)});N)\\[0.5cm]
\hfil\mathcal Q^{(\bar{\mathfrak h}_0)}(\psi^{(\leqslant \bar{\mathfrak h}_0)})+\bar{\mathcal V}^{(\bar{\mathfrak h}_0)}(\psi^{(\leqslant \bar{\mathfrak h}_0)}):=\mathcal V^{(\bar{\mathfrak h}_0)}(\psi^{(\leqslant\bar{\mathfrak h}_0)})
\end{largearray}\label{effpotoh}\end{equation}
and
\begin{equation}\begin{largearray}
-\beta|\Lambda|(F_h-F_{h+1})-\mathcal Q^{(h)}(\psi^{(\leqslant h)})-\bar{\mathcal W}^{(h)}(\psi^{(\leqslant h)},\hat J_{\mathbf k,\underline\alpha})\\
\hfill:=\sum_{N=1}^\infty\frac{(-1)^N}{N!}\bar{\mathcal E}_{h+1}^T(\bar{\mathcal W}^{(h+1)}(\psi^{(h+1)}+\psi^{(\leqslant h)},\hat J_{\mathbf k,\underline\alpha});N)\\[0.5cm]
\hfil\mathcal Q^{(\bar{\mathfrak h}_0)}(\psi^{(\leqslant \bar{\mathfrak h}_0)})+\bar{\mathcal W}^{(\bar{\mathfrak h}_0)}(\psi^{(\leqslant \bar{\mathfrak h}_0)},\hat J_{\mathbf k,\underline\alpha}):=\mathcal W^{(\bar{\mathfrak h}_0)}(\psi^{(\leqslant\bar{\mathfrak h}_0)},\hat J_{\mathbf k,\underline\alpha})
\end{largearray}\label{effpotWoh}\end{equation}
in which $\mathcal Q^{(h)}$ is quadratic in the $\psi$, $\bar{\mathcal V}^{(h)}$ is at least quartic and $\bar{\mathcal W}^{(h)}$ has no terms that are both quadratic in $\psi$ and constant in $\hat J_{\mathbf k,\underline\alpha}$; and $\bar{\mathcal E}^T_{h+1}$ is the truncated expectation defined from the Gaussian measure $P_{\hat{\bar g}_{h+1,+}}(d\psi_+^{(h+1)})P_{\hat{\bar g}_{h+1,-}}(d\psi_-^{(h+1)})$; in which $\hat{\bar g}_{h+1,\omega}$ is the {\it dressed propagator} and is defined as follows. Let $\hat W_2^{(h)}(\mathbf k)$ denote the {\it kernel} of $\mathcal Q^{(h)}$ i.e.
\begin{equation}
\mathcal Q^{(h)}(\psi^{(\leqslant h)})=:\frac{1}{\beta|\Lambda|}\sum_{\omega,(\alpha,\alpha')}\sum_{\mathbf k\in\mathcal B_{\beta,L}^{(\leqslant h,\omega)}}\hat\psi^{(\leqslant h)+}_{\mathbf k,\omega,\alpha}\hat W^{(h)}_{2,(\alpha,\alpha')}(\mathbf k)\hat\psi^{(\leqslant h)-}_{\mathbf k,\omega,\alpha'}
\label{hatWtdeft}\end{equation}
(remark: the $_\omega$ index in $\hat\psi_{\mathbf k,\omega,\alpha}^\pm$ is redundant since given $\mathbf k$, it is defined as the unique $\omega$ that is such that $\mathbf k\in\mathcal B_{\beta,L}^{(\leqslant h,\omega)}$; it will however be needed when defining the $\mathbf x$-space counterpart of $\hat\psi_{\mathbf k,\omega,\alpha}^\pm$ below). We define $\hat{\bar g}_{h,\omega}$ and $\hat{\bar g}_{\leqslant h,\omega}$ by induction: $\hat{\bar g}_{\leqslant \bar{\mathfrak h}_0,\omega}({\bf k}):=(\hat g^{-1}_{\leqslant\bar{\mathfrak h}_0,\omega}(\mathbf k)+\hat W^{(\bar{\mathfrak h}_0)}_{2}(\mathbf k))^{-1}$ and, for $h\in\{\mathfrak h_1+1,\ldots,\bar{\mathfrak h}_0\}$,
\begin{equation}\left\{\begin{array}l
\hat{\bar g}_{h,\omega}(\mathbf k):=f_{h,\omega}(\mathbf k)f_{\leqslant h,\omega}^{-1}(\mathbf k)\hat{\bar g}_{\leqslant h,\omega}(\mathbf k)\\[0.2cm]
\left(\hat{\bar g}_{\leqslant h-1,\omega}(\mathbf k)\right)^{-1}:=f^{-1}_{\leqslant h-1,\omega}(\mathbf k)\left(\hat{\bar g}_{\leqslant h,\omega}(\mathbf k)\right)^{-1}+\hat W_2^{(h-1)}(\mathbf k)
\end{array}\right.\label{inddresspropo}\end{equation}
in which $f_{\leqslant h,\omega}^{-1}(\mathbf k)$ is equal to $1/f_{\leqslant h,\omega}(\mathbf k)$ if $f_{\leqslant h,\omega}(\mathbf k)\neq0$ and to $0$ if not.\par
\bigskip
\indent The dressed propagator is thus defined so that
\begin{equation}\left\{\begin{array}{>{\displaystyle}l}
\int P_M(d\psi)\ e^{-\mathcal V(\psi)}=e^{-\beta|\Lambda|F_h}\int \bar P_{\leqslant h}(\psi^{(\leqslant h)})\ e^{-\bar{\mathcal V}^{(h)}(\psi^{(\leqslant h)})}\\[0.5cm]
\int P_M(d\psi)\ e^{-\mathcal W(\psi,\hat J_{\mathbf k,\underline\alpha})}=e^{-\beta|\Lambda|F_h}\int \bar P_{\leqslant h}(\psi^{(\leqslant h)})\ e^{-\bar{\mathcal W}^{(h)}(\psi^{(\leqslant h)},\hat J_{\mathbf k,\underline\alpha})}
\end{array}\right.\label{freeengrassoh}\end{equation}
in which $\bar P_{\leqslant h}\equiv P_{\hat{\bar g}_{\leqslant h,+}}(d\psi_+^{(\leqslant h)})P_{\hat{\bar g}_{\leqslant h,-}}(d\psi_-^{(\leqslant h)})$. Equation~(\ref{inddresspropo}) can be expanded into a more explicit form: for $h'\in\{\mathfrak h_1+1,\ldots,\bar{\mathfrak h}_0\}$ and $h\in\{\mathfrak h_1,\cdots,\bar{\mathfrak h}_0\}$,
\begin{equation}
\hat{\bar g}_{h',\omega}(\mathbf k)=f_{h',\omega}(\mathbf k)\left(\hat{\bar A}_{h',\omega}(\mathbf k)\right)^{-1},\quad
\hat{\bar g}_{\leqslant h,\omega}(\mathbf k)=f_{\leqslant h,\omega}\left(\hat{\bar A}_{h,\omega}(\mathbf k)\right)^{-1}
\label{inddresspropexpo}\end{equation}
where
\begin{equation}
\hat{\bar A}_{h,\omega}(\mathbf k):=\hat A(\mathbf k)+f_{\leqslant h,\omega}(\mathbf k)\hat W_2^{(h)}(\mathbf k)+\sum_{h'=h+1}^{\bar{\mathfrak h}_0}\hat W_2^{(h')}(\mathbf k)
\label{Adefo}\end{equation}
(in which the sum should be interpreted as zero if $h=\bar{\mathfrak h}_0$).\par
\bigskip
\point{\bf Intermediate regime.} We briefly discuss the intermediate region between regimes~1 and~2. We define
\begin{equation}
f_{\mathfrak h_1,\omega}(\mathbf k):=\chi_0(2^{-\mathfrak h_1}\|\mathbf k-\mathbf p_{F,0}^\omega\|_{\mathrm I})-\chi_0(2^{-\bar{\mathfrak h}_1}\|\mathbf k-\mathbf p_{F,0}^\omega\|_{\mathrm{II}})=:f_{\leqslant\mathfrak h_1,\omega}(\mathbf k)-f_{\leqslant\bar{\mathfrak h}_1,\omega}(\mathbf k)
\label{indicot}\end{equation}
where $\bar{\mathfrak h}_1:=\lfloor\log_2(\bar\kappa_1\epsilon)\rfloor$, from which we define $\hat{\bar g}_{\mathfrak h_1,\omega}(\mathbf k)$ and $\hat{\bar g}_{\leqslant\bar{\mathfrak h}_1,\omega}(\mathbf k)$ in the same way as in~(\ref{inddresspropexpo}) with
\begin{equation}
\hat{\bar A}_{\bar{\mathfrak h}_1,\omega}(\mathbf k):=\hat A(\mathbf k)+f_{\leqslant \bar{\mathfrak h}_1,\omega}(\mathbf k)\hat W_2^{(\bar{\mathfrak h}_1)}(\mathbf k)+\sum_{h'=\mathfrak h_1}^{\bar{\mathfrak h}_0}\hat W_2^{(h')}(\mathbf k).
\label{Adefot}\end{equation}
The analogue of (\ref{freeengrassoh}) holds here as well.\par
\bigskip
\point{\bf Second regime.} We now define a multiscale decomposition for the integration in the second regime. Proceeding as we did in the first regime, we define $\mathfrak h_2:=\lceil\log_2(\kappa_2\epsilon^3)\rceil$, for $h\in\{\mathfrak h_2,\cdots,\bar{\mathfrak h}_1\}$ and $h'\in\{\mathfrak h_2+1,\cdots,\bar{\mathfrak h}_1\}$, we define
\begin{equation}\begin{array}c
f_{\leqslant h,\omega}(\mathbf k):=\chi_0(2^{-h}\|\mathbf k-\mathbf p_{F,0}^\omega\|_{\mathrm{II}}),\quad f_{h',\omega}(\mathbf k):=f_{\leqslant h',\omega}(\mathbf k)-f_{\leqslant h'-1,\omega}(\mathbf k)\\[0.2cm]
\mathcal B_{\beta,L}^{(\leqslant h,\omega)}:=\mathcal B_{\beta,L}\cap\mathrm{supp}f_{\leqslant h,\omega},\quad
\mathcal B_{\beta,L}^{(h',\omega)}:=\mathcal B_{\beta,L}\cap\mathrm{supp}f_{\leqslant h',\omega}.
\end{array}\label{indicth}\end{equation}
The analogues of (\ref{freeengrassoh}), and~(\ref{inddresspropexpo}) hold with
\begin{equation}
\hat{\bar A}_{h-1,\omega}(\mathbf k):=\hat A(\mathbf k)+ f_{\leqslant h-1,\omega}(\mathbf k)\hat W_2^{(h-1)}(\mathbf k)+\sum_{h'=h}^{\bar{\mathfrak h}_1}\hat W_2^{(h')}(\mathbf k)+\sum_{h'=\mathfrak h_1}^{\bar{\mathfrak h}_0}\hat W_2^{(h')}(\mathbf k).
\label{Adeft}\end{equation}
\bigskip
\point{\bf Intermediate regime.} The intermediate region between regimes~2 and~3 is defined in analogy with that between regimes~1 and~2: we let
\begin{equation}\begin{array}{>{\displaystyle}c}
f_{\mathfrak h_2,\omega}(\mathbf k):=\chi_0(2^{-\mathfrak h_2}\|\mathbf k-\mathbf p_{F,0}^\omega\|_{\mathrm{II}})-\sum_{j\in\{0,1,2,3\}}\chi_0(2^{-\bar{\mathfrak h}_2}\|\mathbf k-\mathbf p_{F,j}^\omega\|_{\mathrm{III}})\\[0.5cm]
f_{\leqslant\bar{\mathfrak h}_2,\omega,j}(\mathbf k):=\chi_0(2^{-\bar{\mathfrak h}_2}\|\mathbf k'_{\omega,j}\|_{\mathrm{III}})
\end{array}\label{indictth}\end{equation}
where $\bar{\mathfrak h}_2:=\lfloor\log_2(\bar\kappa_2\epsilon^3)\rfloor$ from which we define $\hat{\bar g}_{\mathfrak h_2,\omega}(\mathbf k)$ and $\hat{\bar g}_{\leqslant\bar{\mathfrak h}_2,\omega}(\mathbf k)$ in the same way as in~(\ref{inddresspropexpo}) with
\begin{equation}
\hat{\bar A}_{\bar{\mathfrak h}_2,\omega}(\mathbf k):=\hat A(\mathbf k)+f_{\leqslant \bar{\mathfrak h}_2,\omega}(\mathbf k)\hat W_2^{(\bar{\mathfrak h}_2)}(\mathbf k)+\sum_{h'=\mathfrak h_2}^{\bar{\mathfrak h}_1}\hat W_2^{(h')}(\mathbf k)+\sum_{h'=\mathfrak h_1}^{\bar{\mathfrak h}_0}\hat W_2^{(h')}(\mathbf k).
\label{Adeftth}\end{equation}
The analogue of (\ref{freeengrassoh}) holds here as well.\par
\bigskip
\point{\bf Third regime.} There is an extra subtlety in the third regime: we will see in section~\ref{thsec} that the singularities of the dressed propagator are slightly different from those of the bare (i.e. non-interacting) propagator: at scale $h$ the effective Fermi points $\mathbf p_{F,j}^{\omega}$ with $j=1,2,3$ are moved to $\tilde{\mathbf p}_{F,j}^{(\omega,h)}$, with
\begin{equation} \label{eq4.20bis}\|\tilde{\mathbf p}_{F,j}^{(\omega,h)}-\mathbf p_{F,j}^{\omega}\|_{\mathrm{III}}\leqslant (\mathrm{const}.)\ |U|\epsilon^3.\end{equation}
The central Fermi points, $j=0$, are left invariant by the interaction. For notational uniformity we set $\tilde{\mathbf p}_{F,0}^{(\omega,h)}\equiv \mathbf p_{F,0}^{\omega}$.
Keeping this in mind, we then proceed in a way reminiscent of the first and second regimes: let $\mathfrak h_\beta:=\lfloor\log_2(\pi/\beta)\rfloor$, for $h\in\{\mathfrak h_\beta,\cdots,\bar{\mathfrak h}_2\}$ and $h'\in\{\mathfrak h_\beta+1,\cdots,\bar{\mathfrak h}_2\}$, we define
\begin{equation}
\begin{array}c
f_{\leqslant h,\omega,j}(\mathbf k):=\chi_0(2^{-h}\|\mathbf k-\tilde{\mathbf p}_{F,j}^{(\omega,h+1)}\|_{\mathrm{III}}),\quad f_{h',\omega,j}(\mathbf k):=f_{\leqslant h',\omega,j}(\mathbf k)-f_{\leqslant h'-1,\omega,j}(\mathbf k)\\[0.2cm]
\mathcal B_{\beta,L}^{(\leqslant h,\omega,j)}:=\mathcal B_{\beta,L}\cap\mathrm{supp}f_{\leqslant h,\omega,j},\quad
\mathcal B_{\beta,L}^{(h',\omega,j)}:=\mathcal B_{\beta,L}\cap\mathrm{supp}f_{\leqslant h',\omega,j}
\end{array}\label{indicthh}\end{equation}
and the analogues of (\ref{freeengrassoh}), and~(\ref{inddresspropexpo}) hold with
\begin{equation}\begin{largearray}
\hat{\bar A}_{h-1,\omega,j}(\mathbf k):=\hat A(\mathbf k)+f_{\leqslant h-1,\omega,j}(\mathbf k)\hat W_2^{(h-1)}(\mathbf k)\\[0.2cm]
\hfill+\sum_{h'=h}^{\bar{\mathfrak h}_2}\hat W_2^{(h')}(\mathbf k)+\sum_{h'=\mathfrak h_2}^{\bar{\mathfrak h}_1}\hat W_2^{(h')}(\mathbf k)+\sum_{h'=\mathfrak h_1}^{\bar{\mathfrak h}_0}\hat W_2^{(h')}(\mathbf k).
\end{largearray}\label{Adefth}\end{equation}
\par
\bigskip
\point{\bf Last scale.} Recalling that $|k_0|\geqslant\pi/\beta$, the smallest possible scale is $\mathfrak h_\beta:=\lfloor\log_2(\pi/\beta)\rfloor$. The last integration is therefore that on scale $h=\mathfrak h_\beta+1$, after which, we are left with
\begin{equation}\left\{\begin{array}{>{\displaystyle}l}
\int P_{\leqslant M}(d\psi)\ e^{-\mathcal V(\psi)}=e^{-\beta|\Lambda|F_{\mathfrak h_\beta}}\\[0.5cm]
\int P_{\leqslant M}(d\psi)\ e^{-\mathcal W(\psi,\hat J_{\mathbf k,\underline\alpha})}=e^{-\beta|\Lambda|F_{\mathfrak h_\beta}}e^{-\mathcal W^{(\mathfrak h_\beta)}(\hat J_{\mathbf k,\underline\alpha})}.
\end{array}\right.\label{freeengrassend}\end{equation}
The subsequent sections are dedicated to the proof of the fact that both $F_{\mathfrak h_\beta}$ and $\mathcal W^{(\mathfrak h_\beta)}$ are analytic in $U$, uniformly in $L$, $\beta$ and $\epsilon$. We will do this by studying each regime, one at a time, performing a {\it tree expansion} in each of them in order to bound the terms of the series (see section~\ref{treeexpsec} and following).\par
\subsection{$\mathbf x$-space representation of the effective potentials}
\label{xspacesec}
\indent We will compute the truncated expectations arising in~(\ref{effpotuv}), (\ref{effpotWuv}), (\ref{effpotoh}) and~(\ref{effpotWoh}) using a determinant expansion (see lemma~\ref{detexplemma}) which, as was mentioned above, is only useful if the propagator and effective potential are expressed in $\mathbf x$-space. We will discuss their definition in this section. We restrict our attention to the effective potentials $\mathcal V^{(h)}$ since, in order to compute the two-point Schwinger function in the regimes we are interested in, we will not need to express the kernels of $\mathcal W^{(h)}$ in $\mathbf x$-space.
\par
\bigskip
\point{\bf Ultraviolet regime.} We first discuss the ultraviolet regime, which differs from the others in that the propagator does not depend on the index $\omega$. We write $\mathcal V^{(h)}$ in terms of its {\it kernels} (anti-symmetric in the exchange of their indices), defined as
\begin{equation}\begin{largearray}
\mathcal V^{(h)}(\psi^{(\leqslant h)})=:\sum_{l=1}^\infty\frac{1}{(\beta|\Lambda|)^{2l-1}}\sum_{\underline\alpha=(\alpha_1,\cdots,\alpha_{2l})}\sum_{\displaystyle\mathop{\scriptstyle (\mathbf k_1,\cdots,\mathbf k_{2l})\in\mathcal B_{\beta,L}^{(\leqslant h)2l}}_{\mathbf k_1-\mathbf k_2+\cdots+\mathbf k_{2l-1}-\mathbf k_{2l}=0}}\hskip-15pt\hat W_{2l,\underline\alpha}^{(h)}(\mathbf k_1,\cdots,\mathbf k_{2l-1})\cdot\\[1.5cm]
\hfill\cdot\hat\psi^{(\leqslant h)+}_{\mathbf k_1,\alpha_1}\hat\psi^{(\leqslant h)-}_{\mathbf k_2,\alpha_2}\cdots\hat\psi^{(\leqslant h)+}_{\mathbf k_{2l-1},\alpha_{2l-1}}\hat\psi^{(\leqslant h)-}_{\mathbf k_{2l},\alpha_{2l}}.
\end{largearray}\label{hatWdefo}\end{equation}
The $\mathbf x$-space expression for $\hat\psi_{\mathbf k,\alpha}^{(\leqslant h)\pm}$ is defined as
\begin{equation}
\psi_{\mathbf x,\alpha}^{(\leqslant h)\pm}:=\frac{1}{\beta|\Lambda|}\sum_{\mathbf k\in\mathcal B_{\beta,L}^{(\leqslant h)}}e^{\pm i\mathbf k\cdot\mathbf x}\hat\psi_{\mathbf k,\alpha}^{(\leqslant h)\pm}
\label{psixdefo}\end{equation}
so that the propagator's formulation in $\mathbf x$-space is
\begin{equation}
g_{h}(\mathbf x-\mathbf y):=\frac{1}{\beta|\Lambda|}\sum_{\mathbf k\in\mathcal B_{\beta,L}^{(\leqslant h)}}e^{i\mathbf k\cdot(\mathbf x-\mathbf y)}\hat g_h(\mathbf k)
\label{propxuv}\end{equation}
and similarly for $g_{\leqslant h}$, and the effective potential~(\ref{hatWdefo}) becomes
\begin{equation}\begin{largearray}
\mathcal V^{(h)}(\psi^{(\leqslant h)})=\sum_{l=1}^\infty\sum_{\underline\alpha}\frac{1}{\beta|\Lambda|}\int d\mathbf x_1\cdots d\mathbf x_{2l}\ W^{(h)}_{2l,\underline\alpha}(\mathbf x_1-\mathbf x_{2l},\cdots,\mathbf x_{2l-1}-\mathbf x_{2l})\cdot\\[0.5cm]
\hfill\cdot\psi^{(\leqslant h)+}_{\mathbf x_1,\alpha_1}\psi^{(\leqslant h)-}_{\mathbf x_2,\alpha_2}\cdots\psi^{(\leqslant h)+}_{\mathbf x_{2l-1},\alpha_{2l-1}}\psi^{(\leqslant h)-}_{\mathbf x_{2l},\alpha_{2l}}
\end{largearray}\label{Wdefo}\end{equation}
with
\begin{equation}
W^{(h)}_{2l,\underline\alpha}(\mathbf u_1,\cdots,\mathbf u_{2l-1})\\[0.2cm]
\hfill:=\frac{1}{(\beta|\Lambda|)^{2l-1}}\sum_{(\mathbf k_1,\cdots,\mathbf k_{2l-1})\in\mathcal B_{\beta,L}^{2l-1}}\hskip-25pt e^{i(\sum_{i=1}^{2l-1}(-1)^i\mathbf k_i\cdot\mathbf u_i)}\hat W^{(h)}_{2l,\underline\alpha}(\mathbf k_1,\cdots,\mathbf k_{2l-1}).
\label{xspaceWuv}\end{equation}
\bigskip
{\bf Remark}: From (\ref{hatWdefo}), $\hat W_{2l,\underline\alpha}^{(h)}(\underline{\mathbf k})$ is not defined for $\mathbf k_i\not\in\mathcal B_{\beta,L}^{(\leqslant h)}$, however, one can easily check that~(\ref{xspaceWuv}) holds for any extension of $\hat W_{2l,\underline\alpha}^{(h)}$ to $\mathcal B_{\beta,L}^{2l-1}$, thanks to the compact support properties of
$\psi^{(\leqslant h)}$ in momentum space. In order to get satisfactory bounds on $W_{2l,\underline\alpha}^{(h)}(\underline{\mathbf x})$, that is in order to avoid Gibbs phenomena, we define the extension of $\hat W_{2l,\underline\alpha}^{(h)}(\underline{\mathbf k})$ similarly to~(\ref{hatWdefo}) by relaxing the condition that $\psi^{(\leqslant h)}$ is supported on $\mathcal B_{\beta,L}^{(\leqslant h)}$ and iterating~(\ref{effpotuv}). In other words, we let $\hat W_{2l,\underline\alpha}^{(h)}(\underline{{\bf k}})$ for $\underline{{\bf k}}\in \mathcal B_{\beta,L}^{2l-1}$ be the kernels of $\mathcal V^{*(h)}$ defined inductively by
\begin{equation}
-\beta|\Lambda|\mathfrak e_h-\mathcal V^{*(h)}(\Psi):=\sum_{N=1}^\infty\frac{(-1)^N}{N!}\mathcal E_{h+1}^T(\mathcal V^{*(h+1)}(\psi^{(h+1)}+\Psi);N)
\label{extendedVuv}\end{equation}
in which $\{\hat\Psi_{\mathbf k,\alpha}\}_{\mathbf k\in\mathcal B_{\beta,L},\alpha\in\mathcal A}$ is a collection of {\it external fields} (in reference to the fact that, contrary to $\psi^{(\leqslant h)}$, they have a non-compact support in momentum space). The use of this specific extension can be justified {\it ab-initio} by re-defining the cutoff function $\chi$ in such a way that its support is $\mathbb{R}$, e.g. using exponential tails that depend on a parameter $\epsilon_\chi$ in such a way that the support tends to be compact as $\epsilon_\chi$ goes to 0. Following this logic, we could first define $\hat W$ using the non-compactly supported cutoff function and then take the $\epsilon_\chi\to0$ limit, thus recovering (\ref{extendedVuv}). Such an approach is dicussed in~[\cite{BM02}]. From now on, with some abuse of notation, we shall identify $\mathcal V^{*(h)}$ with $\mathcal V^{(h)}$ and denote them by the same symbol $\mathcal V^{(h)}$, which is justified by the fact that their kernels are (or can be chosen, from what said above, to be) the same.\par
\bigskip
\point{\bf First and second regimes.} We now discuss the first and second regimes (the third regime is very slightly different in that the index $\omega$ is complemented by an extra index $j$ and the Fermi points are shifted). Similarly to~(\ref{hatWdefo}), we define the {\it kernels} of $\bar{\mathcal V}$:
\begin{equation}\begin{largearray}
\bar{\mathcal V}^{(h)}(\psi^{(\leqslant h)})=:\sum_{l=2}^\infty\frac{1}{(\beta|\Lambda|)^{2l-1}}\sum_{\underline\omega,\underline\alpha}\sum_{\displaystyle\mathop{\scriptstyle (\mathbf k_1,\cdots,\mathbf k_{2l})\in\mathcal B_{\beta,L}^{(\leqslant h,\underline\omega)}}_{\mathbf k_1-\mathbf k_2+\cdots+\mathbf k_{2l-1}-\mathbf k_{2l}=0}}\hskip-15pt\hat W_{2l,\underline\alpha}^{(h)}(\mathbf k_1,\cdots,\mathbf k_{2l-1})\cdot\\[1.5cm]
\hfill\cdot\hat\psi^{(\leqslant h)+}_{\mathbf k_1,\alpha_1,\omega_1}\hat\psi^{(\leqslant h)-}_{\mathbf k_2,\alpha_2,\omega_2}\cdots\hat\psi^{(\leqslant h)+}_{\mathbf k_{2l-1},\alpha_{2l-1},\omega_{2l-1}}\hat\psi^{(\leqslant h)-}_{\mathbf k_{2l},\alpha_{2l},\omega_{2l}}.
\end{largearray}\label{hatWdeft}\end{equation}
where $\mathcal B_{\beta,L}^{(\leqslant h,\underline\omega)}=\mathcal B_{\beta,L}^{(\leqslant h,\omega_1)}\times\cdots\times\mathcal B_{\beta,L}^{(\leqslant h,\omega_{2l})}$. Note that the kernel $\hat W_{2l,\underline\alpha}^{(h)}$ is independent of $\underline\omega$, which can be easily proved using the symmetry $\omega_i\mapsto-\omega_i$. The $\mathbf x$-space expression for $\hat\psi_{\mathbf k,\alpha,\omega}^{(\leqslant h)\pm}$ is
\begin{equation}
\psi_{\mathbf x,\alpha,\omega}^{(\leqslant h)\pm}:=\frac{1}{\beta|\Lambda|}\sum_{\mathbf k\in\mathcal B_{\beta,L}^{(\leqslant h,\omega)}}e^{\pm i(\mathbf k-\mathbf p_{F,0}^\omega)\cdot\mathbf x}\hat\psi_{\mathbf k,\alpha,\omega}^{(\leqslant h)\pm}.
\label{psixdeft}\end{equation}
\bigskip
{\bf Remark}: Unlike $\hat\psi_{\mathbf k,\alpha,\omega}$, the $_\omega$ index in $\psi_{\mathbf x,\alpha,\omega}^{(\leqslant h)\pm}$ is {\it not} redundant. Keeping track of this dependence is required to prove properties of $\hat W_{2l}(\mathbf k)$ and $\hat{\bar g}_h(\mathbf k)$ close to $\mathbf p_{F,0}^\omega$ while working in $\mathbf x$-space. Such considerations were first discussed in~[\cite{BG90}] in which $\psi_{\mathbf x,\alpha,\omega}$ were called {\it quasi-particle} fields.\par
\bigskip
\indent We then define the propagator in $\mathbf x$-space:
\begin{equation}
\hat g_{h,\omega}(\mathbf x-\mathbf y):=\frac{1}{\beta|\Lambda|}\sum_{\mathbf k\in\mathcal B_{\beta,L}^{(\leqslant h,\omega)}}e^{i(\mathbf k-\mathbf p_{F,0}^\omega)\cdot(\mathbf x-\mathbf y)}\hat{\bar g}_{h,\omega}(\mathbf k)
\label{propxoh}\end{equation}
and similarly for $\bar g_{\leqslant h,\omega}$, and the effective potential~(\ref{hatWdeft}) becomes
\begin{equation}\begin{largearray}
\bar{\mathcal V}^{(h)}(\psi^{(\leqslant h)})=\sum_{l=2}^\infty\sum_{\underline\omega,\underline\alpha}\frac{1}{\beta|\Lambda|}\int d\mathbf x_1\cdots d\mathbf x_{2l}\ W^{(h)}_{2l,\underline\alpha,\underline\omega}(\mathbf x_1-\mathbf x_{2l},\cdots,\mathbf x_{2l-1}-\mathbf x_{2l})\cdot\\[0.5cm]
\hfill\cdot\psi^{(\leqslant h)+}_{\mathbf x_1,\alpha_1,\omega_1}\psi^{(\leqslant h)-}_{\mathbf x_2,\alpha_2,\omega_2}\cdots\psi^{(\leqslant h)+}_{\mathbf x_{2l-1},\alpha_{2l-1},\omega_{2l-1}}\psi^{(\leqslant h)-}_{\mathbf x_{2l},\alpha_{2l},\omega_{2l}}
\end{largearray}\label{Wdeft}\end{equation}
and
\begin{equation}
\mathcal Q^{(h)}(\psi^{(\leqslant h)})=\sum_{\omega,(\alpha,\alpha')}\int d\mathbf xd\mathbf y\ \psi^{(\leqslant h)+}_{\mathbf x,\omega,\alpha}W^{(h)}_{2,\omega,(\alpha,\alpha')}(\mathbf x-\mathbf y)\psi^{(\leqslant h)-}_{\mathbf y,\omega,\alpha'}
\label{Wtdeft}\end{equation}
in which
\begin{equation}\begin{largearray}
W^{(h)}_{2l,\underline\alpha,\underline\omega}(\mathbf u_1,\cdots,\mathbf u_{2l-1})\\[0.5cm]
\hfill:=\frac{\delta_{0,\sum_{j=1}^{2l}(-1)^{j}\mathbf p_{F,0}^{\omega_j}}}{(\beta|\Lambda|)^{2l-1}}\sum_{(\mathbf k_1,\cdots,\mathbf k_{2l-1})\in\mathcal B_{\beta,L}^{2l-1}} e^{i(\sum_{j=1}^{2l-1}(-1)^j(\mathbf k_j-\mathbf p_{F,0}^{\omega_j})\cdot\mathbf u_j)}\hat W^{(h)}_{2l,\underline\alpha}(\mathbf k_1,\cdots,\mathbf k_{2l-1}).
\end{largearray}\label{xspaceWo}\end{equation}
As in the ultraviolet, the definition of $\hat W_{2l,\underline\alpha}^{(h)}(\underline{\mathbf k})$ is extended to $\mathcal B_{\beta,L}^{2l-1}$ by defining it as the kernel of $\mathcal V^{*(h)}$:
\begin{equation}
-\beta|\Lambda|\mathfrak e_h-\mathcal V^{*(h)}(\Psi):=\sum_{N=1}^\infty\frac{(-1)^N}{N!}\bar{\mathcal E}_{h+1}^T(\mathcal V^{*(h+1)}(\psi^{(h+1)}+\Psi);N)
\label{extendedVo}\end{equation}
in which $\{\hat\Psi_{\mathbf k,\alpha}\}_{\mathbf k\in\mathcal B_{\beta,L},\alpha\in\mathcal A}$ is a collection of {\it external fields}. The definition~(\ref{xspaceWo}) suggests a definition for $\bar A_{h,\omega}$ (see~(\ref{Adefo}) and~(\ref{Adeft})):
\begin{equation}
\bar A_{h,\omega}(\mathbf x):=\frac{1}{\beta|\Lambda|}\sum_{\mathbf k\in\mathcal B_{\beta,L}}e^{i(\mathbf k-\mathbf p_{F,0}^\omega)\cdot\mathbf x}\hat{\bar A}_{h,\omega}(\mathbf k).
\label{Ahoox}\end{equation}
\bigskip
\point{\bf Third regime.} We now turn our attention to the third regime. As discussed in section~\ref{multiscalesec}, in addition to there being an extra index $j$, the Fermi points are also shifted in the third regime. The kernels of $\bar{\mathcal V}$ and $\mathcal Q$ are defined as in~(\ref{hatWdeft}), but with $\omega$ replaced by $(\omega,j)$. The $\mathbf x$-space representation of $\hat\psi_{\mathbf k,\alpha,\omega,j}^{(\leqslant h)\pm}$ is defined as
\begin{equation}
\psi_{\mathbf x,\alpha,\omega,j}^{(\leqslant h)\pm}:=\frac{1}{\beta|\Lambda|}\sum_{\mathbf k\in\mathcal B_{\beta,L}^{(\leqslant h,\omega,j)}}e^{\pm i(\mathbf k-\tilde{\mathbf p}_{F,j}^{(\omega,h)})\cdot\mathbf x}\hat\psi_{\mathbf k,\alpha,\omega,j}^{(\leqslant h)\pm}
\label{psixdefth}\end{equation}
and the $\mathbf x$-space expression of the propagator and the kernels of $\bar{\mathcal V}$ and $\mathcal Q$ are defined by analogy with the first regime:
\begin{equation}
\hat g_{h,\omega,j}(\mathbf x-\mathbf y):=\frac{1}{\beta|\Lambda|}\sum_{\mathbf k\in\mathcal B_{\beta,L}^{(\leqslant h,\omega,j)}}e^{i(\mathbf k-\tilde{\mathbf p}_{F,j}^{(\omega,h)})\cdot(\mathbf x-\mathbf y)}\hat{\bar g}_{h,\omega,j}(\mathbf k)
\label{propxothh}\end{equation}
and
\begin{equation}\begin{largearray}
W^{(h)}_{2l,\underline\alpha,\underline\omega,\underline j}(\mathbf u_1,\cdots,\mathbf u_{2l-1})
:=\frac{\delta_{0,\sum_{n=1}^{2l}(-1)^{n}\tilde{\mathbf p}_{F,j}^{(h,\omega_n)}}}{(\beta|\Lambda|)^{2l-1}}\cdot\\[0.5cm]
\hfill\cdot\sum_{(\mathbf k_1,\cdots,\mathbf k_{2l-1})\in\mathcal B_{\beta,L}^{2l-1}} e^{i(\sum_{n=1}^{2l-1}(-1)^n(\mathbf k_n-\tilde{\mathbf p}_{F,j}^{(\omega_n,h)})\cdot\mathbf u_j)}\hat W^{(h)}_{2l,\underline\alpha}(\mathbf k_1,\cdots,\mathbf k_{2l-1}).
\end{largearray}\label{xspaceWth}\end{equation}
In addition
\begin{equation}
\bar A_{h,\omega,j}(\mathbf x):=\frac{1}{\beta|\Lambda|}\sum_{\mathbf k\in\mathcal B_{\beta,L}}e^{i(\mathbf k-\tilde{\mathbf p}_{F,j}^{(\omega,h)})\cdot\mathbf x}\hat{\bar A}_{h,\omega,j}(\mathbf k).
\label{Ahojx}\end{equation}
\subsection{Estimates of the free propagator}
\label{estpropsubsec}
\indent Before moving along with the tree expansion, we first compute a bound on $\hat g_{h}$ in the different regimes, which will be used in the following.\par
\bigskip
\point{\bf Ultraviolet regime.} We first study the ultraviolet regime, i.e. $h\in\{1,\cdots,M\}$.\par
\medskip
\subpoint{\bf Fourier space bounds.} We have
$$\hat A(\mathbf k)^{-1}:=-(ik_0\mathds1+H_0(k))^{-1}=-\frac{1}{ik_0}\left(\mathds1+\frac{H_0(k)}{ik_0}\right)^{-1}$$
and
$$|\hat g_{h}(\mathbf k)|=|f_h(\mathbf k)\hat A^{-1}(\mathbf k)|\leqslant (\mathrm{const.})\ 2^{-h},$$
where $|\cdot|$ is the operator norm.
Therefore
\begin{equation}
\frac{1}{\beta|\Lambda|}\sum_{\mathbf k\in\mathcal B_{\beta,L}^*}|\hat g_{h}(\mathbf k)|\leqslant (\mathrm{const}.).
\label{grambounduv}\end{equation}
Furthermore, for all $m_0+m_k\leqslant7$ (we choose the constant $7$ in order to get adequate bounds on the real-space decay of the free propagator,
good enough for performing the localization and renormalization procedure described below; any other larger constant would yield identical results),
\begin{equation}
|2^{hm_0}\partial_{k_0}^{m_0}\partial_k^{m_k}\hat g_{h}(\mathbf k)|\leqslant (\mathrm{const.})\ 2^{-h}
\label{estgkuv}\end{equation}
in which $\partial_{k_0}$ denotes the discrete derivative with respect to $k_0$ and, with a slightly abusive notation, $\partial_k$ the discrete derivative with respect to either $k_x$ or $k_y$. Indeed the derivatives over $k$ land on $ik_0\hat A^{-1}$, which does not change the previous estimate, and the derivatives over $k_0$ either land on $f_h$, $1/(ik_0)$, or $ik_0\hat A^{-1}$, which yields an extra $2^{-h}$ in the estimate.\par
\medskip
{\bf Remark}: The previous argument implicitly uses the Leibnitz rule, which must be used carefully since the derivatives are discrete. However, since the estimate is purely dimensional,
we can replace the discrete discrete derivative with a continuous one without changing the order of magnitude of the resulting bound.\par
\bigskip
\subpoint{\bf Configuration space bounds.} We now prove that the inverse Fourier transform of $\hat g_{h}$
\begin{equation} g_{h}(\mathbf x):=\frac{1}{\beta|\Lambda|}\sum_{\mathbf k\in\mathcal B^*_{\beta,L}}\ e^{-i\mathbf k\cdot\mathbf x}\hat g_{h}(\mathbf k)\label{5.41bis}\end{equation}
satisfies the following estimate: for all $m_0+m_k\leqslant3$,
\begin{equation}
\int d\mathbf x\ x_0^{m_0}x^{m_k}|g_{h}(\mathbf x)|\leqslant (\mathrm{const.})\ 2^{-h-m_0h},
\label{estguv}\end{equation}
where we recall that $\int d\mathbf x$ is a shorthand for $\int_{0}^\beta dt \sum_{x\in \Lambda}$.
Indeed, note that the right side of (\ref{5.41bis}) can be thought of as the Riemann sum approximation of
\begin{equation} \int_{\mathbb R}\frac{dk_0}{2\pi}\int_{\hat\Lambda_\infty} \frac{d{ k}}{|\hat\Lambda_\infty|} e^{-i\mathbf k\cdot\mathbf x}\hat g_{h}(\mathbf k)\label{5.42bis}\end{equation}
where $\hat\Lambda_\infty=\{t_1 G_1+t_2 G_2: t_i\in[0,1)\}$ is the limit as $L\to\infty$ of $\hat \Lambda$, see (\ref{lae}) and following lines. The dimensional estimates one finds using this continuum approximation are the same as those using (\ref{5.41bis}) therefore, integrating (\ref{5.42bis}) 7 times by parts and using (\ref{estgkuv}) we find
$$|g_{h}(\mathbf x)|\leqslant\frac{(\mathrm{const.})}{1+(2^{h}|x_0|+|x|)^7}$$
so that by changing variables in the integral over $x_0$ to $2^hx_0$, and using
$$\int d\mathbf x\ \frac{x_0^{m_0}x^{m_k}}{1+(|x_0|+|x|)^7}<(\mathrm{const.})$$
we find~(\ref{estguv}).\par
\bigskip
\point{\bf First regime.} We now consider the first regime, i.e. $h\in\{\mathfrak h_1+1,\cdots,\bar{\mathfrak h}_0\}$.\par
\medskip
\subpoint{\bf Fourier space bounds.} From~(\ref{freepropzo}) we find
$$|\hat g_{h,\omega}(\mathbf k)|\leqslant (\mathrm{const.})\ 2^{-h}$$
therefore
\begin{equation}
\frac{1}{\beta|\Lambda|}\sum_{\mathbf k\in\mathcal B_{\beta,L}^*}|\hat g_{h,\omega}(\mathbf k)|\leqslant (\mathrm{const}.)\ 2^{2h}
\label{gramboundo}\end{equation}
and for $m\leqslant7$,
\begin{equation}
|2^{mh}\partial_{\mathbf k}^m\hat g_{h,\omega}(\mathbf k)|\leqslant (\mathrm{const.})\ 2^{-h}
\label{estgko}\end{equation}
in which we again used the slightly abusive notation of writing $\partial_{\mathbf k}$ to mean any derivative with respect to $k_0$, $k_x$ or $k_y$. Equation~(\ref{estgko}) then follows from similar considerations as those in the ultraviolet regime.\par
\bigskip
\subpoint{\bf Configuration space bounds.} We estimate the real-space counterpart of $\hat g_{h,\omega}$,
$$g_{h,\omega}(\mathbf x):=\frac{1}{\beta|\Lambda|}\sum_{\mathbf k\in\mathcal B_{\beta,L}^{(h,\omega)}}e^{-i(\mathbf k-\mathbf p_{F,0}^\omega)\cdot\mathbf x}\hat g_{h,\omega}(\mathbf k),$$
and find that for $m\leqslant3$,
\begin{equation}
\int d\mathbf x\ |\mathbf x^mg_{h,\omega}(\mathbf x)|\leqslant (\mathrm{const.})\ 2^{-(1+m)h}
\label{estgo}\end{equation}
which follows from very similar considerations as the ultraviolet estimate.\par
\bigskip
\point{\bf Second regime.} We treat the second regime, i.e. $h\in\{\mathfrak h_2+1,\cdots,\bar{\mathfrak h}_1\}$ in a very similar way (we skip the intermediate regime which can be treated in the same way as either the first or second regimes):
\begin{equation}
\frac{1}{\beta|\Lambda|}\sum_{\mathbf k\in\mathcal B_{\beta,L}^*}|\hat g_{h,\omega}(\mathbf k)|\leqslant (\mathrm{const}.)\ 2^{h+h_\epsilon}
\label{gramboundt}\end{equation}
and for all $m_0+m_k\leqslant7$,
\begin{equation}
\big|2^{m_0h}\partial_{k_0}^{m_0}2^{m_k\frac{h+h_\epsilon}2}\partial_k^{m_k}\hat g_{h,\omega}(\mathbf k)\big|\leqslant (\mathrm{const.})\ 2^{-h}
\label{estgkt}\end{equation}
where $h_\epsilon:=\log_2(\epsilon)$. Therefore for all $m_0+m_k\leqslant 3$,
\begin{equation}
\int d\mathbf x\ |x_0^{m_0}x^{m_k}g_{h,\omega}(\mathbf x)|\leqslant (\mathrm{const.})\ 2^{-h-m_0h-m_k\frac{h+h_\epsilon}2}.
\label{estgt}\end{equation}
\bigskip
\point{\bf Third regime.} Finally, the third regime, i.e. $h\in\{\mathfrak h_3+1,\cdots,\bar{\mathfrak h}_2\}$:
\begin{equation}
\frac{1}{\beta|\Lambda|}\sum_{\mathbf k\in\mathcal B_{\beta,L}^*}|\hat g_{h,\omega}(\mathbf k)|\leqslant (\mathrm{const}.)\ 2^{2h-2h_\epsilon}
\label{gramboundth}\end{equation}
and for all $m_0+m_k\leqslant7$,
\begin{equation}
|2^{m_0h}\partial_{k_0}^{m_0}2^{m_k(h-h_\epsilon)}\hat g_{h,\omega,j}(\mathbf k)|\leqslant (\mathrm{const.})\ 2^{-h}.
\label{estgkth}\end{equation}
Therefore for all $m_0+m_k\leqslant 3$,
\begin{equation}
\int d\mathbf x\ |x_0^{m_0}x^{m_k}g_{h,\omega,j}(\mathbf x)|\leqslant (\mathrm{const.})\ 2^{-h-m_0h-m_k(h-h_\epsilon)}
\label{estgth}\end{equation}
where
$$g_{h,\omega,j}(\mathbf x):=\frac{1}{\beta|\Lambda|}\sum_{\mathbf k\in\mathcal B_{\beta,L}^{(h,\omega,j)}}e^{-i(\mathbf k-\tilde{\mathbf p}_{F,j}^{(\omega,h+1)})\cdot\mathbf x}\hat g_{h,\omega}(\mathbf k).$$
\section{Tree expansion and constructive bounds}
\label{treeexpsec}
\indent In this section, we shall define the Gallavotti-Nicol\`o tree expansion~[\cite{GN85}], and show how it can be used to compute bounds for the $\mathfrak e_h$, $\mathcal V^{(h)}$,
$\mathcal Q^{(h)}$ and $\bar{\mathcal V}^{(h)}$ defined above in~(\ref{effpotuv}) and (\ref{effpotoh}), using the estimates~(\ref{estguv}), (\ref{estgo}), (\ref{estgt}) and~(\ref{estgth}). We follow~[\cite{BG90}, \cite{GM01}, \cite{GM10}].
We conclude the section by showing how to compute the terms in $\bar{\mathcal W}^{(h)}$ that are quadratic in $\hat J_{\mathbf k,\underline\alpha}$ from $\mathcal V^{(h)}$ and $\hat{\bar g}_h$.\par
\bigskip
\indent The discussion in this section is meant to be somewhat general, in order to be applied to the ultraviolet, first, second and third regimes (except for lemma~\ref{powercountinglemma} which does not apply to the ultraviolet regime).\par
\subsection{Gallavotti-Nicol\`o Tree expansion}
\indent In this section, we will define a tree expansion to re-express equations {\it of the type}
\begin{equation}
-v^{(h)}(\psi^{(\leqslant h)})-\mathcal V^{(h)}(\psi^{(\leqslant h)})=\sum_{N=1}^\infty\frac{(-1)^N}{N!}\mathcal E^T_{h+1}\left(\mathcal V^{(h+1)}(\psi^{(\leqslant h)}+\psi^{(h+1)});N\right)
\label{treeindeq}\end{equation}
for $h\in\{h_2^*,\cdots,h_1^*-1\}$ (in the ultraviolet regime $h_2^*=\bar{\mathfrak h}_0$, $h_1^*=M$; in the first $h_2^*=\mathfrak h_1$, $h_1^*=\bar{\mathfrak h}_0$; in the second $h_2^*=\mathfrak h_2$, $h_1^*=\bar{\mathfrak h}_1$;
and in the third, $h_2^*=\mathfrak h_\beta$, $h_1^*=\bar{\mathfrak h}_2$), with
\begin{equation} \left\{\begin{array}{>{\displaystyle}l}
\mathcal V^{(h)}(\psi^{(\leqslant h)})=\sum_{l=q}^\infty\sum_{\underline\varpi}\int d\underline{\mathbf x}\ W_{2l,\underline\varpi}^{(h)}(\underline{\mathbf x})\psi^{(\leqslant h)+}_{\mathbf x_1,\varpi_1}\psi^{(\leqslant h)-}_{\mathbf x_2,\varpi_2}\cdots\psi^{(\leqslant h)+}_{\mathbf x_{2l-1},\varpi_{2l-1}}\psi^{(\leqslant h)-}_{\mathbf x_{2l},\varpi_{2l}}\\[0.5cm]
v^{(h)}(\psi^{(\leqslant h)})=\sum_{l=0}^{q-1}\sum_{\underline\varpi}\int d\underline{\mathbf x}\ W_{2l,\underline\varpi}^{(h)}(\underline{\mathbf x})\psi^{(\leqslant h)+}_{\mathbf x_1,\varpi_1}\psi^{(\leqslant h)-}_{\mathbf x_2,\varpi_2}\cdots\psi^{(\leqslant h)+}_{\mathbf x_{2l-1},\varpi_{2l-1}}\psi^{(\leqslant h)-}_{\mathbf x_{2l},\varpi_{2l}}
\end{array}\right.\label{defW}\end{equation}
($q=1$ in the ultraviolet regime and $q=2$ in the first, second and third) in which $\underline\varpi$ and $\underline{\mathbf x}$ are shorthands for $(\varpi_1,\cdots,\varpi_{2l})$ and $(\mathbf x_1,\cdots,\mathbf x_{2l})$; $\varpi$ denotes a collection of indices: $(\alpha,\omega)$ in the first and second regimes, $(\alpha,\omega,j)$ in the third, and $(\alpha)$ in the ultraviolet; and $W_{2l,\underline\varpi}^{(h)}(\underline{\mathbf x})$ is a function that only depends on the differences $\mathbf x_i-\mathbf x_j$.
The propagator associated with $\mathcal E^T_{h+1}$ will be denoted $g_{(h+1),(\varpi,\varpi')}(