I'm trying to find the explicit solution of the wave equation in three spacial dimensions with the initial condition: $$ u(x,0) = \begin{cases} \sqrt{1-x^2}, & \text{if } |x| \leq 1, \\ 0, & \text{if } |x| \geq 1, \end{cases} $$ and $$ \partial_t u(x,0) = 0,$$ where $x = (x_1,x_2,x_3)$. I guess that using spherical coordinates in $\mathbb{R}^3$ is the wisest choice, and that I should make use of Kirchhoff's formula, $$ u(x,t) = \int_{\partial B(x,t)} t h(y) + g(y) + \nabla g(y) \cdot (y-x) dS(y), $$ where $$g = u(x,0), \quad h = \partial_t u(x,0),$$ and the integral sign should be read as the mean integral (I tried typing \fint but I had no luck).
I'm having a bit of trouble making the change of variables from rectangular to spherical coordinates, and how that would change the integration process...
We're looking for solutions $u(\mathbf{r},t)$ to the wave equation $$\frac{\partial^2 u}{\partial t^2}=c^2\nabla^2u$$ Which, in spherical coordinates, is $$\frac{c^2}{r^2\sin(\phi)}\left(\frac{\partial }{\partial r}\left(r^2\sin(\phi)\frac{\partial u}{\partial r}\right)+\frac{\partial }{\partial \theta}\left(\frac{1}{\sin(\phi)}\frac{\partial u}{\partial \theta}\right)+\frac{\partial }{\partial \phi}\left(\sin(\phi)\frac{\partial u}{\partial \phi}\right)\right)=\frac{\partial^2u}{\partial t^2}$$ Using the $(r,\theta,\phi)$ (radius, azimuthal angle, polar angle) convention. First step: Assume the solution is separable, i.e, $u(\mathbf{r},t)=R(r)\Theta(\theta)\Phi(\phi)T(t)$. Now we plug this into our wave equation. $$R\Theta\Phi T''=\frac{c^2}{r^2\sin\phi}\left(\frac{\partial }{\partial r}\left(r^2\sin(\phi) R'\Theta\Phi T\right)+\frac{\partial }{\partial \theta}\left(\frac{1}{\sin\phi}R\Theta'\Phi T\right)+\frac{\partial }{\partial \phi}\left(\sin(\phi) R\Theta\Phi'T\right)\right)$$ Simplifying, $$R\Theta\Phi T''=c^2\left(\frac{\Theta\Phi T}{r^2}\frac{\partial }{\partial r}\left(r^2R'\right)+\frac{R\Phi T}{r^2\sin^2\phi}\frac{\partial }{\partial \theta}\left(\Theta'\right)+\frac{R\Theta T}{r^2\sin\phi}\frac{\partial }{\partial \phi}\left(\sin(\phi) \Phi'\right)\right)$$ Dividing out, $$ \frac{T''}{T}=\frac{c^2}{r^2}\left(\frac{1 }{R}\frac{\partial }{\partial r}\left(r^2R'\right)+\frac{1 }{\Theta \sin^2\phi}\frac{\partial }{\partial \theta}\left(\Theta'\right)+\frac{1 }{\Phi \sin\phi}\frac{\partial }{\partial \phi}\left(\sin(\phi) \Phi'\right)\right)$$
Now I've illustrated all this so far since this is how you'd approach this problem in general. However, since the initial conditions you've provided are radially symmetric, i.e, independent of $\theta$ and $\phi$, we can assume $\Theta'=\Phi'=0.$ Therefore the above simplifies to $$\frac{T''}{T}=\frac{c^2}{r^2}\frac{1 }{R}\frac{\partial }{\partial r}\left(r^2R'\right)$$ We now assume our solution is "pure frequency", i.e, $$\frac{T''}{T}=-\omega^2 \implies T(t)=a_1\cos(\omega t)+a_2\sin(\omega t)$$ Substituting into the previous, $$-\omega^2=\frac{c^2}{r^2}\frac{1 }{R}\frac{\partial }{\partial r}\left(r^2R'\right)$$
Defining $k=\frac{\omega}{c}$, $$\frac{\partial }{\partial r}\left(r^2R'\right)+k^2r^2R=0$$ This is an Emden-Fowler differential equation, and has the solution $$R(r)=\frac{1}{r}\left(b_1e^{-ikr}+\frac{ib_2}{2k}e^{ikr}\right)$$
Let $b_1=p+iq$, $b_2=p'+iq'$. Then, $$r\cdot R(r)=(p+iq)e^{-ikr}+\frac{i(p'+iq')}{2k}e^{ikr}$$ $$r\cdot R(r)=(p+iq)e^{-ikr}+\frac{1}{2k}(-q'+ip')e^{ikr}$$ Using Euler's formula, $$r\cdot R(r)=p\cos(kr)-ip\sin(kr)+iq\cos(kr)+q\sin(kr)+\frac{1}{2k}(-q'\cos(kr)+iq'\sin(kr)+ip'\cos(kr)-p'\sin(kr))$$
Grouping the cosine and sine terms together, $$r\cdot R(r)=\cos(kr)\left(p+iq-\frac{q'}{2k}+\frac{ip'}{2k}\right)+\sin(kr)\left(-ip+q+\frac{iq'}{2k}-\frac{p'}{2k}\right)$$ In order for $R(r)$ to be real-valued, we require $$q+\frac{p'}{2k}=0 ~\text{ and } -p+\frac{q'}{2k}=0$$ Thus, $$\begin{bmatrix} p'\\ q' \end{bmatrix} =2k\begin{bmatrix} 0 & -1\\ 1 & 0 \end{bmatrix}\begin{bmatrix} p\\ q \end{bmatrix}$$ Therefore given $b_1=p+iq$, we require $b_2=2k(-q+ip)$. Thus, $$R(r)=\frac{1}{r}\left(\cos(kr)(p-\frac{2kp}{2k})+\sin(kr)(q-\frac{-2kq}{2k})\right)$$
Finally, substituting $k=\frac{\omega}{c}$ back in and simplifying, $$R(r)=2q\frac{\omega}{c}\operatorname{sinc}\left(\frac{\omega}{c}r\right)$$ With $\operatorname{sinc}(x):=(1/x)\sin(x)$. Side note: $\operatorname{sinc}$ is the zeroth spherical Bessel function of the first kind, knowledge of which will be needed both for the more general PDE involving $\theta$ and $\phi$, and dealing with the initial conditions later on.
Now, our solution $u(\mathbf{r},t)$(which we can restate as simply $u(r,t)$ because of no $\theta$ or $\phi$ dependence) for this radially symmetric wave equation is $$u(r,t)=\left(a_1\cos(\omega t)+a_2\sin(\omega t)\right)2q\frac{\omega}{c}\operatorname{sinc}\left(\frac{\omega}{c}r\right)$$ Or, more nicely, $$u(r,t)=C\frac{\omega}{c}\operatorname{sinc}\left(\frac{\omega}{c}r\right)(A\cos(\omega t)+B\sin(\omega t))$$ Now let's examine our initial conditions. We want $\frac{\partial u}{\partial t}(r,0)=0$. $$\frac{\partial u}{\partial t}=C\frac{\omega}{c}\operatorname{sinc}\left(\frac{\omega}{c}r\right)(-A\sin(\omega t)+B\cos(\omega t))$$ $$\frac{\partial u}{\partial t}\bigg|_{(r,0)}=C\frac{\omega}{c}\operatorname{sinc}\left(\frac{\omega}{c}r\right)B$$
In order for this to be non trivially $=0 ~ \forall r$, we need $B=0$. As both $A$ and $C$ are constants we can rename $A\cdot C \to A$, therefore $$u(r,t)=A\frac{\omega}{c}\operatorname{sinc}\left(\frac{\omega}{c}r\right)\cos(\omega t)$$ Before we deal with the other initial condition, we need to generalize our solution a bit. Recall that we can state $u(r,t)$ in terms of spherical Bessel functions: $$u(r,t)=A\frac{\omega}{c}\cos(\omega t)j_0\left(\frac{\omega}{c}r\right)$$ Or, better, in terms of ordinary Bessel functions using the identities here: $$u(r,t)=A\frac{\omega}{c}\cos(\omega t)\sqrt{\frac{\pi}{2r}}J_{1/2}\left(\frac{\omega}{c}r\right)$$
Now, one thing we can notice is that the wave equation is linear, i.e if $u_1$ and $u_2$ are both solutions, $a\cdot u_1 +b\cdot u_2$ is also a solution, given real numbers $a$ and $b$. Therefore we can restate our solution more generally as a linear combination of pure frequency solutions (also known as normal modes), $$u(r,t)=\frac{A}{c}\sqrt{\frac{\pi}{2r}}\sum_{n=0}^\infty a_n\cos(\omega_nt)\omega_n J_{1/2}\left(\frac{\omega_n}{c}r\right)$$ Where $a_1,a_2,...$ is any sequence of real numbers (provided the sum converges), and $\omega_n$ is a factor that ensures we have roots at $r=\pm 1$. Specifically, let $\alpha_{1/2,n}$ be the $n$th positive root of $J_{1/2}$. Then $\omega_n =c \alpha_{1/2,n}$. Also let $A_n=A a_n\alpha_n$. We can now state $u(r,0)$ as $$u(r,0)=\sqrt{\frac{\pi}{2r}}\sum_{n=0}^\infty A_n J_{1/2}\left(\alpha_{1/2,n} r\right)$$ Recall that in this example, our initial conditions are $$ u(r,0)=u_0(r) = \begin{cases} \sqrt{1-r^2}, & \text{if } |r| \leq 1, \\ 0 & \text{if } |r| \geq 1 \end{cases} $$ So now the challenge is to find a sequence of real numbers $A_1,A_2,...$ such that $$u_0(r)\sqrt{\frac{2r}{\pi}}=\sum_{n=0}^\infty A_n J_{1/2}\left(\alpha_{1/2,n} r\right)$$ For which you'll need Fourier-Bessel series and/or Hankel Transforms. In this case it turns out that $$A_n=\frac{2}{J_{3/2}(\alpha_{1/2,n})^2}\int_0^1 r\sqrt{\frac{2r}{\pi}}\sqrt{1-r^2} J_{1/2}(\alpha_{1/2,n}r)\mathrm{d}r$$ If you can't be bothered with that, you can use the very rough approximation $$\sqrt{1-x^2}\approx \frac{\sin(2x)}{2x}=j_0(2x)=\sqrt{\frac{\pi}{2x}}J_{1/2}(2x).$$ I know this is not a complete answer, but I hope it's a good start!