B. Mathematical background

B.1. Special functions

B.1.2. Bessel functions


Bessel functions, named after Friedrich Bessel who studied them in detail (and also measured the first astronomical parallax!), are solutions \(y(x)\) of the differential equation

\begin{equation}\label{eq-math-bessel-diffeq} x^2\frac{\mathrm{d}^2 y}{\mathrm{d} x^2} + x\,\frac{\mathrm{d}y}{\mathrm{d}x} + \left(x^2-\alpha^2\right)\,y = 0\,. \end{equation}

They often come up in physics and in gravitational dynamics in particular, because they appear in the solution of the Helmholtz equation in cylindrical coordinates (the Poisson equation with the density proportional to the density which is used to find eigenfunctions of the Laplace operator).

Bessel functions come in various kinds. In this text, we encounter Bessel functions of the first kind \(J_m(x)\) with \(\alpha\) an integer \(m\), which are solutions that are finite at the origin \(x=0\). Bessel functions of the first kind are highly oscillatory, with an amplitude the falls off as \(\approx x^{-1/2}\) and a set of zero-crossings that does not display regular behavior. The following figure shows the first few Bessel functions of the first kind:

from scipy import special
xs= numpy.linspace(0.,10.,101)

Bessel functions of the first kind of different integer orders are related by relations such as

\begin{equation}\label{eq-math-bessel-negm-posm} J_{-m}(x) = (-1)^m\,J_m(x)\, \end{equation}


\begin{equation}\label{eq-math-bessel-dJmx} {\mathrm{d} \over \mathrm{d} x} \Big[ x^m\,J_m(x)\Big] = x^m\, J_{m-1}(x)\,, \end{equation}

Combining these, we can derive that

\begin{equation}\label{eq-math-bessel-dJ0J1} {\mathrm{d} J_0(x) \over \mathrm{d} x} = -J_1(x)\,. \end{equation}

Another useful relation that follows from Equation \eqref{eq-math-bessel-dJmx} is

\begin{equation}\label{eq-math-bessel-dJ1} {\mathrm{d} J_1(x) \over \mathrm{d} x} = J_0(x)-{J_1(x)\over x}\,. \end{equation}

Bessel functions of the first kind satisfy the following orthonormality relation (the Fourier-Bessel theorem):

\begin{equation}\label{eq-math-bessel-FourierBessel} \int_0^\infty\mathrm{d}R\,J_m(kR)\,J_m(k'R)R = \delta(k-k')/k,\quad \forall k, k' > 0\,, \end{equation}

