An integral related to the pdf of the non-central chi-squared distribution

Problem Statement

I was attempting to produce a direct derivation of the pdf of the non-central chi-square distribution and came across the following integral, which I will denote as $g_{n-1}(x_n, \mathbf{\mu}^{(n)})$.

\begin{align} \int_{(0,1)^{n-1}}\left(\ \prod^{n-1}_{k=1}\frac{x_k^{k/2-1}}{\sqrt{1-x_k}}\ \right)\left(\ \prod^n_{k=1}\cosh\left(\mu_k\sqrt{x_n}\sqrt{\ (1-x_{k-1})\prod^{n-k}_{j=1}x_{n-j}}\right)\ \right)\ dx_1dx_2\cdots dx_{n-1} \end{align}

Here $x_n>0$, $x_0:=0$ and $\mathbb{R}^n\ni\mu^{(n)}:=(\mu_1,\mu_2,\cdots,\mu_n)$.

Thus far, I have succeeded in showing that if $f_n$ is the (Lebesgue) pdf of a random variable having a non-central chi-square distribution with $n$ degrees of freedom and non-centrality parameter $\delta_n:=\sum^n_{k=1}\mu_k^2$, then

$$f_n(x_n, \delta_n)=\frac{e^{(-\delta_n+x_n)/2}}{(2\pi)^{n/2}}x_n^{n/2-1}g_{n-1}(x_n,\mu^{(n)})\tag{1}$$

and I am interested to know how one might complete the proof to show that

$$f_n(x_n,\delta_n)=\frac{e^{(-\delta_n+x_n)/2}}{2}\left(\frac{x_n}{\delta_n}\right)^{n/4-1/2}I_{n/2-1}\left(\sqrt{\delta_nx_n}\right)\tag{2}$$

with $I_\nu$ being the modified Bessel function of the first kind.

If one takes the following for granted:

(*) Claim: The value of $g_{n-1}(x_n,\mu^{(n)})$ depends only on $x_n$ and $\delta_n$.

then finding the value of $g_{n-1}(x_n,\mu^{(n)})$ would be a rather easy task as one could simply take $\mu_1=\sqrt{\delta_n}, \mu_2=\cdots=\mu_n=0$ without loss of generality and the integral simplifies to

\begin{align} g_{n-1}(x_n, \mu^{(n)}) &=\int_{(0,1)^{n-1}}\left(\ \prod^{n-1}_{k=1}\frac{x_k^{k/2-1}}{\sqrt{1-x_k}}\ \right)\cosh\left(\sqrt{\delta_n x_n}\sqrt{x_1\cdots x_{n-1}}\right)\ dx_1\cdots dx_{n-1}\\ &=\int_{(0,1)^{n-1}}\sum^\infty_{j=0}\frac{\left(\left(\sqrt{\delta_n x_n}\right)^2\right)^j}{(2j)!}\left(\ \prod^{n-1}_{k=1}\frac{x_k^{j+k/2-1}}{\sqrt{1-x_k}}\ \right)\ dx_1\cdots dx_{n-1}\\ &=\sum^\infty_{j=0}\frac{\left(\left(\sqrt{\delta_n x_n}\right)^2\right)^j}{(2j)!}\int_{(0,1)^{n-1}}\left(\ \prod^{n-1}_{k=1}\frac{x_k^{j+k/2-1}}{\sqrt{1-x_k}}\ \right)\ dx_1\cdots dx_{n-1}\\ &=\sum^\infty_{j=0}\frac{\left(\left(\sqrt{\delta_n x_n}\right)^2\right)^j}{(2j)!}\left(\ \prod^{n-1}_{k=1}\int^1_0\frac{x_k^{j+k/2-1}}{\sqrt{1-x_k}}\ dx_k\ \right)\\ &=\sum^\infty_{j=0}\frac{\left(\left(\sqrt{\delta_n x_n}\right)^2\right)^j}{(2j)!}\left(\ \prod^{n-1}_{k=1}\frac{\Gamma\left(j+\frac{k}{2}\right)\Gamma\left(\frac{1}{2}\right)}{\Gamma\left(j+\frac{k+1}{2}\right)} \right)\\ &=\sum^\infty_{j=0}\frac{\left(\left(\sqrt{\delta_n x_n}\right)^2\right)^j}{(2j)!}\left(\frac{\Gamma\left(j+\frac{1}{2}\right)}{\Gamma\left(j+\frac{n}{2}\right)}\cdot \pi^{(n-1)/2}\right)\\ &=\sum^\infty_{j=0}\frac{\left(\left(\sqrt{\delta_n x_n}\right)^2\right)^j}{(2j)!}\left(\frac{(2j)!}{4^jj!\Gamma\left(j+\frac{n}{2}\right)}\cdot \pi^{n/2}\right)\\ &=\pi^{n/2}\sum^\infty_{j=0}\frac{\left(\frac{\left(\sqrt{\delta_n x_n}\right)^2}{4}\right)^j}{j!\Gamma\left(j+\frac{n}{2}\right)}=\pi^{n/2}\left(\frac{2}{\sqrt{\delta_n x_n}}\right)^{n/2-1}I_{n/2-1}\left(\sqrt{\delta_n x_n}\right)\\ &=\frac{(2\pi)^{n/2}}{2(\delta_nx_n)^{n/4-1/2}}I_{n/2-1}\left(\sqrt{\delta_n x_n}\right) \end{align} which proves $(2)$ after substituting this into equation $(1)$. Hence, here are my two questions:

