Langevin equation and Brownian motion

import numpy as np
import matplotlib.pyplot as plt
from numpy.linalg import norm
from scipy.constants import c,epsilon_0,e,physical_constants

%config InlineBackend.figure_format = 'retina' 

Langevin Equation

Langevin theory provides now an equation of motion for objects whoch are subject to fluctuating forces. In the most general way, Newtons equation of motion is written as

\[ m\frac{d\vec{v}}{dt}=\vec{F}(t) \]

and contains on the right side the sum of all forces \(\vec{F}\) acting on a particle. This total force can be separated into

• a viscous drag force, e.g. \(6\pi \eta R \vec{v}\) for a spherical particle

• a force coming from and external potential \(U\)

• random fluctuating force \(\vec{\zeta}\) due to the collisions with the solvent

With the force arising from an external potential we find the Langevin equation

\[ m\frac{d\vec{v}}{dt}=-6\pi \eta R \vec{v}-\nabla U +\vec{\zeta}(t) \tag{Langevin Equation} \]

The time dependence of the fluctuating force is not known in detail. Thus we can only refer to the statistical properties (e.g. its moments) when solving the Langevin equation. Therefore also not analytical solution for the position can be written down, only average proerties over an ensemble or over time (given ergodicity).

The first moment (the mean) of the fluctuating force gives

\[ \langle\vec{\zeta}(t) \rangle =0 \]

i.e. there is no net force acting on average on the particle considered. Yet the mean velocity of the particle is not zero, \(\langle \vec{v}\rangle\neq 0\). Considering the mean velocity at short times, for example, yields

\[ m\frac{d\langle \vec{v}\rangle}{dt}=-6\pi \eta R \langle \vec{v}\rangle \]

giving

\[ \langle \vec{v}\rangle=\vec{v}(0)\exp\left ( -\frac{t}{\tau_{mr}}\right) \]

with \(\tau_{mr}\) is a characteristic relaxation time for velocity of brownian particle. We can further continue to derive the mean squared position by taking the scalar product of the Langevin equation with \(\vec{r}\) and neglecting the potential. Taking the ensemble average yields

\[ m\frac{d}{dt}\langle \vec{r}\cdot \vec{v}\rangle =-6\pi \eta R \langle \vec{r}\cdot \vec{v}\rangle + m\langle v^2\rangle \]

since

\[ \vec{r}\cdot \frac{d\vec{v}}{dt}=\frac{d}{dt}(\vec{r}\cdot \vec{v})-v^2 \]

and \(\langle \vec{r}\cdot \vec{\zeta}\rangle =0\). Using our previous result for the mean squared velocity, i.e. \(\langle v^2\rangle =3k_B T\) we obtain after integrating over time

\[ \langle \vec{r}\cdot \vec{v}\rangle =C\exp\left ( -\frac{t}{\tau_{mr}}\right)+ \frac{k_B T}{2\pi\eta R} \]

Using the initial condition \(\vec{r}(t=0)=0\) this gives

\[ \langle \vec{r}\cdot \vec{v}\rangle =\frac{k_B T}{2\pi\eta R}\left (1 -\exp\left ( -\frac{t}{\tau_{mr}}\right)\right) \]

A second integration over time of \(\langle \vec{r}\cdot \vec{v}\rangle\) can be carried out using

\[ \langle \vec{r}\cdot \vec{v}\rangle=\frac{1}{2}\frac{d}{dt}\langle r^2\rangle \]

and we finally obtain the mean squared displacement

\[ \langle r^2\rangle =\frac{k_B T}{\pi \eta R}\left [ t-\tau_{mr}\left ( 1-\exp\left (-\frac{t}{\tau_{mr}}\right )\right )\right ] \]

This is the mean squared displacement for an Ornstein Uhlenbeck process, a process, which has a persistence for a certain time and then is randomized. It described the transition from the ballistic to diffusive regime for a Brownian particle except that it misses the hydrodynamic memory in the intermediate regime.

t=np.logspace(-15,0,1000)
k_B=1.381e-23 #J/K
T=300 #K
eta=1e-3 #mPa s
R=200e-9 #m
gamma=6*np.pi*eta*R
k=1e-9 #N/m

def MSD(t,taumr):
    return(k_B*T/(np.pi*R)*(t-taumr*(1-np.exp(-t/taumr))))

fig=plt.figure()
plt.ion()
plt.loglog(t,MSD(t,1e-9),label=r"$\tau_{mr}=1e-9$ s")
plt.axvline(x=1e-9,linestyle="--")
plt.xlabel(r"$t$")
plt.ylabel(r"$<r^2(t)$")
plt.legend()
plt.tight_layout()

We can look at different limits of the process with respect to the momentum relaxation time, ie.

1) long time limits: \(t\gg \tau_{mr}\)

In this regime, the exponential function has already decayed to zero and we obtain

\[ \langle r^{2}\rangle =6Dt \]

with

\[ D=\frac{k_B T}{6\pi\eta R} \]

1) short times: \(t\approx \tau_{mr}\)

For short times, we can expand the exponential function up to second order

\[ \exp\left ( -\frac{t}{\tau_{mr}}\right )\approx 1-\frac{t}{\tau_{mr}}+\frac{t^2}{\tau_{mr}^2} \]

to obtain for

\[ t-\tau_{mr}\left ( 1-\exp\left (-\frac{t}{\tau_{mr}}\right )\right )\approx \frac{t^2}{2\tau_{mr}} \]

and for the mean squared displacement

\[ \langle r^2 \rangle \propto t^2 \]

