%% Copyright 2021 Ian Jauslin
%% 
%% Licensed under the Apache License, Version 2.0 (the "License");
%% you may not use this file except in compliance with the License.
%% You may obtain a copy of the License at
%% 
%%     http://www.apache.org/licenses/LICENSE-2.0
%% 
%% Unless required by applicable law or agreed to in writing, software
%% distributed under the License is distributed on an "AS IS" BASIS,
%% WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
%% See the License for the specific language governing permissions and
%% limitations under the License.

\documentclass[subsection_in_all]{ian}

\usepackage{code}
\usepackage{dsfont}
\usepackage{largearray}
\usepackage{dlmf}

\def\Eta{H}

\begin{document}

\hbox{}
\hfil{\bf\Large
{\tt simplesolv}\par
\medskip
\hfil \large v0.3
}

\vfill

\tableofcontents

\vfill
\eject

\setcounter{page}1
\pagestyle{plain}

\indent
{\tt simplesolv} is a tool to compute the solution of the equations of the ``Simplified approach'' to the repulsive Bose gas introduced in\-~\cite{CJL20,CJL20b,CHe21}.
This approach provides an approximation to various observables of the ground state of the Hamiltonian
\begin{equation}
  H_N=-\frac12\sum_{i=1}^N\Delta_i+\sum_{1\leqslant i\leqslant j\leqslant N}v(|x_i-x_j|)
\end{equation}
in three dimensions, with periodic boundary conditions, in the thermodynamic limit $N\to\infty$ at fixed density $\rho$.
\bigskip

\indent
{\tt simplesolv} is written in {\tt julia}.
The source code is located in the {\tt src} directory of this bundle.
Throughout the documentation, we will refer to the directory containing the {\tt src} directory as the ``installation directory'', and will denote it by the bash variable {\tt\$SIMPLESOLV} (so that the main julia file, for instance, is located at {\tt \$SIMPLESOLV/src/main.jl}).
\vskip20pt

\section{Basic usage}
\indent
Denoting the location of the installation directory by {\tt\$SIMPLESOLV}, {\tt simplesolv} is run by calling
\begin{code}
  julia \$SIMPLESOLV/src/main.jl [args] <command>
\end{code}
where the optional arguments {\tt [args]} take the form {\tt [-p params] [-U potential] [-M method] [-s savefile]}.
\bigskip

\indent
A few commands support multithreaded execution.
To enable {\tt julia} to run on several processors, it should be run with the {\tt -p} option.
For example, to run on 8 CPUs, run
\begin{code}
  julia -p 8 \$SIMPLESOLV/src/main.jl [args] <command>
\end{code}
\bigskip

\indent
{\tt command} specifies which computation is to be carried out, such as {\tt energy} to compute the ground state energy, or {\tt condensate\_fraction} for the uncondensed fraction.
The list of available commands depends on the \refname{sec:methods}{{\tt method}} argument, which specifies one of the available methods to solve the equation at hand.
The available methods are (see section\-~\ref{sec:methods} for further details)
\begin{itemize}
  \item \refname{sec:easyeq}{{\tt easyeq}} ({\bf default}) for the Simple or Medium equation, or any interpolation between them, with a soft potential using the Newton algorithm,
  \item \refname{sec:anyeq}{{\tt anyeq}} for any equation in the ``Simplified approach'' using the Newton algorithm.
  \item \refname{sec:simpleq-Kv}{{\tt simpleq-Kv}} for the Simple equation using explicit expressions involving $\mathfrak Kv$ (see\-~(\ref{Kv})),
  \item \refname{sec:simpleq-hardcore}{{\tt simpleq-hardcore}} for the Simple equation with a hard core potential using the Newton algorithm,
  \item \refname{sec:simpleq-iteration}{{\tt simpleq-iteration}} for the Simple equation with a soft potential using the iteration defined in\-~\cite{CJL20}.
\end{itemize}
Each method is described in detail below, along with the list of commands ({\tt command}) and parameters ({\tt params}) compatible with them.
\bigskip

\indent
{\tt params} should be a `{\tt;}' separated list of entries, each of which is of the form {\tt key=value}.
For example {\tt -p "rho=1e-6;v\_a=2"}.
(Note that you should not end the list of parameters by a `{\tt;}', otherwise {\tt simplesolv} will interpret that as there being an empty parameter entry, which it cannot split into a key and value, and will fail.)
\bigskip

\indent
\refname{sec:potentials}{{\tt potential}} specifies which potential $v$ should be used, from the following list (see section\-~\ref{sec:potentials} for further details).
\begin{itemize}
  \item \refname{sec:exp}{{\tt exp}} ({\bf default}) for $v(|x|)=ae^{-|x|}$,
  \item \refname{sec:tent}{{\tt tent}} for $v(|x|)=\mathds 1_{|x|<b}a\frac{2\pi}{3}(1-\frac{|x|}b)^2(\frac{|x|}b+2)$,
  \item \refname{sec:expcry}{{\tt expcry}} for $v(|x|)=e^{-|x|}-ae^{-b|x|}$,
  \item \refname{sec:npt}{{\tt npt}} for $v(|x|)=x^2e^{-|x|}$,
  \item \refname{sec:alg}{{\tt alg}} for $v(|x|)=\frac1{1+\frac14|x|^4}$,
  \item \refname{sec:algwell}{{\tt algwell}} for $v(|x|)=\frac{1+a|x|^4}{(1+|x|^2)^4}$.
  \item \refname{sec:exact}{{\tt exact}} for $v(|x|)=\frac{12c( |x|^6b^6(2e-b^2) +b^4|x|^4(9e-7b^2) +4b^2|x|^2(3e-2b^2) +(5e+16b^2))}{(1+b^2|x|^2)^2(4+b^2|x|^2)^2((1+b^2|x|^2)^2-c)}$
\end{itemize}
The parameters in the potential can be set using the {\tt params} argument: to set $a$ set {\tt v\_a}, to set $b$ set {\tt v\_b}, to set $c$ set {\tt v\_c}, and to set $e$ set {\tt v\_e}.
\bigskip

\indent
{\tt savefile} can be used to accelerate the computation of observables in the \refname{eq:anyeq_compleq}{{\tt compleq}} equation.
Indeed, as is discussed in section\-~\ref{sec:anyeq}, the computation of \refname{eq:anyeq_compleq}{{\tt compleq}} is based on the computation of a large matrix, which can be pre-computed, saved in a file using the \refname{command:anyeq_save_Abar}{{\tt save\_Abar}} command, and reused by specifying that file in the {\tt savefile} argument.

\section{Methods}\label{sec:methods}
\indent
In this section, we describe the different computation methods.
\bigskip

\subsection{\tt easyeq}\label{sec:easyeq}
\indent
This method is used to solve a family of equations, called {\tt easyeq}, that interpolate between the Simple equation and the Medium equation:
\begin{equation}
  -\Delta u
  =
  v(1-u)-2\rho K+\rho^2 L
  \label{easyeq}
\end{equation}
with
\begin{equation}
  K:=\beta_K u\ast S+(1-\beta_K)\frac{2e}\rho u
  ,\quad
  L:=\beta_Lu\ast u\ast S+(1-\beta_L)\frac{2e}\rho u\ast u
  \label{easyeq_KL}
\end{equation}
\begin{equation}
  S:=(1-u)v
  ,\quad
  e:=\frac\rho2\int dx\ (1-u(|x|))v(|x|)
  .
  \label{easyeq_Se}
\end{equation}
for a soft potential $v$ at density $\rho>0$.
\bigskip

\indent
\makelink{eq:easyeq_simpleq}{}
\makelink{eq:easyeq_medeq}{}
The special choice $\beta_K=\beta_L=0$ is called the Simple equation ({\tt simpleq}), and the choice $\beta_K=\beta_L=1$ is called the Medium equation ({\tt medeq})
\bigskip

