First, we prove the integral formula for $I_{\nu}$. The idea of the proof was originally inspired by Hankel’s representation formula for the reciprocal gamma function. In fact, let $\mathcal{L}\subset\mathbb{C}$ be a contour winding around the negative $Ox$-axis. In our particular case, $\mathcal{L}=\mathcal{L}_{1}+\mathcal{L}_{2}+\mathcal{L}_{3}$, where $\mathcal{L}_{1}$ and $\mathcal{L}_{3}$ are two rays parallel to $Ox$ and $\mathcal{L}_{2}$ is an arc of the unit circle centered at the origin (See Figure A.1). So Hankel’s representation formula gives that

(A.9)

$$\frac{1}{\Gamma(k+\nu+1)}=\frac{1}{2\pi\sqrt{-1}}\int_{\mathcal{L}}e^{w}w^{-(k+\nu+1)}dw,\ w\in\mathbb{C}.$$

By the power series definition of $I_{\nu}$,

(A.10)		$\displaystyle I_{\nu}(y)$	$\displaystyle=$	$\displaystyle\sum\limits_{k=0}^{\infty}\frac{1}{\Gamma(k+1)\Gamma(k+\nu+1)}\Big(\frac{y}{2}\Big)^{2k+\nu}$
			$\displaystyle=$	$\displaystyle(\frac{y}{2})^{\nu}\frac{1}{2\pi\sqrt{-1}}\int_{\mathcal{L}}e^{w}w^{-\nu-1}\sum\limits_{k=0}^{\infty}\frac{(\frac{y^{2}}{4w})^{k}}{k!}dw$
			$\displaystyle=$	$\displaystyle(\frac{y}{2})^{\nu}\frac{1}{2\pi\sqrt{-1}}\int_{\mathcal{L}}e^{w+\frac{y^{2}}{4w}}w^{-\nu-1}dw.$

For every $y>0$, we make change of variables for each $w\in\mathbb{C}$,

(A.11)

$$w=\frac{y\cdot e^{\zeta}}{2}=\frac{ye^{t}}{2}\cdot e^{\sqrt{-1}\theta},\ 0<t<\infty,\ 0\leq\theta\leq 2\pi.$$

Letting $\mathcal{L}_{1}$ and $\mathcal{L}_{3}$ tend to each other, then in terms of the variables $(t,\theta)$,

(A.12)

$$\int_{\mathcal{L}}e^{w+\frac{y^{2}}{4w}}w^{-\nu-1}dw=\frac{1}{\pi}\int_{0}^{\pi}e^{y\cos\theta}\cos(\nu\theta)d\theta-\frac{\sin(\nu\pi)}{\pi}\int_{0}^{\infty}e^{-y\cosh t-\nu t}dt.$$

The integral formula for $K_{\nu}$ follows easily from the above integral representation for $I_{\nu}$ and the definition

(A.13)

$$K_{\nu}(y)=\frac{\pi(I_{-\nu}(y)-I_{\nu}(y))}{2\sin(\nu\pi)}.$$

∎