Update to v0.4v0.4

  feature: compute the 2-point correlation function in easyeq.

  feature: compute the Fourier transform of the 2-point correlation function
           in anyeq and easyeq.

  feature: compute the local maximum of the 2-point correlation function and
           its Fourier transform.

  feature: compute the compressibility for anyeq.

  feature: allow for linear spacing of rho's.

  feature: print the scattering length.

  change: ux and uk now return real numbers.

  fix: error in the computation of the momentum distribution: wrong
       definition of delta functions.

  fix: various minor bugs.

  optimization: assign explicit types to variables.

9 files changed, 1064 insertions, 128 deletions

diff --git a/doc/Makefile b/doc/Makefile
index 4b58b3c..d0c4c64 100644
--- a/doc/Makefile
+++ b/doc/Makefile
@@ -1,4 +1,4 @@
-## Copyright 2021 Ian Jauslin
+## Copyright 2021-2023 Ian Jauslin
 ## 
 ## Licensed under the Apache License, Version 2.0 (the "License");
 ## you may not use this file except in compliance with the License.
diff --git a/doc/bibliography/bibliography.tex b/doc/bibliography/bibliography.tex
index 53454da..e8893ca 100644
--- a/doc/bibliography/bibliography.tex
+++ b/doc/bibliography/bibliography.tex
@@ -4,10 +4,13 @@ doi:{\tt\color{blue}\href{http://dx.doi.org/10.1103/PhysRevA.103.053309}{10.1103
 \bibitem[CJL20]{CJL20}E.A. Carlen, I. Jauslin, E.H. Lieb - {\it Analysis of a simple equation for the ground state energy of the Bose gas}, Pure and Applied Analysis, volume~\-2, issue~\-3, pages~\-659-684, 2020,\par\penalty10000
 doi:{\tt\color{blue}\href{http://dx.doi.org/10.2140/paa.2020.2.659}{10.2140/paa.2020.2.659}}, arxiv:{\tt\color{blue}\href{http://arxiv.org/abs/1912.04987}{1912.04987}}.\par\medskip
  
-\bibitem[CJL20b]{CJL20b}E.A. Carlen, I. Jauslin, E.H. Lieb - {\it Analysis of a simple equation for the ground state of the Bose gas II: Monotonicity, Convexity and Condensate Fraction}, 2020, to appear in the SIAM journal of Mathematical Analysis,\par\penalty10000
-arxiv:{\tt\color{blue}\href{http://arxiv.org/abs/2010.13882}{2010.13882}}.\par\medskip
+\bibitem[CJL21]{CJL21}E.A. Carlen, I. Jauslin, E.H. Lieb - {\it Analysis of a Simple Equation for the Ground State of the Bose Gas II: Monotonicity, Convexity, and Condensate Fraction}, SIAM Journal on Mathematical Analysis, volume~\-53, number~\-5, pages~\-5322-5360, 2021,\par\penalty10000
+doi:{\tt\color{blue}\href{http://dx.doi.org/10.1137/20M1376820}{10.1137/20M1376820}}, arxiv:{\tt\color{blue}\href{http://arxiv.org/abs/2010.13882}{2010.13882}}.\par\medskip
  
 \bibitem[DLMF]{DLMF1.1.3}F.W.J. Olver, A.B. Olde Daalhuis, D.W. Lozier, B.I. Schneider, R.F. Boisvert, C.W. Clark, B.R. Miller, B.V. Saunders, H.S. Cohl, M.A. McClain (editors) - {\it NIST Digital Library of Mathematical Functions}, Release~\-1.1.3 of~\-2021-09-15, 2021.\par\medskip
  
+\bibitem[Ja23]{Ja23}I. Jauslin - {\it The Simplified approach to the Bose gas without translation invariance}, 2023\par\penalty10000
+{\tt\color{blue}\href{http://ian.jauslin.org/publications/23j/}{http://ian.jauslin.org/publications/23j/}}.\par\medskip
+ 
 \bibitem[Ta87]{Ta87}Y. Taguchi - {\it Fourier coefficients of periodic functions of Gevrey classes and ultradistributions}, Yokohama Mathematical Journal, volume~\-35, pages~\-51-60, 1987.\par\medskip
  
diff --git a/doc/libs/code.sty b/doc/libs/code.sty
index d4f4070..e067dfe 100644
--- a/doc/libs/code.sty
+++ b/doc/libs/code.sty
@@ -1,4 +1,4 @@
-%% Copyright 2021 Ian Jauslin
+%% Copyright 2021-2023 Ian Jauslin
 %% 
 %% Licensed under the Apache License, Version 2.0 (the "License");
 %% you may not use this file except in compliance with the License.
diff --git a/doc/libs/dlmf.sty b/doc/libs/dlmf.sty
index f435e96..f66157b 100644
--- a/doc/libs/dlmf.sty
+++ b/doc/libs/dlmf.sty
@@ -1,4 +1,4 @@
-%% Copyright 2021 Ian Jauslin
+%% Copyright 2021-2023 Ian Jauslin
 %% 
 %% Licensed under the Apache License, Version 2.0 (the "License");
 %% you may not use this file except in compliance with the License.
diff --git a/doc/libs/ian.cls b/doc/libs/ian.cls
index f29e6bd..762f679 100644
--- a/doc/libs/ian.cls
+++ b/doc/libs/ian.cls
@@ -219,6 +219,7 @@
   % counters
   \stepcounter{subsectioncount}
   \setcounter{subsubsectioncount}{0}
+  \setcounter{subsubsubsectioncount}{0}
   \ifsubsectionsineq
     \setcounter{seqcount}0
   \fi
@@ -268,6 +269,7 @@
 \def\subsubsection#1{
   % counters
   \stepcounter{subsubsectioncount}
+  \setcounter{subsubsubsectioncount}{0}
 
   % space before subsubsection (if not first)
   \ifnum\thesubsubsectioncount>1
@@ -302,6 +304,49 @@
   \medskip\penalty10000
 }
 
+%% subsubsubsection
+% counter
+\newcounter{subsubsubsectioncount}
+% space before subsubsubsection
+\newlength\subsubsubsecskip
+\setlength\subsubsubsecskip{20pt}
+\def\subsubsubsection#1{
+  % counters
+  \stepcounter{subsubsubsectioncount}
+
+  % space before subsubsubsection (if not first)
+  \ifnum\thesubsubsubsectioncount>1
+    \vskip\subsubsubsecskip
+    \penalty-500
+  \fi
+
+  % hyperref anchor
+  \hrefanchor
+  % define tag (for \label)
+  \xdef\tag{\sectionprefix\thesectioncount.\thesubsectioncount.\thesubsubsectioncount.\thesubsubsubsectioncount}
+
+  % get widths
+  \def\@secnum{{\bf\hskip1.5cm\sectionprefix\thesectioncount.\thesubsectioncount.\thesubsubsectioncount.\thesubsubsubsectioncount.\hskip5pt}}
+  \settowidth\secnumwidth{\@secnum}
+  \setlength\sectitlewidth\textwidth
+  \addtolength\sectitlewidth{-\secnumwidth}
+  % print name
+  \parbox{\textwidth}{
+  \@secnum
+  \parbox[t]{\sectitlewidth}{\bf #1}}
+
+  % write to table of contents
+  \iftoc
+    % save lncount in aux variable which is written to toc
+    \immediate\write\tocoutput{\noexpand\expandafter\noexpand\edef\noexpand\csname toc@subsubsubsec.\thesectioncount.\thesubsectioncount.\thesubsubsectioncount.\thesubsubsubsectioncount\endcsname{\thelncount}}
+    \write\tocoutput{\noexpand\tocsubsubsubsection{#1}{\thepage}}
+  \fi
+
+  % space
+  \par\penalty10000
+  \medskip\penalty10000
+}
+
 %% itemize
 \newlength\itemizeskip
 % left margin for items
@@ -638,6 +683,7 @@
   \stepcounter{tocsectioncount}
   \setcounter{tocsubsectioncount}{0}
   \setcounter{tocsubsubsectioncount}{0}
+  \setcounter{tocsubsubsubsectioncount}{0}
   % write
   \smallskip\hyperlink{ln.\csname toc@sec.\thetocsectioncount\endcsname}{{\bf \tocsectionprefix\thetocsectioncount}.\hskip5pt {\color{blue}#1}\leaderfill#2}\par
 }
@@ -645,6 +691,7 @@
 \def\tocsubsection #1#2{
   \stepcounter{tocsubsectioncount}
   \setcounter{tocsubsubsectioncount}{0}
+  \setcounter{tocsubsubsubsectioncount}{0}
   % write
   {\hskip10pt\hyperlink{ln.\csname toc@subsec.\thetocsectioncount.\thetocsubsectioncount\endcsname}{{\bf \thetocsectioncount.\thetocsubsectioncount}.\hskip5pt {\color{blue}\small #1}\leaderfill#2}}\par
 }
@@ -654,6 +701,12 @@
   % write
   {\hskip20pt\hyperlink{ln.\csname toc@subsubsec.\thetocsectioncount.\thetocsubsectioncount.\thetocsubsubsectioncount\endcsname}{{\bf \thetocsectioncount.\thetocsubsectioncount.\thetocsubsubsectioncount}.\hskip5pt {\color{blue}\small #1}\leaderfill#2}}\par
 }
+\newcounter{tocsubsubsubsectioncount}
+\def\tocsubsubsubsection #1#2{
+  \stepcounter{tocsubsubsubsectioncount}
+  % write
+  {\hskip30pt\hyperlink{ln.\csname toc@subsubsubsec.\thetocsectioncount.\thetocsubsectioncount.\thetocsubsubsectioncount.\thetocsubsubsubsectioncount\endcsname}{{\bf \thetocsectioncount.\thetocsubsectioncount.\thetocsubsubsectioncount.\thetocsubsubsubsectioncount}.\hskip5pt {\color{blue}\small #1}\leaderfill#2}}\par
+}
 \def\tocappendices{
   \medskip
   \setcounter{tocsectioncount}0
diff --git a/doc/libs/iantheo.sty b/doc/libs/iantheo.sty
index 1945a5f..31842a4 100644
--- a/doc/libs/iantheo.sty
+++ b/doc/libs/iantheo.sty
@@ -1,4 +1,4 @@
-%% Copyright 2021 Ian Jauslin
+%% Copyright 2021-2023 Ian Jauslin
 %% 
 %% Licensed under the Apache License, Version 2.0 (the "License");
 %% you may not use this file except in compliance with the License.
diff --git a/doc/libs/largearray.sty b/doc/libs/largearray.sty
index cf9075f..3a8762e 100644
--- a/doc/libs/largearray.sty
+++ b/doc/libs/largearray.sty
@@ -1,4 +1,4 @@
-%% Copyright 2021 Ian Jauslin
+%% Copyright 2021-2023 Ian Jauslin
 %% 
 %% Licensed under the Apache License, Version 2.0 (the "License");
 %% you may not use this file except in compliance with the License.
diff --git a/doc/libs/point.sty b/doc/libs/point.sty
index a396d1c..796e0d2 100644
--- a/doc/libs/point.sty
+++ b/doc/libs/point.sty
@@ -1,4 +1,4 @@
-%% Copyright 2021 Ian Jauslin
+%% Copyright 2021-2023 Ian Jauslin
 %% 
 %% Licensed under the Apache License, Version 2.0 (the "License");
 %% you may not use this file except in compliance with the License.
@@ -33,12 +33,15 @@
 \newif\ifresetatsubsubsection
 \DeclareOption{reset_at_subsubsection}{\resetatsubsubsectiontrue}
 \DeclareOption{no_reset_at_subsubsection}{\resetatsubsubsectionfalse}
+\newif\ifresetatsubsubsubsection
+\DeclareOption{reset_at_subsubsubsection}{\resetatsubsubsubsectiontrue}
+\DeclareOption{no_reset_at_subsubsubsection}{\resetatsubsubsubsectionfalse}
 \newif\ifresetattheo
 \DeclareOption{reset_at_theo}{\resetattheotrue}
 \DeclareOption{no_reset_at_theo}{\resetattheofalse}
 
 \def\point@defaultoptions{
-  \ExecuteOptions{reset_at_section, reset_at_subsection, reset_at_subsubsection, no_reset_at_theo}
+  \ExecuteOptions{reset_at_section, reset_at_subsection, reset_at_subsubsection, reset_at_subsubsubsection, no_reset_at_theo}
   \ProcessOptions
 
   %% reset at every new section
@@ -56,6 +59,11 @@
     \let\point@oldsubsubsection\subsubsection
     \gdef\subsubsection{\resetpointcounter\point@oldsubsubsection}
   \fi
+  %% reset at every new subsubsubsection
+  \ifresetatsubsubsubsection
+    \let\point@oldsubsubsubsection\subsubsubsection
+    \gdef\subsubsubsection{\resetpointcounter\point@oldsubsubsubsection}
+  \fi
 
   %% reset at every new theorem
   \ifresetattheo
diff --git a/doc/simplesolv-doc.tex b/doc/simplesolv-doc.tex
index 722282a..23c7247 100644
--- a/doc/simplesolv-doc.tex
+++ b/doc/simplesolv-doc.tex
@@ -1,4 +1,4 @@
-%% Copyright 2021 Ian Jauslin
+%% Copyright 2021-2023 Ian Jauslin
 %% 
 %% Licensed under the Apache License, Version 2.0 (the "License");
 %% you may not use this file except in compliance with the License.
@@ -27,7 +27,7 @@
 \hfil{\bf\Large
 {\tt simplesolv}\par
 \medskip
-\hfil \large v0.3
+\hfil \large v0.4
 }
 
 \vfill
@@ -41,7 +41,7 @@
 \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}.
+{\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}.
 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|)
@@ -85,6 +85,10 @@ The available methods are (see section\-~\ref{sec:methods} for further details)
   \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.
+In addition, the scattering length can be displayed using the {\tt scattering\_length} command:
+\begin{code}
+  julia \$SIMPLESOLV/src/main.jl scattering\_length
+\end{code}
 \bigskip
 
 \indent
@@ -165,10 +169,10 @@ Unless otherwise noted, this method takes the following parameters (specified vi
   {\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.
+  {\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.
 
   \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.
+  {\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.
 
   \item\makelink{param:easyeq_bK}{}
   {\tt bK}, {\tt bL} ({\tt Float64}, default: 1, 1): the values of $\beta_K$ and $\beta_L$.
@@ -188,7 +192,7 @@ The available {\tt commands} are the following.
   \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{Disabled parameters}: \refname{param:easyeq_rho}{{\tt rho}}.\par
   \underline{Extra parameters}:
   \begin{itemize}
     \item\makelink{param:easyeq_minlrho}{}
@@ -200,8 +204,18 @@ The available {\tt commands} are the following.
     \item\makelink{param:easyeq_nlrho}{}
     {\tt nlrho} ({\tt Int64}, default: $100$): number of values for $\rho$ (spaced logarithmically).
 
+    \item\makelink{param:easyeq_minrho}{}
+    {\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.
+
     \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.
+    {\tt rhos} ({\tt Array\{Float64\}}): list of values for $\rho$, specified as a `{\tt,}' separated list.
+    This parameter takes precedence over \refname{param:easyeq_minlrho}{{\tt minlrho}}, \refname{param:easyeq_maxlrho}{{\tt maxlrho}}, \refname{param:easyeq_nlrho}{{\tt nlrho}}, \refname{param:easyeq_minrho}{{\tt minrho}}, \refname{param:easyeq_maxrho}{{\tt maxrho}}, \refname{param:easyeq_nrho}{{\tt nrho}}.
   \end{itemize}
   \underline{Output} (one line for each value of $\rho$): [$\rho$] [$e$] [Newton error $\epsilon$].
 
@@ -241,15 +255,100 @@ The available {\tt commands} are the following.
   \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|)$]
 
+  \item\makelink{command:easyeq_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: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.
+  \end{itemize}
+  \underline{Output}: [$|x_{\mathrm{max}}|$] [$C_2(|x_{\mathrm{max}}|)$]
+
+  \item\makelink{command:easyeq_2pt_max_rho}{}
+  {\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.
+
+    \item\makelink{param:easyeq_2pt_fourier_window_L}{}
+    {\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}$.
+  \end{itemize}
+  \underline{Output}: [$|k_{\mathrm{max}}|$] [$\hat C_2(|k_{\mathrm{max}}|)$]\par
+
+  \item\makelink{command:easyeq_2pt_fourier_max_rho}{}
+  {\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.
+
+  \item\makelink{command:easyeq_momentum_distribution}{}
+  {\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|)$]
+
 \end{itemize}
 
 \subsubsection{Description}
-\point{\bf Fourier space formulation.}
+\subsubsubsection{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)
+  =\frac{4\pi}{|k|}\int_0^\infty dr\ r\sin(|k|r)f(r)
   .
 \end{equation}
 In Fourier space, (\ref{easyeq}) becomes
@@ -268,6 +367,7 @@ with
   ,\quad
   \hat f\hat\ast\hat g(|k|):=\int_{\mathbb R^3}\frac{dp}{(2\pi)^3}\ \hat f(|k-p|)\hat g(|p|)
   .
+  \label{hatS_easyeq}
 \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}
@@ -287,6 +387,15 @@ with
   \Phi(x):=\frac{2(1-\sqrt{1-x})}x
   .
 \end{equation}
+We can write this as a root finding problem:
+\begin{equation}
+  \Omega(u):=
+  \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)
+  =0
+  .
+  \label{Omega_easyeq}
+\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)
@@ -305,7 +414,7 @@ By a simple change of variables,
 \end{equation}
 \bigskip
 
-\point{\bf Evaluating integrals.}
+\subsubsubsection{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)
@@ -353,6 +462,7 @@ with
   \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}
+  \label{bbSE_easyeq}
 \end{equation}
 \begin{equation}
   \mathbb X_i:=\frac{k_i^2}{2\rho\mathbb A_{K,i}}
@@ -361,12 +471,12 @@ with
 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
+\subsubsubsection{Main algorithm to compute $\mathbb U$}\par\penalty10000
 \medskip\penalty10000
-\subpoint
+\point
 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)
+  \Xi_i(\mathbb U):=\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)
   =0
   \label{root_easyeq}
 \end{equation}
@@ -381,7 +491,7 @@ where $D\Xi$ is the Jacobian of $\Xi$:
 \end{equation}
 \bigskip
 
-\subpoint
+\point
 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}
@@ -397,7 +507,7 @@ 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
+\point
 We are left with computing the Jacobian of $\Xi$:
 \begin{equation}
   \begin{largearray}
@@ -436,7 +546,7 @@ with
 \end{equation}
 \bigskip
 
-\subpoint
+\point
 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}
@@ -466,8 +576,8 @@ and the solution $u$ in real space is obtained by inverting the Fourier transfor
 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})
+\subsubsubsection{Condensate fraction}
+Finally, to compute the uncondensed fraction, we solve the modified {\tt easyeq} (see\-~\cite{CJL21})
 \begin{equation}
   (-\Delta+2\mu)u_\mu=v(1-u_\mu)-2\rho K+\rho^2L
   \label{easyeq_mu}
@@ -518,6 +628,351 @@ We then compute, using\-~(\ref{jacobian_easyeq}),
   .
 \end{equation}
 
+\subsubsubsection{Correlation function (spherical average)}
+The two-point correlation function is
+\begin{equation}
+  c_2(x):=
+  2\rho\frac{\delta e}{\delta v(x)}
+\end{equation}
+and its spherical average is
+\begin{equation}
+  C_2(|x|):=\frac1{4\pi|x|^2}\int dy\ \delta(|x|-|y|)c_2(y)
+  .
+\end{equation}
+In Fourier space,
+\begin{equation}
+  c_2(x)=
+  2\rho
+  \int dk\ e^{ikx}\frac{\delta e}{\delta\hat v(k)}
+\end{equation}
+so
+\begin{equation}
+  C_2(|x|)=
+  2\rho\int dk\ \left(\frac1{4\pi|x|^2}\int dy\ \delta(|x|-|y|)e^{iky}\right)\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{easyeq_C2fourier}
+\end{equation}
+\bigskip
+
+\point
+We can compute this quantity by considering a modified {\tt easyeq} 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|}
+  .
+  \label{2pt_addv}
+\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}{\delta\hat v(k)}\partial_\lambda({\textstyle \hat v(k)+\lambda g(|k|)})
+  =\int dk\ g(|k|)\frac{\delta e}{\delta\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^2g(0)-\rho^2\int\frac{dk}{(2\pi)^3}\ (g(k)\hat u(k)+\hat v(k)\partial_\lambda\hat u_\lambda(k)|_{\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
+
+\point
+We compute $\partial_\lambda\Xi|_{\lambda=0}$:
+\begin{equation}
+  \begin{largearray}
+    \partial_\lambda\Xi_i
+    =
+    -\frac1{2(\mathbb X_i+1)}\left(
+      \partial_\lambda\mathbb T_i-\frac{\mathbb T_i\partial_\lambda\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_\lambda\mathbb T_i+\mathbb T_i\partial_\lambda\mathbb B_i-2\frac{\mathbb B_i\mathbb T_i\partial_\lambda\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{dXi_2pt_easyeq}
+\end{equation}
+with
+\begin{equation}
+  \partial_\lambda\mathbb B_i=
+  (\beta_L(1-\beta_K)-\beta_K(1-\beta_L))
+  \frac{\mathbb E\partial_\lambda\mathbb S_i-\mathbb S_i\partial_\lambda\mathbb E}{(\beta_K\mathbb S_i+(1-\beta_K)\mathbb E)^2}
+\end{equation}
+\begin{equation}
+  \partial_\lambda\mathbb T_i=
+  (1-\beta_K)\frac{\mathbb E\partial_\lambda\mathbb S_i-\mathbb S_i\partial_\lambda\mathbb E}{(\beta_K\mathbb S_i+(1-\beta_K)\mathbb E)^2}
+  ,\quad
+  \partial_\lambda\mathbb A_{K,i}=\beta_K\partial_\lambda\mathbb S_i+(1-\beta_K)\partial_\lambda\mathbb E
+\end{equation}
+\begin{equation}
+  \partial_\lambda\mathbb S_i:=g(k_i)
+  -\frac1{16\pi^3\rho}\sum_{j=1}^N w_j\frac{(1-y_j)\mathbb U_j\partial_\lambda\Eta(k_i,k_j)}{y_j^3}
+\end{equation}
+\begin{equation}
+  \partial_\lambda \mathbb E:=g(0)-\frac1{16\pi^3\rho}\sum_{j=1}^N w_j\frac{(1-y_j)\mathbb U_j\partial_\lambda\Eta(0,k_j)}{y_j^3}
+\end{equation}
+\begin{equation}
+  \partial_\lambda\mathbb X_i:=-\frac{k_i^2}{2\rho\mathbb A_{K,i}^2}\partial_\lambda\mathbb A_{K,i}
+\end{equation}
+where $\partial\lambda\Eta$ is computed similarly to $\Eta$\-~(\ref{Eta}) but with $\hat v$ replaced by $g$:
+\begin{equation}
+  \partial_\lambda\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))g((y+t)s+|y-t|(1-s))
+  .
+  \label{GG}
+\end{equation}
+\bigskip
+
+\point
+In order to invert the Fourier transform in\-~(\ref{easyeq_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:easyeq_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
+To compute the maximum of $C_2$, we use a modified Newton algorithm.
+The initial guess for the maximum is $|x_0|=\rho^{-\frac13}$\refname{param:easyeq_x0}{{\tt x0}}.
+The modified Newton algorithm is an iteration:
+\begin{equation}
+  x_{n+1}=x_n+\frac{\partial C_2(|x_n|)}{|\partial^2C_2(|x_n|)|}
+  \label{easyeq_newton_2pt}
+\end{equation}
+in which the derivatives are approximated using finite differences:
+\begin{equation}
+  \partial C_2(x)\approx \frac{C_2(|x|+dx)-C_2(|x|)}{dx}
+  ,\quad
+  \partial^2 C_2(x)\approx \frac{C_2(|x|+dx)+C_2(|x|-dx)-2C_2(|x|)}{dx^2}
+  .
+  \label{easyeq_dx_2pt}
+\end{equation}
+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
+\begin{equation}
+  \hat c_2(q):=
+  2\rho\frac{\delta e}{\delta v(q)}
+\end{equation}
+and its spherical average is
+\begin{equation}
+  \hat C_2(|q|):=\frac1{4\pi|q|^2}\int dk\ \delta(|q|-|k|)c_2(k)
+  =
+  \frac\rho{2\pi|q|^2}\int dk\ \delta(|q|-|k|)\frac{\delta e}{\delta\hat v(k)}
+  .
+\end{equation}
+\bigskip
+
+\point
+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
+\begin{equation}
+  \Gamma_L(|k|):=\left(\frac{L}{2\pi}\right)^{\frac32}e^{-\frac L2k^2}
+  \label{easyeq_gaussian}
+\end{equation}
+\begin{equation}
+  \hat{\mathfrak C}_2(|q|)
+  :=
+  \int dp\ \hat C_2(|q-p|)\Gamma_L(|p|)
+  =\int dk\ \int dp\ \frac\rho{2\pi|q-p|^2}\delta(|q-p|-|k|)\frac{\delta e}{\delta\hat v(k)}\Gamma_L(|p|)
+\end{equation}
+which by lemma\-~\ref{lemma:bipolar} is
+\begin{equation}
+  \hat{\mathfrak C}_2(|q|)
+  =\int dk\ 
+  \frac\rho{|q|}
+  \frac{\delta e}{\delta\hat v(k)}
+  \int_0^\infty dt\ \int_{||q|-t|}^{|q|+t}ds\ s\frac{\delta(t-|k|)}t\Gamma_L(s)
+\end{equation}
+that is
+\begin{equation}
+  \hat{\mathfrak C}_2(|q|)
+  =\int dk\ 
+  \frac{\delta e}{\delta\hat v(k)}
+  \frac{\rho}{|q||k|}
+  \int_{||q|-|k||}^{|q|+|k|}ds\ s\Gamma_L(s)
+\end{equation}
+which is the directional derivative of $e$ with respect to $\hat v$ in the direction of $2\rho g$ with
+\begin{equation}
+  g(|k|):=\frac1{2|q||k|}\int_{||q|-|k||}^{|q|+|k|}ds\ s\Gamma_L(s)
+  =
+  \frac1{2|k|rL}(\Gamma_L(|k|-r)-\Gamma_L(|k|+r))
+  .
+  \label{easyeq_2pt_fourier_g}
+\end{equation}
+Note that
+\begin{equation}
+  g(0):=\Gamma_L(|q|)
+  .
+\end{equation}
+To compute this derivative, we replace $\hat v$ with
+\begin{equation}
+  \hat v+\lambda g(|k|)
+\end{equation}
+so, denoting the solution of the modified equation by $u_\lambda$, for $q\neq 0$,
+\begin{equation}
+  \hat{\mathfrak C}_2(|q|)=2\rho\partial_\lambda e_\lambda|_{\lambda=0}
+  =\rho^2\left(
+    -\int\frac{dk}{(2\pi)^3}\ g(|k|)\hat u(|k|)
+    -\int\frac{dk}{(2\pi)^3}\ \hat v(|k|)\partial_\lambda\hat u_\lambda(|k|)|_{\lambda=0}
+  \right)
+  .
+\end{equation}
+To compute $\partial_\lambda\hat u_\lambda|_{\lambda=0}$, we differentiate $\Xi(\mathbb U,\lambda)=0$:
+\begin{equation}
+  \partial_\lambda\mathbb U|_{\lambda=0}=-(D\Xi)^{-1}\partial_\lambda\Xi|_{\lambda=0}
+  .
+\end{equation}
+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}:
+\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\frac{dk}{(2\pi)^3}\ \hat v(k)\partial_\lambda\hat u_\lambda(k)|_{\lambda=0}
+  .
+\end{equation}
+\bigskip
+
+\point
+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}
+The presence of delta functions does not play well with the quadratures.
+To get around this, we instead compute a regularization of $\mathcal M(q)$ by convolving it with a peaked spherically symmetric function.
+Let $\Gamma_L$ denote the Gaussian with variance $1/\sqrt L$
+\begin{equation}
+  \Gamma_L(|k|):=\left(\frac{L}{2\pi}\right)^{\frac32}e^{-\frac L2k^2}
+  .
+  \label{easyeq_gaussian_momt}
+\end{equation}
+In fact, we will scale $L$ with $k$, and set $L$ to
+\begin{equation}
+  L=\sqrt{\mathrm{\refname{param:easyeq_momentum_window_L}{{\tt window\_L}}}}/k^2
+  .
+\end{equation}
+To compute
+\begin{equation}
+  \mathfrak M(q):=\mathcal M\ast\Gamma_L(q)
+\end{equation}
+we solve the equation
+\begin{equation}
+  -\Delta u_\lambda(|x|)
+  =
+  (1-u_\lambda(|x|))v(|x|)-2\rho K(|x|)+\rho^2 L(|x|)
+  -2\lambda\int dk\ \hat u_0(k)\cos(k\cdot x)\Gamma_L(q-k)
+  .
+\end{equation}
+Note that the Fourier transform of
+\begin{equation}
+  -2\lambda\int dk\ \hat u_0(k)\cos(k\cdot x)\Gamma_L(q-k)
+\end{equation}
+is
+\begin{equation}
+  -(2\pi)^3\lambda\hat u_0(q)(\Gamma_L(k+q)+\Gamma_L(k-q))
+  .
+\end{equation}
+Since the ground state is unique, $\mathcal M$ is spherically symmetric.
+The term $\Gamma_L(k\pm q)$ is not, so we take its spherical average (which will not change the final result): by lemma\-~\ref{lemma:bipolar},
+\begin{equation}
+  -\frac1{4\pi r^2}\int dq\ \delta(|q|-r)(2\pi)^3\lambda\hat u_0(q)(\Gamma_L(k+q)+\Gamma_L(k-q))
+  =
+  -\frac{(2\pi)^3}{|k|r}\lambda\hat u_0(r)\int_{||k|-r|}^{|k|+r} ds\ s\Gamma_L(s)
+  .
+\end{equation}
+In this setup, the approximation of the delta function is thus
+\begin{equation}
+  \tilde\delta(|k|,r):=
+  \frac1{2|k|r}\int_{||k|-r|}^{|k|+r} ds\ s\Gamma_L(s)
+  =
+  \frac{1}{2|k|rL}(\Gamma_L(|k|-r)-\Gamma_L(|k|+r))
+  .
+\end{equation}
+\bigskip
+
+\point
+To compute the momentum distribution at $q$, we define $\Xi(\mathbb U,\lambda)$ by replacing $\mathbb T$ with
+\begin{equation}
+  \mathbb T_{i}=
+  \frac1{\mathbb A_{K,i}}
+  \left(
+    \mathbb S_i
+    -2(2\pi)^3\lambda \hat u(|q|)\tilde\delta(k_i,|q|)
+  \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_{i}|_{\lambda=0}
+  =
+  -\partial_\lambda\mathbb T_{i}|_{\lambda=0}\left(
+    \frac1{2(\mathbb X_{i}+1)}
+    \Phi\left(\mathbb B_i{\textstyle\frac{\mathbb T_{i}}{(\mathbb X_{i}+1)^2}}\right)
+    +\frac{\mathbb B_i\mathbb T_{i}}{2(\mathbb X_{i}+1)^3}
+    \partial\Phi\left(\mathbb B_i{\textstyle\frac{\mathbb T_{i}}{(\mathbb X_{l,j}+1)^2}}\right)
+  \right)
+\end{equation}
+with
+\begin{equation}
+  \partial_\lambda\mathbb T_{i}|_{\lambda=0}=
+  -\frac{2(2\pi)^3}{\mathbb A_{K,i}} \hat u(|q|)\tilde\delta(k_{i},|q|)
+  .
+\end{equation}
+\bigskip
+
 \subsection{\tt anyeq}\label{sec:anyeq}
 \indent
 This method is used to solve any of the equations in the Simplified approach.
@@ -598,11 +1053,11 @@ Unless otherwise noted, this method takes the following parameters (specified vi
   \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.
+  {\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 $\neq 0$, 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}}.
 
   \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}$.
@@ -623,7 +1078,7 @@ The available {\tt commands} are the following.
   \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{Disabled parameters}: \refname{param:anyeq_rho}{{\tt rho}}.\par
   \underline{Extra parameters}:
   \begin{itemize}
     \item\makelink{param:anyeq_minlrho}{}
@@ -635,8 +1090,18 @@ The available {\tt commands} are the following.
     \item\makelink{param:anyeq_nlrho}{}
     {\tt nlrho} ({\tt Int64}, default: $100$): number of values for $\rho$ (spaced logarithmically).
 
+    \item\makelink{param:anyeq_minrho}{}
+    {\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.
+
     \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.
+    {\tt rhos} ({\tt Array\{Float64\}}): list of values for $\rho$, specified as a `{\tt,}' separated list.
+    This parameter takes precedence over \refname{param:anyeq_minlrho}{{\tt minlrho}}, \refname{param:anyeq_maxlrho}{{\tt maxlrho}}, \refname{param:anyeq_nlrho}{{\tt nlrho}}, \refname{param:anyeq_minrho}{{\tt minrho}}, \refname{param:anyeq_maxrho}{{\tt maxrho}}, \refname{param:anyeq_nrho}{{\tt nrho}}.
   \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.
@@ -686,14 +1151,10 @@ The available {\tt commands} are the following.
     {\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|)$]
+  \underline{Output} (one line for each value of $x$): [$|x|$] [$u(|x|)$]
 
-  \item\makelink{command_anyeq_2pt}{}
-  {\tt 2pt}: compute the two-point correlation function $C_2(|x|)$ at a given $\rho$.\par
+  \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}{}
@@ -701,16 +1162,99 @@ The available {\tt commands} are the following.
   \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.
+  \end{itemize}
+  \underline{Output}: [$|x_{\mathrm{max}}|$] [$C_2(|x_{\mathrm{max}}|)$]
+
+  \item\makelink{command:anyeq_2pt_max_rho}{}
+  {\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
+  \underline{Extra parameters}:
+  \begin{itemize}
+    \item\makelink{param:anyeq_kmin}{}
+    {\tt kmin} ({\tt Float64}, default: 0): minimum of the range of $|k|$ to be printed.
+
+    \item\makelink{param:anyeq_kmax}{}
+    {\tt kmax} ({\tt Float64}, default: 10): maximum of the range of $|k|$ to be printed.
+
+    \item\makelink{param:anyeq_nk}{}
+    {\tt nk} ({\tt Int64}, default: 100): number of $|k|$'s to be printed.
+
+    \item\makelink{param:anyeq_2pt_fourier_window_L}{}
+    {\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}).
+  \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:anyeq_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:anyeq_2pt_fourier_window_L}{{\tt window\_L}}, \refname{param:anyeq_maxstep}{{\tt maxstep}} plus
+  \begin{itemize}
+    \item\makelink{param:anyeq_dk}{}
+    {\tt dk} ({\tt Float64}, default: $10^{-7}$): step used to numerically approximate derivatives.
+
+    \item\makelink{param:anyeq_k0}{}
+    {\tt k0} ({\tt Float64}, default: $1$): initial guess for the maximum is $\rho^{1/3}\mathrm{\tt k0}$.
+  \end{itemize}
+  \underline{Output}: [$|k_{\mathrm{max}}|$] [$\hat C_2(|k_{\mathrm{max}}|)$]\par
+
+  \item\makelink{command:anyeq_2pt_fourier_max_rho}{}
+  {\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.
+
+  \item\makelink{command:anyeq_momentum_distribution}{}
+  {\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:anyeq_momentum_kmin}{}
+    {\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}).
+  \end{itemize}
+  \underline{Output} (one line for each value of $|k|$): [$|k|$] [$\mathcal M(|k|)$]
+
+  \item\makelink{command:anyeq_compressibility_rho}{}
+  {\tt compressibility\_rho}: compute the compressibility $\chi$ as a function of $\rho$.
+  The Newton algorithm is initialized with the solution of \refname{eq:easyeq_medeq}{{\tt medeq}}.
+  \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$] [$\chi$].
+
   \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{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_init}{{\tt nlrho\_init}}.\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.}
+\subsubsubsection{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}
@@ -770,11 +1314,15 @@ with
 \end{equation}
 Therefore,
 \begin{equation}
+  \Omega(u)=0
+\end{equation}
+with
+\begin{equation}
+  \Omega(u):=
   \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
+  \label{Omega}
 \end{equation}
 with
 \begin{equation}
@@ -831,7 +1379,7 @@ with
 \end{equation}
 \bigskip
 
-\point{\bf Evaluating integrals.}
+\subsubsubsection{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}
@@ -861,7 +1409,7 @@ In these terms,
   \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}})
+  \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$:
@@ -879,6 +1427,7 @@ and find
   \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))
   .
+  \label{anyeq_preintegration}
 \end{equation}
 We then approximate the integral using a Gauss-Legendre quadrature of order $N$ (see appendix\-~\ref{appGL}):
 \begin{equation}
@@ -901,14 +1450,14 @@ The choice of the change of variables\-~(\ref{chebyshevvars}) is made so that $U
 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.}
+\subsubsubsection{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
+\point
 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.
@@ -934,7 +1483,7 @@ in which $\mathds 1_{\tau_l<\tau\leqslant\tau_{l+1}}\in\{0,1\}$ is equal to 1 if
 (see appendix\-~\ref{appChebyshev} for a discussion of the Chebyshev polynomial expansion and the error of its truncations).
 \bigskip
 
-\subpoint
+\point
 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}
@@ -952,7 +1501,7 @@ with
 and $(x_j,w_j)$ are the abcissa and weights for Gauss-Legendre quadratures.
 \bigskip
 
-\subpoint
+\point
 All in all, we approximate
 \begin{equation}
   a({\textstyle\frac{1-\tau}{1+\tau}})
@@ -964,14 +1513,15 @@ Furthermore, using the Chebyshev polynomial expansion and Gauss-Legendre quadrat
 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.}
+\subsubsubsection{Convolutions}
 Using the Chebyshev polynomial approximation, we can compute convolutions as follows.
 \medskip
 
-\subpoint
+\point
 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)
+  \label{conv_approx1}
 \end{equation}
 We change variables to compactify the integrals $\tau=\frac{1-t}{1+t}$, $\sigma=\frac{1-s}{1+s}$:
 \begin{equation}
@@ -990,30 +1540,31 @@ Therefore, using the approximation\-~(\ref{approxchebyshev}), if $\mathfrak a_{l
 \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}|)
+  \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}})
+    A_{l',n;l'',m}^{(\nu_a,\nu_b)}(|k|):=
+    \frac1{4\pi^2|k_{l,j}|}
+    \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'}}})
+    \mathds 1_{\alpha_-(|k|,\tau)<\tau_{l''+1}}
+    \mathds 1_{\alpha_+(|k|,\tau)>\tau_{l''}}
+    \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}
+  \label{conv_approx2}
 \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
+\point
 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
@@ -1047,7 +1598,7 @@ 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
+\point
 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}
@@ -1077,11 +1628,11 @@ We can then evaluate convolutions at $k=0$:
 \bigskip
 
 
-\subpoint
+\point
 Let us now compute some choices of $\mathfrak a,\mathfrak b$ more explicitly.
 \medskip
 
-\subsubpoint
+\subpoint
 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}
@@ -1095,7 +1646,7 @@ so
 \end{equation}
 \bigskip
 
-\subsubpoint
+\subpoint
 We now turn to $a\hat\ast\hat v$.
 Let
 \begin{equation}
@@ -1132,7 +1683,7 @@ The integral in\-~(\ref{Upsilon}) is computed using Gauss-Legendre quadratures,
 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
+\subpoint
 Finally,
 \begin{equation}
   (\mathds 1_{l',n}\odot\mathfrak v)_{l,i}=
@@ -1146,7 +1697,7 @@ and
 \end{equation}
 \bigskip
 
-\point{\bf Evaluating $\hat l_3$.}
+\subsubsubsection{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}
@@ -1230,12 +1781,12 @@ Since the tensor $\bar A$ is quite large (it contains $(NJ)^5$ entries), and its
 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$.}
+\subsubsubsection{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
+\point
 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}:=
@@ -1283,7 +1834,7 @@ The equation for $(\mathbb U_{l,j})_{l\in\{0,\cdots,J-1\},j\in\{1,\cdots,N\}}$ i
 (see\-~(\ref{odot}), (\ref{otimes}) and\-~(\ref{avg}) for the definitions of $\odot$, $\otimes$ and $\left<\cdot\right>$).
 \bigskip
 
-\subpoint
+\point
 We rewrite\-~(\ref{bbU}) as a root finding problem:
 \begin{equation}
   \Xi_{l,j}(\mathbb U):=
@@ -1303,14 +1854,14 @@ where $D\Xi$ is the Jacobian of $\Xi$:
 \end{equation}
 \bigskip
 
-\subpoint
+\point
 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}}.
+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}}.
 \bigskip
 
-\subpoint
+\point
 We are left with computing the Jacobian of $\Xi$:
 \begin{equation}
   \begin{largearray}
@@ -1444,7 +1995,7 @@ with
 \end{equation}
 \bigskip
 
-\subpoint
+\point
 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}
@@ -1471,8 +2022,8 @@ where $k_{l,j}$ was defined in\-~(\ref{klj}), and the solution $u$ in real space
 \end{equation}
 \bigskip
 
-\point{\bf Condensate fraction.}
-To compute the condensate fraction, we solve the modified {\tt anyeq} (see\-~\cite{CJL20b}):
+\subsubsubsection{Condensate fraction}
+To compute the condensate fraction, we solve the modified {\tt anyeq} (see\-~\cite{CJL21}):
 \begin{equation}
   (-\Delta+2\mu)u_\mu
   =
@@ -1519,41 +2070,53 @@ We then approximate
 (see\-~(\ref{avg})).
 \bigskip
 
-\point{\bf Correlation function.}
+\subsubsubsection{Correlation function (spherical average)}
 The two-point correlation function is
 \begin{equation}
-  C_2(x):=
+  c_2(x):=
   2\rho\frac{\delta e}{\delta v(x)}
+\end{equation}
+and its spherical average is
+\begin{equation}
+  C_2(|x|):=\frac1{4\pi|x|^2}\int dy\ \delta(|x|-|y|)c_2(y)
   .
 \end{equation}
 In Fourier space,
 \begin{equation}
-  C_2(x)=
-  2\rho\int dk\ e^{ikx}\frac{\delta e}{\delta\hat v(k)}
+  c_2(x)=
+  2\rho
+  \int dk\ e^{ikx}\frac{\delta e}{\delta\hat v(k)}
+\end{equation}
+so
+\begin{equation}
+  C_2(|x|)=
+  2\rho\int dk\ \left(\frac1{4\pi|x|^2}\int dy\ \delta(|x|-|y|)e^{iky}\right)\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
+\point
 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|}
   .
