In the year of 1827, a Scottish botanist named Robert Brown made a series of observations of pollen grains under his microscope [1]. The pollen grains were immersed in water and he discovered that their motion was of a rather erratic character, however he was not able to come up with a theory that could explain these observations. Today we know that this random motion is caused by a vast number of collisions between the grains and the water molecules, and the phenomenon is called *Brownian motion* in honour of Robert Brown. The first person to come up with a proper theoretical treatment of Brownian motion was Albert Einstein, who in 1905 published a paper providing a quantitative theory of the phenomenon (*Über die von der molekularkinetischen Theorie der Wärme geforderte Bewegung von in ruhenden Flüssigkeiten suspendierten Teilchen*).

We consider a very simple probability model of the random motion of a Brownian particle proposed by M. Smoluchowski [2]. Nevertheless, this model gives rise to some interesting physics, as we are about to discover soon. The time evolution of the position of the particle takes the form of discrete steps, and in each step there is an equal probability of the particle moving a distance $h$ to the left or a distance $h$ to the right. The total number of steps taken is denoted by $N = N_l + N_r$, where $N_l$ is the number of steps to the left and $N_r$ the number of steps to the right. The probability that the particle moves $n$ steps to the right is thus given by the binomial distribution

\begin{align*} \mathrm{P}_N(N_r = n) &= \binom{N}{n} \left(\frac{1}{2}\right)^n \left(1 - \frac{1}{2}\right)^{N-n} \\ &= \binom{N}{n} \left(\frac{1}{2}\right)^N \quad. \end{align*}

The position of the particle after $N$ steps is equal to the distance travelled to the right minus the distance travelled to the left, assuming that the particle is initially located at $x = 0$

\begin{align*} x &= (N_r - N_l)\,h \\ &= (2N_r - N)\,h \quad, \end{align*}

where we have used that $N_l = N - N_r$. Let us plot the probability distribution for a few different values of $N$.

In [1]:

```
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from scipy.special import binom
```

In [2]:

```
# Set common figure parameters
newparams = {'font.size': 14, 'lines.linewidth': 1.0, 'figure.figsize': (16, 10),
'legend.frameon': True}
plt.rcParams.update(newparams)
```

In [3]:

```
def calc_pos_distribution(N=10, p=0.5, h=1.0):
""" Creates an array containing the probabilities associated with each position
of the particle after a number N steps.
Parameters:
N: int (non-negative). Number of steps in the time evolution.
p: float (0 <= p <= 1). The probability that the particle moves one step to the right.
h: float. The distance the particle moves during each step.
Returns: (xs, ps)
xs: float array. The possible positions of the particle.
ps: float array. The probabilities of each possible position of the particle.
"""
ns = np.arange(0, N+1)
ps = binom(N, ns) * p**ns * (1 - p)**(N - ns)
xs = 2*ns - N
return (xs, ps)
```

In [4]:

```
fig = plt.figure()
fig.suptitle("Probability distributions of the position of a Brownian particle after $N$ steps with $h = 1$.")
# Handles and data arrays for four subplots
axes = [None] * 4
pos_arrays = [None] * 4
prob_arrays = [None] * 4
# Iterate over each subplot
for i in range(0, 4):
N = 10 * i
pos_arrays[i], prob_arrays[i] = calc_pos_distribution(N)
axes[i] = fig.add_subplot(2, 2, i+1)
axes[i].bar(pos_arrays[i], prob_arrays[i], width=2.0, align='center')
axes[i].set_ylim([0, 1])
axes[i].set_xlim([-20, 20])
axes[i].set_xlabel('$x$')
axes[i].set_ylabel('P($x$)')
axes[i].legend(['$N$ = {0}'.format(N)])
```

The expected value of the position of the particle remains at $x = 0$ for all $N$, but as time evolves ($N$ grows bigger), the probability distribution spreads out more and more. This closely resembles the macroscopic phenomenon of *diffusion*, where particles move from a region of high concentration to a region of lower concentration.
To see this mathematically, we investigate the limit of the binomial distribution where it approaches a continuous probability distribution: Let $N$ increase and let the step size $h$ decrease in such a way that $hN$ is kept constant.
The number of steps $N$ is related to time $t$ as follows: $N = t/\tau$, where $\tau$ is the average time between collisions of the Brownian particle and a water molecule.

