On the Langevin and Fokker-Planck Equations - 6/26

Introduction

Statistical Mechanics is cool. It’s also shockingly useful. It’s one of those areas of physics that ends up being helpful for a bunch of other areas and for lots of flavours of applied maths.

Here we discuss the Langevin (do not ask me to pronounce that word) and then the Fokker-Planck Equations, obtaining the fun result that a system immersed in a fluid that is at equilibrium also tends towards equilibrium.

Movement of a Fluid-Immersed Colloid

Consider a colloid of mass \(M\) with position \(x\) and momentum \(P\) immersed in a fluid composed of particles with positions \(q\) and momenta $p$$. For simplicity, we will treat this problem in one dimension, but our results easily generalise to higher dimensions. In this context, a colloid is essentially just a large particle relative to the particles composing the fluid, where “large” means a length scale of at least three orders of magnitude larger than the fluid particles.

Our Hamiltonian has four sectors. First, the standard description of the colloid motion, an interaction term between the colloid and the fluid, a fluid particle interaction term, and an interaction term between fluid particles, \begin{equation} H=\frac{P^2}{2M}+V_{ext}(x)+\sum_iV_{FC}\left(x-q_i\right)+\sum_i\left[\frac{p_i^2}{2m}+V_{ext}\left(q_i\right)\right]+\sum_{i\neq j}V_{FF}\left(q_i-q_j\right). \end{equation}

This is an extremely general Hamiltonian that is rather unapproachable. We simplify by first assuming that the fluid is equilibrated, such that we may ignore the impact of the external potential on the fluid particles and remove the fluid-fluid interaction terms for the same reason. Here, we have assumed a separation of scales, in that the time it takes the fluid to equilibrate is much shorter than the time describing any of the dynamics of the colloid. We will discuss this equilibrium condition more later.

That leaves the colloid-fluid interactions. In practice, this interaction is rather complicated, but we will assume it can be approximated as a harmonic oscillator. Effectively, we can imagine all the fluid particles attached to the colloid with little springs, interacting with the colloid and oscillating over time. Setting \(m=1\) for convenience and assuming \(M\gg1\), the Hamiltonian becomes, \begin{equation} H=\frac{P^2}{2M}+V_{ext}(x)+\sum_i\left[\frac{p_i^2}{2}+\frac{\omega_i^2}{2}\left(q_i-x\right)^2\right]. \end{equation} We can extract equations of motion from this Hamiltonian for the colloid, \begin{equation} M\dot{x}=P,\quad\dot{P}=-\partial_xV(x)+\sum_i\omega_i^2\left(q_i-x\right), \end{equation} and for the particles, \begin{equation} \dot{q}_i=p_i,\quad\dot{p}_i=-\omega_i^2\left(q_i-x\right). \end{equation} We will solve the motion of the particles and insert this into the motion of the colloid.

The colloid motion can be expressed as a harmonic oscillator plus a constant driving term, \begin{equation}\label{eq: q harmonic oscillator equation} \ddot{q}_i(t)=-\omega_i^2q_i(t)+\omega_i^2 x(t). \end{equation} The homogenous solution is straightforward, with a cosine term and a sine term, \begin{equation} q_i^H(t)=A_i\cos\left(\omega_i t\right)+B_i\sin\left(\omega_i t\right). \end{equation} For the particular solution, we use the Green’s function solution, \begin{equation} q_i^P(t)=\int_0^tds\ f_i(t-s)x(s). \end{equation} Inserting this solution into Eq. \ref{eq: q harmonic oscillator equation}, we find, \begin{equation} \ddot{q}_i^P(t)+\omega_i^2q_i^P(t)=f_i(0)\dot{x}(t)+\partial_xf_i(0)x(t)+\int_0^tds\ \left[ \partial_x^2f_i(t-s)-\omega_i^2f_i(t-s) \right]x(s)=\omega_i^2x(t). \end{equation} Solving for the Green’s function, we see, \begin{equation} f_i(t)=\omega_i\sin\left(\omega_i t\right). \end{equation} And combining the homogeneous and particular solutions to solve for \(q_i(t)\), \begin{equation} q_i(t)=q_i(0)\cos\left(\omega_i t\right)+\frac{p_i(0)}{\omega_i}\sin\left(\omega_i t\right)+\int_0^tds\ \omega_i\sin\left(\omega_i\left(t-s\right)\right)x(s). \end{equation} Via the integral product rule, we note, \begin{equation} \int_0^tds\ \omega_i\sin\left(\omega_i\left(t-s\right)\right)x(s)=x(t)-\cos\left(\omega_i\left(t\right)\right)x(0)-\int_0^tds\ \cos\left(\omega_i\left(t-s\right)\right)\dot{x}(s). \end{equation}

