1.3: Fokker-Planck Equations

Fokker-Planck Equations and Diffusion

We have already generalized the equations governing Markov processes to account for systems that evolve continuously in time, which resulted in the master equations. In this section, we adapt these equations further so that they may be suitable for the description of systems with a continuum of states, rather than a discrete, countable number of states.

Motivation and Derivation

Consider once again the infinite one-dimensional lattice, with lattice spacing \(\Delta x\) and timestep size \(\Delta t\). In the previous section, we wrote down the master equation (discrete sites, continuous time) for this system, but here we will begin with the Markov chain expression (discrete sites, discrete time) for the system,

\[P(n, s+1)=\frac{1}{2}(P(n+1, s)+P(n-1, s))\]

In terms of \(\Delta x\) and \(\Delta t\), this equation is

\[P(x, t+\Delta t)=\frac{1}{2}[P(x+\Delta x, t)+P(x-\Delta x, t)]\]

Rearranging the previous equation as a finite difference, as in

\[\frac{P(x, t+\Delta t)-P(x, t)}{\Delta t}=\frac{(\Delta x)^{2}}{2 \Delta t} \cdot \frac{P(x+\Delta x, t)+P(x-\Delta x, t)-2 P(x, t)}{(\Delta x)^{2}}\]

and taking the limits \(\Delta x \rightarrow 0, \Delta t \rightarrow 0\), we arrive at the following differential equation:

\[\frac{\partial}{\partial t} P(x, t)=D \frac{\partial^{2}}{\partial x^{2}} P(x, t)\]

where \(D=\frac{(\Delta x)^{2}}{2 \Delta t}\). This differential equation is called a diffusion equation with diffusion constant \(D\), and it is a special case of the Fokker-Planck equation, which we will introduce shortly. The most straightforward route to the solution of the diffusion equation is via spatial Fourier transformation,

\[\tilde{P}(k, t)=\int_{-\infty}^{\infty} P(x, t) e^{i k x} d x\]

In Fourier space, the diffusion equation reads

\[\frac{\partial}{\partial t} \tilde{P}(k, t)=-D k^{2} \tilde{P}(k, t)\]

and its solution is

\[\tilde{P}(k, t)=\tilde{P}(k, 0) e^{-D k^{2} t}\]

If we take a delta function \(P(x, 0)=\delta\left(x-x_{0}\right)\) centered at \(x_{0}\) as the initial condition, the solution in \(x\)-space is

\[P(x, t)=\frac{1}{\sqrt{4 \pi D t}} e^{-\frac{\left(x-x_{0}\right)^{2}}{4 D t}}\]

Thus the probability distribution is a Gaussian in \(x\) that spreads with time. Notice that this solution is essentially the same as the long-time solution to the spatially discretized version of the problem presented in the previous example.

We are now in a position to consider a generalization of the diffusion equation known as the Fokker-Planck equation. In addition to the diffusion term \(D \frac{\partial^{2}}{\partial x^{2}}\), we introduce a term linear in the first derivative with respect to \(x\), which accounts for drift of the center of the Gaussian distribution over time.

Consider a diffusion process on a three-dimensional potential energy surface \(U(\mathbf{r})\). Conservation of probability requires that

\[\dot{P}(\mathbf{r}, t)=-\nabla \cdot \mathbf{J}\]

where \(\mathbf{J}\) is the probability current, \(\mathbf{J}=-D \nabla P+\mathbf{J}_{U}\), and \(\mathbf{J}_{U}\) is the current due to the potential \(U(\mathbf{r})\). At equilibrium, we know that the probability current \(\mathbf{J}=0\) and that the probability distribution should be Boltzmann-weighted according to energy, \(P_{\mathrm{eq}}(\mathbf{r}) \propto e^{-\beta U(\mathbf{r})}\). Therefore, at equilibrium,

\[-D \beta \nabla U(\mathbf{r}) P(\mathbf{r})_{\mathrm{eq}}+\mathbf{J}_{U}=0\]

Solving Eq.(1.36) for \(\mathbf{J}_{U}\) and plugging the result into Eq.(1.35) yields the Fokker-Planck equation,

\[\dot{P}(\mathbf{r}, t)=D \nabla[\nabla P(\mathbf{r}, t)+\beta \nabla U(\mathbf{r}) P(\mathbf{r}, t)]\]

Properties of Fokker-Planck Equations

Let’s return to one dimension to discuss some salient features of the Fokker-Planck equation.

• First, the Fokker-Planck equation gives the expected results in the long-time limit:

\[\lim _{t \rightarrow \infty} P=P_{\text {eq }} \text { with } \dot{P}=0\]

• Also, if we define the average position \(\bar{x}=\int_{-\infty}^{\infty} x P(x) d x\), then the differential form of the Fokker-Planck equation can be used to verify that

\[\dot{\bar{x}}=D \beta\left(-\frac{\partial}{\partial x} \overline{U(x)}\right)\]

