I want to prove the following identity calculating the integral of an exponential over an even dimensional sphere in terms of functions $\chi_i(R)$ and $\tilde\psi_i(s)$ (described below) which are essentially modified spherical Bessel functions.

First set up some notation. Let $p>0$ be an integer and $R>0$. Let $S_R^{2p}\subset \mathbb{R}^{2p+1}$ denote the radius $R$ sphere. Let $\mathrm{s}\in \mathbb{R}^{2p+1}$ be a point inside the sphere and let $s=|\mathrm{s}|$ so $0\le s<R$. Finally, let $\omega_n$ be the volume of the unit $n$-ball. I want to show the following. $$ \frac{1}{(2p+1)!\,\omega_{2p+1}} \int_{\mathrm{x}\in S_R^{2p}} e^{-\left|\mathrm{x}-\mathrm{s}\right|}\,\mathrm{d}\mathrm{x} = \frac{(-1)^p e^{-R}}{2^p p!}\sum_{i=0}^p \binom{p}{i}\chi_{p+i}(R)\tilde\psi_i(s). $$

I can get Sage to check this is true up to, say $p=18$, using the integral form of the left hand side given below. I had previously asked for a reference for the integral, but as none was forthcoming I've posted the whole statement in the hope of a proof!

Here $(\chi_i)_{i=0}^\infty$ denote the sequence of 'reverse Bessel polynomials', so that $\chi_i(R)$ is a degree $i$ integer polynomial in $R$. The sequence begins as follows: \begin{align*} \chi_0(R)&=1;\\ \chi_1(R)&=R;\\ \chi_2(R)&=R^2+R;\\ \chi_3(R)&=R^3+3R^2+3R. %;\\ %\chi_4(R)&=R^{4} + 6 R^{3} + 15 R^{2} + 15 R. %\chi_5(R)&=R^{5} + 10 R^{4} + 45 R^{3} + 105 R^{2} + 105 R \end{align*} There are many ways to define this sequence, but we can take the recursion relation as a definition: $$ \chi_{i+2}(R)=R^2\chi_i(R)+(2i+1)\chi_{i+1}(R). $$

These functions are related to the modified spherical Bessel functions by $\chi_i(R)=\frac{2}{\pi} e^R R^{i+1}k_{i-1}(R)$. I will describe an integral form at the bottom.

The sequence $(\tilde\psi_i)_{i=0}^\infty$ denotes the sequence of functions $\mathbb{R}\to \mathbb{R}$ which begins in the following way: \begin{align*} \tilde\psi_0(s)&=\cosh(s);\\ \tilde\psi_1(s) &=-\frac{\sinh\left(s\right)}{s};\\ \tilde\psi_2(s) &=\frac{\cosh\left(s\right)}{s^{2}} - \frac{\sinh\left(s\right)}{s^{3}};\\ \tilde\psi_3(s) &= -\frac{\sinh\left(s\right)}{s^{3}} + \frac{3 \, \cosh\left(s\right)}{s^{4}} - \frac{3 \, \sinh\left(s\right)}{s^{5}}. %;\\ %\tilde\psi_4(s) %&= %\frac{\cosh\left(s\right)}{s^{4}} - \frac{6 \, \sinh\left(s\right)}{s^{5}} + \frac{15 \, \cosh\left(s\right)}{s^{6}} - %\frac{15 \, \sinh\left(s\right)}{s^{7}}. \end{align*} You can take the following recursion relation as a definition. $$ \tilde\psi_{i+1}(s)=-\frac{1}{s}\frac{\mathrm d \tilde\psi_i(s)}{\mathrm{d} s}. $$ These functions are related to the modified spherical Bessel functions by $\tilde\psi_i(s)=(-1)^i i^{(1)}_{i-1}(s)/s^{i-1}$. Again there are integral forms which I will give at the bottom.

The left hand side of the identity I want to prove can be expressed as a line integral as follows which makes it more tractable. \begin{multline*} \frac{1}{(2p+1)!\,\omega_{2p+1}} \int_{\mathrm{x}\in S^{2p}} e^{-\left|\mathrm{x}-\mathrm{s}\right|}\,\mathrm{d}\mathrm{x}\\ =\frac{ R^{2p}}{p!(p-1)!(2)^{2p}}\int_{\theta=0}^{\pi}e^{-\sqrt{R^2+s^2-2Rs\cos \theta}}\sin^{2p-1}\theta\,\mathrm{d}\theta\\ =\frac{ R}{p!(p-1)!(4s)^{2p-1}}\int_{\rho=R-s}^{R+s}e^{-\rho}((R+s)^2-\rho^2)^{p-1} (\rho^2-(R-s)^2)^{p-1}\rho\,\mathrm{d}\rho. \end{multline*}

For $i\ge 1$ and $R, s>0$ we have the following integral forms for the functions. \begin{align*} \chi_i(R) &= \frac{e^R R^{2i}}{2^{i-1} (i-1)!}\int_{t=0}^\infty e^{-R\cosh t}\sinh^{2i-1}t\, \mathrm d t\\ &= \frac{e^R R}{2^{i-1} (i-1)!}\int_{y=R}^\infty e^{-y}(y^2-R^2)^{i-1}\, \mathrm d y. \end{align*} \begin{align*} \tilde\psi_i(s) &= \frac{(-1)^{i}}{2^{i} (i-1)!}\int_{\theta=0}^\pi e^{s\cos \theta}\sin^{2i-1}\theta\, \mathrm d \theta\\ &= \frac{(-1)^{i}}{2^{i} (i-1)!s^{2i-1}}\int_{x=-s}^s e^{x}(s^2-x^2)^{i-1}\, \mathrm d x. \end{align*}