Brownian Motion Simulation

What the simulation shows

In 1827 the botanist Robert Brown saw particles released from pollen grains jitter endlessly under his microscope. The lab reproduces that experiment in a computer: N particles start at the origin and, every Δt seconds, each coordinate receives an independent random kick drawn from a normal distribution with mean 0 and variance 2DΔt.

• Trails. No two paths look alike, and no path settles down. The dashed circle has the theoretical root-mean-square radius √(4Dt); the green circle is the measured one. Both grow like √t, so doubling the distance takes four times as long.
• Mean squared displacement (MSD). Averaged over particles, the squared distance from the start follows the straight line 4Dt. The dotted band shows how far an average over N particles is expected to stray: the squared distance of a single particle is exponentially distributed, so its standard deviation equals its mean and the average has relative spread 1/√N.
• Positions and increments. The single-step increments are normal by construction. The positions are sums of many steps and are normal too, with a variance that grows from 2DΔt to 2DT.

Try a single particle: the MSD curve becomes a wild, jagged line that only follows 4Dt on average. Then try 1,000 particles, or change Δt by a factor of 100 and note that the MSD does not care.

Einstein’s relation

In 1905 Albert Einstein showed that the random kicks of molecules make suspended particles diffuse, and that the spread of their positions grows linearly in time. In d dimensions:

He also connected D to the temperature T, the viscosity η of the fluid and the particle radius a through the Stokes–Einstein relation

For a sphere 1 µm across in water at room temperature this gives D ≈ 0.4 µm²/s, which is why the lab measures lengths in micrometers and times in seconds. Measuring the MSD of real particles and inverting the relation yields Boltzmann's constant, and with it Avogadro's number; Jean Perrin's measurements along these lines helped settle the reality of atoms.

The lab estimates D the same way an experimenter would: a least-squares line through the origin, D̂ = Σt_k·MSD_k / (4Σt_k²). Successive MSD values share the same particles and are strongly correlated, so the fit is not much more precise than the final MSD alone.

From random walk to Wiener process

Take the simple random walk, shrink the steps to Δx and the time per step to Δt with Δx²/Δt fixed, and the path converges to the Wiener process W(t), the mathematical Brownian motion. It is defined by three properties: W(0) = 0, increments over disjoint time intervals are independent, and W(t) − W(s) is normal with variance t − s. A particle with diffusion coefficient D moves as √(2D)·W(t) in each coordinate.

Those properties give an exact simulation method. Unlike an ODE solver, stepping with normal increments has no discretization error: the sampled positions have exactly the right joint distribution at every multiple of Δt. The time step only decides how finely the path is looked at. That is also the only way to look at it, because a Brownian path is continuous but nowhere differentiable. Zoom in on any piece and it is as rough as the whole, since W(ct) has the same distribution as √c·W(t). Its velocity is undefined and its length over any interval is infinite.

Model versus physics

The Wiener process is a model of the overdamped limit. A real particle has mass: after a kick it keeps moving for a velocity-relaxation time m/(6πηa), well under a microsecond for a micrometer bead in water. On shorter time scales its motion is smooth (ballistic), and the Langevin and Ornstein–Uhlenbeck models describe the crossover. For the millisecond-to-second scales of the lab, the overdamped model is excellent.

The simulation also leaves out interactions between particles, walls and hydrodynamic coupling, and it assumes a uniform temperature and viscosity so that D is constant. The two coordinates are independent, which describes motion within a plane, or the 2-D projection of 3-D motion seen through a microscope.

Example result

With 100 particles, D = 1 µm²/s and Δt = 0.01 s for 10 s (seed 42), the mean squared displacement at the end is 38.267 µm² against 4DT = 40.000 µm². The expected spread of that average is 4.000 µm², so the difference is ordinary. A least-squares fit to the whole MSD curve gives D̂ = 1.009 µm²/s, and each step moves a particle 0.2000 µm on average (root mean square).

Implementation

The page runs simulation.js from the package below; the Python file is a line-by-line port with the same random-number generator, so a given seed produces the same result in both languages. The excerpts show the core of each.

