$\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(\boldsymbol{x},t)= \frac{1}{4\pi c ^2 t} \iint _{S(\boldsymbol{x},ct)} h(\boldsymbol{y})\,d\sigma \label{eq-9.1.4} \end{equation} where we integrate along sphere $S(\boldsymbol{x},ct)$ with a center at $\boldsymbol{x}$ and radius $ct$; $d\sigma$ is an area element.
Let us prove (\ref{eq-9.1.4}) first as $h(\boldsymbol{x})=e^{i\boldsymbol{x}\cdot \boldsymbol{\xi}}$ with $\boldsymbol{\xi}\in \mathbb{R}^3\setminus 0$; we use the standard notation $\boldsymbol{x}\cdot \boldsymbol{\xi}=x_1\xi_1+x_2\xi_2+x_3\xi_3$. In this case \begin{equation} u(x,t)=e^{i\boldsymbol{x}\cdot \boldsymbol{\xi}}c^{-1}|\boldsymbol{\xi}|^{-1} \sin (ct|\boldsymbol{\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(\boldsymbol{x},ct)} e^{i\boldsymbol{y}\cdot \boldsymbol{\xi}}\,d\sigma= \frac{1}{4\pi c ^2 t}e^{i\boldsymbol{x}\cdot \xi}\iint _{S(0,ct)} e^{i\boldsymbol{z}\cdot \boldsymbol{\xi}}\,d\sigma \end{equation*} where we changed variable $\boldsymbol{y}=\boldsymbol{x}+\boldsymbol{z}$ with $\boldsymbol{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 $\boldsymbol{\xi}=(0,0,\omega)$ with $\omega =|\boldsymbol{\xi}|$ and introduce corresponding spherical coordinates $(\rho,\phi,\theta)$; then on $S(0,ct)$ $\rho=ct$, $\boldsymbol{z}\cdot \boldsymbol{\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^{i\boldsymbol{x}\cdot\boldsymbol{\xi}}$ and dividing by $4\pi c^2t$ we get $e^{i\boldsymbol{x}\cdot\boldsymbol{\xi}} c^{-1}|\boldsymbol{\xi}|^{-1} \sin (ct|\boldsymbol{\xi}|)$ which is the right-hand expression of (\ref{eq-9.1.5}).
So, for $h(\boldsymbol{x})=e^{i\boldsymbol{x}\cdot\boldsymbol{\xi}}$ (\ref{eq-9.1.4}) has been proven. However the general function $h(\boldsymbol{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{align} h(\boldsymbol{x})=&\iiint \hat{h}(\xi) e^{i\boldsymbol{x}\cdot\boldsymbol{\xi}}\,d\boldsymbol{\xi},\label{eq-9.1.6}\\[3pt] \hat{h} (\boldsymbol{\xi})=&(2\pi)^{-n}\iiint h(\boldsymbol{x}) e^{-i\boldsymbol{x}\cdot\boldsymbol{\xi}}\,d\boldsymbol{x} \label{eq-9.1.7} \end{align} 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(\boldsymbol{x},t)= \frac{1}{4\pi c ^2t} \iint _{S(\boldsymbol{x},c|t|)} h(\boldsymbol{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(\boldsymbol{x},t)$ be given by (\ref{eq-9.1.8}) with $h$ replaced by $g$; 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(\boldsymbol{x},t)= \frac{\partial\ }{\partial t} \Bigl(\frac{1}{4\pi c^2 t} \iint _{S(\boldsymbol{x},c|t|)} g(\boldsymbol{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(\boldsymbol{x},t)= \frac{\partial\ }{\partial t}\Bigl(\frac{1}{4\pi c^2 t} \iint_{S(\boldsymbol{x},c|t|)} g(\boldsymbol{y})\,d\sigma\Bigr)+ \frac{1}{4\pi c ^2 t} \iint _{S(\boldsymbol{x},c|t|)} h(\boldsymbol{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). Namely, consider problem \begin{align*} & U_{tt}-c^2\Delta U=0,\\[3pt] & U|_{t=\tau}=0,\\[3pt] & U_t|_{t=\tau}=f(\boldsymbol{x},\tau). \end{align*} Its solution is given by \begin{equation*} U(\boldsymbol{x},t,\tau) = \frac{1}{4\pi c^2(t-\tau)} \iint_{S(\boldsymbol{x},c|t-\tau|) } f(\boldsymbol{y},\tau)\,d\sigma \end{equation*} and therefore \begin{equation} u(\boldsymbol{x},t) = \int_0^t\frac{1}{4\pi c^2(t-\tau)} \iint_{S(\boldsymbol{x},c|t-\tau|)} f(\boldsymbol{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(\boldsymbol{x},t)= \frac{\partial\ }{\partial t} \Bigl(\frac{1}{4\pi c^2 t} \iint _{S(\boldsymbol{x},c|t|)} g(\boldsymbol{y})\,d\sigma\Bigr)+ \frac{1}{4\pi c^2 t} \iint _{S(\boldsymbol{x},c|t|)} h(\boldsymbol{y})\,d\sigma +\\ \int_0^t \frac{1}{4\pi c^2(t-\tau)} \iint_{S(\boldsymbol{x},c|t-\tau|) } f(\boldsymbol{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(\boldsymbol{x},ct)} \frac{1}{4\pi c^2 |\boldsymbol{x}-\boldsymbol{y}|} f(\boldsymbol{y},t-c^{-1}|\boldsymbol{x}-\boldsymbol{y}|)\,d\boldsymbol{y} \label{eq-9.1.13} \end{equation} where we we integrate over ball $B(\boldsymbol{x},ct)$ of radius $ct$ with the center at $\boldsymbol{x}$. It is called time-delayed potential.
Definition 1. \begin{equation} M_r(h, \boldsymbol{x})=\frac{1}{4\pi r^2} \iint_{S(\boldsymbol{x},r)} h(y)\,d\sigma= \frac{1}{4\pi}\iint_{S(0,1)} h(\boldsymbol{x}+\boldsymbol{\eta} r)\,d\eta, \label{eq-9.1.14} \end{equation} where $d\sigma$ and $d\eta$ are area elements respectively on $S(\boldsymbol{x},r)$ and $S(0,1)$, is a spherical mean of $h$. Recall that $4\pi r^2$ is an area of $S(\boldsymbol{x},r)$.
Therefore (\ref{eq-9.1.8}) is exactly $u(\boldsymbol{x},t)= t M_{c|t|} (h,\boldsymbol{x})$ and all other formulae (\ref{eq-9.1.9})--(\ref{eq-9.1.13}) could be modified similarly.
Remark 3. Another proof Kirchhoff formula is based on spherical means. First, one can prove that if $u(\boldsymbol{x},t)$ satisfies wave equation $u_{tt}-c^2\Delta u=0$, then $v(\boldsymbol{x};r,t)= r M_r(u,\boldsymbol{x})$ satisfies 1D-wave equation $v_{tt}-c^2v_{rr}=0$. Here $\boldsymbol{x}$ is considered as a parameter.
Then $v$ can be written through D'Alembert formula. Plugging it into $\partial _r v |_{r=0}= u(\boldsymbol{x};t)$, we arrive to (\ref{eq-9.1.4}) if $u|_{t=0}=0$, $u_t|_{t=0}=h(\boldsymbol{x})$.
Remark 4. Yet another proof could be done, using Radon transform, in particular, it's intertwining property.
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(\boldsymbol{x},c|t|)} h(\boldsymbol{y})\,d\sigma= \pm \frac{1}{2\pi c} \iint _{D(\boldsymbol{x},c|t|)} \frac{h(\boldsymbol{y})}{\sqrt{c^2t^2-|\boldsymbol{x}-\boldsymbol{y}|^2}}\,dy \label{eq-9.1.15} \end{equation} where $\boldsymbol{y}=(y_1,y_2)$ and we took into account that $S(\boldsymbol{x},ct)$ covers disk $D(\boldsymbol{x},c|t|)$ twice (so factor $2$ appears) and \begin{gather*} d\sigma = \pm \frac{ct}{\sqrt{c^2t^2-|\boldsymbol{x}-\boldsymbol{y}|^2}}\,dy,\qquad dy=dy_1dy_2. \end{gather*}
Thus 3D-formula (\ref{eq-9.1.12}) implies that for $\pm t>0$ the following 2D-formula holds:
\begin{align}
u(\boldsymbol{x},t)= &\pm\frac{\partial\ }{\partial t}
\Bigl(\frac{1}{2\pi c} \iint _{D(\boldsymbol{x},c|t|)} \frac{g(\boldsymbol{y})}{\sqrt{c^2t^2-|\boldsymbol{x}-\boldsymbol{y}|^2}}\,d y\Bigr) \notag\\
&\pm \frac{1}{2\pi c } \iint _{D(\boldsymbol{x},c|t|)} \frac{h(\boldsymbol{y})}{\sqrt{c^2t^2-|\boldsymbol{x}-\boldsymbol{y}|^2}}\,dy \notag\\
&\pm \int_0^t \frac{1}{2\pi c } \iint_{D(\boldsymbol{x},c|t-\tau|) }
\frac{f(\boldsymbol{y},\tau)}{\sqrt{c^2(t-\tau)^2-|\boldsymbol{x}-\boldsymbol{y}|^2}}\,dy d\tau.
\label{eq-9.1.16}
\end{align}
Let $n=3$. Consider solution to inhomogeneous wave equation with a special right-hand expression \begin{equation} \Delta u -c^{-2}u_{tt} = f(\boldsymbol{x})e^{i\omega t} \label{eq-9.1.17} \end{equation} where $\omega\ne 0$ and $f(x)$ does not depend on $t$ and fast decays as $|x|\to \infty$. Assume that $g(\boldsymbol{x})=u(\boldsymbol{x},0)$ and $h(\boldsymbol{x})=u_t (\boldsymbol{x},0)$ also fast decay as $|\boldsymbol{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 the first two terms tend to $0$ for fixed $\boldsymbol{x}$ and $\pm t\to \infty$ while the last term becomes \begin{multline*} \int_0^t \frac{1}{4\pi c^2(t-\tau)} \iint_{S(\boldsymbol{x},c|t-\tau|) } f(\boldsymbol{y})e^{i\omega \tau}\,d\sigma d\tau\\ = -\frac{1}{4\pi}\iiint_{B(\boldsymbol{x},c|t|)} |\boldsymbol{x}-\boldsymbol{y}|^{-1}e^{\mp i\omega c^{-1}|\boldsymbol{x}-\boldsymbol{y}|}\, dy \end{multline*} and therefore \begin{equation} |u (\boldsymbol{x},t) - v^\pm _\omega (\boldsymbol{x})e^{i\omega t}|\to 0 \qquad \text{as }\ t\to \pm \infty \label{eq-9.1.18} \end{equation} with \begin{equation} v^\pm_\omega(\boldsymbol{x})= -\frac{1}{4\pi} \iiint |\boldsymbol{x}-\boldsymbol{y}|^{-1}e^{\mp i\omega c^{-1}|x-y|}\, dy. \label{eq-9.1.19} \end{equation} One can check easily that $v=v^\pm_\omega(\boldsymbol{x})$ satisfies Helmholtz equation \begin{equation} \bigl(\Delta +\frac{\omega^2}{c^2}\bigr)v = f(\boldsymbol{x}) \label{eq-9.1.20} \end{equation} with Sommerfeld radiating conditions \begin{align} &v = o(1) &&\text{as }r\to \infty,\label{eq-9.1.21}\\ &( \partial_r \mp ic^{-1}\omega )v= o( r^{-1}) &&\text{as }r\to \infty \label{eq-9.1.22} \end{align} where $r:=|\boldsymbol{x}|$ and $\partial_r := |\boldsymbol{x}|^{-1} \boldsymbol{x}\cdot \nabla$.
This is called Limiting amplitude principle.
Remark 5.
Let $n=3$. Consider spherically symmetric solution $u(\boldsymbol{x},t)=u(r,t)$ of the wave equation. Since for spherically symmetric $u$ $\displaystyle{\Delta u = u_{rr}+\frac{2}{r}u_r}$ we can rewrite it as \begin{gather*} u_{tt}-c^2\bigl(u_{rr}+\frac{2}{r}u_r\bigr)=0 \end{gather*} or, equivalently} \begin{gather} (ru)_{tt}-c^2(ru)_{rr}=0. \label{eq-9.1.24} \end{gather} Therefore $\displaystyle{ru(r,t)=\phi (r-ct)+\psi(r+ct)}$ and since $(ru)|_{r=0}=0$ we conclude that $\psi(ct)=-\phi(-ct)$ and therefore \begin{gather} u=\frac{1}{r}\bigl(\phi(r-ct)-\phi (-r-ct)\bigr). \label{eq-9.1.25} \end{gather} This is spherical wave. \pause Similar non-trivial solutions $u=r^\alpha \bigl(\phi(r-ct)-\phi (-r-ct)\bigr)$ do not exist for $n=2$.
Remark 6.