Interaction Corrections to Spin-Wave Theory in the Large-S Limit of the Quantum Heisenberg Ferromagnet
Niels Benedikter\vspace{1em
%Understanding the low-temperature properties of the Quantum Heisenberg Ferromagnet is central to a mathematical understanding of magnetism. The Quantum Heisenberg Ferromagnet can be naturally reformulated in terms of interacting bosons (called spin waves or magnons) as an expansion in the inverse spin size. We calculate the first order interaction correction to the free energy, as an upper bound in the limit where the spin size $S \to \infty$ and $\beta S$ is fixed ($\beta$ being the inverse temperature). Our result is valid in two and three spatial dimensions. We extrapolate our result to compare with Dyson's low-temperature expansion. While our first-order correction has the expected temperature dependance, in higher orders of the perturbation theory cancellations are necessary.
heisenberg.tex
\documentclass[a4paper]{article}
\usepackage[a4paper,hmargin=2.6cm,vmargin=2.8cm]{geometry}
\usepackage[utf8x]{inputenc}
\usepackage{amsmath,amsthm,amsfonts,amssymb}
\usepackage{feynmp}
\usepackage[bookmarksnumbered=true]{hyperref}
% \usepackage[color]{showkeys}	% not for final version
 \usepackage{dsfont}
%\usepackage{caption}
\usepackage{graphicx}

\DeclareGraphicsRule{*}{mps}{*}{}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\include{macros}
\renewcommand{\vec}[1]{\textbf{#1}}
\newcommand{\gs}{\lvert\text{gs}\rangle}
%%%%%%%%%%%%%%%%%%%%%%%%% main content %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%\title{A Remark on Interaction Effects in Spin-Wave Theory}
\title{Interaction Corrections to Spin-Wave Theory in the Large-S Limit of the Quantum Heisenberg Ferromagnet}

\author{Niels Benedikter\vspace{1em}\\QMATH\\Department of Mathematical Sciences,  University of Copenhagen\\Universitetsparken 5, 2100 Copenhagen, Denmark\vspace{.25em}}

\begin{document}
\begin{fmffile}{diags}
\maketitle

\begin{abstract}
%Understanding the low-temperature properties of the Quantum Heisenberg Ferromagnet is central to a mathematical understanding of magnetism.
The Quantum Heisenberg Ferromagnet can be naturally reformulated in terms of interacting bosons (called spin waves or magnons) as an expansion in the inverse spin size. We calculate the first order interaction correction to the free energy, as an upper bound in the limit where the spin size $S \to \infty$ and $\beta S$ is fixed ($\beta$ being the inverse temperature). Our result is valid in two and three spatial  dimensions.

We extrapolate our result to compare with Dyson's low-temperature expansion. While our first-order correction has the expected temperature dependance, in higher orders of the perturbation theory cancellations are necessary.
\end{abstract}


\section{Introduction and Main Result}
We consider the ferromagnetic quantum Heisenberg model on a volume $\Lambda_L = [0,L]^d \cap \Zbb^d$ in spatial dimension $d=2$ or $3$, with nearest neighbour interaction. It is described by the Hamilton operator
\[H_{\Lambda_L} := \sum_{\langle x, y\rangle \subset \Lambda_L}\left(S^2 - {{S}}_{x} \cdot {{S}}_{y} \right)\]
acting on the Hilbert space $\bigotimes_{x\in\Lambda_L} \Cbb^{2S+1}$.
The sum is over unordered nearest neighbour pairs in $\Lambda_L$, and ${{S}}_x = (S^1_x,S^2_x,S^3_x)$ is a spin-$S$ operator (where $2S \in \Nbb$), i.\,e.\ satisfying the commutation relations
$[{S}^j_{x},{S}^k_{y}] = \delta_{x,y}\,i \sum_{l=1}^3 \varepsilon_{jkl} {S}^l_{x}$ for $j,k \in \{1,2,3\},\ x,y \in \Lambda_L$,
(with $\varepsilon_{jkl}$ the totally antisymmetric symbol) and the condition $({{S}}_{x})^2 = ( S^1_{x})^2 + ( S^2_{x})^2 + ( S^3_{x})^2 = S(S+1)$. %; this is natural for an upper bound and consistent with our method of proof.
The quantity we are interested in is the free energy in the thermodynamic limit,
\[f(S,\beta) := \lim_{\Lambda_L \to \infty} f(S,\beta,\Lambda_L), \quad f(S,\beta,\Lambda_L) := -\frac{1}{\beta L^d} \log \tr e^{-\beta H_{\Lambda_L}},\]
where $\beta$ is the inverse temperature.
This model is of great importance for the understanding of ferromagnetism, which poses a challenging mathematical problem, namely the spontaneous breaking of continuous symmetries.

\medskip

In this note we are interested in the spin-wave approximation, which goes back to Bloch in 1930 \cite{Bloch1930}. Bloch noticed that the low-energy excitations of the Heisenberg model can be approximately described as independent bosonic modes with energy given by the dispersion relation $\varepsilon(k) := \sum_{j=1}^d 2(1-\cos(k_j))$, and thus was able to predict the behaviour of the thermodynamic quantities at low temperature.

Let us recall this theory in modern language. The ground states of the Heisenberg model at zero temperature are states in which all spins point parallel, while the direction is arbitrary. W.\,l.\,o.\,g.\ we can take the ground state to have all spins pointing in the $-z$-direction, i.\,e.\ 
\[\gs = \bigotimes_{x\in\Lambda_L} \lvert -S\rangle_x,\quad\text{where } \lvert -S\rangle_x \in \Cbb^{2S+1},\quad S^3_x \lvert -S\rangle_x = -S \lvert -S\rangle_x.\]
Spin waves (or magnons) are excitations of the form
\[\lvert k\rangle := (2S L^d)^{-1/2}\sum_{x\in\Lambda_L} e^{ik\cdot x} S^+_x \gs =: (2S)^{-1/2} S^+_k \gs,\]
with a momentum $k \in \frac{2\pi}{L}\Zbb^d$ (taking in this introduction for simplicity periodic boundary conditions). These states are normalized eigenstates of the Hamiltonian, 
\[H_{\Lambda_L}\lvert k \rangle = S \varepsilon(k) \lvert k\rangle.\]
Treating the system as a system of non-interacting bosons with energy $S\varepsilon(k)$, one arrives at the following prediction for the free energy:
\begin{equation}\label{eq:freeenergy}f(S,\beta) \simeq \frac{1}{\beta} \int_{[-\pi,\pi]^d} \frac{\di^d k}{(2\pi)^d} \log(1-e^{-\beta S \varepsilon(k)}) \sim \beta^{-\frac{d+2}{2}} S^{-\frac{d}{2}} \int \frac{\di^d k}{(2\pi)^d} \log( 1 - e^{-k^2}) \quad \text{as } \beta \to \infty.\end{equation}
Unfortunately, the spin-wave excitations do not actually behave like independent bosonic modes; the states with more than one spin wave---i.\,e.\ constructed by applying more than one operators $S^+_k$---are neither eigenstates of the Hamiltonian nor orthogonal.
 This problem is treated more systematically by the Holstein-Primakoff mapping \cite{HP40} of spin operators on bosonic operators: defining the spin raising and lowering operators $ S^+_x :=  S^1_x + i  S^2_x$ and $ S^-_x :=  S^1_x - i  S^2_x$, one introduces the mapping on bosonic creation and annihilation operators $a^*_x$ and $a_x$ by
\begin{equation}\label{eq:hpmapping} S^+_x = \sqrt{2S}a^*_x \left[ 1-\frac{a^*_x a_x}{2S} \right]^{1/2}, \quad S^-_x = \sqrt{2S}\left[ 1-\frac{a^*_x a_x}{2S} \right]^{1/2} a_x, \quad S^3_x = a^*_x a_x - S.\end{equation}
The creation and annihilation operators act on the subspace of bosonic Fock space where the number of bosons per lattice site $x \in \Lambda_L$, $n_x := a^*_x a_x$, is restricted to be at most $2S$. The Hamiltonian becomes
\begin{equation}\label{eq:hphamiltonian}H_{\Lambda_L} = S \sum_{\langle x,y\rangle \subset \Lambda_L} \left( -a^*_x\sqrt{1-\frac{n_x}{2S}}\sqrt{1-\frac{n_y}{2S}} a_y - a^*_y \sqrt{1-\frac{n_y}{2S}}\sqrt{1-\frac{n_x}{2S}} a_x + n_x + n_y - \frac{1}{S} n_x n_y\right).\end{equation}
The Hamiltonian is then formally expanded in $\frac{n_x}{2S}$, a procedure that is expected to lead to good approximations if $S$ is large or if the expected occupation numbers are small, i.\,e.\ at low temperature. The leading term is given by
\[H_{\Lambda_L} \simeq S \sum_{\langle x, y\rangle\subset\Lambda_L} (a^*_x - a^*_y)(a_x - a_y),\]
the second quantized Laplacian on the lattice, giving rise to the free-boson picture and the expressions \eqref{eq:freeenergy}. In this paper, we are interested in the residual interaction between spin waves as given through the higher orders terms of the expansion of the Hamiltonian. The corrections due to interactions have given rise to considerable discussions in the physics community, a topic that we shall discuss further after the statement of our theorem.

\medskip

There are three scaling regimes which are important in the study of the Heisenberg model. Physically most important is the limit of low temperature $\beta \to \infty$ and fixed $S$; it is however very difficult to study. In fact, only recently has the leading order of the free energy been rigorously derived \cite{CGS14} (non-optimal bounds were proved earlier in \cite{CS91,T93}). On the other extreme there is the classical scaling regime, defined by $\beta S^2$ fixed and $S\to \infty$. In this scaling, convergence to the classical Heisenberg model has been proven \cite{Lieb73} (for non-zero magnetic field also in \cite{CS90}). For the classical Heisenberg model, it was proven that it has a critical temperature of order unity \cite{FSS76}. This suggests that by going to the intermediate scaling regime
\begin{equation}\label{eq:intermediateregime}\tb := \beta S \text{ fixed}, \quad S \to \infty\end{equation}
we can study the ferromagnetic phase in a regime that is more accessible (since the attractive interaction between spin waves in this regime is of order $S^{-1}$) than the low temperature regime, and still governed by quantum theory. For $d=3$, the leading order of the free energy in this regime has been obtained in \cite{CG12}. (This regime was introduced and compared to the classical regime in \cite{CS90}, where the leading order in the case of non-zero magnetic field was derived. In this context also the important random walk representation of the Heisenberg model was developed \cite{CS91rw, CS91class}.)

\medskip

In this note we calculate an upper bound on the free energy which includes interaction effects to first order in the intermediate scaling regime \eqref{eq:intermediateregime}. Our result is valid in two and three dimensions.
%
\begin{thm*}\label{thm:mainresult}
Let $d = 2$ or $3$, and $\tb$ fixed and sufficiently large. Then the free energy is bounded above by
 \[\frac{f(S,\beta)}{S} \leq \frac{1}{\tilde \beta} \int_{[-\pi,\pi]^d} \frac{\di^d k}{(2\pi)^d} \log \left( 1-e^{-\tilde \beta \varepsilon(k)} \right) - \frac{1}{S} \frac{1}{4d} \left( \int_{[-\pi,\pi]^d} \frac{\di^d k}{(2\pi)^d} \frac{\varepsilon(k)}{e^{\tilde \beta \varepsilon(k)}-1}\right)^2 + \frac{r_d(S,\tb)}{S^2}. \]
 for $S \to \infty$, where $\varepsilon(k) = \sum_{j=1}^d 2(1-\cos(k_j))$ is the spin-wave dispersion relation.
 
 The error term\footnote{We use $C$ for constants independent of $S$ and $\tb$ (and $\ell$, to be introduced later), with value possibly changing from line to line.} is $r_2(S,\tb) \leq C \tb^{-2} (\log S\tb)^3$ and $r_3(S,\tb) \leq C \tb^{-3}$ (independent of $S$).
 \end{thm*}
% \begin{rem}
%  The biggest error is due to Lemma \ref{lem:II}; in particular it is responsible for the $(\log S)^3$ in $d=2$.
% \end{rem}
% \begin{rem}
% The leading-order \emph{lower} bound of \cite{CG12} also holds for $d=2$. \todo{check again}
% \end{rem}
%
The first summand---the leading order---describes free bosons (the spin waves) on the lattice. The second summand is the first order correction due to the interaction of spin waves. The interaction corrections and their  temperature dependence have long been controversial among physicists, with many contradictory corrections proposed, e.\,g.\ of order $\beta^{-3}$ or $\beta^{-11/4}$ \cite{Kramers36,Opechowski37,Schafroth54,Kranendonk54}. Eventually Dyson mostly settled the issue in his landmark papers \cite{Dyson1,Dyson2}, arguing (for $d=3$) that the correction is very small %\footnote{as an expansion in $\tilde\beta^{-1}$ with $S$ as a parameter, whereas our result expands in $S^{-1}$ with $\tilde\beta$ as a parameter}
 for low temperature\footnote{Here we are rewriting the equation  \cite[(131)]{Dyson2}. The constant summands $-\frac{1}{2}JS^2\gamma_0$ and $-LS$ given in Dyson's paper vanish in our setting, since we shifted the ground state energy by the $S^2$ in the Hamiltonian, and since we do not have an external magnetic field. The next three summands in Dyson's paper are an expansion of the leading integral written here. Dyson's $- k T C_3 S^{-1}[Z_{5/2}(\beta L)]^2 \theta^4$ corresponds to the term written here in square brackets; the $\mathcal{O}(S^{-2})$ is due to the $S$-dependence of $C_3$.}, namely of order $\tilde \beta^{-5}$: % for simple lattice and $S \to \infty$
\begin{equation}\label{eq:star}\begin{split}\frac{f(S,\beta)}{S} & = \frac{1}{\tilde \beta} \int_{[-\pi,\pi]^3} \frac{\di^3 k}{(2\pi)^3} \log \left( 1-e^{-\tilde \beta \varepsilon(k)} \right) - \left[ S^{-1}\frac{3}{128(2\pi)^3}\zeta(5/2)^2 + \mathcal{O}(S^{-2})\right] \tb^{-5} + \mathcal{O}(\tilde\beta^{-11/2}).\nn\end{split}\end{equation}
Nevertheless, the temperature dependence of the interaction corrections is being studied up to recent years, mostly substantiating Dyson's result by other formal methods, e.\,g.\ an effective Lagrangian method \cite{Hofmann02, Hofmann11}. 
One paper that should be highlighted is \cite{Zittartz65}, which lead to Dyson's result by a less cumbersome method of introducing additional bosonic degrees of freedom coupled to the spin system. 
However, there is also work newer than Dyson's papers which contradicts Dyson's result; see \cite{Hewson63} and the list of references therein.

%\todo{i used to read it as $\tilde \beta^{-6}$ before...}
% Before, the coefficient of the $S^{-1}$-term had been expected to be larger than $\tilde\beta^{-5}$.
To compare our result with Dyson's result we now formally think of $\tilde \beta \to \infty$ in our result. Expanding $\varepsilon(k)$ for small $k$ we get
\[-\frac{1}{S}\frac{1}{12} \left( \int_{[-\pi,\pi]^3} \frac{\di^3 k}{(2\pi)^3} \frac{\varepsilon(k)}{e^{\tilde \beta \varepsilon(k)}-1} \right)^2 \simeq -S^{-1}\frac{1}{12}  \frac{1}{(2\pi)^6}\left( \int_{\Rbb^3} \di^3 k' \frac{\lvert k'\rvert^2}{e^{\lvert k'\rvert^2}-1} \right)^2 \tb^{-5}.\]
Using spherical coordinates and $(e^{\lvert k'\rvert^2}-1)^{-1} = \sum_{n=1}^\infty e^{-n\lvert k'\rvert^2}$, we recover Dyson's result at order $S^{-1}$.

Of course this argument is beyond the proven validity of our theorem because $r_3(S,\tb) \sim \tb^{-3}$. Inspection of our proof shows that $r_3$ actually consists of two kinds of errors: corresponding to Dyson's kinematical interaction we have errors controlled by $\tb^{-S}$, and corresponding to Dyson's dynamical interaction we have our main error of order $S^{-2}\tb^{-3}$ (see Lemma $\ref{lem:II}$). However, there is a cancellation mechanism which is supposed (but not proven) to make the latter as small as $\tb^{-5}$ and which we discuss perturbatively in the Appendix.

\begin{rems}
\begin{enumerate}\item Our method also provides a partial result for the case of $\tb$ small, but not too small. As the prime example we choose $\tb = S^{-\alpha}$, $\alpha\in [0,1)$, and  with minor changes (see Remark (ii) after Lemma \ref{lem:oneminusp}) obtain
$r_3(S,\tb)/S^2 = \mathcal{O}(S^{3\alpha-2})$ for $S\to\infty$ (for $d=2$ with logarithmic correction). For comparison: in this case the leading term of $f(S,\beta)/S$ is $\tb^{-1}(c_0 \log \tb + c_1 + c_2\tb + \ldots ) \simeq S^\alpha \log S$ and the first order correction is $S^{-1}\tb^{-2} \simeq S^{2\alpha-1}$.
\item While for $d=3$ the validity of spin-wave theory has long been trusted in by physicists, it remained more disputed in $d=2$. Our result supports the validity of spin-wave theory in $d=2$, as far as the free energy in the intermediate scaling regime is concerned. Notice that also the (leading order) lower bound from \cite{CG12} is easily checked to be valid also for two dimensions.
\item Obtaining the first order correction as a \emph{lower} bound remains open; even in the intermediate scaling regime this is expected to be a very difficult problem.
\end{enumerate}
\end{rems}

%that we can not expect to make the error in our theorem smaller than $S^{-2} \sup_{x} \rho(x)^3$, which even with the improved bound $\rho(x) = \mathcal{O}(\tilde\beta^{-3/2})$ \cite[eq.\ (4.17)]{CGS14} only gives order $\tilde\beta^{-9/2}$ for $\tilde\beta \to \infty$. Since the previous term is already of order $\tilde\beta^{-5}$, by this approach we can not verify Dyson's bound rigorously. 
% for the simple lattice in the large-$S$ limit \cite[Eq.\ (131)]{Dyson2}
% coincides exactly with our result \eqref{eq:finaldyson}.

\section{Proof}
Our proof adapts the methods used recently in \cite{CG12, CGS14}. We use the Gibbs variational principle and the bosonic representation of the Heisenberg model in terms of spin-waves due to Holstein and Primakoff \cite{HP40}. Our trial state is a bosonic quasifree state that we have to supplement with a cutoff on the number of bosons. In our proof we will first remove the cutoff of the particle number. Thereafter we can use Wick's theorem to calculate expectation values, which enables us to bound the error terms and calculate the correction.

\medskip

 To get an upper bound, we use the Gibbs variational principle, which states that
\[f(S,\beta,\Lambda) \leq \frac{1}{\lvert \Lambda \rvert} \tr H_\Lambda \Gamma + \frac{1}{\beta\lvert\Lambda\rvert} \tr \Gamma \log \Gamma\]
for all positive trace class operators $\Gamma$ normalized to $\tr \Gamma = 1$ (i.\,e.\ states).

\medskip

%For definiteness we think of $H_\Lambda$ to have Dirichlet boundary conditions. 
Following a standard procedure, we first use the Gibbs variational principle to break up the system into smaller boxes with Dirichlet boundary conditions. We assume that $L = k(\lcal+1)$ for some integers $k$ and $\lcal$. 
%
On the set $\mathcal{C} := \{x \in \Lambda_L: x_i = n(\ell+1) \text{ for some }i=1,\ldots d \text{ and } n \in \Zbb\}$, we restrict the spins in our states to $S^3_x= -S$.
% This is an imprecise way of saying that they should be the eigenvectors with eigenvalue -S
Now $\Lambda_L \setminus \mathcal{C}$ is a union of translates of boxes $\Lambda_\ell := [1,\ell]^d \cap \Zbb^d$, on which the Hamiltonian on the restricted class of states becomes 
% -- see Figure \ref{fig:tiling}.
% \begin{figure}
%  \centering
%  \includegraphics[scale=0.55]{lattice}
%  \caption{Twodimensional view of a cubic lattice with $L = 9$, tiled with $k=3$ and $l=3$. Dirichlet boundary conditions are imposed on the lattice points located on the bold lines.}
%  \label{fig:tiling}
% \end{figure}\todo{rethink boundaries and the picture}
\begin{equation}H^\text{D}_{\Lambda_\ell} = H_{\Lambda_\ell} + \sum_{x \in \partial\Lambda_\ell} \left( S^2 + S S_x^3 \right).\label{eq:Dirichlet}\end{equation}
Here the boundary $\partial\Lambda_\ell$ consists of the points in $\Lambda_\ell$ having distance $1$ from $\mathcal{C}$. The extra summand is non-negative, and so $H^\text{D}_{\Lambda_\ell} \geq H_{\Lambda_\ell}$.
Due to this extra Dirichlet restriction on the states, the variational principle yields
\[f(S,\beta,\Lambda_L) \leq (1+\ell^{-1})^{-d} f^\text{D}(S,\beta,\Lambda_\ell), \quad \text{with } f^\text{D}(S,\beta,\Lambda_\lcal) := -\frac{1}{\beta \ell^d} \log \tr e^{-\beta H^\text{D}_{\Lambda_\lcal}}.\]
Letting $k\to\infty$, we obtain the following bound for the thermodynamic limit:
\[f(S,\beta) \leq (1+\ell^{-1})^{-d} f^\text{D}(S,\beta,\Lambda_\ell).\]
This bound holds for any integer $\lcal$, and we will later choose $\lcal = \tb^d S^2$ (or more precisely the nearest integer) to optimize the error bounds.

\medskip

The next step is to rewrite the Hamiltonian \eqref{eq:Dirichlet} through the Holstein-Primakoff mapping \eqref{eq:hpmapping}. Leaving aside for a moment the Dirichlet boundary condition, recall the Hamiltonian \eqref{eq:hphamiltonian}. 
We consider the formal Taylor expansion w.\,r.\,t.\ the small parameter $1/S$,
\begin{equation}\label{eq:formalexpansion}\begin{split}
   H_{\Lambda_\ell} & = S \bigg[\sum_{\langle x,y\rangle \subset \Lambda_\ell} \left( -a^*_x a_y - a^*_y a_x + n_x + n_y \right) \\
	     &\quad\quad\  + \frac{1}{S}\sum_{\langle x,y\rangle \subset \Lambda_\ell} \left( a^*_x \left(\frac{n_x}{4}+\frac{n_y}{4}\right) a_y + a^*_y \left(\frac{n_y}{4}+\frac{n_x}{4} \right) a_x - n_x n_y \right) + R_{\Lambda_\ell} \bigg]\\
	     & =: S \left( T_{\Lambda_\ell} + I_{\Lambda_\ell} + R_{\Lambda_\ell} \right) =: H_{0,\Lambda_\ell} + S R_{\Lambda_\ell}.
  \end{split}
\end{equation}
Here $T_{\Lambda_\ell}$ contains the terms which are formally of order unity and $I_{\Lambda_\ell}$ the terms formally of order $S^{-1}$, and $R_{\Lambda_\ell}$ is the remainder. The term $T_{\Lambda_\ell}$ is the second quantization of the discrete Laplacian. We will show that $I_{\Lambda_\ell}$ gives the interaction correction, whereas the contribution of $R_{\Lambda_\ell}$ is estimated to be of order $1/S^2$ and thus negligible.

Including the Dirichlet boundary term from \eqref{eq:Dirichlet} in the Laplace operator, i.\,e.\ setting $T^\text{D}_{\Lambda_\ell} := T_{\Lambda_\ell} + \sum_{x \in \partial\Lambda_\ell} \big( S + S_x^3 \big)$, we similarly get  $H^\text{D}_{\Lambda_\ell} = S \left( T^\text{D}_{\Lambda_\ell} + I_{\Lambda_\ell} + R_{\Lambda_\ell}\right) =: H^\text{D}_{0,\Lambda_\ell} + S R_{\Lambda_\ell}$. In particular, $I_{\Lambda_\ell}$ and $R_{\Lambda_\ell}$ are unchanged by the addition of Dirichlet boundary conditions to the Hamiltonian.

\medskip

In our proof we will obtain an upper bound on $f^\text{D}(S,\beta,\lambdal)$ using the trial state
\begin{equation}\label{eq:trialstate}
\Gamma^\text{D} = \frac{P e^{-\tb T^\text{D}_{\Lambda_\lcal}} P}{ \tr e^{-\tb T^\text{D}_{\Lambda_\lcal}} P}, \quad P = \prod_{x \in \Lambda_\lcal} \id(n_x \leq 2S), 
\end{equation}
in the Gibbs variational principle for $f^\text{D}(S,\beta,\lambdal)$. The projection $P$ ensures that we are in the subspace of bosonic Fock space where there are at most $2S$ particles per site; thus it is valid to use the bosonic formulas for the Hamiltonian equivalently to the formulas in terms of spin operators. (In  \cite{CGS14} a similar trial state was used, projecting on occupation numbers $n_x \leq 1$; projecting on $n_x \leq 2S$ has the advantage of giving exponential decay w.\,r.\,t.\ $S$ in Lemma \ref{lem:oneminusp}.)

\subsection{Bounding the Error Terms}\label{sec:errorterms}
By Gibbs' variational principle
\[f^\text{D}(S,\beta,\Lambda_\lcal) \leq \frac{1}{\lcal^d} \tr H^\text{D}_{\Lambda_\lcal} \Gamma^\text{D} + \frac{1}{\beta \lcal^d} \tr \Gamma^\text{D} \log \Gamma^\text{D}.\]
Inserting the expression \eqref{eq:trialstate} for the trial state we obtain
\begin{align*}f^\text{D}(S,\beta,\Lambda_\lcal) & \leq \frac{S}{\lcal^d} \tr \left(T^\tD_\lambdal + I_\lambdal + R_\lambdal\right) \frac{P e^{-\tb T^\tD_\lambdal}P}{\tr e^{-\tb \tdl} P} + \frac{1}{\beta \lcal^d} \tr \frac{P e^{-\tb \tdl} P}{\tr e^{-\tb \tdl} P} \log P e^{-\tb \tdl} P \\ & \quad - \frac{1}{\beta \lcal^d} \log \tr e^{-\tb \tdl} P.
\end{align*}
From \cite[(4.22)]{CGS14} we have the inequality
\[\begin{split}
   \tr P e^{-\tb \tdl} P \log P e^{-\tb \tdl} P & = \tr e^{-\tb \tdl/2} P e^{-\tb \tdl/2} \log e^{-\tb \tdl/2} P e^{-\tb \tdl/2} \\
   & \leq \tr e^{-\tb \tdl/2} P e^{-\tb \tdl/2} \log e^{-\tb \tdl} = -\tb \tr P e^{-\tb \tdl}\tdl.
  \end{split}
\]
For the expectation values and their normalization we introduce the notation
\[\langle \cdot\rangle = \tr \cdot\,\frac{e^{-\tb \tdl}}{\tr e^{-\tb \tdl}}, \quad \langle \cdot\rangle_P = \tr \cdot\, \frac{P e^{-\tb \tdl} P}{\tr P e^{-\tb \tdl}}, \quad N_P = \frac{\tr e^{-\tb \tdl}}{\tr P e^{-\tb \tdl}}.\]
Dividing by $S$, we thus obtain
\begin{align*}\frac{f^\text{D}(S,\beta,\Lambda_\lcal)}{S} & \leq \frac{1}{\lcal^d} \langle T^\tD_\lambdal + I_\lambdal \rangle_P + \frac{1}{\lcal^d} \langle R_\lambdal \rangle_P  - \frac{1}{\lcal^d} \frac{\tr \tdl e^{-\tb \tdl} P}{\tr e^{-\tb \tdl} P} - \frac{1}{\tb \lcal^d} \log \tr e^{-\tb \tdl} P\\
& =: \text{I} + \text{II} + \text{III} + \text{IV}.
\end{align*}
%%%%%%%%%%%%%%%%%%%%%% LEMMAS %%%%%%%%%%%%%%%%%%%%%

Before analysing terms I through IV, we establish some crucial lemmas. We reprove bounds on the expected number of bosons at site $x \in \lambdal$, $\rho(x) = \langle n_x\rangle = \tr n_x e^{-\tb \tdl}/\tr e^{-\tb \tdl}$ (c.\,f.\ \cite{CGS14} for another proof for $d=3$), and use them to show that $\langle 1-P\rangle$ is exponentially decaying as $S\to \infty$.

\begin{lem}\label{lem:rhobound}
The number of bosons at lattice site $x$, $\rho(x) = \langle n_x\rangle$, is bounded  by
%\[\rho(x) \leq \frac{8}{(2\pi)^{3/2}} \zeta(3/2) \frac{1}{{\tilde \beta}^{3/2}} \quad\text{for } d=3\]
%and 
\begin{equation}\label{eq:numberbound}\sup_{x \in \lambdal} \rho(x) \leq \frac{\pi^{3/2}}{8} \zeta(3/2) \tb^{-3/2} \ (d=3) \quad \text{and}\quad \sup_{x \in \lambdal}\rho(x) \leq 4\pi \tb^{-1} \log(\ell) \ (d=2).\end{equation}
(For $d=2$ we have to assume $2\tb >1 > 2\tb/(\ell+1)$. The constant $4\pi$ is a rather rough estimate.)
\end{lem}
% THERE IS NO ERROR TERM IN THIS ESTIMATE!
\begin{proof}
We use the Fourier transform of the creation and annihilation operators. It is given by
\begin{equation}
\label{eq:deffouriertransform}a^*_x = \sum_{k \in \lambdaldual} \varphi_k(x) a^*_k, \quad a_x = \sum_{k \in \lambdaldual} \varphi_k(x) a_k,\quad \lambdaldual =  \frac{\pi}{\lcal+1}\{1,2,\ldots,l\}^d,\end{equation}
where the $\varphi_k$ are the orthonormal eigenfunctions of the discrete Laplacian with Dirichlet boundary conditions on the box $\lambdal$, i.\,e.\ 
$\varphi_k(x) = 2^{d/2} (\ell+1)^{-d/2} \prod_{j=1}^d \sin\left(x_j k_j\right)$. Their normalization is such that
$\frac{2}{\ell+1}\sum_{x_i=1}^{\ell} \sin\left(x_i k_i\right) \sin\left(x_i k'_i\right) = \delta_{k_i,k'_i}$ and $\sum_{x\in \lambdal} \varphi_k(x) \varphi_{k'}(x) = \delta_{k,k'}$. Thus we find
\begin{equation}\label{eq:rhoboundcalc}\rho(x) = \tr n_x \frac{e^{-\tb \tdl}}{\tr e^{-\tb\tdl}} = \sum_{k \in \lambdaldual} \lvert \varphi_k(x)\rvert^2 \frac{1}{e^{\tilde \beta \varepsilon(k)}-1}\leq \left(\frac{2}{\pi}\right)^d \sum_{k \in\lambdaldual} \left(\frac{\pi}{\ell+1}\right)^d \frac{1}{e^{\tilde\beta \varepsilon(k)}-1}.\end{equation}

Case $d = 3$: The function $f(k):=(e^{\tb \varepsilon(k)}-1)^{-1}$ is monotonically decreasing in $k_1$, $k_2$ and $k_3$; thus  $(\frac{\pi}{\ell+1})^3\sum_{k \in\lambdaldual} f(k)$ is a lower Riemann sum and as such bounded above by the integral of $f$ on $[0,\pi]^3$. The integral is seen to be order $\tb^{-3/2}$ using  $\varepsilon(k) \geq 4 \lvert k\rvert^2/\pi^2$ and substituting $k' = \tb^{1/2} \frac{2}{\pi}k$; the numerical constant is obtained switching to spherical coordinates and using $\int_0^\infty k^2(\exp(k^2)-1)^{-1}\di k = \frac{\sqrt{\pi}}{4} \zeta(3/2)$.

Case $d=2$: Again the sum is a lower Riemann sum. However, the integral diverges at the origin, so we only use it as an upper bound outside the box $[0,\frac{\pi}{\ell+1}]^2$; inside the box we keep the original summand:
%divergent at the singularity $p \to 0$, we use the integral as an upper bound everywhere except on the box $[0,\frac{\pi}{\ell+1}]^2$, where we keep the original summand. Thus
\[\rho(x) \leq \left(\frac{2}{\pi}\right)^2 \int_{[0,\pi]^2 \backslash B_{\pi/(\ell+1)}(0)}\di^2 k \frac{1}{e^{\tb \varepsilon(k)}-1} + \left(\frac{2}{\pi}\right)^2 \frac{1}{e^{\tb \varepsilon\left((\frac{\pi}{\ell+1},\frac{\pi}{\ell+1})\right)}-1} \left(\frac{\pi}{\ell+1} \right)^2.\]
(We have actually enlarged the integral a bit by not excluding $[0,\pi/(\ell+1)]^2$ but only the ball $B_{\pi/(\ell+1)}(0)$ to simplify the further estimates.) Some basic rough estimates yield the constant.
% To conclude, we expand $\varepsilon(k)$ and the exponential around $k=0$.
%The second summand is seen to be bounded uniformly in $\ell$ for $\ell \to \infty$ by expanding $\varepsilon(p)$ and the exponential for small argument. In contrast, the integral is divergent like $\log(\ell)$ as $\ell \to \infty$. 
% 
% We obtain
% \[\rho(x) \leq C\left(\log(\ell)+1\right).\qedhere\]
\end{proof}

\begin{rem}
 For $d=3$, the estimate can be slightly improved to $\sup_{x \in \Lambda_\ell} \rho(x) \leq (2/\pi)^{3/2} \zeta(3/2) \tb^{-3/2}$, estimating the lattice heat kernel through the asymptotics of a modified Bessel function \cite[(4.15)ff]{CGS14}. However, this approach fails for $d=2$, and for this reason we stick with the estimate presented here.
\end{rem}


\begin{lem}[Exponential Decay of $\langle 1-P\rangle$]\label{lem:oneminusp}
With $C$ the respective constant from \eqref{eq:numberbound}, we have
\[\langle 1-P\rangle \leq e \ell^3 (2S+1) \left(C\tb^{-3/2}\right)^{2S} \ (d=3)\quad \text{and}\quad \langle 1-P\rangle \leq e \ell^2 (2S+1) \left( C \tb^{-1} \log (\ell) \right)^{2S}\ (d=2),\]
i.\,e.\ for large enough $\tb$ we have exponential decay as $S \to \infty$.
%, depending on the space dimension $d=2,3$, the constants are 
%$\lambda_3 = \log\left( 1+\tilde\beta^{3/2} \frac{(2\pi)^{3/2}}{16 \zeta(3/2)} \right)$ and $\lambda_2 = C/\log(\ell)$ (for some $C>0$).
\end{lem}
\begin{rems} \begin{enumerate}\item With our later choice $\ell = S^2 \tb^{d}$, this Lemma provides exponential decay of $\langle 1-P\rangle$ as $S \to \infty$. Our method fails for $d=1$ since then $\rho(x)\sim \ell$ and consequently we lose the exponential decay.
%of $\langle 1-P\rangle$.
\item For small $\tb$ and $d=3$, we find the better bound $\rho(x) \leq 8\pi \tb^{-1}$ by expanding the exponential in \eqref{eq:rhoboundcalc}.
%(There is no improvement for $d=2$.)
For the particular case $\tb = S^{-\alpha}$, this bound implies $\langle 1-P\rangle \leq C \ell^3 (2S+1) e^{-S^{1-\alpha}/4\pi}$. In this case the best pick for $\ell$ turns out to be $\ell = S^2$.
\end{enumerate}
\end{rems}
\begin{proof}[Proof of Lemma \ref{lem:oneminusp}]
Recall that $P = \prod_{x\in\lambdal} \id(n_x \leq 2S)$. Thus its expectation value is the probability that on all lattice sites $x$ we have $n_x \leq 2S$, i.\,e.\ $\langle P \rangle = \mathbb{P}(\forall x\in\lambdal :\ n_x\leq 2S)$. Consequently
\begin{align*}\langle 1-P \rangle & = \mathbb{P}(\exists x\in\lambdal :\ n_x > 2S) \leq \sum_{x\in\lambdal} \mathbb{P}(n_x > 2S) = \sum_{x\in\lambdal} \langle \id(n_x>2S) \rangle. \tagg{secondline}
  \end{align*}
We bound the step function by an exponential to see that for all $\lambda >0$ we have $\langle \id(n_x>2S) \rangle \leq \langle e^{\lambda n_x} \rangle e^{-\lambda 2S}$. 
Now denoting by $:\!\cdot\!:$ the normal-ordered product (i.\,e.\ all $a^*$s to the left of the $a$s), we have $\langle e^{\lambda n_x}\rangle = \langle :\!\exp(g(\lambda) a^*_x a_x)\!: \rangle$; this is proven taking the trace in a basis of eigenvectors of $n_x$. Expanding the normal-ordered exponential and using that by Wick's theorem $\langle a^*_x \cdots a^*_x a_x \cdots a_x\rangle = n! \langle a^*_x a_x \rangle^n$, we obtain
\[\langle e^{\lambda n_x} \rangle = \frac{1}{1-g(\lambda)\rho(x)},\]
where $g(\lambda) := e^\lambda-1$. Minimizing $\left(1-g(\lambda)\rho(x)\right)^{-1}e^{-\lambda 2S}$ w.\,r.\,t.\ $\lambda$, we easily find\footnote{For points $x \in \lambdal$ with $\rho(x) = 0$ we can't minimize, but the estimate is trivially true.}
\[\langle 1-P\rangle \leq \sum_{x\in\lambdal} \frac{2S+1}{1+\rho(x)} \left(\frac{2S+1}{2S} \frac{\rho(x)}{1+\rho(x)}\right)^{2S} \leq \sum_{x\in\lambdal} (2S+1)e \rho(x)^{2S}.\]
Finally we use Lemma \ref{lem:rhobound} to bound $\rho(x)$.
\end{proof}
As a corollary we prove that $N_P \simeq 1$, up to an error of order $\tb^{-S}$.
%
\begin{cor}\label{lem:normalizationfactor}
Let $\ell$, $S$, $\tb$ such that $\langle 1-P\rangle \leq 1/2$. Then
\[1 \leq N_P \leq 1 + 2e\ell^d (2S+1) \left(C \tb^{-d/2}\right)^{2S} (\log{\ell})^{(3-d)2S}.\]
\end{cor}
\begin{proof}
 We have $N_P = \left(1-\langle 1-P\rangle\right)^{-1} \leq 1+2\langle 1-P\rangle$ as long as $\langle 1-P\rangle \leq 1/2$.
\end{proof}

% \begin{proof}
% This is a simple consequence of using $P = 1 - (1-P)$ in the denominator of $N_P = \tr e^{-\beta \tdl}/ \tr P e^{-\beta \tdl}$, the previous lemma and the estimate $\frac{1}{1-x} \leq 1 + 2x$ for all $0 \leq x\leq 1/2$ (here $x = \langle 1-P\rangle$; the condition $\ell^{d} e^{-\lambda_d 2 S} \leq 1/4$ stems from demanding $\langle 1-P\rangle \leq 1/2$).
% \end{proof}
We are now ready to analyse terms I through IV.

\medskip
%%%%%%%%%%%%%%%%%%%%%% END OF LEMMAS %%%%%%%%%%%%%%%%%
\emph{Term I.}
%We want to remove the projection $P$ so that we get the expectation of $I_{\lambdal}$ in a quasifree state, which will be the energy correction.
We remove the projection so that we can later calculate the expectation value using Wick's theorem.
\begin{align*}
 \text{I} & = \frac{\langle H^\tD_{0,\lambdal}\rangle_P}{S \lcal^d} = \frac{N_P}{S \lcal^d}  \langle H^\tD_{0,\lambdal}\rangle- \frac{N_P}{S \lcal^d} \langle H^\tD_{0,\lambdal}(1-P) + (1-P)H^\tD_{0,\lambdal}P\rangle.
\end{align*}
Later we are going to show that the kinetic part of $\langle H^\tD_{0,\lambdal}\rangle$
% the expectation of $H^\tD_{0,\lambdal} = S(T^D_\lambdal+I_\lambdal)$
is cancelled by a contribution from term III, and only the expectation value of $I_\lambdal$ remains. Expectation values containing $(1-P)$ are small by the following lemma.
\begin{lem}\label{lem:lem24}
There exists $C > 0$ such that % C is independent of $\tilde\beta$, I checked
\begin{align*}
 &\left\lvert\frac{1}{S\lcal^d} \langle H^\tD_{0,\lambdal}(1-P) + (1-P)H^\tD_{0,\lambdal}P \rangle\right\rvert \leq C \langle 1-P \rangle^{1/2} \sup_{x\in \lambdal} (\rho(x)+1)^{3/2}\rho(x)^{1/2}.
\end{align*}
\end{lem}
\begin{proof}
 By Cauchy-Schwarz
 \begin{equation}\label{eq:oneminuspbound}
 \begin{split}
 &\left\lvert \langle H^\tD_{0,\lambdal}(1-P) + (1-P)H^\tD_{0,\lambdal}P \rangle\right\rvert \leq \langle 1-P \rangle^{1/2} \langle (H^\tD_{0,\lambdal})^2 \rangle^{1/2} +  \langle 1-P \rangle^{1/2} \langle P (H^\tD_{0,\lambdal})^2 P \rangle^{1/2}.
\end{split}
\end{equation}
 %We will bound $(H^\tD_{0,\lambdal})^2$ by powers of the operators $n_x$. Afterwards we can remove the projections $P$ in the last expectation value in \eqref{eq:oneminuspbound}, so that we will be able to use Wick's rule.
 Let us denote by $(x,y) \subset \lambdal$ all \emph{ordered} nearest neighbour pairs. Notice that the boundary term can be written as $\sum_{x \in \partial\lambdal} \big(S^2 + S  S^3_x \big) = \sum_{x \in \partial\lambdal} S n_x$. Using Cauchy-Schwarz
 \begin{equation}\label{eq:xx}\begin{split}
 ( H^\text{D}_{0,\lambdal})^2 & \leq S^2 \sum_{(x,y) \subset \lambdal} \sum_{(x',y')\subset\lambdal} \left( -a^*_x a_y + n_x \right) \left( -a^*_{x'} a_{y'} + n_{x'}\right) +  \sum_{x \in \partial\lambdal} \sum_{x' \in \partial\lambdal} S^2 n_x n_{x'}\\
			      & \quad + \sum_{(x,y) \subset \lambdal} \sum_{(x',y')\subset\lambdal} \left(a^*_x\frac{n_x+n_y}{4}a_y - \frac{n_x n_y}{2}\right) \left( a^*_{x'} \frac{n_{x'}+ n_{y'}}{4} a_{y'} - \frac{n_{x'} n_{y'}}{2} \right).
% \\			      & \leq S^2 \sum_{(x,y) \subset \lambdal} \sum_{(x',y')\subset\lambdal} \left( a^*_x a_y a^*_{x'} a_{y'} - a^*_x a_y n_{x'} - n_x a^*_{x'} a_{y'} + n_{x} n_{x'} \right)  + \sum_{x,x' \in \partial\lambdal} S^2 n_x n_{x'}\\
% 			      & \quad + \sum_{(x,y) \subset \lambdal} \sum_{(x',y')\subset\lambdal} \bigg(a^*_x\frac{n_x+n_y}{4}a_y  a^*_{x'} \frac{n_{x'}+ n_{y'}}{4} a_{y'} -  a^*_x\frac{n_x+n_y}{4}a_y {n_{x'} n_{y'}}
% 			      %\\& \hspace{4cm} + \frac{n_x n_y}{2} a^*_{x'} \frac{n_{x'}+n_{y'}}{4} a_{y'}
% 			      + n_x n_y n_{x'} n_{y'}\bigg).
 \end{split}\end{equation}
%   We will now estimate all these terms using number operators. 
   % We give the estimate for the first term; the other terms can be controlled similarly. For any $\psi$ in Fock space
%  \[
%   \begin{split}
%   & \bigg\lvert S^2 \sum_{(x,y)\subset \lambdal} \sum_{(x',y')\subset\lambdal} \langle \psi, a^*_x a_y a^*_{x'} a_{y'} \psi \rangle \bigg\rvert \leq S^2 \sum_{(x,y)\subset \lambdal} \sum_{(x',y')\subset\lambdal} \norm{a^*_y a_x \psi} \norm{a^*_{x'} a_{y'}\psi}\\
%   & \leq \frac{S^2}{2} \sum_{(x,y)\subset \lambdal} \sum_{(x',y')\subset\lambdal} \big( \langle \psi, (n_x n_y + n_x) \psi\rangle + \langle\psi, (n_{y'}n_{x'} + n_{y'})\psi\rangle \big) \leq S^2 2d \lcal^d \sum_{(x,y)\subset \lambdal} \langle \psi, (n_x n_y + n_x)\psi\rangle.
%   \end{split}
%  \]
%  Thus also
We now focus on the second summand of \eqref{eq:oneminuspbound} (the first summand is simpler) and estimate it using \eqref{eq:xx}. The contribution of the first term of \eqref{eq:xx} can by Cauchy-Schwarz be estimated as
%\footnote{Since in \eqref{eq:oneminuspbound} we have exponential decay w.\,r.\,t.\ $S \to \infty$ from $\langle 1-P\rangle$ and Lemma \ref{lem:oneminusp}, we can afford the $S^2\ell^d$ here.}
 \[
 \begin{split}
 & \langle P S^2 \sum_{(x,y)\subset \lambdal} \sum_{(x',y')\subset\lambdal} a^*_x a_y a^*_{x'} a_{y'} P \rangle  \leq 2d S^2 \lcal^d \sum_{(x,y)\subset \lambdal} \langle P  (n_x n_y + n_x) P \rangle;
 \end{split}
\]
the other contributions of \eqref{eq:xx} similarly. Since $P$ commutes with all the operators $n_x$, we can now drop it for an upper bound. Now let us extend the definition of $\rho$ to
\[\rho(x,y) := \langle a^*_y a_x \rangle = \tr a^*_y a_x \frac{e^{-\tb \tdl}}{\tr e^{-\tb \tdl}} \quad \text{(so that $\rho(x)=\rho(x,x)$)}.\]
Then Wick's theorem, followed by Cauchy-Schwarz $\lvert \rho(x,y)\rvert \leq \rho(x)^{1/2} \rho(y)^{1/2}$, yields
 \[
 \begin{split}
 2d S^2 \lcal^d \sum_{(x,y)\subset \lambdal} \langle n_x n_y + n_x \rangle 
 & \leq 2d S^2 \lcal^d \sum_{(x,y)\subset \lambdal} \left( \rho(x)\rho(y)+\rho(x,y)\rho(y,x)+ \rho(x) \right)\\
 & \leq 8d^2 S^2 \lcal^{2d} \sup_{x \in \lambdal} \left( \rho(x)+1 \right)\rho(x). \qedhere
 \end{split}
\]
% In the last step we used that due to the Cauchy-Schwarz inequality
% % \begin{equation}\label{eq:rhoxy}
% $\lvert \rho(x,y) \rvert \leq \frac{1}{2} \rho(x) + \frac{1}{2}\rho(y)$ and then bounded $\rho(x)$ and $\rho(y)$ by $\sup_{x \in \lambdal} (\rho(x)+1)$.
% , we have lost a volume factor $\lcal^d$ (due to the sum over $x \in \lambdal$; the sum over $y$ is then only over the nearest neighbours of $x$, giving an irrelevant factor $2d$). The other terms work similarly.
\end{proof}
\emph{Term II.} This is an error term of order $S^{-2}$. 
%as proven in the next lemma.
As the only error term that is not exponentially small, it constitutes the biggest error in our main theorem.
%It constitutes the biggest error in our main theorem (the only error that is not exponentially small); in particular it is responsible for the $(\log S)^3$ in $d=2$.
\begin{lem}\label{lem:II}There exists $C>0$ such that
\[ \lvert \text{II}\rvert = \left\lvert \frac{1}{\lcal^d} \langle R_\lambdal \rangle_P \right\rvert \leq \frac{N_P}{S^2} C  \sup_{x\in\lambdal} (\rho(x)+1) \rho(x)^2.\]
\end{lem}
\begin{proof}
According to \eqref{eq:formalexpansion}, on the subspace with at most $2S$ bosons per site
\[R_\lambdal = \sum_{(x,y) \subset \lambdal} a^*_x \left[\left(1-\frac{n_x}{4S}-\frac{n_y}{4S} \right)- \sqrt{1-\frac{n_x}{2S}}\sqrt{1-\frac{n_y}{2S}} \right] a_y =: \sum_{(x,y)\subset \lambdal} a^*_x A_{x,y} a_y.\]
We have again shortened the expression by writing it as a sum over all \emph{ordered} nearest-neighbour pairs $(x,y)$.

As an operator $A_{x,y} \geq 0$; to see this, notice that it depends only on $n_x$ and $n_y$, which can be diagonalized simultaneously.
% We can thus consider $A_{x,y}$ as a function of two variables (taking the place of $n_x$ and $n_y$) in $[0,2S]$. Now it is easily seen to be indeed non-negative.
 Consequently (for $\psi$ a vector with at most $2S$ particles per site)
\begin{align*}
 \lvert \langle \psi, R_\lambdal \psi \rangle\rvert & =  \bigg\lvert \sum_{(x,y) \subset \lambdal} \langle \psi, a^*_x A_{x,y}^{1/2} A_{x,y}^{1/2} a_y \psi\rangle\bigg\rvert
%  \leq S \sum_{(x,y) \subset \lambdal} \norm{A_{x,y}^{1/2} a_x \psi} \norm{A_{x,y}^{1/2} a_y \psi}\\
%  & \leq S \sum_{(x,y) \subset \lambdal} \frac{1}{2}\left(  \langle \psi, a^*_x A_{x,y} a_x \psi\rangle + \langle \psi, a^*_y A_{x,y} a_y \psi \rangle \right)
  \leq \sum_{(x,y) \subset \lambdal} \langle \psi, a^*_x A_{x,y} a_x \psi\rangle,
\end{align*}
and analogously for expectation values $\langle \cdot \rangle_P$.
%Now
% that we have symmetrically $a^*_x$ and $a_x$ on the left and right of $A_{x,y}$,
%we can employ an upper bound for $A_{x,y}$.
Again diagonalizing $n_x$ and $n_y$ simultaneously, and using that $\sqrt{1-t} \geq 1 - \frac{t}{2} - \frac{t^2}{2}$ for all $t \in [0,1]$, we find
\[A_{x,y} \leq \frac{1}{8S^2}\left(n_x^2+n_y^2\right).\]
%\todo{normalordering produces lower powers of $\rho$ which are worse for large $\tilde\beta$}
% Thus we get an expression containing only $n_x$ and $n_y$.
We then write $a^*_x (n_x^2+n_y^2) a_x = n_x(n_x-1)^2+n_x n_y^2$. Now, since $n_x$ and  $n_y$ both commute with $P$, we can drop the $P$s for an upper bound, arriving at
\[
\begin{split}
\frac{1}{\lcal^d} \langle R_\lambdal \rangle_P
% & \leq \frac{1}{\lcal^d} \tr \frac{1}{8S} \sum_{(x,y)\subset \lambdal} \left((n_x+1)^2 n_x + n_y^2 n_x\right) \frac{P e^{-\beta T^\tD_\lambdal}P}{\tr e^{-\beta \tdl} P}\\
 \leq \frac{1}{\lcal^d} \frac{1}{8S^2} \sum_{(x,y)\subset \lambdal} \tr \left( n_x(n_x-1)^2 + n_y^2 n_x\right) \frac{e^{-\tb T^\tD_\lambdal}}{\tr e^{-\tb \tdl}} \frac{\tr e^{-\tb \tdl}}{\tr e^{-\tb \tdl}P}.
\end{split}
\]
Using Wick's theorem this is expanded in terms of $\rho(x)$, $\rho(y)$, and $\rho(x,y)$, and then estimated.
% We use the Wick theorem to evaluate the expectation. The result is
% \begin{equation}\label{eq:boundonR}
% \begin{split}
% & \frac{1}{\lcal^d} \langle R_\lambdal \rangle_P \leq \frac{1}{\lcal^d} \frac{C}{S} N_P \sum_{(x,y)\subset \lambdal} \left(\rho(x)^3 + \rho(x)^2 + \rho(x) + \rho(y)^2\rho(x) + \rho(y) \rho(x,y) \rho(y,x)\right).
% \end{split}\nn
% \end{equation}
% The proof is completed like the proof of Lemma \ref{lem:lem24}.
\end{proof}
\emph{Term III.} This splits into a term which cancels the $\tdl$ of Term I, and an error term, i.\,e.\ 
\begin{align*}
 \text{III} & = - \frac{N_P}{\lcal^d}  \langle \tdl \rangle + \frac{N_P}{\lcal^d}  \langle \tdl(1-P) \rangle, \quad \left\lvert \frac{1}{\lcal^d} \langle \tdl(1-P) \rangle \right \rvert \leq  \langle 1-P\rangle^{1/2} \left(\frac{1}{\lcal^{2d}}\langle (\tdl)^2 \rangle \right)^{1/2}.
\end{align*}
%\begin{lem}
Using the momentum space creation/annihiliation operators and Wick's theorem, we find
\[\frac{1}{\ell^{2d}}\langle (\tdl)^2 \rangle = \left(\frac{1}{\ell^d} \sum_{k \in \lambdaldual} \frac{\varepsilon(k)}{e^{\tb \varepsilon(k)}-1} \right)^2 + \frac{1}{\ell^{2d}}\sum_{k\in\lambdaldual} \frac{\varepsilon(k)^2 e^{\tb \varepsilon(k)}}{\left(e^{\tb\varepsilon(k)}-1\right)^2}
\leq C \tb^{-2},\]
by (crudely) estimating $\varepsilon(k)(e^{\tb \varepsilon(k)}-1)^{-1} \leq \tb^{-1} \sup_{x>0} x(e^x-1)^{-1}$, and similarly for the second sum.
% \leq C \left(\tb^{-(d+2)}+\frac{\tb^{-(d+4)/2}}{\ell^d} \right),\]
%for our choice of $\ell$, depending on powers of $\tb$, actually the first term becomes leading, i.e. we get $\tb^{-5/2}.
% where the estimate can be verified approximating the sums by integrals and controlling the error analogously to Lemma \ref{lem:continuumapproximation} (see also \eqref{eq:result}ff).
% where the estimate can be verified approximating the sums by integrals (on top of the error for replacing a Riemann sum by an integral there is an additional lower order (by $\ell^{-1}$) error because summands with $k_j=0$ (for one or more $j=1,\ldots,d$) are missing in $\lambdaldual$).
% , which is optimal.
%\end{lem}
%\todo{is there a $log$-singularity for two dimensions?}
% \begin{proof}
%  Using the Fourier transformed creation and annihiliation operators, the formula in terms of the sums is a obtained via the Wick theorem. For the bound, notice that
%  \[\frac{1}{\lcal^d}\sum_{p\in \lambdaldual} \frac{\varepsilon(p)}{e^{\tilde \beta \varepsilon(p)}-1} \leq \int_{[-\pi,\pi]^d} \frac{\di^d k}{(2\pi)^d}  \frac{\varepsilon(k)}{e^{\tilde \beta \varepsilon(k)}-1} + \mathcal{O}(\lcal^{-1}).\]
% (On top of the error for replacing a Riemann sum by an integral there is an additional $\mathcal{O}(\lcal^{-1})$-error because summands with $k_j=0$ (for one or more $j=1,\ldots,d$) are missing.) The second summand on the r.\,h.\,s.\ likewise is of order $\ell^{d}$.
% \end{proof}

\medskip

\emph{Term IV.} We have
\[\text{IV} = -\frac{1}{\tb \lcal^d} \log \tr e^{-\tb \tdl} + \frac{1}{\tb \lcal^d} \log N_P
%\left( \frac{\tr e^{-\tb \tdl}}{\tr P e^{-\tb \tdl}} \right)
.\]

\emph{Taking parts I through IV together,} we obtain
\[\begin{split}
   \frac{f^\text{D}(S,\beta,\lambdal)}{S} \leq - \frac{1}{\tb \lcal^d} \log \tr e^{-\tb \tdl} + \frac{1}{\lcal^d} \langle I_{\lambdal} \rangle + \mathcal{E},\\
  \end{split}
\]
% \begin{equation}
% f^\text{D}(S,\beta,\Lambda_\lcal) \leq \frac{1}{\lcal^d} \langle I_\lambdal \rangle N_P - \frac{1}{\beta \lcal^d} \log \tr e^{-\beta \tdl} + \mathcal{E},\nn
% \end{equation}
 where the terms collected in the error $\mathcal{E}$ can be estimated as
\[\begin{split}
   \lvert\mathcal{E}\rvert & \leq (N_P-1) \frac{\langle I_{\lambdal} \rangle}{\ell^d}+\frac{\log N_P}{\tb \lcal^d}
   + C \langle 1-P\rangle^{1/2} \left[\sup_{x \in \lambdal} \left(\rho(x)+1\right)^{3/2} \rho(x)^{1/2}+N_P \tb^{-1} \right]\\ & \quad + \frac{C N_P}{S^2}\sup_{x\in \lambdal} (\rho(x)+1)^2 \rho(x).
%    + C N_P
%    %\left(
%    \left( \frac{(\log \ell)^{3(3-d)}}{S^2} +  \langle 1-P\rangle^{1/2}\right)\sup_{x\in \lambdal} (\rho(x)+1)^2\rho(x). %+\frac{\langle 1-P \rangle^{1/2}  \langle (\tdl)^2\rangle^{1/2}}{\lcal^d} \right).
  \end{split}\]
% 
% We arrive at the following upper bound for the free energy:
% % where
% \[\begin{split}
%    \lvert \mathcal{E}' \rvert &\leq 4 e^{-\lambda_d S} \langle I_{\lambdal} \rangle + \frac{\log\left(1+4\lcal^d e^{-\lambda_d 2S} \right)}{\beta \lcal^d} + (1+4\lcal^d e^{-\lambda_d 2 S})C\left[ S \lcal^{d/2} e^{-\lambda_d S} + \frac{(\log \ell)^{3(3-d)}}{S}
% %   + \lcal^{d/2} e^{-\lambda_d S} S
% %    \cdot \frac{1}{\lcal^d}\sum_{p\in \lambdaldual} \frac{\varepsilon(p)}{e^{\tilde \beta \varepsilon(p)}-1}
%    \right].
%   \end{split}
% \]
% The summand in the error with the logarithm can be simplified using $\log(1+x) \leq x$.
Furthermore, the integral approximation of the leading term is
\[\begin{split}
-\frac{1}{\tb \lcal^d} \log \tr e^{-\tb \tdl} & = \frac{1}{\tb \lcal^d} \sum_{k\in \lambdaldual} \log \left( 1-e^{-\tb \varepsilon(k)}\right) \leq \frac{1}{\tb}\int_{[-\pi,\pi]^d} \frac{\di^d k}{(2\pi)^d} \log \left( 1-e^{-\tilde\beta \varepsilon(k)}\right) + \frac{C}{\tb \lcal}  
  \end{split}
\]
and thus, employing Lemmas \ref{lem:rhobound} through \ref{lem:normalizationfactor}, we arrive at the following proposition:
\begin{prp}[Preliminary upper bound]\label{prp:preliminarybound}
The free energy has the upper bound
\[\frac{f^\text{D}(S,\beta,\lambdal)}{S} \leq \frac{1}{\tilde \beta} \int_{[-\pi,\pi]^d} \frac{\di^d k}{(2\pi)^d} \log \left( 1-e^{-\tilde \beta \varepsilon(k)} \right) + \frac{1}{\lcal^d} \langle I_\lambdal \rangle + \mathcal{E}',\]
where the error term $\mathcal{E}'$ satisfies
\[\begin{split}
   \lvert \mathcal{E}' \rvert &\leq \frac{C}{\tilde \beta \ell}+
   C \left[ \ell^d (2S+1) \left(C \tb^{-d/2} \right)^{2S} \log(\ell)^{(3-d)(2S+4)} \right]^{1/2} \left(\langle I_\lambdal \rangle + 1\right) + \frac{C (\log \ell)^{(3-d)3}}{\tb^d S^2}.
%    4 {e^{-\lambda_d S}} \langle I_{\lambdal} \rangle + \frac{4}{\tilde \beta} e^{-\lambda_d 2S}  + C(1+4\lcal^d e^{-\lambda_d 2 S})\left[ \frac{(\log \ell)^{3(3-d)}}{S^2} + \lcal^{d/2} e^{-\lambda_d S}\right].
% %    + \lcal^{d/2} e^{-\lambda_d S} \cdot \frac{1}{\lcal^d}\sum_{p\in \lambdaldual} \frac{\varepsilon(p)}{e^{\tilde \beta \varepsilon(p)}-1}\right].
  \end{split}
\]
\end{prp}
From the last term we see that the error in the best case will be of order $\tb^{-d}S^{-2}$ (for $d=2$ with a logarithmic correction); to make also the first term and the error from Proposition \ref{prp:twoeight} (see below) that small \emph{we need to choose $\lcal=S^2 \tb^{d}$}. The middle term is exponentially small in $S$ (provided $\tb$ is so large that $C\tb^{-d/2} < 1$, c.\,f.\ Lemma \ref{lem:oneminusp}).
% , i.\,e.\ exponentially small as $S \to \infty$ and at the same time smaller than any power law $\tb^{-m}$ ($m$ fixed).
% The size of the expectation of $I_\lambdal$ in the error estimate will become clear in the next section when we evaluate it explicitly  (providing us at the same time with the next-to-leading order for the main result).

To complete the proof of our theorem, it remains to calculate $\langle I_{\lambdal}\rangle$.

\subsection{Evaluating the Energy Correction}
Here we calculate $\lcal^{-d} \langle I_\lambdal \rangle$. With periodic boundary conditions this is formally simple, but involves infinities for $k=0$. As before, for the rigorous proof we have Dirichlet boundary conditions, making the evaluation somewhat more complicated. In particular, we will find finite-size errors (smaller by $1/\ell$ compared to the leading `bulk' term).
\begin{prp}[First order correction]\label{prp:twoeight} We have
\[\frac{1}{\ell^{d}} \langle I_\lambdal \rangle \leq - \frac{1}{S} \frac{1}{4 d} \left( \int_{[-\pi,\pi]^d} \frac{\di^d k}{(2\pi)^d} \frac{\varepsilon(k)}{e^{\tilde \beta \varepsilon(k)}-1}\right)^2 + \frac{C(\log \ell)^{3-d}}{S\ell}.\]
%In particular, its absolute value is bounded by a constant (independent of $\ell$ and $S$).
\end{prp}
% Combining this result with Proposition \ref{prp:preliminarybound}, the proof of our theorem is complete.

\begin{proof}
 Recall that according to \eqref{eq:formalexpansion} we have
 \begin{equation}\label{eq:termI}
  I_\lambdal = \frac{1}{4S}\sum_{\langle x,y\rangle \subset \lambdal} \left( a^*_x a^*_x a_x a_y + a^*_x a^*_y a_y a_y + a^*_y a^*_x a_x a_x + a^*_y a^*_y a_y a_x - 4 a^*_x a^*_y a_x a_y \right).
 \end{equation}
 By Wick's theorem
\[\begin{split} \langle I_\lambdal\rangle  = \frac{1}{S}\sum_{\langle x,y\rangle \subset \lambdal} \left( \rho(x,x)\rho(x,y) + \rho(x,x)\rho(y,x) - \rho(x,x)\rho(y,y) - \rho(x,y) \rho(y,x)\right).
  \end{split}
\]
We use the Fourier representation \eqref{eq:deffouriertransform} of $a^*_x$ and $a_x$ to calculate $\rho(y,x)$, giving
\[\rho(y,x) = \tr a^*_x a_y \gibbs = \sum_{k,k' \in \lambdaldual} \varphi_k(x) \varphi_{k'}(y) \tr a^*_k a_{k'} \gibbs = \sum_{k \in \lambdaldual} \varphi_k(x) \varphi_k(y) \frac{1}{e^{\tilde \beta \varepsilon(k)}-1}.\]
% We have the Bose-Einstein distribution $\tr a^*_p a_{p'} e^{-\beta \tdl}/\tr e^{-\beta \tdl} = \delta_{p,p'} (e^{\tilde \beta \varepsilon(p)}-1)^{-1}$. Thus
% \[\rho(y,x) \]
Thus we obtain
% \footnote{Notice that for $x_j=\ell+1$ the point $x_j + \delta_{j,j}$ is outside the box---however, in this case in the above expression there is always a factor $\sin((l+1)p_j) = 0$ going along, so contributions of that point actually vanish. We are also free to include or not include the points $x_j=0$, since their contribution is vanishing as well.}
(using the abbreviation $f(k)=(e^{\tilde \beta \varepsilon(k)}-1)^{-1}$)
\[\begin{split}
\langle I_\lambdal\rangle &= - \frac{1}{S}\sum_{\langle x,y\rangle \subset \lambdal}\sum_{k,k' \in \lambdaldual} f(k) f(k')\big[\varphi_k(x)^2 \varphi_{k'}(y)^2 - 2 \varphi_k(x)^2 \varphi_{k'}(x) \varphi_{k'}(y)  + \varphi_k(x)\varphi_k(y)\varphi_{k'}(x)\varphi_{k'}(y)\big]\\
&= - \frac{1}{S}\sum_{\substack{x_j=1,\ldots,\ell\\j=1,\ldots,d}} \sum_{i=1}^d \sum_{k,k'\in\lambdaldual} f(k) f(k')\left(\frac{2}{\ell+1}\right)^{2d}\\
&\quad\times \big[\prod_{j=1}^d \sin(x_j k_j)^2 \sin((x_j+\delta_{i,j})k'_j)^2+2\prod_{j=1}^d \sin(x_j k_j)^2 \sin(x_j k'_j) \sin((x_j+\delta_{i,j})k'_j)\\
&\quad\quad+\prod_{j=1}^d \sin(x_j k_j) \sin((x_j+\delta_{i,j})k_j)\sin(x_j k'_j)\sin((x_j+\delta_{i,j})k'_j)\big].
\end{split}\]
A little regrouping of the last expression yields
\begin{equation}\label{eq:prevstep}
\begin{split}
   \langle I_\lambdal\rangle & = - \frac{1}{S}\sum_{k,k'\in \lambdaldual} f(k)f(k') \left(\frac{2}{\ell+1}\right)^d \sum_{i=1}^d \prod_{j\in\{1,\ldots,d\}\backslash\{i\}} \left( \frac{2}{\ell+1} \sum_{x_j=1}^\ell\sin(x_j k_j)^2\sin(x_j k'_j)^2\right)\\
   &\quad \times \frac{2}{\ell+1}\sum_{x_i=1}^\ell\big[ \sin(x_i k_i)^2 \sin((x_i+1)k'_i)^2 - 2 \sin(x_i k_i)^2 \sin(x_i k'_i)\sin((x_i+1)k'_i)\\
   &\hspace{3cm}+\sin(x_i k_i)\sin((x_i+1)k_i)\sin(x_i k'_i)\sin((x_i+1)k'_i) \big].
  \end{split}
\end{equation}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% If instead of $x_i + 1$ in the previous formula we express it as $x_i+\delta$ and sum over $\delta = \pm1$, then we get a sign in using the trigonometric identity next,
% and then all the terms vanishing in the next equation actually vanish because the $\delta = +1$ and $\delta=-1$ summands cancel each other.
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Let us for the moment look only at the last factor of \eqref{eq:prevstep}, and expand it using $\sin(x_i k_i + k_i)=\sin(x_i k_i)\cos(k_i)+\cos(x_i k_i)\sin(k_i)$. Using Lemma \ref{lem:vanishing}(i) we find
\begin{align*}
& \frac{2}{\ell+1}\sum_{x_i=1}^{\ell}\big[ \sin(x_i k_i)^2 \sin((x_i+1)k'_i)^2 - 2 \sin(x_i k_i)^2 \sin(x_i k'_i)\sin((x_i+1)k'_i)\nn\\
   &\hspace{2.2cm}+\sin(x_i k_i)\sin((x_i+1)k_i)\sin(x_i k'_i)\sin((x_i+1)k'_i) \big]\nn\\
   & = \frac{2}{l+1}\sum_{x_i=1}^{\ell} \big[\sin(x_i k_i)^2 \sin(x_i k'_i)^2 \cos(p'_i)^2
%    & \hspace{2.7cm}+ 2 \sin(x_i p_i)^2 \sin(x_i p'_i) \cos(x_i p'_i) \cos(p'_i) \sin(p'_i)\label{eq:ii}\\
%    & \hspace{2.7cm}
   + \sin(x_i k_i)^2 \cos(x_i k'_i)^2 \sin(k'_i)^2\\
   & \hspace{2.7cm}- 2 \sin(x_i k_i)^2 \sin(x_i k'_i)^2 \cos(k'_i)
%    & \hspace{2.7cm}- 2 \sin(x_i p_i)^2 \sin(x_i p'_i) \cos(x_i p'_i) \sin(p'_i)\label{eq:v}\\
%    & \hspace{2.7cm}
   + \sin(x_i k_i)^2 \sin(x_i k'_i)^2 \cos(k_i) \cos(k'_i)\\
%    &\hspace{2.7cm} + \sin(x_i p_i) \cos(x_i p_i) \sin(x_i p_i)^2 \sin(p_i) \cos(p'_i)\label{eq:vii}\\
%    & \hspace{2.7cm}+ \sin(x_i p_i)^2 \sin(x_i p'_i) \cos(x_i p'_i) \cos(p_i) \sin(p'_i)\label{eq:viii}\\
   & \hspace{2.7cm}+ \sin(x_i k_i) \cos(x_i k_i) \sin(x_i k'_i)\cos(x_i k'_i)\sin(k_i)\sin(k'_i)\big].
\end{align*}
% The contribution of \eqref{eq:ii}, \eqref{eq:v}, \eqref{eq:vii} and \eqref{eq:viii} is vanishing, see . Having observed this,
Thereafter we further evaluate \eqref{eq:prevstep} using Lemma \ref{lem:quartictrigonometry}(ii). Then we use the symmetry between $k$ and $k'$ to write $-2\cos(k'_i)$ as $-\cos(k_i)-\cos(k'_i)$. After these steps we have
\[
% \label{eq:product}
\begin{split}
 \langle I_{\lambdal} \rangle & = - \frac{1}{S}\sum_{k,k' \in \lambdaldual} f(k) f(k') \left(\frac{2}{\ell+1}\right)^d\sum_{i=1}^d\prod_{j\in\{1,\ldots,d\}\backslash \{i\}} \left( \frac{1}{2} + \frac{1}{4}\left(\delta_{k_j,k'_j} + \delta_{k_j+k'_j,\pi} \right) \right)\\
  &\quad\quad\quad \times \frac{1}{2}\bigg\{ \left(1-\cos(k_i)\right)\left(1-\cos(k'_i)\right)
%   \\
%   &\quad\quad\ +
%   %\frac{1}{4}\left(\delta_{p_j,p'_j} + \delta_{p_j+p'_j,\pi} \right)\left(\cos(p'_i)^2-\sin(p'_i)^2 - 2\cos(p'_i)+\cos(p_i)\cos(p'_i)\right)
  +\delta_{k_i,k'_i}\left(\cos(k_i)^2-\cos(k_i)\right) + \delta_{k_i+k'_i,\pi}\cos(k_i)
  \bigg\}.
\end{split}
\]
Next we use symmetry among $k_1, \ldots, k_d$ to replace the sum over $i=1,\ldots,d$ by a factor of $d$. Since deltas eliminate a sum, every factor of a delta\footnote{These deltas constitute the finite-size errors mentioned before, and are a consequence of using Dirichlet boundary conditions. (They do not appear when formally using periodic boundary conditions.)} is effectively of order $1/\ell$. However, for $d=2$, some of the sums appearing are $\log \ell$-divergent at small $k$ and $k'$.
%Writing $\epsilon_d(\tb)$ for (from line to line differing) error terms satisfying $\lvert \epsilon_d(\tb)\rvert \leq C$, we find
The resulting estimate is
%\todo{Can we do without log-singularities in the 2D case?}
% \footnote{One has to be a bit careful here: Expande the product in \eqref{eq:product} and look at all but the leading term. Terms containing a factor $(1-\cos(p_1))(1-\cos(p'_1))$ have a bounded function as the summand, so we use its maximum to bound them; this also works for terms with a $\delta_{p_1+p'_1,\pi(\ell+1)}$. The other terms are singular; however an expansion for small $\lvert p\rvert$ shows that the integrals are finite.}
% Remaining are the terms containing $\delta_{p_1,p_1'}(2\cos(p_1)^2-\sin(p_1)^2-2\cos(p_1))$, e.\,g.\ 
% \[-d \sum_{p,p'\in\lambdaldual} f(p) f(p') \left(\frac{2}{\ell+1} \right)^d \frac{1}{2^{d+1}} \delta_{p_1,p_1'}(2\cos(p_1)^2-\sin(p_1)^2-2\cos(p_1)).\]
% The trigonometric expression in the numerator expands to $-2 p_1^2$ for small $p$. Expanding $f(p)$ and $f(p')$ we get $(p_1^2+p_2^2 + p_3^2)^{-1} (p^2_1+p'^2_2 + p'^2_3)^{-1}$ (or in two dimensions $(p_1^2+p_2^2)^{-1} (p^2_1+p'^2_2)^{-1}$). One of these two factors is compensated by the $p_1^2$ in the numerator, what remains is integrable.)
\begin{equation}\label{eq:result}
\begin{split}
\frac{1}{\ell^{d}}\langle I_{\lambdal} \rangle & \leq - \frac{d}{S\ell^{d}} \sum_{k,k' \in \lambdaldual} f(k) f(k') \left(\frac{2}{\ell+1}\right)^d \frac{1}{2^d}\left(1-\cos(k_1)\right)\left(1-\cos(k'_1)\right) + C S^{-1}(\log \ell)^{3-d}\ell^{-1}\\
&= -\frac{d}{S} \Bigg[\frac{1}{(\ell+1)^d}\sum_{k \in \lambdaldual} f(k) \left(1-\cos(k_1) \right) \Bigg]^2  + C S^{-1}(\log \ell)^{3-d}\ell^{-1}.
\end{split}
\end{equation}
Again we use the symmetry among $k_1,\ldots,k_d$, now to replace $(1-\cos(k_1))$ by $\frac{1}{2d}\varepsilon(k)$.% in \eqref{eq:result}.

It remains to employ the continuum approximation for $g(k) := f(k)\varepsilon(k)$: Lemma \ref{lem:continuumapproximation} with $D_1 = \tb^{-1} \sup_{x\in\Rbb} x(e^x-1)^{-1}$ and $D_2 = \sup_{k\in \lambdaldual} \lvert \nabla g(k)\rvert \leq 2\sqrt{3} \sup_{x\in\Rbb} \lvert x e^x +1 - e^x\rvert(e^x-1)^{-2}$ implies 
\[\frac{1}{\ell^{d}} \langle I_{\lambdal}\rangle \leq -\frac{1}{4 d S}\left[ \int_{[0,\pi]^d} \frac{\di^d k}{\pi^d} \frac{\varepsilon(k)}{e^{\tb \varepsilon(k)}-1} \right]^2 + CS^{-1} (\log \ell)^{3-d}\ell^{-1}.\]
Since $\varepsilon(k)$ depends on $k$ only through $\cos(k_i)$, the integral remains unchanged by any reflection $k_i \mapsto -k_i$.  Thus the integral over $[0,\pi]^d$ is the same as $2^{-d}$ times the integral over $[-\pi,\pi]^d$.
\end{proof}

The next two lemmas were used in the previous proof.
\begin{lem}\label{lem:vanishing}\label{lem:quartictrigonometry}
Let $k, k' \in \lambdaldual$ and $i=1,\ldots d$.
\begin{enumerate}
 \item  We have
\[\sum_{x_i=1}^{\ell} \sin(x_i k_i)^2 \sin(x_i k'_i) \cos(x_i k'_i) = 0.\]
\item Furthermore
 \begin{align*}& \frac{2}{\ell+1}\sum_{x_i=1}^{\ell} \sin(x_i k_i)^2 \sin(x_i k'_i)^2 = \frac{1}{2}+\frac{1}{4}\left( \delta_{k_i,k'_i} + \delta_{k_i+k'_i,\pi}\right),\\
 & \frac{2}{\ell+1}\sum_{x_i=1}^{\ell} \sin(x_i k_i)^2 \cos(x_i k'_i)^2 = \frac{1}{2}-\frac{1}{4}\left( \delta_{k_i,k'_i} + \delta_{k_i+k'_i,\pi} \right),\\
 & \frac{2}{\ell+1}\sum_{x_i=1}^{\ell} \sin(x_i k_i) \cos(x_i k_i) \sin(x_i k'_i) \cos(x_i k'_i) = \frac{1}{4}\left( \delta_{k_i,k'_i} - \delta_{k_i+k'_i,\pi} \right).
 \end{align*}
\end{enumerate}
\end{lem}
\begin{proof}
(i) Expanding the trigonometric functions in terms of exponentials we find
\begin{align*}
&\sum_{x_i=1}^\ell \sin(x_i k_i)^2 \sin(x_i k'_i) \cos(x_i k'_i) = \sum_{x_i=1}^\ell \left[\frac{\sin(2x_i(k_i-k'_i))-\sin(2 x_i(k_i+k'_i))}{4} + \frac{\sin(2 x_i k'_i)}{2} \right].% \tagg{wavefunctions}
\end{align*}
Due to the factor of $2$ in all the arguments on the r.\,h.\,s., the summation range  contains only multiples of full phases of the sine (recall that $k_i = \frac{\pi}{\ell+1}n$ for some $n \in \Zbb$), so negative and positive parts cancel. (ii) We simply expand into exponentials and use the finite geometric sum.
\end{proof}

%
\begin{lem}\label{lem:continuumapproximation}
 Let $g: (0,\pi]^n \to \Rbb$ be bounded above by some $D_1 < \infty$
 %. Assume that there exist $0 < a < b < \pi$ such that $g(p)$ on $[0,b]^n$ is growing w.\,r.\,t.\ all the $p_i$,
 and have Lipschitz constant $D_2 < \infty$.
 %on $[0,\pi]^n \backslash [0,a]^n$.
 Then there exists a constant $C > 0$ (depending only on the dimension $n$) such that
 \[\left(\frac{\pi}{\ell+1}\right)^n \sum_{k \in \lambdaldual} g(k) \geq \int_{[0,\pi]^n} g(k) \di^n k - C (D_1+D_2) \frac{1}{\ell+1} \quad \forall \ell \in \Nbb.\]
\end{lem}
%\todo{$\tilde\beta$-dependence?}
\begin{proof} Due to Lipschitz continuity, for every $k_0 \in \lambdaldual$ we have
 \[\left( \frac{\pi}{\ell+1} \right)^n \left[ g(k_0) + D_2\sqrt{n}\frac{\pi}{\ell+1} \right] \geq \max_{k \in k_0+[0,\pi/(\ell+1)]^n} g(k) \left( \frac{\pi}{\ell+1} \right)^n \geq \int_{k_0+[0,\pi/(\ell+1)]^n} g(k)\di^n k.\]
 Thus (the factor $\ell^n$ in the numerator of the last summand being the number of boxes in the partition, or equivalently, the number of elements of $\lambdaldual$)
 \[\left(\frac{\pi}{\ell+1} \right)^n \sum_{k_0 \in \lambdaldual} g(k_0) \geq \int_{[\pi/(\ell+1),\pi]^n} g(k)\di^n k -  \frac{\pi^n\ell^n}{(\ell+1)^n} D_2 \sqrt{n} \frac{\pi}{\ell+1}.\]
 Since $g$ is bounded above, by extending the integration range from $[\pi/(\ell+1),\pi]^n$ to $[0,\pi]^n$, we make the integral larger by a quantity of at most $C D_1/(\ell+1)$ for some $C < \infty$ depending only on $n$.
\end{proof}
% \begin{proof}
% Choose $\ell$ large enough such that there exists an $m \in \Nbb$ with $a < m\pi/(\ell+1) < b$.
%  Due to Lipschitz continuity, for every $p_0 \in \lambdaldual \backslash [0,m\pi/(\ell+1)]^n$ we have
%  \[\left( \frac{\pi}{\ell+1} \right)^n \left[ g(p_0) + D_2\sqrt{n}\frac{\pi}{\ell+1} \right] \geq \max_{p \in p_0-[0,\pi/(\ell+1)]^n} g(p) \left( \frac{\pi}{\ell+1} \right)^n \geq \int_{p_0-[0,\pi/(\ell+1)]^n} g(p)\di^n p.\]
%  Thus (the factor $\ell^n$ in the numerator of the last summand being the number of boxes in the partition, or equivalently, the number of elements of $\lambdaldual$)
%  \[\left(\frac{\pi}{\ell+1} \right)^n \sum_{p \in \lambdaldual \backslash [0,m\pi/(\ell+1)]^n} g(p) \geq \int_{[0,\ell\pi/(\ell+1)]^n\backslash [0,m\pi/(\ell+1)]^n} g(p)\di^n p -  \frac{\pi^n\ell^n}{(\ell+1)^n} D_2 \sqrt{n} \frac{\pi}{\ell+1}.\]
%  Now, since $g$ is growing on $[0,m\pi/(\ell+1)]^n$ the sum there is actually an upper Riemann sum to $g$, and the integral does a lower bound to the sum. Therefore
%  \[\left(\frac{\pi}{\ell+1} \right)^n \sum_{p \in \lambdaldual} g(p) \geq \int_{[0,\ell\pi/(\ell+1)]^n} g(p)\di^n p - \frac{\pi^n\ell^n}{(\ell+1)^n} D_2 \sqrt{n} \frac{\pi}{\ell+1}.\]
%  Finally, since $g$ is bounded above, by extending the integration range from $[0,\ell\pi/(\ell+1)]^n$ to $[0,\pi]^n$, we make the integral larger by a quantity at most of order $D_1/(\ell+1)$.
% \end{proof}

% //////////////////////
% 
% \begin{lem}[Bound on the expectation of $1-P$ in 2D]
% There exists a $\tilde \lambda_d > 0$ such that
% \[\langle 1-P \rangle \leq 2 \lcal^2 e^{-\tilde\lambda_d S/\log(\ell)}.\]
% \end{lem}
% \begin{proof}
% The proof is the same as in the threedimensional case, just using Lemma \ref{lem:rhoboundtwod} instead of Lemma \ref{lem:rhobound}: the $\lambda_d$ from Lemma \ref{lem:oneminusp} becomes\footnote{The inequality is valid for e.\,g.\ $\ell \geq 2$.} $\lambda_d = \log(1+C/\log(\ell)).$
% For $\ell \to \infty$ this diverges like $1/\log(\ell)$; more precisely we can factor out the $1/\log(\ell)$ by introducing $\tilde \lambda_d$ through
% \[\lambda_d = \log(1+C/\log(\ell)) \geq C/2\log(\ell) =: \tilde \lambda_d/\log(\ell).\qedhere \]
% \end{proof}

\appendix

\section{Appendix: Cancellation at Second Order}
\label{sec:appendix} We now consider the three dimensional case only. 
As explained in the beginning, the dependence of our error bound on $\tb$ (order $\tb^{-3}$) is not in agreement with Dyson's paper, which claims that all corrections are $\tb^{-5}$ and smaller. However, it was pointed out \cite{Oguchi} that at second order of formal perturbation theory there is a cancellation, by which Dyson's result is reproduced to order $1/S^2$ (if one corrects for a trivial numerical imprecision \cite{PD66}).

\medskip

Below we reproduce the calculation of \cite{Oguchi} in the language of modern perturbation theory and in detail. For simplicity we work in periodic boundary conditions, i.\,e.\ the eigenfunctions of the Laplacian are $\ell^{-3/2} e^{i k\cdot x}$ with momenta\footnote{All $k$-sums in this appendix are over this range.} $k \in \frac{2\pi}{\ell}\left(\Zbb^{3} \cap (-\ell,\ell)^3\right)$. We expand the Hamiltonian one order further,
\[H_{\Lambda_\ell} =: S\left( T_{\Lambda_\ell} + I_{\Lambda_\ell} + J_{\lambdal} + \tilde{R}_{\Lambda_\ell} \right), \quad \text{i.\,e.}\quad R_\lambdal = J_{\lambdal} + \tilde{R}_{\Lambda_\ell}.\]
Here $J_\lambdal$ is formally of order $S^{-2}$ and after normal-ordering given by
\[J_{\lambdal} = \frac{1}{32S^2} \sum_{(x,y) \subset \lambdal} \left( a^*_x a^*_y a^*_y a_y a_y a_y + a^*_x a^*_y a_y a_y - 2 a^*_x a^*_x a^*_y a_x a_y a_y + a^*_x a^*_x a^*_x a_x a_x a_y + a^*_x a^*_x a_x a_y \right).\]
% has been double checked
Using Wick's theorem we find the correction to $f(S,\beta,\lambdal)/S$ to be
\begin{equation}\label{eq:firstorderpert}\frac{\langle J_\lambdal \rangle}{\ell^3} = \frac{1}{4 S^2}\sum_{i=1}^3 \bigg[ (\rho+\rho^2) \frac{1}{\ell^3} \sum_{k} f(k) \cos(k_i) - \frac{1}{\ell^9} \sum_{k_1,k_2,k_3} f(k_1) f(k_2) f(k_3)\cos(k_{1,i}-k_{2,i}+k_{3,i})\bigg],\end{equation}
%has beeen double checked
where $\rho = \frac{1}{\ell^3}\sum_{k} f(k)$ and $f(k) = 1/(e^{\tb \varepsilon(k)}-1)$. The biggest contribution\footnote{Strictly speaking this sum is infinite because the contributions of $k_1=0$ and $k_2=0$ are infinite as a remnant of using periodic boundary conditions. The cancellation below also resolves this issue.} is
\begin{equation}\label{eq:biggesterror}\frac{\rho}{4 S^2} \frac{1}{\ell^3} \sum_{k} f(k) \sum_{i=1}^3 \cos(k_i) = \frac{1}{16 S^2 \ell^9} \sum_{k_1,k_2,k_3} f(k_1) f(k_2) \left[12-\varepsilon(k_1)-\varepsilon(k_2) \right] \sim \tb^{-3}.\end{equation}
The remaining part of \eqref{eq:firstorderpert} is finite and of order $\tb^{-11/2}$, as can be seen by expanding the cosines, observing that the lowest terms cancel, replacing the sum by an integral for $\ell \to \infty$ and using the scaling $\int \di^3 k\, ( \exp(\tb \lvert k\rvert^2) - 1 )^{-1} \sim \tb^{-3/2}$. % (by substituting $k' = \tb^{1/2}k$).

The big error \eqref{eq:biggesterror} originates from the summands in $J_{\lambdal}$ with four creation and annihilation operators (which all originate from normal-ordering of the formal expansion). Thus as an operator bound on $R_\lambdal$ we can not expect an estimate better than $\tb^{-3}$. Instead we have to take into account the structure of the interacting Gibbs state. We verify this by showing that in second order of perturbation theory, \eqref{eq:biggesterror} is cancelled up to a $\tb^{-5}$-remainder.

\medskip

\emph{Second order perturbation theory.} The second order perturbation theory is given through the Duhamel formula as
\begin{align}
\frac{f(S,\beta,\Lambda_\ell)}{S}  = - \frac{\log \tr e^{-\beta H_\lambdal}}{\tb \ell^3} & = - \frac{\log \tr e^{-\tb T_\lambdal}}{\tb \ell^3} + \frac{1}{\ell^3} \langle I_\lambdal+J_\lambdal \rangle\nn\\
& \quad - \frac{1}{\tb \ell^3}\int_0^{\tb}\di s \int_0^{s} \di s' \left( \langle I_\lambdal(s) I_\lambdal(s')\rangle -  \langle I_\lambdal\rangle \langle I_\lambdal\rangle \right) + \order(S^{-3}),\label{eq:secondorder} \end{align}
where $I_\lambdal(s) := e^{sT_\lambdal}I_\lambdal e^{-sT_\lambdal}$.
The contractions in momentum space (represented in the Feynman diagrams by lines with arrow pointing from the creation to the annihilation operator) are
\begin{equation}\label{eq:contractions}\begin{split}\langle a^*_k(s)a_{k'}(s')\rangle = \delta_{k,k'} e^{(s-s')\varepsilon(k)}f(k)\quad \text{and}\quad\langle a_k(s)a^*_{k'}(s')\rangle = \delta_{k,k'} e^{-(s-s')\varepsilon(k)}(1+f(k)).\end{split}\end{equation}
Using $a^*_x = \ell^{-3/2} \sum_k e^{ik\cdot x} a^*_k$, we find the interaction operator $I_{\lambdal}$ in momentum space to be
 \[\begin{split}I_\lambdal = 
%  \frac{1}{4}\sum_{\langle x,y\rangle \subset \lambdal} \left(a^*_x-a^*_y\right)^2 a_x a_y + \hc = 
\frac{1}{8 S \ell^3} \!\sum_{k_1,k_2,k_3,k_4}\! &  \delta_{k_1+k_2,k_3+k_4} \nu(k_1,k_2,k_3,k_4)a^*_{k_1} a^*_{k_2} a_{k_3} a_{k_4},\end{split}\]
where\footnote{We abbreviate $\varepsilon(k_1) = \varepsilon_1$, $\varepsilon(k_1-k_2) = \varepsilon_{1-2}$ etc.} $\nu(k_1,k_2,k_3,k_4) := \varepsilon_{4-2}-\varepsilon_4 - \varepsilon_1 + \varepsilon_{4-1} - \varepsilon_2 + \varepsilon_{3-2} + \varepsilon_{3-1} -\varepsilon_3$. Notice that we have symmetrized $\nu$ w.\,r.\,t.\ exchange of $k_1$ with $k_2$ or $k_3$ with $k_4$, so in Feynman diagrams $I_{\lambdal}$ is represented as a vertex with two equivalent outgoing and two equivalent ingoing legs.

The second order contribution \eqref{eq:secondorder} can through Wick's theorem be represented by the sum of all connected Feynman diagrams having two $I_{\lambdal}$-vertices, where the numerical factors count the number of equivalent diagrams:
\vspace{-.5em}
\begin{equation*}\label{eq:diagrams}4\begin{gathered}
\begin{fmfgraph*}(90,90)
\fmfleft{i}
\fmfright{o}
\fmf{phantom,tension=10}{i,v1}
\fmf{phantom,tension=10}{v2,o}
\fmflabel{$s$}{v1}
\fmflabel{$s'$}{v2}
\fmf{fermion,left,tension=0.4,label=$k_1$}{v1,v2}
\fmf{fermion,left,tension=0.4,label=$k_4$}{v2,v1}
\fmf{fermion,left=0.3,label=$k_2$}{v1,v2}
\fmf{fermion,left=0.3,label=$k_3$}{v2,v1}
% \fmf{fermion,right=0.5}{v1,v2}
\fmfdot{v1}
\fmfdot{v2}
\end{fmfgraph*}
\end{gathered}
+\ 16\ \hspace{2em} \begin{gathered}
\begin{fmfgraph*}(90,90)
\fmftop{t0,t1,t2,t3}
\fmfbottom{b0,b1,b2,b3}
\fmf{phantom}{t1,v1,b1}
\fmf{phantom}{t2,v2,b2}
\fmffreeze
\fmflabel{$s$}{v1}
\fmflabel{$s'$}{v2}
\fmf{fermion,right,label=$k_2$}{v1,v2}
\fmf{fermion,right,label=$k_2$}{v2,v1}
\fmf{plain,tension=0.7,right=270,label=$k_1$}{v1,v1}
\fmf{plain,tension=0.7,right,label=$k_3$}{v2,v2}
\fmfdot{v1,v2}
 \end{fmfgraph*}
\end{gathered}\hspace{2.3em}.\vspace{-.8em}\end{equation*}
\emph{Left Feynman diagram.} Using the contractions \eqref{eq:contractions} we write out the left Feynman diagram
\[\begin{split}
& -\frac{4}{\tb\ell^3} \int_0^\tb\di s\int_0^s \di s'\frac{1}{\left(8S\ell^3\right)^2}
\sum_{\substack{k_1,k_2,k_3,k_4\\k_1',k_2',k_3',k_4'}} \delta_{k_1+k_2,k_3+k_4} \delta_{k'_1+k'_2,k'_3+k'_4}
\nu(k_1,k_2,k_3,k_4)\nu(k'_1,k'_2,k'_3,k'_4)\\
&\quad\times\delta_{k_1,k_4'}\delta_{k_2,k_3'}\delta_{k_3,k_2'}\delta_{k_4,k_1'}e^{(s-s')\varepsilon_1}e^{(s-s')\varepsilon_2}e^{-(s-s')\varepsilon_3} e^{-(s-s')\varepsilon_4} f(k_1) f(k_2) (1+f(k_3)) (1+f(k_4)).
\end{split}\]
The integrations over $s'$ and $s$ are done explicitly, yielding
\[\begin{split}
& -\frac{1}{16\tb S^2 \ell^9} \sum_{k_1,k_2,k_3,k_4}\delta_{k_1+k_2,k_3+k_4}\nu(k_1,k_2,k_3,k_4)^2\\
&\quad\times\left[-\frac{\tb}{\varepsilon_1+\varepsilon_2-\varepsilon_3-\varepsilon_4}+ \frac{e^{\tb(\varepsilon_1+\varepsilon_2-\varepsilon_3-\varepsilon_4)}-1}{(\varepsilon_1+\varepsilon_2-\varepsilon_3-\varepsilon_4)^2} \right] \frac{e^{-\tb\varepsilon_1}}{1-e^{-\tb\varepsilon_1}} \frac{e^{-\tb\varepsilon_2}}{1-e^{-\tb\varepsilon_2}} \frac{1}{1-e^{-\tb\varepsilon_3}} \frac{1}{1-e^{-\tb\varepsilon_4}}.\end{split}\]
We observe that
$[e^{\tb(\varepsilon_1+\varepsilon_2-\varepsilon_3-\varepsilon_4)}-1]e^{-\tb\varepsilon_1}e^{-\tb\varepsilon_2}$ has odd sign under simultaneous exchange of $k_1$ with $k_3$ and $k_2$ with $k_4$, and is multiplied with an expression of even sign, thus summing to zero. We are left with the contribution of the first summand from the square brackets,
\begin{equation}\label{eq:leftdiagram}\begin{split}
& \frac{1}{16 S^2 \ell^9} \sum_{k_1,k_2,k_3,k_4}\delta_{k_1+k_2,k_3+k_4}\frac{\nu(k_1,k_2,k_3,k_4)^2}{\varepsilon_1+\varepsilon_2-\varepsilon_3-\varepsilon_4}f(k_1)f(k_2)(1+f(k_3))(1+f(k_4))\\
& = \frac{1}{16 S^2 \ell^9}\sum_{k_1,k_2,k_3} \frac{\left(2\varepsilon_{1-3}+2\varepsilon_{2-3}-\varepsilon_1-\varepsilon_2-\varepsilon_3-\varepsilon_{1+2-3} \right)^2}{\varepsilon_1+\varepsilon_2-\varepsilon_3-\varepsilon_{1+2-3}} \left[f(k_1)f(k_2) + 2 f(k_1) f(k_2) f(k_3) \right].\end{split}\end{equation}
(In the last line the contribution $f(k_1)f(k_2)f(k_3)f(k_4)$ disappeared by an antisymmetry argument.)

\emph{Right Feynman diagram.} The right Feynman diagram is evaluated similarly, and found to be
\begin{equation}\label{eq:rightdiagram}\begin{split}- \frac{\tb}{2 S^2 \ell^9}\sum_{k_1,k_2,k_3} \left(\varepsilon_{2-3} - \varepsilon_3-\varepsilon_2 \right) \left(\varepsilon_{2-1} - \varepsilon_1-\varepsilon_2 \right)
%\\&\quad\times
f(k_1) f(k_2) \left(1+f(k_2) \right) f(k_3).
\end{split}\end{equation} 

 Replacing the sum by an integral as $\ell \to \infty$ and expanding $\varepsilon(k) \simeq \lvert k\rvert^2$, we find that \eqref{eq:rightdiagram} is finite and of order $\tb^{-11/2}$ (so for us negligible). By the same argument, the part of \eqref{eq:leftdiagram} containing $f(k_1)f(k_2)f(k_3)$ is of order $\tb^{-11/2}$ (and thus also negligible).
 
 \medskip
 
 \emph{The cancellation.} As possibly big contributions remain the part of \eqref{eq:leftdiagram} containing only $f(k_1)f(k_2)$ and \eqref{eq:biggesterror}. The crucial observation \cite{Oguchi} is that these two terms cancel up to terms of order $\tb^{-5}$! Let us spell out the argument more carefully.
The sum of the two terms is
 \[\begin{split}& \frac{1}{16 S^2 \ell^9} \sum_{k_1,k_2,k_3} f(k_1)f(k_2) \left[ 12-\varepsilon_1-\varepsilon_2 + \frac{\left[ (\varepsilon_1+\varepsilon_2-\varepsilon_3-\varepsilon_{1+2-3})+2(\varepsilon_3 +  \varepsilon_{1+2-3} - \varepsilon_{2-3}-\varepsilon_{1-3}) \right]^2}{\varepsilon_1+\varepsilon_2 - \varepsilon_3 - \varepsilon_{1+2-3}} \right]\\
 & = \frac{1}{16 S^2 \ell^9} \sum_{k_1,k_2,k_3} f(k_1)f(k_2) \left[ 12 + 3\varepsilon_3 + 3\varepsilon_{1+2-3} - 4\varepsilon_{2-3}- 4\varepsilon_{1-3} + 4 \frac{(\varepsilon_3 +  \varepsilon_{1+2-3} - \varepsilon_{2-3}-\varepsilon_{1-3})^2}{\varepsilon_1+\varepsilon_2 - \varepsilon_3 - \varepsilon_{1+2-3}}\right]
   \end{split}
\]
 Recalling that $\varepsilon(k) = \sum_{j=1}^3 2(1-\cos(k_j))$, the first part vanishes since the sum
 is over all $k$, i.\,e.\ positive and negative parts of the cosine sum to zero:
 %is over the whole Brillouin zone:
 \[\begin{split}& \sum_{k_3} \left[ 12 + 3\varepsilon_3 + 3\varepsilon_{1+2-3} - 4\varepsilon_{2-3}- 4\varepsilon_{1-3} \right] =0.
 %\\ & = - \sum_{k_3} \sum_{i=1}^3 \left[ 6 \cos(k_{3,i}) + 6 \cos(k_{1,i} + k_{2,i} - k_{3,i}) + 8 \cos(k_{2,i}-k_{3,i}) + 8 \cos(k_{1,i}-k_{3,i}) \right] = 0.
 \end{split}\] 
 Concerning the second part, by expanding $\varepsilon(k)$, we find
 \[\begin{split}& \frac{1}{16 S^2 \ell^9} \sum_{k_1,k_2,k_3} f(k_1)f(k_2)\,4 \frac{(\varepsilon_3 +  \varepsilon_{1+2-3} - \varepsilon_{2-3}-\varepsilon_{1-3}))^2}{\varepsilon_1+\varepsilon_2 - \varepsilon_3 - \varepsilon_{1+2-3}}\\&\simeq
 \frac{1}{4S^2 \ell^9} \sum_{k_1,k_2,k_3} f(k_1)f(k_2)\frac{(k_1\cdot k_2)^2}{(k_1-k_3)\cdot (k_2-k_3)} \sim \tb^{-5} S^{-2}.\end{split}\]
 This agrees with Dyson's result at order $1/S^2$. 

Presumably such cancellations appear at all orders in perturbation theory. In Dyson's work the problematic quartic  terms in $J_{\lambdal}$ from normal-ordering do not appear, at the cost of working with a non-selfadjoint Hamiltonian. While this approach is supposed to be equivalent \cite{OtherOguchi} to the Holstein-Primakoff approach followed here, it  has never been made completely rigorous. It remains an interesting problem to rigorously obtain an upper bound with optimal $\tb$-dependence.

\section*{Acknowledgements}
I would like to thank Alessandro Giuliani and Michele Correggi for many discussions and the introduction to the topic. I acknowledge support by ERC grant CoMBoS-239694, ERC Advanced grant 321029, by MIUR through the FIR grant 2013 ``Condensed Matter in Mathematical Physics (Cond-Math)'' (code RBFR13WAET), and by VILLUM FONDEN via the QMATH Centre of Excellence (Grant No. 10059).

\bibliographystyle{plain}
\bibliography{heisenberg}
\end{fmffile}
\end{document}

macros.tex

% Set enumerate environment using roman counters
\renewcommand{\theenumi}{{\rm (\roman{enumi})}}
\renewcommand{\labelenumi}{\theenumi}

% Nicer marginpar with smaller font
\let\oldmarginpar\marginpar
\renewcommand\marginpar[1]{\-\oldmarginpar[\raggedleft\footnotesize #1]%
  {\raggedright\footnotesize #1}}

%%%%%%%%%%%%%%%%%%%%%%%%% Theorems etc %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%\newtheorem{theorem}{Theorem}[section]
\newtheorem{thm}{Theorem}[section]  % use thm for Theorems to keep numbering consistent
\newtheorem*{thm*}{Theorem}
\newtheorem{cor}[thm]{Corollary}
\newtheorem{prp}[thm]{Proposition}
\newtheorem{lem}[thm]{Lemma}
\newtheorem{ex}[thm]{Example}
\newtheorem{dfn}[thm]{Definition}
\newtheorem*{rem}{Remark}
\newtheorem*{rems}{Remarks}
\newtheorem*{hyp}{Hypothesis}
\renewcommand{\theequation}{\thesection.\arabic{equation}}

%%%%%%%%%%%%%%%%%%%%%%%% Niels' Notation %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\newcommand{\tb}{{\tilde\beta}}
\newcommand{\order}{\mathcal{O}}



\newcommand{\slater}{f_1 \wedge \ldots \wedge f_N}
\newcommand{\Scal}{\mathcal{S}}
\newenvironment{provisional}{\vspace{1em}\hrule\footnotesize\vspace{1em}}{\vspace{1em}\hrule\vspace{1em}}

\newcommand{\ad}{\operatorname{ad}}	% abbreviation for iterated commutators
\newcommand{\fock}{\mathcal{F}}		% fock space symbol
\newcommand{\di}{{\textrm{d}}}		% differential (for integrals)
\newcommand{\Lcal}{\mathcal{L}}		% calligraphic L
\newcommand{\Ncal}{\mathcal{N}}		% calligraphic N
\newcommand{\Kcal}{\mathcal{K}}		% calligraphic N
\newcommand{\Vcal}{\mathcal{V}}		% calligraphic V
\newcommand{\Hcal}{\mathcal{H}}		% calligraphic H
\newcommand{\Ocal}{\mathcal{O}}		% big-O, order-of
\newcommand{\lcal}{\ell}
\newcommand{\Zcal}{\mathcal{Z}}
\newcommand{\hc}{\mbox{h.c.}}		%hermitian conjugate
\newcommand{\scal}[2]{\big<#1,#2\big>} % scalar product
\newcommand{\cc}[1]{\overline{#1}}	% complex conjugate
\newcommand{\Rbb}{\mathbb{R}}		% real numbers
\newcommand{\Cbb}{\mathbb{C}}		% complex numbers
\newcommand{\Nbb}{\mathbb{N}}		% natural numbers
\newcommand{\Zbb}{\mathbb{Z}}
\renewcommand{\Re}{\operatorname{Re}\,} 	%RealPart
\renewcommand{\Im}{\operatorname{Im}\,} 	%ImaginaryPart
\newcommand{\id}{\mathds{1}}
\newcommand{\tD}{\text{D}}
\newcommand{\lambdal}{{\Lambda_\lcal}}
\newcommand{\lambdaldual}{\Lambda_\lcal^{*}}
\newcommand{\tdl}{T^\tD_\lambdal}
\newcommand{\norm}[1]{\lVert#1\rVert}	%Norm
\newcommand{\proj}[1]{\lvert #1 \big>\big< #1\rvert}	% orthogonal projection operator
\newcommand{\Tr}{\operatorname{tr}}	% Trace
\newcommand{\tr}{\operatorname{tr}}
\newcommand{\HS}{_{\textrm{HS}}}
\newcommand{\sgn}{\operatorname{sgn}}
\newcommand{\tagg}[1]{ \stepcounter{equation} \tag{\theequation} \label{eq:#1} } % add tag and label in align*-environments

 \newcommand{\todo}[1]{\fbox{\footnotesize #1}\marginpar{TODO}}
% \newcommand{\todo}[1]{}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Benjamin's notation %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\newcommand{\tr}{\mbox{Tr}}
\newcommand{\ov}{\overline}
\newcommand{\const}{\mbox{const}}
\newcommand{\La}{\Lambda}
\newcommand{\e}{\varepsilon}
\newcommand{\eps}{\varepsilon}
\newcommand{\pt}{\partial}
\newcommand{\rd}{{\rm d}}
\newcommand{\bR}{{\mathbb R}}
\newcommand{\bbZ}{{\mathbb Z}}
\newcommand{\bke}[1]{\left( #1 \right)}
\newcommand{\bkt}[1]{\left[ #1 \right]}
\newcommand{\bket}[1]{\left\{ #1 \right\}}
%\newcommand{\norm}[1]{\| #1 \|}
\newcommand{\Norm}[1]{\left\Vert #1 \right\Vert}
\newcommand{\abs}[1]{\left| #1 \right|}
\newcommand{\bka}[1]{\left\langle #1 \right\rangle}
\newcommand{\vect}[1]{\begin{bmatrix} #1 \end{bmatrix}}
\newcommand{\bs}{{\bf{s}}}
\newcommand{\ba}{{\bf{a}}}
\newcommand{\bb}{{\bf{b}}}
\newcommand{\bx}{{\bf{x}}}
\newcommand{\by}{{\bf{y}}}
\newcommand{\bt}{{\bf{t}}}
\newcommand{\bu}{{\bf{u}}}
\newcommand{\bv}{{\bf{v}}}
\newcommand{\Be}{{\bf{e}}}
\newcommand{\bw}{{\bf{w}}}
\newcommand{\bz}{{\bf{z}}}
\newcommand{\bc}{{\bf{c}}}
\newcommand{\bd}{{\bf{d}}}
\newcommand{\bh}{{\bf{h}}}
\newcommand{\bbe}{{\bf{e}}}
\newcommand{\wh}{\widehat}

\newcommand{\ui}{{\underline i}}
\newcommand{\uj}{{\underline j}}
\newcommand{\ual}{{\underline \al}}
\newcommand{\uz}{{\underline z}}


\newcommand{\Z}{{\mathbb Z}}


\newcommand{\bX}{{\bf{X}}}
\newcommand{\bY}{{\bf{Y}}}
\newcommand{\bZ}{{\mathbb{Z}}}


\newcommand{\wG}{{\widehat G}}
\newcommand{\al}{\alpha}
\newcommand{\de}{\delta}

\newcommand{\be}{\begin{equation}}
\newcommand{\ee}{\end{equation}}


\newcommand{\ga}{{\gamma}}
\newcommand{\Ga}{{\Gamma}}
\newcommand{\la}{\lambda}
\newcommand{\Om}{{\Omega}}
\newcommand{\om}{{\omega}}
\newcommand{\si}{\sigma}
\renewcommand{\th}{\theta}
\newcommand{\td}{\tilde}
\newcommand{\ze}{\zeta}

\newcommand{\cA}{{\cal A}}
\newcommand{\cG}{{\cal G}}
\newcommand{\cF}{{\cal F}}
\newcommand{\cH}{{\cal H}}
\newcommand{\cL}{{\cal L}}
\newcommand{\cE}{{\cal E}}
\newcommand{\cN}{{\cal N}}
\newcommand{\cU}{{\cal U}}
\newcommand{\cW}{{\cal W}}

\newcommand{\im}{{\text{Im} }}
\newcommand{\re}{{\text{Re} }}

\newcommand{\E}{{\mathbb E }}
\newcommand{\bQ}{{\mathbb Q }}
\newcommand{\bN}{{\mathbb N}}
\renewcommand{\P}{{\mathbb P}}
\newcommand{\bC}{{\mathbb C}}
\newcommand{\pd}{{\partial}}
\newcommand{\nb}{{\nabla}}
\newcommand{\lec}{\lesssim}
\newcommand{\ind}{{\,\mathrm{d}}}
%\newcommand{\qed}{\hfill\fbox{}\par\vspace{0.3mm}}
\newcommand{\ph}{{\varphi}}

\renewcommand{\div}{\mathop{\mathrm{div}}}
\newcommand{\curl}{\mathop{\mathrm{curl}}}
\newcommand{\spt}{\mathop{\mathrm{spt}}}
\newcommand{\wkto}{\rightharpoonup}
\newenvironment{pf}{{\bf Proof.}} {\hfill\qed}
\newcommand{\wt}{\widetilde}

\newcommand{\lv}{{\bar v}}
\newcommand{\lp}{{\bar p}}

%DEFINITIONS

\def\a{{\alpha}}
\def\b{{\beta}}
\def\g{{\gamma}}
\def\r{{\rho}}
\def\s{{\sigma}}
\def\o{{\omega}}
\def\n{{\nu}}
\def\D{{\Delta}}
\def\O{{\Omega}}
\def\G{{\Gamma}}
%\def \be{\begin{equation}}
%\def \ee{\end{equation}}
\def\rhs{{\it r.h.s. }}
\def\lhs{{\it l.h.s. }}
\def \bea{\begin{eqnarray}}
\def \eea{\end{eqnarray}}
\def\nn{{\nonumber}}
\def\media#1{{\langle#1\rangle}}


\newcommand{\gibbs}{\frac{e^{-\tb\tdl}}{\tr e^{-\tb \tdl}}}
heisenberg.bbl
\begin{thebibliography}{10}

\bibitem{Bloch1930}
F.~Bloch.
\newblock {Zur Theorie des Ferromagnetismus}.
\newblock {\em Zeitschrift f{\"u}r Physik}, 61:206--219, 1930.

\bibitem{CS90}
J.~G. Conlon and J.~P. Solovej.
\newblock {On asymptotic limits for the quantum Heisenberg model}.
\newblock {\em J.\ Phys.\ A}, 23:3199--3213, 1990.

\bibitem{CS91rw}
J.~G. Conlon and J.~P. Solovej.
\newblock {Random walk representations of the Heisenberg model}.
\newblock {\em J.\ Stat.\ Phys.}, 64:251--270, 1991.

\bibitem{CS91class}
J.~G. Conlon and J.~P. Solovej.
\newblock {Uniform convergence of the free energy of the classical Heisenberg
  model to that of the Gaussian model}.
\newblock {\em J.\ Statist.\ Phys.}, 65:235--245, 1991.

\bibitem{CS91}
J.~G. Conlon and J.~P. Solovej.
\newblock {Upper Bound on the Free Energy of the Spin 1/2 Heisenberg
  Ferromagnet}.
\newblock {\em Lett.\ Math.\ Phys.}, 23:223--231, 1991.

\bibitem{CG12}
M.~Correggi and A.~Giuliani.
\newblock {The Free Energy of the Quantum Heisenberg Ferromagnet at Large
  Spin}.
\newblock {\em J.\ Stat.\ Physics}, 149(2):234--245, 2012.

\bibitem{CGS14}
M.~Correggi, A.~Giuliani, and R.~Seiringer.
\newblock {Validity of the Spin-Wave Approximation for the Free Energy of the
  Heisenberg Ferromagnet}.
\newblock {\em Comm.\ Math.\ Phys.}, 339:279--307, 2015.

\bibitem{Dyson1}
F.~Dyson.
\newblock General theory of spin-wave interactions.
\newblock {\em Phys.\ Rev.}, 102:1217--1230, 1956.

\bibitem{Dyson2}
F.~Dyson.
\newblock {Thermodynamic Behavior of an Ideal Ferromagnet}.
\newblock {\em Phys.\ Rev.}, 102:1230--1244, 1956.

\bibitem{FSS76}
J.~Fr\"ohlich, B.~Simon, and T.~Spencer.
\newblock {Phase Transitions and Continuous Symmetry Breaking}.
\newblock {\em Phys.\ Rev.\ Lett.}, 36:804--806, 1976.

\bibitem{Hewson63}
A.~C. Hewson and D.~Ter Haar.
\newblock {On the $T^3$-term in the low-temperature expansion for the
  magnetisation of a spin-1/2 Heisenberg ferromagnet}.
\newblock {\em Physics Letters}, 6(2):136--137, 1963.

\bibitem{Hofmann02}
C.~P. Hofmann.
\newblock {Spontaneous magnetization of the O(3) ferromagnet at low
  temperatures}.
\newblock {\em Phys.\ Rev.\ B}, 65:094430, 2002.

\bibitem{Hofmann11}
C.~P. Hofmann.
\newblock {Spontaneous magnetization of an ideal ferromagnet: Beyond Dyson's
  analysis}.
\newblock {\em Phys.\ Rev.\ B}, 84:064414, 2011.

\bibitem{HP40}
T.~Holstein and H.~Primakoff.
\newblock {Field Dependence of the Intrinsic Domain Magnetization of a
  Ferromagnet}.
\newblock {\em Phys.\ Rev.}, 58:1098--1113, 1940.

\bibitem{Kramers36}
H.~A. Kramers.
\newblock {Zur Theorie des Ferromagnetismus}.
\newblock {\em Commun.\ Kamerlingh Onnes Lab.\ Univ.\ Leiden}, 22:(Suppl.\ 83)
  67,68,72,75--81, 1936.

\bibitem{Lieb73}
E.~H. Lieb.
\newblock {The Classical Limit of Quantum Spin Systems}.
\newblock {\em Comm.\ Math.\ Phys.}, 31:327--340, 1973.

\bibitem{PD66}
P.~D. Loly and S.~Doniach.
\newblock {Removal of an Apparent Discrepancy between Calculations of Dyson and
  of Oguchi for the Heisenberg Ferromagnet}.
\newblock {\em Phys.\ Rev.}, 144:319--321, 1966.

\bibitem{Oguchi}
T.~Oguchi.
\newblock {Theory of Spin-Wave Interactions in Ferro- and Antiferromagnetism}.
\newblock {\em Phys.\ Rev.}, 117:117--123, 1960.

\bibitem{OtherOguchi}
T.~Oguchi.
\newblock {The Equivalence of Hamiltonians of Holstein-Primakoff and Dyson in
  Spin-Wave Theory in Ferromagnetism}.
\newblock {\em Prog.\ Theor.\ Phys.}, 25:721--722, 1961.

\bibitem{Opechowski37}
W.~Opechowski.
\newblock {Über die Temperaturabhängigkeit der Magnetisierung eines
  Ferromagnetikums bei tiefen Temperaturen}.
\newblock {\em Physica}, 4(8):715--722, 1937.

\bibitem{Schafroth54}
M.~R. Schafroth.
\newblock {Self-Consistent Spin-Wave Theory for the Ferromagnetic Exchange
  Problem}.
\newblock {\em Proceedings of the Physical Society. Section A}, 67(1):33, 1954.

\bibitem{T93}
B.~Toth.
\newblock {Improved Lower Bound on the Thermodynamic Pressure of the Spin 1/2
  Heisenberg Ferromagnet}.
\newblock {\em Lett.\ Math.\ Phys.}, 28:75--84, 1993.

\bibitem{Kranendonk54}
J.~van Kranendonk.
\newblock {Theory of the Low-Temperature Properties of Ferromagnetic Crystals}.
\newblock {\em Physica}, 21(1):81--82, 1954.

\bibitem{Zittartz65}
J.~Zittartz.
\newblock {On the Spin Wave Problem in the Heisenberg model of Ferromagnetism}.
\newblock {\em Zeitschrift f{\"u}r Physik}, 184(5):506--520, 1965.

\end{thebibliography}