// Brownian motion sampled every Δt: each coordinate receives an independent
// Gaussian increment with mean 0 and variance 2DΔt. This is exact at the sample
// times (no discretization error), because Brownian increments are Gaussian.
export function step(s, D, dt, rng, onIncrement) {
  const sigma = Math.sqrt(2 * D * dt);
  let sum = 0;
  for (let i = 0; i < s.n; i++) {
    const dx = sigma * rng.normal();
    const dy = sigma * rng.normal();
    s.x[i] += dx;
    s.y[i] += dy;
    sum += s.x[i] * s.x[i] + s.y[i] * s.y[i];
    if (onIncrement) onIncrement(dx, dy);
  }
  s.t++;
  s.msd = sum / s.n; // mean squared displacement from the start
}

// Least-squares fit of MSD(t) = 4Dt through the origin (Einstein, 2-D):
// minimizing Σ (MSD_k − 4D t_k)² gives D = Σ t_k·MSD_k / (4 Σ t_k²).
export function fitD(times, msd) {
  let num = 0, den = 0;
  for (let k = 0; k < times.length; k++) {
    num += times[k] * msd[k];
    den += times[k] * times[k];
  }
  return den > 0 ? num / (4 * den) : 0;
}

def step(s, D, dt, rng):
    """Add an independent N(0, 2*D*dt) increment to each coordinate of every particle."""
    sigma = math.sqrt(2 * D * dt)
    x, y = s["x"], s["y"]
    total = 0.0
    for i in range(s["n"]):
        dx = sigma * rng.normal()
        dy = sigma * rng.normal()
        x[i] += dx
        y[i] += dy
        total += x[i] * x[i] + y[i] * y[i]
    s["t"] += 1
    s["msd"] = total / s["n"]  # mean squared displacement from the start


def fit_d(times, msd):
    """Least-squares fit of MSD(t) = 4*D*t through the origin."""
    num = den = 0.0
    for t, m in zip(times, msd):
        num += t * m
        den += t * t
    return num / (4 * den) if den > 0 else 0

Run the full model from a terminal:

node simulation.js --particles 100 --diffusion 1 --seed 42
python simulation.py --particles 100 --diffusion 1 --seed 42

Download the model

Assumptions & limitations

• Overdamped, no inertia. Positions change by independent Gaussian increments; velocities do not exist in the model. Valid on time scales much longer than the velocity-relaxation time.
• Independent particles. No collisions, no hydrodynamic or electrostatic interactions, no walls: a dilute suspension in an unbounded fluid.
• Constant, isotropic diffusion. One D for both axes, everywhere, for all time.
• Exact sampling, finite resolution. Positions are exact samples of Brownian motion at multiples of Δt; the straight lines between them in the trails are only a drawing.
• Normal variates from a seeded generator. Increments use the Box–Muller transform on xoshiro128** uniforms; see pseudo-random numbers and seeds.

Related simulations

• OS-012Random walkCoin-flip steps on a line or grid: distance grows like √n, and in 1-D and 2-D every walker comes back.
• OS-037DiffusionThe heat equation solved two ways, by finite differences and by random walkers, against the exact Gaussian.
• OS-017Geometric Brownian motionThe standard stock-price model: log-normal paths whose mean and typical outcome drift apart.
• OS-010Galton boardBalls bouncing left or right through rows of pegs pile up into a binomial bell curve.

Background: Pseudo-Random Numbers and Seeds in Simulation

References

1. Einstein, A. (1905). Über die von der molekularkinetischen Theorie der Wärme geforderte Bewegung von in ruhenden Flüssigkeiten suspendierten Teilchen. Annalen der Physik, 322(8), 549–560. doi.org/10.1002/andp.19053220806
2. Uhlenbeck, G. E. & Ornstein, L. S. (1930). On the theory of the Brownian motion. Physical Review, 36(5), 823–841. doi.org/10.1103/PhysRev.36.823
3. Crank, J. (1975). The Mathematics of Diffusion, 2nd ed. Oxford University Press.