+  \label{2pt_addv}
 \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)}
+  =\int dk\ \frac{\delta e}{\delta\hat v(k)}\partial_\lambda({\textstyle \hat v(k)+\lambda g(|k|)})
+  =\int dk\ g(|k|)\frac{\delta e}{\delta\hat v(k)}
+  .
 \end{equation}
-so, denoting the solution of the modified equation by $u_\lambda$,
+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})
+  =\rho^2g(0)-\rho^2\int\frac{dk}{(2\pi)^3}\ (g(k)\hat u(k)+\hat v(k)\partial_\lambda\hat u_\lambda(k)|_{\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:
@@ -1563,7 +2126,7 @@ We compute $\partial_\lambda u_\lambda|_{\lambda=0}$ in the same way as the unco
 \end{equation}
 \bigskip
 
-\subpoint
+\point
 We compute $\partial_\lambda\Xi|_{\lambda=0}$:
 \begin{equation}
   \begin{largearray}
@@ -1581,6 +2144,7 @@ We compute $\partial_\lambda\Xi|_{\lambda=0}$:
       -\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}
+  \label{2pt_dXi}
 \end{equation}
 with
 \begin{equation}
@@ -1610,7 +2174,7 @@ with
 \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}
+    \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)
@@ -1623,7 +2187,7 @@ with
 \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\beta_K\left(\mathbb U\odot\left(\partial_\lambda\mathbb S\mathbb U\right)\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}