Now computing the colloid motion, first we find, \begin{equation} x(t)-q_i(t)=\int_0^tds\ \cos\left(\omega_i\left(t-s\right)\right)\dot{x}(s)+\left(x(0)-q_i(0)\right)\cos\left(\omega_i t\right)-\frac{p_i(0)}{\omega_i}\sin\left(\omega_i t\right). \end{equation} We return to our momentum equation of motion for the colloid, \begin{equation} \dot{P}(t)=-\partial_xV(x)-\int_0^tds\ \left[\sum_i\omega_i^2\cos\left(\omega_i\left(t-s\right)\right)\dot{x}(s)\right] \end{equation} \begin{equation} +\sum_i\left[\omega_ip_i(0)\sin\left(\omega_i t\right)+\omega_i^2\left(q_i(0)-x(0)\right)\cos\left(\omega_i t\right)\right].\nonumber \end{equation} We rewrite in terms of a memory kernel, \(K(t-s)\) and fluctations depending on the initial condition \(\xi(t)\), \begin{equation}\label{eq: kernel form langevin} \dot{P}(t)=-\partial_xV(x)-\int_0^tds\ K(t-s)\dot{x}(s)+\xi(t), \end{equation} with, \begin{equation}\label{eq: memory kernel definition} K(t-s)=\sum_i\omega_i^2\cos\left(\omega_i\left(t-s\right)\right),\quad\xi(t)=\sum_i\left[\omega_ip_i(0)\sin\left(\omega_i t\right)+\omega_i^2\left(q_i(0)-x(0)\right)\cos\left(\omega_i t\right)\right]. \end{equation}

This should begin to look at least a little familiar. Now we take a slight deviation to discuss our assumption of the fluid being at equilibrium.

Simplifying the Fluctuations

We assumed that the fluid is at equilibrium, so at time \(t=0\) the probability distribution of the fluid position and momenta follows the Boltzmann distribution, \begin{equation} \mathcal{P}\left({q_i(0),p_i(0)}\right)\sim\exp\left[-\beta H_{F}\right]=\frac{1}{Z}\exp\left[-\beta\sum_i\left(\frac{p_i(0)^2}{2}+\frac{\omega_i}{2}\left(q_i(0)-x(0)\right)^2\right)\right], \end{equation} where \(\beta\) is the Boltzmann weight. We use \(\sim\) since we exclude the normalization factors, which are not significant in this calculation. Because the exponential is separable, the global probability distribution is a product of the probability distributions for each position and momentum, \begin{equation} \mathcal{P}\left({q_i(0),p_i(0)}\right)\sim\prod_i\mathcal{P}\left(q_i(0)\right)\mathcal{P}\left(p_i(0)\right), \end{equation} where \begin{equation} \mathcal{P}\left(q_i(0)\right)\sim\exp\left[-\beta\frac{\omega_i^2}{2}\left(q_i(0)-x(0)\right)^2\right],\quad \mathcal{P}\left(p_i(0)\right)\sim\exp\left[-\beta\frac{p_i(0)^2}{2}\right]. \end{equation}

