$\renewcommand{\Re}{\operatorname{Re}}$ $\renewcommand{\Im}{\operatorname{Im}}$ $\newcommand{\erf}{\operatorname{erf}}$ $\newcommand{\dag}{\dagger}$ $\newcommand{\const}{\mathrm{const}}$ $\newcommand{\arcsinh}{\operatorname{arcsinh}}$
Consider Cauchy problem for $3$-dimensional wave equation \begin{align} & u_{tt}-c^2\Delta u=f,\label{eq-9.1.1}\\[3pt] & u|_{t=0}=g,\label{eq-9.1.2}\\[3pt] & u_t|_{t=0}=h.\label{eq-9.1.3} \end{align}
Assume first that $f=g=0$. We claim that in this case as $t>0$ \begin{equation} u(x,t)= \frac{1}{4\pi c ^2 t} \iint _{S(x,ct)} h(y)\,d\sigma \label{eq-9.1.4} \end{equation} where we integrate along sphere $S(x,ct)$ with a center at $x$ and radius $ct$; $d\sigma$ is an area element.
Let us prove (\ref{eq-9.1.4}) first as $h(x)=e^{ix\cdot \xi}$ with $\xi\in \mathbb{R}^3\setminus 0$; we use the standard notation $x\cdot \xi=x_1\xi_1+x_2\xi_2+x_3\xi_3$. In this case \begin{equation} u(x,t)=e^{ix\cdot \xi}c^{-1}|\xi|^{-1}\sin (ct|\xi|) \label{eq-9.1.5} \end{equation} is obviously a solution to Cauchy problem (\ref{eq-9.1.1})--(\ref{eq-9.1.3}).
On the other hand, the right-hand expression of (\ref{eq-9.1.4}) becomes \begin{equation*} \frac{1}{4\pi c ^2 t}\iint _{S(x,ct)} e^{iy\cdot \xi}\,d\sigma= \frac{1}{4\pi c ^2 t}e^{ix\cdot \xi}\iint _{S(0,ct)} e^{iz\cdot \xi}\,d\sigma \end{equation*} where we changed variable $y=x+z$ with $z$ running $S(0,ct)$ (sphere with a center at $0$) and we need to calculate integral in the right-hand expression. Let us select coordinate system in which $\xi=(0,0,\omega)$ with $\omega =|\xi|$ and introduce corresponding spherical coordinates $(\rho,\phi,\theta)$; then on $S(0,ct)$ $\rho=ct$, $z\cdot \xi= z_3\omega= ct \omega \cos(\phi)$ and $d\sigma=c^2t^2 \sin(\phi)d\phi d\theta$; so integral becomes \begin{multline*} c^2t^2 \int_0^\pi e^{i ct \omega \cos(\phi)}\sin(\phi) \,d\phi \int_0^{2\pi} d\theta=\\ -2\pi c^2t^2 \int_0^\pi e^{i ct \omega \cos(\phi)} \,d\cos(\phi) = 2\pi c^2t^2 \frac{1}{ic t\omega} \Bigl( e^{i ct \omega }- e^{-i ct \omega } \Bigr)=\\ 4\pi ct \omega^{-1} \sin (ct\omega) \end{multline*} and multiplying by $e^{ix\cdot\xi}$ and dividing by $4\pi c^2t$ we get $e^{ix\cdot\xi} c^{-1}|\xi|^{-1} \sin (ct|\xi|)$ which is the right-hand expression of (\ref{eq-9.1.5}).
So, for $h(x)=e^{ix\cdot \xi}$ (\ref{eq-9.1.4}) has been proven. However the general function $h(x)$ could be decomposed into such special functions using multidimensional Fourier transform and multidimensional Fourier integral which is nothing but repeated $1$-dimensional Fourier transform and Fourier integral: \begin{gather} h(x)=\iiint \hat{h}(\xi) e^{ix\cdot\xi}\,d\xi,\label{eq-9.1.6}\\[3pt] \hat{h}(\xi)=(2\pi)^{-n}\iiint h(x) e^{-ix\cdot\xi}\,dx\label{eq-9.1.7} \end{gather} and therefore (\ref{eq-9.1.4}) extends to general functions as well.
Remark 1. We should deal with the fact that only decaying functions could be decomposed into Fourier integral, but this is easy due to the fact that integral in (\ref{eq-9.1.4}) is taken over bounded domain.
To cover $t<0$ we replace (\ref{eq-9.1.4}) by \begin{equation} u(x,t)= \frac{1}{4\pi c t} \iint _{S(x,c|t|)} h(y)\,d\sigma \label{eq-9.1.8} \end{equation} which is obvious as $u$ must be odd with respect to $t$.
Consider now $g\ne 0$. Let $v$ be given by (\ref{eq-9.1.8}) with $h$ replaced by $h$; then \begin{align*} & v_{tt}-c^2\Delta v=0,\\[3pt] & v|_{t=0}=0,\\[3pt] & v_t|_{t=0}=g. \end{align*} Then $\Delta v|_{t=0}=0$ and therefore $v_{tt}|_{t=0}=0$ and therefore differentiating equation with respect to $t$ we conclude that $u:=v_t$ solves \begin{align*} & u_{tt}-c^2\Delta u=0,\\[3pt] & u|_{t=0}=g,\\[3pt] & u_t|_{t=0}=0. \end{align*} Now \begin{equation} u(x,t)= \frac{\partial\ }{\partial t} \Bigl(\frac{1}{4\pi c^2 t} \iint _{S(x,c|t|)} g(y)\,d\sigma\Bigr). \label{eq-9.1.9} \end{equation} Therefore solving separately (\ref{eq-9.1.1})--(\ref{eq-9.1.3}) for $f=g=0$ (given by (\ref{eq-9.1.8})) and for $f=h=0$ (given by (\ref{eq-9.1.9}) and adding solutions due to linearity we arrive to \begin{equation} u(x,t)= \frac{\partial\ }{\partial t} \Bigl(\frac{1}{4\pi c^2 t} \iint_{S(x,c|t|)} g(y)\,d\sigma\Bigr)+ \frac{1}{4\pi c ^2 t} \iint _{S(x,c|t|)} h(y)\,d\sigma \label{eq-9.1.10} \end{equation} covering case $f=0$.
To cover case of arbitrary $f$ but $g=h=0$ we apply Duhanel integral formula (see Subsection 2.5.2). Consider problem \begin{align*} & U_{tt}-c^2\Delta U=0,\\[3pt] & U|_{t=\tau}=0,\\[3pt] & U_t|_{t=\tau}=f(x,\tau). \end{align*} Its solution is given by \begin{equation*} U(x,t,\tau) = \frac{1}{4\pi c^2(t-\tau)} \iint_{S(x,c|t-\tau|) } f(y,\tau)\,d\sigma \end{equation*} and therefore \begin{equation} u(x,t) = \int_0^t\frac{1}{4\pi c^2(t-\tau)} \iint_{S(x,c|t-\tau|)} f(y,\tau)\,d\sigma d\tau. \label{eq-9.1.11} \end{equation} Assembling (\ref{eq-9.1.10}) and (\ref{eq-9.1.11}) together we arrive to Kirchhoff formula \begin{multline} u(x,t)= \frac{\partial\ }{\partial t} \Bigl(\frac{1}{4\pi c^2 t} \iint _{S(x,c|t|)} g(y)\,d\sigma\Bigr)+ \frac{1}{4\pi c^2 t} \iint _{S(x,c|t|)} h(y)\,d\sigma +\\ \int_0^t \frac{1}{4\pi c^2(t-\tau)} \iint_{S(x,c|t-\tau|) } f(y,\tau)\,d\sigma d\tau. \qquad \label{eq-9.1.12} \end{multline} providing solution to (\ref{eq-9.1.1})--(\ref{eq-9.1.3}).
Remark 2. As $t>0$ one can rewrite the right-hand expression in (\ref{eq-9.1.11}) as \begin{equation} \iiint_{B(x,ct)} \frac{1}{4\pi c^2 |x-y|}f(y,t-c^{-1}|x-y|)\,dy \label{eq-9.1.13} \end{equation} where we we integrate over ball $B(x,ct)$ of radius $ct$ with the center at $x$.
Definition 1. $M_r(h, x)=\frac{1}{4\pi r^2} \int_{S(x,r)} h(y)\,d\sigma$ is a spherical mean of $h$. Recall that $4\pi r^2$ is an area of $S(x,r)$.
Therefore (\ref{eq-9.1.8}) is exactly $ u(x,t)= t M_{c|t|} (h,x)$ and all other formulae (\ref{eq-9.1.9})--(\ref{eq-9.1.13}) could be modified similarly.
Consider now the same problem (\ref{eq-9.1.1})--(\ref{eq-9.1.3}) but in dimension $2$. To apply (\ref{eq-9.1.12}) we introduce a third spatial variable $x_3$ and take $f$, $g$, $h$ not depending on it; then $u$ also does not depend on $x_3$ and solves original $2D$-problem.
So, the right-hand expression in (\ref{eq-9.1.8}) becomes for $\pm t>0$ \begin{equation} \frac{1}{4\pi c^2t} \iint _{S(x,c|t|)} h(y)\,d\sigma= \pm \frac{1}{2\pi c} \iint _{B(x,c|t|)} \frac{h(y)}{\sqrt{c^2t^2-|x-y|^2}}\,dy \label{eq-9.1.14} \end{equation} where $y=(y_1,y_2)$ and we took into account that $S(x,ct)$ covers disk $B(x,c|t|)$ twice (so factor $2$ appears) and $d\sigma = \pm \frac{ct}{\sqrt{c^2t^2-|x-y|^2}}\,dy$, $dy=dy_1dy_2$.
Thus (\ref{eq-9.1.12}) implies that for $\pm t>0$ \begin{multline} u(x,t)= \pm\frac{\partial\ }{\partial t} \Bigl(\frac{1}{2\pi c} \iint _{B(x,c|t|)} \frac{g(y)}{\sqrt{c^2t^2-|x-y|^2}} \,d y\Bigr) \\ \pm \frac{1}{2\pi c } \iint _{B(x,c|t|)} \frac{h(y)}{\sqrt{c^2t^2-|x-y|^2}}\,dy \\ \pm \int_0^t \frac{1}{4\pi c } \iint_{B(x,c|t-\tau|) } \frac{f(y,\tau)}{\sqrt{c^2(t-\tau)^2-|x-y|^2}}\,dy d\tau. \qquad \label{eq-9.1.15} \end{multline}
Let $n=3$. Consider solution to inhomogeneous wave equation with a special right-hand expression \begin{equation} \Delta u -c^{-2}u_{tt} = f(x)e^{i\omega t} \label{eq-9.1.16} \end{equation} where $\omega\ne 0$ and $f(x)$ does not depend on $t$ and fast decays as $|x|\to \infty$. Assume that $g(x)=u(x,0)$ and $h(x)=u_t (x,0)$ also fast decay as $|x|\to \infty$. Plugging all these functions into Kirchhoff formula (\ref{eq-9.1.12}) and considering $|t|\gg 1$ and fixed $x$ we see that \begin{equation} |u (x,t) - v^\pm _\omega (x)e^{i\omega t}|\to 0 \qquad \text{as\ \ }t\to \pm \infty \label{eq-9.1.17} \end{equation} with \begin{equation} v^\pm_\omega= -\frac{1}{4\pi}\iiint |x-y|^{-1}e^{\mp i\omega c^{-1}|x-y|}\, dy. \label{eq-9.1.18} \end{equation} One can check easily that $v=v^\pm_\omega$ satisfies Helmholtz equation \begin{equation} \bigl(\Delta +\frac{\omega^2}{c^2}\bigr)v = f(x) \label{eq-9.1.19} \end{equation} with Sommerfeld radiating conditions \begin{align} &v = o(1) &&\text{as }r\to \infty,\label{eq-9.1.20}\\ &( \partial_r \mp ic^{-1}\omega )v= o( r^{-1}) &&\text{as }r\to \infty \label{eq-9.1.21} \end{align} where $r:=|x|$ and $\partial_r := |x|^{-1} x\cdot \nabla$.
This is called Limiting amplitude principle
Remark 3.
Remark 4. Formulae (\ref{eq-9.1.13}) and (\ref{eq-9.1.15}) could be generalized to the case of odd $n\ge 3$ and even $n\ge 2$ respectively. These formulae imply that $u(x,t)$ does not depend on $g(y),h(y)$ with $|y-x|>ct$ and on $f(y,\tau)$ with $|y-x|>c|t-\tau|$. This could be interpreted as "nothing propagates with a speed exceeding $c$". We will prove it again by completely different method in the next Section 9.2.
Remark 5.
As $n\ge 3$ is odd $u(x,t)$ is given by the following formula \begin{align} u(x,t)= c^{1-n}\kappa_n \frac{\partial\ }{\partial t} &\left(\frac{1}{t}\frac{\partial\ }{\partial t}\right)^{\frac{n-3}{2}} \Biggl(t^{-1}\iint_{S(x,c|t|)}g(y)\,d \sigma\Biggr)\notag\\ +c^{1-n}\kappa_n &\left(\frac{1}{t}\frac{\partial\ }{\partial t}\right)^{\frac{n-3}{2}} \Biggl(t^{-1}\iint_{S(x,c|t|)}h(y)\,d \sigma\Biggr) \label{eq-9.1.23} \end{align} provided $f=0$ (which could be generalized to $f\ne 0$ using Duhamel principle).
As $n\ge 2$ is even $u(x,t)$ is given by the following formula obtained by the method of descent \begin{align} \hskip{-20pt}u(x,t)\qquad&\notag\\= c^{1-n}\kappa_n \frac{\partial\ }{\partial t} &\left(\frac{1}{t}\frac{\partial\ }{\partial t}\right)^{\frac{n-2}{2}} \Biggl(\iiint_{B(x,c|t|)}\frac{g(y)}{(c^2t^2-|x-y|^2)^{\frac{1}{2}}}\,dy \Biggr)\notag\\ +c^{1-n}\kappa_n &\left(\frac{1}{t}\frac{\partial\ }{\partial t}\right)^{\frac{n-2}{2}} \Biggl(\iiint_{B(x,c|t|)}\frac{g(y)}{(c^2t^2-|x-y|^2)^{\frac{1}{2}}}\,dy \Biggr) \label{eq-9.1.24} \end{align} provided $f=0$ (which could be generalized to $f\ne 0$ using Duhamel principle).
Here $\kappa_n$ is a numerical coefficient which could be easily calculated from $g \equiv 1,\ h\equiv 0,\ f\equiv 0 \implies u\equiv 1$.
In particular, for odd $n\ge 3$ solution $u(x,t)$ does not depend on $g(y),h(y)$ with $|y-x|< ct$ and on $f(y,\tau)$ with $|y-x|< c|t-\tau|$. This could be interpreted as "no leftovers after front passed with a speed $c$". In mathematical literature this is called Huygens principle (there is another Huygens principle aka Huygens-Fresnel principle). This property is a rare commodity: adding lower-order terms to the equation breaks it.