@@ -1641,7 +2205,7 @@ To evaluate $(\mathbb U\odot\mathfrak g)$ and $\left<\mathbb U\mathfrak g\right>
 To do so, we replace $\hat v$ with $g$ in the computation of $\Upsilon$.
 \bigskip
 
-\subpoint
+\point
 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}})
@@ -1652,7 +2216,203 @@ 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.}
+\point
+To compute the maximum of $C_2$, we use a modified Newton algorithm.
+The initial guess for the maximum is $|x_0|=\rho^{-\frac13}$\refname{param:anyeq_x0}{{\tt x0}}.
+The modified Newton algorithm is an iteration:
+\begin{equation}
+  x_{n+1}=x_n+\frac{\partial C_2(|x_n|)}{|\partial^2C_2(|x_n|)|}
+  \label{anyeq_newton_2pt}
+\end{equation}
+in which the derivatives are approximated using finite differences:
+\begin{equation}
+  \partial C_2(x)\approx \frac{C_2(|x|+dx)-C_2(|x|)}{dx}
+  ,\quad
+  \partial^2 C_2(x)\approx \frac{C_2(|x|+dx)+C_2(|x|-dx)-2C_2(|x|)}{dx^2}
+  .
+  \label{anyeq_dx_2pt}
+\end{equation}
+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
+\begin{equation}
+  \hat c_2(q):=
+  2\rho\frac{\delta e}{\delta v(q)}
+\end{equation}
+and its spherical average is
+\begin{equation}
+  \hat C_2(|q|):=\frac1{4\pi|q|^2}\int dk\ \delta(|q|-|k|)c_2(k)
+  =
+  \frac\rho{2\pi|q|^2}\int dk\ \delta(|q|-|k|)\frac{\delta e}{\delta\hat v(k)}
+  .
+\end{equation}
+\bigskip
+
+\point
+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
+\begin{equation}
+  \Gamma_L(|k|):=\left(\frac{L}{2\pi}\right)^{\frac32}e^{-\frac L2k^2}
+  \label{anyeq_gaussian}
+\end{equation}
+\begin{equation}
+  \hat{\mathfrak C}_2(|q|)
+  :=
+  \int dp\ \hat C_2(|q-p|)\Gamma_L(|p|)
+  =\int dk\ \int dp\ \frac\rho{2\pi|q-p|^2}\delta(|q-p|-|k|)\frac{\delta e}{\delta\hat v(k)}\Gamma_L(|p|)
+\end{equation}
+which by lemma\-~\ref{lemma:bipolar} is
+\begin{equation}
+  \hat{\mathfrak C}_2(|q|)
+  =\int dk\ 
+  \frac\rho{|q|}
+  \frac{\delta e}{\delta\hat v(k)}
+  \int_0^\infty dt\ \int_{||q|-t|}^{|q|+t}ds\ s\frac{\delta(t-|k|)}t\Gamma_L(s)
+\end{equation}
+that is
+\begin{equation}
+  \hat{\mathfrak C}_2(|q|)
+  =\int dk\ 
+  \frac{\delta e}{\delta\hat v(k)}
+  \frac{\rho}{|q||k|}
+  \int_{||q|-|k||}^{|q|+|k|}ds\ s\Gamma_L(s)
+\end{equation}
+which is the directional derivative of $e$ with respect to $\hat v$ in the direction of $2\rho g$ with
+\begin{equation}
+  g(|k|):=\frac1{2|q||k|}\int_{||q|-|k||}^{|q|+|k|}ds\ s\Gamma_L(s)
+  =
+  \frac1{2|k|rL}(\Gamma_L(|k|-r)-\Gamma_L(|k|+r))
+  .
+  \label{anyeq_2pt_fourier_g}
+\end{equation}
+Note that
+\begin{equation}
+  g(0):=\Gamma_L(|q|)
+  .
+\end{equation}
+To compute this derivative, we replace $\hat v$ with
+\begin{equation}
+  \hat v+\lambda g(|k|)
+\end{equation}
+so, denoting the solution of the modified equation by $u_\lambda$, for $q\neq 0$,
+\begin{equation}
+  \hat{\mathfrak C}_2(|q|)=2\rho\partial_\lambda e_\lambda|_{\lambda=0}
+  =\rho^2\left(
+    -\int\frac{dk}{(2\pi)^3}\ g(|k|)\hat u(|k|)
+    -\int\frac{dk}{(2\pi)^3}\ \hat v(|k|)\partial_\lambda\hat u_\lambda(|k|)|_{\lambda=0}
+  \right)
+  .
+\end{equation}
+To compute $\partial_\lambda\hat u_\lambda|_{\lambda=0}$, we differentiate $\Xi(\mathbb U,\lambda)=0$:
+\begin{equation}
+  \partial_\lambda\mathbb U|_{\lambda=0}=-(D\Xi)^{-1}\partial_\lambda\Xi|_{\lambda=0}
+  .
+\end{equation}
+The computation of $\partial_\lambda\Xi|_{\lambda=0}$ is identical to\-~(\ref{2pt_dXi}), but taking
+\begin{equation}
+  \mathfrak g_{l,j}=g(|k_{l,j}|)
+  .
+\end{equation}
+with $g$ defined in\-~(\ref{anyeq_2pt_fourier_g}).
+\bigskip
+
+\point
+To compute the maximum of $\hat C_2$, we proceed as for $C_2$, see\-~(\ref{anyeq_newton_2pt})-(\ref{anyeq_dx_2pt}).
+The only difference is that the algorithm is initialized with $|k_0|=\rho^{\frac13}$\refname{param:anyeq_k0}{{\tt k0}}.
+
+\subsubsubsection{Correlation function of uncondensed particles (spherical average)}
+To compute the correlation function among uncondensed particles, denoted by $\gamma_2(|\xi|)$, we solve the modified {\tt anyeq} (see\-~\cite{Ja23}):
+\begin{equation}
+  -\Delta u_\mu
+  =
+  (1-u_\mu)v-2\rho K+\rho^2 L
+  -\frac{\rho\mu}{2\pi|\xi|^2} u(|\xi|)\delta(|\xi|-|x|)
+\end{equation}
+where $K$ and $L$ are defined as in\-~(\ref{anyeq_K})-(\ref{anyeq_e}) in which $u$ is replaced with $u_\mu$.
+In Fourier space,
+\begin{equation}
+  k^2 \hat u_\mu
+  =
+  \hat S-2\rho\hat K+\rho^2\hat L
+  -2\rho\mu u(|\xi|)\frac{\sin(|k||\xi|)}{|k||\xi|}
+  .
+\end{equation}
+The uncondensed correlation function is then
+\begin{equation}
+  \gamma_2(|\xi|)=\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 the term
+\begin{equation}
+  \mu g(|k|):=
+  -\mu 2\rho u(|\xi|)\frac{\sin(|k||\xi|)}{|k||\xi|}
+\end{equation}
+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
+\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}+\mu g(k_{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{\partial_\mu\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)
+  -\frac{(1+\mathbb Y_{l,j})\partial_\mu\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)
+  ,\quad
+  \partial_\mu\mathbb Y_{l,j}=\frac{g(k_{l,j})}{\mathbb L_{l,j}}
+  .
+\end{equation}
+We then approximate
+\begin{equation}
+  \gamma_2(\xi)\approx
+  -\frac12\left<\mathfrak v\partial_\mu\mathbb U\right>
+\end{equation}
+(see\-~(\ref{avg})).
+\bigskip
+
+\indent
+In order to avoid numerical oscillations due to the $\sin$ function, 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)g(k)$.
+
+\subsubsubsection{Momentum distribution}
 To compute the momentum distribution (see\-~\cite{CHe21}), we add a parameter $\lambda$ to {\tt anyeq}:
 \begin{equation}
   -\Delta u_\lambda(|x|)
@@ -1664,7 +2424,7 @@ To compute the momentum distribution (see\-~\cite{CHe21}), we add a parameter $\
 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}
