On the Wikipedia page here , it states that the Green's function for 3D relativistic heat conduction (with $c=1$)

$$[\partial_t^2 + 2\gamma\partial_t -\Delta_{3D}] u(t,x) = \delta(t,x) = \delta(t)\delta(x)$$

is given by $$u(t,x) = \frac{e^{-\gamma t}}{20\pi}\bigg( \big[8 - 3e^{-\gamma t} + 2\gamma t + 4\gamma^2t^2\big]\frac{\delta(t-|x|)}{|x|^2} + \gamma^2\Theta(t - |x|)\big[ \frac{I_1\big(\gamma\sqrt{t^2 - |x|^2}\big)}{\sqrt{t^2 - |x|^2}} + t\frac{I_2\big(\gamma\sqrt{t^2 - |x|^2}\big)}{t^2 - |x|^2}\big] \bigg)$$

Where $I_1$ and $I_2$ are modified Bessel functions of the 1st and 2nd kind. Furthermore, $\Theta$ is the Heaviside step function . I am looking for a derivation of this or a reference to one.

**Soft Attempt:**
Taking the Fourier transform with respect to $x$ we obtain
$$[\partial_t^2 + 2\gamma\partial_t + (4\pi^2|\xi|^2)] \hat{u}(t,\xi) = \delta(t)$$

On the same Wikipedia page, this is a 1D damped harmonic oscillator and they claim the Green's function is given by $$\hat{u}(t,\xi) = \Theta(t)e^{-\gamma t} \frac{\sin\big(t\sqrt{4\pi^2|\xi|^2 - \gamma^2}\big)}{\sqrt{4\pi^2|\xi|^2 - \gamma^2}}$$ So I should expect that

$$u(t,x) = \int e^{2\pi i \xi\cdot x}\Theta(t) e^{-\gamma t}\frac{\sin\big(t\sqrt{4\pi^2|\xi|^2 - \gamma^2}\big)}{\sqrt{4\pi^2|\xi|^2 - \gamma^2}}d\xi$$

But I am not sure how we would get the formula above or if this is even the correct approach.