where \(\delta(k-k')\) is the delta function. The Laplace transform of \(J_m(kx)\) is

\begin{equation}\label{eq-math-bessel-laplace} \mathcal{L}\left[J_m(kx)\right](s) = \int_0^\infty \,\mathrm{d} x\,e^{-st}\,J_m(kx) = {k^m \over \sqrt{s^2+k^2}\,\left(s + \sqrt{s^2+k^2}\right)^m}\,, \end{equation}

and in particular (set \(s=0\) and \(k=1\))

\begin{equation}\label{eq-math-bessel-norm} \int_0^\infty \mathrm{d} x\,J_m(x) = 1\,. \end{equation}

Bessel functions of the second kind are denoted as \(Y_\alpha(x)\) and they are, for integer \(\alpha\), the second linearly independent solution of the defining differential equation \eqref{eq-math-bessel-diffeq} (for non-integer alpha, \(J_m\) and \(J_{-m}\) are two linearly independent solutions, because Equation [\ref{eq-math-bessel-negm-posm}] does not hold for non-integer \(\alpha\)). Bessel functions of the second kind have a singularity at the origin.

Modified Bessel functions of the first, \(I_\alpha(x)\), and second, \(K_\alpha(x)\), kind are defined in terms of the Bessel function of the first kind evaluated at purely imaginary values (Bessel functions of the first kind are defined on the entire complex plane) as

\begin{align}\label{eq-math-bessel-modified} I_\alpha(x) & = i^{-\alpha}\,J_\alpha(ix)\,,\\ K_\alpha(x) & = {\pi \over 2}\,{I_{-\alpha}(x)-I_\alpha(x) \over \sin (\alpha \pi)}\,. \end{align}

For integer \(\alpha = m\), \(K_m\) is defined as the limit \(\alpha \rightarrow m\). Similar to how the oscillatory sine and cosine functions become a growing or decaying exponential when evaluated for imaginary arguments, the modified Bessel functions of the first and second kinds are exponentially growing and decaying functions, respectively, without any oscillatory behavior.

Combining their definition with relations for the Bessel functions of the first kind from Equation \eqref{eq-math-bessel-negm-posm} to \eqref{eq-math-bessel-dJ1}, one can show that

\begin{align}\label{eq-math-bessel-dI0} {\mathrm{d} I_0 (x) \over \mathrm{d} x} & = I_1(x)\,,\\ {\mathrm{d} K_0 (x) \over \mathrm{d} x} & = -K_1(x)\,,\label{eq-math-bessel-dK0}\\ {\mathrm{d} I_1 (x) \over \mathrm{d} x} & = I_0(x)-I_1(x)/x\,,\label{eq-math-bessel-dI1}\\ {\mathrm{d} K_1 (x) \over \mathrm{d} x} & = -K_0(x)-K_1(x)/x\,,\label{eq-math-bessel-dK1} \end{align}

which are all relations we need in deriving the rotation curve of a razor-thin exponential disk in Chapter

Spherical Bessel functions are the solutions of the differential equation

\begin{equation}\label{eq-math-spherebessel-diffeq} x^2\frac{\mathrm{d}^2 y}{\mathrm{d} x^2} + 2x\,\frac{\mathrm{d}y}{\mathrm{d}x} + \left(x^2-\alpha\,[\alpha+1]\right)\,y = 0\,. \end{equation}

This equation is similar to the defining equation \eqref{eq-math-bessel-diffeq} of the Bessel functions except for the factor of two in the second term and that we have replaced \(\alpha^2\) with \(\alpha\,(\alpha+1)\). Like the defining equation \eqref{eq-math-bessel-diffeq}, Equation \eqref{eq-math-spherebessel-diffeq} appears in solutions of the Helmholtz equation, but now in spherical coordinates. Equation \eqref{eq-math-spherebessel-diffeq} has two linearly-independent solutions \(j_n(x)\) and \(y_n(x)\) that are given in terms of the Bessel functions of the first and second kind through

\begin{align}\label{eq-math-bessel-jn} j_n(x) & = \sqrt{{\pi\over 2x}}\,J_{\alpha+1/2}(x)\,\\ y_n(x) & = \sqrt{{\pi\over 2x}}\,Y_{\alpha+1/2}(x)\,. \end{align}

Their oscillatory behavior makes integrating over Bessel functions of the first kind difficult. Yet we need such integrals in, for example, the general solution for the potential of a razor-thin disk or to compute the potential and forces of a double-exponential disk. A simple strategy for obtaining a high-precision integral over a Bessel function of the first kind is to numerically integrate between the zeros for a large number of zeros, making use of the fact that the amplitude of the oscillation drops off and that functions that return the zeros of the Bessel functions are readily available in many programming languages.

A different and better way is to use numerical quadrature methods specifically designed to deal with oscillatory functions. For Bessel functions of the first kind, the quadrature formula from Ogata (2005) is very useful. Ogata (2005) demonstrates that integrals of the type

\begin{equation} \int_0^\infty \mathrm{d}x\,J_\nu(x)\,f(x) \end{equation}

can be approximated as

\begin{align} \int_0^\infty & \mathrm{d}x\,J_\nu(x)\,f(x) \approx \label{eq-math-besselint}\\ & \pi\,\sum_{k=1}^{\infty}{w_{\nu k}\,f\left({\pi \over h}\,\psi\left[h\,\xi_{\nu k}\right]\right)\,J_\nu\left({\pi \over h}\,\psi\left[h\,\xi_{\nu k}\right]\right)\,\psi'\left(h\,\xi_{\nu k}\right)}\,,\nonumber \end{align}


\begin{equation} \psi(t) = t\,\tanh\left({\pi \over 2}\,\sinh t\right)\,, \end{equation}

and \(\psi'(t)\) is the derivative of \(\psi(t)\)

\begin{equation} \psi'(t) = {\sinh\left(\pi\,\sinh t\right)+\pi\,t\,\cosh t \over \cosh\left(\pi\,\sinh t\right)+1}\,. \end{equation}

The weights \(w_{\nu k}\) are given by

\begin{equation}\label{eq-math-besselint-weights} w_{\nu k} = {2\over \pi^2\,\xi_{\nu k}\,J^2_{\nu+1}\left(\pi\,\xi_{\nu k}\right)}\,. \end{equation}

The \(\xi_{\nu k}\) are the zeros of the Bessel function \(J_{\nu}(\pi\,x)\) and \(h\) is a numerical step size that should be set as small as possible. Because as \(k \rightarrow \infty\), \({\pi \over h}\,\psi(h\,\xi_{\nu k}) \rightarrow \pi\,\xi_{\nu k}\), which are the zeros of \(J_{\nu}(x)\), the terms in the sum quickly go to zero, such that convergence is obtained with a small number of terms. Mathematically, this approach happens “double exponentially” (that is, fast), but this has nothing to do with and should not be confused for the double exponential disk model common in galactic dynamics (however, this formula is highly useful for working with double exponential disk models!).

To illustrate this quadrature method for Bessel functions, we can apply it to some integrals with known solutions, for example, relation \eqref{eq-math-bessel-norm} above. First, we implement the \(\psi(t)\) and \(\psi'(t)\) functions:

def de_psi(t):
    return t*numpy.tanh(numpy.pi/2.*numpy.sinh(t))
def de_psiprime(t):
    return (numpy.sinh(numpy.pi*numpy.sinh(t))+numpy.pi*t*numpy.cosh(t))\

and then some functions that return the relevant zeros of the Bessel function and the weights from Equation \eqref{eq-math-besselint-weights}

from scipy import special
def de_xinuk(nu,nzeros):
    return special.jn_zeros(nu,nzeros)/numpy.pi
def de_wnuk(nu,nzeros,xinuk=None):
    if xinuk is None: xinuk= de_xinuk(nu,nzeros)
    return 2./(numpy.pi**2.*xinuk*special.jv(nu+1,numpy.pi*xinuk)**2.)

and finally a function that implements the quadrature formula of Equation \eqref{eq-math-besselint} (here, we stop adding terms once they have a relative contribution less than \(10^{-15}\))

def de_approx(fun,nu,h,nzeros):
    xinuk= de_xinuk(nu,nzeros)
    sumthis= fun(numpy.pi/h*de_psi(h*xinuk))\
    sumthis[numpy.fabs(sumthis/numpy.cumsum(sumthis)) < 1e-15]= 0.
    sumthis[numpy.isnan(sumthis)]= 0.
    return numpy.pi*numpy.sum(de_wnuk(nu,nzeros,xinuk=xinuk)*sumthis)

Applying this formula to Equation \eqref{eq-math-bessel-norm} for \(m=0,1,2\) for a range of step sizes \(h\), we get

hs= numpy.geomspace(1e-3,1,100)
loglog(hs,[numpy.fabs(1.-de_approx(lambda x: 1.,0,h,10000)) for h in hs],
loglog(hs,[numpy.fabs(1.-de_approx(lambda x: 1.,1,h,10000)) for h in hs],
loglog(hs,[numpy.fabs(1.-de_approx(lambda x: 1.,2,h,10000)) for h in hs],
ylabel(r'$1-\int \mathrm{d}x\,J_m(x)$')

It is clear that we reach very high precision when using a suitably small step \(h\).

B.1.3. Elliptic integrals


Elliptic integrals are a set of integrals that originally appeared in the problem of finding the arc length of an ellipse, but that more generally pop up in the solution of various mathematical-physics problems. They are notorious for being defined in slightly different ways by different authors and different software packages. We use the following definitions.

The incomplete elliptic integral of the first kind, \(K(\phi,m)\) is defined as the integral

\begin{equation} K(\phi,m) = \int_0^\phi \mathrm{d}\theta\,{1 \over \sqrt{1-m\,\sin^2 \theta}}\,. \end{equation}

The complete elliptic integral of the first kind is \(K(m) \equiv K(\pi/2,m)\) or

\begin{equation}\label{eq-math-ellipint-Km} K(m) = \int_0^{\pi/2} \mathrm{d}\theta\,{1 \over \sqrt{1-m\,\sin^2 \theta}}\,. \end{equation}

In this book, we occasionally need the value of \(K(0)\) for which the integrand is one and therefore \(K(0) = \pi/2\).

Elliptic integrals of the second kind are defined similarly, with the incomplete elliptic integral of the second kind \(E(\phi,m)\)

\begin{equation} E(\phi,m) = \int_0^\phi \mathrm{d}\theta\,\sqrt{1-m\,\sin^2 \theta}\,. \end{equation}

and the complete elliptic integral of the second kind is \(E(m) \equiv E(\pi/2,m)\) or

\begin{equation} E(m) = \int_0^{\pi/2} \mathrm{d}\theta\,\sqrt{1-m\,\sin^2 \theta}\,. \end{equation}

There are elliptic integrals of the third kind as well, but we do not use these.

The derivatives of the complete elliptical integrals of the first and second kind are

\begin{align}\label{eq-math-elliptint-dK} {\mathrm{d} K(m) \over \mathrm{d} m} & = {E(m) - (1-m)\,K(m) \over 2\,(1-m)\,m}\,,\\ {\mathrm{d} E(m) \over \mathrm{d} m} & = {E(m)-K(m) \over 2\,m}\,.\label{eq-math-elliptint-dE} \end{align}