# 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?

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.

• Just so you know, this is (almost) the telegraphers equation (with $G$ or $R = 0$), except you have a $\color{red}{+} \partial_{t} p$ instead of a $\color{orange}{-} \partial_{t} p$ which is a reasonably well studied problem. Maybe a search along those lines will yield something useful (though maybe not). Jun 7, 2021 at 3:57
• @mattos A search along those lines did indeed yield something useful, or as Kac would put it "extremely amusing". See my posted answer. Jun 10, 2021 at 22:48

$$\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)}.$$

• +1 I will have to check my finite difference scheme for a sign error because I really did intend to study $(A)$. But the Hamiltonian argument is quite useful, thank you. I wrote this PDE down after studying Fokker Planck on various manifolds and eventually got curious about Minkowski space (with one space dimension). Using the sign convention, expressed in line element form, $ds^2=c^2 dt^2-dx^2$, I derived the d’Alembert operator in the OP. But I see using the other sign convention $(-,+)$ will yield the operator corresponding to $(B)$. Jun 8, 2021 at 14:46
• With that in mind, I will probably accept this answer and award the bounty. But I will wait a little so I may check some things first for my own edification. Thanks again! Jun 8, 2021 at 14:57
• No problem! I'm surprised different sign conventions yield such different equations, but I'm not very familiar with Minkowski spacetime so I could be quite misguided. In any case, seeing how (A) is indeed the relevant equation, then I'd be very curious to see where it comes from; do you have some references?
– R_B
Jun 8, 2021 at 16:32
• I do not have a reference for the PDE per se but here is more context. Following along Rogers and Williams Volume 2 textbooks, I learned probabilists may define the infinitesimal generator of a Brownian motion on a Riemannian manifold as $\frac 12 \Delta$ where $\Delta$ is the Laplace-Beltrami operator. This can be calculated (in local coordinates) directly from knowledge of the metric tensor of the space. The Fokker-Planck equation is then defined as $p_t = \frac 12 \Delta$. Naturally, after working out this for a few examples like the sphere, disk, etc, (continued) Jun 8, 2021 at 21:59
• Also some people seem to be doing stochastic processes in relativity, arxiv.org/abs/0812.1996 but this is really far away from my area, so I cannot comment. So, in short, it seems your question was a good one, because nobody seems to know what the "right" answer is :)
– R_B
Jun 10, 2021 at 16:32

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...