--- jupytext: formats: ipynb,md:myst,py:light text_representation: extension: .md format_name: myst format_version: 0.13 jupytext_version: 1.13.8 kernelspec: display_name: Python 3 language: python name: python3 --- ```{code-cell} :tags: [hide-cell] 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 ``` (sec:BesselFunctions)= # 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*} :::{admonition} Do It! Derive this from Laplacian in spherical coordinates. :class: dropdown Review {ref}`sec: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 {ref}`sec:A3-5` of {ref}`sec:Assignment3` 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*} :::{admonition} Do It! Show that if $F_{n}(z)$ satisfies Bessel's equation, $f_{\ell}(z)$ solves the ODE above. :class: dropdown 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 ```{code-cell} :tags: [hide-cell] 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)]))) ``` :::{admonition} Do It! Derive these solutions. :class: dropdown 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](sec:SturmLiouville): \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*} :::{admonition} Do It! Check that $\mathcal{L}$ is self-adjoint. :class: dropdown 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 {cite}`Arfken: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: \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} :::{admonition} Do It! Show these from the series expansion. :class: dropdown This is a bit messy, but you should be able to demonstrate that these are correct from the series representations. ::: :::{dropdown} (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. :::{admonition} Do It! Prove this using induction. :class: dropdownn 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*} ::: :::{admonition} Do It! Check that these are consistent with the series representation. :class: dropdown 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. ::: :::{dropdown} (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$. ::: ```{code-cell} :tags: [hide-cell] 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()) ``` ## Generating Functions :::{margin} This partly motivates the sign choice in the definition of $j_{-1}(z) = -y_0(z)$ and $y_{-1}(z) = j_{0}(z)$. ::: From {cite}`Olver: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)* ```{code-cell} ``` ## Integral Representations ```{code-cell} :tags: [margin, hide-input] 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); ``` :::{margin} Contour for the Schläfli integral to compute $j_{n}(z)$ for arbitrary $n$. This avoids the branch cut and branch point. Note: The book {cite}`Arfken:2013` has an error in (14.16): the exponent $e^{z/2)(t - 1/t)}$ here is correct (see code check). ::: 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*} :::{admonition} Do It! Check this! :class: dropdown (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*} ::: ```{code-cell} :tags: [hide-input] 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]) ``` [Frobenius method]: [spherical harmonics]: [Bessel Functions]: