# Solution of an inhomogeneous modified Bessel equation

I'm solving the equation

$$x^2y''+xy'-(x^{2}\lambda^{2}+1)y=-C\frac{I_{1}(\lambda x)}{I_{1}(\lambda)}$$

where $I_{\alpha}$ is the modified Bessel function of the first kind.

The complementary function for this problem is

$$y_{cf}=AI_{1}(\lambda x)+BK_{1}(\lambda x),$$

where $K_{\alpha}$ is the modified Bessel function of the second kind.

Using the method of variation of parameters I find that

$$y_{p}=u_{1}I_{1}(\lambda x)+u_{2}K_{1}(\lambda x),$$

where

\begin{align*} u_{1}&=\frac{C}{I_{1}(\lambda)}\int\frac{I_{1}(\lambda x)K_{1}(\lambda x)} {W[I_{1}(\lambda x),K_{1}(\lambda x)]}\,\textrm{d}x =-\frac{C}{I_{1}(\lambda)}\int xI_{1}(\lambda x)K_{1}(\lambda x)\,\textrm{d}x, \\ u_{2}&=-\frac{C}{I_{1}(\lambda)}\int\frac{I_{1}(\lambda x)I_{1}(\lambda x)} {W[I_{1}(\lambda x),K_{1}(\lambda x)]}\,\textrm{d}x =\frac{C}{I_{1}(\lambda)}\int xI_{1}(\lambda x)I_{1}(\lambda x)\,\textrm{d}x. \end{align*}

Here I have used the fact that

$$W[I_{1}(\lambda x),K_{1}(\lambda x)]=I_{1}(\lambda x)[K_{1}(\lambda x)]'-K_{1}(\lambda x)[I_{1}(\lambda x)]'=-\frac{1}{x}.$$

Can anyone suggest the best way to go about computing these integrals? I've tried 'by parts' (seemed the natural choice) but haven't managed to get anywhere thus far - thanks!

Consider the following differential equation $$\frac{{{d}^{2}}u}{d{{z}^{2}}}-{{k}^{2}}u-\frac{q\left( q+1 \right)}{{{z}^{2}}}u=0$$ Note firstly that if we set $y=\frac{u}{\sqrt{z}}$ then we have $${{z}^{2}}y''+zy'-\left( {{k}^{2}}{{z}^{2}}+{{\left( q+\tfrac{1}{2} \right)}^{2}} \right)y=0$$ Then we set $q+1/2=\nu ,\,\,ikz=x$and so obtain $${{x}^{2}}y''+xy'+\left( {{x}^{2}}-{{\nu }^{2}} \right)y=0$$ Bessel’s DE. Hence $y={{Z}_{\nu }}\left( x \right)$ where Z is the cylindrical function of interest and of order $\nu$. Inverting the transform therefore we see $u={{z}^{1/2}}{{Z}_{q+\tfrac{1}{2}}}\left( ikz \right)$ . In a similar fashion if we set $u=w{{z}^{-q}}$ we obtain $$\frac{{{d}^{2}}w}{d{{z}^{2}}}-\frac{2q}{z}\frac{dw}{dz}-{{k}^{2}}w=0$$ Which has a solution $w={{z}^{q+1/2}}{{Z}_{q+\tfrac{1}{2}}}\left( icz \right)$ . We now get to a result due to Lommel, one which I will simply state. Replace $w=y/\theta \left( z \right)$ and $z\to \psi \left( z \right)$ and $2q$ by $2\nu -1$ and $k=i$ , then \begin{align} & \frac{{{d}^{2}}y}{d{{z}^{2}}}-\left( \frac{\psi ''}{\psi '}+\left( 2\nu -1 \right)\frac{\psi '}{\psi }+2\frac{\theta '}{\theta } \right)\frac{dy}{dz} \\ & +\left\{ \left( \frac{\psi ''}{\psi '}+\left( 2\nu -1 \right)\frac{\psi '}{\psi }+2\frac{\theta '}{\theta } \right)\frac{\theta '}{\theta }-\frac{\theta ''}{\theta }+\psi {{'}^{2}} \right\}y=0 \\ \end{align} The solution of which is $y=\theta \left( z \right)\psi {{\left( z \right)}^{\nu }}{{Z}_{\nu }}\left( \psi \left( z \right) \right)$ . Now introduce a new function defined by $$\frac{\varphi '}{\varphi }=\frac{\psi ''}{\psi '}+\left( 2\nu -1 \right)\frac{\psi '}{\psi }+2\frac{\theta '}{\theta }$$ Using this we eliminate $\theta$ to obtain \begin{align} & \frac{{{d}^{2}}y}{d{{z}^{2}}}-\frac{\varphi '}{\varphi }\frac{dy}{dz} \\ & +\left\{ \frac{3}{4}{{\left( \frac{\varphi '}{\varphi } \right)}^{2}}-\frac{1}{2}\frac{\varphi ''}{\varphi }-\frac{3}{4}{{\left( \frac{\psi ''}{\psi '} \right)}^{2}}+\frac{1}{2}\frac{\psi '''}{\psi '}+\left( {{\psi }^{2}}-\nu +\frac{1}{4} \right){{\left( \frac{\psi '}{\psi } \right)}^{2}} \right\}y=0 \\ \end{align} The solution of which has the form $$y=\sqrt{\frac{\varphi \psi }{\psi '}}{{Z}_{\nu }}\left( \psi \right)$$ Now let $\varphi =1$ and hence
$$\frac{{{d}^{2}}y}{d{{z}^{2}}}+\left\{ \frac{1}{2}\frac{\psi '''}{\psi '}-\frac{3}{4}{{\left( \frac{\psi ''}{\psi '} \right)}^{2}}+\left( {{\psi }^{2}}-\nu +\frac{1}{4} \right){{\left( \frac{\psi '}{\psi } \right)}^{2}} \right\}y=0$$ Now observe that if ${{y}_{\nu }}''+P{{y}_{\nu }}=0$ and ${{y}_{\mu }}''+Q{{y}_{\mu }}=0$ then we have $${{y}_{\nu }}{{y}_{\mu }}''-{{y}_{\nu }}''{{y}_{\mu }}=\left( P-Q \right){{y}_{\nu }}{{y}_{\mu }}$$ or $$\frac{d}{dx}\left( {{y}_{\nu }}{{y}_{\mu }}'-{{y}_{\nu }}'{{y}_{\mu }} \right)=\left( P-Q \right){{y}_{\nu }}{{y}_{\mu }}$$ or $$\int{\left( P-Q \right){{y}_{\nu }}{{y}_{\mu }}dx}={{y}_{\nu }}{{y}_{\mu }}'-{{y}_{\nu }}'{{y}_{\mu }}$$ Here we have $$P,Q=\left\{ \frac{1}{2}\frac{\psi '''}{\psi '}-\frac{3}{4}{{\left( \frac{\psi ''}{\psi '} \right)}^{2}}+\left( {{\psi }^{2}}-\nu +\frac{1}{4} \right){{\left( \frac{\psi '}{\psi } \right)}^{2}} \right\}$$ We now let $\psi =az,\,\,\nu =\nu$ for P, and $\psi =bz,\,\,\nu =\mu$for Q. This yields $$P=\left( {{a}^{2}}{{z}^{2}}-\nu +\frac{1}{4} \right)\frac{1}{{{z}^{2}}},\,\,Q=\left( {{b}^{2}}{{z}^{2}}-\mu +\frac{1}{4} \right)\frac{1}{{{z}^{2}}}$$ Let${{y}_{\nu }}=\sqrt{z}{{Z}_{\nu }}\left( az \right)$and ${{y}_{\mu }}=\sqrt{z}{{\bar{Z}}_{\mu }}\left( bz \right)$ where the bar on Z here indicates that it is possibly a different cylinder function than that in ${{y}_{\nu }}$. We have then \begin{align} & \int{\left( \left( {{a}^{2}}-{{b}^{2}} \right)z+\left( \mu -\nu \right)\frac{1}{z} \right){{Z}_{\nu }}\left( az \right){{{\bar{Z}}}_{\mu }}\left( bz \right)dz} \\ & =z\left( {{Z}_{\nu }}\left( az \right)\frac{d}{dz}{{{\bar{Z}}}_{\mu }}\left( bz \right)-{{{\bar{Z}}}_{\mu }}\left( bz \right)\frac{d}{dz}{{Z}_{\nu }}\left( az \right) \right) \\ \end{align} Now let $\nu =\mu$ and so $$\int{z{{Z}_{\mu }}\left( az \right){{{\bar{Z}}}_{\mu }}\left( bz \right)dz}=\frac{z\left( {{Z}_{\mu }}\left( az \right)\frac{d}{dz}{{{\bar{Z}}}_{\mu }}\left( bz \right)-{{{\bar{Z}}}_{\mu }}\left( bz \right)\frac{d}{dz}{{Z}_{\mu }}\left( az \right) \right)}{\left( {{a}^{2}}-{{b}^{2}} \right)}$$ Using L’Hopital’s rule we have $$\int{z{{Z}_{\mu }}\left( az \right){{{\bar{Z}}}_{\mu }}\left( az \right)dz}=\frac{{{z}^{2}}}{4}\left( 2{{Z}_{\mu }}\left( az \right){{{\bar{Z}}}_{\mu }}\left( az \right)-{{Z}_{\mu -1}}\left( az \right){{{\bar{Z}}}_{\mu +1}}\left( az \right)-{{Z}_{\mu +1}}\left( az \right){{{\bar{Z}}}_{\mu -1}}\left( az \right) \right)$$ Therefore, for example $$\int{z{{I}_{1}}\left( az \right){{I}_{1}}\left( az \right)dz}=\frac{{{z}^{2}}}{2}\left( {{I}_{1}}{{\left( az \right)}^{2}}-{{I}_{0}}\left( az \right){{I}_{2}}\left( az \right) \right)$$

All of this is in Watson's treatise (pg 134-135), although a lot more detail has been added here.

• Thanks for the very thorough response. So am I correct in thinking that $$\int zI_{1}(az)K_{1}(az)\,\textrm{d}z=\frac{z^{2}}{4}(2I_{1}(az)K_{1}(az)-I_{0}(az)K_{2}(az)-I_{2}(az)K_{0}(az))$$ ? When I differentiate this expression I can't seem to return $$zI_{1}(az)K_{1}(az)$$ Unless I am making an error when differentiating these modified Bessel functions, which is wholly possible. Again, thanks for all the help. – Juggler Oct 24 '17 at 11:12
• Hmmm that is interesting. I must admit that i only looked at the I-based integrals because i have used such things in the past. The I-K combination is unkown to me. I think the problem lies with L'Hopital's rule. I believe (at this stage) that the limit doesn't exist for that integral and so the rule is not applicable. I'll test all this and post another "solution" – mathstackuser12 Oct 24 '17 at 22:07

To examine why this method yields one integral but fails on the other let us remove some of the complexity by removing some generality and focus on the specific problem at hand. Therefore consider $${{z}^{2}}y''+zy'-\left( {{k}^{2}}{{z}^{2}}+1 \right)y=0$$ The general solution of which is $$y=A{{I}_{1}}\left( kz \right)+B{{K}_{1}}\left( kz \right)$$ So using $y=\frac{u}{\sqrt{z}}$we have $$\frac{{{d}^{2}}{{u}_{k}}}{d{{z}^{2}}}-\left( {{k}^{2}}+\frac{3}{4{{z}^{2}}} \right){{u}_{k}}=0$$ And thus $u=A\sqrt{z}{{I}_{1}}\left( kz \right)+B\sqrt{z}{{K}_{1}}\left( kz \right)$ which we can confirm directly. Therefore directly from the DE we have $${{\bar{u}}_{p}}\frac{{{d}^{2}}{{u}_{k}}}{d{{z}^{2}}}=\left( {{k}^{2}}+\frac{3}{4{{z}^{2}}} \right){{u}_{k}}{{\bar{u}}_{p}}$$ $${{u}_{k}}\frac{{{d}^{2}}{{{\bar{u}}}_{p}}}{d{{z}^{2}}}=\left( {{p}^{2}}+\frac{3}{4{{z}^{2}}} \right){{u}_{k}}{{\bar{u}}_{p}}$$ So $${{\bar{u}}_{p}}\frac{{{d}^{2}}{{u}_{k}}}{d{{z}^{2}}}-{{u}_{k}}\frac{{{d}^{2}}{{{\bar{u}}}_{p}}}{d{{z}^{2}}}=\left( {{k}^{2}}-{{p}^{2}} \right){{u}_{k}}{{\bar{u}}_{p}}$$ $$\int{\left( {{k}^{2}}-{{p}^{2}} \right){{u}_{k}}{{{\bar{u}}}_{p}}dz}={{\bar{u}}_{p}}\frac{d{{u}_{k}}}{dz}-{{u}_{k}}\frac{d{{{\bar{u}}}_{p}}}{dz}$$ Where $\bar{u}$represents a different linear combination of solutions to the DE, i.e. $${{u}_{k}}=A\sqrt{z}{{I}_{1}}\left( kz \right)+B\sqrt{z}{{K}_{1}}\left( kz \right)$$ $${{\bar{u}}_{p}}=C\sqrt{z}{{I}_{1}}\left( pz \right)+D\sqrt{z}{{K}_{1}}\left( pz \right)$$ $$\int{{{u}_{k}}{{{\bar{u}}}_{p}}dz}=\frac{{{{\bar{u}}}_{p}}\frac{d{{u}_{k}}}{dz}-{{u}_{k}}\frac{d{{{\bar{u}}}_{p}}}{dz}}{\left( {{k}^{2}}-{{p}^{2}} \right)}$$ Let us now confirm the first integral. Consider therefore B=D=0, and A=C=1, hence $${{u}_{k}}=\sqrt{z}{{I}_{1}}\left( kz \right),\,\,\frac{d{{u}_{k}}}{dz}=\frac{{{I}_{1}}\left( kz \right)-2kz{{I}_{0}}\left( kz \right)}{2\sqrt{z}}\,$$ $${{\bar{u}}_{p}}=\sqrt{z}{{I}_{1}}\left( pz \right),\,\,\,\frac{d{{u}_{p}}}{dz}=\frac{{{I}_{1}}\left( pz \right)-2pz{{I}_{0}}\left( pz \right)}{2\sqrt{z}}$$ We obtain $$\int{z{{I}_{1}}\left( kz \right){{I}_{1}}\left( pz \right)dz}=\frac{pz{{I}_{0}}\left( pz \right){{I}_{1}}\left( kz \right)-kz{{I}_{0}}\left( kz \right){{I}_{1}}\left( pz \right)}{{{k}^{2}}-{{p}^{2}}}$$ Notice that in the limit as $p\to k$ both numerator and denominator approach zero, hence invoking L’Hopital’s rule (diff w.r.t to k say), and we have $$\lim \int{z{{I}_{1}}\left( kz \right){{I}_{1}}\left( pz \right)dz}=\lim \frac{z\left( {{I}_{0}}\left( kz \right){{I}_{1}}\left( pz \right)+kz{{I}_{1}}\left( kz \right){{I}_{1}}\left( pz \right)-\tfrac{1}{2}pz{{I}_{0}}\left( pz \right)\left( {{I}_{0}}\left( kz \right)+{{I}_{2}}\left( kz \right) \right) \right)}{2k}$$ Now observe \begin{align} & {{I}_{0}}\left( kz \right){{I}_{1}}\left( pz \right)-\tfrac{1}{2}pz{{I}_{0}}\left( pz \right)\left( {{I}_{0}}\left( kz \right)+{{I}_{2}}\left( kz \right) \right)= \\ & {{I}_{0}}\left( kz \right)\left( {{I}_{1}}\left( pz \right)-\tfrac{1}{2}pz{{I}_{0}}\left( pz \right) \right)-\tfrac{1}{2}pz{{I}_{0}}\left( pz \right){{I}_{2}}\left( kz \right) \\ \end{align} However from the recurrence relations $$\frac{1}{\tfrac{1}{2}pz}{{I}_{1}}\left( pz \right)={{I}_{0}}\left( pz \right)-{{I}_{2}}\left( pz \right)$$ So $${{I}_{0}}\left( kz \right){{I}_{1}}\left( pz \right)-\tfrac{1}{2}pz{{I}_{0}}\left( pz \right)\left( {{I}_{0}}\left( kz \right)+{{I}_{2}}\left( kz \right) \right)=-\tfrac{1}{2}pz\left( {{I}_{0}}\left( kz \right){{I}_{2}}\left( pz \right)+{{I}_{0}}\left( pz \right){{I}_{2}}\left( kz \right) \right)$$ Hence in the limit we obtain $$\int{z{{I}_{1}}\left( kz \right){{I}_{1}}\left( kz \right)dz}=\frac{{{z}^{2}}}{2}\left( {{I}_{1}}{{\left( kz \right)}^{2}}-{{I}_{0}}\left( kz \right){{I}_{2}}\left( kz \right) \right)$$ Now consider the same thing again except we place A=D=1, B=C=0, we have then $${{u}_{k}}=\sqrt{z}{{I}_{1}}\left( kz \right),\,\,\frac{d{{u}_{k}}}{dz}=\frac{{{I}_{1}}\left( kz \right)-2kz{{I}_{0}}\left( kz \right)}{2\sqrt{z}}\,$$ $${{\bar{u}}_{p}}=\sqrt{z}{{K}_{1}}\left( pz \right),\,\,\,\frac{d{{u}_{p}}}{dz}=-\frac{{{K}_{1}}\left( pz \right)+2pz{{K}_{0}}\left( pz \right)}{2\sqrt{z}}$$ Therefore $$\int{z{{I}_{1}}\left( kz \right){{K}_{1}}\left( pz \right)dz}=\frac{pz{{K}_{0}}\left( pz \right){{I}_{1}}\left( kz \right)+kz{{I}_{0}}\left( kz \right){{K}_{1}}\left( pz \right)}{{{k}^{2}}-{{p}^{2}}}$$ However observe that in this instance in the limit as $p\to k$ the numerator does not approach zero, hence the general indefinite integral arrived at from the limit is not applicable to this particular combination of cylinder functions. As to how one would evaluate this particular combination...I'll have to have more of a think.

Update: the results so far do not allow for the indefinite integral of the $z{{I}_{w}}\left( z \right){{K}_{w}}\left( z \right)$. Towards this end consider the Mellin transform representation of the product $$\int{z{{I}_{1}}\left( z \right){{K}_{1}}\left( z \right)dz}=\frac{1}{2\pi i}\int\limits_{c-i\infty }^{c+i\infty }{\frac{1}{4\left( 1-s \right)}\frac{\Gamma \left( s \right)\Gamma \left( s+1 \right)\Gamma \left( \tfrac{1}{2}-s \right)}{\sqrt{\pi }\Gamma \left( 2-s \right)}{{z}^{2-2s}}ds}$$

$$\left| \arg z \right|<\tfrac{1}{2}\pi ,\,\,0<c<\tfrac{1}{2}$$

as given in Titchmarsh’s “Theory of fourier integrals” 7.10.8 pg.199. This has poles at $\operatorname{Re}\left( s \right)={{\mathbb{Z}}^{-}}\cup \{0\}$ driven by $\Gamma \left( s \right)\Gamma \left( s+1 \right)=s\Gamma {{\left( s \right)}^{2}}$. Consider therefore a contour consisting of a straight line segment $\left[ -iR+c,iR+c \right]$ where $0<c<\tfrac{1}{2}$and a half circle of radius R centered at x=c which closes the contour to the left of the origin. The integrand can be shown to approach zero on the circular contour as R becomes infinite.

Now $$\underset{z=-n}{\mathop{\operatorname{res}}}\,\left\{ \Gamma {{\left( z \right)}^{2}}f\left( z \right) \right\}=\frac{1}{{{\left( n! \right)}^{2}}}\left\{ 2\psi \left( n+1 \right)f\left( -n \right)+f'\left( -n \right) \right\}$$ Proof: Consider a function $g\left( z \right)$ that has a simple pole at z=a, and a function $f\left( z \right)$ that is analytic at z=a. Then $g\left( z \right)=\frac{{{c}_{-1}}}{\left( z-{{z}_{0}} \right)}+{{c}_{0}}+{{c}_{1}}\left( z-{{z}_{0}} \right)+...$ and hence $$g{{\left( z \right)}^{2}}={{\left( \frac{{{c}_{-1}}}{\left( z-a \right)}+{{c}_{0}}+{{c}_{1}}\left( z-a \right)+... \right)}^{2}}=\frac{c_{-1}^{2}}{{{\left( z-a \right)}^{2}}}+\frac{2{{c}_{-1}}{{c}_{0}}}{\left( z-a \right)}+O\left( 1 \right)$$

Therefore \begin{align} & \underset{z=a}{\mathop{res}}\,g{{\left( z \right)}^{2}}f\left( z \right)=\underset{z=a}{\mathop{res}}\,\left( f\left( a \right)+f'\left( a \right)\left( z-a \right)+\frac{1}{2!}f''\left( a \right){{\left( z-a \right)}^{2}}+... \right)\left( \frac{c_{-1}^{2}}{{{\left( z-a \right)}^{2}}}+\frac{2{{c}_{-1}}{{c}_{0}}}{\left( z-a \right)}+O\left( 1 \right) \right) \\ & =2f\left( a \right){{c}_{-1}}{{c}_{0}}+f'\left( a \right)c_{-1}^{2} \\ \end{align} It is reasonably well known $$\Gamma \left( z+n \right)=\frac{{{\left( -1 \right)}^{n}}}{n!}\left( \frac{1}{z+n}+\psi \left( 1+n \right)+O\left( z+n \right) \right)$$ so $${{c}_{-1}}=\frac{{{\left( -1 \right)}^{n}}}{n!},\,\,{{c}_{0}}=\frac{{{\left( -1 \right)}^{n}}}{n!}\psi \left( 1+n \right)$$ where $\psi \left( z \right)$is the digamma function. We have therefore the often-useful result (for Mellin transform inversion) $$\underset{z=-n}{\mathop{res}}\,\Gamma {{\left( z \right)}^{2}}f\left( z \right)=\frac{1}{{{\left( n! \right)}^{2}}}\left( 2\psi \left( 1+n \right)f\left( -n \right)+f'\left( -n \right) \right)$$ Here $$f\left( s \right)=\frac{1}{4\left( 1-s \right)}\frac{s\Gamma \left( \tfrac{1}{2}-s \right)}{\sqrt{\pi }\Gamma \left( 2-s \right)}{{z}^{2-2s}}$$ and so after considerable algebra $$\int{z{{I}_{1}}\left( z \right){{K}_{1}}\left( z \right)dz}=\sum\limits_{n=0}^{\infty }{\frac{\Gamma \left( \tfrac{1}{2}+n \right)\left\{ 1-n+2n\left( 1+n \right)\log \left( z \right)+\left( 1+n \right)\left( \psi \left( \tfrac{1}{2}+n \right)-3\psi \left( 1+n \right) \right) \right\}}{4\sqrt{\pi }\Gamma {{\left( 2+n \right)}^{3}}}}{{z}^{2+2n}}$$ Observe that the n=0 term is simply $\tfrac{1}{4}{{z}^{2}}$. Taking this term out and then using $\Gamma \left( z+{\scriptstyle{}^{1}/{}_{2}} \right)=\frac{\Gamma \left( 2z \right){{\left( 2\pi \right)}^{1/2}}}{\Gamma \left( z \right){{2}^{2z-1/2}}}$ yields the indefinite integral as a series $$\int{z{{I}_{1}}\left( z \right){{K}_{1}}\left( z \right)dz}=\frac{{{z}^{2}}}{4}+\frac{{{z}^{2}}}{2}\sum\limits_{n=1}^{\infty }{\frac{\left\{ 1-n+2n\left( 1+n \right)\log \left( z \right)+\left( 1+n \right)\left( \psi \left( \tfrac{1}{2}+n \right)-3\psi \left( 1+n \right) \right) \right\}\Gamma \left( 2n \right)}{\Gamma {{\left( 2+n \right)}^{3}}\Gamma \left( n \right)}}{{\left( \frac{z}{2} \right)}^{2n}}$$ Bessel functions can be extracted from this horrible thing in the following manner. Consider $${{S}_{1}}=\sum\limits_{n=1}^{\infty }{\frac{\Gamma \left( 2n \right)}{\Gamma {{\left( 2+n \right)}^{3}}\Gamma \left( n \right)}}{{\left( \frac{z}{2} \right)}^{2n}}=\frac{1-z-F\left[ -\tfrac{1}{2};1,1;2z \right]}{2z}$$ Where F is the generalised hypergeometric function. This is relatively straightforward to show (see Wilf’s generatingfunctionology for examples of the algorithm). Mathematica then tells me this has some nasty relationship to the confluent hypergeometric functions which lead straight to modified Bessel functions. I’m assuming if you had enough time the identities between all these geometric functions would be reasonably easy to determine. There are two more terms that are derivatives of the above function, call them ${{S}_{2}},{{S}_{3}}$. All these are therefore also related to the modified Bessel function. The summations involving the digamma function, however, requires perhaps another approach.

• I have searched high and low for book references/internet sources but have come up with nothing. It seems strange to me that this problem doesn't appear to have been addressed more widely before! For reference I'm attempting to follow the analytical solution derived in this paper: aip.scitation.org/doi/abs/10.1063/1.1692738 – Juggler Oct 25 '17 at 7:51
• Well i'd need to get my hands on that paper to know where you're trying to go with this. However the reamining integral in this post appears unlikely to have a 'nice' solution strictly in terms of bessel functions. I say this because indefinite integrals of special functions are hard to come by. The trouble with all this is that tables (such as Bateman's Vol II) simply list the results we have above, as though they work with any Bessel function. As we've seen this is not true. This specific type of integral requires the wronkskian to be zero. – mathstackuser12 Oct 25 '17 at 21:32
• I should add that i have a series representation of the integral if that helps? Although it's fairly disgusting. I might post it as an edit to the above post. – mathstackuser12 Oct 26 '17 at 0:50