\subsubsection{Usage}
\indent
Unless otherwise noted, this method takes the following parameters (specified via the {\tt [-p params]} flag, as a `{\tt;}' separated list of entries).
\begin{itemize}
  \item\makelink{param:easyeq_rho}{}
  {\tt rho} ({\tt Float64}, default: $10^{-6}$): density $\rho$.

  \item\makelink{param:easyeq_tolerance}{}
  {\tt tolerance} ({\tt Float64}, default: $10^{-11}$): maximal size of final step in Newton iteration.

  \item\makelink{param:easyeq_maxiter}{}
  {\tt maxiter} ({\tt Int64}, default: 21): maximal number of iterations of the Newton algorithm before giving up.

  \item\makelink{param:easyeq_order}{}
  {\tt order} ({\tt Int64}, default: 100): order used for all Gauss quadratures (denoted by $N$ below).

  \item\makelink{param:easyeq_minlrho_init}{}
  {\tt minlrho\_init} ({\tt Float64}, default: {\tt rho}): to initialize the Newton algorithm, we first compute the solution for a smaller $\rho$, {\tt minlrho} is the minimal value for $\log_{10}\rho$ to start this initialization process.

  \item\makelink{param:easyeq_nlrho_init}{}
  {\tt nlrho\_init} ({\tt Int64}, default: 1): number of steps in the initialization process described above. Set to 1 to disable the incremental initialization process.

  \item\makelink{param:easyeq_bK}{}
  {\tt bK}, {\tt bL} ({\tt Float64}, default: 1, 1): the values of $\beta_K$ and $\beta_L$.

  \item\makelink{param:easyeq_eq}{}
  {\tt eq} ({\tt String}, default: ``{\tt simpleq}'', acceptable values: ``{\tt simpleq}'', ``{\tt medeq}''): A shortcut to select either the \refname{eq:easyeq_simpleq}{Simple equation} ($\beta_K=\beta_L=0$) or the \refname{eq:easyeq_medeq}{Medium equation} ($\beta_L=\beta_K=1$). When this option is set, neither {\tt bK} nor {\tt bL} should be set.
\end{itemize}
\bigskip

The available {\tt commands} are the following.

\begin{itemize}
  \item {\tt energy}\makelink{command:easyeq_energy}{}:
  compute the energy $e$ at a given $\rho$.\par
  \underline{Output}: [$e$] [Newton error $\epsilon$].

  \item {\tt energy\_rho}\makelink{command:easyeq_energy_rho}{}:
  compute the energy $e$ as a function of $\rho$.
  The Newton algorithm is initialized with the hardcore scattering solution (\ref{easyeq_init}) for the lowest $\rho$, and with the previously computed $\rho$ for the larger densities.\par
  \underline{Disabled parameters}: \refname{param:easyeq_rho}{{\tt rho}}, \refname{param:easyeq_minlrho_init}{{\tt minlrho\_init}} and \refname{param:easyeq_nlrho_init}{{\tt nlrho\_init}}.\par
  \underline{Extra parameters}:
  \begin{itemize}
    \item\makelink{param:easyeq_minlrho}{}
    {\tt minlrho} ({\tt Float64}, default: $10^{-6}$): minimal value for $\log_{10}\rho$.

    \item\makelink{param:easyeq_maxlrho}{}
    {\tt maxlrho} ({\tt Float64}, default: $10^{2}$): maximal value for $\log_{10}\rho$.

    \item\makelink{param:easyeq_nlrho}{}
    {\tt nlrho} ({\tt Int64}, default: $100$): number of values for $\rho$ (spaced logarithmically).

    \item\makelink{param:easyeq_rhos}{}
    {\tt rhos} ({\tt Array\{Float64\}}, default: $(10^{{\tt minlrho}+\frac{{\tt maxlrho}-{\tt minlrho}}{{\tt nlrho}}n})_n$: list of values for $\rho$, specified as a `{\tt,}' separated list.
  \end{itemize}
  \underline{Output} (one line for each value of $\rho$): [$\rho$] [$e$] [Newton error $\epsilon$].

  \item\makelink{command:easyeq_condensate_fraction}{}
  {\tt condensate\_fraction}: compute the uncondensed fraction $\eta$ at a given $\rho$.\par
  \underline{Output}: [$\eta$] [Newton error $\epsilon$].

  \item\makelink{command:easyeq_condensate_fraction_rho}{}
  {\tt condensate\_fraction\_rho}: compute the uncondensed fraction $\eta$ as a function of $\rho$.
  The Newton algorithm is initialized with the hardcore scattering solution (\ref{easyeq_init}) for the lowest $\rho$, and with the previously computed $\rho$ for the larger densities.\par
  \underline{Disabled parameters}: same as \refname{command:easyeq_energy_rho}{{\tt energy\_rho}}.\par
  \underline{Extra parameters}: same as \refname{command:easyeq_energy_rho}{{\tt energy\_rho}}.\par
  \underline{Output} (one line for each value of $\rho$): [$\rho$] [$\eta$] [Newton error $\epsilon$].

  \item\makelink{command:easyeq_uk}{}
  {\tt uk}: compute the Fourier transform $\hat u(|k|)$.
  The values $|k|$ at which $\hat u$ is computed are those coming from the Gauss quadratures, and cannot be set.\par
  \underline{Output} (one line for each value of $|k|$): [$|k|$] [$\hat u(|k|)$]

  \item\makelink{command:easyeq_ux}{}
  {\tt ux}: compute $u$ as a function of $|x|$.\par
  \underline{Extra parameters}:
  \begin{itemize}
    \item\makelink{param:easyeq_xmin}{}
    {\tt xmin} ({\tt Float64}, default: 0): minimum of the range of $|x|$ to be printed.

    \item\makelink{param:easyeq_xmax}{}
    {\tt xmax} ({\tt Float64}, default: 100): maximum of the range of $|x|$ to be printed.

    \item\makelink{param:easyeq_nx}{}
    {\tt nx} ({\tt Int64}, default: 100): number of points to print (linearly spaced).
  \end{itemize}
  \underline{Output} (one line for each value of $x$): [$|x|$] [$u(|x|)$]

  \item\makelink{command:easyeq_uux}{}
  {\tt uux}: compute $2u-\rho u\ast u$ as a function of $|x|$.\par
  \underline{Extra parameters}: Same as \refname{command:easyeq_ux}{{\tt ux}}.\par
  \underline{Output} (one line for each value of $x$): [$|x|$] [$2u(|x|)-\rho u\ast u(|x|)$]

\end{itemize}

\subsubsection{Description}
\point{\bf Fourier space formulation.}
The computation is carried out in Fourier space.
We take the convention that the Fourier transform of a function $f(|x|)$ is
\begin{equation}
  \hat f(|k|)=\int_{\mathbb R^3} dx\ e^{ikx}f(|x|)
  =\frac{4\pi}{|k|}\int_0^\infty dr\ \frac{\sin(|k|r)}rf(r)
  .
\end{equation}
In Fourier space, (\ref{easyeq}) becomes
\begin{equation}
  k^2\hat u=\hat S-2\rho\hat A_K\hat u+\rho^2\hat A_L\hat u^2
\end{equation}
\begin{equation}
  \hat A_K:=\beta_K\hat S+(1-\beta_K)\frac{2e}\rho
  ,\quad
  \hat A_L:=\beta_L\hat S+(1-\beta_L)\frac{2e}\rho
\end{equation}
with
\begin{equation}
  \hat S(|k|)=\int_{\mathbb R^3} dx\ e^{ikx}(1-u(|x|))v(|x|)
  =\hat v(k)-\hat u\hat\ast\hat v(k)
  ,\quad
  \hat f\hat\ast\hat g(|k|):=\int_{\mathbb R^3}\frac{dp}{(2\pi)^3}\ \hat f(|k-p|)\hat g(|p|)
  .
\end{equation}
We write this as a quadratic equation for $\hat u$, and solve it, keeping the solution that decays as $|k|\to\infty$:
\begin{equation}
  \rho\hat u(|k|)=
  \frac{A_K(|k|)}{A_L(|k|)}\left(\xi(|k|)+1-\sqrt{(\xi(|k|)+1)^2-\frac{A_L(|k|)}{A_K^2(|k|)}\hat S(|k|)}\right)
\end{equation}
that is,
\begin{equation}
  \rho\hat u(|k|)=
  \frac{\hat S(|k|)}{2A_K(|k|)(\xi(|k|)+1)}\Phi\left({\textstyle\frac{A_L(|k|)}{A_K^2(|k|)}\frac{\hat S(|k|)}{(\xi(|k|)+1)^2}}\right)
  \label{uk_easyeq}
\end{equation}
with
\begin{equation}
  \xi(|k|):=\frac{k^2}{2\rho\hat A_K}
  ,\quad
  \Phi(x):=\frac{2(1-\sqrt{1-x})}x
  .
\end{equation}
Furthermore, using bipolar coordinates (see lemma\-~\ref{lemma:bipolar}), we write $\hat S$ as
\begin{equation}
  \hat S(|k|)=\hat v(|k|)-\frac1{8\pi^3}\int_0^\infty dt\ t\hat u(t)\Eta(|k|,t)
  ,\quad
  \Eta(y,t):=\frac{2\pi}y\int_{|y-t|}^{y+t}ds\ s\hat v(s)
  .
  \label{easyeq_Eta}
\end{equation}
By a simple change of variables,
\begin{equation}
  \Eta(y,t)=
  4\pi\left(\mathds 1_{y>t}\frac{t}y+\mathds 1_{y\leqslant t}\right)
  \int_0^1 ds\ ((y+t)s+|y-t|(1-s))v((y+t)s+|y-t|(1-s))
  .
  \label{Eta}
\end{equation}
\bigskip

\point{\bf Evaluating integrals.}
To compute these integrals numerically, we will use Gauss-Legendre quadratures:
\begin{equation}
  \int_0^1 dt\ f(t)\approx\frac12\sum_{i=1}^N w_if\left({\textstyle\frac{r_i+1}2}\right)
\end{equation}
where $w_i$ and $r_i$ are the Gauss-Legendre {\it weights} and {\it abcissa}.
The order $N$ corresponds to the parameter \refname{param:easyeq_order}{{\tt order}}.
The error made by the quadrature is estimated in appendix\-~\ref{appGL}.
We compactify the integrals to the interval $(0,1)$: $y:=\frac1{t+1}$:
\begin{equation}
  \hat S(|k|)=\hat v(|k|)-\frac1{8\pi^3}\int_0^1 dy\ \frac{(1-y)\hat u(\frac{1-y}y)\Eta(|k|,\frac{1-y}y)}{y^3}
\end{equation}
so, using the Gauss-Legendre quadrature, we approximate
\begin{equation}
  \hat S(|k|)\approx\hat v(|k|)-\frac1{16\pi^3}\sum_{j=1}^N w_j \frac{(1-y_j)\hat u(\frac{1-y_j}{y_j})\Eta(|k|,\frac{1-y_j}{y_j})}{y_j^3}
  ,\quad
  y_j:=\frac{r_j+1}2
  .
\end{equation}
This suggests a natural discretization of Fourier space: let
\begin{equation}
  k_i:=\frac{1-r_i}{1+r_i}\equiv\frac{1-y_i}{y_i}
  .
  \label{easyeq_ki}
\end{equation}
Thus, defining
\begin{equation}
  \mathbb U_i:=\rho\hat u(k_i)
  ,\quad
  \hat v_i:=\hat v(k_i)
\end{equation}
and approximate\-~(\ref{uk_easyeq}):
\begin{equation}
  \mathbb U_i=\frac{\mathbb T_i}{2(\mathbb X_i+1)}\Phi\left({\textstyle\mathbb B_i\frac{\mathbb T_i}{(\mathbb X_i+1)^2}}\right)
  \label{bbU_easyeq}
\end{equation}
with
\begin{equation}
  \mathbb A_{K,i}:=\beta_K\mathbb S_i+(1-\beta_K)\mathbb E
  ,\quad
  \mathbb B_{i}:=\frac{\beta_L\mathbb S_i+(1-\beta_L)\mathbb E}{\mathbb A_{K,i}}
  ,\quad
  \mathbb T_i:=\frac{\mathbb S_i}{\mathbb A_{K,i}}
\end{equation}
\begin{equation}
  \mathbb S_i:=\hat v_i-\frac1{16\pi^3\rho}\sum_{j=1}^N w_j\frac{(1-y_j)\mathbb U_j\Eta(k_i,k_j)}{y_j^3}
  ,\quad
  \mathbb E:=\hat v(0)-\frac1{16\pi^3\rho}\sum_{j=1}^N w_j\frac{(1-y_j)\mathbb U_j\Eta(0,k_j)}{y_j^3}
\end{equation}
\begin{equation}
  \mathbb X_i:=\frac{k_i^2}{2\rho\mathbb A_{K,i}}
  .
\end{equation}
This is a discrete equation for the vector $(\mathbb U_i)_{i=1}^N$.
\bigskip

\point{\bf Main algorithm to compute $\mathbb U$.}\par\penalty10000
\medskip\penalty10000
\subpoint
We rewrite\-~(\ref{bbU_easyeq}) as a root finding problem:
\begin{equation}
  \Xi_i(\mathbb U):=\mathbb U_i-\frac{\mathbb T_i}{2(\mathbb X_i+1)}\Phi\left({\textstyle B_i\frac{\mathbb T_i}{(\mathbb X_i+1)^2}}\right)
  =0
  \label{root_easyeq}
\end{equation}
which we solve using the Newton algorithm, that is, we define a sequence of $\mathbb U$'s:
\begin{equation}
  \mathbb U^{(n+1)}=\mathbb U^{(n)}-(D\Xi(\mathbb U^{(n)}))^{-1}\Xi(\mathbb U^{(n)})
\end{equation}
where $D\Xi$ is the Jacobian of $\Xi$:
\begin{equation}
  (D\Xi)_{i,j}:=\frac{\partial\Xi_i}{\partial\mathbb U_j}
  .
\end{equation}
\bigskip

\subpoint
For small values of $\rho$, we initialize the algorithm with the hardcore scattering solution
\begin{equation}
  \hat u_0(k)=\frac{4\pi a_0}{k^2}
  \label{easyeq_init}
\end{equation}
where $a_0$ is the scattering length of the potential $v$ (or an approximation thereof, which need not be very good).
Thus,
\begin{equation}
  \mathbb U^{(0)}_i=\frac{4\pi a_0}{k_i^2}
  .
\end{equation}
This is a good approximation for small $\rho$.
For larger $\rho$, we choose $\mathbb U^{(0)}$ as the solution of {\tt easyeq} for a slightly smaller $\rho$, and proceed inductively (using the parameters \refname{param:easyeq_minlrho_init}{{\tt minlrho\_init}} and \refname{param:easyeq_nlrho_init}{{\tt nlrho\_init}}).
\bigskip

\subpoint
We are left with computing the Jacobian of $\Xi$:
\begin{equation}
  \begin{largearray}
    \frac{\partial\Xi_i}{\partial\mathbb U_j}
    =\delta_{j,i}
    -\frac1{2(\mathbb X_i+1)}\left(
      \partial_j\mathbb T_i-\frac{\mathbb T_i\partial_j\mathbb X_i}{\mathbb X_i+1}
    \right)\Phi\left({\textstyle\mathbb B_i\frac{\mathbb T_i}{(\mathbb X_i+1)^2}}\right)
    \\[0.5cm]\hfill
    -\frac{\mathbb T_i}{2(\mathbb X_i+1)^3}\left(
      \mathbb B_i\partial_j\mathbb T_i+\mathbb T_i\partial_j\mathbb B_i-2\frac{\mathbb B_i\mathbb T_i\partial_j\mathbb X_i}{\mathbb X_i+1}
    \right)\partial\Phi\left({\textstyle\mathbb B_i\frac{\mathbb T_i}{(\mathbb X_i+1)^2}}\right)
  \end{largearray}
  \label{jacobian_easyeq}
\end{equation}
with
\begin{equation}
  \partial_j\mathbb B_i=
  (\beta_L(1-\beta_K)-\beta_K(1-\beta_L))
  \frac{\mathbb E\partial_j\mathbb S_i-\mathbb S_i\partial_j\mathbb E}{(\beta_K\mathbb S_i+(1-\beta_K)\mathbb E)^2}
\end{equation}
\begin{equation}
  \partial_j\mathbb T_i=
  (1-\beta_K)\frac{\mathbb E\partial_j\mathbb S_i-\mathbb S_i\partial_j\mathbb E}{(\beta_K\mathbb S_i+(1-\beta_K)\mathbb E)^2}
  ,\quad
  \partial_j\mathbb A_{K,i}=\beta_K\partial_j\mathbb S_i+(1-\beta_K)\partial_j\mathbb E
\end{equation}
\begin{equation}
  \partial_j\mathbb S_i:=-\frac1{16\pi^3\rho}w_j\frac{(1-y_j)\Eta(k_i,k_j)}{y_j^3}
  ,\quad
  \partial_j\mathbb E:=-\frac1{16\pi^3\rho}\sum_{j=1}^N w_j\frac{(1-y_j)\Eta(0,k_j)}{y_j^3}
\end{equation}
\begin{equation}
  \partial_j\mathbb X_i:=-\frac{k_i^2}{2\rho\mathbb A_{K,i}^2}\partial_j\mathbb A_{K,i}
  .
\end{equation}
\bigskip

\subpoint
We iterate the Newton algorithm until the Newton relative error $\epsilon$ becomes smaller than the \refname{param:easyeq_tolerance}{{\tt tolerance}} parameter.
The Newton error is defined as
\begin{equation}
  \epsilon=\frac{\|\mathbb U^{(n+1)}-\mathbb U^{(n)}\|_2}{\|\mathbb U^{(n)}\|_2}
\end{equation}
where $\|\cdot\|_2$ is the $l_2$ norm.
The energy thus obtained is
\begin{equation}
  e=\frac\rho2\mathds E
\end{equation}
the Fourier transform $\hat u$ of the solution is
\begin{equation}
  \hat u(k_i)\approx\frac{\mathbb U_i}\rho
\end{equation}
and the solution $u$ in real space is obtained by inverting the Fourier transform
\begin{equation}
  \begin{array}{>\displaystyle r@{\ }>\displaystyle l}
    u(|x|)=\int\frac{dk}{8\pi^3}\ e^{-ikx}\hat u(|k|)
    =&\frac1{2\pi^2|x|}\int_0^1dy\ \frac{(1-y)\sin(|x|\frac{1-y}y)\hat u(\frac{1-y}y)}{y^3}
    \\[0.5cm]
    \approx&
    \frac1{4\pi^2|x|}\sum_{j=1}^Nw_j (1+k_j)^2k_j\sin(|x|k_j)\hat u(k_j)
    .
  \end{array}
  \label{ux}
\end{equation}
To compute $2u-\rho u\ast u$, we replace $\hat u$ with $2\hat u-\rho\hat u^2$ in the previous equation.
\bigskip

\point{\bf Condensate fraction.}
Finally, to compute the uncondensed fraction, we solve the modified {\tt easyeq} (see\-~\cite{CJL20b})
\begin{equation}
  (-\Delta+2\mu)u_\mu=v(1-u_\mu)-2\rho K+\rho^2L
  \label{easyeq_mu}
\end{equation}
where $K$ and $L$ are defined as in\-~(\ref{easyeq_KL})-(\ref{easyeq_Se}) in which $u$ is replaced with $u_\mu$.
The uncondensed fraction is then
\begin{equation}
  \eta=\partial_\mu e|_{\mu=0}
  =-\frac\rho2\int dx\ v(|x|)\partial_\mu u_\mu(|x|)|_{\mu=0}
  .
\end{equation}
To compute the energy in the presence of the parameter $\mu$, we proceed in the same way as for $\mu=0$, the only difference being that $k^2$ should formally be replaced by $k^2+2\mu$.
In other words, we consider $\mathbb U_i=u_\mu(|k_i|)$ and define $\Xi(\mathbb U,\mu)$ in the same way as in\-~(\ref{root_easyeq}), except that $\mathbb X_i$ should be replaced by
\begin{equation}
  \frac{k_i^2+2\mu}{2\rho\mathbb A_{K,i}}
  .
\end{equation}
We then solve
\begin{equation}
  \Xi(\mathbb U,\mu)=0
  .
\end{equation}
By differentiating this identity with respect to $\mu$, we find $\partial_\mu u_\mu$:
\begin{equation}
  D\Xi\partial_\mu \mathbb U=-\partial_\mu\Xi
\end{equation}
and the uncondensed fraction is
\begin{equation}
  \eta
  =\frac12\int dx\ v(x)(D\Xi)^{-1}\partial_\mu\Xi|_{\mu=0}
\end{equation}
thus
\begin{equation}
  \eta=\frac1{16\pi^3}\int_0^1 dy\ \frac{(1-y)\Eta(0,\frac{1-y}y)(D\Xi)^{-1}\partial_\mu\Xi|_{\mu=0}(\frac{1-y}y)}{y^3}
\end{equation}
which we approximate using a Gauss-Legendre quadrature:
\begin{equation}
  \eta\approx\frac1{32\pi^3}\sum_{j=1}^N w_j(1-k_j)^2k_j\Eta(0,k_j)\sum_{i=1}^N(D\Xi)^{-1}_{j,i}\partial_\mu\Xi_i|_{\mu=0}
  .
  \label{eta_easyeq}
\end{equation}
We then compute, using\-~(\ref{jacobian_easyeq}),
\begin{equation}
  \partial_\mu\Xi_i|_{\mu=0}
  =
  \frac{\mathbb T_i}{2\rho\mathbb A_{K,i}(\mathbb X_i+1)^2}\Phi\left({\textstyle\mathbb B_i\frac{\mathbb T_i}{(\mathbb X_i+1)^2}}\right)
  +\frac{\mathbb B_i\mathbb T_i^2}{\rho\mathbb A_{K,i}(\mathbb X_i+1)^4}\partial\Phi\left({\textstyle\mathbb B_i\frac{\mathbb T_i}{(\mathbb X_i+1)^2}}\right)
  .
\end{equation}

\subsection{\tt anyeq}\label{sec:anyeq}
\indent
This method is used to solve any of the equations in the Simplified approach.
Specifically, it solves the equation
\begin{equation}
  -\Delta u
  =
  v(1-u)-2\rho K+\rho^2 L
  \label{anyeq}
\end{equation}
with
\begin{equation}
  K=\gamma_K(1-\alpha_Ku)\left(\beta_K u\ast S+(1-\beta_K)\frac{2e}\rho u\right)
  \label{anyeq_K}
\end{equation}
and
\begin{equation}
  L:=L_1+L_2+L_3
\end{equation}
\begin{equation}
  L_1:=(1-\alpha_{L,1}u)\left(\beta_{L,1}u\ast u\ast S+(1-\beta_{L,1})\frac{2e}\rho u\ast u\right)
\end{equation}
\begin{equation}
  L_2:=-\gamma_{L,2}(1-\alpha_{L,2}u)\left(\beta_{L,2}2u\ast(u(u\ast S))+(1-\beta_{L,2})\frac{4e}\rho u\ast u^2\right)
\end{equation}
\begin{equation}
  L_3:=\gamma_{L,3}(1-\alpha_{L,3}u)\left(\beta_{L,3}\frac12\int dydz\ u(y)u(z-x)u(z)u(y-x)S(z-y)+(1-\beta_{L,3})\frac e\rho u^2\ast u^2\right)
\end{equation}
\begin{equation}
  e=\frac\rho2\int dx\ (1-u(|x|))v(|x|)
  .
  \label{anyeq_e}
\end{equation}
The parameters $\alpha_\cdot$, $\beta_\cdot$ and $\gamma_\cdot$ can be set to turn\-~(\ref{anyeq}) into any of the approximations of the Simplified approach.
For ease of use, there are several predefined equations, given in the following table.
\bigskip

\def\arraystretch{1.5}
\setlength\tabcolsep{5pt}
\makelink{table:anyeq_eqs}{}
\hfil\begin{tabular}{|r|c|c|c|c|c|c|c|c|c|c|c|c|}
  \hline
  &$\alpha_K$ &$\beta_K$ &$\gamma_K$ &$\alpha_{L,1}$ &$\beta_{L,1}$ &$\alpha_{L,2}$ &$\beta_{L,2}$ &$\gamma_{L,2}$ &$\alpha_{L,3}$ &$\beta_{L,3}$ &$\gamma_{L,3}$\\
  \hline
  \makelink{eq:anyeq_compleq}{}{\tt compleq}&1&1&1&1&1&1&1&1&1&1&1\\
  \makelink{eq:anyeq_bigeq}{}  {\tt bigeq  }&1&1&1&1&1&1&1&1&-&-&0\\
  \makelink{eq:anyeq_fulleq}{} {\tt fulleq }&1&1&1&1&1&1&1&1&1&0&1\\
  \makelink{eq:anyeq_medeq}{}  {\tt medeq  }&0&1&1&0&1&-&-&0&-&-&0\\
  \makelink{eq:anyeq_simpleq}{}{\tt simpleq}&0&0&1&0&0&-&-&0&-&-&0\\
  \hline
\end{tabular}
\bigskip

Note that there is no $\gamma_{L,1}$, whose computation would be rather different.
Note, in addition, that {\tt simpleq} and {\tt medeq} coincide with their definitions in\-~(\ref{easyeq}).
The method used to solve this equation is very different from \refname{sec:easyeq}{{\tt easyeq}}, and is significantly longer to run.
\bigskip

\subsubsection{Usage}
\indent
Unless otherwise noted, this method takes the following parameters (specified via the {\tt [-p params]} flag, as a `{\tt;}' separated list of entries).
\begin{itemize}
  \item\makelink{param:anyeq_rho}{}
  {\tt rho} ({\tt Float64}, default: $10^{-6}$): density $\rho$.

  \item\makelink{param:anyeq_tolerance}{}
  {\tt tolerance} ({\tt Float64}, default: $10^{-11}$): maximal size of final step in Newton iteration.

  \item\makelink{param:anyeq_maxiter}{}
  {\tt maxiter} ({\tt Int64}, default: 21): maximal number of iterations of the Newton algorithm before giving up.

  \item\makelink{param:anyeq_P}{}
  {\tt P} ({\tt Int64}, default: 11): order of all Chebyshev polynomial expansions (denoted by $P$ below).

  \item\makelink{param:anyeq_N}{}
  {\tt N} ({\tt Int64}, default: 12): order of all Gauss quadratures (denoted by $N$ below).

  \item\makelink{param:anyeq_J}{}
  {\tt J} ({\tt Int64}, default: 10): number of splines (denoted by $J$ below).

  \item\makelink{param:anyeq_minlrho_init}{}
  {\tt minlrho\_init} ({\tt Float64}, default: \refname{param:anyeq_rho}{{\tt rho}}): we initialize the Newton algorithm using the solution of \refname{eq:easyeq_medeq}{{\tt medeq}}, computed using the methods in \refname{sec:easyeq}{{\tt easyeq}}. This option is passed to the underlying \refname{sec:easyeq}{{\tt easyeq}} routine.

  \item\makelink{param:anyeq_nlrho_init}{}
  {\tt nlrho\_init} ({\tt Int64}, default: 1): this option is passed to the underlying \refname{sec:easyeq}{{\tt easyeq}} routine to initialize the Newton algorithm.

  \item\makelink{param:anyeq_aK}{}
  {\tt aK}, {\tt bK}, {\tt gK}, {\tt aL1}, {\tt bL1}, {\tt aL2}, {\tt bL2}, {\tt gL2}, {\tt aL3}, {\tt bL3}, {\tt gL3} ({\tt Float64}, default: 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0): the values of $\alpha_K$, $\beta_K$, $\gamma_K$, $\alpha_{L,1}$, $\beta_{L,1}$, $\alpha_{L,2}$, $\beta_{L,2}$, $\gamma_{L,2}$, $\alpha_{L,3}$, $\beta_{L,3}$, $\gamma_{L,3}$.

  \item\makelink{param:anyeq_eq}{}
  {\tt eq} ({\tt String}, default: ``{\tt bigeq}'', acceptable values: ``{\tt compleq}'', ``{\tt bigeq}'', ``{\tt fulleq}'', ``{\tt medeq}'', ``{\tt simpleq}''): A shortcut to select any of the equations defined in the \refname{table:anyeq_eqs}{table above}. When this option is set, none of
  \refname{param:anyeq_aK}{\tt aK}, \refname{param:anyeq_aK}{\tt bK}, \refname{param:anyeq_aK}{\tt gK}, \refname{param:anyeq_aK}{\tt aL1}, \refname{param:anyeq_aK}{\tt bL1}, \refname{param:anyeq_aK}{\tt aL2}, \refname{param:anyeq_aK}{\tt bL2}, \refname{param:anyeq_aK}{\tt gL2}, \refname{param:anyeq_aK}{\tt aL3}, \refname{param:anyeq_aK}{\tt bL3}, \refname{param:anyeq_aK}{\tt gL3} should be set.
\end{itemize}
\bigskip

The available {\tt commands} are the following.

\begin{itemize}
  \item\makelink{command:anyeq_energy}{}
  {\tt energy}: compute the energy $e$ at a given $\rho$.\par
  \underline{Output}: [$e$] [Newton error $\epsilon$].

  \item\makelink{command:anyeq_energy_rho}{}
  {\tt energy\_rho}: compute the energy $e$ as a function of $\rho$.
  The Newton algorithm is initialized with the solution of \refname{eq:easyeq_medeq}{{\tt medeq}}.\par
  \underline{Disabled parameters}: \refname{param:anyeq_rho}{{\tt rho}}, \refname{param:anyeq_minlrho_init}{{\tt minlrho\_init}} and \refname{param:anyeq_nlrho_init}{{\tt nlrho\_init}}.\par
  \underline{Extra parameters}:
  \begin{itemize}
    \item\makelink{param:anyeq_minlrho}{}
    {\tt minlrho} ({\tt Float64}, default: $10^{-6}$): minimal value for $\log_{10}\rho$.

    \item\makelink{param:anyeq_maxlrho}{}
    {\tt maxlrho} ({\tt Float64}, default: $10^{2}$): maximal value for $\log_{10}\rho$.

    \item\makelink{param:anyeq_nlrho}{}
    {\tt nlrho} ({\tt Int64}, default: $100$): number of values for $\rho$ (spaced logarithmically).

    \item\makelink{param:anyeq_rhos}{}
    {\tt rhos} ({\tt Array\{Float64\}}, default: $(10^{{\tt minlrho}+\frac{{\tt maxlrho}-{\tt minlrho}}{{\tt nlrho}}n})_n$: list of values for $\rho$, specified as a `{\tt,}' separated list.
  \end{itemize}
  \underline{Output} (one line for each value of $\rho$): [$\rho$] [$e$] [Newton error $\epsilon$].\par
  \underline{Multithread support}: yes, different values of $\rho$ split up among workers.

  \item\makelink{command:anyeq_energy_rho_init_prevrho}{}
  {\tt energy\_rho\_init\_prevrho}: compute the energy $e$ as a function of $\rho$.
  The Newton algorithm is initialized with the solution of \refname{eq:easyeq_medeq}{{\tt medeq}} for the lowest $\rho$, and with the previously computed $\rho$ for the larger densities.\par
  \underline{Disabled parameters}: same as \refname{command:anyeq_energy_rho}{{\tt energy\_rho}}.\par
  \underline{Extra parameters}: same as \refname{command:anyeq_energy_rho}{{\tt energy\_rho}}.\par
  \underline{Output} (one line for each value of $\rho$): [$\rho$] [$e$] [Newton error $\epsilon$].

  \item\makelink{command:anyeq_energy_rho_init_nextrho}{}
  {\tt energy\_rho\_init\_nextrho}: same as \refname{command:anyeq_energy_rho_init_prevrho}{{\tt energy\_rho\_init\_prevrho}} except that the energy is computed for decreasing densities instead of increasing ones.
  The Newton algorithm is initialized with the solution of \refname{eq:easyeq_medeq}{{\tt medeq}} for the largest $\rho$, and with the previously computed $\rho$ for the smaller densities.\par
  \underline{Disabled parameters}: same as \refname{command:anyeq_energy_rho}{{\tt energy\_rho}}.\par
  \underline{Extra parameters}: same as \refname{command:anyeq_energy_rho}{{\tt energy\_rho}}.\par
  \underline{Output} (one line for each value of $\rho$): [$\rho$] [$e$] [Newton error $\epsilon$].

  \item\makelink{command:anyeq_condensate_fraction}{}
  {\tt condensate\_fraction}: compute the uncondensed fraction $\eta$ at a given $\rho$.\par
  \underline{Output}: [$\eta$] [Newton error $\epsilon$].

  \item\makelink{command:anyeq_condensate_fraction_rho}{}
  {\tt condensate\_fraction\_rho}: compute the uncondensed fraction $\eta$ as a function of $\rho$.
  The Newton algorithm is initialized with the solution of \refname{eq:easyeq_medeq}{{\tt medeq}}.\par
  \underline{Disabled parameters}: same as \refname{command:anyeq_energy_rho}{{\tt energy\_rho}}.\par
  \underline{Extra parameters}: same as \refname{command:anyeq_energy_rho}{{\tt energy\_rho}}.\par
  \underline{Output} (one line for each value of $\rho$): [$\rho$] [$\eta$] [Newton error $\epsilon$].\par
  \underline{Multithread support}: yes, different values of $\rho$ split up among workers.

  \item\makelink{command:anyeq_uk}{}
  {\tt uk}: compute the Fourier transform $\hat u(|k|)$.
  The values $|k|$ at which $\hat u$ is computed are those coming from the Gauss quadratures, and cannot be set.\par
  \underline{Output} (one line for each value of $|k|$): [$|k|$] [$\hat u(|k|)$]

  \item\makelink{command:anyeq_ux}{}
  {\tt ux}: compute $u$ as a function of $|x|$.\par
  \underline{Extra parameters}:
  \begin{itemize}
    \item\makelink{param:anyeq_xmin}{}
    {\tt xmin} ({\tt Float64}, default: 0): minimum of the range of $|x|$ to be printed.

    \item\makelink{param:anyeq_xmax}{}
    {\tt xmax} ({\tt Float64}, default: 100): maximum of the range of $|x|$ to be printed.

    \item\makelink{param:anyeq_nx}{}
    {\tt nx} ({\tt Int64}, default: 100): number of points to print (linearly spaced).

  \end{itemize}
  \underline{Output} (one line for each value of $x$): [$|x|$] [$\mathcal Re(u(|x|))$] [$\mathcal Im(u(|x|))$]

  \item\makelink{command:anyeq_momentum_distribution}{}
  {\tt momentum\_distribution}: compute the momentum distribution $\mathfrak M(|k|)$ at a given $\rho$.\par
  \underline{Output} (one line for each value of $|k|$): [$|k|$] [$\mathfrak M(|k|)$]

  \item\makelink{command_anyeq_2pt}{}
  {\tt 2pt}: compute the two-point correlation function $C_2(|x|)$ at a given $\rho$.\par
  \underline{Extra parameters}: same as \refname{command:anyeq_ux}{{\tt ux}}, plus\par
  \begin{itemize}
    \item\makelink{param:anyeq_window_L}{}
    {\tt window\_L} ({\tt Float64}, default: $10^3$): size of the Hann window used to numerically invert the Fourier transform in the computation of the tow-point correlation function, see\-~(\ref{hann}).
  \end{itemize}
  \underline{Output} (one line for each value of $|x|$): [$|x|$] [$C_2(|x|)$]

  \item \makelink{command:anyeq_save_Abar}{}
  {\tt save\_Abar}: compute the matrix $\bar A$. This matrix is used to compute observables for \refname{eq:anyeq_compleq}{{\tt compleq}}. This command is useful to output the value of $\bar A$ to a file once and for all, and use this file to run commands without recomputing $\bar A$.\par
  \underline{Disabled parameters}: \refname{param:anyeq_rho}{{\tt rho}}, \refname{param:anyeq_tolerance}{{\tt tolerance}}, \refname{param:anyeq_maxiter}{{\tt maxiter}}, \refname{param:anyeq_minlrho_init}{{\tt minlrho\_init}}, \refname{param:anyeq_nlrho}{{\tt nlrho}}.\par
  \underline{Output}: [$\bar A$] (the output is not designed to be human-readable; it is obtained through nested {\tt for} loops; for details, see the code).\par
  \underline{Multithread support}: yes, the first indices are split up among workers, which produces $NJ$ jobs.

\end{itemize}

\subsubsection{Description}
\point{\bf Fourier space formulation.}
The computation is carried out in Fourier space.
We take the convention that the Fourier transform of a function $f(|x|)$ is
\begin{equation}
  \hat f(|k|)=\int_{\mathbb R^3} dx\ e^{ikx}f(|x|)
  =\frac{4\pi}{|k|}\int_0^\infty dr\ \frac{\sin(|k|r)}rf(r)
  .
  \label{anyeq_fourier}
\end{equation}
We define a Fourier-space convolution:
\begin{equation}
  \hat f\hat\ast\hat g(|k|):=\int_{\mathbb R^3}\frac{dp}{(2\pi)^3}\ \hat f(|k-p|)\hat g(|p|)
  .
  \label{anyeq_hatast}
\end{equation}
In Fourier space, (\ref{anyeq}) becomes
\begin{equation}
  k^2\hat u(|k|)=\hat S(|k|)-2\rho\hat K(|k|)+\rho^2\hat L(|k|)
\end{equation}
with
\begin{equation}
  \hat S(|k|)=\int_{\mathbb R^3} dx\ e^{ikx}(1-u(|x|))v(|x|)
  =\hat v(k)-\hat u\hat\ast\hat v(k)
\end{equation}
\begin{equation}
  \rho\hat K=
  \gamma_K(\beta_K\rho\hat S+(1-\beta_K)2e)\hat u
  -
  \gamma_K\alpha_K\hat u\hat\ast((\beta_K\rho\hat S+(1-\beta_K)2e)\hat u)
\end{equation}
\begin{equation}
  \rho\hat L_1=
  (\beta_{L,1}\rho\hat S+(1-\beta_{L,1})2e)\hat u^2
  -
  \alpha_{L,1}\hat u\hat\ast((\beta_{L,1}\rho\hat S+(1-\beta_{L,1})2e)\hat u^2)
\end{equation}
\begin{equation}
  \rho\hat L_2=
  -\gamma_{L,2}(\beta_{L,2}2\rho(\hat u\hat\ast(\hat u\hat S))+(1-\beta_{L,2})4e\hat u\hat\ast\hat u)\hat u
  +
  \gamma_{L,2}\alpha_{L,2}\hat u\hat\ast((\beta_{L,2}2\rho(\hat u\hat\ast(\hat u\hat S))+(1-\beta_{L,2})4e\hat u\hat\ast\hat u)\hat u)
\end{equation}
\begin{equation}
  \rho\hat L_3=
  \gamma_{L,3}(1-\beta_{L,3})e(\hat u\hat\ast\hat u)^2
  -
  \gamma_{L,3}\alpha_{L,3}(1-\beta_{L,3})\hat u\hat\ast(e(\hat u\hat\ast\hat u)^2)
  +
  \gamma_{L,3}\beta_{L,3}\hat l_3
  -
  \frac1{2\rho}\gamma_{L,3}\alpha_{L,3}\beta_{L,3}\hat u\hat\ast\hat l_3
\end{equation}
with
\begin{equation}
  \hat l_3(|k|):=
  \frac1{\rho^2}\int\frac{dq}{(2\pi)^3}\frac{dq'}{(2\pi)^3}\ \rho\hat u(|k-q'|)\rho\hat u(|q|)\rho\hat u(|q'|)\rho\hat u(|k-q|)\hat S(|q'-q|)
  \label{l3}
\end{equation}
Therefore,
\begin{equation}
  \rho^2\hat u(|k|)^2\sigma_{L,1}(|k|)
  -\left(\frac{k^2}\rho+2\sigma_K(|k|)+2f_1(|k|)\right)\rho\hat u(|k|)
  +\left(\hat S(|k|)+\frac12f_2(|k|)+G(|k|)\right)
  =
  0
\end{equation}
with
\begin{equation}
  \sigma_K:=\frac1{\rho}\gamma_K(\beta_K\rho\hat S+(1-\beta_K)2e)
  ,\quad
  \sigma_{L,1}:=\frac1{\rho}(\beta_{L,1}\rho\hat S+(1-\beta_{L,1})2e)
\end{equation}
\begin{equation}
  f_1:=\frac1{\rho^2}\gamma_{L,2}(\beta_{L,2}\rho(\rho\hat u\hat\ast(\rho\hat u\hat S))+(1-\beta_{L,2})2e\rho\hat u\hat\ast\rho\hat u)
\end{equation}
\begin{equation}
  f_2:=\gamma_{L,3}(1-\beta_{L,3})\frac{2e}{\rho^3}(\rho\hat u\hat\ast\rho\hat u)^2
  +\gamma_{L,3}\beta_{L,3}\hat l_3
\end{equation}
\begin{equation}
  \begin{largearray}
    G:=
    \frac2{\rho^2}\gamma_K\alpha_K\rho\hat u\hat\ast((\beta_K\rho\hat S+(1-\beta_K)2e)\rho\hat u)
    -
    \frac1{\rho^2}\alpha_{L,1}\rho\hat u\hat\ast((\beta_{L,1}\rho\hat S+(1-\beta_{L,1})2e)\rho\hat u^2)
    \\[0.5cm]\hfill
    +
    \frac2{\rho}\alpha_{L,2}\rho\hat u\hat\ast(f_1\rho\hat u)
    -
    \frac1{2\rho}\alpha_{L,3}\rho\hat u\hat\ast f_2
    .
  \end{largearray}
\end{equation}
Therefore,
\begin{equation}
  \rho\hat u=1+\xi
  -\sqrt{(1+\xi)^2-(1+\zeta)}
  \label{u}
\end{equation}
with
\begin{equation}
  \xi:=\frac1{\sigma_{L,1}}\left(\sigma_K-\sigma_{L,1}+\frac{k^2}{2\rho}+f_1\right)
  ,\quad
  \zeta:=\frac1{\sigma_{L,1}}\left(\hat S-\sigma_{L,1}+\frac12f_2+G\right)
\end{equation}
We rewrite\-~(\ref{u}) as
\begin{equation}
  U=
  \frac{1+\zeta}{2(\xi+1)}\ \Phi\left({\textstyle\frac{1+\zeta}{(\xi+1)^2}}\right)
  \label{anyeq_hatu}
\end{equation}
with
\begin{equation}
  U:=\rho\hat u
\end{equation}
\begin{equation}
  \Phi(x):=\frac{2(1-\sqrt{1-x})}x
  .
\end{equation}
\bigskip

\point{\bf Evaluating integrals.}
To evaluate integrals numerically, we will split integration intervals over splines and use Gauss-Legendre quadratures.
More specifically, to compute integrals of the form
\begin{equation}
  \int \frac{dk}{(2\pi)^3}\ f(U(k),k)
  =
  \frac1{2\pi^2}\int_0^\infty dt\ t^2f(U(t),t)
\end{equation}
we first compactify space to $(-1,1]$ by changing variables to $\tau=\frac{1-t}{1+t}$:
\begin{equation}
  \int \frac{dk}{(2\pi)^3}\ f(U(k),k)
  =
  \frac1{\pi^2}\int_{-1}^1 d\tau\ \frac{(1-\tau)^2}{(1+\tau)^4}f({\textstyle U(\frac{1-\tau}{1+\tau}),\frac{1-\tau}{1+\tau}})
  .
\end{equation}
We then split $(-1,1]$ into $J$ sub-intervals (given by the parameter \refname{param:anyeq_J}{{\tt J}}), called {\it splines}:
\begin{equation}
  (-1,1]=\bigcup_{j=0}^{J-1}(\tau_j,\tau_{j+1}]
  .
\end{equation}
The $\tau$ are taken to be equally spaced, but the code is designed in such a way that this could be changed easily in the future:
\begin{equation}
  \tau_j=-1+\frac{2j}J
  .
\end{equation}
In these terms,
\begin{equation}
  \int \frac{dk}{(2\pi)^3}\ f(U(k),k)
  =
  \sum_{l=0}^{J-1}
  \int_{\tau_l}^{\tau_{l+1}} d\tau\ \frac{(1-\tau)^2}{(1+\tau^4}g(U({\textstyle\frac{1-\tau}{1+\tau},\frac{1-\tau}{1+\tau}})
  .
\end{equation}
We then change variables to $r$:
\begin{equation}
  \tau=
  -\frac{\tau_{l+1}-\tau_l}2\sin({\textstyle\frac{\pi r}2})+\frac{\tau_{l+1}+\tau_l}2
  =:\frac{1-\vartheta_l(r)}{1+\vartheta_l(r)}
  \label{chebyshevvars}
\end{equation}
(the reason for this specific change of variables will become clear at the end of this paragraph and in the next one)
and find
\begin{equation}
  \int \frac{dk}{(2\pi)^3}\ f(U(k),k)
  =
  \sum_{l=0}^{J-1}
  \frac{\tau_{l+1}-\tau_l}{16\pi}\int_{-1}^1 dr\ \cos({\textstyle\frac{\pi r}2})(1+\vartheta_l(r))^2\vartheta_l^2(r)f(U(\vartheta_l(r)),\vartheta_l(r))
  .
\end{equation}
We then approximate the integral using a Gauss-Legendre quadrature of order $N$ (see appendix\-~\ref{appGL}):
\begin{equation}
  \int \frac{dk}{(2\pi)^3}\ f(k)
  \approx
  \sum_{l=0}^{J-1}
  \frac{\tau_{l+1}-\tau_l}{16\pi}
  \sum_{j=1}^Nw_j
  \cos({\textstyle\frac{\pi x_j}2})(1+\vartheta_l(x_j))^2\vartheta_l^2(x_j)f(\mathbb U_{l,j},\vartheta_l(x_j))
  \label{anyeq_integration}
\end{equation}
with
\begin{equation}
  \mathbb U_{l,j}:=U(k_{l,j})
  ,\quad
  k_{l,j}:=\frac{2+(\tau_{l+1}-\tau_l)\sin(\frac{\pi x_j}2)-(\tau_{l+1}+\tau_\alpha)}{2-(\tau_{l+1}-\tau_l)\sin(\frac{\pi x_j}2)+(\tau_{l+1}+\tau_l)}
  .
\end{equation}
The choice of the change of variables\-~(\ref{chebyshevvars}) is made so that $U$ is evaluated at $k_{l,j}$, which are the momenta that appear naturally in the Chebyshev polynomial expansion below\-~(\ref{klj}).
In this way, we can compute arbitrary integrals of functions of $U$ by just computing the values of $\mathbb U_{l,j}$.
\bigskip

\point{\bf Chebyshev polynomial expansion.}
In the case of \refname{sec:easyeq}{{\tt easyeq}}, we saw that using Gauss quadratures reduced the computation to evaluating $U$ at a fixed set of discrete momenta.
This is not the case here, due to the presence of more complicated convolutions of $U$.
Instead, we will approximate $U$ using polynomial approximations.
We will now discuss this approximation for an arbitrary function $a(|k|)$, as we will need it for other functions than simply $U$.
\bigskip

\subpoint
As in the previous paragraph, we compactify $[0,\infty)$ to $(-1,1]$ via the change of variables $r\mapsto\frac{1-r}{1+r}$, and split the interval into splines.
Inside each spline, we use a Chebyshev polynomial expansion.
However, we need to be careful in the first spline, which encapsulates the behavior at $|k|\to\infty$: $U$ decays to 0 as $|k|\to\infty$, and this behavior cannot be reproduced by a polynomial expansion.
To avoid this, if $a\sim|k|^{-\nu_a}$, we expand $|k|^{\nu_a}a$ instead of $a$.
Actually, to simplify expressions in the presence of the compactification, we will expand $2^{-\nu_a}(1+|k|)^{\nu_a}a$, which makes no conceptual difference.
In the case of $a=U$, we will use $\nu=2$, which we know holds for the Simple equation\-~\cite{CJL20}.
Putting all this together, we write, for $\tau\in(-1,1]$, if $|k|\equiv\frac{1-\tau}{1+\tau}$,
\begin{equation}
  2^{-\nu_a}(1+|k|)^{\nu_a}a(|k|)\equiv
  \frac1{(1+\tau)^{\nu_a}}a({\textstyle\frac{1-\tau}{1+\tau}})
  =\sum_{l=0}^{J-1}\mathds 1_{\tau_l<\tau\leqslant\tau_{l+1}}\sum_{n=0}^\infty \mathbf F_{l,n}^{(\nu_a)}(a) T_n({\textstyle\frac{2\tau-(\tau_l+\tau_{l+1})}{\tau_{l+1}-\tau_l}})
  \label{a_chebyshev_spline}
\end{equation}
in which $\mathds 1_{\tau_l<\tau\leqslant\tau_{l+1}}\in\{0,1\}$ is equal to 1 if and only if $\tau_l<\tau\leqslant\tau_{l+1}$, and
\begin{equation}
  \mathbf F_{l,n}^{(\nu_a)}(a):=
  \frac{2-\delta_{n,0}}{\pi}
  \int_0^{\pi} d\theta\ 
  \frac{\cos(n\theta)}{(1+\frac{\tau_{l+1}-\tau_l}2\cos(\theta)+\frac{\tau_{l+1}+\tau_l}2)^{\nu_a}}
  a({\textstyle\frac{2-(\tau_{l+1}-\tau_l)\cos(\theta)-(\tau_{l+1}+\tau_\alpha)}{2+(\tau_{l+1}-\tau_l)\cos(\theta)+(\tau_{l+1}+\tau_l)}})
  \label{bfF}
\end{equation}
(see appendix\-~\ref{appChebyshev} for a discussion of the Chebyshev polynomial expansion and the error of its truncations).
\bigskip

\subpoint
In order to compute an approximation for $a$ using\-~(\ref{a_chebyshev_spline}), we will truncate the sum over $n$ to a finite value $P$ (given by the parameter \refname{param:anyeq_P}{{\tt P}}).
In addition, to compute the integral in\-~(\ref{bfF}), we will use a Gauss-Legendre quadrature of order $N$ (given by the parameter \refname{param:anyeq_N}{{\tt N}}), see appendix\-~\ref{appGL}:
\begin{equation}
  \mathbf F_{l,n}^{(\nu_a)}(a)\approx\mathfrak F_{l,n}^{(\nu_a)}(\mathfrak a):=
  \frac{2-\delta_{n,0}}2\sum_{j=1}^N w_j \cos({\textstyle\frac{n\pi(1+x_j)}2})\frac{\mathfrak a_{l,j}}{(1-\frac{\tau_{l+1}-\tau_l}2\sin(\frac{\pi x_j}2)+\frac{\tau_{l+1}+\tau_l}2)^{\nu_a}}
  \label{frakF}
\end{equation}
with
\begin{equation}
  \mathfrak a_{l,j}:=a(k_{l,j})
  ,\quad
  k_{l,j}:=\frac{2+(\tau_{l+1}-\tau_l)\sin(\frac{\pi x_j}2)-(\tau_{l+1}+\tau_\alpha)}{2-(\tau_{l+1}-\tau_l)\sin(\frac{\pi x_j}2)+(\tau_{l+1}+\tau_l)}
  \label{klj}
\end{equation}
and $(x_j,w_j)$ are the abcissa and weights for Gauss-Legendre quadratures.
\bigskip

\subpoint
All in all, we approximate
\begin{equation}
  a({\textstyle\frac{1-\tau}{1+\tau}})
  \approx(1+\tau)^{\nu_a}\sum_{l=0}^{J-1}\mathds 1_{\tau_l<\tau\leqslant\tau_{l+1}}\sum_{n=0}^P \mathfrak F_{l,n}^{(\nu_a)}(\mathfrak a) T_n({\textstyle\frac{2\tau-(\tau_l+\tau_{l+1})}{\tau_{l+1}-\tau_l}})
  \label{approxchebyshev}
\end{equation}
with $\mathfrak F$ defined in\-~(\ref{frakF}).
Furthermore, using the Chebyshev polynomial expansion and Gauss-Legendre quadratures, we can compute all the observables we are interested in by computing $\mathbb U_{l,j}\equiv U(k_{l,j})$.
With this in mind, we will represent the function $U$ as a vector of dimension $NJ$ whose components are $\mathbb U_{l,j}$.
\bigskip

\point{\bf Convolutions.}
Using the Chebyshev polynomial approximation, we can compute convolutions as follows.
\medskip

\subpoint
First, we rewrite the convolution as a two-dimensional integral, using bipolar coordinates (see lemma\-~\ref{lemma:bipolar}):
\begin{equation}
  (a\hat\ast b)(|k|)=\frac1{4\pi^2|k|}\int_0^\infty dt\ ta(t)\int_{||k|-t|}^{|k|+t}ds\ sb(s)
\end{equation}
We change variables to compactify the integrals $\tau=\frac{1-t}{1+t}$, $\sigma=\frac{1-s}{1+s}$:
\begin{equation}
  (a\hat\ast b)(|k|)=\frac1{4\pi^2|k|}\int_{-1}^1 d\tau\ \frac{2(1-\tau)}{(1+\tau)^3}a({\textstyle\frac{1-\tau}{1+\tau}})\int_{\alpha_-(|k|,\tau)}^{\alpha_+(|k|,\tau)}d\sigma\ \frac{2(1-\sigma)}{(1+\sigma)^3}b({\textstyle\frac{1-\sigma}{1+\sigma}})
  \label{convvar}
\end{equation}
with
\begin{equation}
  \alpha_-(|k|,\tau):=\frac{1-|k|-\frac{1-\tau}{1+\tau}}{1+|k|+\frac{1-\tau}{1+\tau}}
  ,\quad
  \alpha_+(|k|,\tau):=\frac{1-||k|-\frac{1-\tau}{1+\tau}|}{1+||k|-\frac{1-\tau}{1+\tau}|}
  .
  \label{anyeq_alpha}
\end{equation}
Therefore, using the approximation\-~(\ref{approxchebyshev}), if $\mathfrak a_{l,j}:=a(k_{l,j})$ and $\mathfrak b_{l,j}:=b(k_{l,j})$ (\ref{klj}),
\begin{equation}
  (a\hat\ast b)(|k_{l,j}|)\approx
  (\mathfrak a\odot \mathfrak b)_{l,j}:=
  \frac1{4\pi^2|k_{l,j}|}
  \sum_{n,m=0}^P
  \sum_{l,l'=0}^{J-1}
  \mathfrak F_{l,n}^{(\nu_a)}(\mathfrak a)\mathfrak F_{l',m}^{(\nu_b)}(\mathfrak b)
  A_{l,n;l',m}^{(\nu_a,\nu_b)}(|k_{l,j}|)
  \label{odot}
\end{equation}
with
\begin{equation}
  \begin{largearray}
    A_{l,n;l',m}^{(\nu_a,\nu_b)}(|k|):=
    \int_{\tau_l}^{\tau_{l+1}} d\tau\ \frac{2(1-\tau)}{(1+\tau)^{3-\nu_a}}T_n({\textstyle\frac{2\tau-(\tau_l+\tau_{l+1})}{\tau_{l+1}-\tau_l}})
    \cdot\\[0.5cm]\hfill\cdot
    \mathds 1_{\alpha_-(|\tilde k|,\tau)<\tau_{\zeta'+1}}
    \mathds 1_{\alpha_+(|\tilde k|,\tau)>\tau_{\zeta'}}
    \int_{\max(\tau_{l'},\alpha_-(|k|,\tau))}^{\min(\tau_{l'+1},\alpha_+(|k|,\tau))}d\sigma\ \frac{2(1-\sigma)}{(1+\sigma)^{3-\nu_b}}T_m({\textstyle\frac{2\sigma-(\tau_{l'}+\tau_{l'+1})}{\tau_{l'+1}-\tau_{l'}}})
  \end{largearray}
\end{equation}
(we need the indicator functions to ensure that the bounds of the integral are correct).
Note that $A$ is independent of $a$ and $b$, and can be computed once and for all at the beginning of execution for all values of $k_{l,j}$ (\ref{klj}).
We then compute the integrals using Gauss-Legendre quadratures (as is proved in the next paragraph, the integrand is non singular provided $\nu_a,\nu_b\geqslant 2$).
\bigskip

\subpoint
Note that these integrals are not singular as long as $\nu_a,\nu_b\geqslant 2$: indeed (since the only possible problems occur at $-1$, it suffices to consider the case with only one spline), $\alpha_-,\alpha_+>-1$ for $\tau>-1$, and
\begin{equation}
  \alpha_\pm(k,\tau)=-1+(1+\tau)\pm\frac k2(1+\tau)^2+O(1+\tau)^3
\end{equation}
and
\begin{equation}
  \int_{\alpha_-(|k|,\tau)}^{\alpha_+(|k|,\tau)}d\sigma\ \left|\frac{2(1-\sigma)}{(1+\sigma)^{3-\nu_b}}T_m(\sigma)\right|
  \leqslant
  4\int_{\alpha_-(|k|,\tau)}^{\alpha_+(|k|,\tau)}d\sigma\ \frac1{(1+\sigma)^{3-\nu_b}}
\end{equation}
which, if $\nu_b=2$, yields
\begin{equation}
  4\log\left(\frac{1+\alpha_+(|k|,\tau)}{1+\alpha_-(|k|,\tau)}\right)
  =4|k|(1+\tau)+O(1+\tau)^2
\end{equation}
and if $\nu_b>2$, yields
\begin{equation}
  \frac4{\nu_b-2}\left((1+\alpha_+(|k|,\tau))^{\nu_b-2}-(1+\alpha_-(|k|,\tau))^{\nu_b-2}\right)
  =
  \frac4{\nu_b-2}|k|(1+\tau)^{\nu_b-1}(1+O(1+\tau))
  .
\end{equation}
Therefore,
\begin{equation}
  \left|\frac{2(1-\tau)}{(1+\tau)^{3-\nu_a}}T_n(\tau)\int_{\alpha_-(|k|,\tau)}^{\alpha_+(|k|,\tau)}d\sigma\ \frac{2(1-\sigma)}{(1+\sigma)^{3-\nu_b}}T_m(\sigma)\right|
  \leqslant
  \frac8{c_{\nu_b}}
  |k|(1-\tau)(1+\tau)^{\nu_a+\nu_b-4}(1+O(1+\tau))
\end{equation}
where $c_{2}:=1$ and $c_{\nu_b}:=\nu_b-2$ if $\nu_b>2$.
The integrand is, therefore, not singular as long as $\nu_a,\nu_b\geqslant 2$.
\bigskip

\subpoint
Evaluating convolutions at $k=0$ is not immediate, as the formula for $A(0)$ involves a bit of a computation.
To compute $A(0)$, we expand
\begin{equation}
  \alpha_-(|k|,\tau)=\tau-\frac{|k|(1+\tau)^2}2+O(k^2)
  ,\quad
  \alpha_+(|k|,\tau)=\tau+\frac{|k|(1+\tau)^2}2+O(k^2)
\end{equation}
so
\begin{equation}
  A_{\zeta,n;\zeta',m}^{(\nu_a,\nu_b)}(0)=\mathds 1_{\zeta=\zeta'}
  \frac1{\pi^2}
  \int_{\tau_\zeta}^{\tau_{\zeta+1}} d\tau\ \frac{(1-\tau)^2}{(1+\tau)^{4-\nu_a-\nu_b}}T_n({\textstyle\frac{2\tau-(\tau_\zeta+\tau_{\zeta+1})}{\tau_{\zeta+1}-\tau_\zeta}})
  T_m({\textstyle\frac{2\tau-(\tau_{\zeta}+\tau_{\zeta+1})}{\tau_{\zeta+1}-\tau_{\zeta}}})
\end{equation}
which we then compute using a Gauss-Legendre quadrature.
We can then evaluate convolutions at $k=0$:
\begin{equation}
  (a\hat\ast b)(0)\approx\left<\mathfrak a\mathfrak b\right>
  :=
  \frac1{4\pi^2|k_{l,j}|}
  \sum_{n,m=0}^P
  \sum_{l,l'=0}^{J-1}
  \mathfrak F_{l,n}^{(\nu_a)}(\mathfrak a)\mathfrak F_{l',m}^{(\nu_b)}(\mathfrak b)
  A_{l,n;l',m}^{(\nu_a,\nu_b)}(0)
  \label{avg}
\end{equation}
\bigskip


\subpoint
Let us now compute some choices of $\mathfrak a,\mathfrak b$ more explicitly.
\medskip

\subsubpoint
Let us start with $\mathds 1_{l,n}\hat\ast a$ where $\mathds 1_{l,n}$ is the vector which has 0's everywhere except at position $(l,n)$.
We have
\begin{equation}
  \mathfrak F_{l',m}^{(\nu_1)}(\mathds 1_{l,n})=\delta_{l',l}\frac{2-\delta_{m,0}}2w_n\frac{\cos({\textstyle\frac{m\pi(1+x_n)}2})}{(1-\frac{\tau_{l+1}-\tau_l}2\sin(\frac{\pi x_j}2)+\frac{\tau_{l+1}+\tau_l}2)^{\nu_1}}
\end{equation}
so
\begin{equation}
  (\mathds 1_{n,l''}\odot\mathfrak a)_{m,l}=
  \sum_{l,p}\sum_{l'}\frac{2-\delta_{l,0}}2w_n\frac{\cos({\textstyle\frac{l\pi(1+x_n)}2})}{(1-\frac{\tau_{l''+1}-\tau_{l''}}2\sin(\frac{\pi x_n}2)+\frac{\tau_{l''+1}+\tau_{l''}}2)^{\nu_1}}A_{l'',l;l',p}^{(\nu_1,\nu_a)}(k_{l,m})\mathfrak F_{l',p}^{(\nu_a)}(\mathfrak a)
  .
\end{equation}
\bigskip

\subsubpoint
We now turn to $a\hat\ast\hat v$.
Let
\begin{equation}
  \Upsilon(t,|k|):=
  \int_{||k|-t|}^{|k|+t}ds\ \frac{s\hat v(s)}{|k|}
  \label{Upsilon}
\end{equation}
We have
\begin{equation}
  a\hat\ast\hat v=\frac 1{4\pi^2}\int_0^\infty dt\ ta(t)\Upsilon(t,|k|)
  .
  \label{avUpsilon}
\end{equation}
Following the integration scheme\-~(\ref{anyeq_integration}),
\begin{equation}
  (\mathfrak a\odot\mathfrak v)_{l,i}=
  \frac1{32\pi}\sum_{l'=0}^{J-1}(\tau_{l'+1}-\tau_{l'})\sum_{j=1}^Nw_j \cos({\textstyle\frac{\pi x_j}2})(1+\vartheta_{l'}(x_j))^2\vartheta_{l'}(x_j)\mathfrak a_{l',j}\Upsilon(\vartheta(x_j),k_{l,i})
  .
  \label{adotv}
\end{equation}
At $k=0$,
\begin{equation}
  \left<\mathfrak a\mathfrak v\right>=
  \frac1{32\pi}\sum_{l'=0}^{J-1}(\tau_{l'+1}-\tau_{l'})\sum_{j=1}^Nw_j \cos({\textstyle\frac{\pi x_j}2})(1+\vartheta_{l'}(x_j))^2\vartheta_{l'}(x_j)\mathfrak a_{l',j}\Upsilon(\vartheta(x_j),0)
  \label{<av>}
\end{equation}
with
\begin{equation}
  \Upsilon(t,0)=2tv(t)
  .
\end{equation}
The quantities $\Upsilon(\vartheta(x_n),k_{l,i})$ and $\Upsilon(\vartheta(x_n),0)$ are independent of $a$ and can be computed once and for all at execution.
The integral in\-~(\ref{Upsilon}) is computed using Gauss-Legendre quadratures, but without splitting into splines.
To maintain a high precision, we set the order of the integration to \refname{param:anyeq_J}{{\tt J}}$\times$\refname{param:anyeq_N}{{\tt N}}.
\bigskip

\subsubpoint
Finally,
\begin{equation}
  (\mathds 1_{l',n}\odot\mathfrak v)_{l,i}=
  \frac{\tau_{l'+1}-\tau_{l'}}{32\pi}w_n\cos({\textstyle\frac{\pi x_n}2})(1+\vartheta_{l'}(x_n))^2\vartheta_{l'}(x_n)\Upsilon(\vartheta(x_n),k_{l,i})
\end{equation}
and
\begin{equation}
  \left<\mathds 1_{l',n}\mathfrak v\right>=
  \frac{\tau_{l'+1}-\tau_{l'}}{32\pi}w_n\cos({\textstyle\frac{\pi x_n}2})(1+\vartheta_{l'}(x_n))^2\vartheta_{l'}(x_n)\Upsilon(\vartheta(x_n),0)
  .
\end{equation}
\bigskip

\point{\bf Evaluating $\hat l_3$.}
The only term in\-~(\ref{anyeq}) that does not involve just convolutions (whose computation was described above) is $\hat l_3$ (\ref{l3}).
To evaluate it, we first change variables, using a generalization of bipolar coordinates (see lemma\-~\ref{lemma:bipolar2}):
\begin{equation}
  \begin{largearray}
    \hat l_3(|k|)
    =
    \frac1{(2\pi)^5\rho^2| k|^2}\int_0^\infty dt\int_{|| k|-t|}^{| k|+t}ds\int_0^\infty dt'\int_{|| k|-t'|}^{| k|+t'}ds'\int_0^{2\pi}d\theta\ sts't'U(s)U(t)U(s')U(t')
    \cdot\\\hfill\cdot
    S(\mathfrak R(s,t,s',t',\theta,| k|))
  \end{largearray}
\end{equation}
with
\begin{equation}
  \begin{largearray}
    \mathfrak R(s,t,s',t',\theta,| k|)
    :=
    \frac1{\sqrt 2| k|}
    \Big(
      | k|^2(s^2+t^2+(s')^2+(t')^2)-| k|^4-(s^2-t^2)((s')^2-(t')^2)
      \\[0.5cm]\hfill
      -
      \sqrt{(4| k|^2s^2-(| k|^2+s^2-t^2)^2)(4| k|^2(s')^2-(| k|^2+(s')^2-(t')^2)^2)}\cos\theta
    \Big)^{\frac12}
    .
  \end{largearray}
\end{equation}
We change variables as in\-~(\ref{convvar}):
\begin{equation}
  \begin{largearray}
    \hat l_3(|k|)
    =
    \frac1{(2\pi)^5\rho^2| k|^2}
    \int_{-1}^1 d\tau\ \frac{2(1-\tau)}{(1+\tau)^3}U({\textstyle\frac{1-\tau}{1+\tau}})
    \int_{\alpha_-(| k|,\tau)}^{\alpha_+(| k|,\tau)} d\sigma\ \frac{2(1-\sigma)}{(1+\sigma)^3}U({\textstyle\frac{1-\sigma}{1+\sigma}})
    \int_{-1}^1 d\tau'\ \frac{2(1-\tau')}{(1+\tau')^3}
    \cdot\\[0.5cm]\hfill\cdot
    U({\textstyle\frac{1-\tau'}{1+\tau'}})
    \int_{\alpha_-(| k|,\tau')}^{\alpha_+(| k|,\tau')} d\sigma'\ \frac{2(1-\sigma')}{(1+\sigma')^3}U({\textstyle\frac{1-\sigma'}{1+\sigma'}})
    \int_0^{2\pi}d\theta\  S(\mathfrak R({\textstyle \frac{1-\sigma}{1+\sigma},\frac{1-\tau}{1+\tau},\frac{1-\sigma'}{1+\sigma'},\frac{1-\tau'}{1+\tau'},\theta,| k|}))
    .
  \end{largearray}
\end{equation}
We expand $U$ and $ S$ into Chebyshev polynomials as in\-~(\ref{a_chebyshev_spline}), and split the integrals into splines:
\begin{equation}
  \hat l_3(|k_{l,j}|)
  \approx
  (\mathbb U^{\otimes4}\otimes\mathbb S)_{l,j}:=
  \sum_{\lambda_1,\cdots,\lambda_5=0}^{J-1}
  \sum_{n_1,\cdots,n_5=1}^\infty
  \bar A_{\lambda_1,n_1;\cdots;\lambda_5,n_5}^{(\nu)}(|k_{l,j}|)
  \prod_{i=1}^4
  \left(\mathfrak F_{\lambda_i,n_i}^{(\nu)}(\mathbb U)\right)
  \mathfrak F_{\lambda_5,n_5}^{(\nu)}(\mathbb S)
  \label{otimes}
\end{equation}
with $\mathbb U_{l,j}:=U(k_{l,j})$, $\mathbb S_{l,j}:=S(k_{l,j})$ and
\begin{equation}
  \begin{largearray}
    \bar A_{\lambda_1,n_1;\cdots;\lambda_5,n_5}^{(\nu)}(| k|)
    =
    \frac1{(2\pi)^5\rho^2| k|^2}
    \int_{\tau_{\lambda_1}}^{\tau_{\lambda_1+1}} d\tau\ \frac{2(1-\tau)}{(1+\tau)^{3-\nu}}T_{n_1}({\textstyle\frac{2\tau-(\tau_{\lambda_1}+\tau_{\lambda_1+1})}{\tau_{\lambda_1+1}-\tau_{\lambda_1}}})
    \cdot\\[0.5cm]\hskip1cm\cdot
    \mathds 1_{\alpha_-(| k|,\tau)<\tau_{\lambda_2+1}}\mathds 1_{\alpha_+(| k|,\tau)>\tau_{\lambda_2}}
    \int_{\max(\tau_{\lambda_2},\alpha_-(| k|,\tau))}^{\min(\tau_{\lambda_2+1},\alpha_+(| k|,\tau))} d\sigma\ \frac{2(1-\sigma)}{(1+\sigma)^{3-\nu}}T_{n_2}({\textstyle\frac{2\sigma-(\tau_{\lambda_2}+\tau_{\lambda_2+1})}{\tau_{\lambda_2+1}-\tau_{\lambda_2}}})
    \cdot\\[0.5cm]\hskip1cm\cdot
    \int_{\tau_{\lambda_3}}^{\tau_{\lambda_3+1}} d\tau'\ \frac{2(1-\tau')}{(1+\tau')^{3-\nu}}T_{n_3}({\textstyle\frac{2\tau'-(\tau_{\lambda_3}+\tau_{\lambda_3+1})}{\tau_{\lambda_3+1}-\tau_{\lambda_3}}})
    \cdot\\[0.5cm]\hskip1cm\cdot
    \mathds 1_{\alpha_-(| k|,\tau')<\tau_{\lambda_4+1}}\mathds 1_{\alpha_+(| k|,\tau')>\tau_{\lambda_4}}
    \int_{\max(\tau_{\lambda_4},\alpha_-(| k|,\tau'))}^{\min(\tau_{\lambda_4+1},\alpha_+(| k|,\tau'))} d\sigma'\ \frac{2(1-\sigma')}{(1+\sigma')^{3-\nu}}T_{n_4}({\textstyle\frac{2\sigma'-(\tau_{\lambda_4}+\tau_{\lambda_4+1})}{\tau_{\lambda_4+1}-\tau_{\lambda_4}}})
    \cdot\\[0.5cm]\hskip1cm\cdot
    \int_0^{2\pi}d\theta\ 
    \frac{2^\nu}{(2+\mathfrak R)^\nu}
    T_{n_5}({\textstyle\frac{2\frac{1-\mathfrak R}{1+\mathfrak R}-(\tau_{\lambda_5}+\tau_{\lambda_5+1})}{\tau_{\lambda_5+1}-\tau_{\lambda_5}}})
    \mathds 1_{\tau_{\lambda_5}<\frac{1-\mathfrak R}{1+\mathfrak R}<\tau_{\lambda_5+1}}
  \end{largearray}
\end{equation}
in which $\mathfrak R$ is a shorthand for $\mathfrak R({\textstyle \frac{1-\sigma}{1+\sigma},\frac{1-\tau}{1+\tau},\frac{1-\sigma'}{1+\sigma'},\frac{1-\tau'}{1+\tau'},\theta,| k|})$.
Note that $\bar A$ is independent of $U$, and can be computed once and for all at the beginning of execution.
Since the tensor $\bar A$ is quite large (it contains $(NJ)^5$ entries), and its computation can be rather long it can be inconvenient to compute $\bar A$ at every execution.
Instead, one can use the \refname{command:anyeq_save_Abar}{{\tt save\_Abar}} method to compute $\bar A$ and save it to a file, which can then be recalled via the {\tt -s} option on the command line.
\bigskip

\point{\bf Main algorithm to compute $U$.}
We are now ready to detail how $U\equiv\rho\hat u$ is computed.
All of the observables we are interested in are approximated from the values $\mathbb U_{l,j}:=U(k_{l,j})$ (\ref{klj}).
\bigskip

\subpoint
The equation for $(\mathbb U_{l,j})_{l\in\{0,\cdots,J-1\},j\in\{1,\cdots,N\}}$ is obtained by approximating\-~(\ref{anyeq_hatu}) according to the prescriptions detailed above:
\begin{equation}
  \mathbb U_{l,j}:=
  \frac{1+\mathbb Y_{l,j}}{2(\mathbb X_{l,j}+1)}\ \Phi\left({\textstyle\frac{1+\mathbb Y_{l,j}}{(\mathbb X_{l,j}+1)^2}}\right)
  \label{bbU}
\end{equation}
\begin{equation}
  \mathbb X_{l,j}:=\frac1{\mathbb L_{l,j}}\left(\mathbb K_{l,j}-\mathbb L_{l,j}+\frac{k_{l,j}^2}{2\rho}+\mathbb I_{l,j}\right)
  ,\quad
  \mathbb Y_{l,j}:=\frac1{\mathbb L_{l,j}}\left(\mathbb S_{l,j}-\mathbb L_{l,j}+\frac12\mathbb J_{l,j}+\mathbb G_{l,j}\right)
\end{equation}
\begin{equation}
  \mathbb S_{l,j}:=
  \mathfrak v_{l,j}-\frac1{\rho}(\mathbb U\odot\mathfrak v)_{l,j}
  ,\quad
  \mathfrak v_{l,j}:=\hat v(k_{l,j})
  ,\quad
  \mathbb E:=\hat v(0)-\frac1{\rho}\left<\mathbb U\mathfrak v\right>
\end{equation}
\begin{equation}
  \mathbb I_{l,j}:=\frac1{\rho}\gamma_{L,2}(\beta_{L,2}(\mathbb U\odot(\mathbb U\mathbb S))_{l,j}+(1-\beta_{L,2})\mathbb E(\mathbb U\odot\mathbb U)_{l,j})
\end{equation}
\begin{equation}
  \mathbb K_{l,j}:=\gamma_K(\beta_K\mathbb  S_{l,j}+(1-\beta_K)\mathbb E)
  ,\quad
  \mathbb L_{l,j}:=\gamma_{L,1}(\beta_{L,1}\mathbb  S_{l,j}+(1-\beta_{L,1})\mathbb E)
\end{equation}
\begin{equation}
  \mathbb J_{l,j}:=\gamma_{L,3}(1-\beta_{L,3})\frac{\mathbb E}{\rho^2}(\mathbb U\odot\mathbb U)^2_{l,j}
  +\gamma_{L,3}\beta_{L,3}(\mathbb U^{\otimes4}\otimes\mathbb S)_{l,j}
\end{equation}
\begin{equation}
  \begin{largearray}
    \mathbb G_{l,j}:=
    \frac2{\rho}\gamma_K\alpha_K(\mathbb U\odot((\beta_K\mathbb S+(1-\beta_K)\mathbb E)\mathbb U))_{l,j}
    \\[0.5cm]\hfill
    -
    \frac1{\rho}\alpha_{L,1}(\mathbb U\odot((\beta_{L,1}\mathbb S+(1-\beta_{L,1})\mathbb E)\mathbb U^2))_{l,j}
    +
    \frac2{\rho}\alpha_{L,2}(\mathbb U\odot(\mathbb I\mathbb U))_{l,j}
    -
    \frac1{2\rho}\alpha_{L,3}(\mathbb U\odot \mathbb J)_{l,j}
  \end{largearray}
\end{equation}
(see\-~(\ref{odot}), (\ref{otimes}) and\-~(\ref{avg}) for the definitions of $\odot$, $\otimes$ and $\left<\cdot\right>$).
\bigskip

\subpoint
We rewrite\-~(\ref{bbU}) as a root finding problem:
\begin{equation}
  \Xi_{l,j}(\mathbb U):=
  \mathbb U_{l,j}-
  \frac{1+\mathbb Y_{l,j}}{2(\mathbb X_{l,j}+1)}\ \Phi\left({\textstyle\frac{1+\mathbb Y_{l,j}}{(\mathbb X_{l,j}+1)^2}}\right)
  =0
  \label{root_anyeq}
\end{equation}
which we solve using the Newton algorithm, that is, we define a sequence of $\mathbb U$'s:
\begin{equation}
  \mathbb U^{(n+1)}=\mathbb U^{(n)}-(D\Xi(\mathbb U^{(n)}))^{-1}\Xi(\mathbb U^{(n)})
\end{equation}
where $D\Xi$ is the Jacobian of $\Xi$:
\begin{equation}
  (D\Xi)_{l,j;l',i}:=\frac{\partial\Xi_{l,j}}{\partial\mathbb U_{l',i}}
  .
\end{equation}
\bigskip

\subpoint
We initialize the algorithm with the solution of the \refname{eq:easyeq_medeq}{Medium equation}, which is computed using the \refname{sec:easyeq}{{\tt easyeq}} method.
However, \refname{sec:easyeq}{{\tt easyeq}} only computes $\hat u$ at the momenta given by the Gauss-Legendre quadratures\-~(\ref{easyeq_ki}).
To obtain a value for $\hat u$ at $k_{l,j}$, we use a linear interpolation (code is provided for a polynomial interpolation, which has not performed as well).
The parameters \refname{param:anyeq_tolerance}{{\tt tolerance}}, \refname{param:anyeq_maxiter}{{\tt maxiter}}, \refname{param:anyeq_minlrho_init}{{\tt minlrho\_init}}, \refname{param:anyeq_nlrho_init}{{\tt nlrho\_init}} are passed on to \refname{sec:easyeq}{easyeq} as is, and the order of the Gauss-Legendre quadratures is set to \refname{param:anyeq_J}{{\tt J}}$\times$\refname{param:anyeq_N}{{\tt N}}.
\bigskip

\subpoint
We are left with computing the Jacobian of $\Xi$:
\begin{equation}
  \begin{largearray}
    \partial_{l',i}\Xi_{l,j}\equiv
    \frac{\partial\Xi_{l,j}}{\partial\mathbb U_{l',i}}=
    \delta_{l,l'}\delta_{i,j}
    +\frac1{2(\mathbb X_{l,j}+1)}\left(
      \frac{(1+\mathbb Y_{l,j})\partial_{l',i}\mathbb X_{l,j}}{\mathbb X_{l,j}+1}
      -
      \partial_{l',i}\mathbb Y_{l,j}
    \right)
    \Phi\left({\textstyle\frac{1+\mathbb Y_{l,j}}{(\mathbb X_{l,j}+1)^2}}\right)
    \\[0.5cm]\hfill
    +\frac{(1+\mathbb Y_{l,j})}{2(\mathbb X_{l,j}+1)^3}\left(
      \frac{2(1+\mathbb Y_{l,j})}{\mathbb X_{l,j}+1}\partial_{l',i}\mathbb X_{l,j}
      -\partial_{l',i}\mathbb Y_{l,j}
    \right)\ \partial\Phi\left({\textstyle\frac{1+\mathbb Y_{l,j}}{(\mathbb X_{l,j}+1)^2}}\right)
  \end{largearray}
\end{equation}
with
\begin{equation}
  \partial_{l',i}\mathbb X_{l,j}=
  \frac1{\mathbb L_{l,j}}\left(\partial_{l',i}\mathbb K_{l,j}-\partial_{l',i}\mathbb L_{l,j}+\partial_{l',i}\mathbb I_{l,j}\right)
  -
  \frac{\partial_{l',i}\mathbb L_{l,j}}{\mathbb L_{l,j}}\mathbb X_{l,j}
\end{equation}
\begin{equation}
  \partial_{l',i}\mathbb Y_{l,j}=
  \frac1{\mathbb L_{l,j}}\left(\partial_{l',i}\mathbb S_{l,j}-\partial_{l',i}\mathbb L_{l,j}+\frac12\partial_{l',i}\mathbb J_{l,j}+\partial_{l',i}\mathbb G_{l,j}\right)
  -
  \frac{\partial_{l',i}\mathbb L_{l,j}}{\mathbb L_{l,j}}\mathbb Y_{l,j}
\end{equation}
\begin{equation}
  \partial_{l',i}\mathbb S_{l,j}=
  -\frac1{\rho}(\mathds 1_{l',i}\odot\mathfrak v)_{l,j}
  ,\quad
  \partial_{l',i}\mathbb E=-\frac1{\rho}\left<\mathds 1_{l',i}\mathfrak v\right>
\end{equation}
\begin{equation}
  \partial_{l',i}\mathbb K_{l,j}:=\gamma_K(\beta_K\partial_{l',i}\mathbb  S_{l,j}+(1-\beta_K)\partial_{l',i}\mathbb E)
  ,\quad
  \partial_{l',i}\mathbb L_{l,j}:=\gamma_{L,1}(\beta_{L,1}\partial_{l',i}\mathbb  S_{l,j}+(1-\beta_{L,1})\partial_{l',i}\mathbb E)
\end{equation}
\begin{equation}
  \begin{largearray}
    \partial_{l',i}\mathbb I_{l,j}=
    \frac1{\rho}\gamma_{L,2}\left(\beta_{L,2}(
      \mathds 1_{l',i}\odot(\mathbb U\mathbb S)
      +
      \mathbb U\odot(\mathds 1_{l',i}\mathbb S+\mathbb U\partial_{l',i}\mathbb S)
    )_{l,j}
    \right.\\[0.5cm]\hfill\left.
    +
    (1-\beta_{L,2})(
      \partial_{l',i}\mathbb E(\mathbb U\odot\mathbb U)_{l,j}
      +
      2\mathbb E(\mathbb U\odot\mathds 1_{l',i})_{l,j}
    )\right)
  \end{largearray}
\end{equation}
\begin{equation}
  \begin{largearray}
    \partial_{l',i}\mathbb J_{l,j}=
    \frac1{\rho^2}\gamma_{L,3}(1-\beta_{L,3})\left(
      \partial_{l',i}\mathbb E(\mathbb U\odot\mathbb U)^2_{l,j}
      +
      4\mathbb E(\mathbb U\odot\mathbb U)_{l,j}(\mathbb U\odot\mathds 1_{l',i})_{l,j}
    \right)
    \\[0.5cm]\hfill
    +\gamma_{L,3}\beta_{L,3}(
      4\mathds 1_{l',i}\otimes\mathbb U^{\otimes3}\otimes\mathbb S
      +
      \mathbb U^{\otimes4}\otimes\partial_{l',i}\mathbb S
    )_{l,j}
  \end{largearray}
\end{equation}
\begin{equation}
  \begin{largearray}
    \partial_{l',i}\mathbb G_{l,j}=
    \frac2{\rho}\gamma_K\alpha_K\beta_K\left(
      \mathds 1_{l',i}\odot(\mathbb S\mathbb U)
      +
      \mathbb U\odot(\mathbb S\mathds 1_{l',i}
      +\partial_{l',i}\mathbb S\mathbb U)
    \right)_{l,j}
    \\[0.5cm]\hfill
    +\frac2{\rho}\gamma_K\alpha_K(1-\beta_K)\left(
      \partial_{l',i}\mathbb E(\mathbb U\odot\mathbb U)_{l,j}
      +2\mathbb E(\mathds 1_{l',i}\odot\mathbb U)_{l,j}
    \right)
    \\[0.5cm]\indent
    -
    \frac1{\rho}\alpha_{L,1}\beta_{L,1}\left(
      \mathds 1_{l',i}\odot(\mathbb S\mathbb U^2)
      +
      \mathbb U\odot(
	\partial_{l',i}\mathbb S\mathbb U^2
	+
	2\mathbb S\mathbb U\mathds 1_{l',i}
      )
    \right)_{l,j}
    \\[0.5cm]\hfill
    -
    \frac1{\rho}\alpha_{L,1}(1-\beta_{L,1})\left(
      \mathbb E\mathds 1_{l',i}\odot\mathbb U^2
      +
      \mathbb U\odot(
	\partial_{l',i}\mathbb E\mathbb U^2
	+
	2\mathbb E\mathbb U\mathds 1_{l',i}
      )
    \right)_{l,j}
    \\[0.5cm]\hfill
    +
    \frac2{\rho}\alpha_{L,2}(
      (\mathds 1_{l',i}\odot(\mathbb I\mathbb U))_{l,j}
      +
      \mathbb U\odot(
	\partial_{l',i}\mathbb I\mathbb U
	+
	\mathbb I\mathds 1_{l',i})
    )_{l,j}
    -
    \frac1{2\rho}\alpha_{L,3}(
      \mathds 1_{l',i}\odot \mathbb J
      +
      \mathbb U\odot\partial_{l',i}\mathbb J
    )_{l,j}
    .
  \end{largearray}
\end{equation}
\bigskip

\subpoint
We iterate the Newton algorithm until the Newton relative error $\epsilon$ becomes smaller than the \refname{param:easyeq_tolerance}{{\tt tolerance}} parameter.
The Newton error is defined as
\begin{equation}
  \epsilon=\frac{\|\mathbb U^{(n+1)}-\mathbb U^{(n)}\|_2}{\|\mathbb U^{(n)}\|_2}
\end{equation}
where $\|\cdot\|_2$ is the $l_2$ norm.
The energy thus obtained is
\begin{equation}
  e=\frac\rho2\mathds E
\end{equation}
the Fourier transform $\hat u$ of the solution is
\begin{equation}
  \hat u(k_{l,j})\approx\frac{\mathbb U_{l,j}}\rho
\end{equation}
where $k_{l,j}$ was defined in\-~(\ref{klj}), and the solution $u$ in real space is obtained by inverting the Fourier transform, following the prescription of\-~(\ref{anyeq_integration}):
\begin{equation}
  u(|x|)=\int \frac{dk}{(2\pi)^3}\ e^{-ikx}\hat u(|k|)
  \approx
  \sum_{l=0}^{J-1}
  \frac{\tau_{l+1}-\tau_l}{16\pi}
  \sum_{j=1}^Nw_j
  \cos({\textstyle\frac{\pi x_j}2})(1+\vartheta_l(x_j))^2\vartheta_l(x_j)\mathbb U_{l,j}\sin(\vartheta_l(x_j)|x|)
  .
\end{equation}
\bigskip

\point{\bf Condensate fraction.}
To compute the condensate fraction, we solve the modified {\tt anyeq} (see\-~\cite{CJL20b}):
\begin{equation}
  (-\Delta+2\mu)u_\mu
  =
  (1-u_\mu)v-2\rho K+\rho^2 L
  .
\end{equation}
where $K$ and $L$ are defined as in\-~(\ref{anyeq_K})-(\ref{anyeq_e}) in which $u$ is replaced with $u_\mu$.
The uncondensed fraction is then
\begin{equation}
  \eta=\partial_\mu e|_{\mu=0}
  =-\frac\rho2\int dx\ v(x)\partial_\mu u_\mu(x)|_{\mu=0}
  .
\end{equation}
To compute the energy in the presence of the parameter $\mu$, we proceed in the same way as for $\mu=0$, the only difference being that $k^2$ should formally be replaced by $k^2+2\mu$.
In other words, we consider $\mathbb U_{j,l}=u_\mu(|k_{j,l}|)$ and define $\Xi(\mathbb U,\mu)$ in the same way as in\-~(\ref{root_anyeq}), except that $\mathbb X_i$ should be replaced by
\begin{equation}
  \mathbb X_{l,j}=\frac1{\mathbb L_{l,j}}\left(\mathbb K_{l,j}-\mathbb L_{l,j}+\frac{k_{l,j}^2+2\mu}{2\rho}+\mathbb I_{l,j}\right)
  .
\end{equation}
We then solve
\begin{equation}
  \Xi(\mathbb U_\mu,\mu)=0
\end{equation}
By differentiating this identity with respect to $\mu$, we find $\partial_\mu u_\mu$:
\begin{equation}
  \partial_\mu \mathbb U|_{\mu=0}=-(D\Xi)^{-1}\partial_\mu\Xi|_{\mu=0}
\end{equation}
and
\begin{equation}
  \partial_\mu\Xi|_{\mu=0}=
  \frac{(1+\mathbb Y_{l,j})\partial_\mu\mathbb X_{l,j}}{2(\mathbb X_{l,j}+1)^2}
  \Phi\left({\textstyle\frac{1+\mathbb Y_{l,j}}{(\mathbb X_{l,j}+1)^2}}\right)
  +\frac{(1+\mathbb Y_{l,j})^2}{(\mathbb X_{l,j}+1)^4}\partial_\mu\mathbb X_{l,j}
  \partial\Phi\left({\textstyle\frac{1+\mathbb Y_{l,j}}{(\mathbb X_{l,j}+1)^2}}\right)
  ,\quad
  \partial_\mu\mathbb X_{l,j}=\frac1{\rho\mathbb L_{l,j}}
  .
\end{equation}
We then approximate
\begin{equation}
  \eta\approx
  -\frac12\left<\mathfrak v\partial_\mu\mathbb U\right>
\end{equation}
(see\-~(\ref{avg})).
\bigskip

\point{\bf Correlation function.}
The two-point correlation function is
\begin{equation}
  C_2(x):=
  2\rho\frac{\delta e}{\delta v(x)}
  .
\end{equation}
In Fourier space,
\begin{equation}
  C_2(x)=
  2\rho\int dk\ e^{ikx}\frac{\delta e}{\delta\hat v(k)}
  =2\rho\int dk\ \frac{\sin(|k||x|)}{|k||x|}\frac{\delta e}{\delta\hat v(k)}
  .
  \label{C2fourier}
\end{equation}
\bigskip

\subpoint
We can compute this quantity by considering a modified {\tt anyeq} in Fourier space, by formally replacing $\hat v$ with
\begin{equation}
  \hat v+\lambda g(|k|)
  ,\quad
  g(|k|):=\frac{\sin(|k||x|)}{|k||x|}
  .
\end{equation}
Indeed, if $e_\lambda$ denotes the energy of this modified equation,
\begin{equation}
  \partial_\lambda e_\lambda|_{\lambda=0}
  =\int dk\ \frac{\delta e}{\hat v(k)}\partial_\lambda({\textstyle \hat v(k)+\lambda g(|k|)})
  =\int dk\ g(|k|)\frac{\delta e}{\hat v(k)}
\end{equation}
so, denoting the solution of the modified equation by $u_\lambda$,
\begin{equation}
  C_2(x)=2\rho\partial_\lambda e_\lambda|_{\lambda=0}
  =-\rho^2\int dx\ (g(x)u_0(x)+v(x)\partial_\lambda u_\lambda(x)|_{\lambda=0})
  .
\end{equation}
We compute $\partial_\lambda u_\lambda|_{\lambda=0}$ in the same way as the uncondensed fraction: we define $\Xi(\mathbb U,\lambda)$ by formally adding $\lambda g(|k|)$ to $\hat v$, solve $\Xi(\mathbb U,\lambda)=0$, and differentiate:
\begin{equation}
  \partial_\lambda\mathbb U|_{\lambda=0}=-(D\Xi)^{-1}\partial_\lambda\Xi|_{\lambda=0}
  .
\end{equation}
\bigskip

\subpoint
We compute $\partial_\lambda\Xi|_{\lambda=0}$:
\begin{equation}
  \begin{largearray}
    \partial_\lambda\Xi_{l,j}=
    \delta_{l,l'}\delta_{i,j}
    +\frac1{2(\mathbb X_{l,j}+1)}\left(
      \frac{(1+\mathbb Y_{l,j})\partial_\lambda\mathbb X_{l,j}}{\mathbb X_{l,j}+1}
      -
      \partial_\lambda\mathbb Y_{l,j}
    \right)
    \Phi\left({\textstyle\frac{1+\mathbb Y_{l,j}}{(\mathbb X_{l,j}+1)^2}}\right)
    \\[0.5cm]\hfill
    +\frac{(1+\mathbb Y_{l,j})}{2(\mathbb X_{l,j}+1)^3}\left(
      \frac{2(1+\mathbb Y_{l,j})}{\mathbb X_{l,j}+1}\partial_\lambda\mathbb X_{l,j}
      -\partial_\lambda\mathbb Y_{l,j}
    \right)\ \partial\Phi\left({\textstyle\frac{1+\mathbb Y_{l,j}}{(\mathbb X_{l,j}+1)^2}}\right)
  \end{largearray}
\end{equation}
with
\begin{equation}
  \partial_\lambda\mathbb S_{l,j}=
  \mathfrak g_{l,j}
  -\frac1{\rho}(\mathbb U\odot\mathfrak g)_{l,j}
  ,\quad
  \partial_\lambda\mathbb E=g(0)-\frac1{\rho}\left<\mathbb U\mathfrak g\right>
\end{equation}
\begin{equation}
  \partial_\lambda\mathbb X_{l,j}=
  \frac1{\mathbb L_{l,j}}\left(\partial_\lambda\mathbb K_{l,j}-\partial_\lambda\mathbb L_{l,j}+\partial_\lambda\mathbb I_{l,j}\right)
  -
  \frac{\partial_\lambda\mathbb L_{l,j}}{\mathbb L_{l,j}}\mathbb X_{l,j}
\end{equation}
\begin{equation}
  \partial_\lambda\mathbb Y_{l,j}=
  \frac1{\mathbb L_{l,j}}\left(\partial_\lambda\mathbb S_{l,j}-\partial_\lambda\mathbb L_{l,j}+\frac12\partial_\lambda\mathbb J_{l,j}+\partial_\lambda\mathbb G_{l,j}\right)
  -
  \frac{\partial_\lambda\mathbb L_{l,j}}{\mathbb L_{l,j}}\mathbb Y_{l,j}
\end{equation}
\begin{equation}
  \partial_\lambda\mathbb K_{l,j}:=\gamma_K(\beta_K\partial_\lambda\mathbb  S_{l,j}+(1-\beta_K)\partial_\lambda\mathbb E)
  ,\quad
  \partial_\lambda\mathbb L_{l,j}:=\gamma_{L,1}(\beta_{L,1}\partial_\lambda\mathbb  S_{l,j}+(1-\beta_{L,1})\partial_\lambda\mathbb E)
\end{equation}
\begin{equation}
  \partial_\lambda\mathbb I_{l,j}=
  \frac1{\rho}\gamma_{L,2}\left(
    \beta_{L,2}(\mathbb U\odot(\mathbb U\partial_\lambda\mathbb S)_{l,j}
    +
    (1-\beta_{L,2})\partial_\lambda\mathbb E(\mathbb U\odot\mathbb U)_{l,j}
  \right)
\end{equation}
\begin{equation}
  \partial_\lambda\mathbb J_{l,j}=
  \frac1{\rho^2}\gamma_{L,3}(1-\beta_{L,3})\partial_\lambda\mathbb E(\mathbb U\odot\mathbb U)^2_{l,j}
  +\gamma_{L,3}\beta_{L,3}(\mathbb U^{\otimes4}\otimes\partial_\lambda\mathbb S)_{l,j}
\end{equation}
\begin{equation}
  \begin{largearray}
    \partial_\lambda\mathbb G_{l,j}=
    \frac2{\rho}\gamma_K\alpha_K\beta_K\left(\partial_\lambda\mathbb S\mathbb U\right)_{l,j}
    +\frac2{\rho}\gamma_K\alpha_K(1-\beta_K)\partial_\lambda\mathbb E(\mathbb U\odot\mathbb U)_{l,j}
    -
    \frac1{\rho}\alpha_{L,1}\beta_{L,1}\left(\mathbb U\odot\partial_\lambda\mathbb S\mathbb U^2\right)_{l,j}
    \\[0.5cm]\hfill
    -
    \frac1{\rho}\alpha_{L,1}(1-\beta_{L,1})\partial_\lambda\mathbb E\left(\mathbb U\odot(\mathbb U^2)\right)_{l,j}
    +
    \frac2{\rho}\alpha_{L,2}\mathbb U\odot(\partial_\lambda\mathbb I\mathbb U)_{l,j}
    -
    \frac1{2\rho}\alpha_{L,3}(\mathbb U\odot\partial_\lambda\mathbb J)_{l,j}
    .
  \end{largearray}
\end{equation}
To evaluate $(\mathbb U\odot\mathfrak g)$ and $\left<\mathbb U\mathfrak g\right>$, we proceed as in\-~(\ref{adotv}) and\-~(\ref{<av>}).
To do so, we replace $\hat v$ with $g$ in the computation of $\Upsilon$.
\bigskip

\subpoint
In order to invert the Fourier transform in\-~(\ref{C2fourier}) numerically, we will use a Hann window (see appendix\-~\ref{appendix:hann})
\begin{equation}
  H_L(k):=\mathds 1_{|k|<\frac L2}\cos^2({\textstyle\frac{\pi|k|}{L}})
  .
  \label{hann}
\end{equation}
The parameter $L$ is set using \refname{param:anyeq_window_L}{{\tt window\_L}}.
The computation is changed only in that $g$ is changed to $H_L(k)\frac{\sin(|k||x|)}{|k||x|}$.
\bigskip

\point {\bf Momentum distribution.}
To compute the momentum distribution (see\-~\cite{CHe21}), we add a parameter $\lambda$ to {\tt anyeq}:
\begin{equation}
  -\Delta u_\lambda(|x|)
  =
  (1-u_\lambda(|x|))v(|x|)-2\rho K(|x|)+\rho^2 L(|x|)
  -2\lambda \hat u_0(q)\cos(q\cdot x)
\end{equation}
($\hat u_0\equiv\hat u_\lambda|_{\lambda=0}$).
The momentum distribution is then
\begin{equation}
  \mathcal M(q)=\partial_\lambda e|_{\lambda=0}
  =-\frac\rho2\int dx\ v(x)\partial_\lambda u_\lambda(x)|_{\lambda=0}
  .
\end{equation}
Note that the Fourier transform of $2\lambda\hat u_0(q)\cos(q\cdot x)$ is
\begin{equation}
  -(2\pi)^3\lambda\hat u_0(q)(\delta(q+k)+\delta(q-k))
  .
  \label{fouriermomentum}
\end{equation}
We compute $\partial_\lambda u_\lambda|_{\lambda=0}$ in the same way as the uncondensed fraction.
\bigskip

\subpoint
In order to do so we will discretize momentum space, see\-~(\ref{klj}), and so it is necessary to construct a discrete analog of the delta-functions in\-~(\ref{fouriermomentum}).
The starting point we take is
\begin{equation}
  \int dk f(k)\delta(k-q)
  =f(q)
\end{equation}
so, when approximating the integral according to\-~(\ref{anyeq_integration}), we find
\begin{equation}
  \frac{\pi^2}2\sum_{l=0}^{J-1}
  (\tau_{l+1}-\tau_l)
  \sum_{j=1}^Nw_j
  \cos({\textstyle\frac{\pi x_j}2})(1+\vartheta_l(x_j))^2\vartheta_l^2(x_j)f(\vartheta_l(x_j))\tilde\delta(\vartheta_l(x_j)-k_{l',i})
  \approx
  f(k_{l',i})
\end{equation}
where $\tilde\delta$ is the approximation of the delta-function.
Since
\begin{equation}
  \vartheta_l(x_j)\equiv k_{l,j}
\end{equation}
(see\-~(\ref{chebyshevvars})), we find that the definition of $\tilde\delta$ must be
\begin{equation}
  \tilde\delta_{l,j;l',i}
  :=
  \delta_{l,l'}\delta_{j,i}
  \frac2{\pi^2}
  \left(
    (\tau_{l+1}-\tau_l)
    w_j
    \cos({\textstyle\frac{\pi x_j}2})
    (1+k_{l,j})^2k_{l,j}^2
  \right)^{-1}
  .
\end{equation}
Note that, due to the invariance under rotation, the approximation for $\delta(q+k)$ is equal to that for $\delta(q-k)$.
\bigskip

\subpoint
To compute the momentum distribution at $q=k_{l',i}$, we define $\Xi(\mathbb U,\lambda)$ by formally adding $-2(2\pi)^3\lambda \hat u_0(k_{l',i})\tilde\delta_{l,j;l',i}$ to $\mathbb G_{l,j}$, which corresponds to replacing $\mathbb Y$ with
\begin{equation}
  \mathbb Y_{l,j}=
  \frac1{\mathbb L_{l,j}}\left(\mathbb S_{l,j}-\mathbb L_{l,j}+\frac12\mathbb J_{l,j}+\mathbb G_{l,j}-2(2\pi)^3\lambda\frac{\mathbb U_{l',i}}\rho\tilde\delta_{l,j;l',i}\right)
\end{equation}
Then we solve $\Xi(\mathbb U,\lambda)=0$, and differentiate:
\begin{equation}
  \partial_\lambda\mathbb U|_{\lambda=0}=-(D\Xi)^{-1}\partial_\lambda\Xi|_{\lambda=0}
  .
\end{equation}
Finally,
\begin{equation}
  \partial_\lambda\Xi_{l,j}|_{\lambda=0}
  =
  -\partial_\lambda\mathbb Y_{l,j}|_{\lambda=0}\left(
    \frac1{2(\mathbb X_{l,j}+1)}
    \Phi\left({\textstyle\frac{1+\mathbb Y_{l,j}}{(\mathbb X_{l,j}+1)^2}}\right)
    +\frac{(1+\mathbb Y_{l,j})}{2(\mathbb X_{l,j}+1)^3}
    \partial\Phi\left({\textstyle\frac{1+\mathbb Y_{l,j}}{(\mathbb X_{l,j}+1)^2}}\right)
  \right)
\end{equation}
with
\begin{equation}
  \partial_\lambda\mathbb Y_{l,j}|_{\lambda=0}=
  \frac{2(2\pi)^3\mathbb U_{l',i}\tilde\delta_{l,j;l',i}}{\rho\mathbb L_{l,j}}
  .
\end{equation}

\subsection{{\tt simpleq-Kv}}\label{sec:simpleq-Kv}
\indent
The method is used to compute observables for the simple equation
\begin{equation}
  -\Delta u
  =
  v(1-u)-4eu+2e\rho u\ast u
  ,\quad
  e=\frac\rho2\int dx\ (1-u(|x|))v(|x|)
  .
  \label{simpleq}
\end{equation}
One can show\-~\cite[Theorem\-~1.6]{CJL20b} that the condensate fraction is
\begin{equation}
  \eta=\frac{\rho\int v(x)\mathfrak Ku(x)\ dx}{1-\rho\int v(x)\mathfrak K(2u(x)-\rho u\ast u(x))\ dx}
  \label{simpleq-Kv_eta}
\end{equation}
with
\begin{equation}\label{bg6}
  \mathfrak K  = (-\Delta + v +  4e(1 - \rho u\ast))^{-1}
  .
  \label{Kv}
\end{equation}
Similarly, the two-point correlation function is\-~\cite[(45)]{CHe21}
\begin{equation}
  C_2=
  \rho^2(1-u)+
  \rho^2\frac{\mathfrak Kv(1-u)-2\rho u\ast \mathfrak Kv+\rho^2u\ast u\ast \mathfrak Kv}{1-\rho\int dx\ v(x)\mathfrak K(2u(x)-\rho u\ast u(x))}
  .
  \label{simpleq-Kv_C2}
\end{equation}
Thus, using the fact that $\mathfrak K$ is self-adjoint, we can compute these observables of the simple equation directly from the knowledge of $\mathfrak Kv$.
\bigskip

\subsubsection{Usage}
\indent
The computation uses the same approximation scheme as \refname{sec:anyeq}{{\tt anyeq}}, as well as using the solution of \refname{sec:anyeq}{{\tt anyeq}}.
As such, it takes a similar list of parameters:
\refname{param:anyeq_rho}{{\tt rho}},
\refname{param:anyeq_tolerance}{{\tt tolerance}},
\refname{param:anyeq_maxiter}{{\tt maxiter}},
\refname{param:anyeq_P}{{\tt P}},
\refname{param:anyeq_N}{{\tt N}},
\refname{param:anyeq_J}{{\tt J}},
\refname{param:anyeq_minlrho_init}{{\tt minlrho\_init}},
\refname{param:anyeq_nlrho_init}{{\tt nlrho\_init}}.
\bigskip

The available {\tt commands} are the following.
\bigskip

\begin{itemize}
  \item\makelink{command:simpleq-Kv_Kv}{}
  {\tt Kv}: compute $\mathfrak Kv$ as a function of $|x|$.\par
  \underline{Extra parameters}:
  \begin{itemize}
    \item\makelink{param:simpleq-Kv_xmin}{}
    {\tt xmin} ({\tt Float64}, default: 0): minimum of the range of $|x|$ to be printed.

    \item\makelink{param:simpleq-Kv_xmax}{}
    {\tt xmax} ({\tt Float64}, default: 100): maximum of the range of $|x|$ to be printed.

    \item\makelink{param:simpleq-Kv_nx}{}
    {\tt nx} ({\tt Int64}, default: 100): number of points to print (linearly spaced).

  \end{itemize}
  \underline{Output} (one line for each value of $|x|$): [$|x|$] [$\mathfrak K v$]

  \item\makelink{command:simpleq-Kv_condensate-fraction}{}
  {\tt condensate-fraction}: compute the uncondensed fraction as a function of $\rho$\par
  \underline{Extra parameters}:
  \begin{itemize}
    \item\makelink{param:simpleq-Kv_minlrho}{}
    {\tt minlrho} ({\tt Float64}, default: $10^{-6}$): minimal value for $\log_{10}\rho$.

    \item\makelink{param:simpleq-Kv_maxlrho}{}
    {\tt maxlrho} ({\tt Float64}, default: $10^{2}$): maximal value for $\log_{10}\rho$.

    \item\makelink{param:simpleq-Kv_nlrho}{}
    {\tt nlrho} ({\tt Int64}, default: $100$): number of values for $\rho$ (spaced logarithmically).

    \item\makelink{param:simpleq-Kv_rhos}{}
    {\tt rhos} ({\tt Array\{Float64\}}, default: $(10^{{\tt minlrho}+\frac{{\tt maxlrho}-{\tt minlrho}}{{\tt nlrho}}n})_n$: list of values for $\rho$, specified as a `{\tt,}' separated list.
  \end{itemize}
  \underline{Output} (one line for each value of $\rho$): [$\rho$] [$\eta$]

  \item\makelink{command:simpleq-Kv_2pt}{}
  {\tt 2pt}: compute the two-point correlation as a function of $|x|$.\par
  \underline{Extra parameters}: same as \refname{command:simpleq-Kv_Kv}{{\tt Kv}}.\par
  \underline{Output} (one line for each value of $|x|$): [$|x|$] [$C_2$]
\end{itemize}

\subsubsection{Description}
\indent
In Fourier space\-~(\ref{anyeq_fourier}),
\begin{equation}
  \hat{\mathfrak K}^{-1}f
  :=
  \int dx\ e^{ikx}
  (-\Delta+4e(1-\rho u\ast)+v)f
  =(k^2+4e(1-\rho\hat u(|k|))+\hat v\hat\ast)\hat f
\end{equation}
where $\hat\ast$ is defined in\-~(\ref{anyeq_hatast}).
We follow the same approximation scheme as \refname{sec:anyeq}{{\tt anyeq}}:
\begin{equation}
  \hat{\mathfrak K}^{-1}f(k_{l,j})
  \approx(k_{l,j}^2+4e(1-\mathbb U_{l,j}))\mathfrak f_{l,j}
  +
  (\mathfrak v\odot\mathfrak f)_{l,j}
\end{equation}
with $\mathfrak f_{l,j}:=f(k_{l,j})$, $\mathbb U_{l,j}:=\rho u(|k_{l,j}|)$, $\odot$ is defined in\-~(\ref{odot}), and $k_{l,j}$ is defined in\-~(\ref{klj}).
Therefore, we approximate the operator $\hat{\mathfrak K}^{-1}$ by a matrix:
\begin{equation}
  \hat{\mathfrak K}^{-1}f(k_{l,j})\approx
  \sum_{l'=0}^{J-1}\sum_{i=1}^N M_{l,j;l',i}\mathfrak f_{l',i}
\end{equation}
with, by\-~(\ref{odot}) and\-~(\ref{frakF}),
\begin{equation}
  \begin{largearray}
    M_{l,j;l',i}:=
    \delta_{l',l}\delta_{j,i}(k_{l,j}^2+4e(1-\mathbb U_{l,j}))
    \\\hfill
    +
    \frac1{4\pi^2|k_{l,j}|}\sum_{n,m=0}^P\sum_{l''=0}^{J-1}\mathfrak F_{l'',n}^{(\nu_v)}(\mathfrak v)
    A_{l'',n;l',m}^{(\nu_v,\nu_f)}(|k_{l,j}|)
    \frac{\frac{2-\delta_{m,0}}2w_i\cos(\frac{m\pi(1+x_i)}2)}{(1-\frac{\tau_{l'+1}-\tau_{l'}}2\sin(\frac{\pi x_i}2)+\frac{\tau_{l'+1}+\tau_{l'}}2)^{\nu_f}}
    .
  \end{largearray}
\end{equation}
Defining $\mathds 1_{l',i}$ as the vector whose only non-vanishing component is that indexed by $l',i$ which is equal to 1, we can rewrite
\begin{equation}
  M_{l,j;l',i}:=
  \delta_{l',l}\delta_{j,i}(k_{l,j}^2+4e(1-\mathbb U_{l,j}))
  +
  (\mathfrak v\odot\mathds 1_{l',i})_{l,j}
  .
\end{equation}
Thus
\begin{equation}
  \mathfrak Kv\approx M^{-1}\mathfrak v
  .
\end{equation}
\bigskip

\indent
To compute\-~(\ref{simpleq-Kv_eta}), we write
\begin{equation}
  \eta=\frac{\rho\int\frac{dk}{(2\pi)^3} \hat u(k)\hat{\mathfrak K}\hat v(k)}{1-\rho\int \frac{dk}{(2\pi)^3)}(2\hat u(k)-\rho\hat u^2(k))\hat{\mathfrak K}\hat v(k)}
\end{equation}
which we approximate as
\begin{equation}
  \eta\approx\frac{\left<\mathbb UM^{-1}\mathfrak v\right>}{1-\left<(2\mathbb U-\mathbb U^2)M^{-1}\mathfrak v\right>}
\end{equation}
where $\left<\cdot\right>$ is defined in\-~(\ref{avg}).
We can thus compute $\eta$ using the solution $\mathbb U$ computed by \refname{sec:anyeq}{{\tt anyeq}}.
\bigskip

\indent
To compute\-~(\ref{simpleq-Kv_C2}), we write
\begin{equation}
  C_2=
  \rho^2(1-u(x))+
  \rho^2\frac{\mathfrak Kv(x)(1-u(x))+\int\frac{dk}{(2\pi)^3}(-2\rho\hat u\hat{\mathfrak K}\hat v+\rho^2\hat u^2\hat{\mathfrak K}\hat v)}{1-\rho\int\frac{dk}{(2\pi)^3}\ (2\hat u(k)-\rho \hat u^2(k))\hat{\mathfrak K}\hat v(k)}
\end{equation}
which we approximate as
\begin{equation}
  C_2\approx
  \rho^2-\rho\left<e^{-ik_{l,j}x}\mathbb U\right>+
  \rho^2\frac{\left<M^{-1}\mathfrak v\right>(1-\rho^{-1}\left<e^{-ik_{l,j}x}\mathbb U\right>)+\left<(-2\mathbb U M^{-1}\mathfrak v+\mathbb UM^{-1}\mathfrak v)\right>}{1-\left<(2\mathbb U-\mathbb U^2)M^{-1}\mathfrak v\right>}
  .
\end{equation}
We can thus compute $C_2$ using the solution $\mathbb U$ computed by \refname{sec:anyeq}{{\tt anyeq}}.

\subsection{\tt simpleq-hardcore}\label{sec:simpleq-hardcore}
\indent
This method is used to solve the Simple equation with a hardcore potential:
\begin{equation}
  \left\{\begin{array}{>\displaystyle ll}
    (-\Delta+4e) u(x)=2e\rho u\ast u(x)
    &\mathrm{for\ }|x|>1\\
    u(x)=1
    &\mathrm{for\ }|x|\leqslant1
  \end{array}\right.
  \label{linearhc}
\end{equation}
with
\begin{equation}
  e=-\frac{4\pi\rho\partial u|_{|x|\searrow 1}}{2(1-\frac83\pi\rho+\rho^2\int_{|x|<1}dx\ u\ast u(x))}
  .
  \label{hardcore_energy}
\end{equation}
\bigskip

This equation is solved in $x$-space, and as such is very different from \refname{sec:easyeq}{{\tt easyeq}}, and significantly longer to run.
\bigskip

\subsubsection{Usage}
\indent
Unless otherwise noted, this method takes the following parameters (specified via the {\tt [-p params]} flag, as a `{\tt;}' separated list of entries).
\begin{itemize}
  \item\makelink{param:simpleq-hardcore_rho}{}
  {\tt rho} ({\tt Float64}, default: $10^{-6}$): density $\rho$.

  \item\makelink{param:simpleq-hardcore_tolerance}{}
  {\tt tolerance} ({\tt Float64}, default: $10^{-11}$): maximal size of final step in Newton iteration.

  \item\makelink{param:simpleq-hardcore_maxiter}{}
  {\tt maxiter} ({\tt Int64}, default: 21): maximal number of iterations of the Newton algorithm before giving up.

  \item\makelink{param:simpleq-hardcore_P}{}
  {\tt P} ({\tt Int64}, default: 11): order of all Chebyshev polynomial expansions (denoted by $P$ below).

  \item\makelink{param:simpleq-hardcore_N}{}
  {\tt N} ({\tt Int64}, default: 12): order of all Gauss quadratures (denoted by $N$ below).

  \item\makelink{param:simpleq-hardcore_J}{}
  {\tt J} ({\tt Int64}, default: 10): number of splines (denoted by $J$ below).
\end{itemize}
\bigskip

The available {\tt commands} are the following.

\begin{itemize}
  \item\makelink{command:simpleq-hardcore_energy_rho}{}
  {\tt energy\_rho}: compute the energy $e$ as a function of $\rho$.\par
  \underline{Disabled parameters}: \refname{param:simpleq-hardcore_rho}{{\tt rho}}.\par
  \underline{Extra parameters}:
  \begin{itemize}
    \item\makelink{param:simpleq-hardcore_minlrho}{}
    {\tt minlrho} ({\tt Float64}, default: $10^{-6}$): minimal value for $\log_{10}\rho$.

    \item\makelink{param:simpleq-hardcore_maxlrho}{}
    {\tt maxlrho} ({\tt Float64}, default: $10^{2}$): maximal value for $\log_{10}\rho$.

    \item\makelink{param:simpleq-hardcore_nlrho}{}
    {\tt nlrho} ({\tt Int64}, default: $100$): number of values for $\rho$ (spaced logarithmically).

    \item\makelink{param:simpleq-hardcore_rhos}{}
    {\tt rhos} ({\tt Array\{Float64\}}, default: $(10^{{\tt minlrho}+\frac{{\tt maxlrho}-{\tt minlrho}}{{\tt nlrho}}n})_n$: list of values for $\rho$, specified as a `{\tt,}' separated list.
  \end{itemize}
  \underline{Output} (one line for each value of $\rho$): [$\rho$] [$e$] [Newton error $\epsilon$].\par
  \underline{Multithread support}: yes, different values of $\rho$ split up among workers.

  \item\makelink{command:simpleq-hardcore_condensate_fraction_rho}{}
  {\tt condensate\_fraction\_rho}: compute the uncondensed fraction $\eta$ as a function of $\rho$.\par
  \underline{Disabled parameters}: same as \refname{command:simpleq-hardcore_energy_rho}{{\tt energy\_rho}}.\par
  \underline{Extra parameters}: same as \refname{command:simpleq-hardcore_energy_rho}{{\tt energy\_rho}}.\par
  \underline{Output} (one line for each value of $\rho$): [$\rho$] [$\eta$] [Newton error $\epsilon$].\par
  \underline{Multithread support}: yes, different values of $\rho$ split up among workers.

  \item\makelink{command:simpleq-hardcore_ux}{}
  {\tt ux}: compute $u$ as a function of $|x|$.\par
  \underline{Extra parameters}:
  \begin{itemize}
    \item\makelink{param:simpleq-hardcore_xmin}{}
    {\tt xmin} ({\tt Float64}, default: 0): minimum of the range of $|x|$ to be printed.

    \item\makelink{param:simpleq-hardcore_xmax}{}
    {\tt xmax} ({\tt Float64}, default: 100): maximum of the range of $|x|$ to be printed.

    \item\makelink{param:simpleq-hardcore_nx}{}
    {\tt nx} ({\tt Int64}, default: 100): number of points to print (linearly spaced).

  \end{itemize}
  \underline{Output} (one line for each value of $x$): [$|x|$] [$u(|x|)$]

\end{itemize}

\subsubsection{Description}
\indent
In order to carry out the computation of the solution of\-~(\ref{linearhc}) and compute the condensate fraction at the same time, we will consider the equation with an added parameter $\mu>0$:
\begin{equation}
  (-\Delta+4\epsilon)u=2e\rho u\ast u
  ,\quad
  \epsilon:=e+\frac\mu2
  \label{hardcore_simple}
\end{equation}
for $|x|>1$.
\bigskip

\point{\bf Energy.}
To compute the energy $e$ of this equation, with the extra parameter $\mu$, we consider the limit of the soft sphere potential $\lambda\mathds 1_{|x|<1}$ (see\-~(\ref{easyeq}) with $\beta_K=\beta_L=0$):
\begin{equation}
  (-\Delta+2\mu+4e)u(x)=s_\lambda(x)+2e\rho u\ast u
  ,\quad
  s_\lambda(x):=\lambda(1-u(x))\mathds 1_{|x|\leqslant 1}
  ,\quad
  \frac{2e}\rho=\int dx\ s_\lambda(x)
  .
\end{equation}
Furthermore, since $\partial u$ need not be continuous at $|x|=1$, by integrating $-\Delta u$ over a thin spherical shell of radius 1, we find that, for $|x|\leqslant 1$,
\begin{equation}
  -\Delta u(x)=-\delta(|x|-1)\partial u|_{|x|\searrow1}
\end{equation}
so, formally,
\begin{equation}
  s_\infty(x)=\mathds 1_{|x|\leqslant 1}(2e(2-\rho u\ast u(x))+2\mu)-\delta(|x|-1)\partial u|_{|x|\searrow1}
\end{equation}
and
\begin{equation}
  \frac{2e}\rho=\int dx\ s_\infty(x)
  =
  2e\left(\frac{8\pi}3-\rho\int_{|x|<1}u\ast u(x)\right)+\frac{8\pi}3\mu-4\pi\partial u|_{|x|\searrow1}
  .
\end{equation}
Therefore,
\begin{equation}
  e=\frac{2\pi\rho(\frac23\mu-\partial u|_{|x|\searrow1})}{1-\frac{8\pi}3\rho+\rho^2\int_{|x|<1}dx\ u\ast u(x)}
  .
  \label{emu}
\end{equation}
\bigskip

\point{\bf Integral equation.}
We turn the differential equation in\-~(\ref{hardcore_simple}) into an integral equation.
Let
\begin{equation}
  w(|x|):=|x|u(x)
\end{equation}
we have, for $r>1$,
\begin{equation}
  (-\partial_r^2+4\epsilon)w(r)=2e\rho ru\ast u(r)
  .
\end{equation}
Furthermore, the bounded solution of
\begin{equation}
  (-\partial_r^2+\alpha^2)w(r)=F(r)
  ,\quad
  w(1)=1
\end{equation}
is
\begin{equation}
  w(r)=e^{-\alpha(r-1)}+\int_0^\infty ds\ \frac{F(s+1)}{2\alpha}\left(e^{-\alpha|(r-1)-s|}-e^{-\alpha((r-1)+s)}\right)
  \label{solw}
\end{equation}
so, for $r>1$,
\begin{equation}
  u(r)=
  \frac 1re^{-2\sqrt\epsilon(r-1)}
  +
  \frac{\rho e}{2r\sqrt\epsilon}\int_0^\infty ds\ (s+1)(u\ast u(s+1))\left(e^{-2\sqrt\epsilon|(r-1)-s|}-e^{-2\sqrt\epsilon((r-1)+s)}\right)
  .
\end{equation}
In order to compute the integral more easily, we split it:
\begin{equation}
  \begin{largearray}
    u(r)=
    \frac 1re^{-2\sqrt\epsilon(r-1)}
    +
    \frac{\rho e}{r\sqrt\epsilon}\int_0^{r-1} ds\ (s+1)(u\ast u(s+1))\sinh(2\sqrt\epsilon s)e^{-2\sqrt\epsilon(r-1)}
    \\[0.5cm]\hfill
    +
    \frac{\rho e}{r\sqrt\epsilon}\sinh(2\sqrt\epsilon(r-1))\int_{r-1}^\infty ds\ (s+1)(u\ast u(s+1))e^{-2\sqrt\epsilon s}
    .
  \end{largearray}
\end{equation}
We change variables in the last integral:
\begin{equation}
  \begin{largearray}
    u(r)=
    \frac 1re^{-2\sqrt\epsilon(r-1)}
    +
    \frac{\rho e}{r\sqrt\epsilon}\int_0^{r-1} ds\ (s+1)(u\ast u(s+1))\sinh(2\sqrt\epsilon s)e^{-2\sqrt\epsilon(r-1)}
    \\[0.5cm]\hfill
    +
    \frac{\rho e}{2r\sqrt\epsilon}\left(1-e^{-4\sqrt\epsilon(r-1)}\right)\int_0^\infty d\sigma\ (\sigma+r)(u\ast u(\sigma+r))e^{-2\sqrt\epsilon\sigma}
    .
  \end{largearray}
  \label{uinteq}
\end{equation}
\bigskip

\point{\bf The auto-convolution term.} We split
\begin{equation}
  u(r)=\mathds 1_{r>1}u_+(r-1)+\mathds 1_{r\leqslant 1}
  .
\end{equation}
in terms of which
\begin{equation}
  u\ast u=
  \mathds 1_{r\leqslant 1}\ast\mathds 1_{r\leqslant 1}
  +2\mathds 1_{r\leqslant 1}\ast(u_+(r-1)\mathds 1_{r>1})
  +(\mathds 1_{r>1}u_+(r-1))\ast(u_+(r-1)\mathds 1_{r>1})
  \label{u*u_dcmp}
\end{equation}
In bipolar coordinates (see lemma\-~\ref{lemma:bipolar}),
\begin{equation}
  \mathds 1_{r\leqslant 1}\ast\mathds 1_{r\leqslant 1}(r)
  =
  \frac{2\pi}{r}\int_0^1 dt\ t\int_{|r-t|}^{r+t}ds\ s\mathds 1_{s\leqslant 1}
  \label{1*1}
\end{equation}
and, if $r\geqslant 0$ and $t\geqslant 0$, then
\begin{equation}
  \int_{|r-t|}^{r+t}ds\ s\mathds 1_{s\leqslant 1}
  =
  \mathds 1_{r-1\leqslant t\leqslant r+1}\int_{|r-t|}^{\min(1,r+t)}ds\ s
  =
  \mathds 1_{r-1\leqslant t\leqslant r+1}
  \left(
    \mathds 1_{t\leqslant 1-r}
    2rt
    +
    \mathds 1_{t>1-r}
    \frac{1-(r-t)^2}2
  \right)
  .
  \label{int1}
\end{equation}
Therefore, if $r\geqslant 1$
\begin{equation}
  \mathds 1_{r\leqslant 1}\ast\mathds 1_{r\leqslant 1}(r)
  =
  \mathds 1_{r\leqslant 2}
  \frac{2\pi}{r}\int_{r-1}^1 dt\ t
  \frac{1-(r-t)^2}2
  =
  \mathds 1_{r\leqslant 2}\frac\pi{12}(r-2)^2(r+4)
  \label{u*u_1}
\end{equation}
and if $r<1$,
\begin{equation}
  \mathds 1_{r\leqslant 1}\ast\mathds 1_{r\leqslant 1}(r)
  =
  \frac{2\pi}r\int_0^1dt\ t
  \left(
    \mathds 1_{t\leqslant 1-r}
    2rt
    +
    \mathds 1_{t>1-r}
    \frac{1-(r-t)^2}2
  \right)
\end{equation}
so
\begin{equation}
  \mathds 1_{r\leqslant 1}\ast\mathds 1_{r\leqslant 1}(r)
  =\frac43\pi(1-r)^3
  +\frac\pi{12}r(36-48r+17r^2)
  =
  \frac\pi{12}(r-2)^2(r+4)
  .
\end{equation}
Thus, (\ref{u*u_1}) holds for all $r$.
Furthermore,
\begin{equation}
  2\mathds 1_{r\leqslant 1}\ast(u_+(r-1)\mathds 1_{r>1})
  =
  \frac{4\pi}{r}\int_1^\infty dt\ tu_+(t-1)\int_{|r-t|}^{r+t}ds\ s\mathds 1_{s\leqslant 1}
\end{equation}
so, if $r\geqslant 0$ then, by\-~(\ref{int1}),
\begin{equation}
  2\mathds 1_{r\leqslant 1}\ast(u_+(r-1)\mathds 1_{r>1})(r)
  =
  \frac{2\pi}{r}\int_{\max(1,r-1)}^{r+1} dt\ tu_+(t-1)(1-(r-t)^2)
  .
  \label{u*u_cross}
\end{equation}
Finally, if $r\geqslant 0$ then
\begin{equation}
  (\mathds 1_{r>1}u_+(r-1))\ast(u_+(r-1)\mathds 1_{r>1})(r)
  =
  \frac{2\pi}{r}\int_1^\infty dt\ tu_+(t-1)\int_{\max{(1,|r-t|)}}^{r+t}ds\ su_+(s-1)
  .
  \label{u*u_u}
\end{equation}
Thus, by~\-(\ref{u*u_1}), (\ref{u*u_cross}) and\-~(\ref{u*u_u}), for $r\geqslant -1$,
\begin{equation}
  \begin{largearray}
    u\ast u(1+r)=
    \mathds 1_{r\leqslant 1}\frac\pi{12}(r-1)^2(r+5)
    +
    \frac{2\pi}{r+1}\int_{\max(0,r-1)}^{r+1} dt\ (t+1)u_+(t)(1-(r-t)^2)
    \\[0.5cm]\hfill
    +
    \frac{2\pi}{r+1}\int_0^\infty dt\ (t+1)u_+(t)\int_{\max{(0,|r-t|-1)}}^{r+t+1}ds\ (s+1)u_+(s)
    .
  \end{largearray}
\end{equation}
We then compactify the integrals by changing variables to $\tau=\frac{1-t}{1+t}$ and $\sigma=\frac{1-s}{1+s}$:
\begin{equation}
  \begin{largearray}
    u\ast u(1+r)=
    \mathds 1_{r\leqslant 1}\frac\pi{12}(r-1)^2(r+5)
    +
    \frac{8\pi}{r+1}\int_{-\frac r{2+r}}^{\min(1,\frac{2-r}r)} d\tau\ \frac1{(1+\tau)^3}u_+({\textstyle\frac{1-\tau}{1+\tau}})(1-(r-{\textstyle\frac{1-\tau}{1+\tau}})^2)
    \\[0.5cm]\hfill
    +
    \frac{32\pi}{r+1}\int_{-1}^1 d\tau\ \frac1{(1+\tau)^3}u_+({\textstyle\frac{1-\tau}{1+\tau}})\int_{\alpha_-(1+r,\tau)}^{\min(1,\beta_+(r,\tau))}d\sigma\ \frac1{(1+\sigma)^3}u_+({\textstyle\frac{1-\sigma}{1+\sigma}})
    .
  \end{largearray}
  \label{u*u}
\end{equation}
with
\begin{equation}
  \alpha_-(1+r,\tau):=\frac{-r-\frac{1-\tau}{1+\tau}}{2+r+\frac{1-\tau}{1+\tau}}
  ,\quad
  \beta_+(r,\tau):=\frac{2-|r-\frac{1-\tau}{1+\tau}|}{|r-\frac{1-\tau}{1+\tau}|}
\end{equation}
(note that $\alpha_-$ is the same as defined for fulleq).
Finally, note that, if $\beta_+<1$, that is, if $|r-\frac{1-\tau}{1+\tau}|>1$, then
\begin{equation}
  \beta_+(r,\tau)=\alpha_+({\textstyle |r-\frac{1-\tau}{1+\tau}|-1+\frac{1-\tau}{1+\tau}},\tau)
\end{equation}
where
\begin{equation}
  \alpha_+(r,\tau):=\frac{1-|r-\frac{1-\tau}{1+\tau}|}{1+|r-\frac{1-\tau}{1+\tau}|}
\end{equation}
is the same as is defined for \refname{sec:anyeq}{{\tt anyeq}} (\ref{anyeq_alpha}).
\bigskip

\point{\bf Chebyshev polynomial expansion.} We use the same interpolation as we used in \refname{sec:anyeq}{{\tt anyeq}}: (\ref{a_chebyshev_spline})
\begin{equation}
  \frac1{(1+\tau)^{\nu_u}}u_+({\textstyle\frac{1-\tau}{1+\tau}})=\sum_{l=0}^{J-1}\mathds 1_{\tau_l\leqslant\tau<\tau_{l+1}}\sum_{n=0}^\infty \mathbf F_{l,n}^{(\nu_u)}(u_+) T_n({\textstyle\frac{2\tau-(\tau_l+\tau_{l+1})}{\tau_{l+1}-\tau_l}})
  \label{a_chebyshev_spline}
\end{equation}
($J$ is set by the parameter \refname{param:simpleq-hardcore_J}{{\tt J}})
with
\begin{equation}
  \mathbf F_{l,n}^{(\nu_u)}(u_+):=
  \frac{2-\delta_{n,0}}{\pi}
  \int_0^{\pi} d\theta\ 
  \frac{\cos(n\theta)}{(1+\frac{\tau_{l+1}-\tau_l}2\cos(\theta)+\frac{\tau_{l+1}+\tau_l}2)^{\nu_u}}
  u_+({\textstyle\frac{2-(\tau_{l+1}-\tau_l)\cos(\theta)-(\tau_{l+1}+\tau_\alpha)}{2+(\tau_{l+1}-\tau_l)\cos(\theta)+(\tau_{l+1}+\tau_l)}})
  \label{bfF}
\end{equation}
and we take $\nu_u$ to be the decay exponent of $u$, which we will assume is $\nu_u=4$.
In particular, by\-~(\ref{u*u}),
\begin{equation}
  u\ast u(1+r)=\mathds 1_{r\leqslant 1}B^{(0)}(r)+\sum_l\sum_n\mathbf F_{l,n}^{(\nu_u)}(u_+)B^{(1)}_{l,n}(r)+\sum_{l,l'}\sum_{n,m}\mathbf F_{l,n}^{(\nu_u)}(u_+)\mathbf F_{l',m}^{(\nu_u)}(u_+)B^{(2)}_{l,n;l',m}(r)
  \label{u*u}
\end{equation}
with
\begin{equation}
  B^{(0)}(r):=
  \frac\pi{12}(r-1)^2(r+5)
\end{equation}
\begin{equation}
  B^{(1)}_{l,n}(r):=
  \mathds 1_{\tau_l<\frac{2-r}r}\mathds 1_{\tau_{l+1}>-\frac r{2+r}}
  \frac{8\pi}{r+1}\int_{\max(\tau_l,-\frac r{2+r})}^{\min(\tau_{l+1},\frac{2-r}r)} d\tau\ \frac1{(1+\tau)^{3-\nu_u}}(1-(r-{\textstyle\frac{1-\tau}{1+\tau}})^2)T_n({\textstyle\frac{2\tau-(\tau_l+\tau_{l+1})}{\tau_{l+1}-\tau_l}})
\end{equation}
\begin{equation}
  \begin{largearray}
    B^{(2)}_{l,n;l',m}(r):=
    \frac{32\pi}{r+1}\int_{\tau_l}^{\tau_{l+1}} d\tau\ \frac1{(1+\tau)^{3-\nu_u}}T_n({\textstyle\frac{2\tau-(\tau_l+\tau_{l+1})}{\tau_{l+1}-\tau_l}})
    \mathds 1_{\tau_{l'}<\alpha_+(|r-\frac{1-\tau}{1+\tau}|-\frac{2\tau}{1+\tau},\tau)}\mathds 1_{\tau_{l'+1}>\alpha_-(1+r,\tau)}
    \cdot\\[0.5cm]\hfill\cdot
    \int_{\max(\tau_{l'},\alpha_-(1+r,\tau))}^{\min(\tau_{l'+1},\alpha_+(|r-\frac{1-\tau}{1+\tau}|-\frac{2\tau}{1+\tau},\tau))}d\sigma\ \frac1{(1+\sigma)^{3-\nu_u}}T_m({\textstyle\frac{2\sigma-(\tau_{l'}+\tau_{l'+1})}{\tau_{l'+1}-\tau_{l'}}})
    .
  \end{largearray}
\end{equation}
Thus, by\-~(\ref{uinteq}), for $r>0$,
\begin{equation}
  u_+(r,e,\epsilon)=
  D^{(0)}(r,e,\epsilon)
  +
  \sum_l\sum_n\mathbf F_{l,n}^{(\nu_u)}(u_+)D^{(1)}_{l,n}(r,e,\epsilon)
  +
  \sum_{l,l'}\sum_{n,m}\mathbf F_{l,n}^{(\nu_u)}(u_+)\mathbf F_{l',m}^{(\nu_u)}(u_+)D^{(2)}_{l,n;l',m}(r,e,\epsilon)
  \label{hardcore_u+}
\end{equation}
with
\begin{equation}
  \begin{largearray}
    D^{(0)}(r,e,\epsilon):=
    \frac 1{r+1}e^{-2\sqrt\epsilon r}
    +
    \frac{\rho e}{(r+1)\sqrt\epsilon}\int_0^{\min(1,r)} ds\ (s+1)B^{(0)}(s)\sinh(2\sqrt\epsilon s)e^{-2\sqrt\epsilon r}
    \\[0.5cm]\hfill
    +
    \mathds 1_{r\leqslant 1}
    \frac{\rho e}{2(r+1)\sqrt\epsilon}\left(1-e^{-4\sqrt\epsilon r}\right)\int_0^{1-r} d\sigma\ (\sigma+r+1)B^{(0)}(\sigma+r)e^{-2\sqrt\epsilon\sigma}
  \end{largearray}
  \label{D0}
\end{equation}
\begin{equation}
  \begin{largearray}
    D^{(1)}_{l,n}(r,e,\epsilon):=
    \frac{\rho e}{(r+1)\sqrt\epsilon}\int_0^r ds\ (s+1)B^{(1)}_{l,n}(s)\sinh(2\sqrt\epsilon s)e^{-2\sqrt\epsilon r}
    \\[0.5cm]\hfill
    +
    \frac{\rho e}{2(r+1)\sqrt\epsilon}\left(1-e^{-4\sqrt\epsilon r}\right)\int_0^\infty d\sigma\ (\sigma+r+1)B_{l,n}^{(1)}(\sigma+r)e^{-2\sqrt\epsilon\sigma}
    .
  \end{largearray}
  \label{D1}
\end{equation}
and
\begin{equation}
  \begin{largearray}
    D^{(2)}_{l,n;l',m}(r,e,\epsilon):=
    \frac{\rho e}{(r+1)\sqrt\epsilon}\int_0^r ds\ (s+1)B^{(2)}_{l,n;l',m}(s)\sinh(2\sqrt\epsilon s)e^{-2\sqrt\epsilon r}
    \\[0.5cm]\hfill
    +
    \frac{\rho e}{2(r+1)\sqrt\epsilon}\left(1-e^{-4\sqrt\epsilon r}\right)\int_0^\infty d\sigma\ (\sigma+r+1)B_{l,n;l',m}^{(2)}(\sigma+r)e^{-2\sqrt\epsilon\sigma}
    .
  \end{largearray}
  \label{D2}
\end{equation}
\bigskip

\point{\bf Energy.} We now compute the approximation for the energy, using\-~(\ref{hardcore_energy}).
\bigskip

\subpoint We start with $\partial u|_{|x|\searrow 1}$.
By\-~(\ref{uinteq}),
\begin{equation}
  \partial u|_{|x|\searrow 1}
  =
  -(1+2\sqrt\epsilon)
  +
  2\rho e\int_0^\infty d\sigma\ (\sigma+1)(u\ast u(\sigma+1))e^{-2\sqrt\epsilon\sigma}
\end{equation}
\bigskip
which, by\-~(\ref{u*u}), becomes
\begin{equation}
  \begin{largearray}
    \partial u|_{|x|\searrow 1}
    =
    -(1+2\sqrt\epsilon)
    +
    \gamma^{(0)}(e,\epsilon)
    +
    \sum_l\sum_n\mathbf F_{l,n}^{(\nu_u)}(u_+)\gamma^{(1)}_{l,n}(e,\epsilon)
    +\\\hfill+
    \sum_{l,l'}\sum_{n,m}\mathbf F_{l,n}^{(\nu_u)}(u_+)\mathbf F_{l',m}^{(\nu_u)}(u_+)\gamma^{(2)}_{l,n;l',m}(e,\epsilon)
  \end{largearray}
\end{equation}
with
\begin{equation}
  \gamma^{(0)}(e,\epsilon):=
  2\rho e\int_0^1 d\sigma\ (\sigma+1)B^{(0)}(\sigma)e^{-2\sqrt\epsilon\sigma}
\end{equation}
\begin{equation}
  \gamma^{(1)}_{l,n}(e,\epsilon):=
  2\rho e\int_0^\infty d\sigma\ (\sigma+1)B_{l,n}^{(1)}(\sigma)e^{-2\sqrt\epsilon\sigma}
\end{equation}
and
\begin{equation}
  \gamma^{(2)}_{l,n;l',m}(e,\epsilon):=
  2\rho e\int_0^\infty d\sigma\ (\sigma+1)B_{l,n;l',m}^{(2)}(\sigma)e^{-2\sqrt\epsilon\sigma}
  .
\end{equation}

\subpoint Let us now turn to $\int_{|x|<1}dx\ u\ast u(x)$.
We have
\begin{equation}
  \int_{|x|<1}dx\ u\ast u(x)
  =
  4\pi\int_0^1 dr\ r^2 u\ast u(r)
\end{equation}
so, by\-~(\ref{u*u}),
\begin{equation}
  \int_{|x|<1}dx\ u\ast u(x)
  =
  \bar\gamma^{(0)}(r)
  +
  \sum_l\sum_n\mathbf F_{l,n}^{(\nu_u)}(u_+)\bar\gamma^{(1)}_{l,n}(r)
  +
  \sum_{l,l'}\sum_{n,m}\mathbf F_{l,n}^{(\nu_u)}(u_+)\mathbf F_{l',m}^{(\nu_u)}(u_+)\bar\gamma^{(2)}_{l,n;l',m}(r)
  \label{intusu}
\end{equation}
with
\begin{equation}
  \bar\gamma^{(0)}:=
  4\pi\int_0^1 d\sigma\ \sigma^2B^{(0)}(\sigma-1)
\end{equation}
\begin{equation}
  \bar\gamma^{(1)}_{l,n}:=
  4\pi\int_0^1 d\sigma\ \sigma^2B_{l,n}^{(1)}(\sigma-1)
\end{equation}
and
\begin{equation}
  \bar\gamma^{(2)}_{l,n;l',m}:=
  4\pi\int_0^1 d\sigma\ \sigma^2B_{l,n;l',m}^{(2)}(\sigma-1)
  .
\end{equation}
\bigskip

\subpoint Thus,
\begin{equation}
  \begin{largearray}
    e=2\pi\rho
    \cdot\\\hfill\cdot
    \frac
    {
      C(\epsilon)
      -
      \gamma^{(0)}(e,\epsilon)
      -
      \sum_l\sum_n\mathbf F_{l,n}^{(\nu_u)}(u_+)\gamma^{(1)}_{l,n}(e,\epsilon)
      -
      \sum_{l,l'}\sum_{n,m}\mathbf F_{l,n}^{(\nu_u)}(u_+)\mathbf F_{l',m}^{(\nu_u)}(u_+)\gamma^{(2)}_{l,n;l',m}(e,\epsilon)
    }
    {1-\frac83\pi\rho+\rho^2\left(
      \bar\gamma^{(0)}
      +
      \sum_l\sum_n\mathbf F_{l,n}^{(\nu_u)}(u_+)\bar\gamma^{(1)}_{l,n}
      +
      \sum_{l,l'}\sum_{n,m}\mathbf F_{l,n}^{(\nu_u)}(u_+)\mathbf F_{l',m}^{(\nu_u)}(u_+)\bar\gamma^{(2)}_{l,n;l',m}
    \right)}
    .
  \end{largearray}
\end{equation}
\begin{equation}
  C(\epsilon):=\frac23\mu+(1+2\sqrt\epsilon)
  =\frac43(\epsilon-e)+1+2\sqrt\epsilon
\end{equation}
\bigskip

\point{\bf Newton algorithm.}
In this paragraph, we set $\epsilon=e$, that is, $\mu=0$.
\bigskip

\subpoint
As we did for \refname{sec:anyeq}{{\tt anyeq}}, we discretize the integral in\-~(\ref{bfF}) by using a Gauss-Legendre quadrature, and truncate the sum over Chebyshev polynomials to order \refname{param:simpleq-hardcore_P}{{\tt P}}.
We then reduce the computation to a finite system of equations, whose variables are
\begin{equation}
  e
  ,\quad
  \mathfrak u_{l,j}:=u_+(r_{l,j})
  ,\quad
  r_{l,j}:=
  \frac{2+(\tau_{l+1}-\tau_l)\sin(\frac{\pi x_j}2)-(\tau_{l+1}+\tau_\alpha)}{2-(\tau_{l+1}-\tau_l)\sin(\frac{\pi x_j}2)+(\tau_{l+1}+\tau_l)}
\end{equation}
where $x_j$ are the abcissa for Gauss-Legendre quadratures (see\-~(\ref{klj})).
In other words, we define a vector $\mathbb U$ in dimension \refname{param:simpleq-hardcore_N}{{\tt N}}$\times$\refname{param:simpleq-hardcore_J}{{\tt J}}$+1$, where the first \refname{param:simpleq-hardcore_N}{{\tt N}}$\times$\refname{param:simpleq-hardcore_J}{{\tt J}} terms are $\mathfrak u_{l,j}$ and the last component is $e$.
We then write\-~(\ref{hardcore_u+}) as, for $l\in\{0,\cdots,J-1\}$, $j\in\{1,\cdots,N\}$,
\begin{equation}
  \Xi_{lN+j}(\mathbb U)=0
  ,\quad
  \Xi_{NJ+1}(\mathbb U)=0
\end{equation}
(note that $\Xi_{lN+j}$ corresponds to the pair $(l,j)$)
with
\begin{equation}
  \begin{largearray}
    \Xi_{lN+j}(\mathbb U):=
    -\mathfrak u_{l,j}
    +
    \mathbb D^{(0)}(r_{l,j},e)
    +
    \sum_{l'}\sum_n\mathfrak F_{l',n}^{(\nu_u)}(\mathfrak u)\mathbb D^{(1)}_{l',n}(r_{l,j},e)
    +\\\hfill+
    \sum_{l',l''}\sum_{n,m}\mathfrak F_{l',n}^{(\nu_u)}(\mathfrak u)\mathfrak F_{l'',m}^{(\nu_u)}(\mathfrak u)\mathbb D^{(2)}_{l',n;l'',m}(r_{l,j},e)
  \end{largearray}
\end{equation}
\begin{equation}
  \begin{largearray}
    \Xi_{NJ+1}(\mathbb U):=-e
    +\\\hfill+
    2\pi\rho
    \frac
    {
      C(e)
      -
      \mathfrak g^{(0)}(e)
      -
      \sum_l\sum_n\mathfrak F_{l,n}^{(\nu_u)}(\mathfrak u)\mathfrak g^{(1)}_{l,n}(e)
      -
      \sum_{l,l'}\sum_{n,m}\mathfrak F_{l,n}^{(\nu_u)}(\mathfrak u)\mathfrak F_{l',m}^{(\nu_u)}(\mathfrak u)\mathfrak g^{(2)}_{l,n;l',m}(e)
    }
    {1-\frac83\pi\rho+\rho^2\left(
      \bar{\mathfrak g}^{(0)}
      +
      \sum_l\sum_n\mathfrak F_{l,n}^{(\nu_u)}(\mathfrak u)\bar{\mathfrak g}^{(1)}_{l,n}
      +
      \sum_{l,l'}\sum_{n,m}\mathfrak F_{l,n}^{(\nu_u)}(\mathfrak u)\mathfrak F_{l',m}^{(\nu_u)}(\mathfrak u)\bar{\mathfrak g}^{(2)}_{l,n;l',m}
    \right)}
    .
  \end{largearray}
\end{equation}
where $\mathfrak F$ is defined in\-~(\ref{frakF}), and $\mathbb D^{(i)}$, $\mathfrak g^{(i)}$, $\bar{\mathfrak g}^{(i)}$ are defined like $D^{(i)}$, $\gamma^{(i)}$ and $\bar\gamma^{(i)}$ except that the integrals over bounded intervals are approximated using Gauss-Legendre quadratures, and the integrals from $0$ to $\infty$ are approximated using Gauss-Laguerre quadratures.
Gauss-Legendre and Gauss-Laguerre quadratures and their errors are discussed in appendix\-~\ref{appGL}.
The orders of the quadratures are given by the variable \refname{param:simpleq-hardcore_N}{{\tt N}}.
\bigskip

\subpoint
We then solve $\Xi=0$ using the Newton algorithm, that is, we define a sequence of $\mathbb U$'s:
\begin{equation}
  \mathbb U^{(n+1)}=\mathbb U^{(n)}-(D\Xi(\mathbb U^{(n)}))^{-1}\Xi(\mathbb U^{(n)})
\end{equation}
where $D\Xi$ is the Jacobian of $\Xi$:
\begin{equation}
  (D\Xi)_{\alpha,\beta}:=\frac{\partial\Xi_\alpha}{\partial\mathbb U_\beta}
  .
\end{equation}
We initialize the algorithm with
\begin{equation}
  \mathbb U^{(0)}_{lN+j}=\frac1{(1+r_{l,j}^2)^2}
  ,\quad
  \mathbb U^{(0)}_{JN+1}=\pi\rho
  .
\end{equation}
\bigskip

\subpoint
We are left with computing the Jacobian of $\Xi$:
\begin{equation}
  \begin{largearray}
    \frac{\partial\Xi_{lN+j}(\mathbb U)}{\partial\mathfrak u_{l',i}}=
    -\delta_{l,l'}\delta_{j,i}
    +
    \sum_{l''}\sum_n\mathfrak F_{l'',n}^{(\nu_u)}(\mathds 1_{l',i})\mathbb D^{(1)}_{l'',n}(r_{l,j},e)
    +\\\hfill+
    2\sum_{l'',l'''}\sum_{n,m}\mathfrak F_{l'',n}^{(\nu_u)}(\mathfrak u)\mathfrak F_{l''',m}^{(\nu_u)}(\mathds 1_{l',j})\mathbb D^{(2)}_{l'',n;l''',m}(r_{l,j},e)
  \end{largearray}
\end{equation}
where $\mathds 1_{l',i}$ is a vector whose only non-vanishing entry is the $(l',i)$-th, which is 1,
\begin{equation}
  \begin{largearray}
    \frac{\partial\Xi_{JN+1}(\mathbb U)}{\partial\mathfrak u_{l',i}}
    =\\\hfill=
    2\pi\rho
    \frac
    {
      -
      \sum_l\sum_n\mathfrak F_{l,n}^{(\nu_u)}(\mathds 1_{l',i})\mathfrak g^{(1)}_{l,n}(e)
      -
      2\sum_{l,l''}\sum_{n,m}\mathfrak F_{l,n}^{(\nu_u)}(\mathfrak u)\mathfrak F_{l'',m}^{(\nu_u)}(\mathds 1_{l',i})\mathfrak g^{(2)}_{l,n;l'',m}(e)
    }
    {1-\frac83\pi\rho+\rho^2\left(
      \bar{\mathfrak g}^{(0)}
      +
      \sum_l\sum_n\mathfrak F_{l,n}^{(\nu_u)}(\mathfrak u)\bar{\mathfrak g}^{(1)}_{l,n}
      +
      \sum_{l,l''}\sum_{n,m}\mathfrak F_{l,n}^{(\nu_u)}(\mathfrak u)\mathfrak F_{l'',m}^{(\nu_u)}(\mathfrak u)\bar{\mathfrak g}^{(2)}_{l,n;l'',m}
    \right)}
    \\[1cm]\hfill
    -(\Xi_{NJ+1}(\mathbb U)+e)
    \frac
    {\rho^2\left(
      \sum_l\sum_n\mathfrak F_{l,n}^{(\nu_u)}(\mathds 1_{l'.i})\bar{\mathfrak g}^{(1)}_{l,n}
      +
      2\sum_{l,l''}\sum_{n,m}\mathfrak F_{l,n}^{(\nu_u)}(\mathfrak u)\mathfrak F_{l'',m}^{(\nu_u)}(\mathds 1_{l',i})\bar{\mathfrak g}^{(2)}_{l,n;l'',m}
    \right)}{1-\frac83\pi\rho+\rho^2\left(
      \bar{\mathfrak g}^{(0)}
      +
      \sum_l\sum_n\mathfrak F_{l,n}^{(\nu_u)}(\mathfrak u)\bar{\mathfrak g}^{(1)}_{l,n}
      +
      \sum_{l,l''}\sum_{n,m}\mathfrak F_{l,n}^{(\nu_u)}(\mathfrak u)\mathfrak F_{l'',m}^{(\nu_u)}(\mathfrak u)\bar{\mathfrak g}^{(2)}_{l,n;l'',m}
    \right)}
    .
  \end{largearray}
\end{equation}
\begin{equation}
  \begin{largearray}
    \frac{\partial\Xi_{lN+j}(\mathbb U)}{\partial e}=
    \partial_e\mathbb D^{(0)}(r_{l,j},e)
    +
    \sum_{l'}\sum_n\mathfrak F_{l',n}^{(\nu_u)}(\mathfrak u)\partial_e\mathbb D^{(1)}_{l',n}(r_{l,j},e)
    +\\\hfill+
    \sum_{l',l''}\sum_{n,m}\mathfrak F_{l',n}^{(\nu_u)}(\mathfrak u)\mathfrak F_{l'',m}^{(\nu_u)}(\mathfrak u)\partial_e\mathbb D^{(2)}_{l',n;l'',m}(r_{l,j},e)
  \end{largearray}
\end{equation}
\begin{equation}
  \begin{largearray}
    \frac{\partial\Xi_{NJ+1}(\mathbb U)}{\partial e}
    =
    -1
    +\\\hfill+
    2\pi\rho
    \frac
    {
      \partial_e C(e)
      -
      \partial_e\mathfrak g^{(0)}(e)
      -
      \sum_l\sum_n\mathfrak F_{l,n}^{(\nu_u)}(\mathfrak u)\partial_e\mathfrak g^{(1)}_{l,n}(e)
      -
      \sum_{l,l'}\sum_{n,m}\mathfrak F_{l,n}^{(\nu_u)}(\mathfrak u)\mathfrak F_{l',m}^{(\nu_u)}(\mathfrak u)\partial_e\mathfrak g^{(2)}_{l,n;l',m}(e)
    }
    {1-\frac83\pi\rho+\rho^2\left(
      \bar{\mathfrak g}^{(0)}
      +
      \sum_l\sum_n\mathfrak F_{l,n}^{(\nu_u)}(\mathfrak u)\bar{\mathfrak g}^{(1)}_{l,n}
      +
      \sum_{l,l'}\sum_{n,m}\mathfrak F_{l,n}^{(\nu_u)}(\mathfrak u)\mathfrak F_{l',m}^{(\nu_u)}(\mathfrak u)\bar{\mathfrak g}^{(2)}_{l,n;l',m}
    \right)}
    .
  \end{largearray}
\end{equation}
To compute $\partial_e\mathbb D^{(i)}$ and $\partial_e\mathfrak g^{(i)}$, we use $\partial_e=\frac1{2\sqrt e}\frac\partial{\partial\sqrt e}$ and
\begin{equation}
  \frac{\partial C(e)}{\partial\sqrt e}=2
\end{equation}
\begin{equation}
  \begin{largearray}
    \frac{\partial D^{(0)}(r,e)}{\partial\sqrt e}
    =
    -\frac{2r}{r+1}e^{-2\sqrt er}
    \\[0.5cm]\hfill
    +
    \frac{\rho}{r+1}\int_0^{\min(1,r)} ds\ (s+1)B^{(0)}(s)e^{-2\sqrt er}\left(
      (1-2\sqrt er)\sinh(2\sqrt e s)
      +2\sqrt es\cosh(2\sqrt e s)
    \right)
    \\[0.5cm]
    +
    \mathds 1_{r\leqslant 1}
    \frac{\rho}{2(r+1)}\int_0^{1-r} d\sigma\ (\sigma+r+1)B^{(0)}(\sigma+r)
    \cdot\\[0.5cm]\hfill\cdot
    \left(
      (1-2\sqrt e\sigma)\left(1-e^{-4\sqrt er}\right)e^{-2\sqrt e\sigma}
      +
      4\sqrt e re^{-2\sqrt e(2r+\sigma)}
    \right)
  \end{largearray}
\end{equation}
\begin{equation}
  \begin{largearray}
    \frac{\partial D_{l,n}^{(1)}(r,e)}{\partial\sqrt e}
    =
    \frac{\rho}{r+1}\int_0^{r} ds\ (s+1)B_{l,n}^{(1)}(s)e^{-2\sqrt er}\left(
      (1-2\sqrt er)\sinh(2\sqrt e s)
      +2\sqrt es\cosh(2\sqrt e s)
    \right)
    \\[0.5cm]
    +
    \frac{\rho}{2(r+1)}\int_0^\infty d\sigma\ (\sigma+r+1)B_{l,n}^{(1)}(\sigma+r)
    \cdot\\[0.5cm]\hfill\cdot
    \left(
      (1-2\sqrt e\sigma)\left(1-e^{-4\sqrt er}\right)e^{-2\sqrt e\sigma}
      +
      4\sqrt e re^{-2\sqrt e(2r+\sigma)}
    \right)
  \end{largearray}
\end{equation}
\begin{equation}
  \begin{largearray}
    \frac{\partial D_{l,n;l',m}^{(2)}(r,e)}{\partial\sqrt e}
    \\[0.5cm]
    =
    \frac{\rho}{r+1}\int_0^{r} ds\ (s+1)B_{l,n;l',m}^{(2)}(s)e^{-2\sqrt er}\left(
      (1-2\sqrt er)\sinh(2\sqrt e s)
      +2\sqrt es\cosh(2\sqrt e s)
    \right)
    \\[0.5cm]
    +
    \frac{\rho}{2(r+1)}\int_0^\infty d\sigma\ (\sigma+r+1)B_{l,n;l',m}^{(2)}(\sigma+r)
    \cdot\\[0.5cm]\hfill\cdot
    \left(
      (1-2\sqrt e\sigma)\left(1-e^{-4\sqrt er}\right)e^{-2\sqrt e\sigma}
      +
      4\sqrt e re^{-2\sqrt e(2r+\sigma)}
    \right)
    .
  \end{largearray}
\end{equation}
Furthermore,
\begin{equation}
  \frac{\partial\gamma^{(0)}(e)}{\partial\sqrt e}=
  4\rho\sqrt e\int_0^1 d\sigma\ (\sigma+1)B^{(0)}(\sigma)(1-\sqrt e\sigma)e^{-2\sqrt e\sigma}
\end{equation}
\begin{equation}
  \frac{\partial\gamma^{(1)}_{l,n}(e)}{\partial\sqrt e}=
  4\rho e\int_0^\infty d\sigma\ (\sigma+1)B_{l,n}^{(1)}(\sigma)(1-\sqrt e\sigma)e^{-2\sqrt e\sigma}
\end{equation}
and
\begin{equation}
  \frac{\partial\gamma^{(2)}_{l,n;l',m}(e)}{\partial\sqrt e}=
  4\rho e\int_0^\infty d\sigma\ (\sigma+1)B_{l,n;l',m}^{(2)}(\sigma)(1-\sqrt e\sigma)e^{-2\sqrt e\sigma}
  .
\end{equation}
Finally, to get from $D$ to $\mathbb D$ and $\gamma$ to $\mathfrak g$, we approximate the integrate using Gauss-Legendre and Gauss-Laguerre quadratures (see appendix\-~\ref{appGL}), as described above.
\bigskip

\subpoint
We iterate the Newton algorithm until the Newton relative error $\epsilon$ becomes smaller than the \refname{param:easyeq_tolerance}{{\tt tolerance}} parameter.
The Newton error is defined as
\begin{equation}
  \epsilon=\frac{\|\mathbb U^{(n+1)}-\mathbb U^{(n)}\|_2}{\|\mathbb U^{(n)}\|_2}
\end{equation}
where $\|\cdot\|_2$ is the $l_2$ norm.
The energy thus obtained is
\begin{equation}
  e=\mathbb U_{JN+1}
  .
\end{equation}
\bigskip

\point{\bf Condensate fraction.}
To compute the condensate fraction, we use the parameter $\mu$ in\-~(\ref{hardcore_simple}).
The uncondensed fraction is
\begin{equation}
  \eta=\partial_\mu e|_{\mu=0}
  .
\end{equation}
To compute $\partial_\mu e$, we use
\begin{equation}
  \Xi(\mathbb U)=0
\end{equation}
which we differentiate with respect to $\mu$:
\begin{equation}
  \partial_\mu\mathbb U=-(D\Xi)^{-1}\partial_\mu\Xi
  .
\end{equation}
We are left with computing $\partial_\mu\Xi$:
\begin{equation}
  \begin{largearray}
    \frac{\partial\Xi_{lN+j}(\mathbb U)}{\partial\mu}=
    \partial_\mu\mathbb D^{(0)}(r_{l,j},e,\epsilon)
    +
    \sum_{l'}\sum_n\mathfrak F_{l',n}^{(\nu_u)}(\mathfrak u)\partial_\mu\mathbb D^{(1)}_{l',n}(r_{l,j},e,\epsilon)
    +\\\hfill+
    \sum_{l',l''}\sum_{n,m}\mathfrak F_{l',n}^{(\nu_u)}(\mathfrak u)\mathfrak F_{l'',m}^{(\nu_u)}(\mathfrak u)\partial_\mu\mathbb D^{(2)}_{l',n;l'',m}(r_{l,j},e,\epsilon)
  \end{largearray}
\end{equation}
\begin{equation}
  \begin{largearray}
    \frac{\partial\Xi_{NJ+1}(\mathbb U)}{\partial \mu}
    =\\\hfill=
    2\pi\rho
    \frac
    {
      \partial_\mu C(\epsilon)
      -
      \partial_\mu\mathfrak g^{(0)}(e,\epsilon)
      -
      \sum_l\sum_n\mathfrak F_{l,n}^{(\nu_u)}(\mathfrak u)\partial_\mu\mathfrak g^{(1)}_{l,n}(e,\epsilon)
      -
      \sum_{l,l'}\sum_{n,m}\mathfrak F_{l,n}^{(\nu_u)}(\mathfrak u)\mathfrak F_{l',m}^{(\nu_u)}(\mathfrak u)\partial_\mu\mathfrak g^{(2)}_{l,n;l',m}(e,\epsilon)
    }
    {1-\frac83\pi\rho+\rho^2\left(
      \bar{\mathfrak g}^{(0)}
      +
      \sum_l\sum_n\mathfrak F_{l,n}^{(\nu_u)}(\mathfrak u)\bar{\mathfrak g}^{(1)}_{l,n}
      +
      \sum_{l,l'}\sum_{n,m}\mathfrak F_{l,n}^{(\nu_u)}(\mathfrak u)\mathfrak F_{l',m}^{(\nu_u)}(\mathfrak u)\bar{\mathfrak g}^{(2)}_{l,n;l',m}
    \right)}
    .
  \end{largearray}
\end{equation}
We then use $\partial_\mu=\frac1{4\sqrt\epsilon}\partial_{\sqrt\epsilon}$ and
\begin{equation}
  \left.\frac{\partial C(\epsilon)}{\partial\sqrt\epsilon}\right|_{\mu=0}
  =
  \frac83\sqrt\epsilon+2
\end{equation}
\begin{equation}
  \begin{largearray}
    \left.\frac{\partial D^{(0)}(r)}{\partial\sqrt\epsilon}\right|_{\mu=0}=
    -\frac{2r}{r+1}e^{-2\sqrt{e}r}
    \\[0.5cm]\hfill
    +
    \frac{\rho}{(r+1)}\int_0^{\min(1,r)} ds\ (s+1)B^{(0)}(s)e^{-2\sqrt{e}r}((-1-2\sqrt er)\sinh(2\sqrt{e} s)+2\sqrt es\cosh(2\sqrt{e} s))
    \\[0.5cm]
    +
    \mathds 1_{r\leqslant 1}
    \frac{\rho}{2(r+1)}\int_0^{1-r} d\sigma\ (\sigma+r+1)B^{(0)}(\sigma+r)
    \cdot\\[0.5cm]\hfill\cdot
    \left((-1-2\sqrt e\sigma)\left(1-e^{-4\sqrt er}\right)e^{-2\sqrt{e}\sigma}+4\sqrt ere^{-2\sqrt{e}(2r+\sigma)}\right)
  \end{largearray}
\end{equation}
\begin{equation}
  \begin{largearray}
    \left.\frac{\partial D_{l,n}^{(1)}(r)}{\partial\sqrt\epsilon}\right|_{\mu=0}
    =
    \frac{\rho}{r+1}\int_0^{r} ds\ (s+1)B_{l,n}^{(1)}(s)e^{-2\sqrt er}\left(
      (-1-2\sqrt er)\sinh(2\sqrt e s)
      +2\sqrt es\cosh(2\sqrt e s)
    \right)
    \\[0.5cm]
    +
    \frac{\rho}{2(r+1)}\int_0^\infty d\sigma\ (\sigma+r+1)B_{l,n}^{(1)}(\sigma+r)
    \cdot\\[0.5cm]\hfill\cdot
    \left(
      (-1-2\sqrt e\sigma)\left(1-e^{-4\sqrt er}\right)e^{-2\sqrt e\sigma}
      +
      4\sqrt e re^{-2\sqrt e(2r+\sigma)}
    \right)
  \end{largearray}
\end{equation}
\begin{equation}
  \begin{largearray}
    \left.\frac{\partial D_{l,n;l',m}^{(2)}(r)}{\partial\sqrt\epsilon}\right|_{\mu=0}
    \\[0.5cm]
    =
    \frac{\rho}{r+1}\int_0^{r} ds\ (s+1)B_{l,n;l',m}^{(2)}(s)e^{-2\sqrt er}\left(
      (-1-2\sqrt er)\sinh(2\sqrt e s)
      +2\sqrt es\cosh(2\sqrt e s)
    \right)
    \\[0.5cm]
    +
    \frac{\rho}{2(r+1)}\int_0^\infty d\sigma\ (\sigma+r+1)B_{l,n;l',m}^{(2)}(\sigma+r)
    \cdot\\[0.5cm]\hfill\cdot
    \left(
      (-1-2\sqrt e\sigma)\left(1-e^{-4\sqrt er}\right)e^{-2\sqrt e\sigma}
      +
      4\sqrt e re^{-2\sqrt e(2r+\sigma)}
    \right)
    .
  \end{largearray}
\end{equation}
\begin{equation}
  \left.\frac{\partial \gamma^{(0)}}{\partial\sqrt\epsilon}\right|_{\mu=0}
  =
  -4\rho e\int_0^1 d\sigma\ \sigma(\sigma+1)B^{(0)}(\sigma)e^{-2\sqrt e\sigma}
\end{equation}
\begin{equation}
  \left.\frac{\partial\gamma^{(1)}_{l,n}}{\partial\sqrt\epsilon}\right|_{\mu=0}
  =
  -4\rho e\int_0^\infty d\sigma\ (\sigma+1)\sigma B_{l,n}^{(1)}(\sigma)e^{-2\sqrt e\sigma}
\end{equation}
and
\begin{equation}
  \left.\frac{\partial\gamma^{(2)}_{l,n;l',m}}{\partial\sqrt\epsilon}\right|_{\mu=0}
  =
  -4\rho e\int_0^\infty d\sigma\ (\sigma+1)\sigma B_{l,n;l',m}^{(2)}(\sigma)e^{-2\sqrt e\sigma}
  .
\end{equation}

\subsection{\tt simpleq-iteration}\label{sec:simpleq-iteration}
\indent
This method is used to solve the Simple equation using the iteration described in\-~\cite{CJL20}.
The Simple equation is
\begin{equation}
  -\Delta u
  =
  S-4e u+2e\rho u\ast u
  \label{simpleq_iteration}
\end{equation}
\begin{equation}
  S:=(1-u)v
  ,\quad
  \rho:=\frac{2e}{\int dx\ (1-u(|x|))v(|x|)}
  .
  \label{simpleq_Se}
\end{equation}
for a soft potential $v$ at fixed energy $e>0$.
The iteration is defined as
\begin{equation}
  -\Delta u_n=S_n-4eu_n+2e\rho_{n-1}u_{n-1}\ast u_{n-1}
  ,\quad
  u_0=0
  ,\quad
  \rho_{n-1}=
  \frac{2e}{\int dx\ (1-u_{n-1}(|x|))v(|x|)}
  .
  \label{iteration}
\end{equation}
\bigskip

\subsubsection{Usage}
\indent
Unless otherwise noted, this method takes the following parameters (specified via the {\tt [-p params]} flag, as a `{\tt;}' separated list of entries).
\begin{itemize}
  \item\makelink{param:simpleq-iteration_e}{}
  {\tt e} ({\tt Float64}, default: $10^{-4}$): energy $e$.

  \item\makelink{param:simpleq-iteration_maxiter}{}
  {\tt maxiter} ({\tt Int64}, default: 21): maximal number of iterations.

  \item\makelink{param:simpleq-iteration_order}{}
  {\tt order} ({\tt Int64}, default: 100): order used for all Gauss quadratures (denoted by $N$ below).

\end{itemize}
\bigskip

The available {\tt commands} are the following.

\begin{itemize}
  \item {\tt rho\_e}\makelink{command:simpleq-iteration_rho_e}{}:
  compute the density $\rho$ as a function of $e$.\par
  \underline{Disabled parameters}: \refname{param:simpleq-iteration_e}{{\tt e}}.\par
  \underline{Extra parameters}:
  \begin{itemize}
    \item\makelink{param:simpleq-iteration_minle}{}
    {\tt minle} ({\tt Float64}, default: $10^{-6}$): minimal value for $\log_{10}e$.

    \item\makelink{param:simpleq-iteration_maxle}{}
    {\tt maxle} ({\tt Float64}, default: $10^{2}$): maximal value for $\log_{10}e$.

    \item\makelink{param:simpleq-iteration_nle}{}
    {\tt nle} ({\tt Int64}, default: $100$): number of values for $e$ (spaced logarithmically).

    \item\makelink{param:simpleq-iteration_es}{}
    {\tt es} ({\tt Array\{Float64\}}, default: $(10^{{\tt minle}+\frac{{\tt maxle}-{\tt minle}}{{\tt nle}}n})_n$: list of values for $e$, specified as a `{\tt,}' separated list.
  \end{itemize}
  \underline{Output} (one line for each value of $e$): [$e$] [$\rho$].

  \item\makelink{command:simpleq-iteration_ux}{}
  {\tt ux}: compute $u$ as a function of $|x|$.\par
  \underline{Extra parameters}:
  \begin{itemize}
    \item\makelink{param:simpleq-iteration_xmin}{}
    {\tt xmin} ({\tt Float64}, default: 0): minimum of the range of $|x|$ to be printed.

    \item\makelink{param:simpleq-iteration_xmax}{}
    {\tt xmax} ({\tt Float64}, default: 100): maximum of the range of $|x|$ to be printed.

    \item\makelink{param:simpleq-iteration_nx}{}
    {\tt nx} ({\tt Int64}, default: 100): number of points to print (linearly spaced).
  \end{itemize}
  \underline{Output} (one line for each value of $x$): [$|x|$] [$u_1(|x|)$] [$u_2(|x|)$] [$u_3(|x|)$] $\cdots$

\end{itemize}

\subsubsection{Description}
\indent
In Fourier space
\begin{equation}
  \hat u_n(|k|):=\int dx\ e^{ikx}u_n(|x|)
\end{equation}
(\ref{iteration}) becomes
\begin{equation}
  (k^2+4e)\hat u_n(|k|)=\hat S_n(|k|)+2e\rho_{n-1}\hat u_{n-1}(k)^2
  ,\quad
  \rho_n:=\frac{2e}{\hat S_n(0)}
  \label{iteration_uk}
\end{equation}
with
\begin{equation}
  \hat S_n(|k|)
  =\hat v(|k|)-\frac1{(2\pi)^3}\int dp\ \hat u_n(|p|)\hat v(|k-p|)
  .
\end{equation}
We write $\hat S_n$ in bipolar coordinates (see lemma\-~\ref{lemma:bipolar}):
\begin{equation}
  \hat S_n(|k|)=\hat v(|k|)-\frac1{(2\pi)^3}\int_0^\infty dt\ t\hat u_n(t)\Eta(|k|,t)
\end{equation}
with
\begin{equation}
  \Eta(y,t):=\frac{2\pi}y\int_{|y-t|}^{y+t}ds\ s\hat v(s)
\end{equation}
(note that this agrees with the function\-~(\ref{easyeq_Eta}) defined for \refname{sec:easyeq}{{\tt easyeq}}).
We also change variables to
\begin{equation}
  y=\frac1{t+1}
  ,\quad
  t=\frac{1-y}y
\end{equation}
\begin{equation}
  \hat S_n(|k|)=\hat v(|k|)-\frac1{(2\pi)^3}\int_0^1 dy\ \frac{(1-y)\hat u_n(\frac{1-y}y)\Eta(|k|,\frac{1-y}y)}{y^3}
  .
\end{equation}
We approximate this integral using a Gauss-Legendre quadrature (see appendix\-~\ref{appGL}) and discretize Fourier space:
\begin{equation}
  k_i:=\frac{1-x_i}{1+x_i}
  ,\quad
  y_i:=\frac{x_i+1}2
\end{equation}
where $x_i$ are the Gauss-Legendre abscissa, and
\begin{equation}
  \hat S_n(k_i)\approx\hat v(k_i)-\frac1{2(2\pi)^3}\sum_{j=1}^Nw_j\frac{(1-y_j)\hat u_n(\frac{1-y_j}{y_j})\Eta(k_i,\frac{1-y_j}{y_j})}{y_j^3}
\end{equation}
so if
\begin{equation}
  \mathbb U_i(n):=\hat u_n(k_i)
  ,\quad
  \mathbf v_i:=\hat v(k_i)
\end{equation}
we have
\begin{equation}
  \sum_{j=1}^NA_{i,j}\mathbb U_j(n)=b_i^{(n)}
\end{equation}
with
\begin{equation}
  A_{i,j}:=(k_i^2+4e)\delta_{i,j}+\frac{w_j(1-y_j)\Eta(k_i,\frac{1-y_j}{y_j})}{2(2\pi)^3y_j^3}
\end{equation}
and
\begin{equation}
  b_i^{(n)}:=\mathbf v_i+2e\rho_{n-1}\mathbb U_i(n-1)^2
\end{equation}
in terms of which
\begin{equation}
  \mathbb U=A^{-1}b^{(n)}
  .
\end{equation}
Finally, we compute $\rho_n$ using the second of\-~(\ref{iteration_uk}):
\begin{equation}
  \rho_n=\frac{2e}{\hat S_n(0)}
  ,\quad
  S_n(0)\approx\hat v(0)
  -\frac1{2(2\pi)^3}\sum_{j=1}^Nw_j\frac{(1-y_j)\hat u_n(\frac{1-y_j}{y_j})\Eta(0,\frac{1-y_j}{y_j})}{y_j^3}
  .
\end{equation}

\section{Potentials}\label{sec:potentials}
\indent
In this section, we describe the potentials available in {\tt simplesolv}, and provide documentation to add custom potentials to {\tt simplesolv}.
\bigskip

\subsection{Built-in potentials}
\subsubsection{{\tt exp}}\label{sec:exp}
\indent
In $x$-space,
\begin{equation}
  v(|x|)=ae^{-|x|}
  .
\end{equation}
The constant $a$ is specified through the {\tt v\_a} parameter, and can be any real number.
Note that $v\geqslant 0$ if and only if $a\geqslant 0$.
This is the potential that is selected by default.
\bigskip

\point
In Fourier space,
\begin{equation}
  \hat v(|k|)=\int dx\ e^{ikx}v(|x|)
  =\frac{8\pi a}{(1+k^2)^2}
  .
\end{equation}
In particular, $v$ is of positive type (that is, $\hat v\geqslant 0$) if and only if $a\geqslant 0$.
\bigskip

\point
The zero energy scattering solution, that is, the solution of
\begin{equation}
  (-\Delta+v)\psi=0
  ,\quad
  \lim_{|x|\to\infty}\psi=1
\end{equation}
is
\begin{equation}
  \psi(r)=\frac{cI_0(2\sqrt a e^{-\frac r2})+2K_0(2\sqrt a e^{-\frac r2})}r
  ,\quad
  c=-\frac{2K_0(2\sqrt a)}{I_0(2\sqrt a)}
\end{equation}
where $I_0$ and $K_0$ are the modified Bessel functions, and the square root is taken with a branch cut on the negative imaginary axis.
In other words, if $a<0$, $\sqrt a=i\sqrt{|a|}$ and\-~\dlmfcite{(10.27.6),(10.27.9),(10.27.10)}{1.1.3}
\begin{equation}
  I_0(ix)=J_0(x)
  ,\quad
  K_0(ix)=-\frac\pi2(Y_0(x)+iJ_0(x))
\end{equation}
where $J_0$ and $Y_0$ are the Bessel functions, so
\begin{equation}
  \psi(r)=\pi\frac{c'J_0(2\sqrt{|a|} e^{-\frac r2})-Y_0(2\sqrt{|a|} e^{-\frac r2})}r
  ,\quad
  c'=\frac{Y_0(2\sqrt{|a|})}{J_0(2\sqrt{|a|})}
  .
\end{equation}
The scattering length is\-~\dlmfcite{(10.25.2),(10.31.2)}{1.1.3}
\begin{equation}
  \mathfrak a_0=\lim_{r\to\infty}r(1-\psi(r))
  =\log(a)+2\gamma+\frac{2K_0(2\sqrt a)}{I_0(2\sqrt a)}
\end{equation}
where $\gamma$ is the Euler constant, which, for $a<0$, is
\begin{equation}
  \mathfrak a_0=\log|a|+2\gamma-\frac{\pi Y_0(2\sqrt{|a|})}{J_0(2\sqrt{|a|})}
  .
\end{equation}

\subsubsection{{\tt tent}}\label{sec:tent}
\indent
In $x$-space,
\begin{equation}
  v(|x|)=\mathds 1_{|x|<b}a\frac{2\pi}3\left(1-\frac{|x|}b\right)^2\left(\frac{|x|}b+2\right)
  .
\end{equation}
The constants $a$ and $b$ are specified through the {\tt v\_a} and {\tt v\_b} parameters, and $a\in\mathbb R$, $b>0$.
Note that $v\geqslant 0$ if and only if $a\geqslant 0$.
Note that
\begin{equation}
  v(|x|)=\frac a{b^3}\mathds 1_{|x|<\frac b2}\ast\mathds 1_{|x|<\frac b2}
\end{equation}
(this is more easily checked in Fourier space than explicitly computing the convolution).
\bigskip

\indent
In Fourier space,
\begin{equation}
  \hat v(|k|)=\int dx\ e^{ikx}v(|x|)=a\frac{b^3}8\left(4\pi\frac{\sin(\frac{|k|b}2)-\frac{|k|b}2\cos(\frac{|k|b}2)}{\frac{(|k|b)^3}8}\right)^2
  .
\end{equation}
Note that this is $\frac a{b^3}(\hat\mathds 1_{|x|<\frac b2})^2$.
In particular, $v$ is of positive type (that is, $\hat v\geqslant 0$) if and only if $a\geqslant 0$.

\subsubsection{{\tt expcry}}\label{sec:expcry}
\indent
In $x$-space
\begin{equation}
  v(|x|)=e^{-|x|}-ae^{-b|x|}
  .
\end{equation}
The constants $a$ and $b$ are specified through the {\tt v\_a} and {\tt v\_b} parameters, and $a\in\mathbb R$, $b>0$.
Note that $v\geqslant 0$ if and only if $a\leqslant 1$ and $b\geqslant 1$.
\bigskip

\indent
In Fourier space
\begin{equation}
  \hat v(|k|)=\int dx\ e^{ikx}v(|x|)
  =8\pi\left(\frac1{(1+k^2)^2}-\frac{ab}{(b^2+k^2)^2}\right)
\end{equation}
In particular, $\hat v$ is of positive type (that is, $\hat v\geqslant 0$) if and only if $ab\leqslant 1$, $a\leqslant b$ and $a\leqslant b^3$.
If $a\leqslant 1$, $b\geqslant 1$ and $ab>1$, then $\hat v$ has a unique minimum at
\begin{equation}
  |k_*|=\sqrt{\frac{b^2-(ab)^{\frac13}}{(ab)^{\frac13}-1}}
  ,\quad
  \hat v(|k_*|)=
  -8\pi\frac{((ab)^{\frac13}-1)^3}{b^2-1}
  .
\end{equation}

\subsubsection{{\tt npt}}\label{sec:npt}
\indent
In $x$-space
\begin{equation}
  v(|x|)=x^2e^{-|x|}
  .
\end{equation}
Note that $v\geqslant 0$.
\bigskip

\indent
In Fourier space
\begin{equation}
  \hat v(|k|)=\int dx\ e^{ikx}v(|x|)
  =96\pi\frac{1-k^2}{(1+k^2)^4}
  .
\end{equation}
In particular, $v$ is not of positive type (that is, $\hat v$ is not $\geqslant 0$).

\subsubsection{{\tt alg}}\label{sec:alg}
\indent
In $x$-space
\begin{equation}
  v(|x|)=\frac1{1+\frac14|x|^4}
  .
\end{equation}
Note that $v\geqslant 0$.
\bigskip

\indent
In Fourier space
\begin{equation}
  \hat v(|k|)=\int dx\ e^{ikx}v(|x|)
  =4\pi^2 e^{-|k|}\frac{\sin(|k|)}{|k|}
  .
\end{equation}
In particular, $v$ is not of positive type (that is, $\hat v$ is not $\geqslant 0$).

\subsubsection{{\tt algwell}}\label{sec:algwell}
\indent
In $x$-space
\begin{equation}
  v(|x|)=\frac{1+a|x|^4}{(1+|x|^2)^4}
  .
\end{equation}
The constant $a$ can be set using the parameter {\tt v\_a} and can be any real number.
Note that $v\geqslant 0$ if and only if $a\geqslant 0$.
If $a>8$, this potential has a local minimum at $|x_-|$ and a local maximum at $|x_+|$:
\begin{equation}
  |x_\pm|=\sqrt{\frac12(1\pm\sqrt{1-8a^{-1}})}
  .
\end{equation}
\bigskip

\indent
In Fourier space
\begin{equation}
  \hat v(|k|)=\int dx\ e^{ikx}v(|x|)
  =\frac{\pi^2}{24}e^{-|k|}(a(k^2-9|k|+15)+k^2+3|k|+3)
  .
\end{equation}
In particular, $v$ is of positive type (that is, $\hat v$ is not $\geqslant 0$) if and only if
\begin{equation}
  -\frac15\leqslant a\leqslant 3+\frac87\sqrt7\approx6.02
  .
\end{equation}

\subsubsection{{\tt exact}}\label{sec:exact}
\indent
In $x$-space
\begin{equation}
  v(|x|)=
  \frac{12c( x^6b^6(2e-b^2) +b^4x^4(9e-7b^2) +4b^2x^2(3e-2b^2) +(5e+16b^2))}{(1+b^2x^2)^2(4+b^2x^2)^2((1+b^2x^2)^2-c)}
\end{equation}
The constants $a,b,c,e$ can be set using the parameters {\tt v\_a}, {\tt v\_b}, {\tt v\_c}, {\tt v\_e}, and $a,e\in\mathbb R$, $b\neq0$, $c>0$, $c\neq9$.
Note that $v\geqslant 0$ if and only if\-~\cite{CJL20b}
\begin{equation}
  b>0
  ,\quad
  0<c<1
  ,\quad
  \frac e{b^2}\geqslant\frac{-263+23\sqrt{161}}{48}\approx 0.601
  .
\end{equation}
With this potential, the Simple equation has an exact solution:
\begin{equation}
  u=\frac c{(1+b^2x^2)^2}
  ,\quad
  \rho=\frac{b^3}{c\pi^2}
  .
\end{equation}
\bigskip

\indent
In Fourier space $\hat v(|k|)=\int dx\ e^{ikx}v(x)$ is
\begin{equation}
  \begin{largearray}
    \hat v(|k|)=\frac{48 c\pi^2}{|k|}
    \left(
      \frac{18+3\sqrt c-(4-3\alpha)c-(1-2\alpha)c^{\frac32}}{4(3+\sqrt c)^2c^{\frac32}}
      e^{-\sqrt{1-\sqrt c}\frac{|k|}b}
      \right.\\[0.5cm]\hfill\left.
      +
      \frac{-18+3\sqrt c+(4-3\alpha)c-(1-2\alpha)c^{\frac32}}{4(3-\sqrt c)^2c^{\frac32}}
      e^{-\sqrt{1+\sqrt c}\frac{|k|}b}
      +\frac{1-\frac{|k|}b}{2c}
      e^{-\frac{|k|}b}
      -\frac{\alpha(3(9-c)\frac{|k|}b+8c)}{8(9-c)^2}
      e^{-2\frac{|k|}b}
    \right)
  \end{largearray}
\end{equation}
with $\alpha:=\frac e{b^2}$.

\subsection{Programming custom potentials}
\indent
In this section, we provide documentation for programming custom potentials.
\bigskip

\indent
The potentials are implemented in the file `{\tt \$SIMPLESOLV/src/potentials.jl}', and consist of two functions, one specifying the potential in Fourier space (the ``potential function''), and the other returning an approximate value for the scattering length (the ``scatterin glength function'') (as is explained below, a precise value of the scattering length is not actually needed).
For instance, the potential \refname{sec:exp}{{\tt exp}} has two functions: {\tt v\_exp} and {\tt a0\_exp}.
The potential function should take the following arguments:
\begin{itemize}
  \item {\tt k} ({\tt Float64}): the Fourier momentum
  \item and any parameters that the potential depends on, such as $a$ in \refname{sec:exp}{{\tt exp}} (can be of any type, provided the appropriate changes are made to {\tt main.jl} as explained below)
\end{itemize}
and it must return a {\tt Float64}: the value of $\hat v$ at {\tt k}.
The scattering length function takes the same parameters as an input, and returns a {\tt Float64}: the approximate value for the scattering length.
\bigskip

\indent
In addition, the potential must be linked in `{\tt \$SIMPLESOLV/src/main.jl}'.
In that file, the potential is read from the command line option {\tt U}.
The relevant code is in lines 197-222.
To add a new potential, add
\begin{code}
  elseif potential=="{\rm \{name of potential\}}"\par
  \indent  v=k->{\rm \{potential function\}}(k,v\_param\_a,v\_param\_b,{\rm...})\par
  \indent  a0={\rm \{scattering length function\}}(v\_param\_a,v\_param\_b,{\rm...})\par
\end{code}
where the number of {\tt v\_param} entries should be the number of parameters of the potential.
The parameters that are currently read from the parameters list are {\tt a}, {\tt b}, {\tt c} and {\tt e}.
To add a parameter, it must first be declared and initialized after line 35, and code to read it should be added after line 172:
\begin{code}
  elseif lhs="v\_{\rm\{name of parameter\}}"\par
  \indent  v\_param\_{\rm\{name of parameter\}}=parse(Float64,rhs)
\end{code}
If the new parameter has a type other than {\tt Float64}, this should be changed in the {\tt parse} function, and in the initialization.
\bigskip

\indent
The approximation of the scattering length is only used to initialize the Newton algorithm for \refname{sec:easyeq}{{\tt easyeq}}, so it is not important that it be exact.
In fact, some of the built-in potentials set the scattering length to 1, when it has proved too difficult to compute it exactly.

\appendix

\section{Chebyshev polynomial expansion}\label{appChebyshev}
\indent
In this appendix, we compute the error of Chebyshev polynomial expansions of regular functions.
Specifically, we will consider class-$s$ Gevrey functions (which is a generalization of the notion of analyticity: analytic functions are Gevrey functions with $s=1$).
A class-$s$ Gevrey function on $[-1,1]$ is a $\mathcal C^\infty$ function that satisfies, $\forall n\in\mathbb N$,
\begin{equation}
  \mathop{\mathrm{sup}}_{x\in[-1,1]}\left|\frac{d^nf(x)}{dx^n}\right|\leqslant C_0C^n(n!)^s.
\end{equation}
Formally, the Chebyshev polynomial expansion of $f$ is
\begin{equation}
  f(x)=\frac{c_0}2+\sum_{j=1}^\infty c_j T_j(x)
  \label{eqcheby}
\end{equation}
where $T_j$ is the $j$-th Chebyshev polynomial:
\begin{equation}
  T_j(x):=\cos(j\arccos(x))
\end{equation}
and
\begin{equation}
  c_j:=\frac2\pi\int_0^\pi d\theta\ f(\cos\theta)\cos(j\theta)
  .
\end{equation}

\bigskip

\theo{Lemma}\label{lemmaChebyshev}
  Let $f$ be a class-$s$ Gevrey function on $[-1,1]$ with $s\in\mathbb N\setminus\{0\}$.
  There exist $b_0,b>0$, which are independent of $s$, such that the coefficients $c_j$ of the Chebyshev polynomial expansion are bounded as
  \begin{equation}
    c_j\leqslant b_0e^{-bj^{\frac1s}}
    .
    \label{boundchebyshevcoef}
  \end{equation}
  In particular, the Chebyshev polynomial expansion is absolutely convergent pointwise (and, therefore in any $L_p$ norm), and, for every $N\geqslant 1$,
  \begin{equation}
    \left|f(x)-\frac{c_0}2-\sum_{j=1}^N c_j T_j(x)\right|
    \leqslant
    \frac{b_0}{1-e^{-b}}(s-1)!(N^{\frac1s}+b^{-1})^{s-1}e^{-bN^{\frac1s}}
    .
    \label{boundchebyshev}
  \end{equation}
\endtheo
\bigskip

\indent\underline{Proof}:
  Note that~\-(\ref{eqcheby}) is nothing other than the Fourier cosine series expansion of $F(\theta):=f(\cos(\theta))$, which is an even, periodic, class-$s$ Gevrey function on $[-\pi,\pi]$, whose $j$-th Fourier coefficient for $j\in\mathbb Z$ is equal to $\frac12c_{|j|}$.
  The bound\-~(\ref{boundchebyshevcoef}) follows from a well-known estimate of the decay of Fourier coefficients of class-$s$ Gevrey functions (see e.g.~\-\cite[Theorem~\-3.3]{Ta87}).
  The bound\-~(\ref{boundchebyshev}) then follows from $|T_j(x)|\leqslant 1$ and lemma\-~\ref{lemma:ineqgevrey} below.
\qed
\bigskip

\theo{Lemma}\label{lemma:ineqgevrey}
  Given $b>0$ and two integers $N,s>0$,
  \nopagebreakaftereq
  \begin{equation}
    \sum_{j=N}^\infty e^{-bj^{\frac1s}}
    \leqslant
    (s-1)!\left(N^{\frac1s}+\frac1{1-e^{-b}}\right)^{s-1}\frac{e^{-bN^{\frac1s}}}{1-e^{-b}}
  \end{equation}
\endtheo
\restorepagebreakaftereq
\bigskip

\indent\underline{Proof}:
  If $\nu_{N,s}^s:=\lfloor N^{\frac1s}\rfloor^s$ denotes the largest integer that is $\leqslant N$ and has an integer $s$-th root, then
  \begin{equation}
    \sum_{j=N}^\infty e^{-bj^{\frac1s}}\leqslant
    \sum_{j=\nu_{N,s}^s}^\infty e^{-bj^{\frac1s}}\leqslant
    \sum_{k=\nu_{N,s}}^\infty(k^s-(k-1)^s)e^{-bk}\leqslant
    s\sum_{k=\nu_{N,s}}^\infty k^{s-1}e^{-bk}
    .
  \end{equation}
  We then estimate
  \begin{equation}
    \sum_{k=\nu_{N,s}}^\infty k^{s-1}e^{-bk}=\frac{d^{s-1}}{d(-b)^{s-1}}\sum_{k=\nu_{N,s}}^\infty e^{-bk}
    \leqslant (s-1)!\left(\nu_{N,s}+\frac1{1-e^{-b}}\right)^{s-1}\frac{e^{-b\nu_{N,s}}}{1-e^{-b}}
  \end{equation}
  which concludes the proof of the lemma.
\qed

\section{Gauss quadratures}\label{appGL}
\indent
Gauss quadratures are approximation schemes to compute integrals of the form
\begin{equation}
  \int_a^bdt\ \omega(t)f(t)
\end{equation}
where $\omega(x)\geqslant 0$ is one of several functions that Gauss quadratures can treat.
The possible choices of $\omega$ and $(a,b)$ are
\begin{itemize}
  \item $(a,b)=(-1,1)$, $\omega(t)=1$: Gauss-Legendre quadratures.
  \item $(a,b)=(-1,1)$, $\omega(t)=(1-t)^\alpha(1+t)^\beta$, $\alpha,\beta>-1$: Gauss-Jacobi quadratures.
  \item $(a,b)=(-1,1)$, $\omega(t)=(1-t^2)^{-\frac12}$: Gauss-Chebyshev quadratures.
  \item $(a,b)=(0,\infty)$, $\omega(t)=e^{-t}$: Gauss-Laguerre quadratures.
  \item $(a,b)=(0,\infty)$, $\omega(t)=t^\alpha e^{-t}$, $\alpha>-1$: generalized Gauss-Laguerre quadratures.
  \item $(a,b)=(0,\infty)$, $\omega(t)=e^{-t^2}$: Gauss-Hermite quadratures.
\end{itemize}
It is not our goal here to discuss Gauss quadratures in detail, or their relation to orthogonal polynomials.
Instead, we will compute the error made when approximating such an integral by a Gauss quadrature.
\bigskip

\indent
For each Gauss quadrature, the integral is approximated in the form
\begin{equation}
  \int_a^b dt\ \omega(t) f(t)
  \approx
  \sum_{i=1}^Nw_if(r_i)
\end{equation}
where $w_i$ are called the {\it weights}, $r_i$ are the {\it abcissa}, and $N$ is the {\it order}.
The weights and abcissa depend on both $\omega$ and the order $N$.
The crucial property of Gauss quadratures is that they are {\it exact} when $f$ is a polynomial of order $\leqslant 2N-1$.
\bigskip

\indent
In this appendix, we compute the error of Gauss quadratures when used to integrate regular functions.
Specifically, we will consider class-$s$ Gevrey functions (which is a generalization of the notion of analyticity: analytic functions are Gevrey functions with $s=1$).
A class-$s$ Gevrey function on is a $\mathcal C^\infty$ function that satisfies, $\forall n\in\mathbb N$,
\begin{equation}
  \mathop{\mathrm{sup}}_{x\in[-1,1]}\left|\frac{d^nf(x)}{dx^n}\right|\leqslant C_0C^n(n!)^s.
\end{equation}

\bigskip

\theo{Lemma}\label{lemmaGL}
  Let $f$ be a class-$s$ Gevrey function with $s\in\mathbb N\setminus\{0\}$.
  There exist $b_0,b>0$, which are independent of $s$, and $N_0>0$, which is independent of $s$ and $f$, such that, if $N\geqslant N_0$, then, denoting the Gauss weights and abcissa by $w_i$ and $r_i$,
  \begin{equation}
    \left|\int_a^bdt\ \omega(t)f(t)-\sum_{i=1}^Nw_if(r_i)\right|
    \leqslant
    b_0(s-1)!((2N-1)^{\frac1s}+b^{-1})^{s-1}e^{-b(2N-1)^{\frac1s}}\left(\int_a^b\omega(t)+\sum_{i=1}^Nw_i\right)
    .
  \end{equation}
  In particular, if $f$ is analytic (i.e. $s=1$), then
  \begin{equation}
    \left|\int_a^bdt\ \omega(t)f(t)-\sum_{i=1}^Nw_if(r_i)\right|
    \leqslant b_0e^{-b(2N-1)}\left(\int_a^b\omega(t)+\sum_{i=1}^Nw_i\right)
    .
  \end{equation}
\endtheo
\bigskip

\indent\underline{Proof}:
  We expand $f$ into Chebyshev polynomials as in\-~(\ref{eqcheby}):
  \begin{equation}
    f(x)=\frac{c_0}2+\sum_{j=1}^{\infty}c_jT_j(x)
  \end{equation}
  Let
  \begin{equation}
    g(x):=f(x)-\frac{c_0}2-\sum_{j=1}^{2N-1} c_j T_j(x)
    .
  \end{equation}
  Since order-$N$ Gauss quadratures are exact on polynomials of order $\leqslant 2N-1$, we have
  \begin{equation}
    \int_a^bdt\ \omega(t)f(t)-\sum_{i=1}^Nw_if(r_i)
    =
    \int_a^bdt\ \omega(t)g(t)-\sum_{i=1}^Nw_ig(r_i)
  \end{equation}
  and, by lemma\-~\ref{lemmaChebyshev},
  \begin{equation}
    |g(x)|\leqslant
    (\mathrm{const}.)(s-1)!((2N-1)^{\frac1s}+b^{-1})^{s-1}e^{-b(2N-1)^{\frac1s}}
    .
    \label{boundchebyshev}
  \end{equation}
\qed

\section{Bipolar coordinates}\label{appendix:bipolar}
\indent
Bipolar coordinates are very useful for computing convolutions of radial functions in three dimensions.
\bigskip

\theo{Lemma}\label{lemma:bipolar}
  For $y\in\mathbb R^3$,
  \nopagebreakaftereq
  \begin{equation}
    \int_{\mathbb R^3}dx\ f(|x|,|x-y|)
    =
    \frac{2\pi}{|y|}\int_0^\infty dt\int_{||y|-t|}^{|y|+t}ds\ stf(s,t)
  \end{equation}
\endtheo
\restorepagebreakaftereq

\indent\underline{Proof}:
  Without loss of generality, we assume that $y=(0,0,a)$ with $a>0$. We first change to cylindrical coordinates: $(\rho,\theta,x_3)$:
  \begin{equation}
    \int_{\mathbb R^3}dx\ f(|x|,|x-y|)
    =
    2\pi\int_0^\infty d\rho\int_{-\infty}^\infty dx_3\ \rho f(|(\rho,0,x_3)|,|(\rho,0,x_3-a)|)
    .
  \end{equation}
  Next, we change variables to
  \begin{equation}
    s=|(\rho,0,x_3)|
    ,\quad
    t=|(\rho,0,x_3-a)|
    .
  \end{equation}
  The inverse of this transformation is
  \begin{equation}
    x_3=\pm\frac{a^2+s^2-t^2}{2a}
    ,\quad
    \rho=\frac1{2a}\sqrt{4a^2s^2-(a^2+s^2-t^2)^2}
    \label{inverse}
  \end{equation}
  and its Jacobian is
  \begin{equation}
    \frac{2ts}{\sqrt{-2a^4+2a^2(t^2+s^2)-(t^2-s^2)^2}}
    =\frac{ts}{\rho a}
    .
    \label{jacobian}
  \end{equation}
\qed
\bigskip

\indent
The following is a generalization of the previous lemma to functions of four variables.
\bigskip

\theo{Lemma}\label{lemma:bipolar2}
  For $y\in\mathbb R^3$,
  \begin{equation}
    \begin{largearray}
      \int_{\mathbb R^6}dxdx'\ f(|x|,|x-y|,|x'|,|x'-y|,|x-x'|)
      \\[0.5cm]\hfill
      =
      \frac{2\pi}{|y|^2}\int_0^\infty dt\int_{||y|-t|}^{|y|+t}ds\int_0^\infty dt'\int_{||y|-t'|}^{|y|+t'}ds'\int_0^{2\pi}d\theta\ sts't'f(s,t,s',t',\xi(s,t,s',t',\theta,|y|))
    \end{largearray}
  \end{equation}
  with
  \nopagebreakaftereq
  \begin{equation}
    \begin{largearray}
      \xi(s,t,s',t',\theta,|y|)
      :=
      \frac1{\sqrt 2|y|}
      \Big(
	|y|^2(s^2+t^2+(s')^2+(t')^2)-|y|^4-(s^2-t^2)((s')^2-(t')^2)
	\\[0.5cm]\hfill
	-
	\sqrt{(4|y|^2s^2-(|y|^2+s^2-t^2)^2)(4|y|^2(s')^2-(|y|^2+(s')^2-(t')^2)^2)}\cos\theta
      \Big)^{\frac12}
      .
    \end{largearray}
  \end{equation}
\endtheo
\restorepagebreakaftereq

\indent\underline{Proof}:
  Without loss of generality, we assume that $y=(0,0,a)$ with $a>0$. We first change to cylindrical coordinates: $(\rho,\theta,x_3;\rho',\theta',x_3')$:
  \begin{equation}
    \begin{largearray}
      \int_{\mathbb R^6}dxdx'\ f(|x|,|x-y|,|x'|,|x'-y|,|x-x'|)
      \\[0.5cm]\hfill
      =
      2\pi\int_0^\infty d\rho\int_{-\infty}^\infty dx_3
      \int_0^\infty d\rho'\int_{-\infty}^\infty dx_3'\int_0^{2\pi}d\theta'
      \ \rho\rho'
      f(
	s,t,s',t',
	|(\rho-\rho'\cos\theta',\rho'\sin\theta',x_3-x_3')|
      )
      .
    \end{largearray}
  \end{equation}
  where
  \begin{equation}
    s:=|(\rho,0,x_3)|
    ,\quad
    t:=|(\rho,0,x_3-a)|
    ,\quad
    s':=|(\rho'\cos\theta',\rho'\sin\theta',x_3')|,
    t':=|(\rho'\cos\theta',\rho'\sin\theta',x_3'-a)|
  \end{equation}
  Next, we change variables to $(s,t,s',t',\theta')$.
  The Jacobian of this transformation is, by\-~(\ref{jacobian}),
  \begin{equation}
    \frac{tst's'}{\rho\rho' a^2}
    .
  \end{equation}
  Furthermore, by\-~(\ref{inverse}),
  \begin{equation}
    x_3-x_3'=\frac{s^2-t^2-(s')^2+(t')^2}{2a}
  \end{equation}
  and
  \begin{equation}
    \rho=\frac{\sqrt{4a^2s^2-(a^2+s^2-t^2)^2}}{2a}
    ,\quad
    \rho'=\frac{\sqrt{4a^2(s')^2-(a^2+(s')^2-(t')^2)^2}}{2a}
  \end{equation}
  so
  \begin{equation}
    \begin{largearray}
      |(\rho-\rho'\cos\theta',\rho'\sin\theta',x_3-x_3')|
      =
      \frac1{2a}
      \left(
	4a^2s^2-(a^2+s^2-t^2)^2
	+
	4a^2(s')^2-(a^2+(s')^2-(t')^2)^2
	\right.\\[0.5cm]\hfill\left.
	-
	2\sqrt{(4a^2s^2-(a^2+s^2-t^2)^2)(4a^2(s')^2-(a^2+(s')^2-(t')^2)^2)}\cos\theta'
	+
	(s^2-t^2-(s')^2+(t')^2)^2
      \right)^{\frac12}
    \end{largearray}
  \end{equation}
  and
  \begin{equation}
    \begin{largearray}
      |(\rho-\rho'\cos\theta',\rho'\sin\theta',x_3-x_3')|
      =
      \frac1{\sqrt 2a}
      \Big(
	a^2(s^2+t^2+(s')^2+(t')^2)-a^4-(s^2-t^2)((s')^2-(t')^2)
	\\[0.5cm]\hfill
	-
	\sqrt{(4a^2s^2-(a^2+s^2-t^2)^2)(4a^2(s')^2-(a^2+(s')^2-(t')^2)^2)}\cos\theta'
      \Big)^{\frac12}
      .
    \end{largearray}
  \end{equation}
\qed

\section{Hann windowing}\label{appendix:hann}
\indent
In this appendix, we discuss the use of Hann windows to compute Fourier transforms.
Consider the Fourier transform
\begin{equation}
  \hat f(k)=\int dx\ e^{ikx}f(x)
  .
\end{equation}
Evaluating this integral numerically can be tricky, especially at high $|k|$, because of the rapid oscillations at large $|x|$.
A trick to palliate such a problem is to multiply $f$ by a {\it window function} $h_L$, which cuts off distances of order $L$.
We then compute, instead of $\hat f$,
\begin{equation}
  \tilde f(k)=\int dx\ e^{ikx}h_L(x)f(x)
  .
\end{equation}
We can then evaluate $\tilde f$ using standard numerical techniques, such as Gauss quadratures (see appendix\-~\ref{appGL}), without issues at large $|x|$.
However, in doing so, we will make an error in the Fourier transform.
To quantify this error, note that
\begin{equation}
  \tilde f(k)
  =\hat h_L\hat\ast\hat f(k)
  \equiv\int\frac{dq}{(2\pi)^d}\ \hat h_L(q)\hat f(k-q)
\end{equation}
so if we choose $h_L$ in such a way that $\hat h_L$ is peaked around the origin, then $\tilde f$ will not differ too much from $\hat f$:
\begin{equation}
  \hat f(k)-\tilde f(k)
  =((2\pi)^d\delta(k)-\hat h_L)\hat\ast\hat f(k)
  .
\end{equation}
\bigskip

\indent
The {\it Hann} window is defined as
\begin{equation}
  h_L(x)=\cos^2({\textstyle\frac{\pi|x|}L})\mathds 1_{|x|<\frac L2}
\end{equation}
whose Fourier transform is, in $d=3$,
\begin{equation}
  \hat h_L(k)=
  \frac{4\pi^3 L^2}{|k|}\frac{((|k|L)^3-4|k|L\pi^2)\cos(\frac{|k|L}2)-2(3(|k|L)^2-4\pi^2)\sin(\frac{|k|L}2)}{((|k|L)^3-4|k|L\pi^2)^2}
\end{equation}
which decays at large $|k|L$ as
\begin{equation}
  \hat h_L(k)\sim\frac{4\pi^3}{|k|^4L}\cos({\textstyle\frac{|k|L}2})
  .
\end{equation}



\vfill
\eject

\begin{thebibliography}{WWW99}
\small
\IfFileExists{bibliography/bibliography.tex}{\input bibliography/bibliography.tex}{}
\end{thebibliography}


\end{document}