$\text{a)}$ How might one rigorously justify the claim (*) above? I am most likely missing something but the veracity of the claim is not immediately obvious to me.

$\text{b)}$ How might one evalaute the integral $g_{n-1}$ without exploiting the claim?

Solution to Question b) for $n=2$

If $n=2$, we have

$$g_1(x_2, \mu^{(2)})=\int^1_0\frac{\cosh(\mu_1\sqrt{x_2}\sqrt{x_1})\cosh(\mu_2\sqrt{x_2}\sqrt{1-x_1})}{\sqrt{x_1(1-x_1)}}\ dx_1$$

Let $(a,b)=(\mu_1\sqrt{x_2}, \mu_2\sqrt{x_2})$ and $2\omega=b+ia$. With the substitution $x_1=\sin^2\theta$ and then $z=e^{i\theta}$, we get

\begin{align} g_1(x_2, \mu^{(2)}) &=\int^\frac{\pi}{2}_02\cosh(a\sin{\theta})\cosh(b\cos\theta)\ d\theta\\ &=\frac{1}{2}\int^{2\pi}_0\cosh(a\sin{\theta})\cosh(b\cos\theta)\ d\theta\\ &=\frac{1}{2}\int_{|z|=1}\cosh\left(a\cdot\frac{z-z^{-1}}{2i}\right)\cosh\left(b\cdot\frac{z+z^{-1}}{2}\right)\frac{dz}{iz}\\ &=\frac{1}{4}\int_{|z|=1}\left[\cosh\left(\bar{\omega}z+\omega z^{-1}\right)+\cosh\left(\omega z+\bar{\omega}z^{-1}\right)\right]\frac{dz}{iz} \end{align}

Now deform the contour along the essential singularity $z=0$ to obtain

\begin{align} \int_{|z|=1}\cosh\left(\bar{\omega}z+\omega z^{-1}\right)\frac{dz}{iz} &=\lim_{\epsilon\to 0^+}\int^{2\pi}_0\cosh\left(\bar{\omega}\epsilon e^{i\theta}+\omega\epsilon^{-1}e^{-i\theta}\right)\ d\theta\\ &=\lim_{\epsilon\to 0^+}\int^{2\pi}_0\sum^\infty_{k=0}\frac{1}{(2k)!}\sum^{2k}_{j=0}\binom{2k}{j}\bar{\omega}^{2k-j}\omega^j(\epsilon e^{i\theta})^{2k-2j}\ d\theta\\ &=\lim_{\epsilon\to 0^+}\sum^\infty_{k=0}\frac{1}{(2k)!}\sum^{2k}_{j=0}\binom{2k}{j}\bar{\omega}^{2k-j}\omega^j\int^{2\pi}_0(\epsilon e^{i\theta})^{2k-2j}\ d\theta\\ &=\lim_{\epsilon\to 0^+}\sum^\infty_{k=0}\frac{1}{(2k)!}\sum^{2k}_{j=0}\binom{2k}{j}\bar{\omega}^{2k-j}\omega^j(2\pi\mathbf{1}_{\{k\}}(j))\\ &=\sum^\infty_{k=0}\frac{1}{(2k)!}\binom{2k}{k}|\omega|^{2k}\cdot 2\pi\\ &=\sum^\infty_{k=0}\frac{\left(|\omega|^2\right)^k}{k!\Gamma(k+1)}\cdot 2\pi\\ &=2\pi I_0(2|\omega|) \end{align}