Since the quantity in parentheses is just the average force \(\bar{F}\), Eq. \((1.39)\) can be combined with the Einstein relation \(D \beta \zeta=1\) (see section 1.4) to justify that \(\zeta \bar{v}=\bar{F}\); the meaning and significance of this equation, including the definition of \(\zeta\), will be discussed in section 1.4.

• The Fokker-Planck equation is linear in the first and second derivatives of \(P\) with respect to \(x\); it turns out that any spatial operator that is a linear combination of \(\frac{\partial}{\partial x}, x \frac{\partial}{\partial x}\), and \(\frac{\partial^{2}}{\partial x^{2}}\) will define a Gaussian process when used to describe the time evolution of a probability density. Thus, both the diffusion equation and the more general Fokker-Planck equation will generally always describe a Gaussian process. - One final observation about the Fokker-Planck equation is that it is only analytically solvable in a small number of special cases. This situation is exacerbated by the fact that it is not of Hermitian (self-adjoint) form. However, we can introduce the change of variable \(P=e^{-\frac{\beta U}{2}} \Phi\); in terms of \(\Phi\), the Fokker-Planck equation is Hermitian,

\[\frac{\partial \Phi}{\partial t}=D\left[\nabla^{2} \Phi-U_{\mathrm{eff}} \Phi\right]\]

where \(U_{\text {eff }}=\frac{(\beta \nabla U)^{2}}{4}-\frac{\beta \nabla^{2} U}{2}\). This transformed Fokker-Planck equation now bears the same functional form as the time-dependent Schrödinger equation, so all of the techniques associated with its solution can likewise be applied to Eq.(1.40).

Example: One of the simplest, yet most useful, applications of the Fokker-Planck equation is the description of the diffusive harmonic oscillator, which can be treated analytically. Here we solve the Fokker-Planck equation for the one-dimensional diffusive oscillator with frequency \(\omega\). The differential equation is

\[\frac{\partial P}{\partial t}=D \frac{\partial^{2}}{\partial x^{2}} P+\gamma \frac{\partial}{\partial x}(x P)\]

where \(\gamma=m \omega^{2} D \beta\). We can solve this equation in two steps: first, solve for the average position using Eq. \((1.39)\),

\[\dot{\bar{x}}=-\gamma \bar{x}\]

Given the usual delta function initial condition \(P(x, 0)=\delta\left(x-x_{0}\right)\), the average position is given by

\[\bar{x}(t)=x_{0} e^{-\gamma t}\]

Thus, memory of the initial conditions decays exponentially for the diffusive oscillator.

Then, since the Fokker-Planck equation is linear in \(P\) and bilinear in \(x\) and \(\frac{\partial}{\partial x}\), the full solution must take the form of a Gaussian, so we can write

\[P\left(x_{0}, x, t\right)=\frac{1}{\sqrt{2 \pi \alpha(t)}} \exp \left[-\frac{(x-\bar{x}(t))^{2}}{2 \alpha(t)}\right]\]

where \(\bar{x}(t)\) is the time-dependent mean position and \(\alpha(t)\) is the time-dependent standard deviation of the distribution. But we’ve already found \(\bar{x}(t)\), so we can substitute it into the solution,

\[P\left(x_{0}, x, t\right)=\frac{1}{\sqrt{2 \pi \alpha(t)}} \exp \left[-\frac{\left(x-x_{0} e^{-\gamma t}\right)^{2}}{2 \alpha(t)}\right]\]

Finally, from knowledge that the equilibrium distribution must satisfy the stationary condition

\[P_{\mathrm{eq}}(x)=\int_{-\infty}^{\infty} P\left(x_{0}, x, t\right) P_{\mathrm{eq}}\left(x_{0}\right) d x_{0}\]

we can determine that

\[\alpha(t)=\frac{1-e^{-2 \gamma t}}{m \omega^{2} \beta}\]

Thus the motion of the diffusive oscillator is fully described.

The long and short-time limits of \(P\left(x_{0}, x, t\right)\) are both of interest to us. At short times,

\[\lim _{t \rightarrow 0} P\left(x_{0}, x, t\right)=\sqrt{\frac{1}{4 \pi D t}} \exp \left[-\frac{\left(x-x_{0}\right)^{2}}{4 D t}\right]\]

and the evolution of the probability looks like that of a random walk. In the long-time limit, on the other hand, we find the equilibrium probability distribution

\[\lim _{t \rightarrow \infty} P\left(x_{0}, x, t\right)=\sqrt{\frac{m \omega^{2} \beta}{2 \pi}} \exp \left[-\frac{1}{2} m \omega^{2} \beta x^{2}\right]\]

which is Gaussian with no mean displacement and with variance determined by a thermal parameter and a parameter describing the shape of the potential. A Gaussian, Markovian process that exhibits exponential memory decay, such as this diffusive oscillator, is called an Ornstein-Uhlenbeck process.