Complete Algorithm

Let us record here the complete algorithm, with all the details. Recall that the integrals in the $\gamma$ variable should not be used in the cases $p=0$ and $p=2$, in which cases we should use the integrals in terms of the variable $\chi$. For all other values of $p$ we should use the integrals in terms of the variable $\gamma$ for ease of numerical integration.

1. Choose a value for $\xi$ in $(-\infty,\infty)$. The method for this is: choose some typically small $\varepsilon$ and integers $n=1,2,3,\ldots$, then use the values

   
   

   $$\xi_{n}=\frac{\mathfrak{N}}{4d}\sqrt{n\varepsilon}-\frac{d}{\sqrt{n\varepsilon}},$$
   
   

   where $n\varepsilon$ is a very rough estimate of the corresponding value of $\lambda$. For example, for $\mathfrak{N}=1$ using $\varepsilon=0.1$ and $n$ up to $1000$ results on $\lambda$ up to approximately $40$, and using $\varepsilon=1.2$ and $n$ up to $1000$ results on $\lambda$ up to approximately $1000$.

2. Calculate the corresponding value of $\sqrt{\lambda}$ by means of two numerical integrations, in one of the two forms,

   
   

   $$\sqrt{\lambda} =\frac{4d}{\mathfrak{N}}\;\frac{\displaystyle... ...ma\;\gamma^{\frac{\mathfrak{N}-2}{2}}\;e^{-(\gamma-\xi)^{2}}},$$
   
   

   where the integrals in $\chi$ and in $\gamma$ are related by

   
   

   $$I_{p}=\int_{0}^{\infty}{\rm d}\chi\;\chi^{p}\;e^{-(\chi^{2}-... ...{\rm d}\gamma \;\gamma^{\frac{p-1}{2}}\;e^{-(\gamma-\xi)^{2}}.$$
   
   

   The form of the integral to be used depends on the value of $p$:

   $p=0$:

   
   

   $$I_{0}(\xi)=\int_{0}^{\infty}{\rm d}\chi\;e^{-(\chi^{2}-\xi)^... ... \mbox{    }\chi_{\rm max}=\sqrt{\frac{\xi+\vert\xi\vert}{2}}.$$

   
   

   $p=1$:

   
   

   $$I_{1}(\xi)=\frac{1}{2}\int_{0}^{\infty}{\rm d}\gamma\;e^{-(\... ...{2}}, \mbox{    }\gamma_{\rm max}=\frac{\xi+\vert\xi\vert}{2}.$$

   
   

   $p=2$:

   
   

   $$I_{2}(\xi)=\int_{0}^{\infty}{\rm d}\chi\;\chi^{2}\;e^{-(\chi... ...box{    }\chi_{\rm max}=\sqrt{\frac{\xi+\sqrt{\xi^{2}+2}}{2}}.$$

   
   

   $p\geq 3$:

   
   

   $$I_{p}(\xi)=\frac{1}{2}\int_{0}^{\infty}{\rm d}\gamma \;\gamm... ...mbox{    }\gamma_{\rm max}=\frac{\xi+\sqrt{\xi^{2}+(p-1)}}{2}.$$

   
   

   The integration is to be by cubic splines, in which the integration increment $\Delta I$ for a function $y(x)$ with derivative $y'(x)$, from a point $(x_{1},y_{1})$ to a point $(x_{2},y_{2})$, is given by

   
   

   $$\Delta I=\frac{1}{2}(y_{1}+y_{2})\Delta x-\frac{1}{12}\Delta y'(\Delta x)^{2},$$
   
   

   where $\Delta x=x_{2}-x_{1}$, $\Delta y'=y'_{2}-y'_{1}$, and so on. The functions involved and their derivatives are

   $$\begin{eqnarray*} f_{p}(\chi) & = & \chi^{p}\;e^{-(\chi^{2}-\xi)^{2}}, g_{p}(... ...rac{p-1}{2}+2\xi\gamma -2\gamma^{2}\right]e^{-(\gamma-\xi)^{2}}. \end{eqnarray*}$$
   
   
   
   

3. Calculate the corresponding value of $\alpha$ using the formula

   
   

   $$\alpha=-2\left(\xi\sqrt{\lambda}+d\right).$$
   
   

This completes the algorithm. One should also recall that for even $\mathfrak{N}$ the relevant integrals, with odd $p$, can be written in terms of the error function. In particular, for $\mathfrak{N}=2$ one can use

$$\begin{eqnarray*} I_{1}(\xi) & = & \frac{\sqrt{\pi}}{4}\;[1+\Phi(\xi)], I_{3}... ...ac{\sqrt{\pi}}{4}\;\xi\;[1+\Phi(\xi)]+\frac{1}{4}\;e^{-\xi^{2}}. \end{eqnarray*}$$

This can be used to solve the problem in this case, probably in a more efficient way. For larger even values of $\mathfrak{N}$ one has to derive the corresponding relations for each case in order to do this.