+  =-\frac\rho2\int\frac{dk}{(2\pi)^3}\ \hat v(k)\partial_\lambda\hat u_\lambda(k)|_{\lambda=0}
   .
 \end{equation}
 Note that the Fourier transform of $2\lambda\hat u_0(q)\cos(q\cdot x)$ is
@@ -1673,52 +2433,118 @@ Note that the Fourier transform of $2\lambda\hat u_0(q)\cos(q\cdot x)$ is
   .
   \label{fouriermomentum}
 \end{equation}
-We compute $\partial_\lambda u_\lambda|_{\lambda=0}$ in the same way as the uncondensed fraction.
+The presence of delta functions does not play well with the Chebyshev polynomial expansion and the quadratures.
+\bigskip
+
+\point
+We will consider two different ways of getting around this.
 \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
+One idea is to compute a regularization of $\mathcal M(q)$ by convolving it with a peaked spherically symmetric function.
+Let $\Gamma_L$ denote the Gaussian with variance $1/\sqrt L$:
 \begin{equation}
-  \int dk f(k)\delta(k-q)
-  =f(q)
+  \Gamma_L(|k|):=\left(\frac{L}{2\pi}\right)^{\frac32}e^{-\frac L2k^2}
+  .
+  \label{anyeq_gaussian_momt}
+\end{equation}
+In fact, we will scale $L$ with $k$, and set $L$ to
+\begin{equation}
+  L=\sqrt{\mathrm{\refname{param:anyeq_momentum_window_L}{{\tt window\_L}}}}/k^2
+  .
+\end{equation}
+To compute
+\begin{equation}
+  \mathfrak M(q):=\mathcal M\ast\Gamma_L(q)
+\end{equation}
+we solve the equation
+\begin{equation}
+  -\Delta u_\lambda(|x|)
+  =
+  (1-u_\lambda(|x|))v(|x|)-2\rho K(|x|)+\rho^2 L(|x|)
+  -2\lambda\int dk\ \hat u_0(k)\cos(k\cdot x)\Gamma_L(q-k)
+  .
+\end{equation}
+Note that the Fourier transform of
+\begin{equation}
+  -2\lambda\int dk\ \hat u_0(k)\cos(k\cdot x)\Gamma_L(q-k)
+\end{equation}
+is
+\begin{equation}
+  -(2\pi)^3\lambda\hat u_0(q)(\Gamma_L(k+q)+\Gamma_L(k-q))
+  .
+\end{equation}
+Since the ground state is unique, $\mathcal M$ is spherically symmetric.
+The term $\Gamma_L(k\pm q)$ is not, so we take its spherical average (which will not change the final result): by lemma\-~\ref{lemma:bipolar},
+\begin{equation}
+  -\frac1{4\pi r^2}\int dq\ \delta(|q|-r)(2\pi)^3\lambda\hat u_0(q)(\Gamma_L(k+q)+\Gamma_L(k-q))
+  =
+  -\frac{(2\pi)^3}{|k|r}\lambda\hat u_0(r)\int_{||k|-r|}^{|k|+r} ds\ s\Gamma_L(s)
+  .
+\end{equation}
+In this setup, the approximation of the delta function is thus
+\begin{equation}
+  \tilde\delta(|k|,r):=
+  \frac1{2|k|r}\int_{||k|-r|}^{|k|+r} ds\ s\Gamma_L(s)
+  =
+  \frac{1}{2|k|rL}(\Gamma_L(|k|-r)-\Gamma_L(|k|+r))
+  .
+\end{equation}
+To choose this method, set \refname{param:anyeq_window_L}{{\tt window\_L}} to a finite value.
+\bigskip
+
+\subpoint
+Another approach is to contruct a discrete analog of the delta-functions in\-~(\ref{fouriermomentum}).
+The starting point we take is, for $q\neq0$,
+\begin{equation}
+  \int dk\ f(|k|)\delta(k-q)
+  \equiv\frac1{4\pi|q|^2}\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}
+  \frac{\pi}{8|q|^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})
+  \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),|q|)
+  =
+  f(|q|)
 \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
