B. Mathematical background

B.1. Special functions

B.1.1. The Dirac delta function and related functions

\label{sec-math-delta}

We encounter the Dirac delta function in many situations in (astro)physics. As physicists, we heuristically define the Dirac delta function, which we will generally simply refer to as the delta function, as the function

\begin{equation}\label{eq-math-delta} \delta(x) = \begin{cases} +\infty\,,\quad &x = 0\,,\\ 0\,, &x\neq 0\,.\end{cases} \end{equation}

However, this is only a useful way to think about the delta function, not a mathematically rigorous way to define it. Mathematically, the delta funtion isn’t technically a function, but something that is rather described as a measure or a distribution. A measure is a technical, mathematical specification of what we denote as \(\mathrm{d} x\) in an integral and a loose way of defining the delta function as a measure is through the equation

\begin{equation}\label{eq-math-delta-asint} \int_{-\infty}^{+\infty}\mathrm{d}x\,\delta(x)\,f(x) = f(0)\,, \end{equation}

where the combination \(\mathrm{d}x\,\delta(x)\) is the Dirac measure (mathematicians would be horrified at this abuse of notation and insist on a form like \(\int_{-\infty}^{+\infty}\delta\{\mathrm{d}x\}\,f(x) = f(0)\); but it will do for us). Thus, the delta function only leads to physical quantities through integration over it.

While one is in general quite happy when the solution of a physics problem requires an integral over a delta function (because Equation \ref{eq-math-delta-asint} makes these integrals simple compared to just about any other integral you encounter), we will occasionally need to perform integrals that involve a delta function that do not quite have the simple form of Equation \eqref{eq-math-delta-asint} and are therefore trickier to obtain. One way to perform such integrals is to view the delta function as the limit of a family of well-defined functions that approach the delta function as some parameter goes to, say, zero. There are many families that approach a delta function in this way, but one particularly simple one is that of the rectangular function \(\Pi(x)\) defined as

\begin{equation}\label{eq-math-rect1} \Pi(x) = \begin{cases} 0\,,\quad & |x| > \frac{1}{2}\,,\\ \frac{1}{2}\,,\quad & |x| = \frac{1}{2}\,,\\ 1\,,\quad & |x| < \frac{1}{2}\,.\end{cases} \end{equation}

In terms of the rectangular function, the delta function can be defined as

\begin{equation} \delta(x) = \lim_{h\rightarrow 0} {1 \over h}\,\Pi\left({x\over h}\right)\,. \end{equation}

As an example of how we can use this limit approximation, we can perform the integral of a delta function multiplied by a function \(f(x)\) over half the real line rather than the full real line

\begin{align}\label{eq-math-deltainthalfline-1} \int_{0}^{+\infty}\mathrm{d}x\,\delta(x)\,f(x) & = \lim_{h\rightarrow 0}\int_{0}^{+\infty}\mathrm{d}x\,{1 \over h}\,\Pi\left({x\over h}\right)\,f(x)\,,\\ & = \lim_{h\rightarrow 0}{1 \over h}\int_{0}^{+\infty}\mathrm{d}x\,\,\Pi\left({x\over h}\right)\,\left[f(0) + f'(0)\,x+ \mathcal{O}(x^2)\right]\,,\\ & = \lim_{h\rightarrow 0}\left[ {f(0) \over 2} + h\,{f'(0) \over 4}+ \mathcal{O}(h^2)\right]\,,\\ & = {f(0) \over 2}\,.\label{eq-math-intdeltahalfline-2} \end{align}

This result makes intuitive sense, but the limit approach of the delta function makes it rigorous. We could use a similar approach to show that \(\int_{-a}^{+b}\mathrm{d}x\,\delta(x)\,f(x) = f(0)\) for any \(a,b > 0\).

A function that is related to the delta function is the Heaviside function \(\Theta(x)\) defined as

\begin{equation}\label{eq-math-heaviside} \Theta(x) = \begin{cases} 1\,,\quad & x > 0\,,\\ \frac{1}{2}\,, \quad & x = 0\,,\\ 0\,,\quad &x < 0\,.\end{cases} \end{equation}

The Heaviside function is also known as the step function. The Heaviside can be thought of as the integral of the delta function \(\Theta(x) = \int_{-\infty}^{x} \mathrm{d}s\,\delta(s)\), which makes sense if we define the delta function as the limit of rectangular functions (but for other limit definitions of the delta function, this equality does not necessarily hold at \(x=0\)). In terms of the Heaviside function, we can write the integral of the delta function over half the real line from Equation \eqref{eq-math-deltainthalfline-1} as

\begin{equation}\label{eq-math-intdeltaheaviside} \int_{0}^{+\infty}\mathrm{d}x\,\delta(x)\,f(x) = \int_{-\infty}^{+\infty}\mathrm{d}x\,\delta(x)\,\Theta(x)\,f(x) = {f(0)\over 2}\,. \end{equation}

Written in this way, it becomes clear that the integral \(\int_{-\infty}^{+\infty}\mathrm{d}r\,r\,\delta(r-R)\,\Theta(r-R)\) that we need in Chapter 15.1.3 is a special case of this integral and that \(\int_{-\infty}^{+\infty}\mathrm{d}r\,r\,\delta(r-R)\,\Theta(r-R) = R/2\).

B.1.2. Bessel functions

\label{sec-math-bessel}

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)
figsize(7,5)
plot(xs,special.jv(0,xs),label=r'$J_0$')
plot(xs,special.jv(1,xs),label=r'$J_1$')
plot(xs,special.jv(2,xs),label=r'$J_2$')
plot(xs,special.jv(3,xs),label=r'$J_3$')
xlabel(r'$x$')
ylabel(r'$J_m(x)$')
legend(frameon=False,fontsize=18.);

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}

and

\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 8.3.3.2.

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}

where

\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))\
        /(numpy.cosh(numpy.pi*numpy.sinh(t))+1)

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))\
                *special.jv(nu,numpy.pi/h*de_psi(h*xinuk))\
                *de_psiprime(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],
      label=r'$m=0$')
loglog(hs,[numpy.fabs(1.-de_approx(lambda x: 1.,1,h,10000)) for h in hs],
      label=r'$m=1$')
loglog(hs,[numpy.fabs(1.-de_approx(lambda x: 1.,2,h,10000)) for h in hs],
      label=r'$m=2$')
xlabel(r'$h$')
ylabel(r'$1-\int \mathrm{d}x\,J_m(x)$')
legend(frameon=False,fontsize=18.);

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

B.1.3. Elliptic integrals

\label{sec-math-ellipticint}

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}