Solution of a differentiation in integral form

First of all, taking the variables $$ t = \frac{t_0}{m}, \qquad r = \frac{r_0}{m}, $$ the PDE becomes

$$ \phi_{tt} - \frac{1}{r^{d-1}}\frac{\partial}{\partial r} \left(r^{d-1} \frac{\partial \phi}{\partial r}\right) + 2 \phi = 0 $$ where we have lost the sub-indexes in the notation.

Now, using separation of variables $\phi(t,r) = A(t) \Phi(r)$, we end up with $$ A''(t) \Phi(r) - \frac{A(t)}{r^{d-1}}\frac{d}{d r} \left(r^{d-1} \frac{d \Phi(r)}{d r}\right) + 2 A(t)\Phi(r) = 0. $$ Dividing by $A(t)\Phi(r)$ and rearranging, $$ \frac{A''(t)}{A(t)} + 2 = \frac{1}{r^{d-1}\Phi(r)}\frac{d}{d r} \left(r^{d-1} \frac{d \Phi(r)}{d r}\right). $$ The left hand size is a function of $t$, while the right hand is a function of $r$. The only way for the equality to hold is that both functions are equal to a constant (which I'll call -$k^2$), hence \begin{align} A''(t) + (2 + k^2)A(t) &= 0,\\ \\ \frac{1}{r^{d-1}}\frac{d}{d r} \left(r^{d-1} \frac{d \Phi(r)}{d r}\right) + k^2 \Phi(r) &=0. \end{align} The first equation is easily solved, yielding $$ A(t) = A_0 \cos \omega t + B_0 \sin \omega t $$ where $\omega = \sqrt{2 + k^2}$ is called the dispersion relation. For the radial equation, taking $$ \Phi(r) = r^\beta \chi(r), $$ we have $$ \chi''(r) + \frac{d - 1 + 2\beta}{r}\chi'(r) + \left(k^2 + \frac{\beta(\beta + d - 2)}{r^2}\right)\chi(r) = 0. $$ Let $\beta = 1 - \frac{d}{2}$, $$ \chi''(r) + \frac{1}{r} \chi'(r) + \left(k^2 - \frac{(1-\frac{d}{2})^2}{r^2}\right)\chi(r) = 0. $$ Taking $r = \frac{\rho}{k}$ and $\chi(r) = \hat{\chi}(\rho)$, $$ \hat{\chi}''(\rho) + \frac{1}{\rho} \hat{\chi}'(\rho) + \left(1 - \frac{(1-\frac{d}{2})^2}{\rho^2}\right)\hat{\chi}(\rho) = 0, $$ which is the Bessel differential equation of order $1-\frac{d}{2}$, and has as solution $$ \hat{\chi}(\rho) = c_1 J_{1-\frac{d}{2}}(\rho) + c_2 Y_{-1+\frac{d}{2}}(\rho), $$ where $J_\alpha$ and $Y_\alpha$ are the Bessel function of first and second kind, of order $\alpha$.

The spatial solution for every $k$ is $$ \Phi_k(r) = r^{1-\frac{d}{2}} \left(c_1 J_{1-\frac{d}{2}}(k r) + c_2 Y_{-1+\frac{d}{2}}(kr)\right). $$ In the case $d = 3$, the solutions can be written as $$ \Phi_k(r) = \hat{c}_1 j_0(k r) + \hat{c}_2 y_0(k r) = \hat{c}_1(k) \frac{\sin k r}{r} + \hat{c}_2(k) \frac{\cos k r}{r}, $$ where $j_0$ and $y_0$ are the Spherical Bessel Functions.

To determine the constants $\hat{c}_1(k)$ and $\hat{c}_2(k)$, first we note that $r = 0$ is part of the domain of the field, and $\Phi_k(r)$ must be well behaved at the origin. This means that $\hat{c}_2(k) = 0$. On the other hand, there is a perfectly valid solution for every $k >0$, then $$ \phi(t,r) = \int_0^\infty \hat{c}_1(k) \frac{\sin k r}{r} [A_0(k) \cos \omega t + B_0(k) \sin \omega t] d k. $$ To determine $b(k)$, $A_0(k)$ and $B_0(k)$, we use the conditions \begin{align} \phi(0,r) &= a_0 e^{-r^2/R^2}, \\ \phi_r(0,r) &= 0, \\ \phi_t(r,0) &=0, \end{align} and then $$ a_0 e^{-r^2/R^2} = \int_0^\infty \hat{c}_1(k) \frac{\sin(k r)}{r} dk, $$ where $\hat{c}_1(k)$ can be obtained using the Fourier Transform, yielding $$ \hat{c}_1(k) = \frac{2 a_0}{\pi} \Im\left\{\int_0^\infty r e^{-r^2/R^2 + ikr}dr\right\}. $$ This can be easily integrated, and then the solution is $$ \phi(t,r) = \frac{a_0 R^3}{2 \sqrt{\pi}} \int_0^\infty k e^{-R^2 k ^2/4} \frac{\sin k r}{r} \cos \omega t \,dk. $$