Initial commitv0.0

18 files changed, 3076 insertions, 0 deletions

diff --git a/Carlen_Jauslin_Lieb_2019.tex b/Carlen_Jauslin_Lieb_2019.tex
new file mode 100644
index 0000000..9e020a5
--- /dev/null
+++ b/Carlen_Jauslin_Lieb_2019.tex
@@ -0,0 +1,1437 @@
+\documentclass{ian}
+
+\usepackage{array}
+\usepackage{graphicx}
+\usepackage{largearray}
+\usepackage{dsfont}
+\usepackage{etoolbox}
+
+\begin{document}
+
+\pagestyle{empty}
+
+\hbox{}
+\hfil{\bf\LARGE Analysis of a simple equation for the
+\par
+\vskip10pt
+\hfil ground state energy of the Bose gas
+}
+\vfill
+
+\hfil{\bf\large Eric Carlen}\par
+\hfil{\it Department of Mathematics, Rutgers University}\par
+\hfil{\tt\color{blue}\href{mailto:carlen@rutgers.edu}{carlen@rutgers.edu}}\par
+\vskip20pt
+
+\hfil{\bf\large Ian Jauslin}\par
+\hfil{\it Department of Physics, Princeton University}\par
+\hfil{\tt\color{blue}\href{mailto:ijauslin@princeton.edu}{ijauslin@princeton.edu}}\par
+\vskip20pt
+
+\hfil{\bf\large Elliott H. Lieb}\par
+\hfil{\it Departments of Mathematics and Physics, Princeton University}\par
+\hfil{\tt\color{blue}\href{mailto:lieb@princeton.edu}{lieb@princeton.edu}}\par
+
+\vfill
+
+
+\hfil{\bf Abstract}\par
+\medskip
+{\small
+In 1963 a partial differential equation with a convolution non-linearity was introduced in connection with a quantum mechanical many-body problem, namely the gas of bosonic particles.
+This equation is mathematically interesting for several reasons.
+(1) Further investigation showed that predictions based on the equation agree well over the {\it entire range} of parameters with what is expected to be true for the solution of the true many-body problem.
+(2) The novel nonlinearity is easy to state but seems to have almost no literature up to now.
+(3) The earlier work did not prove existence and uniqueness of a solution, which we provide here along with properties of the solution such as decay at infinity. 
+}
+\vfill
+
+\tableofcontents
+ 
+\vfill
+\hfil{\scriptsize\copyright\ by the authors. This paper may be reproduced in its entirety for non-commercial purposes.}
+\eject
+
+\setcounter{page}1
+\pagestyle{plain}
+
+\section{Introduction}
+\indent This paper is devoted to the study of an integro-differential equation introduced in\-~\cite{Li63} in connection with the study of the Bose gas, a many body problem in quantum mechanics.
+The equation is
+\begin{equation}
+ (-\Delta+4e+\mathcal V(x))u(x)=\mathcal V(x)+2e\rho(u\ast u)(x)\ ,
+ \label{simpleq}
+\end{equation}
+with $x\in\mathbb R^d$, and $*$ denoting convolution: $u\ast u(x):=\int u(x-y)u(y)\ dy$.
+Here, $\mathcal V$ is a given function, (called the {\it potential}), in $L^1 (\mathbb{R}^d ) \cap L^p(\mathbb{R}^d )$,
+with $p >d/2$ for $d\geqslant 2$ and $p>1 $for $d=1$. We assume $\mathcal{V}$ to be non-negative. (This corresponds to a repulsive interaction between the particles in the underlying quantum system). The two parameters
+$e$ and $\rho$ are non-negative numbers, and they are related by a constraint, namely
+\begin{equation}
+ e=\frac\rho2\int (1-u(x))\mathcal V(x)\ dx\ .
+ \label{energy}
+\end{equation}
+
+\indent We are interested in solutions of\-~(\ref{simpleq}) that satisfy the constraint\-~(\ref{energy}), or, in other words, solutions of the system\-~(\ref{simpleq}) and\-~(\ref{energy}). We are particularly interested in the case $d=3$, though other dimensions are also of interest. 
+As explained in \-~\cite{Li63}, the parameter $\rho$ corresponds to the particle density $\frac NV$ of the underlying Bose gas in the large volume and large particle number limit, and $e=\frac EN$ stands for the energy per particle. 
+\bigskip
+
+\indent One would like to fix a value $\rho$ for the density, and then one expects, on the basis of the arguments in \-~\cite{Li63}, that there will be a unique value of $e = e(\rho)$ such that there is a solution of\-~(\ref{simpleq}) and\-~(\ref{energy}) with $u$ taking values in $[0,1]$.
+This value of $e$ is then the energy per particle of the Bose gas in its ground state.
+
+\indent The problem of determining this ground state energy per particle, as a function of the density, has attracted the attention of a great many researchers since the pioneering work of Lenz in 1929\-~\cite{Le29}.
+In that paper and subsequent work\-~\cite{Bo47,LHY57}, an asymptotic expansion of $e(\rho)$, for $d=3$ and small $\rho$ was obtained:
+\begin{equation}
+ e=2\pi\rho a\left(1+\frac{128}{15\sqrt\pi}\sqrt{\rho a^3}+o(\sqrt\rho)\right)\ \qquad \rho \to 0
+ \label{lhy}
+\end{equation}
+where $a$, called the {\it scattering length}, is a property of the pair interaction potential $\mathcal V(x)$, and is defined in\-~(\ref{scattering})-(\ref{ephi}) below.
+Here, we set both the mass $m$ of the particle and Planck's constant $\hbar$ to 1. This early work was not mathematically rigorous, and it was not until
+ 1998\-~\cite{LY98} that the validity of the first term $2\pi\rho a$ was proved, and not until 2019\-~\cite{FS19} that the validity of second term was also proved, utilizing upper bounds proved earlier in\-~\cite{Dy57,YY09}.
+
+\indent This timeline gives some idea of the complexity of the problem of directly studying the Bose gas ground state as a many body problem. The complexity makes it very attractive to try to show that the system\-~(\ref{simpleq}) and 
+ (\ref{energy})
+provides a useful and illuminating route to the computation of the properties of the ground state for a Bose gas. Interest is piqued further by the fact that numerical studies show that the function $e(\rho)$ computed using the system\-~(\ref{simpleq}) and\-~(\ref{energy}) is surprisingly accurate for {\em all } densities, not only low densities, as we discuss later in this paper. Until now, however, there has been no mathematically rigorous study of this system, and even the most basic questions concerning existence and uniqueness of solutions had remained open. 
+\bigskip
+
+\indent In this paper, we settle some of these basic questions and raise others. It may at first appear surprising that the equation\-~(\ref{simpleq}) poses any serious mathematical challenges. After all, if one replaced the convolution nonlinearity $u*u$ in\-~(\ref{simpleq}) by a power non-linearity, say $u^2$, one would have a familiar sort of local elliptic equation:
+\begin{equation}
+ (-\Delta+4e+\mathcal V(x))u(x)=\mathcal V(x)+2e\rho u^2(x)
+ .
+ \label{simpleq2}
+\end{equation}
+However, the convolution nonlinearity in\-~(\ref{simpleq}) makes it non-local, and very different from\-~(\ref{simpleq2}).
+\bigskip
+
+\indent As explained in\-~\cite{Li63} the solutions of physical interest are integrable, and satisfy
+\begin{equation}
+ 0 \leqslant u(x) \leqslant 1 \quad{\rm for \ all}\ x \ .
+ \label{con1}
+\end{equation}
+We shall see that any non-negative solution automatically satisfies this upper bound.
+\bigskip
+
+\indent Before stating our main theorems, we make a few observations.
+\medskip
+
+\point The system\-~(\ref{simpleq})-(\ref{energy}) is actually equivalent to\-~(\ref{simpleq}) and the constraint
+\begin{equation}
+ \int u(x)\ dx=\frac1\rho.
+ \label{con2}
+\end{equation}
+To prove this, consider the operator
+\begin{equation}
+ G_e := [-\Delta + 4e]^{-1}
+ \label{Ge}
+\end{equation}
+which is given by
+\begin{equation}
+ G_e f = Y_{4e}*f
+ \label{yukawa}
+\end{equation}
+where $Y_{4e}$ is the {\em Yukawa potential}\-~\cite[section 6.23]{LL01}, which is non-negative and $\int Y_{4e} dx = (4e)^{-1}$.
+When $d=3$,
+\begin{equation}
+ Y_{4e}(x)=\frac{e^{-2\sqrt e|x|}}{4\pi|x|}
+ .
+\end{equation}
+Equation\-~(\ref{simpleq}) can be rewritten as
+\begin{equation}\label{simpleq3}
+u(x) = Y_{4e}*(\mathcal V (1- u(x)) +2e\rho Y_{4e}*u*u \ .
+\end{equation}
+Since $u$ and $\mathcal V$ are assumed to be integrable, and $u(x)$ is assumed to satisfy\-~(\ref{con1}), all terms in\-~(\ref{simpleq3}) are integrable, and integrating yields 
+\begin{equation}
+\int u(x)\ dx = \frac{1}{4e} \int \mathcal V(x)(1- u(x)) dx +\frac{\rho}{2} \left(\int u(x)\ dx\right)^2\ .
+\end{equation}
+Thus, for integrable solutions $u$ of\-~(\ref{simpleq}) satisfying\-~(\ref{con1}), the constraint\-~(\ref{energy}) is equivalent to\-~(\ref{con2}).
+\bigskip
+
+\point
+There is another useful way to write the system \-~(\ref{simpleq})-(\ref{energy}). The damped heat semigroup $e^{-t(-\Delta + 4e)}$ is a strongly continuous contraction semigroup on $L^p(\mathbb{R}^d)$, and the domain of its generator is 
+$\mathcal{D}(-\Delta + 4e)= W^{2,p}(\mathbb{R}^d)$. 
+By the Sobolev embedding theorem\-~\cite[Theorem 10.2]{LL01}, since $p>d/2$, all functions $f\in \mathcal{D}(-\Delta + 4e)$ are continuous and vanish at infinity. Since $\mathcal{V} \geqslant 0$, $e^{-t\mathcal{V}}$ is also a
+strongly continuous contraction semigroup on $L^p(\mathbb{R}^d)$, and since $\mathcal{V} \in L^p(\mathbb R^d)$, the domain of its generator, $\mathcal{D}(\mathcal{V})$, contains all bounded functions, and in particular
+$W^{2,p}(\mathbb{R}^d)$. Writing $\mathcal{V}$ as the sum of a piece with a small norm in $L^p(\mathbb R^d)$ and another piece that is bounded, it is easy to see that there are numbers $a,b>0$ with $a< 1/2$ such that for all $f\in 
+W^{2,p}(\mathbb{R}^d)$, 
+\begin{equation}
+ \|\mathcal{V}f\|_p \leqslant a\|(-\Delta + 4e)f\|_p + b\|f\|_p\ .
+\end{equation}
+Then by the Banach space version
+of the Kato-Rellich theorem, \cite[p.~244]{RS75b}, the operator $-\Delta + 4e + \mathcal V(x)$ maps 
+$W^{2,p}(\mathbb{R}^d)$ invertibly onto $L^p(\mathbb R^d)$. Define $K_e$ to be the inverse operator
+\begin{equation}
+ K_e := [-\Delta + 4e + \mathcal V(x)]^{-1}\ .
+ \label{Kedef}
+\end{equation}
+By the Trotter product formula, the operator $K_e$ has a positive kernel that we denote by $K_e(x,y)$; in particular, $K_e$ preserves positivity.
+By the resolvent identity
+\begin{equation}
+ K_e=G_e-G_e\mathcal VK_e
+ \label{resolvent}
+\end{equation}
+we conclude that
+\begin{equation}
+ 0\leqslant K_e(x,y)\leqslant G_e(x,y)
+ \label{ineqKe}
+\end{equation}
+for all $x,y$. Thus, the operator $K_e$ extends to a bounded operator on $L^1(\mathbb R^d)$ and all terms in the equation 
+\begin{equation}\label{simpleq4}
+u(x) = K_e \mathcal V(x) +2e\rho K_eu*u(x) \ .
+\end{equation}
+are well-defined whenever $u$ is integrable. Moreover, since
+$\mathcal V \in L^p(\mathbb R^d)$, and since $u\ast u\in L^p(\mathbb R^d)$ when $u$ is integrable and satisfies\-~(\ref{con1}), every integrable solution $u$ of\-~(\ref{simpleq4}) that satisfies\-~(\ref{con1}) actually belongs to
+$W^{2,p}(\mathbb{R}^d)$ and satisfies\-~(\ref{simpleq}).
+\bigskip
+
+\indent Several simple bounds follow almost immediately from this form of the equation.
+First of all, since the last term on the right of\-~(\ref{simpleq4}) is non-negative, we have an {\em a-priori} lower bound on $u(x)$, namely
+\begin{equation}
+ u(x) \geqslant u_1(x) := K_e\mathcal V(x)\ .
+ \label{u1def}
+\end{equation}
+Integrating both sides of\-~(\ref{u1def}), and using\-~(\ref{con2}) yields an upper bound on $\rho$ depending only on $e$, namely, $\rho \leqslant \left(\int K_e\mathcal V(x)\ dx\right)^{-1}$. 
+By\-~(\ref{energy}) and\-~(\ref{u1def}),
+\begin{equation}
+\rho = {2e}\left(\int\mathcal{V}(1- u)(x) dx\right)^{-1} \geqslant {2e}\left(\int\mathcal{V}(1- K_e\mathcal V)(x) dx\right)^{-1} \ .
+\end{equation}
+Altogether,
+\begin{equation}
+{2e}\left(\int\mathcal{V}(1- K_e\mathcal V)(x) dx\right)^{-1} \leqslant \rho \leqslant \left(\int K_e\mathcal V(x)\ dx\right)^{-1} \ .
+ \label{con4}
+\end{equation}
+In fact, the left side of\-~(\ref{con4}) is equal to one half the right side. To see this observe that $u_1 = K_d \mathcal{V}$ satisfies $(-\Delta +4e + \mathcal{V})u_1 = \mathcal{V}$, and hence $u_1 = G_e(\mathcal{V}(1-u_1))$.
+Integrating both side yields $\int u_1 dx = \frac{1}{4e} \int \mathcal{V}(1-u_1) dx$. By\-~(\ref{u1def}), we obtain the following simpler (albeit less sharp) bounds:
+\begin{equation}
+{2e}\left(\int\mathcal{V}dx \right)^{-1} \leqslant \rho \leqslant{4e}\left(\int\mathcal{V}dx \right)^{-1} \ ,
+ \label{con4B}
+\end{equation}
+or equivalently
+\begin{equation}
+\left( \frac{1}{4}\int\mathcal{V}dx \right)\rho \leqslant e \leqslant \left(\frac{1}{2}\int\mathcal{V}dx\right)\rho \ .
+ \label{con4C}
+\end{equation}
+In particular, this shows that the system\-~(\ref{simpleq})-(\ref{energy}) does not have a solution for arbitrary values of $\rho$ and $e$: when either is small, a solution of the type we seek can only exist if the other is correspondingly small, as specified by\-~(\ref{con4B}) and\-~(\ref{con4C}).
+In fact, as is stated in the following theorem, $\rho$ and $e$ are constrained to be related by a functional equation.
+\bigskip
+
+\theo{Theorem}[existence and uniqueness]\label{theorem:existence} Let $\mathcal V\in L^1(\mathbb{R}^d)\cap L^p(\mathbb{R}^d)$, $p>\max\{\frac d2,1\}$, be non-negative. 
+Then there is a constructively defined continuous function $\rho(e)$ on $(0,\infty)$ such that 
+ $\lim_{e\to 0}\rho(e) = 0$ and $\lim_{e\to \infty} \rho(e) = \infty$ and such that for any $e\geqslant 0$ and $\rho = \rho(e)$, 
+ the system\-~(\ref{simpleq}) and\-~(\ref{energy}) has a unique integrable solution $u(x)$ satisfying\-~(\ref{con1}). Moreover, if $\rho \neq \rho(e)$, the system\-~(\ref{simpleq}) and\-~(\ref{energy}) has {\em no} integrable solution $u(x)$ satisfying\-~(\ref{con1}).
+\endtheo
+\bigskip
+
+{\bf Remarks}:
+\begin{itemize}
+ \item We do not assume here that the potential is radially symmetric.
+ However, we shall see from our proof that if $\mathcal V$ is radially symmetric, then so is $u$.
+
+ \item The function $\rho(e)$ is the {\em density function}, which specifies the density as a function of the energy. Thus, our system together with\-~(\ref{con1}) constrains the parameters $e$ and $\rho$ to be related by a strict functional relation $\rho = \rho(e)$. In most of the early literature on the Bose gas, $\rho$ is taken as the independent parameter, as suggested by\-~(\ref{lhy}): One puts $N$ particles in a box of volume $N/\rho$, and 
+seeks to find the ground state energy per particle, $e$, as a function of $\rho$. Our theorem goes in the other direction, with $\rho$ specified as a function of $e$. We prove that $e\mapsto\rho(e)$ is continuous, and we conjecture that $\rho(e)$ is a strictly 
+monotone increasing function. In that case, the functional relation could be inverted, and we would have a well-defined function $e(\rho)$.
+
+ \item Since $ \lim_{e\to 0}\rho(e)=0$ and $\lim_{e\to \infty}\rho(e)=\infty$, the continuity of $e\to \rho(e) $ implies that for each $\rho \in (0,\infty)$ there is {\it at least one} $e$ such that $\rho(e) =\rho$.
+\end{itemize}
+\bigskip
+
+\indent Having proved that the solution to the simple equation is unique, our second main result is an asymptotic expression for $e(\rho)$, both for low and for high density.
+\bigskip
+
+\theo{Theorem}[asymptotics of the energy for $d=3$]\label{theorem:asymptotics}
+ Consider the case $d=3$.
+ Let $\mathcal V$ be non-negative, integrable and square-integrable. Then, for each 
+ $\rho>0$ there is at least one $e>0$ such that $\rho = \rho(e)$. For any such $\rho $ and $e$ we have the following bounds for low and high density (i.e., small and large $\rho$).
+ For low density, 
+ \begin{equation}
+ e=2\pi\rho a\left(1+\frac{128}{15\sqrt\pi}\sqrt{\rho a^3}+o(\sqrt\rho)\right)
+ \label{3d}
+ \end{equation}
+ where $a$ is the {\it scattering length} of the potential, which is defined in~(\ref{sl}).
+ For high density, in any dimension $d\geqslant 1$,
+ \begin{equation}
+ e=\frac\rho2\int\mathcal V(x)\ dx+o(\rho)
+ .
+ \label{largerho}
+ \end{equation}
+\endtheo
+\bigskip
+
+{\bf Remarks}:
+\begin{itemize}
+ \item For low densities in $d=3$, the energy $e$ predicted by the simple equation\-~(\ref{simpleq})-(\ref{energy}) is asymptotically equal to the ground state energy of the Bose gas\-~\cite{LHY57,YY09,FS19}.
+ For high densities, when the potential has a non-negative Fourier transform, the asymptotic formula for the ground state energy of the Bose gas coincides with\-~(\ref{largerho})\-~\cite[appendix]{Li63}.
+ Thus, the simple equation yields the same asymptotes for both low and high densities as the Bose gas does (at least when the potential has a non-negative Fourier transform, as in the example $\mathcal V(x)=e^{-|x|}$ discussed in section\-~\ref{subsec:numerics}).
+\end{itemize}
+\bigskip
+
+\theo{Theorem}[decay of $u$ in $d=3$]\label{theorem:decay}
+ For $d=3$, if $\mathcal V$ is non-negative, square-integrable, spherically symmetric (that is, $\mathcal V(x)=\mathcal V(|x|)$), and, for $|x|>R$,
+ \begin{equation}
+ \mathcal V(|x|)\leqslant Ae^{-B|x|}
+ \label{expdecay}
+ \end{equation}
+ for some $A,B>0$ then there exists $\alpha>0$ such that
+ \begin{equation}
+ u(|x|)\mathop\sim_{|x|\to\infty}\frac\alpha{|x|^4}
+ .
+ \end{equation}
+\endtheo
+\bigskip
+
+{\bf Remarks}:
+\begin{itemize}
+ \item This result is consistent with a prediction in\-~\cite{LHY57} that the truncated 2-point correlation function in the ground state of the Bose gas decays like $|x|^{-4}$.
+
+ \item To prove this theorem, we will use analytical properties of the Fourier transform $\widehat{\mathcal V}$ of $\mathcal V$, which is why we assume that $\mathcal V$ decays exponentially at infinity.
+ For potentials with slower decay, it seems that the decay of $u$ should still be $|x|^{-4}$, except if $\mathcal V$ itself decays slower than $|x|^{-4}$, in which case $u$ should decay like $\mathcal V$.
+
+ \item It is presumably not too difficult to extend this result to cases with potentials that are not spherically symmetric.
+\end{itemize}
+\bigskip
+
+{\bf Remark}:
+The simple equation\-~(\ref{simpleq}) is actually an approximation of a richer equation for $u$\-~\cite{Li63}, which should more accurately depict the Bose gas, see\-~(\ref{fulleq}).
+ Little is known about this richer equation.
+\bigskip
+
+\indent The paper is organized as follows.
+We prove theorems\-~\ref{theorem:existence}, \ref{theorem:asymptotics} and\-~\ref{theorem:decay} in sections\-~\ref{sec:existence}, \ref{sec:asymptotics} and\-~\ref{sec:decay} respectively.
+In section\-~\ref{sec:bosegas}, we explain how the simple equation is related to the Bose gas, and present some numerical evidence that it is very good at predicting the ground state energy.
+In section\-~\ref{sec:open} we discuss a few open problems and extensions.
+
+\section{Proof of Theorem\-~\expandonce{\ref{theorem:existence}}}\label{sec:existence}
+
+\indent As was shown in\-~(\ref{simpleq3}) and\-~(\ref{simpleq4}), there are at least two ways to write\-~(\ref{simpleq}) as a fixed point equation.
+As it turns out, only the latter one
+\begin{equation}
+u(x) = \Phi(u)(x) := K_e (\mathcal V(x) +2e\rho u*u(x))
+\label{iteration0}
+\end{equation}
+is adapted to solution by iteration, because of its monotonicity properties.
+Starting with $u_0(x) = 0$, define
+\begin{equation}
+ u_n(x) = \Phi(u_{n-1})(x)
+\end{equation}
+for $n\geq1$. It is easy to see that for arbitrary $e,\rho \geqslant 0$, this produces a monotone increasing sequence of non-negative integrable functions. Thus, $u(x) := \lim_{n\to \infty}u_n(x)$ will exist, but it need not be integrable and it need not satisfy\-~(\ref{energy}) or\-~(\ref{con1}). 
+
+\indent To bring\-~(\ref{energy}) into the iteration scheme, we take $e$ as the independent parameter, and define a sequence $\{\rho_n\}$ along with the sequence $\{u_n(x)\}$, both depending on $e$, through
+\begin{equation}
+ u_n(x)=K_e \mathcal V(x)+2e\rho_{n-1}K_e u_{n-1}\ast u_{n-1}(x) \ .
+ \quad
+ u_0(x)=0
+ \label{iteration}
+\end{equation}
+and
+\begin{equation}
+ \rho_n:=\frac{2e}{\int(1-u_n(x))\mathcal V(x)} .
+ \label{rhon}
+\end{equation} 
+Comparing\-~(\ref{iteration}) to\-~(\ref{iteration0}), note that the analog of $\Phi$ now depends on $n$.
+
+\theo{Lemma}\label{lem1} Let $\mathcal{V} \in L^1(\mathbb{R}^d)\cap L^p(\mathbb{R}^d)$, $p>\max\{\frac d2,1\}$. Both sequences $\{\rho_n\}$ and $\{u_n\}$ are well defined and increasing, and for all $n$,
+\begin{equation}\label{simple17}
+\int_{\mathbb{R}^d} u_n{\rm d}x < \frac{1}{2e}\int_{\mathbb{R}^d}\mathcal{V}(1-u_n)dx \ .
+\end{equation}
+\endtheo
+\bigskip
+
+{\bf Proof:}
+We proceed by induction. By definition, $u_0 =0$ and $\rho_0 = 2e\left(\int_{\mathbb{R}^d}\mathcal V(x)dx\right)^{-1}$. Also by definition
+$u_1 = K_e\mathcal V \geqslant u_0$ and $\rho_1 = 2e\left( \int \mathcal V(1- K_e\mathcal{V})dx \right)^{-1}$.
+ As noted in the discussion between\-~(\ref{con4}) and\-~(\ref{con4B}), 
+\begin{equation}
+2\int_{\mathbb{R}^d} u_1dx = \frac{1}{e}\int_{\mathbb{R}^d}\mathcal{V}(1-u_1)dx \leq \ \frac{1}{e}\int_{\mathbb{R}^d}\mathcal{V}dx.
+\end{equation}
+Since $t\mapsto t^{-1}$ is monotone decreasing on $(0,\infty)$, this shows that $\rho_1 > \rho_0$, and that\-~(\ref{simple17}) holds for $n=1$. 
+\bigskip
+
+
+\indent Now suppose that $u_{n} \geqslant u_{n-1} \geqslant 0$, $\rho_{n} \geqslant \rho_{n-1} \geqslant 0$, and $\int_{\mathbb{R}^d} u_{n}{\rm d}x < \frac{1}{2e}\int_{\mathbb{R}^d}\mathcal V(1-u_{n})$, all of which we have just verified for $n=1$. Then
+\begin{equation}
+ u_{n+1} = K_e\mathcal V + 2e \rho_n K_e u_{n}*u_{n}(x) \geqslant K_e\mathcal V + 2e\rho_{n-1} K_e u_{n-1}*u_{n-1}(x) = u_n(x)\ ,
+\end{equation}
+and then 
+\begin{equation}
+\int_{\mathbb{R}^d} \mathcal V(1-u_{n+1}){\rm d}x < \int_{\mathbb{R}^d}\mathcal V(1-u_{n}){\rm d}x \ .
+\end{equation}
+Integrating both sides of $u_{n+1} = G_e \mathcal {V}(1- u_{n+1}) + 2e \rho_n G_e u_{n}*u_{n}$ yields,
+\begin{equation}\label{simple15}
+2\int_{\mathbb{R}^d} u_{n+1}{\rm d}x = \frac{1}{2e}\int_{\mathbb{R}^d}\mathcal V(1-u_{n+1}) + \rho_n \left( \int_{\mathbb{R}^d}u_{n}{\rm d}x \right)^2
+\end{equation}
+Then since $\int_{\mathbb{R}^d} u_{n}{\rm d}x < \frac{1}{2e}\int_{\mathbb{R}^d}\mathcal V(1-u_{n}) = \frac{1}{\rho_n}$, (\ref{simple15}) implies 
+\begin{equation}
+ 2\int_{\mathbb{R}^d} u_n{\rm d}x \leq \frac{1}{2e}\int_{\mathbb{R}^d}\mathcal{V}(1-u_n) + \int_{\mathbb{R}^d}u_{n-1}{\rm d}x\ .
+\end{equation}
+Then because $\int_{\mathbb{R}^d} u_{n}{\rm d}x < \int_{\mathbb{R}^d} u_{n+1} {\rm d}x $, we have that $\int_{\mathbb{R}^d} u_{n+1}{\rm d}x < \frac{1}{2e}\int_{v}\mathcal{V}(1-u_{n+1}) $.
+This proves\-~(\ref{simple17}) for $n+1$, and shows that 
+\begin{equation}
+ 0 \leqslant \frac{1}{2e}\int_{\mathbb{R}^d}\mathcal{V}(1-u_{n+1})dx \leqslant \frac{1}{2e}\int_{\mathbb{R}^d}\mathcal{V}(1-u_n)dx \ ,
+\end{equation}
+and then, as before, $\rho_{n+1}\geqslant \rho_n$. 
+\qed
+\bigskip
+
+\theo{Lemma}\label{lem2} Let $\mathcal{V} \in L^1(\mathbb{R}^d)\cap L^p(\mathbb{R}^d)$, $p>\max\{\frac d2,1\}$. Then for all $n$ and $x$, $u_n(x)$ is continuous, vanishing at infinity, and $0 \leqslant u_n(x) \leqslant 1$. 
+\endtheo
+\bigskip
+
+{\bf Proof:} First consider $n=1$. Since $u_n = K_e\mathcal{V}$ with $\mathcal{V}\in L^p(\mathbb{R}^d)$, $u_1\in W^{2,p}(\mathbb{R}^d)$ and 
+\begin{equation}
+\Delta u_1(x) = \mathcal{V}(x)(u_1(x) - 1) + 4e u_1(x)\ .
+\end{equation}
+Since $K_e$ maps $L^p(\mathbb R^d)$ into $W^{2,p}(\mathbb R^d)$, $u_1$ is continuous and vanishes at infinity. Let $A := \{x\ :\ u_1(x) > 1\}$. Then $A$ is open. If $A$ is non-empty, then $u_1$ is subharmonic on $A$, and hence takes on its maximum on the boundary of $A$. Since $u_1$ would equal $1$ on the boundary, this is impossible, and $A$ is empty. This proves the assertion for $n=1$. 
+\bigskip
+
+\indent Now make the inductive hypothesis that $0 \leqslant u_n(x) \leq 1$ for all $x$. Then $\|u_n\|_p^p \leq \|u_n\|_1 \leq\frac{1}{2e}\int_{\mathbb{R}^d}\mathcal{V}dx$. By Young's inequality, $\|u_n\ast u_n\|_p \leq \|u_n\|_p\|u_1\|_1$, and hence 
+$\mathcal{V} + 2e\rho_n u_n\ast u_n \in L^p(\mathbb{R}^d)$. Therefore, 
+$u_{n+1}= K_e(\mathcal{V} + 2e\rho_n u_n\ast u_n) \in W^{2,p}(\mathbb{R}^d)$. It follows as before that $u_{n+1}$ is continuous and vanishing at infinity, and in particular, bounded, and 
+\begin{eqnarray*}
+\Delta u_{n+1}(x) &=& \mathcal{V}(x)(u_n(x) - 1) + 4e u_n(x) - 2e\rho_n u_n\ast u_n\\
+&\geqslant& \mathcal{V}(x)(u_n(x) - 1) + 4e u_n(x) -2e\rho_n\|u_n\|_1 \|u_n\|_\infty\\
+& \geqslant& 
+ \mathcal{V}(x)(u_n(x) - 1) + 4e u_n(x) -2e
+ \end{eqnarray*}
+ where we have used $\rho_n\|u_n\|_1 \leqslant 1$, which is valid on account of\-~(\ref{simple17}). Define $A := \{x:\ u_{n+1}(x) > 1\}$. Then $u_{n+1}$ is subharmonic on $A$, and maximal on the boundary of $A$, where $u_n(x)$ would equal 
+ $1$. This contradiction shows that $\|u_{n+1}\|_\infty \leqslant 1$.
+\qed
+\bigskip
+
+\theo{Lemma}\label{lem3} Let $\mathcal{V} \in L^1(\mathbb{R}^d)\cap L^p(\mathbb{R}^d)$, $p>\max\{\frac d2,1\}$. 
+Now let
+\begin{equation}
+ u(x): = \lim_{n\to\infty}u_n(x) \quad{\rm and}\quad \rho(e) = \lim_{n\to\infty}\rho_n(e)\ .
+\end{equation}
+Then both limits exist, $u\in W^{2,p}(\mathbb{R}^d)$ and $u$ satisfies\-~(\ref{simpleq}), (\ref{energy}) and\-~(\ref{con1}). 
+\endtheo
+\bigskip
+
+{\bf Proof:} 
+By Lemma~\ref{lem1}, both limits exist, and by\-~(\ref{simple17}), $\rho(e) \leqslant \left(\int_{\mathbb{R}^d} K_e\mathcal{V}dx\right)^{-1}$.
+Also by Lemma~\ref{lem1}, $\int_{\mathbb{R}^d} \leqslant \frac{1}{2e}\int_{\mathbb{R}^d}\mathcal{V}(x)dx$,
+$u$ is integrable and $\lim_{n\to\infty}\|u_n - u\|_1 = 0$. Moreover, by Lemma~\ref{lem2}, $0 \leq u \leq 1$, and then $\|u\|_p^p \leq \|u\|_1$ and $\|u_n- u\|_p^p \leq (p+1)\|u_n-u\|_1$, and then by Young's Inequality
+\begin{equation}
+ \|u\ast u - u_n\ast u_n\|_p \leq \|u_n\|_1\|u_n -u\|_p + \leq \|u\|_1\|u_n -u\|_p\ .
+\end{equation}
+Therefore, $\lim_{n\to\infty}(\mathcal{V} +2e\rho_n(e) u_n\ast u_n) = (\mathcal{V} +2e\rho(e) u\ast u)$
+with convergence in $L^p(\mathbb{R}^d)$. Then $\lim_{n\to\infty}K_e (\mathcal{V} +2e\rho_n(e) u_n\ast u_n) = K_e(\mathcal{V} +2e\rho(e) u\ast u)$ with convergence in $W^{2,p}(\mathbb{R}^d)$, and in particular, in
+$L^p(\mathbb{R}^d)$. It now follows that $u = K_e(\mathcal{V} +2e\rho(e) u\ast u) $, and by the Dominated Convergence Theorem, the constraint $\rho = \frac{1}{2e}\int_{\mathbb{R}^d}\mathcal{V}(1-u)dx$ is satisfied.
+By remarks made above, this means that $u$ satisfies\-~(\ref{simpleq}) and\-~(\ref{energy}). \qed
+\bigskip
+
+\theo{Lemma}\label{lem3}
+ For all $e\in (0,\infty)$, the solution $u$ of the system\-~(\ref{simpleq}) and\-~(\ref{energy}) that we have constructed by iteration on Lemma~\ref{lem3} is the unique non-negative integrable solution for $\rho = \rho(e)$.
+ Moreover, there does not exist such any such solution when $\rho \neq \rho(e)$. 
+\endtheo
+\bigskip
+
+{\bf Proof:} 
+ Consider any non-negative solution integrable $\tilde u$, with
+ \begin{equation}
+ \tilde\rho=\frac{2e}{\int(1-\tilde u(x))\mathcal V(x)\ dx}
+ .
+ \end{equation}
+ We first show that $\tilde u\geqslant u_n$ by induction.
+ We have
+ \begin{equation}
+ \tilde u(x)-u_n(x)=2eK_e(\tilde\rho\tilde u\ast\tilde u(x)-\rho_{n-1}u_{n-1}\ast u_{n-1}(x))
+ \end{equation}
+ Since $u_{n-1} =0$, the positivity of $\tilde u$ implies the positivity of $\tilde u(x) - u_1(x)$.
+ If $\tilde u\geqslant u_{n-1}$, then, by\-~(\ref{rhon}), $\tilde\rho\geqslant\rho_{n-1}$, from which $\tilde u\geqslant u_n$ follows easily.
+ This proves that both $\tilde\rho\geqslant\rho$ and $\tilde u\geqslant u$. However, integrating both sides of the latter inequality yields 
+ \begin{equation}
+ \frac{1}{\tilde \rho(e)} = \int\tilde u(x)\ dx \geqslant\int u(x)\ dx=\frac1{\rho(e)}\ .
+ \end{equation}
+ Since $\tilde\rho\geqslant\rho$, equality must hold, and then since $\tilde u\geqslant u$, it must be that 
+ so $u=\tilde u$.
+\qed
+\bigskip 
+ 
+\theo{Lemma}\label{lem4}
+ The function $\rho(e)$ is continuous on $(0,\infty)$, with
+ \begin{equation}
+ \lim_{e\to 0}\rho(e) = 0
+ ,\quad
+ \lim_{e\to\infty}\rho(e) = \infty.
+ \end{equation}
+ In particular, for each $\rho\in (0,\infty)$, there is at least one $e\in (0,\infty)$ such that $\rho = \rho(e)$. 
+\endtheo
+\bigskip
+
+{\bf Proof:}
+ We now turn to the continuity of $e\to \rho(e)$.
+ For $n\in \mathbb{N}$, define functions $a_n(e)$ and $b_n(e)$ by
+ \begin{equation}
+ a_n:=\int u_n(x,e)\ dx \quad{\rm and}\quad b_n(e) = \frac1{2e}\int (1-u_n(x,e))\mathcal V(x)\ dx\ .
+ \end{equation} 
+ where we have temporarily made the dependence of $u_n$ on $e$ explicit. 
+ Note that $b_n(e) = 1/\rho_n(e)$. 
+ $u_1(x,e) = K_e\mathcal{V}$ is continuous in $e$ (and monotone decreasing) for each $x$. A simple induction shows that $u_n(x,e)$ is continuous in $e$ for each $x$.
+ Then since $(1-u_n(x,e))\mathcal V(x) \leqslant \mathcal{V}(x)$, the Dominated Convergence Theorem yields the continuity of $\rho_n(e)$ for each $n$.
+ Writing our iteration in the equivalent form (as in\-~(\ref{simpleq3})):
+ \begin{equation}\label{simpleq3B}
+u_n(x,e) = Y_{4e}*(\mathcal V (1- u_n(x,e)) +2e \frac{1}{b_{n-1}(e)} Y_{4e}\ast u_{n-1}\ast u_{n-1}(x,e) \ ,
+\end{equation}
+and integrating, we obtain
+ \begin{equation}\label{simpleq3C} 
+2a_n(x) = b_n(e) + \frac{1}{b_{n-1}(e)} a_{n-1}^2(e) \ ,
+\end{equation}
+Now an easy induction shows that $a_n(e)$ is continuous for each $n$. By\-~(\ref{simple17}), for each $n$,
+ \begin{equation}
+ a_n(e)\leqslant\frac{1}{\rho(e)}\leqslant b_n(e) \ .
+ \end{equation}
+ By Lemma\-~\ref{lem1}, as $n$ increases to infinity, $a_n(e)$ increases to $1/\rho(e)$, while $b_n(e) $ decreases to $1/\rho(e)$. It remains to show that this convergence is uniform on any compact interval in $(0,\infty)$.
+ By\-~(\ref{simpleq3C}), 
+ \begin{equation}\label{telescope}
+ \frac{1}{b_n(e)}(a_n(e) - b_n(e))^2 = \frac{a_n^2(e)}{b_n(e)} -( 2a_n(e) - b_n(e)) = \frac{a_n^2(e)}{b_n(e)} - \frac{a_{n-1}^2(e)}{b_{n-1}(e)} \ .
+ \end{equation}
+ Sum both sides over $n\in \mathbb{N}$. The sum on the right telescopes, and since for all $e$, $a_0^2/b_0 = 0$ while $\lim_{n\to\infty}a_n^2(e)/b_n(e) = 1/\rho_n(e)$, 
+ \begin{equation}
+ \sum_{n=1}^\infty \frac{1}{b_n(e)}(a_n(e) - b_n(e))^2 = \frac{1}{\rho(e)} \ .
+ \end{equation}
+ By the bounds on $b_(e) = 1/\rho_n(e)$ and $\rho(e)$ provided by Lemma~\ref{lem1}, for all $e>0$,
+ \begin{equation}
+ \sum_{n=1}^\infty (a_n(e) - b_n(e))^2 \leqslant\frac{ \int \mathcal{V}dx}{\int K_e\mathcal{V}dx}\ ,
+ \end{equation}
+ and on any compact interval $[e_1,e_2]$, the right hand side is uniformly bounded by $C$, its value at $e_2$. Then since the summand on the left is monotone decreasing in $n$, we obtain for each $n$ that 
+ \begin{equation}\label{rate}
+ (a_n(e) - b_n(e))^2 \leqslant \frac{C}{n}
+ \end{equation} uniformly on $[e_1,e_2]$. This proves the desired uniform convergence, and hence the continuity of $\rho(e)$. The final statement now follows from 
+ (\ref{con4B}). \qed
+\bigskip
+ 
+{\bf Remark}: Note that $\|u - u_n\|_1 = \frac{1}{\rho} - a_n$, and hence by By\-~(\ref{rate}), $\|u - u_n\|_1 \leq Cn^{-1/2}$. 
+In fact, numerically, we find that the rate is significantly faster than this. For example, with $\mathcal V(x)=e^{-|x|}$ and $e=10^{-4}$, $\|u - u_n\|_1$ decays at least as fast as $n^{-3.5}$.
+\bigskip
+
+{\bf Proof of Theorem\-~\ref{theorem:existence}} Every statement in the theorem has been established in Lemma~\ref{lem1} through Lemma~\ref{lem4}. \qed
+
+\bigskip
+
+We close this section by remarking that if $\mathcal{V}$ is radially symmetric, then so is each $u_1 = K_e\mathcal{V}$, and then by a simple induction, so is each $u_n$, hence also $u$, the unique solution $u$ provided by Theorem\-~\ref{theorem:existence}.
+This justifies the first remark following Theorem\-~\ref{theorem:existence}.
+
+\section{Asymptotics}\label{sec:asymptotics}
+\indent In this section, we prove Theorem\-~\ref{theorem:asymptotics}.
+We will first prove the high density asymptote\-~(\ref{largerho}), and then proceed to the low density\-~(\ref{3d}).
+\bigskip
+
+\indent
+By Theorem\-~\ref{theorem:existence}, for each $\rho>0$ there exists at least one $e$ such that $\rho(e)=\rho$.
+If there is more than one, the theorems proved in this section apply to every such solution.
+Throughout this section, let $u_\rho$ denote the solution provided by Theorem\-~\ref{theorem:existence} and any such choice of $e$.
+\bigskip
+
+\subsection{High density $\rho$}
+
+\theo{Lemma}[high density asymptotics]\label{lemma:large}
+ If $\mathcal V$ is integrable, then as $\rho\to\infty$,
+ \begin{equation}
+ e=\frac\rho2\left(\int \mathcal V(x)\ dx\right)(1+o(1))
+ .
+ \end{equation}
+\endtheo
+\bigskip
+
+{\bf Remark}: From\-~(\ref{energy}),
+\begin{equation}
+ e\leqslant\frac\rho2\int\mathcal V(x)\ dx
+ .
+\end{equation}
+Note that this is not an optimal bound, as follows from\-~(\ref{con4}).
+\bigskip
+
+{\bf Proof}:
+ By\-~(\ref{energy}), it suffices to prove that
+ \begin{equation}
+ \lim_{\rho\to\infty}\int u_\rho(x)\mathcal V(x)\ dx
+ =0
+ .
+ \label{limintuv}
+ \end{equation}
+ Let
+ \begin{equation}
+ \chi_\gamma:=\{x:\ \mathcal V(x)\geqslant \gamma\}
+ \end{equation}
+ and decompose
+ \begin{equation}
+ \int u_\rho(x)\mathcal V(x)\ dx
+ =
+ \int_{\chi_\gamma} u_\rho(x)\mathcal V(x)\ dx
+ +
+ \int_{\mathbb R^d\setminus\chi_\gamma} u_\rho(x)\mathcal V(x)\ dx
+ \end{equation}
+ which, by\-~(\ref{con2}), is bounded as follows
+ \begin{equation}
+ \int u_\rho(x)\mathcal V(x)\ dx
+ \leqslant
+ \int_{\chi_\gamma}\mathcal V(x)\ dx
+ +
+ \frac \gamma\rho
+ .
+ \end{equation}
+ Since $\mathcal V$ is integrable, $\int_{\chi_\gamma}\mathcal V(x)\ dx\to0$ as $\gamma\to\infty$.
+ Therefore,
+ \begin{equation}
+ \mathop{\mathrm{inf}}_{\gamma>0}
+ \left(
+ \int_{\chi_\gamma}\mathcal V(x)\ dx
+ +
+ \frac \gamma\rho
+ \right)
+ \mathop{\longrightarrow}_{\rho\to\infty}0
+ .
+ \end{equation}
+\qed
+
+\subsection{Low density $\rho$}\label{subsec:lowrho} In this section, we only consider the dimension $d=3$. As before, we suppose that $\mathcal{V}\in L^1(\mathbb{R}^3)\cap L^p(\mathbb{R}^3)$, $p> 3/2$, and $\mathcal{V} \geqslant 0$.
+
+\indent We first recall the definition of the {\it scattering length} of the potential $\mathcal{V}$, and relate it to the solution of the system\-~(\ref{simpleq})-(\ref{energy}) The {\it scattering equation} is defined as
+\begin{equation}
+ -\Delta\varphi(x)=(1-\varphi(x))\mathcal V(x)
+ ,\quad
+ \lim_{|x|\to\infty}\varphi(x)=0\ .
+ \label{scattering}
+\end{equation}
+Note that\-~(\ref{scattering}) can be written as $(-\Delta + \mathcal V)\varphi = \mathcal{V}$, and hence the solution is
+\begin{equation}\label{scattering2}
+ \varphi(x) = \lim_{e\downarrow 0} K_e \mathcal{V}(x) = \lim_{e\downarrow 0}u_1(x,e)\ ,
+ \end{equation}
+ where $u_1$ is the first term of the iteration introduced in the previous section. It follows from Lemma~\ref{lem2} that $0 \leq \varphi(x) \leq1$ for all $x$. 
+ 
+\indent We now impose a mild localization hypothesis on $\mathcal{V}$: For 
+$R>0$ define $\mathcal{V}_R(x) = \mathcal{V}(x)$ for $|x| > R$ and otherwise $\mathcal{V}_R(x) =0$. We require that for some $q>1$ and all sufficiently large $R$,
+\begin{equation}\label{local}
+\|\mathcal{V}_R\|_1 < R^{-q} \quad{\rm and}\quad \|\mathcal{V}_R \|_p < R^{-q} \ .
+\end{equation}
+By the lemma below, $\lim_{x|\to\infty}|x|\varphi(x)$ exists. The {\em scattering length} $a$ is
+ defined to be (in dimension $d=3$).
+ \begin{equation}\label{sl}
+ a= \lim_{|x|\to \infty} |x| \phi(x).
+ \end{equation}
+For more information on the scattering length, see\-~\cite[appendix A]{LY01}.
+
+ 
+
+\theo{Lemma}\label{sctlem} Let $\mathcal{V}\in L^1(\mathbb{R}^3)\cap L^p(\mathbb{R}^3)$, $p> 3/2$, and suppose that
+the localization condition\-~(\ref{local}) is satisfied. Let $\varphi$ be the corresponding scattering solution given by\-~(\ref{scattering2}). Then the scattering length $a := \lim_{|x|\to\infty}\varphi(x)$ exists and satisfies
+\begin{equation}\label{ephi}
+4\pi a = \int \mathcal{V}(x)(1-\varphi(x))dx\ .
+\end{equation}. 
+\endtheo 
+
+{\bf Proof:} By the resolvent identity, 
+$ \varphi(x) = G\ast (\mathcal{V}(1- \varphi))(x)$
+where $G(x) = \frac{1}{4\pi|x|}$. Since $p> 3/2$. $p' < 3$, and it is easy to decompose $G$ into the sum of two pieces,
+$G = G_1+G_2$ where $G_1 \in L^{p'}(\mathbb{R}^d)$ and $G_2 \in L^{4}(\mathbb{R}^d)$.
+Then for all $R$ sufficiently large,
+\begin{equation}
+ 0 \leq G\ast (\mathcal{V_R}(1- \varphi))(x) \leq (\|G_1\|_{p'} + \|G_2\|_4) R^{-q}\ .
+\end{equation}
+For $0 < r < 1$, then for $|y| < r|x|$, 
+${\displaystyle \frac{1}{1+r} \leqslant \frac{|x|}{|x-y|} \leqslant \frac{1}{1-r}}$.
+It follows that for all sufficiently large $|x|$, 
+\begin{equation}
+\frac{1}{1+r}\int_{|y|< r|x|}\mathcal{V}(y)(1-\varphi(y)) dy+ o(1) \leqslant 4\pi|x|\varphi(x) \leqslant 
+ \frac{1}{1-r}\int_{|y|< r|x|}\mathcal{V}(y)(1-\varphi(y) dx + o(1)\ .
+\end{equation}
+Taking $|x|\to \infty$, and then $r \to 0$ proves\-~(\ref{local}).\qed
+\bigskip
+
+{\bf Remark}: The following lemma is valid if the scattering length $a$ were {\em defined} by (\ref{ephi}).
+For this reason, we do not impose the additional  condition (\ref{local}) in the statement of Theorem\-~\ref{theorem:asymptotics}: Lemma\-~\ref{sctlem} reconciles the stated definition with the formula (\ref{ephi}).
+\bigskip
+
+\theo{Lemma}[low density asymptotics]\label{lemma:small}
+ If $\mathcal V$ is non-negative and integrable and $d=3$, then
+ \begin{equation}
+ e=2\pi\rho a\left(1+\frac{128}{15\sqrt\pi}\sqrt{\rho a^3}+o(\sqrt\rho)\right)
+ .
+ \label{3dlemma}
+ \end{equation}
+\endtheo
+\bigskip
+
+{\bf Proof}:
+ The scheme of the proof is as follows.
+ We first approximate the solution $u$ by $w$, which is defined as the decaying solution of
+ \begin{equation}
+ -\Delta w_\rho(x)=(1-u_\rho(x))\mathcal V(x)
+ .
+ \label{u1}
+ \end{equation}
+ The energy of $w_\rho$ is defined to be
+ \begin{equation}
+ e_w:=\frac\rho2\int(1-w_\rho(x))\mathcal V(x)\ dx
+ \label{ew}
+ \end{equation}
+ and, as we will show, it is {\it close} to $e$, more precisely,
+ \begin{equation}
+ e-e_w=\frac{16\sqrt 2 e^{\frac32}}{15\pi^2}\int\mathcal V(x)\ dx+o(\rho^{\frac32})
+ .
+ \label{eew}
+ \end{equation}
+ In addition, (\ref{u1}) is quite similar to the scattering equation\-~(\ref{scattering}).
+ In fact we will show that $e_w$ is {\it close} to the energy $2\pi\rho a$ of the scattering equation
+ \begin{equation}
+ e_w-2\pi\rho a
+ =
+ -\frac{16\sqrt 2e^{\frac32}}{15\pi^2}\int \varphi(x) \mathcal V(x)\ dx+o(\rho^{\frac32})
+ .
+ \label{ewephi}
+ \end{equation}
+ Summing\-~(\ref{eew}) and\-~(\ref{ewephi}), we find
+ \begin{equation}
+ e=2\pi\rho a\left(1+\frac{32\sqrt 2 e^{\frac32}}{15\pi^2\rho}+o(\sqrt\rho)\right)
+ ,
+ \end{equation}
+ from which\-~(\ref{3dlemma}) follows.
+ We are thus left with proving\-~(\ref{eew}) and\-~(\ref{ewephi}).
+ \bigskip
+
+ \point{\bf Proof of\-~(\ref{eew})}.
+ By\-~(\ref{energy}) and\-~(\ref{ew}),
+ \begin{equation}
+ e-e_w=\frac\rho2\int(w_\rho(x)-u_\rho(x))\mathcal V(x)\ dx
+ .
+ \end{equation}
+ We will work in Fourier space
+ \begin{equation}
+ \hat u_\rho(k):=\int e^{ikx}u_\rho(x)\ dx
+ \label{fourieru}
+ \end{equation}
+ which satisfies, by\-~(\ref{simpleq}),
+ \begin{equation}
+ (k^2+4e)\hat u_\rho(k)
+ =\frac{2e}\rho S(k)
+ +2e\rho \hat u^2(k)
+ \end{equation}
+ with
+ \begin{equation}
+ S(k):=\frac\rho{2e}\int e^{ikx}(1-u_\rho(x))\mathcal V(x)\ dx
+ .
+ \label{S}
+ \end{equation}
+ Therefore,
+ \begin{equation}
+ \hat u_\rho(k)=\frac1\rho\left(\frac{k^2}{4e}+1-\sqrt{\left(\frac{k^2}{4e}+1\right)^2-S(k)}\right)
+ .
+ \label{hatu}
+ \end{equation}
+ Similarly, the Fourier transform of $w_\rho$ is
+ \begin{equation}
+ \hat w_\rho(k):=\int e^{ikx}w_\rho(x)\ dx=\frac{2e S(k)}{\rho k^2}
+ .
+ \label{u1hat}
+ \end{equation}
+ Note that, as $|k|\to\infty$, $\hat u\sim \frac{2eS(k)}{\rho k^2}$, so, while $\hat u_\rho$ is not integrable, $\hat u_\rho-\hat w_\rho$ is.
+ We invert the Fourier transform:
+ \begin{equation}
+ u_\rho(x)-w_\rho(x)=\frac1{8\pi^3\rho}\int e^{-ikx}
+ \left(\frac{k^2}{4e}+1-\sqrt{\left(\frac{k^2}{4e}+1\right)^2-S(k)}-\frac{2eS(k)}{k^2}\right)\ dk
+ .
+ \end{equation}
+ We change variables to $\tilde k:=\frac k{2\sqrt e}$:
+ \begin{equation}
+ u_\rho(x)-w_\rho(x)
+ =
+ \frac{e^{\frac32}}{\rho\pi^3}\int e^{-i2\sqrt e\tilde kx}
+ \left(\tilde k^2+1-\sqrt{(\tilde k^2+1)^2-S(2\tilde k\sqrt e)}-\frac{S(2\tilde k\sqrt e)}{2\tilde k^2}\right)\ d\tilde k
+ .
+ \end{equation}
+ Furthermore,
+ \begin{equation}
+ s\mapsto\left|\tilde k^2+1-\sqrt{(\tilde k^2+1)^2-s}-\frac{s}{2\tilde k^2}\right|
+ \end{equation}
+ is monotone increasing.
+ In addition, by~\-(\ref{S}) and~\-(\ref{simpleq}), and using the fact that $u_\rho(x)\leqslant 1$ (see Lemma\-~\ref{lem2}) and $\mathcal V(x)\geqslant 0$,
+ \begin{equation}
+ |S(k)|\leqslant\frac\rho{2e}\int|(1-u_\rho(x))\mathcal V(x)|\ dx=1
+ .
+ \end{equation}
+ Therefore
+ \begin{equation}
+ \left|
+ \tilde k^2+1-\sqrt{(\tilde k^2+1)^2-S(2\tilde k\sqrt e)}-\frac{S(2\tilde k\sqrt e)}{2\tilde k^2}
+ \right|
+ \leqslant
+ \left|
+ \tilde k^2+1-\sqrt{(\tilde k^2+1)^2-1}-\frac{1}{2\tilde k^2}
+ \right|
+ .
+ \end{equation}
+ Therefore,
+ \begin{equation}
+ |u_\rho(x)-w_\rho(x)|\leqslant
+ \frac{e^{\frac32}}{\rho\pi^3}\int\left|\tilde k^2+1-\sqrt{(\tilde k^2+1)^2-1}-\frac{1}{2\tilde k^2}\right|\ d\tilde k
+ =\frac{32\sqrt 2 e^{\frac32}}{15\pi^2\rho}
+ .
+ \end{equation}
+ By dominated convergence, and using the fact that $S(0)=1$,
+ \begin{equation}
+ \begin{largearray}
+ \lim_{e\to 0}\frac1{e^{\frac32}}(e-e_w)
+ =-\lim_{e\to 0}\frac\rho{2e^{\frac32}}\int (u_\rho(x)-w_\rho(x))\mathcal V(x)\ dx
+ \\[0.5cm]\hfill
+ =
+ -\frac12\int\mathcal V(x)\left(
+	\frac1{\pi^3}\int\left(\tilde k^2+1-\sqrt{(\tilde k^2+1)^2-1}-\frac1{2\tilde k^2}\right)d \tilde k
+ \right)\ dx
+ =\frac{16\sqrt 2}{15\pi^2}\int\mathcal V(x)\ dx
+ .
+ \end{largearray}
+ \end{equation}
+ Using\-~(\ref{con4C}), this proves\-~(\ref{eew}).
+ Incidentally, again by dominated convergence,
+ \begin{equation}
+ u_\rho(x)-w_\rho(x)=
+ \frac{e^{\frac32}}{\rho\pi^3}\int\left(\tilde k^2+1-\sqrt{(\tilde k^2+1)^2-1}-\frac{1}{2\tilde k^2}\right)\ d\tilde k
+ =-\frac{32\sqrt 2 e^{\frac32}}{15\pi^2\rho}+\sqrt\rho f_\rho(x)
+ \label{uapproxw}
+ \end{equation}
+ with
+ \begin{equation}
+ 0\leqslant f_\rho(x)\leqslant\frac{32\sqrt2 e^{\frac32}}{15\pi^2\rho}
+ ,\quad
+ f_\rho(x)\mathop{\longrightarrow}_{\rho\to0}0
+ \end{equation}
+ pointwise in $x$.
+
+ \bigskip
+
+ \point{\bf Proof of\-~(\ref{ewephi})}.
+ Let
+ \begin{equation}
+ \xi(r):=w_\rho(r)-\varphi(r).
+ \label{xi}
+ \end{equation}
+ By\-~(\ref{u1}), (\ref{scattering}) and\-~(\ref{simpleq}),
+ \begin{equation}
+ (-\Delta+\mathcal V(x))\xi(x)=-(u_\rho(x)-w_\rho(x))\mathcal V(x)
+ .
+ \end{equation}
+ Therefore, by\-~(\ref{ephi}),
+ \begin{equation}
+ e_w-2\pi\rho a=-\frac\rho2\int \xi(x)\mathcal V(x)\ dx
+ =-\frac\rho2\int\mathcal V(x)(-\Delta+\mathcal V)^{-1}((u-w)\mathcal V)(x)\ dx
+ \end{equation}
+ and
+ \begin{equation}
+ (-\Delta+\mathcal V)^{-1}\mathcal V(x)=\varphi(x)
+ \end{equation}
+ so
+ \begin{equation}
+ e_w-2\pi\rho a
+ =-\frac\rho2\int\varphi(x)(u_\rho(x)-w_\rho(x))\mathcal V(x)\ dx
+ .
+ \end{equation}
+ By\-~(\ref{uapproxw}),
+ \begin{equation}
+ e_w-2\pi\rho a
+ =\frac{16\sqrt 2e^{\frac32}}{15\pi^2}\int\varphi(x)\mathcal V(x)\ dx
+ -\frac{\rho^{\frac32}}2\int\varphi(x)f_\rho(x)\mathcal V(x)\ dx
+ .
+ \end{equation}
+ Since $x\mapsto f_\rho(x)$ is bounded, we can use dominated convergence to show\-~(\ref{ewephi}).
+\qed
+
+\section{Decay of $u$}\label{sec:decay}
+\indent In this section, we prove Theorem\-~\ref{theorem:decay}.
+Our proof assumes that $\mathcal V$ decays exponentially, because we will use analyticity properties of the Fourier transform of the potential $\mathcal V$.
+In particular, the theorem holds if $\mathcal V$ has compact support.
+We expect the result to hold for any potential that decays faster than $|x|^{-4}$.
+Algebraic decay for $u$ seems natural: by\-~(\ref{simpleq}), $u\ast u$ must decay at infinity in the same way as $u$.
+This is the case if $u$ decays algebraically, but would not be so if, say, it decayed exponentially.
+\bigskip
+
+\indent\underline{Proof of Theorem\-~\ref{theorem:decay}}:
+ We recall that the Fourier transform of $u$ (\ref{fourieru}) satisfies\-~(\ref{hatu}):
+ \begin{equation}
+ \hat u(|k|)=\frac1\rho\left(\frac{k^2}{4e}+1-\sqrt{\left(\frac{k^2}{4e}+1\right)^2-S(|k|)}\right)
+ \end{equation}
+ where $S$ was defined in\-~(\ref{S}):
+ \begin{equation}
+ S(|k|):=\frac\rho{2e}\int e^{ikx}(1-u(|x|))\mathcal V(|x|)\ dx
+ .
+ \end{equation}
+ We split
+ \begin{equation}
+ \hat u(|k|)=\widehat{\mathcal U}_1(|k|)+\widehat{\mathcal U}_2(|k|)
+ \label{hatusplit}
+ \end{equation}
+ with
+ \begin{equation}
+ \widehat{\mathcal U}_1(|k|):=\frac{2eS(|k|)}{\rho(1+k^2)}
+ \end{equation}
+ so that, taking the large $|k|$ limit in\-~(\ref{hatu}),
+ \begin{equation}
+ \widehat{\mathcal U}_2(|k|)=O(|k|^{-4}S^2(|k|))
+ \end{equation}
+ so $\widehat{\mathcal U}_2$ is integrable.
+ \bigskip
+
+ \point{\bf Decay of $\mathcal U_1$}. We first show that
+ \begin{equation}
+ \mathcal U_1(|x|):=\frac1{(2\pi)^3}\int e^{-ikx}\widehat{\mathcal U}_1(|k|)\ dk
+ \end{equation}
+ decays exponentially in $|x|$. We have
+ \begin{equation}
+ \mathcal U_1(|x|)=(-\Delta+1)^{-1}(1-u(|x|))\mathcal V(|x|)
+ =Y_1\ast((1-u)\mathcal V)(|x|)
+ \end{equation}
+ with
+ \begin{equation}
+ Y_1(|x|):=\frac{e^{-|x|}}{4\pi|x|}
+ .
+ \end{equation}
+ Therefore, by\-~(\ref{expdecay}),
+ \begin{equation}
+ \mathcal U_1(|x|)
+ \leqslant
+ \frac A{4\pi}\int_{|y|>R} \frac{e^{-|x-y|-B|y|}}{|x-y|}\ dy
+ +
+ \frac 1{4\pi}\int_{|y|<R} \frac{e^{-|x-y|}}{|x-y|}\mathcal V(|y|)\ dy
+ \end{equation}
+ so, denoting $b:=\min(B,1)$,
+ \begin{equation}
+ \mathcal U_1(|x|)\leqslant
+ \frac A{4\pi}\int \frac{e^{-b(|x-y|-|y|)}}{|x-y|}\ dy
+ +
+ \frac{e^{-(|x|-R)}}{4\pi(|x|-R)}\int\mathcal V(|y|)\ dy
+ \end{equation}
+ and since
+ \begin{equation}
+ \frac A{4\pi}\int \frac{e^{-b(|x-y|-|y|)}}{|x-y|}\ dy
+ =
+ \frac{Ae^{-b|x|}}{4b^2}(b|x|+1)
+ \end{equation}
+ we have
+ \begin{equation}
+ \mathcal U_1(|x|)\leqslant
+ \frac{Ae^{-b|x|}}{4b^2}(b|x|+1)
+ +
+ \frac{e^{-(|x|-R)}}{4\pi(|x|-R)}\int\mathcal V(|y|)\ dy
+ .
+ \label{U1}
+ \end{equation}
+ \bigskip
+
+ \point{\bf Analyticity of $\mathcal U_2$}.
+ We now turn to
+ \begin{equation}
+ \mathcal U_2(|x|):=\frac1{(2\pi)^3}\int e^{-ikx}\widehat{\mathcal U}_2(|k|)\ dk
+ =\frac1{4i\pi^2|x|}\sum_{\eta=\pm}\eta\int_0^\infty e^{i\eta\kappa|x|}\kappa\widehat{\mathcal U}_2(\kappa)\ d\kappa
+ .
+ \label{U2}
+ \end{equation}
+ We start by proving some analytic properties of $\widehat{\mathcal U}_2$, which, we recall from\-~(\ref{hatu}) and\-~(\ref{hatusplit}), is
+ \begin{equation}
+ \widehat{\mathcal U}_2(|k|)=\frac1\rho\left(\frac{k^2}{4e}+1-\sqrt{\left(\frac{k^2}{4e}+1\right)^2-S(|k|)}-\frac{2eS(|k|)}{1+k^2}\right)
+ .
+ \end{equation}
+ \bigskip
+
+ \subpoint First of all, $S$ is analytic in a strip about the real axis:
+ \begin{equation}
+ S(\kappa)=4\pi\int_0^\infty \mathrm{sinc}(\kappa r)r^2\mathcal V(r)(1-u(r))\ dr
+ ,\quad
+ \mathrm{sinc}(\xi):=\frac{\sin(\xi)}\xi
+ \end{equation}
+ so
+ \begin{equation}
+ \partial^nS(\kappa)=4\pi\int_0^\infty \partial^n\mathrm{sinc}(\kappa r)r^{n+2}\mathcal V(r)(1-u(r))\ dr
+ .
+ \end{equation}
+ We will show that if $\mathcal Im(\kappa)\leqslant \frac B2$ (the factor $\frac12$ can be improved to any factor that is $<1$, but this does not matter here), then there exists $C>0$ which only depends on $A$ and $B$ such that
+ \begin{equation}
+ |\partial^nS(\kappa)|\leqslant n!C^n
+ .
+ \label{decaydS}
+ \end{equation}
+ As a consequence, $S$ is analytic in a strip around the real line of height $\frac B2$.
+ In particular, if we define the strip
+ \begin{equation}
+ H_\tau:=\{z:\ |\mathcal Im(z)|\leqslant r^{-\tau},\ \mathcal Re(z)>0\}
+ \end{equation}
+ with $0<\tau<1$, and take
+ \begin{equation}
+ r>\left(\frac B2\right)^{-\frac1\tau}
+ \end{equation}
+ then $S$ is analytic in $H_\tau$.
+ \bigskip
+
+ \bigskip
+
+ \subsubpoint We now prove\-~(\ref{decaydS}).
+ We first treat the case $|\kappa|\leqslant\frac B2$. We have
+ \begin{equation}
+ \mathrm{sinc}(\xi)=\sum_{p=0}^\infty\frac{(-1)^p\xi^{2p}}{(2p+1)!}
+ \end{equation}
+ so
+ \begin{equation}
+ \partial^n\mathrm{sinc}(\xi)=\sum_{p=\lceil\frac n2\rceil}^\infty\frac{(-1)^p\xi^{2p-n}}{(2p+1)(2p-n)!}
+ .
+ \end{equation}
+ Therefore
+ \begin{equation}
+ |\partial^n\mathrm{sinc}(\xi)|\leqslant
+ \sum_{p=\lceil\frac n2\rceil}^\infty\frac{|\xi|^{2p-n}}{(2p-n)!}
+ \leqslant \cosh(|\xi|)
+ .
+ \end{equation}
+ Thus,
+ \begin{equation}
+ |\partial^n S(\kappa)|\leqslant 4\pi\int_0^\infty \cosh(|\kappa|r)r^{n+2}\mathcal V(r)(1-u(r))\ dr
+ \end{equation}
+ so, by\-~(\ref{expdecay}),
+ \begin{equation}
+ |\partial^n S(\kappa)|\leqslant
+ 4A\pi\int_R^\infty \cosh(|\kappa|r)r^{n+2}e^{-Br}\ dr
+ +
+ 4\pi\int_0^R \cosh(|\kappa|r)r^{n+2}\mathcal V(r)\ dr
+ \end{equation}
+ and
+ \begin{equation}
+ |\partial^n S(\kappa)|\leqslant
+ 8A\pi\int_0^\infty r^{n+2}e^{-(B-|\kappa|)r}\ dr
+ +
+ 8\pi e^{|\kappa|R}R^{n}\int r^2\mathcal V(r)\ dr
+ \end{equation}
+ which, if $|\kappa|\leqslant \frac B2$, implies that
+ \begin{equation}
+ 8A\pi\int_0^\infty r^{n+2}e^{-(B-|\kappa|)r}\ dr
+ \leqslant
+ 8A\pi\int_0^\infty r^{n+2}e^{-\frac B2r}\ dr
+ =
+ \frac{2^{n+6}A\pi}{B^{n+3}}(n+2)!
+ \end{equation}
+ and
+ \begin{equation}
+ 8\pi e^{|\kappa|R}R^{n+2}\int\mathcal V(r)\ dr
+ \leqslant
+ 8\pi e^{\frac B2R}R^{n}\int r^2\mathcal V(r)\ dr
+ \end{equation}
+ which implies\-~(\ref{decaydS}) in this case.
+ \bigskip
+
+ \subsubpoint We now turn to $|\kappa|\geqslant \frac B2$:
+ \begin{equation}
+ \partial^n\mathrm{sinc}(\xi)=\sum_{p=0}^n{n\choose p}\partial^p\sin(\xi)\frac{(n-p)!(-1)^{n-p}}{\xi^{n-p+1}}
+ \end{equation}
+ so
+ \begin{equation}
+ |\partial^n\mathrm{sinc}(\xi)|\leqslant 2e^{\mathcal Im(\xi)}\sum_{p=0}^n\frac{n!}{p!}|\xi|^{-(n-p+1)}
+ .
+ \end{equation}
+ Therefore,
+ \begin{equation}
+ |\partial^nS(\kappa)|\leqslant
+ 8\pi \sum_{p=0}^n\frac{n!}{p!|\kappa|^{n-p+1}}
+ \int_0^\infty e^{\mathcal Im(\kappa)r}r^{p+1}\mathcal V(r)(1-u(r))\ dr
+ \end{equation}
+ so, by\-~(\ref{expdecay}),
+ \begin{equation}
+ |\partial^nS(\kappa)|\leqslant
+ \sigma_1+\sigma_2
+ \end{equation}
+ with
+ \begin{equation}
+ \sigma_1:=8A\pi\sum_{p=0}^n\frac{n!}{p!|\kappa|^{n-p+1}}
+ \int_R^\infty r^{p+1}e^{-(B-\mathcal Im(\kappa))r}\ dr
+ \end{equation}
+ and
+ \begin{equation}
+ \sigma_2:=8\pi\sum_{p=0}^n\frac{n!}{p!|\kappa|^{n-p+1}}
+ \int_0^R r^{p+1}e^{\mathcal Im(\kappa)r}\mathcal V(r)\ dr
+ .
+ \end{equation}
+ Furthermore,
+ \begin{equation}
+ \sigma_1
+ =
+ 8A\pi n!
+ \sum_{p=0}^n\frac{p+1}{(B-\mathcal Im(\kappa))^{p+2}|\kappa|^{n-p+1}}
+ \end{equation}
+ so, as long as $|\kappa|\geqslant \frac12B$ and $\mathcal Im(\kappa)\leqslant\frac12B$,
+ \begin{equation}
+ \sigma_1
+ \leqslant
+ \frac{2^{n+6}A\pi}{B^{n+3}}n!\sum_{p=0}^n(p+1)
+ =
+ \frac{2^{n+5}A\pi}{B^{n+3}}(n+2)!
+ .
+ \end{equation}
+ In addition,
+ \begin{equation}
+ \sigma_2
+ \leqslant
+ 8\pi\sum_{p=0}^n\frac{n!}{p!|\kappa|^{n-p+1}}R^{p-1}e^{\mathcal Im(\kappa)R}\int_0^R r^2\mathcal V(r)\ dr
+ \end{equation}
+ so
+ \begin{equation}
+ \sigma_2
+ \leqslant
+ 8\pi\sum_{p=0}^n\frac{n!2^{n-p+1}}{p!B^{n-p+1}}R^{p-1}e^{\mathcal Im(\kappa)R}\int_0^R r^2\mathcal V(r)\ dr
+ \leqslant
+ \frac{2^{n+4}\pi}{RB^{n+1}} n!e^{RB}\int_0^R r^2\mathcal V(r)\ dr
+ \end{equation}
+ which implies\-~(\ref{decaydS}) in this case.
+ \bigskip
+
+ \subpoint We have thus proved that $S$ is analytic in $H_\tau$, which implies that the singularities of $\widehat{\mathcal U}_2$ in $H_\tau$ all come from the branch points of $\sqrt F$.
+ For $\kappa\in\mathbb R$,
+ \begin{equation}
+ |S(\kappa)|\leqslant 1
+ \end{equation}
+ so, for $\kappa\in\mathbb R$,
+ \begin{equation}
+ F(\kappa)\geqslant \frac{\kappa^2}{2e}
+ .
+ \end{equation}
+ Therefore, since $F$ is analytic in a finite strip around the real axis, $F$ cannot have any roots in the vicinity of the real axis, except at $0$, so the only branch point of $\sqrt F$ near the real axis is $0$.
+ Thus, $\widehat{\mathcal U}_2$ is analytic in $H_\tau$.
+ \bigskip
+
+ \point{\bf Decay of $\mathcal U_2$}. We deform the integral to the path
+ \begin{equation}
+ \{i\eta y,\ 0<y<|x|^{-\tau}\}\cup\{i\eta |x|^{-\tau}+y,\ y>0\}
+ \end{equation}
+ and find
+ \begin{equation}
+ \int_0^\infty e^{i\eta\kappa|x|}\kappa\widehat{\mathcal U}_2(\kappa)\ d\kappa
+ =
+ I_1+I_2
+ \label{I12}
+ \end{equation}
+ with
+ \begin{equation}
+ I_1:=
+ -\int_0^{|x|^{-\tau}}e^{-y|x|}y\widehat{\mathcal U}_2(i\eta y)\ dy
+ \end{equation}
+ and
+ \begin{equation}
+ I_2:=
+ e^{-|x|^{1-\tau}}
+ \int_0^\infty e^{i\eta y|x|}(i\eta|x|^{-\tau}+y)\widehat{\mathcal U}_2(i\eta|x|^{-\tau}+y)\ dy
+ .
+ \end{equation}
+ \bigskip
+
+ \subpoint We first estimate $I_1$.
+ We expand $S$:
+ \begin{equation}
+ S(\kappa)=1-\beta \kappa^2+O(|\kappa|^4)
+ \end{equation}
+ with $\beta>0$ (since $S$ is analytic and symmetric, and $|S(|k|)|\leqslant 1$).
+ Therefore, $y\mapsto\widehat{\mathcal U}_2(iy)$ is $\mathcal C^2$ for $y\neq0$, and
+ \begin{equation}
+ \widehat{\mathcal U}_2(i\eta y)
+ =
+ \frac1\rho-\frac{i\eta y}\rho\sqrt{\frac1{2e}+\beta}
+ +O(y^2)
+ .
+ \end{equation}
+ Furthermore,
+ \begin{equation}
+ -\int_0^{|x|^{-\tau}}e^{-y|x|}y\ dy
+ =
+ -\frac1{|x|^2}+\frac{1+|x|^{1-\tau}}{|x|^2}e^{-|x|^{1-\tau}}
+ \end{equation}
+ \begin{equation}
+ -\int_0^{|x|^{-\tau}}e^{-y|x|}y^2\ dy
+ =
+ -\frac2{|x|^3}+\frac{1+|x|^{1-\tau}(2+x^{1-\tau})}{|x|^3}e^{-|x|^{1-\tau}}
+ \end{equation}
+ and
+ \begin{equation}
+ -\int_0^{|x|^{-\tau}}e^{-y|x|}y^3\ dy
+ =
+ O(|x|^{-4})
+ \end{equation}
+ \begin{equation}
+ I_1=
+ -\frac1{\rho|x|^2}
+ +\frac{2i\eta}{\rho|x|^3}\sqrt{\frac1{2e}+\beta}
+ +O(|x|^{-4})
+ \end{equation}
+ so
+ \begin{equation}
+ \frac1{4i\pi^2|x|}\sum_{\eta=\pm}\eta I_1=\frac1{\pi^2\rho|x|^4}\sqrt{\frac1{2e}+\beta}
+ +O(|x|^{-5})
+ .
+ \label{ineqI1}
+ \end{equation}
+ \bigskip
+
+ \subpoint We now bound $I_2$.
+ Recall that, for $\kappa\in\mathbb R$, $|S(\kappa)|\leqslant 1$.
+ Recalling\-~(\ref{decaydS}),
+ \begin{equation}
+ |S(\kappa+i\eta|x|^{-\tau})|
+ \leqslant
+ \sum_{n=0}^\infty\frac1{n!}|\partial^nS(\kappa)|^n|x|^{-n\tau}
+ \leqslant
+ \frac1{1-C|x|^{-\tau}}
+ \leqslant 2
+ \end{equation}
+ provided $|x|^\tau>2C$.
+ Therefore, for large $\kappa$,
+ \begin{equation}
+ |\widehat{\mathcal U}_2(\kappa+i\eta)|=O(\kappa^{-4})
+ \end{equation}
+ so
+ \begin{equation}
+ I_2\leqslant C'e^{-|x|^{1-\tau}}
+ \label{ineqI2}
+ \end{equation}
+ for some constant $C'>0$.
+ \bigskip
+
+ \subpoint Inserting\-~(\ref{ineqI1}) and\-~(\ref{ineqI2}) into\-~(\ref{I12}) and\-~(\ref{U2}), we find that
+ \begin{equation}
+ \mathcal U_2(|x|)=\frac1{\pi^2\rho|x|^4}\sqrt{\frac1{2e}+\beta}+O(|x|^{-5})
+ \end{equation}
+ which, using\-~(\ref{U1}), concludes the proof of the lemma.
+ \bigskip
+\qed
+
+
+\section{Comparison with the Bose gas}\label{sec:bosegas}
+\subsection{Sketch of the derivation of the simple equation}
+\indent The simple equation\-~(\ref{simpleq})-(\ref{energy}) was originally derived\-~\cite{Li63} to approximate the ground state energy $E_0$ of a repulsive Bose gas, which is a system of $N$ quantum particles interacting via the repulsive potential $\mathcal V$.
+The ground state energy of this system is the lowest eigenvalue of the Hamiltonian operator
+\begin{equation}
+ H_N:=-\frac12\sum_{i=1}^N\Delta_i+\sum_{1\leqslant i<j\leqslant N}\mathcal V(x_i-x_j)
+\end{equation}
+acting on the space of $L_2$ functions on the torus $\mathbb T_V$ of volume $V$.
+The corresponding eigenfunction, which we will denote by $\psi_N$, satisfies
+\begin{equation}
+ H_N\psi_N(x_1,\cdots,x_N)=E_0\psi_N(x_1,\cdots,x_N)
+ \label{eigenvalue}
+\end{equation}
+with $x_i\in\mathbb T_V$.
+As is well known, by a Perron-Frobenius argument, $\psi_N$ is unique, non-negative, and hence symmetric under exchanges $x_i\leftrightarrow x_j$, and under translations.
+\bigskip
+
+\indent We can write $E_0$ by integrating both sides of\-~(\ref{eigenvalue}):
+\begin{equation}
+ E_0=\frac {N(N-1)}{2V}\int g_N^{(2)}(x)\mathcal V(x)\ dx
+ \label{energyg}
+\end{equation}
+with
+\begin{equation}
+ g_N^{(p)}(x_1,\cdots,x_p):=\frac{V^j}{\int\psi_N(x_1,\cdots,x_N)\ dx_1\cdots dx_N}\int\psi_N(x_1,\cdots,x_N)\ dx_{p+1}\cdots dx_N
+ \label{g}
+\end{equation}
+and $g_N^{(2)}(x_1,x_2)\equiv g_N^{(2)}(x_1-x_2)$.
+The computation of $E_0$ thus reduces to that of $g_N^{(2)}$.
+Note that the kinetic energy does not appear explicitly in\-~(\ref{energyg}).
+\bigskip
+
+\indent To compute $g_N^{(2)}$, integrate both sides of\-~(\ref{eigenvalue}) with respect to $x_3,\cdots,x_N$.
+This yields an equation relating $g_N^{(2)}$, $g_N^{(3)}$ and $g_N^{(4)}$.
+The main approximation made in\-~\cite{Li63}, is to write $g_N^{(3)}$ and $g_N^{(4)}$ products of $g_N^{(2)}$ factors: roughly,
+\begin{equation}
+ g_N^{(p)}(x_1,\cdots,x_p)\approx\prod_{1\leqslant i<j\leqslant p}g_N^{(2)}(x_i-x_j)
+ .
+ \label{approxprod}
+\end{equation}
+This is a sensible approximation in the case of low density $\rho=\frac NV\ll 1$.
+Indeed, in this regime, one might expect $\psi_N$ to be approximately a {\it Bijl-Dingle-Jastrow function}:
+\begin{equation}
+ \psi_N(x_1,\cdots,x_N)\approx\prod_{1\leqslant i<j\leqslant N}e^{-\phi(x_i-x_j)}
+\end{equation}
+for some appropriately chosen real function $\phi$.
+Thus, $\psi_N$ is approximated by the {\it partition function} of a classical statistical mechanical model of particles interacting via the pair-potential $\phi$.
+In this setting, $g_N^{(p)}$ is the {\it $p$-point correlation function} of the canonical Gibbs distribution of this model.
+When\-~(\ref{approxprod}) holds asymptotically as the particles move away from each other (remember, the density is low), the statistical mechanics system is said to satisfy the {\it clustering property}.
+There is a long literature on proving the clustering property for a large class of potentials $\phi$, see, among many others, \cite{Ru99,Ga99,PT12}.
+\bigskip
+
+\indent Assuming the clustering property for the potential $\phi$, the assumption\-~(\ref{approxprod}) does not seem far fetched.
+This product structure leads to an equation for $g_N^{(2)}$.
+At this stage, one takes the thermodynamic limit: $N\to\infty$ and $\rho=\frac NV$ fixed.
+There are some subtleties to taking this limit, which are explained in\-~\cite{Li63}.
+Defining $u:=1-g_\infty^{(2)}$, the equation for $u$ is\-~\cite[(3.29)]{Li63}.
+After a few extra reasonable approximations, this equation reduces to\-~(\ref{simpleq}).
+The equation for the energy\-~(\ref{energy}) is simply the $N\to\infty$ limit of\-~(\ref{energyg}).
+
+\subsection{Numerical comparison}\label{subsec:numerics}
+\indent One of the motivations for studying the simple equation is that it provides a simple tool to approximate the ground state energy of the Bose gas.
+In\-~\cite{LL64}, it was found that in one dimension the simple equation gives a value for the energy that differs from the Bose gas ground state energy by at most 69\% (a more complete form of the equation yields even better, result, with a maximal error of 19\%).
+In one dimension, the difference is larger at high density.
+\bigskip
+
+\indent In three dimensions, by Theorem\-~\ref{theorem:asymptotics}, the simple equation predicts the correct low density asymptote as the Bose gas.
+This is a not so surprising, since the derivation of the simple equation from the ground state equation of the Bose gas sketched above seems somewhat sensible when the density is low.
+However, when the density is high, at least in the case in which the potential has a non-negative Fourier transform, the simple equation also yields the same asymptote as the Bose gas.
+In fact, considering the case
+\begin{equation}
+ \mathcal V(x)=e^{-|x|}
+\end{equation}
+(which has a positive Fourier transform), we compared the ground state energy of the simple equation with values from a Monte Carlo simulation of the Bose gas computed by M.\-~Holzmann\-~\cite{CHe}, to whom we are most grateful for sharing his unpublished work.
+The comparison is in figure\-~\ref{fig:bosegas}, in which we found that the maximal error made by the simple equation, over the {\it entire range of densities}, is 5\%!
+This is a promising result, which we will investigate in more depth and with more rigor in a later publication.
+\bigskip
+
+\begin{figure}
+ \hfil\includegraphics[width=10cm]{bosegas.pdf}
+ \caption{
+ Plot of $\frac e{4\pi\rho}$ as a function of $\rho$ on a log scale.
+ The potential is $\mathcal V(r)=e^{-r}$, in which case the scattering length is $a\approx 1.25$.
+ The solid curve is the energy computed from the simple equation\-~(\ref{simpleq})-(\ref{energy}), and the discrete points are the values of the energy of the Bose gas computed by M.\-~Holzmann\-~\cite{CHe} using a Monte Carlo algorithm.
+ The gray area corresponds to a 5\% error on the value of the energy.
+ At low densities, we recover the Lenz asymptote $\frac e{4\pi\rho}\sim\frac a2$ and at high densities, we recover $\frac e{4\pi\rho}\sim1$.
+ The difference between the Monte Carlo simulation and the solution of the simple equation is smaller than 5\%.
+ }
+ \label{fig:bosegas}
+\end{figure}
+
+
+\section{Open problems and conjectures}\label{sec:open}
+
+\point{\bf Monotonicity}.
+An important open problem is to show that $e\mapsto\rho(e)$ is an increasing function.
+If the solution of the simple equation is in any way related to the ground state wave function of the Bose gas, then this should hold: if the density increases, the energy should increase.
+In addition, it would enable us to prove the uniqueness of the solution of the simple equation with fixed $\rho$, and might even allow us to generalize our result to potentials with hard core components, as well as to relax the constraint that $\mathcal V$ decays exponentially in Theorem\-~\ref{theorem:decay}.
+By running a few numerical computations, it seems clear that $\rho(e)$ should be increasing, see figure\-~\ref{fig:bosegas}.
+Using a modified iteration in which $\rho$ is fixed, we have proved that $e\rho(e)$ is strictly monotone increasing in $e$, but the proof that $\rho(e)$ is as well has eluded us thus far.
+\bigskip
+
+\point{\bf Convexity}.
+Another open problem is to prove that $\rho e(\rho)$ {\it is a convex function}, or, equivalently, that $\frac1{\rho(e)}$ is concave.
+In a physical setting, one expects $\rho e(\rho)$ to be convex.
+Indeed if $\rho e=:e_{\mathrm v}$ were {\it not convex}, there would exist $\rho_1<\rho<\rho_2$ such that $\frac{\rho_1+\rho_2}2=\rho$ and $e_{\mathrm v}(\rho_1)+e_{\mathrm v}(\rho_2)<2e_{\mathrm v}(\rho)$.
+Furthermore, $e_{\mathrm v}$ is the energy per unit volume, and, considering a volume $V$ that is split into two equal halves, we find that a configuration in which one half of the volume holds a density $\rho_1$ of particles, whereas the other holds $\rho_2$ would have energy
+\begin{equation}
+ \frac V2(e_{\mathrm v}(\rho_1)+e_{\mathrm v}(\rho_2))<Ve_{\mathrm v}(\rho)
+ .
+\end{equation}
+Therefore, it would pay to have more particles in one half than in the other, which is unstable.
+Numerically, it seems quite clear that $\rho e(\rho)$ is convex, see figure\-~\ref{fig:convexity}.
+\bigskip
+
+\begin{figure}
+ \hfil\includegraphics[width=10cm]{convexity.pdf}
+ \caption{
+ Numerical evaluation of $\frac1{4\pi}\partial_\rho^2(\rho e)$ for $\mathcal V(r)=e^{-r}$.
+ The asymptotic values are, for $\rho\to0$, $a\approx 1.25$ and for $\rho\to\infty$, $2$.
+ This second derivative seems to be clearly positive, so $\rho e$ appears to be convex.
+ }
+ \label{fig:convexity}
+\end{figure}
+
+\point{\bf Solutions with negative values}.
+In this paper, we solved the simple equation for functions $u$ that satisfy\-~(\ref{con1}).
+The condition that $u(x)\leqslant 1$ comes from physical considerations, and we are gratified that our simple equation has this property automatically, see\-~(\ref{con1}).
+The correlation function $g_N^{(2)}$ defined in\-~(\ref{g}) is non-negative, which means that $u(x)\leqslant 1$.
+However, one may wonder whether the condition $u(x)\geqslant0$ must be imposed, or whether it may follow from the simple equation.
+In principle, Theorem\-~\ref{theorem:existence} does not exclude the existence of other solutions of\-~(\ref{simpleq})-(\ref{energy}) in which $u(x)<0$ for some $x\in\mathbb R^d$.
+Proving that there are no such solutions is another interesting open problem.
+It seems rather unlikely that such solutions exist: defining
+\begin{equation}
+ \omega(x):=G_eu(x)
+\end{equation}
+(\ref{simpleq}) becomes
+\begin{equation}
+ \omega(x)=
+ G_e^2(1-u)\mathcal V(x)+2e\rho\omega\ast\omega(x)
+\end{equation}
+but, $u(x)\leqslant 1$, so $G_e^2(1-u)\mathcal V(x)\geqslant 0$ and
+\begin{equation}
+ \omega(x)\geqslant 2e\rho\omega\ast\omega(x)
+ .
+\end{equation}
+In particular,
+\begin{equation}
+ \{x:\ \omega(x)<0\}\subset
+ \{x:\ \omega\ast\omega(x)<0\}
+\end{equation}
+which does not seem to be possible, although a proof that it is not so has eluded us.
+\bigskip
+
+\point {\bf Solution of the full equation}.
+The simple equation\-~(\ref{simpleq}) is actually a simplified version of an equation that should approximate the Bose gas more accurately\-~\cite{Li63}:
+\begin{equation}
+  (-\Delta+\mathcal V(x))u(x)
+  =\mathcal V(x)
+  -\rho(1-u(x))(2K(x)-\rho L(x))
+  \label{fulleq}
+\end{equation}
+with
+\begin{equation}
+  K(x):=u\ast S(x)
+  ,\quad
+  S(x):=(1-u(x))\mathcal V(x)
+\end{equation}
+\begin{equation}
+  L(x):=
+  \int u(y)u(z-x)\left(
+    1-u(z)-u(y-x)+\frac12u(z)u(y-x)
+    \ dydz
+  \right)S(y)
+  .
+  \label{L}
+\end{equation}
+Note that $e$ appears only as the integral of $S$, see\-~(\ref{energy}).
+Little is known about the solutions of this equation: even proving the existence of a solution of\-~(\ref{fulleq}) is open.
+We hope to study this equation numerically in a later publication.
+
+\vfill
+
+{\bf Acknowledgements}:
+  {\it
+  Thanks go to Markus Holzmann for his hard work and skill computing the ground state energy of the Bose gas, and for sharing his results with us.
+  E.C. gratefully acknowledges support through NSF grant DMS-174625.
+  I.J. gratefully acknowledges support through NSF grant DMS-1802170.
+  }
+\eject
+
+\begin{thebibliography}{WWW99}
+\small
+ 
+\bibitem[Bo47]{Bo47}N. Bogolubov - {\it On the theory of superfluidity}, Journal of Physics (USSR), volume\-~11, number\-~1 , pages\-~23-32 (translated from the Russian Izv.Akad.Nauk Ser.Fiz, volume\-~11, pages\-~77-90), 1947.\par\medskip
+ 
+\bibitem[CHe]{CHe}E.\-~Carlen, M.\-~Holzmann, I.\-~Jauslin, E.H.\-~Lieb, in preparation.\par\medskip
+ 
+\bibitem[Dy57]{Dy57}F.J. Dyson - {\it Ground-State Energy of a Hard-Sphere Gas}, Physical Review, volume\-~106, issue\-~1, pages\-~20-26, 1957,\par\penalty10000
+doi:{\tt\color{blue}\href{http://dx.doi.org/10.1103/PhysRev.106.20}{10.1103/PhysRev.106.20}}.\par\medskip
+ 
+\bibitem[FS19]{FS19}S.\-~Fournais, J.P.\-~Solovej - {\it The energy of dilute Bose gases}, 2019,\par\penalty10000
+arxiv:{\tt\href{http://arxiv.org/abs/1904.06164}{1904.06164}}.\par\medskip
+ 
+\bibitem[Ga99]{Ga99}G. Gallavotti - {\it Statistical mechanics, a short treatise}, Springer, 1999.\par\medskip
+ 
+\bibitem[LHY57]{LHY57}T.D. Lee, K. Huang, C.N. Yang - {\it Eigenvalues and Eigenfunctions of a Bose System of Hard Spheres and Its Low-Temperature Properties}, Physical Review, volume\-~106, issue\-~6, pages\-~1135-1145, 1957,\par\penalty10000
+doi:{\tt\color{blue}\href{http://dx.doi.org/10.1103/PhysRev.106.1135}{10.1103/PhysRev.106.1135}}.\par\medskip
+ 
+\bibitem[Le29]{Le29}W. Lenz - {\it Die Wellenfunktion und Geschwindigkeitsverteilung des entarteten Gases}, Zeitschrift f\"ur Physik, volume\-~56, issue\-~11-12, pages\-~778-789, 1929,\par\penalty10000
+doi:{\tt\color{blue}\href{http://dx.doi.org/10.1007/BF01340138}{10.1007/BF01340138}}.\par\medskip
+ 
+\bibitem[Li63]{Li63}E.H. Lieb - {\it Simplified Approach to the Ground-State Energy of an Imperfect Bose Gas}, Physical Review, volume\-~130, issue\-~6, pages\-~2518-2528, 1963,\par\penalty10000
+doi:{\tt\color{blue}\href{http://dx.doi.org/10.1103/PhysRev.130.2518}{10.1103/PhysRev.130.2518}}.\par\medskip
+ 
+\bibitem[LL64]{LL64}E.H. Lieb, W. Liniger - {\it Simplified Approach to the Ground-State Energy of an Imperfect Bose Gas. III. Application to the One-Dimensional Model}, Physical Review, volume\-~134, issue\-~2A, pages A312-A315, 1964,\par\penalty10000
+doi:{\tt\color{blue}\href{http://dx.doi.org/10.1103/PhysRev.134.A312}{10.1103/PhysRev.134.A312}}.\par\medskip
+ 
+\bibitem[LL01]{LL01}E.H. Lieb, M. Loss - {\it Analysis}, Second edition, Graduate studies in mathematics, Americal Mathematical Society, 2001.\par\medskip
+ 
+\bibitem[LY98]{LY98}E.H. Lieb, J. Yngvason - {\it Ground State Energy of the Low Density Bose Gas}, Physical Review Letters, volume\-~80, issue\-~12, pages\-~2504-2507, 1998,\par\penalty10000
+doi:{\tt\color{blue}\href{http://dx.doi.org/10.1103/PhysRevLett.80.2504}{10.1103/PhysRevLett.80.2504}}, arxiv:{\tt\color{blue}\href{http://arxiv.org/abs/cond-mat/9712138}{cond-mat/9712138}}.\par\medskip
+ 
+\bibitem[LY01]{LY01}E.H. Lieb, J. Yngvason - {\it The Ground State Energy of a Dilute Two-Dimensional Bose Gas}, Journal of Statistical Physics, volume\-~103, issue\-~3-4, pages\-~509-526, 2001,\par\penalty10000
+doi:{\tt\color{blue}\href{http://dx.doi.org/10.1023/A:1010337215241}{10.1023/A:1010337215241}}, arxiv:{\tt\color{blue}\href{http://arxiv.org/abs/math-ph/0002014}{math-ph/0002014}}.\par\medskip
+ 
+\bibitem[PT12]{PT12}E. Pulvirenti, D. Tsagkarogiannis - {\it Cluster Expansion in the Canonical Ensemble}, Communications in Mathematical Physics, volume\-~316, issue\-~2, pages\-~289-306, 2012,\par\penalty10000
+doi:{\tt\color{blue}\href{http://dx.doi.org/10.1007/s00220-012-1576-y}{10.1007/s00220-012-1576-y}}, arxiv:{\tt\color{blue}\href{http://arxiv.org/abs/1105.1022}{1105.1022}}.\par\medskip
+ 
+\bibitem[RS75]{RS75b}M. Reed, B. Simon - {\it Methods of Modern Mathematical Physics II: Fourier Analysis, Self-Adjointness}, second edition, Academic Press, New York, 1975.\par\medskip
+ 
+\bibitem[Ru99]{Ru99}D. Ruelle - {\it Statistical mechanics: rigorous results}, Imperial College Press, World Scientific, (first edition: Benjamin, 1969), 1999.\par\medskip
+ 
+\bibitem[YY09]{YY09}H. Yau, J. Yin - {\it The Second Order Upper Bound for the Ground Energy of a Bose Gas}, Journal of Statistical Physics, volume\-~136, issue\-~3, pages\-~453-503, 2009,\par\penalty10000
+doi:{\tt\color{blue}\href{http://dx.doi.org/10.1007/s10955-009-9792-3}{10.1007/s10955-009-9792-3}}, arxiv:{\tt\color{blue}\href{http://arxiv.org/abs/0903.5347}{0903.5347}}.\par\medskip
+ 
+\end{thebibliography}
+
+
+
+\end{document}
diff --git a/Makefile b/Makefile
new file mode 100644
index 0000000..77c40ae
--- /dev/null
+++ b/Makefile
@@ -0,0 +1,54 @@
+PROJECTNAME=$(basename $(wildcard *.tex))
+LIBS=$(notdir $(wildcard libs/*))
+FIGS=$(notdir $(wildcard figs/*.fig))
+
+PDFS=$(addsuffix .pdf, $(PROJECTNAME))
+SYNCTEXS=$(addsuffix .synctex.gz, $(PROJECTNAME))
+
+all: $(PROJECTNAME)
+
+$(PROJECTNAME): $(LIBS) $(FIGS)
+	pdflatex -file-line-error $@.tex
+	pdflatex -file-line-error $@.tex
+	pdflatex -synctex=1 $@.tex
+
+$(PROJECTNAME).aux: $(LIBS) $(FIGS)
+	pdflatex -file-line-error -draftmode $(PROJECTNAME).tex
+
+
+$(SYNCTEXS): $(LIBS) $(FIGS)
+	pdflatex -synctex=1 $(patsubst %.synctex.gz, %.tex, $@)
+
+libs: $(LIBS)
+
+$(LIBS):
+	ln -fs libs/$@ ./
+
+figs: $(FIGS)
+
+$(FIGS):
+	make -C figs/$@
+	for pdf in $$(find figs/$@/ -name '*.pdf'); do ln -fs "$$pdf" ./ ; done
+
+
+clean-aux: clean-figs-aux
+	rm -f $(addsuffix .aux, $(PROJECTNAME))
+	rm -f $(addsuffix .log, $(PROJECTNAME))
+	rm -f $(addsuffix .out, $(PROJECTNAME))
+	rm -f $(addsuffix .toc, $(PROJECTNAME))
+
+clean-libs:
+	rm -f $(LIBS)
+
+clean-figs:
+	$(foreach fig,$(addprefix figs/, $(FIGS)), make -C $(fig) clean; )
+	rm -f $(notdir $(wildcard figs/*.fig/*.pdf))
+
+clean-figs-aux:
+	$(foreach fig,$(addprefix figs/, $(FIGS)), make -C $(fig) clean-aux; )
+
+
+clean-tex:
+	rm -f $(PDFS) $(SYNCTEXS)
+
+clean: clean-aux clean-tex clean-libs clean-figs
diff --git a/README b/README
new file mode 100644
index 0000000..cffdc99
--- /dev/null
+++ b/README
@@ -0,0 +1,47 @@
+This directory contains the source files to typeset the article, and generate
+the figures. This can be accomplished by running
+  make
+
+A program is bundled with this document to compute solutions to the simple
+equation numerically. It is located in figs/numerical.figs/simpleq. We will not
+give a full documentation of the program, which will be provided with a more
+complete version at some later time. The program is run when running 'make' to
+compute the data for the figures.
+
+This document uses a custom class file, located in the 'libs' directory, which
+defines a number of commands. Most of these are drop-in replacements for those
+defined in the 'article' class.
+
+Some extra functionality is provided in custom style files, located in the
+'libs' directory.
+
+
+* Dependencies:
+
+  pdflatex
+  TeXlive packages:
+    amsfonts
+    array
+    color
+    dsfont
+    etoolbox
+    graphics
+    hyperref
+    latex
+    marginnote
+    pgf
+    standalone
+  GNU make
+  julia
+
+* Files:
+
+  Carlen_Jauslin_Lieb_2019.tex:
+    main LaTeX file
+
+  libs:
+    custom LaTeX class file
+
+  figs:
+    source code for the figures
+
diff --git a/bibliography/bibliography.tex b/bibliography/bibliography.tex
new file mode 100644
index 0000000..7236c3e
--- /dev/null
+++ b/bibliography/bibliography.tex
@@ -0,0 +1,42 @@
+\bibitem[Bo47]{Bo47}N. Bogolubov - {\it On the theory of superfluidity}, Journal of Physics (USSR), volume\-~11, number\-~1 , pages\-~23-32 (translated from the Russian Izv.Akad.Nauk Ser.Fiz, volume\-~11, pages\-~77-90), 1947.\par\medskip
+ 
+\bibitem[CHe]{CHe}E.\-~Carlen, M.\-~Holzmann, I.\-~Jauslin, E.H.\-~Lieb, in preparation.\par\medskip
+ 
+\bibitem[Dy57]{Dy57}F.J. Dyson - {\it Ground-State Energy of a Hard-Sphere Gas}, Physical Review, volume\-~106, issue\-~1, pages\-~20-26, 1957,\par\penalty10000
+doi:{\tt\href{http://dx.doi.org/10.1103/PhysRev.106.20}{10.1103/PhysRev.106.20}}.\par\medskip
+ 
+\bibitem[FS19]{FS19}S.\-~Fournais, J.P.\-~Solovej - {\it The energy of dilute Bose gases}, 2019,\par\penalty10000
+arxiv:{\tt\href{http://arxiv.org/abs/1904.06164}{1904.06164}}.\par\medskip
+ 
+\bibitem[Ga99]{Ga99}G. Gallavotti - {\it Statistical mechanics, a short treatise}, Springer, 1999.\par\medskip
+ 
+\bibitem[LHY57]{LHY57}T.D. Lee, K. Huang, C.N. Yang - {\it Eigenvalues and Eigenfunctions of a Bose System of Hard Spheres and Its Low-Temperature Properties}, Physical Review, volume\-~106, issue\-~6, pages\-~1135-1145, 1957,\par\penalty10000
+doi:{\tt\href{http://dx.doi.org/10.1103/PhysRev.106.1135}{10.1103/PhysRev.106.1135}}.\par\medskip
+ 
+\bibitem[Le29]{Le29}W. Lenz - {\it Die Wellenfunktion und Geschwindigkeitsverteilung des entarteten Gases}, Zeitschrift f\"ur Physik, volume\-~56, issue\-~11-12, pages\-~778-789, 1929,\par\penalty10000
+doi:{\tt\href{http://dx.doi.org/10.1007/BF01340138}{10.1007/BF01340138}}.\par\medskip
+ 
+\bibitem[Li63]{Li63}E.H. Lieb - {\it Simplified Approach to the Ground-State Energy of an Imperfect Bose Gas}, Physical Review, volume\-~130, issue\-~6, pages\-~2518-2528, 1963,\par\penalty10000
+doi:{\tt\href{http://dx.doi.org/10.1103/PhysRev.130.2518}{10.1103/PhysRev.130.2518}}.\par\medskip
+ 
+\bibitem[LL64]{LL64}E.H. Lieb, W. Liniger - {\it Simplified Approach to the Ground-State Energy of an Imperfect Bose Gas. III. Application to the One-Dimensional Model}, Physical Review, volume\-~134, issue\-~2A, pages A312-A315, 1964,\par\penalty10000
+doi:{\tt\href{http://dx.doi.org/10.1103/PhysRev.134.A312}{10.1103/PhysRev.134.A312}}.\par\medskip
+ 
+\bibitem[LL01]{LL01}E.H. Lieb, M. Loss - {\it Analysis}, Second edition, Graduate studies in mathematics, Americal Mathematical Society, 2001.\par\medskip
+ 
+\bibitem[LY98]{LY98}E.H. Lieb, J. Yngvason - {\it Ground State Energy of the Low Density Bose Gas}, Physical Review Letters, volume\-~80, issue\-~12, pages\-~2504-2507, 1998,\par\penalty10000
+doi:{\tt\href{http://dx.doi.org/10.1103/PhysRevLett.80.2504}{10.1103/PhysRevLett.80.2504}}, arxiv:{\tt\href{http://arxiv.org/abs/cond-mat/9712138}{cond-mat/9712138}}.\par\medskip
+ 
+\bibitem[LY01]{LY01}E.H. Lieb, J. Yngvason - {\it The Ground State Energy of a Dilute Two-Dimensional Bose Gas}, Journal of Statistical Physics, volume\-~103, issue\-~3-4, pages\-~509-526, 2001,\par\penalty10000
+doi:{\tt\href{http://dx.doi.org/10.1023/A:1010337215241}{10.1023/A:1010337215241}}, arxiv:{\tt\href{http://arxiv.org/abs/math-ph/0002014}{math-ph/0002014}}.\par\medskip
+ 
+\bibitem[PT12]{PT12}E. Pulvirenti, D. Tsagkarogiannis - {\it Cluster Expansion in the Canonical Ensemble}, Communications in Mathematical Physics, volume\-~316, issue\-~2, pages\-~289-306, 2012,\par\penalty10000
+doi:{\tt\href{http://dx.doi.org/10.1007/s00220-012-1576-y}{10.1007/s00220-012-1576-y}}, arxiv:{\tt\href{http://arxiv.org/abs/1105.1022}{1105.1022}}.\par\medskip
+ 
+\bibitem[RS75]{RS75b}M. Reed, B. Simon - {\it Methods of Modern Mathematical Physics II: Fourier Analysis, Self-Adjointness}, second edition, Academic Press, New York, 1975.\par\medskip
+ 
+\bibitem[Ru99]{Ru99}D. Ruelle - {\it Statistical mechanics: rigorous results}, Imperial College Press, World Scientific, (first edition: Benjamin, 1969), 1999.\par\medskip
+ 
+\bibitem[YY09]{YY09}H. Yau, J. Yin - {\it The Second Order Upper Bound for the Ground Energy of a Bose Gas}, Journal of Statistical Physics, volume\-~136, issue\-~3, pages\-~453-503, 2009,\par\penalty10000
+doi:{\tt\href{http://dx.doi.org/10.1007/s10955-009-9792-3}{10.1007/s10955-009-9792-3}}, arxiv:{\tt\href{http://arxiv.org/abs/0903.5347}{0903.5347}}.\par\medskip
+ 
diff --git a/figs/numerical.fig/Makefile b/figs/numerical.fig/Makefile
new file mode 100644
index 0000000..8fb4411
--- /dev/null
+++ b/figs/numerical.fig/Makefile
@@ -0,0 +1,32 @@
+PROJECTNAME=bosegas convexity
+
+DATS=erho.dat ddrhoe.dat
+PDFS=$(addsuffix .pdf, $(PROJECTNAME))
+TEXS=$(addsuffix .tikz.tex, $(PROJECTNAME))
+
+all: $(PDFS)
+
+$(PDFS): $(DATS)
+	gnuplot $(patsubst %.pdf, %.gnuplot, $@) > $(patsubst %.pdf, %.tikz.tex, $@)
+	pdflatex -jobname $(basename $@) -file-line-error $(patsubst %.pdf, %.tikz.tex, $@)
+
+erho.dat:
+	julia ./simpleq/main.jl -p "tolerance=1e-14;order=100" energy_rho > $@
+ddrhoe.dat:
+	julia ./simpleq/main.jl -p "tolerance=1e-14;order=100" convexity > $@
+
+install: $(PDFS)
+	cp $^ $(INSTALLDIR)/
+
+clean-aux:
+	rm -f $(addsuffix .tikz.tex, $(PROJECTNAME))
+	rm -f $(addsuffix .aux, $(PROJECTNAME))
+	rm -f $(addsuffix .log, $(PROJECTNAME))
+
+clean-dat:
+	rm -f $(DATS)
+
+clean-tex:
+	rm -f $(PDFS)
+
+clean: clean-aux clean-tex
diff --git a/figs/numerical.fig/bosegas.gnuplot b/figs/numerical.fig/bosegas.gnuplot
new file mode 100644
index 0000000..bfcffac
--- /dev/null
+++ b/figs/numerical.fig/bosegas.gnuplot
@@ -0,0 +1,34 @@
+set ylabel "$\\displaystyle\\frac{e}{4\\pi\\rho}$" norotate offset -1,0
+set xlabel "$\\rho$"
+
+set xtics 1e-6, 100, 100
+set xtics add ("$10^{-6}$" 0.000001, "$10^{-4}$" 0.0001, "$10^{-2}$" 0.01, "$1$" 1.0, "$10^2$" 100)
+unset mxtics
+
+set ytics 0.6, 0.1
+set mytics 2
+
+set xrange [0.000001:100]
+set yrange [0.6:1.05]
+
+# default output canvas size: 12.5cm x 8.75cm
+set term lua tikz size 8,6 standalone
+
+set key off
+
+
+# set linestyle
+set style line 1 linetype rgbcolor "#4169E1" linewidth 3
+set style line 2 linetype rgbcolor "#DC143C" linewidth 3
+set style line 3 linetype rgbcolor "#32CD32" linewidth 3
+set style line 4 linetype rgbcolor "#4B0082" linewidth 3
+set style line 5 linetype rgbcolor "#DAA520" linewidth 3
+
+set pointsize 1
+
+set logscale x
+
+plot "erho.dat" using 1:(0.95*$2/$1/(4*pi)):(1.05*$2/$1/(4*pi)) with filledcurves linetype rgbcolor "#DDDDDD", \
+ "erho.dat" using 1:($2/$1/(4*pi)) with lines linestyle 1, \
+ "holzmann_2019-09-28.dat" using 1:($2/$1/(4*pi)) with points linestyle 2
+
diff --git a/figs/numerical.fig/convexity.gnuplot b/figs/numerical.fig/convexity.gnuplot
new file mode 100644
index 0000000..f6f606b
--- /dev/null
+++ b/figs/numerical.fig/convexity.gnuplot
@@ -0,0 +1,30 @@
+set ylabel "$\\frac1{4\\pi}\\partial_\\rho^2(e\\rho)$" offset -4,0 norotate
+set xlabel "$\\rho$"
+
+set xtics 1e-6, 100, 100
+set xtics add ("$10^{-6}$" 0.000001, "$10^{-4}$" 0.0001, "$10^{-2}$" 0.01, "$1$" 1.0, "$10^2$" 100)
+unset mxtics
+set xrange [0.000001:100]
+
+set ytics 1.2, 0.2
+set mytics 2
+set yrange [1.2:2.1]
+
+# default output canvas size: 12.5cm x 8.75cm
+set term lua tikz size 8,6 standalone
+
+# no key
+set key off
+
+# set linestyle
+set style line 1 linetype rgbcolor "#4169E1" linewidth 3
+set style line 2 linetype rgbcolor "#DC143C" linewidth 3
+set style line 3 linetype rgbcolor "#000000" linewidth 1
+set style line 4 linetype rgbcolor "#4B0082" linewidth 3
+set style line 5 linetype rgbcolor "#DAA520" linewidth 3
+
+set pointsize 0.6
+
+set logscale x
+
+plot "ddrhoe.dat" using 1:($2/(4*pi)) with lines linestyle 1
diff --git a/figs/numerical.fig/holzmann_2019-09-28.dat b/figs/numerical.fig/holzmann_2019-09-28.dat
new file mode 100644
index 0000000..fe9c0e5
--- /dev/null
+++ b/figs/numerical.fig/holzmann_2019-09-28.dat
@@ -0,0 +1,9 @@
+## data from M. Holzmann, 2019-09-22
+# rho        E0    E0+dE0
+1e-4  8.3500e-4 8.3600e-4
+1e-3  9.1340e-3 9.1350e-3
+1e-2  1.0609e-1 1.0610e-1
+1e-1  1.1930e+0 1.1940e+0
+1e-0  1.2445e+1 1.2446e+1
+1e+1  1.2544e+2 1.2545e+2
+
diff --git a/figs/numerical.fig/simpleq/integration.jl b/figs/numerical.fig/simpleq/integration.jl
new file mode 100644
index 0000000..0ad456a
--- /dev/null
+++ b/figs/numerical.fig/simpleq/integration.jl
@@ -0,0 +1,28 @@
+# approximate \int_a^b f using Gauss-Legendre quadratures
+function integrate_legendre(f,a,b,weights)
+  out=0
+  for i in 1:length(weights[1])
+    out+=(b-a)/2*weights[2][i]*f((b-a)/2*weights[1][i]+(b+a)/2)
+  end
+  return out
+end
+# \int f*g where g is sampled at the Legendre nodes
+function integrate_legendre_sampled(f,g,a,b,weights)
+  out=0
+  for i in 1:length(weights[1])
+    out+=(b-a)/2*weights[2][i]*f((b-a)/2*weights[1][i]+(b+a)/2)*g[i]
+  end
+  return out
+end
+
+
+
+# approximate \int_a^b f/sqrt((b-x)(x-a)) using Gauss-Chebyshev quadratures
+function integrate_chebyshev(f,a,b,N)
+  out=0
+  for i in 1:N
+    out=out+pi/N*f((b-a)/2*cos((2*i-1)/(2*N)*pi)+(b+a)/2)
+  end
+  return out
+end
+
diff --git a/figs/numerical.fig/simpleq/iteration.jl b/figs/numerical.fig/simpleq/iteration.jl
new file mode 100644
index 0000000..f6f044a
--- /dev/null
+++ b/figs/numerical.fig/simpleq/iteration.jl
@@ -0,0 +1,55 @@
+# \hat u_n
+function hatun_iteration(e,order,d,v,maxiter)
+  # gauss legendre weights
+  weights=gausslegendre(order)
+
+  # init V and Eta
+  (V,V0,Eta,Eta0)=init_veta(weights,d,v)
+
+  # init u and rho
+  u=Array{Array{Complex{Float64}}}(undef,maxiter+1)
+  u[1]=zeros(Complex{Float64},order)
+  rho=zeros(Complex{Float64},maxiter+1)
+
+  # iterate
+  for n in 1:maxiter
+    u[n+1]=A(e,weights,Eta)\b(u[n],e,rho[n],V)
+    rho[n+1]=rhon(u[n+1],e,weights,V0,Eta0)
+  end
+
+  return (u,rho)
+end
+
+# A
+function A(e,weights,Eta)
+  N=length(weights[1])
+  out=zeros(Complex{Float64},N,N)
+  for i in 1:N
+    k=(1-weights[1][i])/(1+weights[1][i])
+    out[i,i]=k^2+4*e
+    for j in 1:N
+      y=(weights[1][j]+1)/2
+      out[i,j]+=weights[2][j]*(1-y)*Eta[i][j]/(2*(2*pi)^3*y^3)
+    end
+  end
+  return out
+end
+
+# b
+function b(u,e,rho,V)
+  out=zeros(Complex{Float64},length(V))
+  for i in 1:length(V)
+    out[i]=V[i]+2*e*rho*u[i]^2
+  end
+  return out
+end
+
+# rho_n
+function rhon(u,e,weights,V0,Eta0)
+  S=V0
+  for i in 1:length(weights[1])
+    y=(weights[1][i]+1)/2
+    S+=-weights[2][i]*(1-y)*u[i]*Eta0[i]/(2*(2*pi)^3*y^3)
+  end
+  return 2*e/S
+end
diff --git a/figs/numerical.fig/simpleq/main.jl b/figs/numerical.fig/simpleq/main.jl
new file mode 100644
index 0000000..ece0703
--- /dev/null
+++ b/figs/numerical.fig/simpleq/main.jl
@@ -0,0 +1,192 @@
+using FastGaussQuadrature
+using Printf
+using LinearAlgebra
+include("tools.jl")
+include("integration.jl")
+include("simpleq.jl")
+include("simpleq-newton.jl")
+include("iteration.jl")
+
+
+function main()
+  ## defaults
+
+  tolerance=1e-14
+  order=100
+  maxiter=21
+
+  rho=1e-6
+  e=1e-4
+
+  d=3
+  v=v_exp3d
+  a0=a0_exp3d
+
+  # plot range when plotting in x
+  xmin=0
+  xmax=100
+  nx=100
+
+  # read cli arguments
+  (params,command)=read_args(ARGS)
+
+  # read params
+  if params!=""
+    for param in split(params,";")
+      terms=split(param,"=")
+      if length(terms)!=2
+	print(stderr,"error: could not read parameter '",param,"'.\n")
+	exit(-1)
+      end
+      lhs=terms[1]
+      rhs=terms[2]
+      if lhs=="rho"
+	rho=parse(Float64,rhs)
+      elseif lhs=="e"
+	e=parse(Float64,rhs)
+      elseif lhs=="tolerance"
+	tolerance=parse(Float64,rhs)
+      elseif lhs=="order"
+	order=parse(Int64,rhs)
+      elseif lhs=="maxiter"
+	maxiter=parse(Int64,rhs)
+      elseif lhs=="xmin"
+	xmin=parse(Float64,rhs)
+      elseif lhs=="xmax"
+	xmax=parse(Float64,rhs)
+      elseif lhs=="nx"
+	nx=parse(Int64,rhs)
+      else
+	print(stderr,"unrecognized parameter '",lhs,"'.\n")
+	exit(-1)
+      end
+    end
+  end
+
+  ## run command
+  # e(rho)
+  if command=="energy_rho"
+    energy_rho(order,a0,d,v,maxiter,tolerance)
+  # d^2(rho*e(rho))
+  elseif command=="convexity"
+    ddrhoe(order,a0,d,v,maxiter,tolerance)
+  # u_n(x)
+  elseif command=="iteration"
+    iteration_ux(order,e,a0,d,v,maxiter)
+  else
+    print(stderr,"unrecognized command '",command,"'.\n")
+    exit(-1)
+  end
+
+end
+
+# read cli arguments
+function read_args(ARGS)
+  # flag
+  flag=""
+
+  # output strings
+  params=""
+  command=""
+
+  # loop over arguments
+  for arg in ARGS
+    # argument that starts with a dash
+    if arg[1]=='-'
+      # go through characters after dash
+      for char in arg[2:length(arg)]
+	
+	# set params
+	if char=='p'
+	  # raise flag
+	  flag="params"
+	else
+	  print_usage()
+	  exit(-1)
+	end
+      end
+    # arguments that do not start with a dash
+    else
+      if flag=="params"
+	params=arg
+      else
+	command=arg
+      end
+      # reset flag
+      flag=""
+    end
+  end
+
+  if command==""
+    print_usage()
+    exit(-1)
+  end
+
+  return (params,command)
+end
+
+# usage
+function print_usage()
+  print(stderr,"usage:  simpleq [-p params] <command>\n")
+end
+
+# exponential potential in 3 dimensions
+function v_exp3d(k)
+  return 8*pi/(1+k^2)^2
+end
+a0_exp3d=1.254356410591064819505409291110046864031181245635821944528
+
+# compute energy as a function of rho
+function energy_rho(order,a0,d,v,maxiter,tolerance)
+  minlrho=-6
+  maxlrho=2
+  nlrho=100
+
+  for j in 0:nlrho-1
+    rho=10^(minlrho+(maxlrho-minlrho)*j/nlrho)
+    # linear spacing
+    #rho=10.0^minlrho+(10.0^maxlrho-10.0^minlrho)*j/nlrho
+
+    (u,E)=hatu_newton(order,rho,a0,d,v,maxiter,tolerance)
+
+    @printf("% .8e % .8e\n",rho,real(E))
+
+  end
+end
+
+# compute \partial^2(e\rho) as a function of rho
+function ddrhoe(order,a0,d,v,maxiter,tolerance)
+  minlrho=-6
+  maxlrho=2
+  nlrho=100
+
+  for j in 0:nlrho-1
+    rho=10^(minlrho+(maxlrho-minlrho)*j/nlrho)
+
+    # interval
+    drho=rho*1.01
+
+    (u,E)=hatu_newton(order,rho,a0,d,v,maxiter,tolerance)
+    (up,Ep)=hatu_newton(order,rho+drho,a0,d,v,maxiter,tolerance)
+    (um,Em)=hatu_newton(order,rho-drho,a0,d,v,maxiter,tolerance)
+
+    @printf("% .8e % .8e\n",rho,real((rho+drho)*Ep+(rho-drho)*Em-2*rho*E)/drho^2)
+
+  end
+end
+
+# compute \int u_n(x) at every step
+function iteration_ux(order,e,a0,d,v,maxiter)
+  (u,rho)=hatun_iteration(e,order,d,v,maxiter)
+
+  weights=gausslegendre(order)
+
+  intun=0.
+  for n in 2:maxiter+1
+    # compute \hat u_n(0)=1/(2*rho_n)+rho_{n-1}/2*\hat u_{n-1}(0)^2
+    intun=real(1/(2*rho[n])+rho[n-1]/2*intun^2)
+    @printf("%2d % .15e\n",n-1,abs(intun-1/rho[maxiter+1]))
+  end
+end
+
+main()
diff --git a/figs/numerical.fig/simpleq/simpleq-newton.jl b/figs/numerical.fig/simpleq/simpleq-newton.jl
new file mode 100644
index 0000000..2c70162
--- /dev/null
+++ b/figs/numerical.fig/simpleq/simpleq-newton.jl
@@ -0,0 +1,95 @@
+# \hat u(k) computed using Newton algorithm
+function hatu_newton(order,rho,a0,d,v,maxiter,tolerance)
+
+  # compute gaussian quadrature weights
+  weights=gausslegendre(order)
+
+  # initialize V and Eta
+  (V,V0,Eta,Eta0)=init_veta(weights,d,v)
+
+  # initialize u, V and Eta
+  u=zeros(Complex{Float64},order)
+  for j in 1:order
+    # transformed k
+    k=(1-weights[1][j])/(1+weights[1][j])
+    if d==2
+      u[j]=2*pi*a0*rho/k
+    elseif d==3
+      u[j]=4*pi*a0*rho/k^2
+    end
+  end
+
+  # iterate
+  for i in 1:maxiter-1
+    new=u-inv(dXi(u,V,V0,Eta,Eta0,weights,rho,d))*Xi(u,V,V0,Eta,Eta0,weights,rho,d)
+
+    if(norm(new-u)/norm(u)<tolerance)
+      u=new
+      break
+    end
+    u=new
+  end
+
+  return (u,en(u,V0,Eta0,rho,weights,d)*rho/2)
+end
+
+# Xi(u)=0 is equivalent to the linear equation
+function Xi(u,V,V0,Eta,Eta0,weights,rho,d)
+  order=length(weights[1])
+
+  # init
+  out=zeros(Complex{Float64},order)
+
+  # compute E before running the loop
+  E=en(u,V0,Eta0,rho,weights,d)
+
+  for i in 1:order
+    # k_i
+    k=(1-weights[1][i])/(1+weights[1][i])
+    # S_i
+    S=V[i]-1/(rho*(2*pi)^d)*integrate_legendre_sampled(y->(1-y)/y^3,Eta[i].*u,0,1,weights)
+    # X_i
+    X=k^2/(2*E*rho)
+
+    # U_i
+    out[i]=u[i]-S/(2*E*(X+1))*Phi(S/(E*(X+1)^2))
+  end
+
+  return out
+end
+
+# derivative of Xi (for Newton)
+function dXi(u,V,V0,Eta,Eta0,weights,rho,d)
+  order=length(weights[1])
+
+  # init
+  out=zeros(Complex{Float64},order,order)
+
+  # compute E before the loop
+  E=en(u,V0,Eta0,rho,weights,d)
+
+  for i in 1:order
+    # k_i
+    k=(1-weights[1][i])/(1+weights[1][i])
+    # S_i
+    S=V[i]-1/(rho*(2*pi)^d)*integrate_legendre_sampled(y->(1-y)/y^3,Eta[i].*u,0,1,weights)
+    # X_i
+    X=k^2/(2*E*rho)
+
+    for j in 1:order
+      y=(weights[1][j]+1)/2
+      dS=-1/rho*(1-y)*Eta[i][j]/(2*(2*pi)^d*y^3)*weights[2][j]
+      dE=-1/rho*(1-y)*Eta0[j]/(2*(2*pi)^d*y^3)*weights[2][j]
+      dX=-k^2/(2*E^2*rho)*dE
+      out[i,j]=(i==j ? 1 : 0)-(dS-S*dE/E-S*dX/(X+1))/(2*E*(X+1))*Phi(S/(E*(X+1)^2))-S/(2*E^2*(X+1)^3)*(dS-S*dE/E-2*S*dX/(X+1))*dPhi(S/(E*(X+1)^2))
+    end
+  end
+
+  return out
+end
+
+# energy
+function en(u,V0,Eta0,rho,weights,d)
+  return V0-1/(rho*(2*pi)^d)*integrate_legendre_sampled(y->(1-y)/y^3,Eta0.*u,0,1,weights)
+end
+
diff --git a/figs/numerical.fig/simpleq/simpleq.jl b/figs/numerical.fig/simpleq/simpleq.jl
new file mode 100644
index 0000000..88209b0
--- /dev/null
+++ b/figs/numerical.fig/simpleq/simpleq.jl
@@ -0,0 +1,40 @@
+# \eta
+function eta(x,t,weights,d,v)
+  if d==2
+    return integrate_chebyshev(y->4*((x+t)*y+abs(x-t)*(1-y))*v((x+t)*y+abs(x-t)*(1-y))/sqrt(((x+t)*y+abs(x-t)*(2-y))*((x+t)*(1+y)+abs(x-t)*(1-y))),0,1,length(weights))
+  elseif d==3
+    return (x>t ? 2*t/x : 2)* integrate_legendre(y->2*pi*((x+t)*y+abs(x-t)*(1-y))*v((x+t)*y+abs(x-t)*(1-y)),0,1,weights)
+  end
+end
+
+# initialize V and Eta
+function init_veta(weights,d,v)
+  order=length(weights[1])
+  V=Array{Complex{Float64}}(undef,order)
+  Eta=Array{Array{Complex{Float64}}}(undef,order)
+  Eta0=Array{Complex{Float64}}(undef,order)
+  V0=v(0)
+  for i in 1:order
+    k=(1-weights[1][i])/(1+weights[1][i])
+    V[i]=v(k)
+    Eta[i]=Array{Complex{Float64}}(undef,order)
+    for j in 1:order
+      y=(weights[1][j]+1)/2
+      Eta[i][j]=eta(k,(1-y)/y,weights,d,v)
+    end
+    y=(weights[1][i]+1)/2
+    Eta0[i]=eta(0,(1-y)/y,weights,d,v)
+  end
+  return(V,V0,Eta,Eta0)
+end
+
+# inverse Fourier transform
+function u_x(x,u,weights,d)
+  order=length(weights[1])
+  if d==2
+    out=integrate_legendre_sampled(y->(1-y)/y^3*besselj(0,x*(1-y)/y)/(2*pi),u,0,1,weights)
+  elseif d==3
+    out=integrate_legendre_sampled(y->(1-y)/y^3*sin(x*(1-y)/y)/x/(2*pi^2),u,0,1,weights)
+  end
+  return out
+end
diff --git a/figs/numerical.fig/simpleq/tools.jl b/figs/numerical.fig/simpleq/tools.jl
new file mode 100644
index 0000000..1635a14
--- /dev/null
+++ b/figs/numerical.fig/simpleq/tools.jl
@@ -0,0 +1,20 @@
+# \Phi(x)=2*(1-sqrt(1-x))/x
+function Phi(x)
+  if abs(x)>1e-5
+    return 2*(1-sqrt(1-x))/x
+  else
+    return 1+x/4+x^2/8+5*x^3/64+7*x^4/128+21*x^5/512
+  end
+end
+# \partial\Phi
+function dPhi(x)
+  #if abs(x-1)<1e-5
+  #  @printf(stderr,"warning: dPhi is singular at 1, and evaluating it at (% .8e+i% .8e)\n",real(x),imag(x))
+  #end
+  if abs(x)>1e-5
+    return 1/(sqrt(1-x)*x)-2*(1-sqrt(1-x))/x^2
+  else
+    return 1/4+x/4+15*x^2/64+7*x^3/32+105*x^4/512+99*x^5/512
+  end
+end
+
diff --git a/libs/ian.cls b/libs/ian.cls
new file mode 100644
index 0000000..aec99ad
--- /dev/null
+++ b/libs/ian.cls
@@ -0,0 +1,673 @@
+%%
+%% Ian's class file
+%%
+
+%% TeX format
+\NeedsTeXFormat{LaTeX2e}[1995/12/01]
+
+%% class name
+\ProvidesClass{ian}[2017/09/29]
+
+%% boolean to signal that this class is being used
+\newif\ifianclass
+
+%% options
+% no section numbering in equations
+\DeclareOption{section_in_eq}{\sectionsineqtrue}
+\DeclareOption{section_in_fig}{\sectionsinfigtrue}
+\DeclareOption{section_in_theo}{\PassOptionsToPackage{\CurrentOption}{iantheo}}
+\DeclareOption{section_in_all}{\sectionsineqtrue\sectionsinfigtrue\PassOptionsToPackage{section_in_theo}{iantheo}}
+\DeclareOption{subsection_in_eq}{\subsectionsineqtrue}
+\DeclareOption{subsection_in_fig}{\subsectionsinfigtrue}
+\DeclareOption{subsection_in_theo}{\PassOptionsToPackage{\CurrentOption}{iantheo}}
+\DeclareOption{subsection_in_all}{\subsectionsineqtrue\subsectionsinfigtrue\PassOptionsToPackage{subsection_in_theo}{iantheo}}
+\DeclareOption{no_section_in_eq}{\sectionsineqfalse}
+\DeclareOption{no_section_in_fig}{\sectionsinfigfalse}
+\DeclareOption{no_section_in_theo}{\PassOptionsToPackage{\CurrentOption}{iantheo}}
+\DeclareOption{no_section_in_all}{\sectionsineqfalse\sectionsinfigfalse\PassOptionsToPackage{no_section_in_theo}{iantheo}}
+\DeclareOption{no_subsection_in_eq}{\subsectionsineqfalse}
+\DeclareOption{no_subsection_in_fig}{\subsectionsinfigfalse}
+\DeclareOption{no_subsection_in_theo}{\PassOptionsToPackage{\CurrentOption}{iantheo}}
+\DeclareOption{no_subsection_in_all}{\subsectionsineqfalse\subsectionsinfigfalse\PassOptionsToPackage{no_subsection_in_theo}{iantheo}}
+% reset point
+\DeclareOption{point_reset_at_section}{\PassOptionsToPackage{reset_at_section}{point}}
+\DeclareOption{point_no_reset_at_section}{\PassOptionsToPackage{no_reset_at_section}{point}}
+\DeclareOption{point_reset_at_theo}{\PassOptionsToPackage{reset_at_theo}{point}}
+\DeclareOption{point_no_reset_at_theo}{\PassOptionsToPackage{no_reset_at_theo}{point}}
+
+\def\ian@defaultoptions{
+  \ExecuteOptions{section_in_all, no_subsection_in_all}
+  \ProcessOptions
+
+  %% required packages
+  \RequirePackage{iantheo}
+  \RequirePackage{point}
+  \RequirePackage{color}
+  \RequirePackage{marginnote}
+  \RequirePackage{amssymb}
+  \PassOptionsToPackage{hidelinks}{hyperref}
+  \RequirePackage{hyperref}
+}
+
+%% paper dimensions
+\setlength\paperheight{297mm}
+\setlength\paperwidth{210mm}
+
+%% fonts
+\input{size11.clo}
+\DeclareOldFontCommand{\rm}{\normalfont\rmfamily}{\mathrm}
+\DeclareOldFontCommand{\sf}{\normalfont\sffamily}{\mathsf}
+\DeclareOldFontCommand{\tt}{\normalfont\ttfamily}{\mathtt}
+\DeclareOldFontCommand{\bf}{\normalfont\bfseries}{\mathbf}
+\DeclareOldFontCommand{\it}{\normalfont\itshape}{\mathit}
+
+%% text dimensions
+\hoffset=-50pt
+\voffset=-72pt
+\textwidth=460pt
+\textheight=704pt
+
+
+%% remove default indentation
+\parindent=0pt
+%% indent command
+\def\indent{\hskip20pt}
+
+%% something is wrong with \thepage, redefine it
+\gdef\thepage{\the\c@page}
+
+%% array lines (to use the array environment)
+\setlength\arraycolsep{5\p@}
+\setlength\arrayrulewidth{.4\p@}
+
+
+%% correct vertical alignment at the end of a document
+\AtEndDocument{
+  \vfill
+  \eject
+}
+
+
+%% hyperlinks
+% hyperlinkcounter
+\newcounter{lncount}
+% hyperref anchor
+\def\hrefanchor{%
+\stepcounter{lncount}%
+\hypertarget{ln.\thelncount}{}%
+}
+
+%% define a command and write it to aux file
+\def\outdef#1#2{%
+  % define command%
+  \expandafter\xdef\csname #1\endcsname{#2}%
+  % hyperlink number%
+  \expandafter\xdef\csname #1@hl\endcsname{\thelncount}%
+  % write command to aux%
+  \immediate\write\@auxout{\noexpand\expandafter\noexpand\gdef\noexpand\csname #1\endcsname{\csname #1\endcsname}}%
+  \immediate\write\@auxout{\noexpand\expandafter\noexpand\gdef\noexpand\csname #1@hl\endcsname{\thelncount}}%
+}
+
+%% can call commands even when they are not defined
+\def\safe#1{%
+  \ifdefined#1%
+    #1%
+  \else%
+    {\color{red}\bf?}%
+  \fi%
+}
+
+%% define a label for the latest tag
+%% label defines a command containing the string stored in \tag
+\def\deflabel{
+  \def\label##1{\expandafter\outdef{label@##1}{\safe\tag}}
+
+  \def\ref##1{%
+    % check whether the label is defined (hyperlink runs into errors if this check is omitted)
+    \ifcsname label@##1@hl\endcsname%
+      \hyperlink{ln.\csname label@##1@hl\endcsname}{{\color{blue}\safe\csname label@##1\endcsname}}%
+    \else%
+      \ifcsname label@##1\endcsname%
+	{\color{blue}\csname ##1\endcsname}%
+	\else%
+	{\bf ??}%
+      \fi%
+    \fi%
+  }
+}
+
+
+%% make a custom link at any given location in the document
+\def\makelink#1#2{%
+  \hrefanchor%
+  \outdef{label@#1}{#2}%
+}
+
+
+%% section command
+% counter
+\newcounter{sectioncount}
+% space before section
+\newlength\secskip
+\setlength\secskip{40pt}
+% a prefix to put before the section number, e.g. A for appendices
+\def\sectionprefix{}
+% define some lengths
+\newlength\secnumwidth
+\newlength\sectitlewidth
+\def\section#1{
+  % reset counters
+  \stepcounter{sectioncount}
+  \setcounter{subsectioncount}{0}
+  \ifsectionsineq
+    \setcounter{seqcount}0
+  \fi
+  \ifsectionsinfig
+    \setcounter{figcount}0
+  \fi
+
+  % space before section (if not first)
+  \ifnum\thesectioncount>1
+    \vskip\secskip
+    \penalty-1000
+  \fi
+
+  % hyperref anchor
+  \hrefanchor
+  % define tag (for \label)
+  \xdef\tag{\sectionprefix\thesectioncount}
+
+  % get widths
+  \def\@secnum{{\bf\Large\sectionprefix\thesectioncount.\hskip10pt}}
+  \settowidth\secnumwidth{\@secnum}
+  \setlength\sectitlewidth\textwidth
+  \addtolength\sectitlewidth{-\secnumwidth}
+
+  % print name
+  \parbox{\textwidth}{
+  \@secnum
+  \parbox[t]{\sectitlewidth}{\Large\bf #1}}
+
+  % write to table of contents
+  \iftoc
+    % save lncount in aux variable which is written to toc
+    \immediate\write\tocoutput{\noexpand\expandafter\noexpand\edef\noexpand\csname toc@sec.\thesectioncount\endcsname{\thelncount}}
+    \write\tocoutput{\noexpand\tocsection{#1}{\thepage}}
+  \fi
+
+  %space
+  \par\penalty10000
+  \bigskip\penalty10000
+}
+
+%% subsection
+% counter
+\newcounter{subsectioncount}
+% space before subsection
+\newlength\subsecskip
+\setlength\subsecskip{30pt}
+\def\subsection#1{
+  % counters
+  \stepcounter{subsectioncount}
+  \setcounter{subsubsectioncount}{0}
+  \ifsubsectionsineq
+    \setcounter{seqcount}0
+  \fi
+  \ifsubsectionsinfig
+    \setcounter{figcount}0
+  \fi
+
+  % space before subsection (if not first)
+  \ifnum\thesubsectioncount>1
+    \vskip\subsecskip
+    \penalty-500
+  \fi
+
+  % hyperref anchor
+  \hrefanchor
+  % define tag (for \label)
+  \xdef\tag{\sectionprefix\thesectioncount.\thesubsectioncount}
+
+  % get widths
+  \def\@secnum{{\bf\large\hskip.5cm\sectionprefix\thesectioncount.\thesubsectioncount.\hskip5pt}}
+  \settowidth\secnumwidth{\@secnum}
+  \setlength\sectitlewidth\textwidth
+  \addtolength\sectitlewidth{-\secnumwidth}
+  % print name
+  \parbox{\textwidth}{
+  \@secnum
+  \parbox[t]{\sectitlewidth}{\large\bf #1}}
+
+  % write to table of contents
+  \iftoc
+    % save lncount in aux variable which is written to toc
+    \immediate\write\tocoutput{\noexpand\expandafter\noexpand\edef\noexpand\csname toc@subsec.\thesectioncount.\thesubsectioncount\endcsname{\thelncount}}
+    \write\tocoutput{\noexpand\tocsubsection{#1}{\thepage}}
+  \fi
+
+  % space
+  \par\penalty10000
+  \medskip\penalty10000
+}
+
+%% subsubsection
+% counter
+\newcounter{subsubsectioncount}
+% space before subsubsection
+\newlength\subsubsecskip
+\setlength\subsubsecskip{20pt}
+\def\subsubsection#1{
+  % counters
+  \stepcounter{subsubsectioncount}
+
+  % space before subsubsection (if not first)
+  \ifnum\thesubsubsectioncount>1
+    \vskip\subsubsecskip
+    \penalty-500
+  \fi
+
+  % hyperref anchor
+  \hrefanchor
+  % define tag (for \label)
+  \xdef\tag{\sectionprefix\thesectioncount.\thesubsectioncount.\thesubsubsectioncount}
+
+  % get widths
+  \def\@secnum{{\bf\hskip1.cm\sectionprefix\thesectioncount.\thesubsectioncount.\thesubsubsectioncount.\hskip5pt}}
+  \settowidth\secnumwidth{\@secnum}
+  \setlength\sectitlewidth\textwidth
+  \addtolength\sectitlewidth{-\secnumwidth}
+  % print name
+  \parbox{\textwidth}{
+  \@secnum
+  \parbox[t]{\sectitlewidth}{\large\bf #1}}
+
+  % write to table of contents
+  \iftoc
+    % save lncount in aux variable which is written to toc
+    \immediate\write\tocoutput{\noexpand\expandafter\noexpand\edef\noexpand\csname toc@subsubsec.\thesectioncount.\thesubsectioncount.\thesubsubsectioncount\endcsname{\thelncount}}
+    \write\tocoutput{\noexpand\tocsubsubsection{#1}{\thepage}}
+  \fi
+
+  % space
+  \par\penalty10000
+  \medskip\penalty10000
+}
+
+%% itemize
+\newlength\itemizeskip
+% left margin for items
+\setlength\itemizeskip{20pt}
+\newlength\itemizeseparator
+% space between the item symbol and the text
+\setlength\itemizeseparator{5pt}
+% penalty preceding an itemize
+\newcount\itemizepenalty
+\itemizepenalty=0
+% counter counting the itemize level
+\newcounter{itemizecount}
+
+% item symbol
+\def\itemizept#1{
+  \ifnum#1=1
+    \textbullet
+  \else
+    $\scriptstyle\blacktriangleright$
+  \fi
+}
+
+
+\newlength\current@itemizeskip
+\setlength\current@itemizeskip{0pt}
+\def\itemize{%
+  \par\expandafter\penalty\the\itemizepenalty\medskip\expandafter\penalty\the\itemizepenalty%
+  \addtocounter{itemizecount}{1}%
+  \addtolength\current@itemizeskip{\itemizeskip}%
+  \leftskip\current@itemizeskip%
+}
+\def\enditemize{%
+  \addtocounter{itemizecount}{-1}%
+  \addtolength\current@itemizeskip{-\itemizeskip}%
+  \par\expandafter\penalty\the\itemizepenalty\leftskip\current@itemizeskip%
+  \medskip\expandafter\penalty\the\itemizepenalty%
+}
+
+% item, with optional argument to specify the item point
+% @itemarg is set to true when there is an optional argument
+\newif\if@itemarg
+\def\item{%
+  % check whether there is an optional argument (if there is none, add on empty '[]')
+  \@ifnextchar [{\@itemargtrue\@itemx}{\@itemargfalse\@itemx[]}%
+}
+\newlength\itempt@total
+\def\@itemx[#1]{
+  \if@itemarg
+    \settowidth\itempt@total{#1}
+  \else
+    \settowidth\itempt@total{\itemizept\theitemizecount}
+  \fi
+  \addtolength\itempt@total{\itemizeseparator}
+  \par
+  \medskip
+  \if@itemarg
+    \hskip-\itempt@total#1\hskip\itemizeseparator
+  \else
+    \hskip-\itempt@total\itemizept\theitemizecount\hskip\itemizeseparator
+  \fi
+}
+
+%% prevent page breaks after itemize
+\newcount\previtemizepenalty
+\def\nopagebreakafteritemize{
+  \previtemizepenalty=\itemizepenalty
+  \itemizepenalty=10000
+}
+%% back to previous value
+\def\restorepagebreakafteritemize{
+  \itemizepenalty=\previtemizepenalty
+}
+
+%% enumerate
+\newcounter{enumerate@count}
+\def\enumerate{
+  \setcounter{enumerate@count}0
+  \let\olditem\item
+  \let\olditemizept\itemizept
+  \def\item{
+    % counter
+    \stepcounter{enumerate@count}
+    % set header
+    \def\itemizept{\theenumerate@count.}
+    % hyperref anchor
+    \hrefanchor
+    % define tag (for \label)
+    \xdef\tag{\theenumerate@count}
+    \olditem
+  }
+  \itemize
+}
+\def\endenumerate{
+  \enditemize
+  \let\item\olditem
+  \let\itemizept\olditemizept
+}
+
+
+%% equation numbering
+% counter
+\newcounter{seqcount}
+% booleans (write section or subsection in equation number)
+\newif\ifsectionsineq
+\newif\ifsubsectionsineq
+\def\seqcount{
+  \stepcounter{seqcount}
+  % the output
+  \edef\seqformat{\theseqcount}
+  % add subsection number
+  \ifsubsectionsineq
+    \let\tmp\seqformat
+    \edef\seqformat{\thesubsectioncount.\tmp}
+  \fi
+  % add section number
+  \ifsectionsineq
+    \let\tmp\seqformat
+    \edef\seqformat{\sectionprefix\thesectioncount.\tmp}
+  \fi
+  % define tag (for \label)
+  \xdef\tag{\seqformat}
+  % write number
+  \marginnote{\hfill(\seqformat)}
+}
+%% equation environment compatibility
+\def\equation{\hrefanchor$$\seqcount}
+\def\endequation{$$\@ignoretrue}
+
+
+%% figures
+% counter
+\newcounter{figcount}
+% booleans (write section or subsection in equation number)
+\newif\ifsectionsinfig
+\newif\ifsubsectionsinfig
+% width of figures
+\newlength\figwidth
+\setlength\figwidth\textwidth
+\addtolength\figwidth{-2.5cm}
+% caption
+\def\defcaption{
+  \long\def\caption##1{
+    \stepcounter{figcount}
+
+    % hyperref anchor
+    \hrefanchor
+
+    % the number of the figure
+    \edef\figformat{\thefigcount}
+    % add subsection number
+    \ifsubsectionsinfig
+      \let\tmp\figformat
+      \edef\figformat{\thesubsectioncount.\tmp}
+    \fi
+    % add section number
+    \ifsectionsinfig
+      \let\tmp\figformat
+      \edef\figformat{\sectionprefix\thesectioncount.\tmp}
+    \fi
+
+    % define tag (for \label)
+    \xdef\tag{\figformat}
+
+    % write
+    \hfil fig \figformat: \parbox[t]{\figwidth}{\leavevmode\small##1}
+
+    % space
+    \par\bigskip
+  }
+}
+%% short caption: centered
+\def\captionshort#1{
+  \stepcounter{figcount}
+
+  % hyperref anchor
+  \hrefanchor
+
+  % the number of the figure
+  \edef\figformat{\thefigcount}
+  % add section number
+  \ifsectionsinfig
+  \let\tmp\figformat
+  \edef\figformat{\sectionprefix\thesectioncount.\tmp}
+  \fi
+
+  % define tag (for \label)
+  \xdef\tag{\figformat}
+
+  % write
+  \hfil fig \figformat: {\small#1}
+
+  %space
+  \par\bigskip
+}
+
+%% environment
+\def\figure{
+  \par
+  \vfil\penalty100\vfilneg
+  \bigskip
+}
+\def\endfigure{
+  \par
+  \bigskip
+}
+
+
+%% start appendices
+\def\appendix{
+  \vfill
+  \pagebreak
+
+  % counter
+  \setcounter{sectioncount}0
+
+  % prefix
+  \def\sectionprefix{A}
+
+  % write
+  {\bf \LARGE Appendices}\par\penalty10000\bigskip\penalty10000
+
+  % add a mention in the table of contents
+  \iftoc
+    \immediate\write\tocoutput{\noexpand\tocappendices}\penalty10000
+  \fi
+
+  %% uncomment for new page for each appendix
+  %\def\seqskip{\vfill\pagebreak}
+}
+
+
+%% bibliography
+% size of header
+\newlength\bibheader
+\def\thebibliography#1{
+  \hrefanchor
+
+  % add a mention in the table of contents
+  \iftoc
+    % save lncount in aux variable which is written to toc
+    \immediate\write\tocoutput{\noexpand\expandafter\noexpand\edef\noexpand\csname toc@references\endcsname{\thelncount}}
+    \write\tocoutput{\noexpand\tocreferences{\thepage}}\penalty10000
+  \fi
+
+  % write
+  {\bf \LARGE References}\par\penalty10000\bigskip\penalty10000
+  % width of header
+  \settowidth\bibheader{[#1]}
+  \leftskip\bibheader
+}
+% end environment
+\def\endthebibliography{
+  \par\leftskip0pt
+} 
+
+%% bibitem command
+\def\bibitem[#1]#2{%
+  \hrefanchor%
+  \outdef{label@cite#2}{#1}%
+  \hskip-\bibheader%
+  \makebox[\bibheader]{\cite{#2}\hfill}%
+}
+
+%% cite command (adapted from latex.ltx)
+% @tempswa is set to true when there is an optional argument
+\newif\@tempswa
+\def\cite{%
+  % check whether there is an optional argument (if there is none, add on empty '[]')
+  \@ifnextchar [{\@tempswatrue\@citex}{\@tempswafalse\@citex[]}%
+}
+% command with optional argument
+\def\@citex[#1]#2{\leavevmode%
+  % initialize loop
+  \let\@citea\@empty%
+  % format
+  \@cite{%
+    % loop over ',' separated list
+    \@for\@citeb:=#2\do{%
+      % text to add at each iteration of the loop (separator between citations)
+      \@citea\def\@citea{,\ }%
+      % add entry to citelist
+      \@writecitation{\@citeb}%
+      \ref{cite\@citeb}%
+    }%
+  }%
+  % add optional argument text (as an argument to '\@cite')
+  {#1}%
+}
+\def\@cite#1#2{%
+  [#1\if@tempswa , #2\fi]%
+}
+%% add entry to citelist after checking it has not already been added
+\def\@writecitation#1{%
+  \ifcsname if#1cited\endcsname%
+  \else%
+    \expandafter\newif\csname if#1cited\endcsname%
+    \immediate\write\@auxout{\string\citation{#1}}%
+  \fi%
+}
+
+%% table of contents
+% boolean
+\newif\iftoc
+\def\tableofcontents{
+  {\bf \large Table of contents:}\par\penalty10000\bigskip\penalty10000
+
+  % copy content from file
+  \IfFileExists{\jobname.toc}{\input{\jobname.toc}}{{\tt error: table of contents missing}}
+
+  % open new toc
+  \newwrite\tocoutput
+  \immediate\openout\tocoutput=\jobname.toc
+
+  \toctrue
+}
+%% close file
+\AtEndDocument{
+  % close toc
+  \iftoc
+    \immediate\closeout\tocoutput
+  \fi
+}
+
+
+%% fill line with dots
+\def\leaderfill{\leaders\hbox to 1em {\hss. \hss}\hfill}
+
+%% same as sectionprefix
+\def\tocsectionprefix{}
+
+%% toc formats
+\newcounter{tocsectioncount}
+\def\tocsection #1#2{
+  \stepcounter{tocsectioncount}
+  \setcounter{tocsubsectioncount}{0}
+  \setcounter{tocsubsubsectioncount}{0}
+  % write
+  \smallskip\hyperlink{ln.\csname toc@sec.\thetocsectioncount\endcsname}{{\bf \tocsectionprefix\thetocsectioncount}.\hskip5pt {\color{blue}#1}\leaderfill#2}\par
+}
+\newcounter{tocsubsectioncount}
+\def\tocsubsection #1#2{
+  \stepcounter{tocsubsectioncount}
+  \setcounter{tocsubsubsectioncount}{0}
+  % write
+  {\hskip10pt\hyperlink{ln.\csname toc@subsec.\thetocsectioncount.\thetocsubsectioncount\endcsname}{{\bf \thetocsubsectioncount}.\hskip5pt {\color{blue}\small #1}\leaderfill#2}}\par
+}
+\newcounter{tocsubsubsectioncount}
+\def\tocsubsubsection #1#2{
+  \stepcounter{tocsubsubsectioncount}
+  % write
+  {\hskip20pt\hyperlink{ln.\csname toc@subsubsec.\thetocsectioncount.\thetocsubsectioncount.\thetocsubsubsectioncount\endcsname}{{\bf \thetocsubsubsectioncount}.\hskip5pt {\color{blue}\small #1}\leaderfill#2}}\par
+}
+\def\tocappendices{
+  \medskip
+  \setcounter{tocsectioncount}0
+  {\bf Appendices}\par
+  \smallskip
+  \def\tocsectionprefix{A}
+}
+\def\tocreferences#1{
+  \medskip
+  {\hyperlink{ln.\csname toc@references\endcsname}{{\color{blue}\bf References}\leaderfill#1}}\par
+  \smallskip
+}
+
+
+%% definitions that must be loaded at begin document
+\let\ian@olddocument\document
+\def\document{
+  \ian@olddocument
+
+  \deflabel
+  \defcaption
+}
+
+%% end
+\ian@defaultoptions
+\endinput
diff --git a/libs/iantheo.sty b/libs/iantheo.sty
new file mode 100644
index 0000000..d33a93d
--- /dev/null
+++ b/libs/iantheo.sty
@@ -0,0 +1,162 @@
+%%
+%% iantheorem package:
+%%   Ian's customized theorem command
+%%
+
+%% boolean to signal that this package was loaded
+\newif\ifiantheo
+
+%% TeX format
+\NeedsTeXFormat{LaTeX2e}[1995/12/01]
+
+%% package name
+\ProvidesPackage{iantheo}[2016/11/10]
+
+%% options
+\newif\ifsectionintheo
+\DeclareOption{section_in_theo}{\sectionintheotrue}
+\DeclareOption{no_section_in_theo}{\sectionintheofalse}
+\newif\ifsubsectionintheo
+\DeclareOption{subsection_in_theo}{\subsectionintheotrue}
+\DeclareOption{no_subsection_in_theo}{\subsectionintheofalse}
+
+\def\iantheo@defaultoptions{
+  \ExecuteOptions{section_in_theo, no_subsection_in_theo}
+  \ProcessOptions
+
+  %%% reset at every new section
+  \ifsectionintheo
+    \let\iantheo@oldsection\section
+    \gdef\section{\setcounter{theocount}{0}\iantheo@oldsection}
+  \fi
+
+  %% reset at every new subsection
+  \ifsubsectionintheo
+    \let\iantheo@oldsubsection\subsection
+    \gdef\subsection{\setcounter{theocount}{0}\iantheo@oldsubsection}
+  \fi
+}
+
+
+%% delimiters
+\def\delimtitle#1{
+  \par%
+  \leavevmode%
+  \raise.3em\hbox to\hsize{%
+    \lower0.3em\hbox{\vrule height0.3em}%
+    \hrulefill%
+    \ \lower.3em\hbox{#1}\ %
+    \hrulefill%
+    \lower0.3em\hbox{\vrule height0.3em}%
+  }%
+  \par\penalty10000%
+}
+
+%% callable by ref
+\def\delimtitleref#1{
+  \par%
+%
+  \ifdefined\ianclass%
+    % hyperref anchor%
+    \hrefanchor%
+    % define tag (for \label)%
+    \xdef\tag{#1}%
+  \fi%
+%
+  \leavevmode%
+  \raise.3em\hbox to\hsize{%
+    \lower0.3em\hbox{\vrule height0.3em}%
+    \hrulefill%
+    \ \lower.3em\hbox{\bf #1}\ %
+    \hrulefill%
+    \lower0.3em\hbox{\vrule height0.3em}%
+  }%
+  \par\penalty10000%
+}
+
+%% no title
+\def\delim{
+  \par%
+  \leavevmode\raise.3em\hbox to\hsize{%
+    \lower0.3em\hbox{\vrule height0.3em}%
+    \hrulefill%
+    \lower0.3em\hbox{\vrule height0.3em}%
+  }%
+  \par\penalty10000%
+}
+
+%% end delim
+\def\enddelim{
+  \par\penalty10000%
+  \leavevmode%
+  \raise.3em\hbox to\hsize{%
+    \vrule height0.3em\hrulefill\vrule height0.3em%
+  }%
+  \par%
+}
+
+
+%% theorem
+% counter
+\newcounter{theocount}
+% booleans (write section or subsection in equation number)
+\def\theo#1{
+  \stepcounter{theocount}
+  \ifdefined\ianclass
+    % hyperref anchor
+    \hrefanchor
+  \fi
+  % the number
+  \def\formattheo{\thetheocount}
+  % add subsection number
+  \ifsubsectionintheo
+    \let\tmp\formattheo
+    \edef\formattheo{\thesubsectioncount.\tmp}
+  \fi
+  % add section number
+  \ifsectionintheo
+    \let\tmp\formattheo
+    \edef\formattheo{\sectionprefix\thesectioncount.\tmp}
+  \fi
+  % define tag (for \label)
+  \xdef\tag{\formattheo}
+  % write
+  \delimtitle{\bf #1 \formattheo}
+}
+\let\endtheo\enddelim
+%% theorem headers with name
+\def\theoname#1#2{
+  \theo{#1}\hfil({\it #2})\par\penalty10000\medskip%
+}
+
+
+%% qed symbol
+\def\qedsymbol{$\square$}
+\def\qed{\penalty10000\hfill\penalty10000\qedsymbol}
+
+
+%% compatibility with article class
+\ifdefined\ianclasstrue
+  \relax
+\else
+  \def\thesectioncount{\thesection}
+  \def\thesubsectioncount{\thesubsection}
+  \def\sectionprefix{}
+\fi
+
+
+%% prevent page breaks after displayed equations
+\newcount\prevpostdisplaypenalty
+\def\nopagebreakaftereq{
+  \prevpostdisplaypenalty=\postdisplaypenalty
+  \postdisplaypenalty=10000
+}
+%% back to previous value
+\def\restorepagebreakaftereq{
+  \postdisplaypenalty=\prevpostdisplaypenalty
+}
+
+
+%% end
+\iantheo@defaultoptions
+\endinput
diff --git a/libs/largearray.sty b/libs/largearray.sty
new file mode 100644
index 0000000..ad5753b
--- /dev/null
+++ b/libs/largearray.sty
@@ -0,0 +1,19 @@
+%%
+%% largearray package:
+%%   Array spanning the entire line
+%%
+
+%% TeX format
+\NeedsTeXFormat{LaTeX2e}[1995/12/01]
+
+%% package name
+\ProvidesPackage{largearray}[2016/11/10]
+
+\RequirePackage{array}
+
+%% array spanning the entire line
+\newlength\largearray@width
+\setlength\largearray@width\textwidth
+\addtolength\largearray@width{-10pt}
+\def\largearray{\begin{array}{@{}>{\displaystyle}l@{}}\hphantom{\hspace{\largearray@width}}\\[-.5cm]}
+\def\endlargearray{\end{array}}
diff --git a/libs/point.sty b/libs/point.sty
new file mode 100644
index 0000000..da06361
--- /dev/null
+++ b/libs/point.sty
@@ -0,0 +1,107 @@
+%%
+%% Points package:
+%%   \point commands
+%%
+
+%% TeX format
+\NeedsTeXFormat{LaTeX2e}[1995/12/01]
+
+%% package name
+\ProvidesPackage{point}[2017/06/13]
+
+%% options
+\newif\ifresetatsection
+\DeclareOption{reset_at_section}{\resetatsectiontrue}
+\DeclareOption{no_reset_at_section}{\resetatsectionfalse}
+\newif\ifresetatsubsection
+\DeclareOption{reset_at_subsection}{\resetatsubsectiontrue}
+\DeclareOption{no_reset_at_subsection}{\resetatsubsectionfalse}
+\newif\ifresetattheo
+\DeclareOption{reset_at_theo}{\resetattheotrue}
+\DeclareOption{no_reset_at_theo}{\resetattheofalse}
+
+\def\point@defaultoptions{
+  \ExecuteOptions{reset_at_section, reset_at_subsection, no_reset_at_theo}
+  \ProcessOptions
+
+  %% reset at every new section
+  \ifresetatsection
+    \let\point@oldsection\section
+    \gdef\section{\resetpointcounter\point@oldsection}
+  \fi
+  %% reset at every new subsection
+  \ifresetatsubsection
+    \let\point@oldsubsection\subsection
+    \gdef\subsection{\resetpointcounter\point@oldsubsection}
+  \fi
+
+  %% reset at every new theorem
+  \ifresetattheo
+    \ifdefined\iantheotrue
+      \let\point@oldtheo\theo
+      \gdef\theo{\resetpointcounter\point@oldtheo}
+    \fi
+  \fi
+}
+
+
+%% point
+% counter
+\newcounter{pointcount}
+\def\point{
+  \stepcounter{pointcount}
+  \setcounter{subpointcount}{0}
+  % hyperref anchor (only if the class is 'ian')
+  \ifdefined\ifianclass
+    \hrefanchor
+    % define tag (for \label)
+    \xdef\tag{\arabic{pointcount}}
+  \fi
+  % header
+  \indent{\bf \arabic{pointcount}\ - }
+}
+
+%% subpoint
+% counter
+\newcounter{subpointcount}
+\def\subpoint{
+  \stepcounter{subpointcount}
+  \setcounter{subsubpointcount}0
+  % hyperref anchor (only if the class is 'ian')
+  \ifdefined\ifianclass
+    \hrefanchor
+    % define tag (for \label)
+    \xdef\tag{\arabic{pointcount}-\arabic{subpointcount}}
+  \fi
+  % header
+  \indent\hskip.5cm{\bf \arabic{pointcount}-\arabic{subpointcount}\ - }
+}
+
+%% subsubpoint
+% counter
+\newcounter{subsubpointcount}
+\def\subsubpoint{
+  \stepcounter{subsubpointcount}
+  % hyperref anchor (only if the class is 'ian')
+  \ifdefined\ifianclass
+    \hrefanchor
+    % define tag (for \label)
+    \xdef\tag{\arabic{pointcount}-\arabic{subpointcount}-\arabic{subsubpointcount}}
+  \fi
+  \indent\hskip1cm{\bf
+\arabic{pointcount}-\arabic{subpointcount}-\arabic{subsubpointcount}\ - }
+}
+
+
+%% reset point counters
+\def\resetpointcounter{
+  \setcounter{pointcount}{0}
+  \setcounter{subpointcount}{0}
+  \setcounter{subsubpointcount}{0}
+}
+
+
+
+%% end
+\point@defaultoptions
+\endinput