+(see\-~(\ref{chebyshevvars})), we define the approximation of the delta function as
 \begin{equation}
-  \tilde\delta_{l,j;l',i}
+  \tilde\delta(k_{l,j},k_{l',i})
   :=
   \delta_{l,l'}\delta_{j,i}
-  \frac2{\pi^2}
+  \frac{8}{\pi}
   \left(
     (\tau_{l+1}-\tau_l)
     w_j
     \cos({\textstyle\frac{\pi x_j}2})
-    (1+k_{l,j})^2k_{l,j}^2
+    (1+k_{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)$.
+To choose this method, set \refname{param:anyeq_window_L}{{\tt window\_L}} to {\tt Inf}.
+{\color{red}This method seems to yield some fairly poor results!}
 \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
+\point
+To compute the momentum distribution at $q$, we define $\Xi(\mathbb U,\lambda)$ by formally adding
+\begin{equation}
+  -2(2\pi)^3\lambda \hat u_0(|q|)\tilde\delta(k_{l,j},|q|)
+\end{equation}
+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)
+  \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 \hat u(|q|)\tilde\delta(k_{l,j},|q|)
+  \right)
+  .
 \end{equation}
 Then we solve $\Xi(\mathbb U,\lambda)=0$, and differentiate:
 \begin{equation}
@@ -1739,7 +2565,33 @@ Finally,
 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}}
+  -\frac{2(2\pi)^3}{\mathbb L_{l,j}} \hat u(|q|)\tilde\delta(k_{l,j},|q|)
+  .
+\end{equation}
+\bigskip
+
+\subsubsubsection{Compressibility}
+The compressibility is defined as
+\begin{equation}
+  \chi:=\frac1{\rho^2\partial_\rho^2(\rho e)}
+  =\frac1{\partial_{\log\rho}^2(\rho e)-\partial_{\log\rho}(\rho e)}
+  .
+\end{equation}
+We approximate these derivatives by finite differences:
+\begin{equation}
+  \partial_{\log\rho}^2 f(\rho)
+  \approx
+  \frac{f(\rho_{j+1})+f(\rho_{j-1})-2f(\rho_j)}{(\log\rho_{j+2}-\log\rho_{j+1})(\log\rho_{j+1}-\log\rho_j)}
+\end{equation}
+and
+\begin{equation}
+  \partial_{\log\rho} f(\rho)
+  \approx
+  \frac12\left(
+    \frac{f(\rho_{j+1})-f(\rho_j)}{\log\rho_{j+1}-\log\rho_{j}}
+    +
+    \frac{f(\rho_{j})-f(\rho_{j-1})}{\log\rho_{j}-\log\rho_{j-1}}
+  \right)
   .
 \end{equation}
 
@@ -1755,7 +2607,7 @@ The method is used to compute observables for the simple equation
   .
   \label{simpleq}
 \end{equation}
-One can show\-~\cite[Theorem\-~1.6]{CJL20b} that the condensate fraction is
+One can show\-~\cite[Theorem\-~1.6]{CJL21} 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}
@@ -1824,8 +2676,18 @@ The available {\tt commands} are the following.
     \item\makelink{param:simpleq-Kv_nlrho}{}
     {\tt nlrho} ({\tt Int64}, default: $100$): number of values for $\rho$ (spaced logarithmically).
 
+    \item\makelink{param:simpleq-Kv_minrho}{}
+    {\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.
+
     \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.
+    {\tt rhos} ({\tt Array\{Float64\}}): list of values for $\rho$, specified as a `{\tt,}' separated list.
+    This parameter takes precedence over \refname{param:simpleq-Kv_minlrho}{{\tt minlrho}}, \refname{param:simpleq-Kv_maxlrho}{{\tt maxlrho}}, \refname{param:simpleq-Kv_nlrho}{{\tt nlrho}}, \refname{param:simpleq-Kv_minrho}{{\tt minrho}}, \refname{param:simpleq-Kv_maxrho}{{\tt maxrho}}, \refname{param:simpleq-Kv_nrho}{{\tt nrho}}.
   \end{itemize}
   \underline{Output} (one line for each value of $\rho$): [$\rho$] [$\eta$]
 
@@ -1980,8 +2842,18 @@ The available {\tt commands} are the following.
     \item\makelink{param:simpleq-hardcore_nlrho}{}
     {\tt nlrho} ({\tt Int64}, default: $100$): number of values for $\rho$ (spaced logarithmically).
 
+    \item\makelink{param:simpleq-hardcore_minrho}{}
+    {\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.
+
     \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.
+    {\tt rhos} ({\tt Array\{Float64\}}): list of values for $\rho$, specified as a `{\tt,}' separated list.
+    This parameter takes precedence over \refname{param:simpleq-hardcore_minlrho}{{\tt minlrho}}, \refname{param:simpleq-hardcore_maxlrho}{{\tt maxlrho}}, \refname{param:simpleq-hardcore_nlrho}{{\tt nlrho}}, \refname{param:simpleq-hardcore_minrho}{{\tt minrho}}, \refname{param:simpleq-hardcore_maxrho}{{\tt maxrho}}, \refname{param:simpleq-hardcore_nrho}{{\tt nrho}}.
   \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.
@@ -2023,7 +2895,7 @@ In order to carry out the computation of the solution of\-~(\ref{linearhc}) and
 for $|x|>1$.
 \bigskip
 
-\point{\bf Energy.}
+\subsubsubsection{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
@@ -2056,7 +2928,7 @@ Therefore,
 \end{equation}
 \bigskip
 
-\point{\bf Integral equation.}
+\subsubsubsection{Integral equation}
 We turn the differential equation in\-~(\ref{hardcore_simple}) into an integral equation.
 Let
 \begin{equation}
@@ -2115,7 +2987,7 @@ We change variables in the last integral:
 \end{equation}
 \bigskip
 
-\point{\bf The auto-convolution term.} We split
+\subsubsubsection{The auto-convolution term} We split
 \begin{equation}
   u(r)=\mathds 1_{r>1}u_+(r-1)+\mathds 1_{r\leqslant 1}
   .
@@ -2253,7 +3125,7 @@ where
 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})
+\subsubsubsection{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}
@@ -2343,10 +3215,10 @@ and
 \end{equation}
 \bigskip
 
-\point{\bf Energy.} We now compute the approximation for the energy, using\-~(\ref{hardcore_energy}).
+\subsubsubsection{Energy} We now compute the approximation for the energy, using\-~(\ref{hardcore_energy}).
 \bigskip
 
-\subpoint We start with $\partial u|_{|x|\searrow 1}$.
+\point We start with $\partial u|_{|x|\searrow 1}$.
 By\-~(\ref{uinteq}),
 \begin{equation}
   \partial u|_{|x|\searrow 1}
@@ -2386,7 +3258,7 @@ and
   .
 \end{equation}
 
-\subpoint Let us now turn to $\int_{|x|<1}dx\ u\ast u(x)$.
+\point 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)
@@ -2421,7 +3293,7 @@ and
 \end{equation}
 \bigskip
 
