Solution (source code)

= Solution

Split the inverse transform into the intervals $|\omega|<\mu$ and $|\omega|>\mu$, and in each case pair the positive-frequency <integral> with its <complex conjugate> at negative <frequency>. The low-frequency pair becomes
$$
\frac1\pi\operatorname{Re}\int_0^\mu A(\omega)e^{ik(\omega)r}e^{i\omega t}\,d\omega,
$$
since $B=0$ and the decaying spatial factor is real. In the high-frequency positive half, replace $A(\omega)$ by $B(-\omega)^*$. The negative half is again its conjugate. Combining the two pieces gives
$$
\boxed{\sigma(t,r)=\frac1\pi\operatorname{Re}\left[
\int_0^\mu A(\omega)e^{ikr}e^{i\omega t}\,d\omega+
\int_\mu^\infty\left\{B(\omega)e^{-ikr}e^{i\omega t}
+B(-\omega)^*e^{ikr}e^{i\omega t}\right\}d\omega\right],\quad\lambda=\frac1\pi.}
$$
Every integrand satisfies the massive <wave equation> because $\omega^2=k^2+\mu^2$. Linearity, and taking the real part, preserve the equation whenever differentiation under the <integrals> is justified, or distributionally under the corresponding spectral assumptions. The construction is exactly the <real bounded Fourier modes of the Klein-Gordon equation> decomposition.

The displayed expression uses ordinary spectral amplitudes. If distributions supported at $\omega=\pm\mu$ are admitted, a further bounded threshold contribution $c_1\cos\mu t+c_2\sin\mu t$, independent of $r$, is possible. It is not a missing propagating or evanescent branch of the open <frequency> intervals.