12.3: Diffusion in a Potential

Fokker–Planck Equation

Diffusion with drift or diffusion in a velocity field is closely related to diffusion of a particle under the influence of an external force f or potential U.

\[ f(x) = - \dfrac{\partial U}{\partial x} \nonumber \]

When random forces on a particle dominate the inertial ones, we can equate the drift velocity and external force through the friction coefficient

\[ \begin{aligned} &\cancel{m\ddot{x}} = f_d +\cancel{f_r(t)} +f_{ext} \\ &f_d = -\zeta v_x \\ &f_{ext} = \zeta v_x \end{aligned}\]

\[ f= \zeta v_x \]

and therefore the contribution of the force or potential to the total flux is

\[ J_U = v_xC = \dfrac{f}{\zeta} C = -\dfrac{C}{\zeta} \dfrac{\partial U}{\partial x} \]

The Fokker–Planck equation refers to stochastic equations of motion for the continuous probability density \(\rho (x,t)\) with units of m⁻¹. The corresponding continuity expression for the probability density is

\[ \dfrac{\partial \rho}{\partial t} = -\dfrac{\partial j}{\partial x} \nonumber \]

where j is the flux, or probability current, with units of s^–1, rather than the flux density we used for continuum diffusion J (m⁻² s⁻¹). If the concentration flux is instead expressed in terms of a probability density eq. (12.1.3) becomes

\[ j = -D \dfrac{\partial \rho}{\partial x} + \dfrac{f(x)}{\zeta}\rho \]

and the continuity expression is used to obtain the time-evolution of the probability density:

\[\dfrac{\partial \rho}{\partial x} = D\dfrac{\partial^2 \rho}{\partial x^2} - \dfrac{\partial}{\partial x} \left[ \dfrac{f(x)}{\zeta} \rho \right] \]

This is known as a Fokker–Planck equation.

Smoluchowski Equation

Similarly, we can express diffusion in the presence of an internal interaction potential U(x) using eq. (12.3.2) and the Einstein relation

\[ \zeta = \dfrac{k_BT}{D} \]

Then the total flux with contributions from the diffusive flux and potential flux can be written as

\[ J=-D\dfrac{\partial C}{\partial x}- \dfrac{DC}{k_BT} \left( \dfrac{\partial U}{\partial x} \right) \]

and the corresponding diffusion equation is

\[ \dfrac{\partial C}{\partial t} = D \left[ \dfrac{\partial^2 C}{\partial x^2} - \dfrac{\partial}{\partial x} \left[ \dfrac{C}{k_BT} \left( \dfrac{\partial U}{\partial x} \right) \right] \right] \]

This is known as the Smoluchowski Equation.

Linear Potential

For the case of a linear external potential, we can write the potential in terms of a constant external force \(U=-f_{ext}x\). This makes eq. (12.3.7) identical to eq. (12.1.3), if we use eqs. (12.3.1) and (12.3.5) to define the drift velocity as

\[ v_x = \dfrac{f_{ext}D}{k_BT} \equiv \underset{sim}{f} D \nonumber \]

\[ J = -D \dfrac{\partial C}{\partial x} + \underset{\sim}{f} DC \nonumber \]

Here I defined \(\underset{\sim}{f}\) as the constant external force expressed in units of k_BT.

The probability distribution that describes the position of particles released at x₀ after a time t is

\[ P(x,t) = \dfrac{1}{\sqrt{4\pi Dt}} \exp \left[ -\dfrac{(x-x_0-\underset{\sim}{f}Dt)^2}{4Dt} \right] \nonumber \]

As expected, the mean position of the diffusing particle is given by ⟨x(t)⟩ = x₀ + v_xt.

To make use of this, let’s calculate the time it takes a monovalent ion to diffuse freely across the width of a membrane (d) under the influence of a linear electrostatic potential of Φ = 0.3V. With U = eΦ

\[ t = \dfrac{d}{v_x}= \dfrac{k_BTd}{f_{ext}D} = \dfrac{k_BTd^2}{e\Phi D} \nonumber \]

Using d = 4 nm, D = 10−5 cm²/s, and e = 1.6×10⁻¹⁹ C, we obtain t = 1.4 ns.

Steady‐State Solutions

For steady-state solutions to the Fokker–Planck or Smoluchowski equations, we can make use of a commonly used mathematical manipulation. As an example, let’s work with eq. (12.3.3), re-writing it as

\[ j = -D \left[ \dfrac{\partial \rho}{\partial x} -\dfrac{\rho}{k_BT} \left( \dfrac{\partial U}{\partial x} \right) \right] \]

We can rewrite the quantity in brackets as:

\[ e^{-U(x) /k_BT} \dfrac{d}{dx} \left[ \rho e^{U(x)/k_BT} \right] \nonumber \]

Separating variables, we obtain

\[ - \dfrac{j}{D} e^{U(x) /k_BT} dx = d(\rho e^{U(x)/k_BT} \nonumber \]

This is an expression that can be manipulated in various ways and integrated over different boundary conditions.¹ For instance, recognizing that j is a constant under steady state conditions, and integrating from x to a boundary b:

\[ \begin{aligned} -\dfrac{j}{D} \int^b_x e^{U(x)/k_BT} dx &= \int^b_x d(\rho e^{U(x)/k_BT}) \\ &= \rho (b) e^{U(b)/k_BT} - \rho (x)e^{U(x)/k_BT} \end{aligned} \]

This leads one to an important expression for the steady state flux in the diffusive limit:

\[ j = \dfrac{-D\left[ \rho (b) e^{U(b)/k_BT}-\rho (x) e^{U(x)/k_BT} \right]}{\int^b_x e^{U(x)/k_BT}dx} \nonumber \]

The boundary chosen depends on the problem, for instance b is set to infinity in diffusion to capture problems or set as a fixed boundary for first-passage time problems.
For problems involving an absorbing boundary condition, ρ(b) = 0, and we can solve for the probability density as

\[ \rho (x) = \dfrac{j}{D} e^{-U(x)/k_BT} \left[ \int^b_x e^{U(x')/k_BT} dx' \right] \nonumber \]

If we integrate both sides of this expression over the entire space, the left hand side is just unity, so we can express the steady-state flux as

\[ j = D^{-1} \left[ \int^b_0 e^{-U(x)/k_BT} \left[ \int^b_x e^{U(x')/k_BT}dx' \right] dx \right]^{-1} \nonumber \]

____________________________________________

1. The general three-dimensional expression is \( \textbf{J}(\textbf{r},t)= -De^{-U(\textbf{r})/k_BT}\nabla \cdot [ e^{U(\textbf{r})/k_BT}\rho (\textbf{r},t) ] \).