{\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,CJL21,CHe21}.
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.
{\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|)=\mathds1_{|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).
{\tt minlrho\_init} ({\tt Float64}, default: $-6$): 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.
{\tt nlrho\_init} ({\tt Int64}, default: 0): number of steps in the initialization process described above. Set to 0 to disable the incremental initialization process.
{\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.
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
{\tt minrho} ({\tt Float64}, default: $10^{-6}$): minimal value for $\rho$.
\item\makelink{param:easyeq_maxrho}{}
{\tt maxrho} ({\tt Float64}, default: $10^{2}$): maximal value for $\rho$.
\item\makelink{param:easyeq_nrho}{}
{\tt nrho} ({\tt Int64}, default: $0$): number of values for $\rho$ (spaced linearly). If {\tt nrho} is $\neq0$, then the linear spacing will be used, and \refname{param:easyeq_minlrho}{{\tt minlrho}}, \refname{param:easyeq_maxlrho}{{\tt maxlrho}}, \refname{param:easyeq_nlrho}{{\tt nlrho}} will be ignored. Otherwise, the logarithmic spacing will be used and \refname{param:easyeq_minrho}{{\tt minrho}}, \refname{param:easyeq_maxrho}{{\tt maxrho}} will be ignored.
{\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|)$]
{\tt 2pt}: compute the spherically averaged two-point correlation function $C_2(|x|)$ at a given $\rho$.\par
\underline{Extra parameters}: same as \refname{command:easyeq_ux}{{\tt ux}}, plus\par
\begin{itemize}
\item\makelink{param:easyeq_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 two-point correlation function, see\-~(\ref{hann}).
\end{itemize}
\underline{Output} (one line for each value of $|x|$): [$|x|$] [$C_2(|x|)$]
\item\makelink{command:easyeq_2pt_max}{}
{\tt 2pt\_max}: compute the maximum of the spherically averaged two-point correlation function $C_2(|x|)$ at a given $\rho$.\par
\underline{Extra parameters}: \refname{param:easyeq_window_L}{{\tt window\_L}} plus
\begin{itemize}
\item\makelink{param:easyeq_dx}{}
{\tt dx} ({\tt Float64}, default: $10^{-7}$): step used to numerically approximate derivatives.
\item\makelink{param:easyeq_x0}{}
{\tt x0} ({\tt Float64}, default: $1$): initial guess for the maximum is $\rho^{-1/3}\mathrm{\tt x0}$.
\item\makelink{param:easyeq_maxstep}{}
{\tt maxstep} ({\tt Float64}, default: $\infty$): maximal size of single step in maximization algorithm.
\item\makelink{param:easyeq_tolerance_max}{}
{\tt tolerance\_max} ({\tt Float64}, default: \refname{param:easyeq_tolerance}{{\tt tolerance}}): same as \refname{param:easyeq_tolerance}{{\tt tolerance}}, used for the Newton algorithm underlying the maximization algorithm.
{\tt 2pt\_max\_rho}: compute the maximum of the spherically averaged two-point correlation function $C_2(|x|)$ for a range of $\rho$.\par
\underline{Extra parameters}: same as \refname{command:easyeq_2pt_max}{{\tt easyeq\_2pt\_max}} plus those of \refname{command:easyeq_energy_rho}{{\tt energy\_rho}}.\par
\underline{Output} (one line for each value of $\rho$): [$\rho$] [$|x_{\mathrm{max}}|$] [$C_2(|x_{\mathrm{max}}|)$]\par
\underline{Multithread support}: yes, different values of $\rho$ split up among workers.
\item\makelink{command:easyeq_2pt_fourier}{}
{\tt 2pt\_fourier}: compute the spherically averaged Fourier transform of the two-point correlation function $\hat C_2(|k|)$ at a given $\rho$.\par
\underline{Extra parameters}:
\begin{itemize}
\item\makelink{param:easyeq_kmin}{}
{\tt kmin} ({\tt Float64}, default: 0): minimum of the range of $|k|$ to be printed.
\item\makelink{param:easyeq_kmax}{}
{\tt kmax} ({\tt Float64}, default: 10): maximum of the range of $|k|$ to be printed.
\item\makelink{param:easyeq_nk}{}
{\tt nk} ({\tt Int64}, default: 100): number of $|k|$'s to be printed.
{\tt window\_L} ({\tt Float64}, default: 1000): what is actually computed is $\hat C_2$ convolved with a Gaussian of variance $1/\sqrt L$ with $L=\mathrm{\tt window\_L}$, see\-~(\ref{easyeq_gaussian}).
\end{itemize}
\underline{Output} (one line for each value of $|k|$): [$|k|$] [$\hat C_2(|k|)$].\par
\underline{Multithread support}: yes, different values of $k$ are split up among workers.
\item\makelink{command:easyeq_2pt_fourier_max}{}
{\tt 2pt\_fourier\_max}: compute the maximum of the spherically averaged Fourier transformed two-point correlation function $\hat C_2(|k|)$.\par
\underline{Extra parameters}: \refname{param:easyeq_2pt_fourier_window_L}{{\tt window\_L}}, \refname{param:easyeq_maxstep}{{\tt maxstep}} plus
\begin{itemize}
\item\makelink{param:easyeq_dk}{}
{\tt dk} ({\tt Float64}, default: $10^{-7}$): step used to numerically approximate derivatives.
\item\makelink{param:easyeq_k0}{}
{\tt k0} ({\tt Float64}, default: $1$): initial guess for the maximum is $\rho^{1/3}\mathrm{\tt k0}$.
{\tt 2pt\_fourier\_max\_rho}: compute the maximum of the spherically averaged Fourier transformed two-point correlation function $\hat C_2(|k|)$ for a range of $\rho$.\par
\underline{Extra parameters}: same as those of \refname{command:easyeq_2pt_fourier_max}{{\tt 2pt\_fourier\_max}} plus those of \refname{command:easyeq_energy_rho}{{\tt energy\_rho}}.\par
\underline{Output} (one line for each value of $\rho$): [$\rho$] [$|k_{\mathrm{max}}|$] [$\hat C_2(|k_{\mathrm{max}}|)$]\par
\underline{Multithread support}: yes, different values of $\rho$ split up among workers.
{\tt momentum\_distribution}: compute the momentum distribution $\mathcal M(|k|)$ at a given $\rho$. The momentum distribution is computed for $k=k_{l,j}$ (see\-~(\ref{klj})).\par
\underline{Extra parameters}:
\begin{itemize}
\item\makelink{param:easyeq_momentum_kmin}{}
{\tt kmin} ({\tt Float64}, default: 0): minimum of the range of $|k|$ to be printed.
\item\makelink{param:easyeq_momentum_kmax}{}
{\tt kmax} ({\tt Float64}, default: 10): maximum of the range of $|k|$ to be printed.
\item\makelink{param:easyeq_momentum_window_L}{}
{\tt window\_L} ({\tt Float64}, default: 1000): what is actually computed is $\mathfrak M_k$ convolved with a Gaussian of variance $1/\sqrt L$ where $L=\sqrt{\mathrm{\tt window\_L}}/k^2$, see\-~(\ref{easyeq_gaussian_momt}).
\end{itemize}
\underline{Output} (one line for each value of $|k|$): [$|k|$] [$\mathcal M(|k|)$]
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}}).
We iterate the Newton algorithm until the Newton relative error $\epsilon$ becomes smaller than the \refname{param:easyeq_tolerance}{{\tt tolerance}} parameter.
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$:
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:
This is a modification of the usual Newton iteration $x_n+\partial C_2/\partial^2C_2$ which is designed to follow the direction of the gradient, and thus to move toward a local maximum.
In addition, if $|\partial C_2|/|\partial^2 C_2|$ is larger than \refname{param:easyeq_maxstep}{{\tt maxstep}}, then the step is replaced with $\pm$\refname{param:easyeq_maxstep}{{\tt maxstep}}.
This prevents the algorithm from stepping over a maximum and land on another, further away.
This is useful if one has a good idea of where the global maximum is, and does not want to get trapped in a smaller local maximum.
\bigskip
\indent
The algorithm is run for a maximum of \refname{param:easyeq_maxiter}{{\tt maxiter}} iterations, or until $|x_{n+1}-x_n|$ is smaller than \refname{param:easyeq_tolerance}{{\tt tolerance}}.
If the maximal number of iterations is reached, or if the solution found is not a local maximum, then the algorithm fails, and returns $+\infty$.
The point thus computed is therefore a local maximum, but it is not guaranteed to be the global maximum.
\subsubsubsection{Fourier transform of two-point correlation (spherical average)}
The Fourier transform of the two-point correlation function is
To compute $\frac{\delta e}{\delta\hat v(q)}$, one idea would be to proceed in the same way as for the two-point correlation function, by replacing $\hat v$ with
\begin{equation}
\hat v+\lambda g(|k|)
,\quad
g(|k|):=\frac1{4\pi|q|^2}\delta(|q|-|k|)
\end{equation}
where $\delta$ is the Dirac-delta function distribution (compare this with\-~(\ref{2pt_addv})).
However, the $\delta$ function causes all sorts of problems with the quadratures.
\bigskip
\point
Instead, we approximate $\hat C_2$ by convolving it with a normalized Gaussian: let
The computation of $\partial_\lambda\Xi|_{\lambda=0}$ is identical to\-~(\ref{dXi_2pt_easyeq}), but with the $g$ defined in\-~(\ref{easyeq_2pt_fourier_g}).
\bigskip
\point
To compute the maximum of $\hat C_2$, we proceed as for $C_2$, see\-~(\ref{easyeq_newton_2pt})-(\ref{easyeq_dx_2pt}).
The only difference is that the algorithm is initialized with $|k_0|=\rho^{\frac13}$\refname{param:easyeq_k0}{{\tt k0}}.
\subsubsubsection{Momentum distribution}
To compute the momentum distribution (see\-~\cite{CHe21}), we add a parameter $\lambda$ to {\tt easyeq}:
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.
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).
{\tt nlrho\_init} ({\tt Int64}, default: 0): 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}}. If {\tt nlrho\_init} is $\neq0$, then the solution of \refname{eq:easyeq_medeq}{{\tt medeq}} is first computed for {\tt nlrho} smaller values of $\rho$ starting from $10^{\mathrm{\refname{param:anyeq_minlrho_init}{{\tt minlrho\_init}}}}$. This is useful when $\rho$ is too large for the solution of \refname{eq:easyeq_medeq}{{\tt medeq}} to be computed directly.
\item\makelink{param:anyeq_minlrho_init}{}
{\tt minlrho\_init} ({\tt Float64}, default: $-6$): see \refname{param:anyeq_nlrho_init}{{\tt nlrho\_init}}.
{\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
{\tt minrho} ({\tt Float64}, default: $10^{-6}$): minimal value for $\rho$.
\item\makelink{param:anyeq_maxrho}{}
{\tt maxrho} ({\tt Float64}, default: $10^{2}$): maximal value for $\rho$.
\item\makelink{param:anyeq_nrho}{}
{\tt nrho} ({\tt Int64}, default: $0$): number of values for $\rho$ (spaced linearly). If {\tt nrho} is $\neq0$, then the linear spacing will be used, and \refname{param:anyeq_minlrho}{{\tt minlrho}}, \refname{param:anyeq_maxlrho}{{\tt maxlrho}}, \refname{param:anyeq_nlrho}{{\tt nlrho}} will be ignored. Otherwise, the logarithmic spacing will be used and \refname{param:anyeq_minrho}{{\tt minrho}}, \refname{param:anyeq_maxrho}{{\tt maxrho}} will be ignored.
{\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$].
{\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$].
{\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|$] [$u(|x|)$]
\item\makelink{command:anyeq_2pt}{}
{\tt 2pt}: compute the spherically averaged 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_2pt_max}{}
{\tt 2pt\_max}: compute the maximum of the spherically averaged two-point correlation function $C_2(|x|)$ at a given $\rho$.\par
\underline{Extra parameters}: \refname{param:anyeq_window_L}{{\tt window\_L}} plus
\begin{itemize}
\item\makelink{param:anyeq_dx}{}
{\tt dx} ({\tt Float64}, default: $10^{-7}$): step used to numerically approximate derivatives.
\item\makelink{param:anyeq_x0}{}
{\tt x0} ({\tt Float64}, default: $1$): initial guess for the maximum is $\rho^{-1/3}\mathrm{\tt x0}$.
\item\makelink{param:anyeq_maxstep}{}
{\tt maxstep} ({\tt Float64}, default: $\infty$): maximal size of single step in maximization algorithm.
\item\makelink{param:easyeq_tolerance_max}{}
{\tt tolerance\_max} ({\tt Float64}, default: \refname{param:easyeq_tolerance}{{\tt tolerance}}): same as \refname{param:easyeq_tolerance}{{\tt tolerance}}, used for the Newton algorithm underlying the maximization algorithm.
{\tt 2pt\_max\_rho}: compute the maximum of the spherically averaged two-point correlation function $C_2(|x|)$ for a range of $\rho$.\par
\underline{Extra parameters}: same as \refname{command:anyeq_2pt_max}{{\tt anyeq\_2pt\_max}} plus those of \refname{command:anyeq_energy_rho}{{\tt energy\_rho}}.\par
\underline{Output} (one line for each value of $\rho$): [$\rho$] [$|x_{\mathrm{max}}|$] [$C_2(|x_{\mathrm{max}}|)$]\par
\underline{Multithread support}: yes, different values of $\rho$ split up among workers.
\item\makelink{command:anyeq_2pt_fourier}{}
{\tt 2pt\_fourier}: compute the spherically averaged Fourier transform of the two-point correlation function $\hat C_2(|k|)$ at a given $\rho$.\par
{\tt window\_L} ({\tt Float64}, default: 1000): what is actually computed is $\hat C_2$ convolved with a Gaussian of variance $1/\sqrt L$ where $L=\sqrt{\mathrm{\tt window\_L}}$, see\-~(\ref{anyeq_gaussian}).
{\tt 2pt\_fourier\_max\_rho}: compute the maximum of the spherically averaged Fourier transformed two-point correlation function $\hat C_2(|k|)$ for a range of $\rho$.\par
\underline{Extra parameters}: same as those of \refname{command:anyeq_2pt_fourier_max}{{\tt 2pt\_fourier\_max}} plus those of \refname{command:anyeq_energy_rho}{{\tt energy\_rho}}.\par
\underline{Output} (one line for each value of $\rho$): [$\rho$] [$|k_{\mathrm{max}}|$] [$\hat C_2(|k_{\mathrm{max}}|)$]\par
\underline{Multithread support}: yes, different values of $\rho$ split up among workers.
{\tt momentum\_distribution}: compute the momentum distribution $\mathcal M(|k|)$ at a given $\rho$. The momentum distribution is computed for $k=k_{l,j}$ (see\-~(\ref{klj})).\par
{\tt kmin} ({\tt Float64}, default: 0): minimum of the range of $|k|$ to be printed.
\item\makelink{param:anyeq_momentum_kmax}{}
{\tt kmax} ({\tt Float64}, default: 10): maximum of the range of $|k|$ to be printed.
\item\makelink{param:anyeq_momentum_window_L}{}
{\tt window\_L} ({\tt Float64}, default: 1000): what is actually computed is $\mathfrak M_k$ convolved with a Gaussian of variance $1/\sqrt L$ where $L=\sqrt{\mathrm{\tt window\_L}}/k^2$, see\-~(\ref{anyeq_gaussian_momt}).
{\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{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.
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}$.
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$.
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}$,
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}:
\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}$.
(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\geqslant2$).
Note that these integrals are not singular as long as $\nu_a,\nu_b\geqslant2$: 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
\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)
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}}.
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.
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:
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}} are passed on to \refname{sec:easyeq}{{\tt 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}}.
We iterate the Newton algorithm until the Newton relative error $\epsilon$ becomes smaller than the \refname{param:easyeq_tolerance}{{\tt tolerance}} parameter.
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}):
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
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:
This is a modification of the usual Newton iteration $x_n+\partial C_2/\partial^2C_2$ which is designed to follow the direction of the gradient, and thus to move toward a local maximum.
In addition, if $|\partial C_2|/|\partial^2 C_2|$ is larger than \refname{param:anyeq_maxstep}{{\tt maxstep}}, then the step is replaced with $\pm$\refname{param:anyeq_maxstep}{{\tt maxstep}}.
This prevents the algorithm from stepping over a maximum and land on another, further away.
This is useful if one has a good idea of where the global maximum is, and does not want to get trapped in a smaller local maximum.
\bigskip
\indent
The algorithm is run for a maximum of \refname{param:anyeq_maxiter}{{\tt maxiter}} iterations, or until $|x_{n+1}-x_n|$ is smaller than \refname{param:anyeq_tolerance}{{\tt tolerance}}.
If the maximal number of iterations is reached, or if the solution found is not a local maximum, then the algorithm fails, and returns $+\infty$.
The point thus computed is therefore a local maximum, but it is not guaranteed to be the global maximum.
\subsubsubsection{Fourier transform of two-point correlation (spherical average)}
The Fourier transform of the two-point correlation function is
To compute $\frac{\delta e}{\delta\hat v(q)}$, one idea would be to proceed in the same way as for the two-point correlation function, by replacing $\hat v$ with
\begin{equation}
\hat v+\lambda g(|k|)
,\quad
g(|k|):=\frac1{4\pi|q|^2}\delta(|q|-|k|)
\end{equation}
where $\delta$ is the Dirac-delta function distribution (compare this with\-~(\ref{2pt_addv})).
However, the $\delta$ function causes all sorts of problems with the Chebyshev polynomial exansion and the quadratures.
\bigskip
\point
Instead, we approximate $\hat C_2$ by convolving it with a normalized Gaussian: let
should formally be added to the right side of\-~(\ref{Omega}).
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 Y_{l,j}$ should be replaced by
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}}.
{\tt minrho} ({\tt Float64}, default: $10^{-6}$): minimal value for $\rho$.
\item\makelink{param:simpleq-Kv_maxrho}{}
{\tt maxrho} ({\tt Float64}, default: $10^{2}$): maximal value for $\rho$.
\item\makelink{param:simpleq-Kv_nrho}{}
{\tt nrho} ({\tt Int64}, default: $0$): number of values for $\rho$ (spaced linearly). If {\tt nrho} is $\neq0$, then the linear spacing will be used, and \refname{param:simpleq-Kv_minlrho}{{\tt minlrho}}, \refname{param:simpleq-Kv_maxlrho}{{\tt maxlrho}}, \refname{param:simpleq-Kv_nlrho}{{\tt nlrho}} will be ignored. Otherwise, the logarithmic spacing will be used and \refname{param:simpleq-Kv_minrho}{{\tt minrho}}, \refname{param:simpleq-Kv_maxrho}{{\tt maxrho}} will be ignored.
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:
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$.
{\tt minrho} ({\tt Float64}, default: $10^{-6}$): minimal value for $\rho$.
\item\makelink{param:simpleq-hardcore_maxrho}{}
{\tt maxrho} ({\tt Float64}, default: $10^{2}$): maximal value for $\rho$.
\item\makelink{param:simpleq-hardcore_nrho}{}
{\tt nrho} ({\tt Int64}, default: $0$): number of values for $\rho$ (spaced linearly). If {\tt nrho} is $\neq0$, then the linear spacing will be used, and \refname{param:simpleq-hardcore_minlrho}{{\tt minlrho}}, \refname{param:simpleq-hardcore_maxlrho}{{\tt maxlrho}}, \refname{param:simpleq-hardcore_nlrho}{{\tt nlrho}} will be ignored. Otherwise, the logarithmic spacing will be used and \refname{param:simpleq-hardcore_minrho}{{\tt minrho}}, \refname{param:simpleq-hardcore_maxrho}{{\tt maxrho}} will be ignored.
{\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$:
To compute the energy $e$ of this equation, with the extra parameter $\mu$, we consider the limit of the soft sphere potential $\lambda\mathds1_{|x|<1}$ (see\-~(\ref{easyeq}) with $\beta_K=\beta_L=0$):
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|\leqslant1$,
\subsubsubsection{Chebyshev polynomial expansion} We use the same interpolation as we used in \refname{sec:anyeq}{{\tt anyeq}}: (\ref{a_chebyshev_spline})
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
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)$)
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}}.
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.
We iterate the Newton algorithm until the Newton relative error $\epsilon$ becomes smaller than the \refname{param:easyeq_tolerance}{{\tt tolerance}} parameter.
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).
{\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$
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$.
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
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.
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$,
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)|\leqslant1$ and lemma\-~\ref{lemma:ineqgevrey} below.
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 $\leqslant2N-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$,
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$,
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.