The result is a normal distribution when $N$ gets large (a derivation is given in the appendix)

\begin{equation*} \mathrm {P}(x) = \frac{1}{\sqrt{2\pi h^2 t/\tau}} \exp\left({-\frac{1}{2} \frac{x^2}{h^2 t/\tau}}\right) \quad . \end{equation*}

The expected value of the position is $x = 0$ at all times, and we have a time-dependent variance $\sigma^2 = h^2 t/\tau$. If we introduce the constant $D \equiv h^2/\tau$, the normal distribution can be written in the following form

\begin{equation*} \mathrm{P}(x) = \frac{1}{\sqrt{4\pi Dt}} \exp\left({-\frac{1}{2} \frac{x^2}{2Dt}}\right) \quad . \end{equation*}

Suppose we have a macroscopic system of $\mathcal{N}$ independent particles. The concentration of the particles can then be written as

\begin{equation*} n(x, t) = \mathcal{N} \, \mathrm{P}(x) \quad . \end{equation*}

A straight-forward calculation shows that $n(x, t)$ satisfies the one-dimensional diffusion equation

\begin{equation*} \frac{\partial}{\partial t} n(x, t) = D \frac{\partial^2}{\partial x^2} n(x, t) \quad , \end{equation*}

hence our simple microscopic model gives rise to the macroscopic phenomenon of diffusion!

In order to further acknowledge the connection between Brownian motion and diffusion, it would be helpful to visualise the time development of a system of independent Brownian particles. By independent we mean that the motion of a Brownian particle does not influence the motion of the others. Suppose we have a collection of $\mathcal{N}$ Brownian particles that can move in two dimensions as described in the model above. We are going to plot the positions of the particles after each time step and indeed observe that they tend to move from regions of high concentration to regions of lower concentration.

During the simulation we are also going to compute and display the arithmetic mean and the sample variance of the positions

\begin{align*} \overline{x} &= \frac{1}{\mathcal{N}} \sum_{k = 1}^\mathcal{N} x_k \\ s_x^2 &= \frac{1}{\mathcal{N}} \sum_{k = 1}^\mathcal{N} (x_k - \overline{x})^2 \quad . \end{align*}

The expected values of these quantities are

\begin{align*} \langle\overline{x}\rangle &= \frac{1}{\mathcal{N}} \sum_{k = 1}^\mathcal{N} \langle x_k \rangle \\ &= 0 \\ \langle s_x^2 \rangle &= \frac{1}{\mathcal{N}} \sum_{k = 1}^\mathcal{N} \langle(x_k - \overline{x})^2 \rangle \\ &= \frac{1}{\mathcal{N}} \sum_{k = 1}^\mathcal{N} \left(\langle x_k^2 \rangle - \langle \overline{x}\rangle^2\right) \\ &= \frac{1}{\mathcal{N}} \sum_{k = 1}^\mathcal{N} \langle x_k^2 \rangle \\ &= \frac{1}{\mathcal{N}} \sum_{k = 1}^\mathcal{N} h^2 N \\ &= h^2 N \quad . \end{align*}

In [5]:

