I've recently come across an integral of the form $$\int_{0}^{\infty} \sinh(\eta)^2\int_{-1}^{1}f_1(a\cosh(\eta)+b\sinh(\eta)x)f_2(a\cosh(\eta)-b\sinh(\eta)x)dxd\eta$$ Where $a>b$, and $f_1$ and $f_2$ are some real valued functions that decay very rapidly (safe to assume $f_i(x) \sim O(e^{-\alpha_i x})$ for $\alpha_i>0$). I am trying to compute this in terms of their Fourier transforms $\hat{f_i}$, for which I use: $$f_i(x)=\int_{-\infty}^{\infty}\hat{f_i}(k)e^{ikx}dk$$

I also am only interested in $f_i(x)$ for $x\geq 0$, so we can assume for the sake of computing the Fourier transform $f_i(-x)=f_i(x)$

Plugging this in to the above, we get:

$$\int_{0}^{\infty} \sinh(\eta)^2\int_{-1}^{1}\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}\hat{f_1}(k_1)\hat{f_2}(k_2)e^{ia\cosh(\eta)(k_1+k_2)}e^{ib\sinh(\eta)(k_1-k_2)x}dk_1dk_2dxd\eta$$

I want to reduce this to a double integral, and have come very close but am encountering some issues. I will show why I believe we can reduce this to a double integral--the following is not rigorous in any way. Some of my steps are unjustified, with the hope that the resulting integral will converge to the same value.

Changing the order of integration (unjustified):

$$\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}\hat{f_1}(k_1)\hat{f_2}(k_2)\int_{0}^{\infty} \sinh(\eta)^2e^{ia\cosh(\eta)(k_1+k_2)}\int_{-1}^{1}e^{ib\sinh(\eta)(k_1-k_2)x}dxd\eta dk_1dk_2$$

The inntermost integral is easy to compute: $$\int_{-1}^{1}e^{ib\sinh(\eta)(k_1-k_2)x}dx=\frac{e^{ib\sinh(\eta)(k_1-k_2)}-e^{-ib\sinh(\eta)(k_1-k_2)}}{ib\sinh(\eta)(k_1-k_2)}$$

This allows us to complete the integral over $\eta$ as well: $$\int_{0}^{\infty} \sinh(\eta)^2e^{ia\cosh(\eta)(k_1+k_2)}\int_{-1}^{1}e^{ib\sinh(\eta)(k_1-k_2)x}dxd\eta =\\ \frac{1}{ib(k_1-k_2)}\int_{0}^{\infty} \sinh(\eta)e^{ia\cosh(\eta)(k_1+k_2)}[e^{ib\sinh(\eta)(k_1-k_2)}-e^{-ib\sinh(\eta)(k_1-k_2)}]d\eta= \\ \frac{1}{ib(k_1-k_2)}\int_{-\infty}^{\infty} \sinh(\eta)e^{ia\cosh(\eta)(k_1+k_2)}e^{ib\sinh(\eta)(k_1-k_2)}d\eta =\\ \frac{-1}{b(k_1-k_2)^2} \frac{\partial}{\partial b}\int_{-\infty}^{\infty} e^{ia\cosh(\eta)(k_1+k_2)}e^{ib\sinh(\eta)(k_1-k_2)}d\eta$$ To my luck, I was able to find a remarkable identity to compute this: $$H_0^{1}(\sqrt{a^2-b^2})=\frac{1}{i\pi}\int_{-\infty}^{\infty} e^{ia\cosh(\eta)}e^{ib\sinh(\eta)}d\eta$$ $$\frac{\partial}{\partial b}H_0^{1}(\sqrt{a^2-b^2})=\frac{b}{\sqrt{a^2-b^2}}H_1^{1}(\sqrt{a^2-b^2})$$ Where $H_{\nu}^{1}$ is the Hankel function of the first type.

Unfortunately, the identity only holds for $\Im(a \pm b)>0$. However, $H_{\nu}^{1}(\sqrt{a^2-b^2})=J_{\nu}(\sqrt{a^2-b^2})+iY_{\nu}(\sqrt{a^2-b^2})$ is still defined for $a>b,a,b\in \mathbb{R}$. So, let me just pretend that we are justified to use this here. Using this formula, we observe:

$$\frac{-1}{b(k_1-k_2)^2} \frac{\partial}{\partial b}\int_{-\infty}^{\infty} e^{ia\cosh(\eta)(k_1+k_2)}e^{ib\sinh(\eta)(k_1-k_2)}d\eta= \\ \frac{-i\pi}{b(k_1-k_2)^2} \frac{\partial}{\partial b}H_0^{1}(\sqrt{a^2(k_1+k_2)^2-b^2(k_1-k_2)^2})= \\ \frac{-i\pi}{\sqrt{a^2(k_1+k_2)^2-b^2(k_1-k_2)^2}} H_1^{1}(\sqrt{a^2(k_1+k_2)^2-b^2(k_1-k_2)^2})$$

Finally, plugging back into the initial expression, we obtain:

$$\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}\hat{f_1}(k_1)\hat{f_2}(k_2) \frac{-i\pi}{\sqrt{a^2(k_1+k_2)^2-b^2(k_1-k_2)^2}} H_1^{1}(\sqrt{a^2(k_1+k_2)^2-b^2(k_1-k_2)^2}) dk_1dk_2$$

The issue with this is the argument of $H_1^{1}$ can be negative, landing us on a branch cut, and numerical tests show this result does not hold. At this point, I am unsure how to deal with this.

I would greatly appreciate if anyone with more knowledge than me could guide me on how to procede/modify my efforts, or tell me if this route likely will not work.