Fokker Planck with d'Alembert operator Question:
Consider the d'Alembert operator with only one space-variable: $\Box:=\frac1{c^2}∂_t^2-∂_x^2$ and then the PDE
$$\frac{∂p}{∂t}=\frac12\Box p,\quad\forall(t,x)\in[0,∞)×\mathbb R$$
with IC $p(0, x)=\delta_{x_0}(x)$ and $p_t(0,x)=0$ (for limiting BC, see below). Does this PDE have a (semi)-closed-form solution?
Some comments:
If we discretize this using a basic implicit scheme, we see solutions that evolve like this:

So, analogous to the classic Fokker-Planck evolution $p_t = \frac 12 p_{xx},$ except that instead of one Gaussian that disperses over time around it's initial point $x_0$, we see two traveling Gaussian-like curves dispersing, where the peaks occur roughly at $\pm ct$ and move away from the initial point $x_0=0$, in this case. This, I think, intuitively captures wave-like properties and diffusion-properties. Additionally, judging by the numerical solutions, we can further impose some regularity conditions/limiting behavior: $\forall t >0$, $\lim_{x\to \pm \infty} p(t, x)=0$ and $\forall x$, $\lim_{t\to \infty} p(t, x)=0$. (The label price should be space, by the way, the time-span was $T=3/365$, $c\approx 0.307$, and $100$ time subintervals and $500$, space sub-intervals, please comment if you want me to include the R-code that implements the solver; it is not too long).
However, if we guess $p(t, x)$ is a $1/2$-mixture of two Gaussians with opposite means $\pm ct$, which we write as $u^+$ and $u^-$ with same variance $\sigma^2 t$, then $p$, it can be seen, satisfies the standard Fokker-Planck,
$$\frac{\partial p}{\partial t} = -\frac{\partial}{\partial_x} \left[\frac{c}{2}(u^+(t, x)-u^-(t,x))\right]+\frac12 \sigma^2 \frac{\partial^2 p}{\partial x^2}.$$
Now we might set the RHS side equal to $\frac 12 \Box p$ and see what this implies on $p=\frac12 u^+ +\frac 12 u^-$, but I have not been able to carry out anything fruitful with moderate effort.
I have tried separation of variables and transform techniques (Laplace followed by Fourier) and I do not quite have the fortitude (yet) to brute-force check whether this Gaussian mixture guess satisfies the PDE by computing the derivatives $p_t, p_{tt}, p_{xx}$. Naturally, this mixture guess (or something like it) seems to match characteristics of the basic heat-diffusion and the d'Alembert formula for wave-equations simultaneously (since the Gaussians depend on $(x-ct)$ and $(x+ct)$), at least informally. To be honest, judging from many numerical computations, the densities look slightly skewed, with the density moving house-left of $x_0$ left-skewed, and the density moving house-right of $x_0$, right-skewed, so perhaps then, a mixture of Gaussians is not entirely accurate, but is an okay first-guess.
Another appealing idea, but for different purposes, is to try to discretize the PDE into a difference equation that we could build a random walk $S_n$ whose PMF satisfies, and then do a scaling limit of some sort. This I am working on today. Of course this won't give us the density form but it'll at least give us a stochastic object (the RW) that is associated with it and this is useful/interesting to me.
Please comment for clarifications of the question, context, or corrections on typos/mistakes I missed, thank you.
 A: $ \newcommand{\pt}{ \partial_{t} } $$ \newcommand{\ptt}{ \partial_{tt} } $$ \newcommand{\px}{ \partial_{x} } $$ \newcommand{\pxx}{ \partial_{xx} } $$ \newcommand{\ptx}{ \partial_{tx} } $$ \newcommand{\pxt}{ \partial_{xt} } $$ \newcommand{\prt}[1]{ \left(#1\right) } $$ \newcommand{\abs}[1]{ \left|#1\right| } $$ \newcommand{\hatp}{ \hat{p} } $$ \newcommand\diff{\mathop{}\!\mathrm{d}} $I wonder if instead of
$$ \pt p = \frac{1}{2}\prt{\frac{1}{c^2}\ptt p - \pxx p}, \tag{A}$$
you mean
$$ \pt p = - \frac{1}{2}\prt{\frac{1}{c^2}\ptt p - \pxx p}, \tag{B}$$
which is more consistent with your simulation. I say this because both of these can be though of as wave equations with a "damping" term, though in the case of (A) the $\pt p$ term has the wrong sign for a damping term.
If you know your Hamiltonians, you can make physical sense of this. The wave equation has an energy which arises from physical principles (kinetic energy plus potential energy),
$$ H = \int_\mathbb{R} \frac{1}{2}\prt{\pt p}^2 + \frac{c^2}{2}\prt{\px p}^2 \diff x.$$
For the homogeneous wave equation, $\ptt p = c^2\pxx p$, you have that energy is conserved:
$$
\frac{\diff H}{\diff t} =
\int_\mathbb{R} \pt p \ptt p + c^2\px p \ptx p  \diff x = 
\int_\mathbb{R} \pt p \ptt p - c^2\pxx p \pt p  \diff x = 
\int_\mathbb{R} \pt p \prt{ \ptt p - c^2 \pxx p } \diff x = 
0,
$$
using integration by parts (and assuming $p$ and its derivatives decay at infinity). For (A), however, we would find
$$
\frac{\diff H}{\diff t} =
\int_\mathbb{R} \pt p \prt{ \ptt p - c^2 \pxx p } \diff x = 
\int_\mathbb{R} \pt p \prt{ 2c^2 \pt p } \diff x = 
2c^2 \int_\mathbb{R} \prt{ \pt p }^2 \diff x \geq 0,
$$
so the energy of the system increases in time, and you would expect solutions to blow up. The same calculation for (B) yields
$$
\frac{\diff H}{\diff t} =
-2c^2 \int_\mathbb{R} \prt{ \pt p }^2 \diff x \leq 0,
$$
so the energy of solutions decays in time, which is more consistent with the simulation you show.
Given that (B) is the right equation (though a similar treatment will work for (A)): it's a linear equation, so we try a Fourier method first. You are solving the equation on the whole real line, so we take Fourier transforms of (B), to find
$$ \pt \hatp = - \frac{1}{2}\prt{\frac{1}{c^2}\ptt \hatp + k^2 \hatp}, \tag{C}$$
where
$$
\hatp(t,k) = \int_\mathbb{R} p(t,x) e^{-ikx} \diff x.
$$
Equation (C) is now a linear ODE for each wavenumber $k$, so we make the usual ansatz $\hatp\propto e^{\lambda t}$, leading to
$$ \lambda = - \frac{1}{2}\prt{\frac{1}{c^2}\lambda^2 + k^2}. \tag{D}$$
Rewriting as
$$ \lambda^2 + 2c^2\lambda + c^2k^2 = 0,$$
we have
$$
\lambda = -c^2 \pm c\sqrt{c^2-k^2}.
$$
Whenever $\abs{k}\leq \abs{c}$, $\lambda$ is real and negative, so you have frequency modes which decay exponentially fast. We write
$$
\lambda =
-c^2 \pm c\sqrt{c^2-k^2} = 
-c^2 \pm \xi(k),
$$
and obtain the general solution
$$
\hatp(t,k) = e^{-c^2t} \prt{\hat f(k) e^{\xi(k)t} + \hat g(k) e^{-\xi(k)t}} \quad \text{for} \abs{k}< \abs{c},
$$
where $f(k)$ and $g(k)$ are constants for each value of $k$, to be determined from the datum of your problem. In the limiting case, $\abs{k}=\abs{c}$, you obtain instead
$$
\hatp(t,k) = e^{-c^2t} \prt{\hat f(k) + t \hat g(k) } \quad \text{for} \abs{k} = \abs{c},
$$
We now concentrate on $\abs{k} > \abs{c}$, which yields complex values of $\lambda$ (these correspond to the travelling pulses in your picture). For convenience, we rewrite
$$
\lambda =
-c^2 \pm c\sqrt{c^2-k^2} = 
-c^2 \mp ic\sqrt{k^2-c^2} =
-c^2 \mp i\omega(k),
$$
which yields a general solution to (C),
$$
\hatp(t,k) = e^{-c^2t} \prt{\hat f(k) e^{-i\omega(k)t} + \hat g(k) e^{+i\omega(k)t}} \quad \text{for} \abs{k}> \abs{c}.
$$
To obtain a particular solution for your datum, we take the transforms
$
\widehat{p(0,x)} = \hat \delta(x) = 1, 
$
and
$
\widehat{\pt p(0,x)} = \hat 0 = 0,
$
and use them to find $\hat f$ and $\hat g$ on each regime. In the $\abs{k}< \abs{c}$ case, for instance, setting $1=\hatp(0,k)$ and $0=\pt\hatp(0,k)$ yields $f(k)+g(k)=1$ and $(f(k)-g(k))=\frac{c^2}{\xi(k)}$, leading to
$$
f(k) = 1+\frac{c^2}{2\xi(k)}, \quad 
g(k) = 1-\frac{c^2}{2\xi(k)}.
$$
A: Thanks to the comment by @mattos and the answer and comment thread by @R_B I have found a very satisfying answer, originally proven by Kac in 1970. It is recalled and discussed further in section $5.2$ (starting at page $90$) of this excellent paper on relativistic Brownian motion. I will restate the result here for completeness, in the latter paper's notation. Also before this result is stated, the solution is explicitly written down, using Bessel functions but I am not going to copy that here. This explicit solution was first obtained by Goldstein in 1939 by considering the continuum limit of a persistent random walk model.
If $q(t,x)$ solves
$$\left(\tau_v \frac{\partial^2}{\partial t^2} + \frac{\partial}{\partial t}\right) q=D \nabla^2 q,$$
where $\tau_v$ and $D$ are a time-relaxation constant and diffusion constant, respectively, with initial conditions $q(0, x) = q_0(x)$ and zero initial velocity, then
$$q(t, x) = \frac 12 \mathbb{E} \left(q_0 \left(x+v \int_0^t (-1)^{N_s} ds\right)+q_0 \left(x-v \int_0^t (-1)^{N_s} ds\right) \right),$$
where $N_s$ is a Poisson process with intensity rate $\frac{1}{2 \tau_v}$ and $v = (D/\tau_v)^{1/2}$. The correspondence to the PDE (after you flip the signs on the RHS; see the @R_B's answer and the ensuing comment thread) in the OP is, I believe, $D = \frac{1}{2}$ and $\tau_v = \frac{1}{2c^2}$. The PDE is equation $(206)$ in the second linked paper, the explicit solution $(209a)$ both on page $91$ and Kac's stochastic solution is equation $(211)$ on page 94.
Kac remarks "[this]  is  really  extremely  amusing  and  shows  that  if  one  hits  upon  the  right  formulation  one  always  gets  more  than  one  has  bargained  for.  Our  solution,  in the  form  [above], is  extremely  reminiscent  of the  solution  of the  equation  of the  vibrating  string.   Remember   that  for  the  vibrating  string  I  had  simply  $(1/2)[q(x  +  vt)  +  q(x  —  vt)]$  and  there  was  no  average.   Now  these  two  differ  in  only  one  respect.  Time   $t$   is  replaced   by   this  "randomized   time"  $\int_0^t (-1)^{N_s}  ds$.    And    then,  because  you  don't  know  what  it  is  going  to be  exactly, you  must  average."
His short paper on this topic is a highly refreshing read. He laments that his proof is too full of computations and is inelegant (he verifies that the Laplace transform of the candidate solution solves the Laplace transform of the telegrapher PDE). "One  feels  that  one  shouldn't  have  to  compute  anything  at  all." Oh, don't we all...
