Spherical Bessel and outgoing Hankel recurrences solve the radial Helmholtz equation

Write $f_0(z)=\sin z\mathrm{/}z$ and $g_0(z)=-i{e}^{iz}\mathrm{/}z$. By [F2]–[F3], $f_0^{\mathrm{\prime }}=\cos z\mathrm{/}z-\sin z\mathrm{/}{z}^2$ and $f_0^{\mathrm{\prime }\mathrm{\prime }}=-\sin z\mathrm{/}z-2\cos z\mathrm{/}{z}^2+2\sin z\mathrm{/}{z}^3$, so $L_0{f}_0=0$. Also $g_0^{\mathrm{\prime }}=e^{iz}(1\mathrm{/}z+i\mathrm{/}{z}^2)$ and $g_0^{\mathrm{\prime }\mathrm{\prime }}=e^{iz}(i\mathrm{/}z-2\mathrm{/}{z}^2-2i\mathrm{/}{z}^3)$, so $L_0{g}_0=0$. For any function with the required finite derivatives on $(0,\mathrm{\infty })$, put $y_{\mathrm{\ell }}=(-z{)}^{\mathrm{\ell }}{D}^{\mathrm{\ell }}f$ with $D=z^{-1}d\mathrm{/}dz$. Since $(D^{\mathrm{\ell }}f{)}^{\mathrm{\prime }}=z{D}^{\mathrm{\ell }+1}f$, differentiating the product gives $y_{\mathrm{\ell }}^{\mathrm{\prime }}=\mathrm{\ell }{y}_{\mathrm{\ell }}\mathrm{/}z-y_{\mathrm{\ell }+1}$, hence $y_{\mathrm{\ell }+1}=\mathrm{\ell }{y}_{\mathrm{\ell }}\mathrm{/}z-y_{\mathrm{\ell }}^{\mathrm{\prime }}$. This recurrence applies to both functions in [F1].

The sine series gives $\sin z\mathrm{/}z={\sum }_{m\ge 0}(-1{)}^m{z}^{2m}\mathrm{/}(2m+1)!$. After $r$ applications of $D$, the formal series is ${\sum }_{m\ge r}(-1{)}^m{2}^{r}m!z^{2m-2r}\mathrm{/}((m-r)!(2m+1)!)$. For $z\mathrm{\ne }0$ its successive nonzero absolute-term ratio is $\mathrm{\mid }z{\mathrm{\mid }}^2(m+1)\mathrm{/}((m+1-r)(2m+2)(2m+3))\to 0$; at $z=0$ only its constant term remains. Thus each such series has infinite radius by [F4]. Termwise differentiation and division by $z$ give the next series: the $m=r$ term is constant and differentiates to zero, so no negative power appears. The operation at $r=\mathrm{\ell }$ kills $m<\mathrm{\ell }$; its first remaining coefficient is $2^{\mathrm{\ell }}\mathrm{\ell }!\mathrm{/}(2\mathrm{\ell }+1)!=1\mathrm{/}(2\mathrm{\ell }+1)!!$ by [F9]. The later powers are $z^{\mathrm{\ell }+2},z^{\mathrm{\ell }+4},\dots$; their coefficient series is bounded on $\mathrm{\mid }z\mathrm{\mid }\le 1$ by compact uniform convergence in [F4], so their sum is $O(z^{\mathrm{\ell }+2})$. This proves the stated origin expansion for $j_{\mathrm{\ell }}$.

Put $\lambda =\mathrm{\ell }(\mathrm{\ell }+1)$ and $\mu =(\mathrm{\ell }+1)(\mathrm{\ell }+2)$. Direct differentiation gives $L_{\mathrm{\ell }+1}(\mathrm{\ell }y\mathrm{/}z-y^{\mathrm{\prime }})=-y^{\mathrm{\prime }\mathrm{\prime }\mathrm{\prime }}+(\mathrm{\ell }-2)y^{\mathrm{\prime }\mathrm{\prime }}\mathrm{/}z+(-1+\mu \mathrm{/}{z}^2)y^{\mathrm{\prime }}+(\mathrm{\ell }\mathrm{/}z-\mathrm{\ell }\mu \mathrm{/}{z}^3)y=\mathrm{\ell }{L}_{\mathrm{\ell }}y\mathrm{/}z-(L_{\mathrm{\ell }}y{)}^{\mathrm{\prime }}$. Thus if $L_{\mathrm{\ell }}y=0$, the recurrence in step 1.1 gives $L_{\mathrm{\ell }+1}{y}_{\mathrm{\ell }+1}=0$. Starting from step 1.1 proves both radial equations by induction. For $\mathrm{\ell }\ge 1$, differentiate $y_{\mathrm{\ell }}=(\mathrm{\ell }-1)y_{\mathrm{\ell }-1}\mathrm{/}z-y_{\mathrm{\ell }-1}^{\mathrm{\prime }}$ and substitute $L_{\mathrm{\ell }-1}{y}_{\mathrm{\ell }-1}=0$; this yields $y_{\mathrm{\ell }}^{\mathrm{\prime }}=(\mathrm{\ell }+1)y_{\mathrm{\ell }-1}^{\mathrm{\prime }}\mathrm{/}z+(1-({\mathrm{\ell }}^2-1)\mathrm{/}{z}^2)y_{\mathrm{\ell }-1}$. Substituting $y_{\mathrm{\ell }-1}^{\mathrm{\prime }}=(\mathrm{\ell }-1)y_{\mathrm{\ell }-1}\mathrm{/}z-y_{\mathrm{\ell }}$ from step 1.1 and collecting terms gives $y_{\mathrm{\ell }}^{\mathrm{\prime }}=y_{\mathrm{\ell }-1}-(\mathrm{\ell }+1)y_{\mathrm{\ell }}\mathrm{/}z$.

