14. Bessel Functions

Hide code cell content

 
import mmf_setup; mmf_setup.nbinit()
import os
from pathlib import Path
FIG_DIR = Path(mmf_setup.ROOT) / '../Docs/_build/figures/'
os.makedirs(FIG_DIR, exist_ok=True)
import logging; logging.getLogger("matplotlib").setLevel(logging.CRITICAL)
%matplotlib inline
import numpy as np, matplotlib.pyplot as plt
try: from myst_nb import glue
except: glue = None

This cell adds /home/docs/checkouts/readthedocs.org/user_builds/physics-571-math-methods/checkouts/latest/src to your path, and contains some definitions for equations and some CSS for styling the notebook. If things look a bit strange, please try the following:

  • Choose "Trust Notebook" from the "File" menu.
  • Re-execute this cell.
  • Reload the notebook.

14. Bessel Functions#

Spherical Bessel Functions#

The spherical Bessel functions \(j_{\ell}(z)\) and \(y_{\ell}(z)\) arise when solving the Helmholtz equation in 3D:

\[\begin{gather*} (\nabla^2 - k^2)R_{\ell}(r)Y^{\ell}_{m}(\theta,\phi) = 0 \end{gather*}\]

where \(Y^{\ell}_{m}(\theta,\phi)\) are the spherical harmonics. Letting \(z=kr\), we have the solution \(R_{\ell}(r) = f_{\ell}(z)\)

