1 of 40

Stochastic Dynamics

Undergraduate Students: Giles El-Assal, Mateo Price-Otero, Rachel Rakushkin

Post Doctorate: Debangana Mukhopadhyay

Advisor: Ken Duffy

1

2 of 40

Presentation Overview

Introduction
Diffusion Limited Aggregation
Monte Carlo Simulations

Metropolis Monte Carlo
Kinetic Monte Carlo

Comparison and Takeaways

2

3 of 40

Introduction

Introduction
Diffusion Limited Aggregation
Monte Carlo Simulations
Comparison and Takeaways

3

4 of 40

Project Goal

To study, model, and analyze generative stochastic processes from which patterns emerge

4

To study, model, and analyze generative stochastic processes from which fractal-like patterns emerge

---------------------------------------------------------------------------------------------------------

Our goal is to study stochastic processes that generate fractal-like structure, and to organize them along one key axis: memory. A memory-dependent process carries a permanent record of its past — each step depends on everything that came before. DLA is our example: once a particle sticks it is frozen, and it shapes where every future particle can land. A memoryless, or Markov, process forgets — the next state depends only on the current state, not the path taken to reach it. Metropolis and Kinetic Monte Carlo are our examples there. The two panels show these generative processes building structure from simple local rules. Everything in the talk is a comparison between these two families.

5 of 40

Motivation

5

Despite much work being done, certain features remain unexplained (no clear analytical explanations, no rigorous proofs)

The same rate-limited growth recurs across unrelated systems

Biological Behaviors

Bacteria colony growth
Antibody to antigen attachments

Dielectric Breakdown

Dielectrics of semiconductor devices
Lightning

Formation in Nature

Snowflakes

Fabrication Techniques

Physical Vapor Deposition

Physical Vapor Deposition

Snowflakes

Dielectric Breakdowns

Bacteria Colony

Why study this? First, despite decades of work, key features still lack rigorous explanation — there is no first-principles analytical derivation of the fractal dimension, and no proof that these clusters are truly scale-invariant.

Second, and more striking, the same rate-limited growth recurs across completely unrelated systems. That is not a coincidence: in each case growth is controlled by a field obeying Laplace's equation at the boundary — concentration for bacteria and for deposition, electric potential for dielectric breakdown and lightning. Same governing equation, same emergent geometry. So we see it in bacterial colonies and blood-cell attachment, in dielectric breakdown and lightning, in snowflake growth, and in thin-film deposition for nanofabrication. One simple stochastic rule reproduces the morphology of all of them, and that universality is what makes the model worth understanding.

6 of 40

Diffusion Limited Aggregation

Introduction
Diffusion Limited Aggregation
Markov Chain Monte Carlo
Comparison and Takeaways

6

7 of 40

Diffusion Limited Aggregation

A Repetitive Game

Initial Conditions: seed at some position and a defined boundary

Randomly select position from boundary
Begin a random walk from that position
When walker reaches seed/cluster it may stick and become part of the cluster; if it does, repeat steps 1-3 with new walker

Increasing Complexity

Vary the sticking probability
Bias the walk
Vary spawn positions
Increase number of particles
Allow for cluster movement

7

8 of 40

8

d=2

Modeling Diffusion

Brownian Motion

d=1