which is the result we expect for ballistic motion. Between both regimes, the MSD could be complicated as the short time motion starts a hydrodynamic flux in the environment which is influencing the particle motion again. These hydrodynamic memory effects cause the so called long time tails in the mean squared displacement.

Fluctuation-dissipation theorem

While we so far only know that \(\langle\vec{\zeta}(t) \rangle =0\), we also know that there must be a certain magnitude of the random force that is connected to the temperature of the sample. The short time autocorrelation of the noise should be of the type

\[ \langle\vec{\zeta}(t_1)\cdot \vec{\zeta}(t_2) \rangle =A\delta (t_1-t_2) \]

stating the the forces are truly random and only correlated, when the times conincide. We would like to determine the prefactor \(A\) now. We can write down the Langevin equation with

\[ \frac{d}{dt}\vec{v}(t)=-\frac{1}{\tau_{mr}}\vec{v}(t)+\frac{1}{m}\vec{\zeta}(t) \]

and multiply by \(\exp(t/\tau_{mr})\), which results in

\[ \exp\left (\frac{t}{\tau_{mr}}\right )\frac{d}{dt}\vec{v}(t)+\frac{1}{\tau_{mr}}\exp\left (\frac{t}{\tau_{mr}}\right )\vec{v}(t)=\frac{d}{dt}\left (\exp\left (\frac{t}{\tau_{mr}}\right )\vec{v}(t) \right )=\exp\left(\frac{t}{\tau_{mr}}\right )\frac{1}{m}\vec{\zeta}(t) \]

Integrating both sides of the two right parts over time results in

\[ \vec{v}(t)=\vec{v}(0)\exp \left (-\frac{t}{\tau_{mr}} \right )+\exp \left (-\frac{t}{\tau_{mr}} \right )\frac{1}{m}\int_0^t \exp \left (\frac{t'}{\tau_{mr}} \right )\vec{\zeta}(t')dt' \]

which is the time evolution of the velocity. To obtain \(\langle\vec{\zeta}(t_1)\cdot \vec{\zeta}(t_2) \rangle\), we need to calculate the product of the velocity at two times and take the ensemble average. The calculation is essentially a writing exercise which at the end yields three different terms.

The first term contains

\[ \langle v_0^2 \rangle \exp\left (-\frac{2t}{\tau_{mr}} \right )\rightarrow 0 \]

which decays to zero for long times. The second term contains two mixed terms of the type

\[ \frac{v_0}{m}\exp\left (-\frac{t}{\tau_{mr}} \right ) \int \ldots \rightarrow 0 \]

which also decay to zero for long times. The last term contains

\[ \frac{1}{m^2}\int_{-\infty}^{t}dt' \int_{-\infty}^{t}dt''\exp\left ( -\frac{t-t'}{\tau_{mr}} \right )\exp\left ( -\frac{t-t''}{\tau_{mr}} \right )\langle \vec{\zeta}(t')\vec{\zeta}(t'')\rangle \]

whcih is the only nonzero term. Inserting our initial assumption \(\langle\vec{\zeta}(t_1)\cdot \vec{\zeta}(t_2) \rangle =A\delta (t_1-t_2)\) results in the mean squared velocity

\[ \langle v^2\rangle =\frac{A}{2m}\int_{-\infty}^{t}dt' \exp\left ( -\frac{2(t-t')}{\tau_{mr}}\right)=\frac{A\tau_{mr}}{2m^2} \]

Since the velocity also has to comply with equipartition, i.e.

\[ \langle v^2\rangle=\frac{k_B T}{m}=\frac{A\tau_{mr}}{2m^2} \]

we find

\[ A=2\gamma k_B T \]

with \(\gamma=6\pi\eta R\) or

\[ \langle\vec{\zeta}(t_1)\cdot \vec{\zeta}(t_2) \rangle =2\gamma k_B T\delta (t_1-t_2) \]

which is a very fundamental result. It states that the fluctuating force \(\vec{\zeta}(t)\) are related to the viscous forces expressed by the friction coefficient \(\gamma\). This makes sensse, since the friction forces must dissipate energy into heat. The motion would thus come to a rest if heat and dynamics would be decoupled. The law we obtained is thus stateing the dissipated energy goes back again into kinetic energy. This is the statement of the fluctuation dissipation theorem,

Simulating Langevin dynamics

• Wiki on langevin dynamics

• Wiki on Brownian dynamics

import numpy as np

def langevin(z, v, tau, eta, dt):
    # Implements numerical solution by means of the Euler-Maruyama method of the Langevin which is known as an SDE (Stochastic Differential equation)
    # dz = v * dt
    # dv = -1/tau * v * dt + eta dW
    
    dW = np.random.normal(loc = 0, scale = np.sqrt(dt), size = z.shape)

    z_new = z + v*dt
    v_new = v - (1/tau)*v*dt + eta*dW
    
    return z_new, v_new

# Numerical parameters and system coefficients
Np   = 10000
Tmax = 1000
tau  = 10
eta  = 0.01
dt   = 0.1
Nt   = int(Tmax/dt) + 1
# Output interval
Nskip = 10

# Initial values
Z   = np.zeros(Np)
V   = np.zeros(Np)
time, var = [], []

for i in range(1, Nt):

    Z, V = langevin(Z, V, tau, eta, dt)
    # write time and variance every Nskip timesteps
    if i % Nskip == 0:
        time.append(i*dt)
        var.append(np.var(Z))

plt.loglog(time, var)
plt.xlabel(r"$t$")
plt.ylabel(r"$<r^2(t)>$")
plt.legend()
plt.tight_layout()