```
from matplotlib import animation
from IPython.display import HTML
plt.rcParams.update({'animation.html':'html5', 'savefig.dpi': 60})
def populate_particles(nParticles, R=0):
""" Creates two arrays containing randomised positions of the particles.
The positions are randomised within a circle of radius R about the origin.
Parameters:
nParticles: int (positive). The number of Brownian particles.
R: float (positive). The radius within which the randomised positions are confined.
Returns: (xs, ys)
xs: float array. The x-coordinates of the particles.
ys: float array. The y-coordinates of the particles.
"""
angles = 2*np.pi*np.random.rand(nParticles)
radii = R*np.random.rand(nParticles)
xs = radii*np.cos(angles)
ys = radii*np.sin(angles)
return (xs, ys)
def evolve_particles(xs, ys, h=0.1, pr=0.5, pu=0.5, vx=0, vy=0):
""" Evolves the positions of the particles by one time step.
Parameters:
xs: float array. The x-coordinates of the particles.
ys: float array. The y-coordinates of the particles.
h: float. The distance a particle moves during a time step.
pr: unsigned int (0 <= pr <= 1). The probability that a particle moves to the right along the x-axis.
pu: unsigned int (0 <= pr <= 1). The probability that a particle moves upwards along the y-axis.
vx: float. The x-component of the average drift velocity (length h per unit time).
vy: float. The y-component of the average drift velocity (length h per unit time).
Returns: None. The updated positions are stored in xs and ys.
"""
if xs.size != ys.size:
print("Error: The sizes of xs and ys do not match.")
return -1
randnumbers = np.random.random(xs.size)
for i, step_right in enumerate(randnumbers < pr):
if step_right:
xs[i] += h
else:
xs[i] -= h
randnumbers = np.random.random(ys.size)
for i, step_up in enumerate(randnumbers < pu):
if step_up:
ys[i] += h
else:
ys[i] -= h
# Particle drift
xs += vx*h
ys += vy*h
def init_anim():
""" Initialises the animation.
"""
global ax, line, textbox_sim, textbox_theory
line, = ax.plot([], [], 'bo')
ax.set_xlim([-20, 20])
ax.set_ylim([-20, 20])
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_title('Numerical simulation of Brownian motion')
props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
# A text box that will display the average and the variance of the positions of the simulated particles.
textbox_sim = ax.text(0.05, 0.95, '', transform=ax.transAxes, fontsize=12,
verticalalignment='top', bbox=props)
# A text box that will display the average and the variance of the theoretical normal distribution of the positions.
textbox_theory = ax.text(0.775, 0.95, '', transform=ax.transAxes, fontsize=12,
verticalalignment='top', bbox=props)
return line, textbox_sim, textbox_theory
def animate(N, xs, ys, h, nParticles, drift_vx=0, drift_vy=0):
""" Runs the simulation and updates the plot.
Parameters:
N: int. The number of time steps taken so far.
xs: float array. The x-coordinates of the particles.
ys: float array. The y-coordinates of the particles.
h: float. The distance a particle moves during a time step.
nParticles: int (positive). The number of Brownian particles.
Returns: Handles to the objects that are going to be updated in the plot.
"""
global ax, line, textbox_sim, textbox_theory
# Evolve the system by one time step.
evolve_particles(xs, ys, h, vx=drift_vx, vy=drift_vy)
N += 1
# Update the plot of the positions.
line.set_data(xs, ys)
# Compute and display quantities associated with the simulated particles.
x_avg, y_avg = np.average(xs), np.average(ys)
x_var = np.sum((xs - x_avg)**2) / nParticles
y_var = np.sum((ys - y_avg)**2) / nParticles
x_avg_str = '$\\langle x \\rangle = %+.2f$' % (x_avg)
y_avg_str = '$\\langle y \\rangle = %+.2f$' % (y_avg)
x_var_str = '$\mathrm{Var}(x) = %.2f$' % (x_var)
y_var_str = '$\mathrm{Var}(y) = %.2f$' % (y_var)
textbox_sim.set_text('Simulation\n' + x_avg_str + '\n' + y_avg_str + '\n' + x_var_str + '\n' + y_var_str)
# Compute and display the theoretical quantities.
x_var = h**2 * N
y_var = h**2 * N
x_avg_str = '$\\langle x \\rangle = %.2f$' % (drift_vx * h * N)
y_avg_str = '$\\langle y \\rangle = %.2f$' % (drift_vy * h * N)
x_var_str = '$\mathrm{Var}(x) = %.2f$' % (x_var)
y_var_str = '$\mathrm{Var}(y) = %.2f$' % (y_var)
textbox_theory.set_text('Theoretical\n' + x_avg_str + '\n' + y_avg_str + '\n' + x_var_str + '\n' + y_var_str)
return line, textbox_sim, textbox_theory
```

In [6]:

```
# The number of Brownian particles.
nParticles = 500
# The total number of steps taken during the simulation.
nSteps = 1000
# Create arrays containing the x and y positions of the particles.
xs, ys = populate_particles(nParticles)
# The distance that the particles move during each time step.
h = 0.1
# Run the simulation and visualise the system as an animation.
fig, ax = plt.subplots()
h_anim = animation.FuncAnimation(fig, animate, init_func=init_anim, frames=nSteps, interval=50, blit=True,
fargs=(xs, ys, h, nParticles))
plt.close(h_anim._fig)
HTML(h_anim.to_html5_video())
```

Out[6]: