Ian Jauslin
summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorIan Jauslin <ian.jauslin@rutgers.edu>2023-02-26 18:36:05 -0500
committerIan Jauslin <ian.jauslin@rutgers.edu>2023-02-26 18:36:05 -0500
commitc1b477a1b2b796617c4e345a7296a8d429d7a067 (patch)
tree8a8a2fc0fb6e7da5f4b0b271382740f858ee4372 /doc/simplesolv-doc.tex
parente72af82c3ed16b81cdb5043c58abbdbb3cf02102 (diff)
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.
Diffstat (limited to 'doc/simplesolv-doc.tex')
-rw-r--r--doc/simplesolv-doc.tex1110
1 files changed, 991 insertions, 119 deletions
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