diff options
Diffstat (limited to 'doc/simplesolv-doc.tex')
| -rw-r--r-- | doc/simplesolv-doc.tex | 1110 |
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 |