Write $h_{\mathrm{\ell }}^{(1)}=C_{\mathrm{\ell }}{e}^{iz}$, starting with $C_0=-i\mathrm{/}z$. Applying the recurrence of step 1.1 and the product rule in [F2]–[F3] gives $C_{\mathrm{\ell }+1}=\mathrm{\ell }{C}_{\mathrm{\ell }}\mathrm{/}z-C_{\mathrm{\ell }}^{\mathrm{\prime }}-i{C}_{\mathrm{\ell }}$. Induction shows $C_{\mathrm{\ell }}$ is a finite Laurent polynomial with exponents $-1,\dots ,-\mathrm{\ell }-1$. If its most singular term is $a_{\mathrm{\ell }}{z}^{-\mathrm{\ell }-1}$, the coefficient of $z^{-\mathrm{\ell }-2}$ in $C_{\mathrm{\ell }+1}$ is $(2\mathrm{\ell }+1)a_{\mathrm{\ell }}$; the term $-i{C}_{\mathrm{\ell }}$ is less singular. Since $a_0=-i$, this gives $a_{\mathrm{\ell }}=-i(2\mathrm{\ell }-1)!!$. Also $e^{iz}=1+O(z)$ as $z\downarrow 0$ by [F3] and differentiability at $0$, so the outgoing origin expansion follows. For infinity, induct on $\mathrm{\ell }$ using step 1.1. The base terms are $j_0=\sin z\mathrm{/}z$ and $h_0^{(1)}=-i{e}^{iz}\mathrm{/}z$. If $j_{\mathrm{\ell }}=\sin \theta \mathrm{/}z+R_{\mathrm{\ell }}$, $\theta =z-\mathrm{\ell }\pi \mathrm{/}2$, with $R_{\mathrm{\ell }}$ a finite sum of $z^{-k}\sin z$ and $z^{-k}\cos z$ terms for $k\ge 2$, then $R_{\mathrm{\ell }}^{\mathrm{\prime }}=O(z^{-2})$ by [F3],[F7],[F8], and $\mathrm{\ell }{j}_{\mathrm{\ell }}\mathrm{/}z=O(z^{-2})$; hence $j_{\mathrm{\ell }+1}=-\cos \theta \mathrm{/}z+O(z^{-2})=\sin (z-(\mathrm{\ell }+1)\pi \mathrm{/}2)\mathrm{/}z+O(z^{-2})$ by the addition identities [F8]. For $h_{\mathrm{\ell }}^{(1)}$, the coefficient recurrence gives a finite form $e^{iz}(b_{\mathrm{\ell }}\mathrm{/}z+{\sum }_{k=2}^{\mathrm{\ell }+1}{b}_{\mathrm{\ell }k}{z}^{-k})$. Its leading coefficient obeys $b_{\mathrm{\ell }+1}=-i{b}_{\mathrm{\ell }}$ because the terms $\mathrm{\ell }{C}_{\mathrm{\ell }}\mathrm{/}z$ and $C_{\mathrm{\ell }}^{\mathrm{\prime }}$ have powers at most $z^{-2}$; thus $b_{\mathrm{\ell }}=(-i{)}^{\mathrm{\ell }+1}$. The remainder and its derivative are $O(z^{-2})$ by [F3],[F7],[F8], so the recurrence gives $h_{\mathrm{\ell }+1}^{(1)}=(-i{)}^{\mathrm{\ell }+2}{e}^{iz}\mathrm{/}z+O(z^{-2})$. Differentiating the leading terms and finite remainders also gives $j_{\mathrm{\ell }}^{\mathrm{\prime }}=\cos (z-\mathrm{\ell }\pi \mathrm{/}2)\mathrm{/}z+O(z^{-2})$ and $(h_{\mathrm{\ell }}^{(1)}{)}^{\mathrm{\prime }}=i(-i{)}^{\mathrm{\ell }+1}{e}^{iz}\mathrm{/}z+O(z^{-2})$.

For two solutions $y,v$ of $L_{\mathrm{\ell }}y=0$, differentiation and substitution of their equations gives $(y{v}^{\mathrm{\prime }}-y^{\mathrm{\prime }}v{)}^{\mathrm{\prime }}=-2(y{v}^{\mathrm{\prime }}-y^{\mathrm{\prime }}v)\mathrm{/}z$, hence $(z^{2}W[y,v]{)}^{\mathrm{\prime }}=0$. Applying [F6] separately to the real and imaginary parts on $(0,\mathrm{\infty })$ shows $z^{2}W[j_{\mathrm{\ell }},h_{\mathrm{\ell }}^{(1)}]$ is a constant complex number. With $\theta =z-\mathrm{\ell }\pi \mathrm{/}2$, the expansions and derivative expansions of step 2.2 give $z^{2}W[j_{\mathrm{\ell }},h_{\mathrm{\ell }}^{(1)}]⟶(-i{)}^{\mathrm{\ell }+1}{e}^{iz}(i\sin \theta -\cos \theta )=i$. Indeed, Euler's formula and the addition identities give $e^{iz}(\cos \theta -i\sin \theta )=e^{i(z-\theta )}=e^{i\mathrm{\ell }\pi \mathrm{/}2}$; repeated use of $\sin (\pi \mathrm{/}2)=1$ and $\cos (\pi \mathrm{/}2)=0$ gives $e^{i\mathrm{\ell }\pi \mathrm{/}2}=i^{\mathrm{\ell }}$. Since $i\sin \theta -\cos \theta =-(\cos \theta -i\sin \theta )$ and $(-i{)}^{\mathrm{\ell }+1}{i}^{\mathrm{\ell }}=-i$, the displayed limit is $i$. The constant is therefore $i$, which proves the Wronskian formula.