\[\begin{gather*} z^2 f_{\ell}''(z) + 2zf_{\ell}'(z) + \bigl[z^2 - \ell(\ell+1)\bigr]f_{\ell}(z) = 0. \end{gather*}\]

Do It! Derive this from Laplacian in spherical coordinates.

Review Laplacian and derive the Laplacian in spherical coordinates \(X^{\alpha} = (r, \theta, \phi)\). Hint: justify and use the line element

\[\begin{gather*} \d{s}^2 = \d{x}^2 + \d{y}^2 + \d{z}^2 = \d{r}^2 + r^2 \d{\theta}^2 + r^2\sin^2\theta \d{\phi}^2 \end{gather*}\]

to determine the components of the metric. Then use the note in 5. Angular Momentum of Assignment 3: Vector Calculus to argue that, with the spherical harmonics, the angular parts just give the term \(\ell(\ell+1)\):

\[\begin{gather*} -\nabla^2 = \frac{L^2}{r^2} - \frac{1}{r^2}\pdiff{}{r}\left(r^2\pdiff{}{r}\right), \qquad L^2Y^{m}_{l}(\theta, \phi) = \ell(\ell+1)Y^{m}_{l}(\theta, \phi). \end{gather*}\]

Relationship to Bessel Functions#

The solutions \(f_{\ell} \in \{j_{\ell}, y_{\ell}\}\) can be expressed in terms of the standard Bessel functions \(F_{n} \in \{J_{n}, Y_{n}\}\) at half-integer values of the index \(n=\ell +1/2\):

\[\begin{align*} f_{\ell}(z) &= \sqrt{\frac{\pi}{2z}}F_{\ell+1/2}(z), & & \text{i.e.} & j_{n}(z) &= \sqrt{\frac{\pi}{2z}}J_{n+1/2}(z), & y_{n}(z) &= \sqrt{\frac{\pi}{2z}}Y_{n+1/2}(z). \end{align*}\]

Do It! Show that if \(F_{n}(z)\) satisfies Bessel’s equation, \(f_{\ell}(z)\) solves the ODE above.

Start from Bessel’s equation

\[\begin{gather*} z^2 F''_{n} + z F_{n}' + (z^2 - n^2) F_{n} = 0. \end{gather*}\]

Series Representation#

Using the Frobenius method, one can find a series solution:

\[\begin{align*} (-1)^{\ell}y_{-\ell-1}(z) = j_{\ell}(z) &= \frac{z^{\ell}}{(2\ell+1)!!} \left( 1 - \frac{\tfrac{1}{2}z^2}{1!(2\ell+3)} + \frac{(\tfrac{1}{2}z^2)^2}{2!(2\ell+3)(2\ell+5)} - \dots \right),\\ (-1)^{\ell+1}j_{-\ell-1}(z) = y_{\ell}(z) &= \frac{(2\ell-1)!!}{-z^{\ell+1}} \left( 1 - \frac{\tfrac{1}{2}z^2}{1!(1-2\ell)} + \frac{(\tfrac{1}{2}z^2)^2}{2!(1-2\ell)(3-2\ell)} - \dots \right). \end{align*}\]

By inspection, we can affirm the following

\[\begin{align*} j_{0}(z) &= \frac{\sin z}{z}, & j_{1}(z) &= \frac{\sin{z} - z \cos{z}}{z^2}, & j_{2}(z) &= \frac{(3-z^2)\sin{z} - 3 z \cos{z}}{z^3}, \\ y_{0}(z) &= \frac{\cos{z}}{-z}, & y_{1}(z) &= \frac{\cos{z} + z \sin{z}}{-z^{2}}, & y_{2}(z) & =\frac{(3-z^{2}) \cos{z} + 3 z \sin{z}}{-z^{3}}. \end{align*}\]

Note that

Hide code cell content

 
from sympy import *
from IPython.display import Latex

z = var('z')

def get_j(l, f0=sin(z)/z):
    res = f0
    for n in range(l):
        res = -res.diff(z)/z
    return (z**l * res).simplify()

def get_y(l):
    return get_j(l, f0=cos(z)/(-z))

N = 3
display(Latex(r', \qquad '.join([
    fr'j_{{{l}}}(z)={latex(get_j(l))}'
    for l in range(N)])))
display(Latex(r', \qquad '.join([
    fr'y_{{{l}}}(z)={latex(get_y(l))}'
    for l in range(N)])))
\[j_{0}(z)=\frac{\sin{\left(z \right)}}{z}, \qquad j_{1}(z)=\frac{- z \cos{\left(z \right)} + \sin{\left(z \right)}}{z^{2}}, \qquad j_{2}(z)=\frac{- z^{2} \sin{\left(z \right)} - 3 z \cos{\left(z \right)} + 3 \sin{\left(z \right)}}{z^{3}}\]
\[y_{0}(z)=- \frac{\cos{\left(z \right)}}{z}, \qquad y_{1}(z)=- \frac{z \sin{\left(z \right)} + \cos{\left(z \right)}}{z^{2}}, \qquad y_{2}(z)=\frac{z^{2} \cos{\left(z \right)} - 3 z \sin{\left(z \right)} - 3 \cos{\left(z \right)}}{z^{3}}\]

Do It! Derive these solutions.

Using the Frobenius method, find the indicial equation and solve it for the. You should get the series for \(j_{\ell}\) for \(s = l \geq 0\) and for \(y_{\ell}\) for negative \(s = -\ell-1 < 0\). Note: this method will not fix the overall coefficients which are arbitrary. The normalization listed above is conventional.

Partial result:

\[\begin{gather*} f_{s}(z) = \sum_{j=0}^{\infty} a_{j} z^{s+j}, \qquad a_{1} = 0, \qquad a_{2m} = \frac{-a_{2m-2}}{2m(2s+2m+1)},\\ a_{2} = \frac{-a_{0}}{2(2s+3)}, \qquad a_{4} = \frac{a_{0}}{2^{2} 2! (2s+3)(2s+5)}, %\qquad a_{6} = \frac{-a_{0}}{(2^{3} 3! (2s+3)(2s+5)(2s+7)}, \\ f_{s}(z) = a_0z^{s}\left( 1 - \frac{z^2}{2(2s+3)} + \frac{z^4}{2^{2} 2! (2s+3)(2s+5)} - \cdots \right)\\ = a_0z^{s}\left( 1 - \frac{\tfrac{1}{2}z^2}{1!(2s+3)} + \frac{(\tfrac{1}{2}z^2)^2}{2!(2s+3)(2s+5)} - \cdots \right). \end{gather*}\]

Note that the coefficient \(a_{0}(s)\) is just a normalization and differs from each series. To get \(j_{l}(z)\) we set \(s=l\), and to get \(y_{l}(z)\), we set \(s=-\ell-1\):

\[\begin{align*} j_{\ell}(z) &\propto z^{\ell}\left( 1 - \frac{\tfrac{1}{2}z^2}{1!(2\ell+3)} + \frac{(\tfrac{1}{2}z^2)^2}{2!(2\ell+3)(2\ell+5)} - \cdots \right),\\ y_{\ell}(z) &\propto \frac{1}{z^{1+\ell}}\left( 1 - \frac{\tfrac{1}{2}z^2}{1!(1-2\ell)} + \frac{(\tfrac{1}{2}z^2)^2}{2!(1-2\ell)(3-2\ell)} - \cdots \right). \end{align*}\]

The ODE above can be re-expressed so they can be analyzed with Sturm-Liouville Theory:

\[\begin{gather*} \overbrace{\left(\diff[2]{}{z} + 2z\diff{}{z} + z^2\right)}^{\mathcal{L}(z)}f_{\ell}(z) = \overbrace{\ell(\ell+1)}^{E_{\ell}}f_{\ell}(z). \end{gather*}\]

Do It! Check that \(\mathcal{L}\) is self-adjoint.

We start by checking whether or not \(\mathcal{L}(z)\) is Hermitian (here we check that it is symmetric since our functions are real):

\[\begin{gather*} \braket{f|\mathcal{L}(z)|g} \stackrel{?}{=} \braket{g|\mathcal{L}(z)|f}\\ \int_0^{\infty}(z^2 fg'' + 2zfg' + z^2fg)\d{z} \stackrel{?}{=} \int_0^{\infty}(z^2gf'' + 2zgf' + z^2gf)\d{z}\\ % \int_0^{\infty}(-(z^2 f)'g' + 2zfg')\d{z} \stackrel{?}{=} \int_0^{\infty}(-(z^2g)'f' + 2zgf')\d{z},\\ \int_0^{\infty}(-z^2f'g')\d{z} \stackrel{!}{=} \int_0^{\infty}(-z^2g'f')\d{z}. \end{gather*}\]

Hence, the equation is self-adjoint as long as the boundary terms vanish.

Keeping only the boundary terms, we have

\[\begin{gather*} \left. z^2 fg' \right|_{z=0}^{\infty} \stackrel{?}{=} \left. z^2 f'g \right|_{z=0}^{\infty}. \end{gather*}\]

This is not obviously true, and it will impact the orthogonality conditions. One way to ensure that this condition is satisfied is to choose boundaries that are roots of the Bessel functions. This is the basis of the somewhat complicated formula (14.183) in [Arfken et al., 2013].

Since \(\mathcal{L}\) is Hermitian, we might expect an orthogonality condition of the form

\[\begin{gather*} \braket{f_{n}|f_{m}} = \int_0^{\infty} f_{n}(z)f_{m}(z)\d{z} \propto \delta_{mn} \end{gather*}\]

however, this also requires self-adjointness which is spoiled by the boundary conditions. It turns out that, since \(j_{\ell}(z)\) are regular at \(z=0\), we can extend the region of integration to \(z\in (-\infty, \infty)\) and orthogonality holds here:

\[\begin{gather*} \int_{-\infty}^{\infty} j_{m}(z)j_{n}(z) \d{z} = \frac{\pi \delta_{mn}}{2m+1}. \end{gather*}\]

If we restrict \(z\in [0, \infty)\), then we have

\[\begin{gather*} \braket{j_{m}|j_{n}} = \frac{\pi}{2(m+n+1)} \frac{\sin\Bigl(\tfrac{\pi}{2}(m-n)\Bigr)}{\tfrac{\pi}{2}(m-n)}. \end{gather*}\]

Note the orthogonality for \(m-n\) even: this means that the even functions \(j_{2n}\) are orthogonal, as are the odd functions \(j_{2n+1}\), but the even functions are not orthogonal to the odd functions.

Recurrence Relations#

For \(f_{\ell} \in \{j_{\ell}, y_{\ell}\}\) we have the following recurrence relations:

(1)#\[\begin{align} f_{\ell-1}(z) + f_{\ell+1}(z) &= (2\ell + 1)\frac{f_{\ell}(z)}{z}, \\ \ell f_{\ell-1}(z) - (\ell +1)f_{\ell+1}(z) &= (2\ell + 1)f'_{\ell}(z). \end{align}\]

Do It! Show these from the series expansion.

This is a bit messy, but you should be able to demonstrate that these are correct from the series representations.

(INCOMPLETE) Here is some preliminary work.

Note that for \(j_{\ell}\), \(\ell \rightarrow \ell \pm 1\) is equivalent to taking \(s \rightarrow s \pm 1\), while for \(y_{\ell}\), \(\ell \rightarrow \ell \pm 1\) is equivalent to taking \(s \rightarrow s \mp 1\) or \(\abs{s} \rightarrow \abs{s} \pm 1\).

From the series solution above, including the normalization, let

\[\begin{gather*} w_{s}(z) = (2s + 1)!!j_{s}(z) = \frac{-1}{(1-2s)!!}y_{-s-1}(z)\\ = z^{s}\left(1 - \frac{z^2}{2^1 1!(2s+3)} + \frac{z^{4}}{2^2 2!(2s+3)(2s+5)} - \dots \right),\\ w_{s+ 1} = z^{s + 1}\left(1 - \frac{z^2}{2^1 1!(2s+5)} + \frac{z^{4}}{2^2 2!(2s+5)(2s+7)} - \dots \right),\\ w_{s - 1} = z^{s-1}\left(1 - \frac{z^2}{2^1 1!(2s+1)} + \frac{z^{4}}{2^2 2!(2s+1)(2s+3)} - \dots \right). \end{gather*}\]
\[\begin{gather*} j_{l\pm 1} = \frac{w_{s\pm 1}}{(2s +1 \pm 2)!!}\\ y_{l\pm 1} = \frac{w_{s\pm 1}}{(2s +1 \pm 2)!!} \end{gather*}\]

Massaging these, we can recover the original form of the series:

\[\begin{gather*} w_{s + 1} = z^{s + 1}\left(1 - \frac{z^2}{2^1 1!(2s+5)} + \frac{z^{4}}{2^2 2!(2s+5)(2s+7)} - \dots \right),\\ w_{s - 1} = z^{s-1}\left(1 - \frac{z^2}{2^1 1!(2s+1)} + \frac{z^{4}}{2^2 2!(2s+1)(2s+3)} - \dots \right). \end{gather*}\]
\[\begin{gather*} f_{s}(z) = z^{s}\sum_{j=0}^{\infty} a_{2j}(s) z^{2j}, \qquad a_{2j}(s) = \frac{-a_{2j-2}(s)}{2j(2s+2j+1)},\\ f_{s\pm 1} = z^{s\pm 1}\sum_{j=0}^{\infty} a_{2j}(s\pm 1) z^{2j},\\ a_{2j}(s\pm 1) = \frac{-a_{2j-2}(s\pm 1)}{2j(2(s\pm 1)+2j+1)},\\ = \frac{-a_{2j-2}(s\pm 1)}{2j(2(s\pm 1)+2j+1)},\\ \end{gather*}\]

These relationships imply the recurrence

\[\begin{gather*} \frac{f_{\ell+1}}{z^{\ell+1}} = -\frac{1}{z}\diff{}{z}\frac{f_{\ell}}{z^{\ell}}. \end{gather*}\]

Adding the normalization

\[\begin{align*} j_{\ell}(z) &\propto z^{\ell}\left( 1 - \frac{\tfrac{1}{2}z^2}{1!(2\ell+3)} + \frac{(\tfrac{1}{2}z^2)^2}{2!(2\ell+3)(2\ell+5)} - \cdots \right),\\ y_{\ell}(z) &\propto \frac{1}{z^{1+\ell}}\left( 1 - \frac{\tfrac{1}{2}z^2}{1!(1-2\ell)} + \frac{(\tfrac{1}{2}z^2)^2}{2!(1-2\ell)(3-2\ell)} - \cdots \right). \end{align*}\]

Note that

\[\begin{gather*} f_{s}'(z) = a_0z^{s-1}\left( 0 - \frac{(s+2)z^{2}}{2^{1}1!(2s+3)} + \frac{(s+4)z^{4}}{2^{2}2!(2s+3)(2s+5)} - \cdots \right). \end{gather*}\]

Here we have used the properties

\[\begin{gather*} a_{2m} = \frac{(-1)^{m}a_0 (2s-1)!!}{(2m)!!(2s+2m+1)!!} = \frac{(-1)^{m}a_0 (2s-1)!!}{(2m)!!(2s+2m+1)!!} = 2\frac{(-1)^{m}a_0 (2s)! (m+s+1)!}{m!(2m+2s+2)! s!}\\ (2m)!! = 2^{m}m!, \qquad (2m-1)!! = \frac{(2m)!}{(2m)!!} = \frac{(2m)!}{2^{m}m!}, \\ (2m+1)!! = \frac{(2m+2)!}{(2m+2)!!} = \frac{(2m+2)!}{2^{m+1}(m+1)!} \end{gather*}\]
\[\begin{gather*} a_{2m} = \frac{(-1)^{m}a_0}{(2m)!!(2s+2m+1)(2s+2(m-1)+1)\dots(2s+1)} . \end{gather*}\]
\[\begin{align*} j_{s}(z) &= \frac{z^{s}}{(2s+1)!!} \left( 1 - \frac{\tfrac{1}{2}z^2}{1!(2s+3)} + \frac{(\tfrac{1}{2}z^2)^2}{2!(2s+3)(2s+5)} - \dots \right),\\ y_{\ell}(z) &= \frac{(-3-2s)!!z^{s}}{-1} \left( 1 - \frac{\tfrac{1}{2}z^2}{1!(2s+3)} + \frac{(\tfrac{1}{2}z^2)^2}{2!(2s+3)(2s+5)} - \dots \right). \end{align*}\]

Rayleigh’s Formula#

For \(f_{\ell} \in \{j_{\ell}, y_{\ell}\}\) we have:

\[\begin{gather*} f_{\ell}(z) = z^{\ell}\left(-\frac{1}{z}\diff{}{z}\right)^{\ell}f_0(z). \end{gather*}\]

Specifically:

\[\begin{align*} j_{\ell}(z) &= z^{\ell}\left(-\frac{1}{z}\diff{}{z}\right)^{\ell}\frac{\sin z}{z}, & y_{\ell}(z) &= z^{\ell}\left(-\frac{1}{z}\diff{}{z}\right)^{\ell}\frac{\cos z}{-z}. \end{align*}\]

These are very useful for symbolically constructing the functions.

Do It! Prove this using induction.

First we define the operator

\[\begin{gather*} \op{A} = -\frac{1}{z}\diff{}{z}, \qquad f_{n} = z^{n}\op{A}^nf_{0}. \end{gather*}\]

Now we can express this a a recurrence relation

\[\begin{gather*} f_{n+1} = z^{n+1}\op{A}^{n+1}f_0 = z^{n+1}\op{A}z^{-n}\underbrace{z^{n}\op{A}^{n}f_0}_{f_{n}} = z^{n+1}\op{A}z^{-n}f_{n}. \end{gather*}\]

We could just act with \(\op{A}\), but it will be convenient to first work the commutation relation:

\[\begin{gather*} [\op{A}, z^{n}] = -nz^{n-2}. \end{gather*}\]

Proof:

\[\begin{gather*} [\op{A}, z^{n}]f = \op{A}z^{n}f - z^{n}\op{A}f = \frac{-1}{z}(nz^{n-1}f + z^{n}f') - z^{n}\op{A}f\\ = -nz^{n-2}f + \underbrace{z^{n}\underbrace{\frac{-1}{z}f'}_{\op{A}f} - z^{n}\op{A}f}_{0}. \end{gather*}\]

Thus, Rayleigh’s formula is equivalent to the recursion

\[\begin{align*} f_{n+1} &= z^{n+1}\bigl([\op{A}, z^{-n}] + z^{-n}\op{A}\bigr)f_{n} = \frac{n}{z}f_{n} + z\op{A}f_{n}\\ &= \frac{n}{z}f_{n} - f'_{n}. \end{align*}\]

To derive the inductive step, apply \(\op{L}\) to \(f_{n+1}\):

\[\begin{gather*} \mathcal{L}f_{n+1} = z^2\left(\frac{n}{z}f_{n} - f_{n}'\right)'' + 2z\left(\frac{n}{z}f_{n} - f_{n}'\right)' + z^2\left(\frac{n}{z}f_{n} - f_{n}'\right) \\ = z^2\left(-\frac{n}{z^2}f_{n} + \frac{n}{z}f_{n}' - f_{n}''\right)' + 2z\left(\frac{n}{z}f_{n} - f_{n}'\right)' + z^2\left(\frac{n}{z}f_{n} - f_{n}'\right). \end{gather*}\]

We are looking for occurrences of \(\mathcal{L}f_{n} = z^2f_{n}'' + 2zf'_{n} + z^2f_{n}\), hence the \(f'''_{n}\) looks potentially troublesome.

\[\begin{gather*} \mathcal{L}f_{n+1} = (n + 2)\Bigl( z(f''_{n} + f_{n}) + (1-n)f'_{n}\Bigr). \end{gather*}\]

This suggests using the ODE to remove the \(f_{n}''\) term:

\[\begin{gather*} f_{n}'' = \frac{n(n+1) - z^2}{z^2}f_{n} - \frac{2}{z}f'_{n}. \end{gather*}\]

This gives

\[\begin{gather*} \mathcal{L}f_{n+1} = z^2\left(\frac{n^2 - z^2}{z^2}f_{n} + \frac{3n}{z}f_{n}'\right)' + 2z\left(\frac{n}{z}f_{n} - f_{n}'\right)' + z^2\left(\frac{n}{z}f_{n} - f_{n}'\right)\\ = \frac{3n-2z^2}{z}f''_{n} + \frac{2nz^2 -2z^4 - 3n}{z^2}f'_{n} + \frac{nz^3-2nz-3n^2z}{z^2}f_{n}. \end{gather*}\]

To get this, we can compute

\[\begin{gather*} (\mathcal{L}f_{n})' = z^2f_{n}'' + 2zf'_{n} + z^2f_{n}, \end{gather*}\]
\[\begin{gather*} = z^2\left(\frac{2n}{z^3}f_{n} - \frac{2n}{z^2}f_{n}' + \frac{n}{z}f_{n}'' - f_{n}'''\right) + 2z\left(-\frac{n}{z^2}f_{n} + \frac{n}{z}f_{n}' - f_{n}''\right) + z^2\left(\frac{n}{z}f_{n} - f_{n}'\right)\\ = -nzf_{n}'' - z^2f_{n}''' + z^2\left(\frac{n}{z}f_{n} - f_{n}'\right) \end{gather*}\]

Now we manipulate this to find $\mathcal{L}f_{n} = z^2f_{n}‘’ + 2zf’

\[\begin{gather*} = -\Bigl(nzf_{n}' + z^2f_{n}''\Bigr)' + nf_{n}' + 2zf_{n}'' \end{gather*}\]

Now suppose that \(f_{n}\) satisfies the ODE:

\[\begin{gather*} z^2f_{n}'' + 2zf'_{n} + \bigl[z^2-n(n+1)\bigr]f_{n} = 0 \end{gather*}\]

Do It! Check that these are consistent with the series representation.

This is somewhat messy, but fairly easy to check. I would like a way to derive these generally, but have not yet found a good approach.

(INCOMPLETE) Here is some preliminary work.

Suppose that \(f(z)\) satisfies the ODE. We can look for functions \(a(z)\) and \(b(z)\) such that

\[\begin{gather*} g(z) = a(z)\diff{b(z)f(z)}{z} = a(z)b'(z)f(z) + a(z)b(z)f'(z) = ab'f + abf' \end{gather*}\]

also satisfies the ODE.

\[\begin{gather*} g = ab'f + abf',\\ p_2 g'' + p_1 g' + p_0 g = E_g g\\ p_2 f'' + p_1 f' + p_0 f = E_f f\\ \end{gather*}\]

This will include terms

\[\begin{gather*} g' = Af + Bf' + abf'', \qquad g'' = A'f + [A+B']f' + [B+(ab)']f'' + abf''',\\ \mathcal{L}f = z^2 f'' + 2z f' + z^2 f = Ef. \end{gather*}\]

Differentiating this gives:

\[\begin{gather*} z^2 f'' = - 2z f' + (E-z^2)f, \\ z^2 f''' = - 4z f'' + (E - 2 - z^2)f' - 2zf \end{gather*}\]

We now substiute this into

\[\begin{gather*} \mathcal{L}g = z^2 (A'f + [A+B']f' + [B+(ab)']f'' + abf''') + 2z (Af + Bf' + abf'') + z^2 (ab'f + abf') = \tilde{E}(ab'f + abf')\\ z^2 (A'f + [A+B']f') + (z^2[B+(ab)'] - 2abz)f'' + ab(E-2-z^2)f' -2abzf + 2z (Af + Bf') + z^2 (ab'f + abf') = \tilde{E}(ab'f + abf')\\ \end{gather*}\]

To do this, note that \(g''\) will contain a term \(abf'''\). We can get this from the ODE:

\[\begin{gather*} g' = (ab')'f + (2ab' + a'b)f' + abf'', \qquad g'' = (ab')''f + [3(ab')' + (a'b)']f' + [(2ab' + a'b) + (ab)'] f'' + ab f''', \qquad \end{gather*}\]

also satisfies the ODE. This requires

\[\begin{gather*} \mathcal{L}g \mathcal{L}g = z^2 g'' \end{gather*}\]

, which we found about to require that \(z^2f'g\) vanish at both \(z=0\) and \(z=\infty\).

Hide code cell content

 
from sympy import *

z, E_f, E_g = var('z, E_f, E_g')
f, a, b, p_0, p_1, p_2 = [Function(f)(z) for f in ['f', 'a', 'b', 'p_0', 'p_1', 'p_2']]
g = a*(b*f).diff(z)
G = p_2 * g.diff(z,z) + p_1 * g.diff(z) + (p_0 - E_g) * g
F = p_2 * f.diff(z,z) + p_1 * f.diff(z) + (p_0 - E_f) * f

d2f = ((E-z**2)*f - 2*z*f)/z**2
d3f = d2f.diff(z).subs([(f.diff(z, z), d2f)])
res = z**2 * g.diff(z,z) + 2*z*g.diff(z) + z**2*g
res = res.subs([(f.diff(z,z,z), d3f), (f.diff(z, z), d2f)]).simplify()
f0, f1 = var('f0, f1')
res = res.subs([(f.diff(z), f1), (f, f0)])
assert(res.subs([(f0, 0), (f1, 0)]) == 0)
display(res.subs([(f0, 0), (f1, 1)]).simplify())
display(res.subs([(f0, 1), (f1, 0)]).simplify())
\[\displaystyle 3 z^{2} a{\left(z \right)} \frac{d^{2}}{d z^{2}} b{\left(z \right)} + z^{2} b{\left(z \right)} \frac{d^{2}}{d z^{2}} a{\left(z \right)} + 4 z^{2} \frac{d}{d z} a{\left(z \right)} \frac{d}{d z} b{\left(z \right)} - 2 z a{\left(z \right)} b{\left(z \right)} + 4 z a{\left(z \right)} \frac{d}{d z} b{\left(z \right)} + 2 z b{\left(z \right)} \frac{d}{d z} a{\left(z \right)} + e a{\left(z \right)} b{\left(z \right)}\]
\[\displaystyle \frac{z^{3} a{\left(z \right)} \frac{d}{d z} b{\left(z \right)} + z^{3} \frac{d^{2}}{d z^{2}} a{\left(z \right)} \frac{d}{d z} b{\left(z \right)} + 2 z^{2} \frac{d}{d z} a{\left(z \right)} \frac{d}{d z} b{\left(z \right)} + 2 z \left(z^{2} \frac{d^{2}}{d z^{2}} b{\left(z \right)} - \left(z^{2} + 2 z - e\right) b{\left(z \right)}\right) \frac{d}{d z} a{\left(z \right)} + 2 \left(z^{2} \frac{d^{2}}{d z^{2}} b{\left(z \right)} - \left(z^{2} + 2 z - e\right) b{\left(z \right)}\right) a{\left(z \right)} + \left(z^{3} \frac{d^{3}}{d z^{3}} b{\left(z \right)} - 3 z \left(z^{2} + 2 z - e\right) \frac{d}{d z} b{\left(z \right)} + 2 \left(z^{2} - z \left(z + 1\right) + 2 z - e\right) b{\left(z \right)}\right) a{\left(z \right)}}{z}\]

Generating Functions#

From [Olver et al., 2010] we have the following generating functions:

\[\begin{gather*} \frac{\cos\sqrt{z^2 - 2zt}}{z} = \sum_{n=0}^{\infty} \frac{t^n}{n!}j_{n-1}(z) = \frac{\cos{z}}{z} + \sum_{n=1}^{\infty} \frac{t^n}{n!}j_{n-1}(z), \\ \frac{\sin\sqrt{z^2 - 2zt}}{z} = \sum_{n=0}^{\infty} \frac{t^n}{n!}y_{n-1}(z) = \frac{\sin{z}}{z} + \sum_{n=1}^{\infty} \frac{t^n}{n!}y_{n-1}(z). \end{gather*}\]

(I have not checked these yet)

Integral Representations#

Hide code cell source

 
from scipy.integrate import quad

def get_t(s):
    """Return (t(s), t'(s)) for the Schlaefli contour for s in [-oo, oo]."""
    y = s*np.exp(-s**2/2)
    dy_ds = (1-s**2)*np.exp(-s**2/2)
    x = 1-s**2/2
    dx_ds = -s
    t = x + 1j*y
    dt_ds = dx_ds + 1j*dy_ds
    return t, dt_ds

fig, ax = plt.subplots(figsize=(4, 3))
s = np.linspace(-4, 4)
t, dt_ds = get_t(s)
ax.plot(t.real, t.imag)
s0 = -2
t0, dt0_ds = get_t(s0)
dt0 = 0.1*dt0_ds
plt.plot([-4, 0], [0, 0], '--r')
plt.plot([0], [0], 'xr')
ax.arrow(t0.real, t0.imag, dt0.real, dt0.imag, width=0.05)
ax.set(xlabel=r"$\Re t$", ylabel=r"$\Im t$", 
       xlim=(-2,1), aspect=1);
_images/e4ff398091aba9bc356a58562536cd202595d008fe03058f9066eb6f86573a39.png

The spherical Bessel function can be defined in terms of a Schläfli integral with the contour shown in the margin:

\[\begin{gather*} j_{n}(z) = \sqrt{\frac{\pi}{2z}}\oint \frac{e^{(z/2)(t-1/t)}}{t^{n+3/2}}\frac{\d{t}}{2\pi\I}. \end{gather*}\]

This can be used to derive the following:

\[\begin{gather*} j_{n}(z) = \frac{z^{n}}{2^{n+1}n!}\int_0^{\pi}\cos(z\cos\theta)(\sin{\theta})^{2n+1}\d{\theta}. \end{gather*}\]

Do It! Check this!

(INCOMPLETE)

Let \(t = R e^{\I\theta}\), \(\d{t} = \I t \d{\theta}\):

\[\begin{gather*} j_{n}(z) = \sqrt{\frac{\pi}{2z}}\oint_C \frac{e^{(z/2)(t-1/t)}}{t^{n+3/2}}\frac{\d{t}}{2\pi\I} \\ = \sqrt{\frac{\pi}{2z}}\int_{-\pi}^{\pi} \frac{e^{\I z R (e^{\I\theta}-e^{-\I\theta})/2\I}}{R^{n+3/2}e^{\I\theta(n+3/2)}} \frac{\d{\theta}}{2\pi}\\ = \sqrt{\frac{\pi}{2z}}\int_{-\pi}^{\pi} \frac{e^{\I z R \sin{\theta}}}{R^{n+3/2}e^{\I\theta n}e^{3\I\theta/2}} \frac{\d{\theta}}{2\pi} \end{gather*}\]

Hide code cell source

 
from functools import partial
from scipy.integrate import quad
from scipy.special import factorial

n = 5
j = get_j(n)
z_ = 1.2

def integrand1(th, z, n):
    return np.cos(z*np.cos(th))*(np.sin(th))**(2*n+1)*z**n/2**(n+1)/factorial(n)

assert np.allclose(j.subs([(z, z_)]),  quad(lambda th: integrand1(th, z_, n), 0, np.pi)[0])

def integrand2(s, z, n):
    t, dt_ds = get_t(s)
    J = np.exp(z/2 * (t-1/t)) / t**(n+3/2)/2j/np.pi
    j = J * np.sqrt(np.pi/2/z)
    return j * dt_ds

assert np.allclose(j.subs([(z, z_)]), 
                   quad(lambda s: integrand2(s, z_, n), 
                   -np.inf, np.inf, complex_func=True)[0])
Contents