-\subpoint Thus,
+\point Thus,
 \begin{equation}
   \begin{largearray}
     e=2\pi\rho
@@ -2452,11 +3324,11 @@ and
 \end{equation}
 \bigskip
 
-\point{\bf Newton algorithm.}
+\subsubsubsection{Newton algorithm}
 In this paragraph, we set $\epsilon=e$, that is, $\mu=0$.
 \bigskip
 
-\subpoint
+\point
 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}
@@ -2519,7 +3391,7 @@ Gauss-Legendre and Gauss-Laguerre quadratures and their errors are discussed in
 The orders of the quadratures are given by the variable \refname{param:simpleq-hardcore_N}{{\tt N}}.
 \bigskip
 
-\subpoint
+\point
 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)})
@@ -2538,7 +3410,7 @@ We initialize the algorithm with
 \end{equation}
 \bigskip
 
-\subpoint
+\point
 We are left with computing the Jacobian of $\Xi$:
 \begin{equation}
   \begin{largearray}
@@ -2709,7 +3581,7 @@ and
 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
+\point
 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}
@@ -2723,7 +3595,7 @@ The energy thus obtained is
 \end{equation}
 \bigskip
 
-\point{\bf Condensate fraction.}
+\subsubsubsection{Condensate fraction}
 To compute the condensate fraction, we use the parameter $\mu$ in\-~(\ref{hardcore_simple}).
 The uncondensed fraction is
 \begin{equation}
@@ -3213,7 +4085,7 @@ In $x$-space
   \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}
+Note that $v\geqslant 0$ if and only if\-~\cite{CJL21}
 \begin{equation}
   b>0
   ,\quad