In 4-dimensional Euclidean QFT, the propagator can be expressed in terms of a modified Bessel function of the second kind:
\begin{align} \Delta_E (x_1-x_2)&\equiv\int\frac{d^4p}{(2\pi)^4}\frac{e^{ip\cdot(x_1-x_2)}}{p^2+m^2}\\ &=\frac{m}{(2\pi)^2|x_1-x_2|}K_1(m|x_1-x_2|), \tag{*}\label{B} \end{align} where $p\cdot x\equiv\sum^{4}_{i=1} p^i x^i$ and $|x_1-x_2|\equiv \sqrt{\sum^{4}_{i=1}(x_1^i-x_2^i)^2}$.
In the context of some manipulations I am doing, it is actually more convenient to work with the 4-dimensional integral representation of this Bessel function. However, I also need to manipulate the modified Bessel function of the first kind $I_1 (z)$, and I don't know how to put this into a similar 4-dimensional integral form.
I am familiar with how to, e.g., derive \eqref{B} by first reducing (via a contour integral, etc.) the 4-dimensional integral in the first line into a 1-dimensional integral like
\begin{align} \int\frac{d^4p}{(2\pi)^4}\frac{e^{ip\cdot(x_1-x_2)}}{p^2+m^2}=\frac{1}{(2\pi)^2}\int^\infty_m d\omega_\vec{p} \sqrt{\omega^2_\vec{p}-m^2}e^{-\omega_\vec{p}|x_1-x_2|},\\ \end{align} where $\omega^2_\vec{p}\equiv \vec{p}^2+m^2$ with $\vec{p}\equiv (p^1,p^2,p^3).$ Rescaling the integration variable so that the bounds of the integral run instead from $1\to\infty$ will then put this into a well-known representation of the modified Bessel function of the second kind.
Therefore, it seems that one should be able to work in reverse: i.e., start with \begin{equation} \frac{m}{|x_1-x_2|}I_1(m|x_1-x_2|), \end{equation} and use a 1-dimensional representation of the modified Bessel function of the first kind to put this into the form of a 4-dimensional integral.
Edit: As a first step, we have the 1-dimensional integral representation,
\begin{align} I_1[z]&=\frac{z}{\pi}\int^1_{-1}\sqrt{1-t^2} e^{\pm zt} dt\\ &=\frac{z}{\pi m^2}\int^m_{-m}\sqrt{m^2-t^2} e^{\pm zt/m} dt, \end{align} where in going to the second line I rescaled the integration variable $t\to t/m$. Now, this implies,
\begin{align} \frac{m}{|x_1-x_2|}I_1(m|x_1-x_2|)&=\frac{1}{\pi}\int^m_{-m} d\omega\sqrt{m^2-\omega^2} e^{\pm \omega|x_1-x_2|} \\ &=\frac{i}{\pi}\int^m_{-m} d\omega\sqrt{\omega^2-m^2} e^{\pm \omega|x_1-x_2|} \tag{**}\label{I}, \end{align} where at this stage "$\omega$" is just an integration variable with no specified relation to coordinates in higher dimensions like $\vec{p}$.
Choosing the negative exponential in \eqref{I} will yield an integrand of the same form as \eqref{B}. However, I'm unsure how to convert this 1-dimensional integral into a 4-dimensional integral in this case (i.e., with integration bounds $-m\to m$ rather than $m\to\infty$).