Modeling DLA begins with modeling the diffusion process which we naturally think of as a continuous process. In fact there’s this famous diffusion equation which can be used to solve for the probability function for a given particle. This function is actually a solution to the continuous limit of a discrete random walk. The catch is that this Gaussian is the free-space solution: it assumes no boundaries. In DLA we have an absorbing aggregate, so that assumption breaks immediately — and worse, the boundary evolves every time a particle sticks. So instead of the diffusion equation with its clean Gaussian solution, we end up needing to solve a slightly different equation (Laplace's equation ∇^2 𝑢=0) on a constantly changing and irregular shaped unpredictable domain, and this has no closed form.

Therefore for our simulation work we use the discrete version of the random walk. In the simplest case we can use this equation as our model, and what this equation is telling us is that the probability u of a random walker being on the green dot in the next time step is simply the sum of the probabilities of it being on each of the two possible adjacent positions at the previous time step. We can generalize this to 2 or more dimensions.

9 of 40

9

DLA Aggregate growth with 10,000 particles and P = 1

DLA Aggregate growth with 10,000 particles and P = 0.1

Example Implementation and Analysis

Density-Density Correlation

r

C(r)

10 of 40

Cluster Diffusion 1/N

10

This is the most complex variant of DLA we looked at

Instead of single particles sticking to one growing seed, we let many clusters (initially particles) move and merge.

Each cluster diffuses as a whole, but with the mobility scaling as one over its size — bigger clusters move more slowly, exactly as real diffusion predicts. When two clusters collide they fuse into one. The left panel is a snapshot of all the clusters in the box; the right is the live distribution of cluster sizes.

Early on, as shown here, there are hundreds of small clusters — mean size near one, the largest only seven. As the run proceeds, small clusters coalesce into fewer, larger ones and the distribution broadens. This is cluster–cluster aggregation, a better model for systems whose building blocks are themselves mobile.

---------------------------------------------------------------------------------------------------------

We start with N many particles, we randomly pick a particle and have it perform a step as part of a random walk, we then decide whether to accept or not accept this move. If we don’t accept the move it return to where it was before, and we retry. We continue sampling and accepting or rejecting moves until the particles come together into one cluster. The catch here however, is that we have the probability of accepting a move be dependent on size; here for example we used a 1/N rule, thus clusters with more particles had lower probability of moving.

11 of 40

Monte Carlo Simulations

Introduction
Diffusion Limited Aggregation
Monte Carlo Simulations
Comparison and Takeaways

11

12 of 40

Context

12

Now we shift to the Monte Carlo half.

Suppose we want to quantify some observable f of a system

We can do this by performing an average over all possible states where we weight each state by its probability p — the sum, or integral, at the top.

The problem with this method is it becomes intractable for cases with many particles

Take one mole of gas. Each particle needs six coordinates, three position and three momentum, and if we discretize each axis into m values, a single particle already has m-to-the-sixth states — and the whole mole has that raised to the power of Avogadro's number. We cannot enumerate that; gridding the space is hopeless.

The alternative is sampling: rather than summing over every state, we draw representative states from the distribution p and average f over just those samples. By the law of large numbers, that converges to the true average.

So everything reduces to one question which is what is our distribution?

---------------------------------------------------------------------------------------------------------

Shifting gears a bit to consider how we can analyze more complex systems so in general terms, suppose we would like to model and analyze a system as it evolves as well as when it is in its final state meaning we want to quantify certain observable features of the system

One way to do this would be to estimate the average value of whatever metric we are interested in by just multiplying the probability of seeing a state by the specific value the metric of our interest takes on for that state and do this for all states (do a weighted average of values)

However this becomes computationally intractable when there are many states. For example, a simple diatomic gas, each particle can be moving in x y or z direction and has momentum in each as well thus each molecule has 6 different states that need to be tracked. Consider a mol of gas that means there are a total of ___ 6^10^23 many different states we must account for numerically or you can think of x as a 6^10^23 dimensional vector

Alternatively we can draw samples from a distribution, note the value of our metric of interest, and then estimate the average value. Law or large numbers means that as we do this for more samples our estimate will converge to the actual true average value.

So it seems all we need to know is our distribution function

13 of 40

13

1

Probability proportional to multiplicity

Probability of system being in microstate i, is proportional to the number of reservoir microstates holding the rest of the energy

2

Taylor expansion

Taylor-expand about total energy

Substitute

Exponentiate

Re-exponentiating gives the reservoir count as a Boltzmann factor

Example: Deriving a Distribution

1.

2.

3.

4.

14 of 40

Markov Chain Monte Carlo

14

Okay great so if we know the distribution or we can derive it from first principles we can easily obtain samples and perform an estimation to get whatever desired metric we cared to obtain :

However we don’t always know our distribution. This is a good case for the Markov Chain Monte Carlo process

A Markov chain is a memoryless sequence of states, each drawn from transition probabilities T.

Markov Chains of a particular sort have a useful feature. Specifically, if p is stationary under T (running one step of the chain leaves p unchanged) and if the chain is ergodic, meaning irreducible and aperiodic, then the chain converges to p from any starting point.

To guarantee stationarity we impose a stronger, sufficient condition called detailed balance, the balanced-flux equation shown here. The Metropolis algorithm is an example of this which we will get into in the next slides

---------------------------------------------------------------------------------------------------------

Irreducible = for any two states there is a non zero probability for the transition to occur over some number of steps

Aperiodic = never gets stuck toggling between states in a periodic fashion forever

Stationary Probability of the next state is a particular state x’ is the same as the probability that that is the current state

Irreducible means for any two states in your set of states, there exists a non zero probability that you transition between them

15 of 40

Metropolis Monte Carlo: A Lattice Gas and Its Energy

The system

• 2D lattice; each site is empty (σ = 0) or holds a particle (σ ∈ {1,…,k})

• Our runs: 150×150 lattice, density 0.15 (3375 particles), k = 1

• Conserved (Kawasaki) dynamics: a particle hops to an adjacent empty site — particle number is fixed

The energy (Hamiltonian)

H = −J Σ δ(σ_r , σ_r′) · [ 1 − δ(σ_r , 0) ]

• Sum runs over nearest-neighbor pairs in 4 lattice directions

• A pair counts only if same-state AND occupied

• Each such “bond” lowers the energy by J

• Lower energy ⇔ more same-type contacts ⇔ larger, denser clusters

(units: J = 1, k_B = 1, so T is in units of J / k_B)

Temperature & the point

• Acceptance Critera: P(config) ∝ e^(−ΔH / k_B * T)

• Low T → energy wins → condensation / clustering

• High T → entropy wins → disordered “gas”

• MMC samples this distribution to measure the equilibrium structures (cluster size, fractal dimension) that emerge vs. T

15

16 of 40

Metropolis Monte Carlo (MMC)

Propose move

Randomly select a particle; propose a new adjacent lattice site.

Compute ΔH

Calculate energy difference: ΔH = H(new) − H(old).

Accept / Reject

If ΔH ≤ 0: accept. If ΔH > 0: accept with prob. exp(−ΔH/T).

Update & record

Apply move if accepted; increment step counter; record observables.

Equilibration

Discard early sweeps (transient); measure only after equilibration.

Sweep = N moves

One MC sweep ≈ N attempted moves (one per particle on average).

What MMC does

- Samples equilibrium distribution correctly��- Efficient configuration-space exploration��- Clean energy / temperature control��

16

1.

2.

3.

4.

5.

6.

17 of 40

Metropolis Monte Carlo – Single Particle Movement

17

18 of 40

Review – Radius of Gyration and Fractal Dimension

Radius of gyration: the root-mean-square distance of a cluster's particles from its own center of mass

18

Dimension is the exponent linking size to mass: Area = Radius/Side-Length^(Dimension)

Example: Take a square of side length 1 and double the length. Area increases from 1 to 4. Take a cube of side length 1 and double the length. Area increases from 1 to 8.

Size ∝ R_g^(D_f) with 1 < D_f < 2 in 2D. The fractional exponent D_f measures how efficiently the cluster fills space — 2 = compact blob, 1 = thin chain.

19 of 40

Metropolis Monte Carlo – Single Particle Movement

19

Panel 1 — Equilibration: Total Energy

System Energy vs. sweeps, averaged over 3 trials
Lower temperature leads to lower energy as same-state bonds form
Hot runs (T = 1.0, 2.5) barely relax due to constant noise breaking bonds

Panel 2 — Phase State: Cluster Count

Number of distinct clusters ⟨N_c⟩ vs. sweeps
Cooling drives coalescence into fewer, larger domains
Structural mirror of the energy drop in Panel 1

Panel 3 — Morphology: Fractal Dimension

R_g vs. cluster size (log–log); slope gives fractal dimension: 1 / D_f
Most compact at T = 0.5 (D_f ≈ 1.70) — same run that hit lowest energy
T = 0.1 most ramified (D_f ≈ 1.30) → frozen before it can compact

20 of 40

Metropolis Monte Carlo – Single Particle Movement

20

21 of 40

Metropolis Monte Carlo – Single Particle Movement

21

T=0.5 ; Sweeps = 500 ; Particles = 3375 ; 150 x 150 lattice

run_simulation 7.73 s — the whole run; everything nests inside it
perform_sweep 6.27 s — the gap to the right (~1.5 s) is measurement: per-sweep measure_clusters_advanced (full-lattice DFS), R_g snapshots
attempt_single_move 5.16 s — the gap to the right (~1.1 s) is per-sweep RNG and looping
Space below attempt_single_move — per-attempt Python logic (boundary/exclusion checks, ΔH arithmetic)
Only ~⅔ of runtime the simulation physics. ~20% is measurement, ~14% is RNG/loop scaffolding

22 of 40

Why Add Cluster Moves?

Single-particle moves

• At low T, any move that breaks a bond has ΔH > 0 and is almost always rejected (e^(−ΔH/T) ≈ 0)

• A compact cluster can only rearrange at its surface — whole domains can’t translate or merge

• System freezes in a metastable, partly-coarsened state; never reaches equilibrium in feasible sweeps

Rigid cluster moves

• Propose translating an entire connected cluster by one lattice step

• Restores large-scale mobility, allowing domains to diffuse and coalesce toward true equilibrium

• Size-dependent acceptance (mobility ∝ 1/√size): big clusters move less often

21

Identifying and moving clusters is expensive, and the cost grows as clusters coarsen

23 of 40

Metropolis Monte Carlo – Single+Cluster Particle Movement

23

24 of 40

Metropolis Monte Carlo – Single+Cluster Particle Movement

24

Panel 1 — Equilibration: Total Energy

• System energy vs. sweeps, single+cluster, averaged over [3] trials

• Cluster moves merge whole domains, meaning energy drops faster and reaches lower values than single moves alone

• Cold runs no longer freeze: T = 0.1 now relaxes toward equilibrium

Panel 2 — Phase State: Cluster Count

• Distinct clusters ⟨N_c⟩ vs. sweeps

• Domain diffusion speeds coalescence → fewer, larger clusters than single-mode

• At low T the system can condense to [a single dominant cluster]

Panel 3 — Morphology: Fractal Dimension

• R_g vs. cluster size (log–log); slope gives 1 / D_f

• All temperatures converge to D_f ≈ 1.8 — close to compact 2D aggregates

• Reaches a more compact equilibrium than single moves at the same T

25 of 40

Metropolis Monte Carlo – Single+Cluster Particle Movement

25

Runtime: the cost of cluster dynamics

• Cumulative runtime vs. sweep count, single+cluster, density 0.15

• Single-particle cost is ≈ linear; single+cluster grows super-linearly at low T. This is because as domains coarsen, each move must traverse a larger cluster

26 of 40

Metropolis Monte Carlo – Single+Cluster Particle Movement

26

• Cluster identification (DFS) dominates the cost

• The self-runtime in attempt_cluster_move covers the visited-array allocation, the boundary/collision checks over every cluster cell, the energy (ΔH) recalculation, and the coordinate bookkeeping that rewrites particle positions on an accepted move.

• Optimizations: kinetic pre-check from a once-per-sweep cluster-size snapshot (scipy labeling) + O(cluster-size) coordinate updates

• Previous 500-sweep runs were taking hours. Now only takes 195 s allowing full data collection

T=0.5 ; Sweeps = 500 ; Particles = 3375 ; 150 x 150 Lattice

27 of 40

Metropolis Monte Carlo – Single vs. Single+Cluster

27

Energy (left): single-particle freezes at ⟨H⟩ ≈ −3,000; single+cluster keeps relaxing to ≈ −4,700 — the true energy floor
Cluster count (right): single stays fragmented at ~600 clusters; single+cluster coalesces to a handful

28 of 40

Why we need another model

Metropolis Monte Carlo has no meaningful time

No rates, no diffusion coefficient, no timescales

Code stalls when temperature is low

28

29 of 40

Kinetic Monte Carlo (KMC)

Program Steps

Enumerate all possible events from the current state, each with rate rₖ
Compute the total rate R
Select event k with vector entry rₖ / R (cumulative sum + uniform draw)
Execute event k and update the configuration
Advance time stochastically by Δt

29

30 of 40

Example Implementation

30

31 of 40

KMC Cluster Analysis

31

32 of 40

Comparison and Takeaways

Introduction
Diffusion Limited Aggregation
Monte Carlo Simulations
Comparison and Takeaways

32

33 of 40

KMC vs MMC Comparison (Cluster Size)

33

34 of 40

KMC vs MMC Comparison (Cluster Number)

34

35 of 40

KMC vs MMC (Final-State Comparison)

35

36 of 40

Takeaways

Each model has its use

DLA for diffusion-limited growth
MMC for correct equilibrium structure
KMC for true kinetics

Equilibrium vs Non-equilibrium

36

Model	Purpose	“Time” is…
DLA	a specific non-equilibrium growth process	Physical
MMC / MCMC	sample an equilibrium distribution	Fictitious (MC steps)
KMC	simulate real-time evolution	Physical

37 of 40

Future Work

Higher dimensions and off-lattice movement
Longer KMC runs
Runtime reduction
Applying these models to real-world experimental results

37

38 of 40

Sources

Tronnolone, Hayden, et al. "Diffusion-limited growth of microbial colonies." Scientific Reports 8.1 (2018): 5992.

Witten, Thomas A., and Leonard M. Sander. "Diffusion-limited aggregation, a kinetic critical phenomenon." Physical review letters 47.19 (1981): 1400.

Witten, Thomas A., and Leonard M. Sander. "Diffusion-limited aggregation." Physical review B 27.9 (1983): 5686.

Newman, Mark. Computational Physics. Revised ed., CreateSpace Independent Publishing Platform, 2013.

38

39 of 40

Thank you

Q & A

39

40 of 40

40