Repeating the same argument, we clearly also have

$$\int_{|z|=1}\cosh\left({\omega}z+\bar{\omega} z^{-1}\right)\frac{dz}{iz}=\int_{|z|=1}\cosh\left(\bar{\omega}z+\omega z^{-1}\right)\frac{dz}{iz}=2\pi I_0(2|\omega|)$$

and therefore,

$$g_1(x_2, \mu^{(2)})=\pi I_0(2|\omega|)=\pi I_0\left(2.\sqrt{\frac{a^2+b^2}{4}}\right)=\pi I_0\left(\sqrt{x_2(\mu_1^2+\mu_2^2)}\right)=\pi I_0\left(\sqrt{\delta_2 x_2}\right)$$

Brief Derivation of Equation (1)

In case anyone is interested, this is the sequence of transformations I used to arrive at equation $(1)$. Let $X_k$, $k=1,\cdots,n$, be independent random variables each distributed as $N(\mu_k,1)$. First, it is trivial to show that each $X_k^2$ has Lebesgue pdf $\frac{1}{\sqrt{2\pi x}}e^{-(x+\mu_k^2)/2}\cosh\left(\mu_k\sqrt{x}\right)\mathbf{1}_{(0,\infty)}(x)$. Integrating the joint density of $(X_1^2,\cdots,X_{n-1}^2, X_1^2+\cdots+X_n^2)$ gives us, for $x_n>0$,

\begin{align} \small{f_n(x_n,\delta_n)} &\small{=\frac{e^{-(x_n+\delta_n)/2}}{(2\pi)^{n/2}}\int_{\mathbb{R}^{n-1}}\frac{\prod^{n-1}_{k=1}\cosh\left(\mu_k\sqrt{x_k-S_{n-1}\mathbf{1}_{\{n\}}(k)}\right)\mathbf{1}_{(0,\infty)}(x_k-S_{n-1}\mathbf{1}_{\{n\}}(k))}{\sqrt{x_1\cdots x_{n-1}(x_n-S_{n-1})}}\ dx_1\cdots dx_{n-1}}\\ &\small{=\frac{e^{-(x_n+\delta_n)/2}}{(2\pi)^{n/2}}\int^{x_n}_0\int^{x_{n-1}}_0\cdots\int^{x_2}_0\prod^{n}_{k=1}\frac{\cosh\left(\mu_k\sqrt{x_k-x_{k-1}}\right)}{\sqrt{x_k-x_{k-1}}}}\ dx_1\ dx_2\cdots dx_{n-1}\tag{i}\\ &\small{=\frac{e^{(-\delta_n+x_n)/2}}{(2\pi)^{n/2}}x_n^{n/2-1}g_{n-1}(x_n,\mu^{(n)})}\tag{ii} \end{align}

where $x_0:=0$ and $S_{n-1}:=\sum^{n-1}_{k=1}x_k$.

Explanation:

$\text{(i)}$: Apply the transformation $(x_1,x_2,\cdots,x_{n-1})\mapsto (x_1, x_2-x_1,\cdots,x_{n-1}-x_{n-2})$. The determinant of the Jacobian for this transformation is $1$.

$\text{(ii)}$: Apply the transformation $(x_1,x_2,\cdots,x_{n-1})\mapsto \left(\prod^n_{k=1}x_k, \prod^n_{k=2}x_k,\cdots, \prod^n_{k=n-1}x_k\right)$. The Jacobian for this transformation is an upper triangular matrix with diagonal elements $\prod^n_{k=2}x_k, \prod^n_{k=3}x_k,\cdots, \prod^n_{k=n}x_k$ and its determinant is thus $\prod^n_{k=2}x_k^{k-1}$. Further simplification gives us the stated result.