Now the cool bit. Return to our fluctuation term, \(\xi(t)\). From our probability distributions, we see that \(\xi(t)\) is the sum of Gaussian Random Variables (GRVs), so \(\xi(t)\) must also be a GRV. Since \(\xi(t)\) is linear in \(p_i(0)\) and \(q_i(0)-x(0)\), and the mean of those distributions is zero, the mean of \(\xi(t)\) is identically zero, \begin{equation} \langle\xi(t)\rangle=0. \end{equation} The covariance is more interesting, consider \begin{equation}\label{eq: xi covariance} \langle\xi(t)\xi(t’)\rangle=\left\langle\left(\sum_i\left[\omega_ip_i(0)\sin\left(\omega_i t\right)+\omega_i^2\left(q_i(0)-x(0)\right)\cos\left(\omega_i t\right)\right]\right)\left(\sum_j\left[t\rightarrow t’\right]\right)\right\rangle, \end{equation} where \(\left[t\rightarrow t'\right]\) schematically refers the same bracketed term with \(t\) swapped to \(t'\). This looks spooky but is really not too bad. Firstly, because the momenta for each particle are independent, \begin{equation} \langle p_i(0)p_j(0)\rangle=\langle p_i(0)\rangle\langle p_j(0)\rangle\quad\text{if }i\neq j, \end{equation} \begin{equation} \langle p_i(0)p_i(0)\rangle=\frac{1}{\beta}. \end{equation} Similarly, for the positions, \begin{equation} \langle \left(q_i(0)-x(0)\right)\left(q_j(0)-x(0)\right)\rangle=\langle q_i(0)-x(0)\rangle\langle q_j(0)-x(0)\rangle\quad\text{if }i\neq j, \end{equation} \begin{equation} \langle \left(q_i(0)-x(0)\right)\left(q_i(0)-x(0)\right)\rangle=\frac{1}{\omega_i^2\beta}. \end{equation} So, Eq. \ref{eq: xi covariance} simplifies to, \begin{equation} \langle\xi(t)\xi(t’)\rangle=\sum_i\frac{\omega_i^2}{\beta}\left[\sin\left(\omega_i t\right)\sin\left(\omega_i t’\right)+\cos\left(\omega_i t\right)\cos\left(\omega_i t’\right)\right]. \end{equation} This reduces nicely to our memory kernel, Eq. \ref{eq: memory kernel definition}, \begin{equation} \langle\xi(t)\xi(t’)\rangle=\sum_i\frac{\omega_i^2}{\beta}\cos\left(\omega_i\left(t-t’\right)\right)=\frac{1}{\beta}K(t-t’). \end{equation} We have seen that we can express the fluctuations \(\xi(t)\) as a Gaussian Random Variable with a covariance given by the memory kernel \(K(t-t')\).

The Langevin Equation

The final key is in recalling our assumption that the bath relaxation time is much shorter than the timescale associated with the colloid. This allows us to simplify the memory kernel into a delta function, \(K(t-s)\sim\delta(t-s)\). Physically, we are saying that because the fluid relaxes quickly to equilibrium relative to the colloid, the “memory” contributions of the fluid are negligible, and we can take its impacts to be instantaneous. Letting, \begin{equation} K(t-t’)=2\gamma\delta\left(t-t’\right), \end{equation} the memory term becomes, \begin{equation} \int_0^tds\ K(t-s)\dot{x}(s)=\gamma\dot{x}(t). \end{equation} After consulting Eq. \ref{eq: kernel form langevin}, we arrive at the Langevin equation, \begin{equation}\label{eq: langevin equation} M\ddot{x}=-\partial_xV(x)-\gamma\dot{x}+\sqrt{\frac{2\gamma}{\beta}}\eta(t), \end{equation} where \(\eta(t)\) is a Gaussian random variable with zero mean and delta covariance, \begin{equation} \langle\eta(t)\rangle=0,\quad\langle\eta(t)\eta(t’)\rangle=\delta(t-t’). \end{equation} This is neat! We have shown that the motion of a colloid immersed in a fluid follows the usual potential gradient plus drag from moving through the fluid plus random thermal fluctuations. However, we will soon see that Eq. \ref{eq: langevin equation} is actually disastrous for doing classical calculus, and we will need a whole new machinery to operate on variables that evolve stochastically.

Itô Calculus

Let’s begin by considering the motion of a variable driven only by random fluctuations, \begin{equation}\label{eq: ito example variable} \dot{x}(t)=\eta(t), \end{equation} where \(\eta(t)\) is a GRV as described earlier. We can solve Eq. \ref{eq: ito example variable} formally to find, \begin{equation} x(t)=x(0)+\int_0^tds\ \eta(s). \end{equation} Then, taking the time derivative of \(x(t)\), \begin{equation} \dot{x}(t)=\lim_{\delta t\rightarrow0}\left[\frac{x\left(t+\delta t\right)-x(t)}{\delta t}\right]=\lim_{\delta t\rightarrow0}\left[\frac{1}{\delta t}\int_{t}^{t+\delta t}ds\ \eta(s)\right]. \end{equation} This is not ideal; consider the covaraiance of \(\frac{1}{\delta t}\int_{t}^{t+\delta t}ds\ \eta(s)\), \begin{equation}\label{eq: eta covar problem} \left\langle\frac{1}{\delta t^2}\int_{t}^{t+\delta t}\int_{t}^{t+\delta t}duds\ \eta(s)\eta(u)\right\rangle=\frac{1}{\delta t^2}\int_{t}^{t+\delta t}\int_{t}^{t+\delta t}duds\ \delta\left(s-u\right)=\frac{1}{\delta t^2}\int_{t}^{t+\delta t}du=\frac{1}{\delta t}. \end{equation} When we take the limit of \(\delta t\rightarrow0\), we see the variance of \(\dot{x}\) is infinite! Formally speaking, \(x(t)\) is not differentiable, which is the first problem.

For the second problem, consider, \begin{equation}\label{eq: chain rule result naive} \left\langle x^2(t)\right\rangle=\int_0^t\int_0^tduds\ \langle\eta(s)\eta(u)\rangle=t. \end{equation} Now via the chain rule, \begin{equation} \frac{d}{dt}\left[x^2(t)\right]=2x(t)\dot{x}(t)\ni2\int_0^tds\ \eta(s)\eta(t)\rightarrow\frac{d}{dt}\left[\langle x^2(t)\rangle\right]\sim\langle x(t)\eta(t)\rangle. \end{equation} Let’s think about this. We know \(x(t)\) is the sum of GRV “pushes” since \(\dot{x}(t)=\eta(t)\). And because \(\eta(t)\) is the current Gaussian “push”, it should be independent of \(x(t)\), yielding \(\langle x(t)\eta(t)\rangle=0\). But that’s in contradiction with our chain rule result, Eq. \ref{eq: chain rule result naive}. The second problem with naively applying the rules of calculus to a stochastic variable is that the chain rule is incompatible with causality!

There are two paths to proceed. Either we modify causality and keep the chain rule, or keep causality and modify the chain rule. Both approaches are valid, but we adopt the second as it is a bit more straightforward. This is known as the Itô prescription.

Consider now a forced stochastic variable, \begin{equation} \dot{x}(t)=F(x)+\eta(t), \end{equation} where \(\eta(t)\) is the usual GRV forcing. We are interested in evaluating the time derivative of \(f\left(x(t)\right)\), a function of \(x(t)\). Formally, \begin{equation}\label{eq: def delta x} dx = x(t+\delta t)-x(t)=\int_t^{t+\delta t}ds\ F(x(s))+\int_t^{t+\delta t}ds\ \eta(s) \end{equation} \begin{equation} = F(x(t))dt + \mathcal{O}\left(\delta t^2\right) + d\eta(t).\nonumber \end{equation} We have compressed the second term into \(d\eta(t)\). From Eq. \ref{eq: eta covar problem}, we refer to \(d\eta(t)\) as proportational to the square root of \(\delta t\), \(d\eta(t)\sim\mathcal{O}\left(\sqrt{\delta t}\right)\). This will turn out to be the correction necessary to save the chain rule.

Now taking the derivative of \(f\left(x(t)\right)\), \begin{equation} \frac{d}{dt}\left[f\left(x(t)\right)\right]=\lim_{\delta t\rightarrow0}\left[\frac{f\left(x+\delta x\right)-f(x)}{\delta t}\right]. \end{equation} Now expanding \(f\) via Taylor series, \begin{equation} \frac{d}{dt}\left[f\left(x(t)\right)\right]=\lim_{\delta t\rightarrow0}\left[\frac{1}{\delta t}\sum_{k=1}^\infty\frac{(dx)^k}{k!}f^{(k)}(x)\right]. \end{equation} At order \(k=1\), we recover the usual chain rule, \begin{equation} \frac{1}{\delta t}\left(dx\partial_xf(x)\right)\rightarrow\dot{x}\partial_x f, \end{equation}

At order \(k=2\), something interesting happens when we insert our expression for \(dx\), Eq. \ref{eq: def delta x}, \begin{equation} \frac{1}{2\delta t}\left(dx^2\partial^2_xf(x)\right)=\frac{\partial^2_xf}{2\delta t}\left)(F(x)^2\delta t^2+2F(x)d\eta(t)\delta t+d\eta(t)^2+\mathcal{O}\left(\delta t^{5/2}\right)\right). \end{equation}

In the limit of \(\delta t\rightarrow 0\), all the terms vanish except for the \(d\eta(t)^2\) since this term is \(\mathcal{O}\left(\delta t\right)\). This leaves, \begin{equation} \frac{1}{2\delta t}\left(dx^2\partial^2_xf(x)\right)\rightarrow\frac{\partial_x^2f}{2}\left(\frac{d^2\eta}{dt^2}\right). \end{equation} Now assume the covariance of the GRV \(\eta(t)\) follows \(\langle\eta(t)\eta(t')\rangle=\sigma\delta(t-t')\), we say \(\frac{d^2\eta}{dt^2}=\sigma\) to represent the dispersion caused by the stochastic forcing. So, we obtain the Itô chain rule, \begin{equation} \frac{d}{dt}\left[f\left(x(t)\right)\right]=\dot{x}(t)f(x)+\frac{\sigma}{2}\partial_x^2f(x). \end{equation} We have “fixed” the chain rule and can make sense of the derivative of a function of a stochastic variable. Now we extend this perspective to a probability distribution.

The Fokker-Planck Equation

It turns out to be much more practical to compute not the path of a single stochastic variable, but the probability distribution of a stochastically forced variable in phase space. Consider \(\mathcal{P}\left(x,t|x_0,0\right)\), the probability to encounter our particle at location \(x\) at time \(t\) given that at time \(t=0\) the particle inhabited \(x_0\). This is a time-dependent probability distribution, so we can take a time derivative. First, notice we can rewrite \(\mathcal{P}\left(x,t|x_0,0\right)\) in terms of a delta function, \begin{equation} \mathcal{P}\left(x,t|x_0,0\right)=\int dy\ \mathcal{P}\left(y,t|x_0,0\right)\delta(y-x)=\left\langle\delta(y(t)-x)\right\rangle. \end{equation} Notice that \(y(t)\) carries the time dependence because we integrate over the variable \(y\). Taking the time derivative and using our Itô chain rule, \begin{equation} \frac{d}{dt}\left[\mathcal{P}\left(x,t|x_0,0\right)\right]=\left\langle\frac{d}{dt}\delta(y(t)-x)\right\rangle=\left\langle\dot{y}\partial_y\delta\left(y(t)-x\right)+\frac{\sigma}{2}\partial^2_y\delta\left(y(t)-x\right)\right\rangle. \end{equation} Using the evolution of the stochastic variable \(\dot{y}=F(y)+\eta(t)\), \begin{equation} \frac{d}{dt}\left[\mathcal{P}\left(x,t|x_0,0\right)\right]=\left\langle F(y)\partial_y\delta\left(y(t)-x\right)\right\rangle+\left\langle\partial_y\delta\left(y(t)-x\right)\eta(t)\right\rangle+\left\langle\frac{\sigma}{2}\partial^2_y\delta\left(y(t)-x\right)\right\rangle. \end{equation} The second term vanishes since it is proportional to \(\langle\eta\rangle=0\). The other two terms leave us with, \begin{equation} \frac{d}{dt}\left[\mathcal{P}\left(x,t|x_0,0\right)\right]=\int dy\ \left[\mathcal{P}\left(y,t|x_0,0\right)F(y)\partial_y\delta\left(y-x\right)\right]+\frac{\sigma}{2}\int dy\ \left[\mathcal{P}\left(y,t|x_0,0\right)\partial^2_y\delta\left(y-x\right)\right]. \end{equation} Because the probability distribution is properly normalized, we can freely apply the chain rule to remove the derivatives off the \(\delta\) functions, \begin{equation} \frac{d}{dt}\left[\mathcal{P}\left(x,t|x_0,0\right)\right]=-\int dy\ \left[\partial_y\left(\mathcal{P}\left(y,t|x_0,0\right)F(y)\right)\delta\left(y-x\right)\right]+\frac{\sigma}{2}\int dy\ \left[\partial^2_y\left(\mathcal{P}\left(y,t|x_0,0\right)\right)\delta\left(y-x\right)\right]. \end{equation} And evaluating the integrals, we obtain the Fokker-Planck Equation, \begin{equation}\label{eq: Fokker Planck} \frac{d}{dt}\left[\mathcal{P}\left(x,t|x_0,0\right)\right]=\partial_x\left(-F(x)+\frac{\sigma}{2}\partial_x\right)\mathcal{P}\left(x,t|x_0,0\right). \end{equation} This is just ordinary advection plus a spatial diffusion term arising from the stochastic forcing. Finally, we can return to our result of the fluid-immersed colloid.

The Spread of Equilibrium

Let’s return to our Langevin equation written in terms of position and momenta, \begin{equation}\label{eq: langevin two variable} \dot{x}=\frac{P}{M},q\quad\dot{P}=-\partial_xV(x)-\frac{\gamma}{M}P+\sqrt{\frac{2\gamma}{\beta}}\eta(t). \end{equation} Armed with the Fokker-Planck Equation, we wonder what the steady-state probability distribution looks like for a Langevin particle.

Our probability distribution is a function of the position, \(x\), momenta, \(P\), and time \(t\), \(\mathcal{P}\left(x,P,t\right)\). In higher dimensions, the Fokker-Planck Equation is simply summed over all the relevant variables, so applying Eq. \ref{eq: Fokker Planck} to Eq. \ref{eq: langevin two variable}, \begin{equation} \frac{d}{dt}\mathcal{P}=-\partial_x\left(\frac{P}{M}\mathcal{P}\right)-\partial_P\left(-(\partial_xV)\mathcal{P}-\frac{\gamma P}{M}\mathcal{P}-\frac{\gamma}{\beta}\partial_P\mathcal{P}\right) \end{equation} \begin{equation} =-\frac{P}{M}\partial_x\mathcal{P}+\partial_P\left((\partial_xV)\mathcal{P}+\frac{\gamma P}{M}\mathcal{P}+\frac{\gamma}{\beta}\partial_P\mathcal{P}\right). \end{equation} This is known as Kramer’s equation. We are looking for a steady state solution satisfying \(\frac{d}{dt}\mathcal{P}=0\). To solve, we note that the differential equation looks fairly separable. Substitute the ansatz, \(\mathcal{P}=\exp\left[\mathcal{A}(x)+\mathcal{B}(P)\right]\), \begin{equation}\label{eq: kramers eq} -\frac{P}{M}\partial_x\mathcal{A}+(\partial_xV)\partial_P\mathcal{B}+\frac{\gamma}{M}+\frac{\gamma P}{M}\partial_P\mathcal{B}+\frac{\gamma}{\beta}\partial_P^2\mathcal{B}+\frac{\gamma}{\beta}\left(\partial_P\mathcal{B}\right)^2=0. \end{equation} Comparing the first two terms yields, \begin{equation} \frac{P}{M}\partial_x\mathcal{A}=(\partial_xV)\partial_P\mathcal{B}\rightarrow\mathcal{A}(x)=CV(x),\quad\mathcal{B}(P)=C\frac{P^2}{2M}, \end{equation} where \(C\) is an undetermined constant. Using this solution for \(\mathcal{B}\) in the remaining terms of Eq. \ref{eq: kramers eq}, \begin{equation} \frac{\gamma}{M}+C\frac{\gamma P^2}{M^2}+C\frac{\gamma}{M\beta}+C^2\frac{\gamma P^2}{\beta M^2}=0, \end{equation} and we see \(C=-\beta\). Therefore, the solution for the probability distribution is, \begin{equation} \mathcal{P}\left(x,P,t\right)=\exp\left[-\beta\left(\frac{P^2}{2M}+V(x)\right)\right]=\exp\left[-\beta H_C\right], \end{equation} where \(H_C\) is the Hamiltonian of the undisturbed colloid. The steady state probability distribution of the colloid in a fluid is the Boltzmann distribution for the colloid Hamiltonian! This is a very interesting result, as it means that the equilibrium of the fluid is “contagious.” Because the fluid is at equilibrium, the colloid will move towards equilibrium as well, towards the steady-state Boltzmann distribution.

Sources

This derivation largely follows from the lecture notes of my Statistical Dynamics II, 8.08, class taught by Professor Julien Tailleur during the 2026 IAP period.