Random Surfaces from Stacking Cubes:
A Visual Journey

Leonid Petrov

University of Virginia

Loading...
NSF Simons Foundation

How Does Nature Make Shape?

Salt — same substance with cubic crystal structure, at different scales and growth conditions

Salt crystals microscopic
Microscopic · Public domain
Hopper salt crystal
Mesoscopic · Doronenko, CC BY 3.0
Halite crystal
Mesoscopic · Lavinsky, CC BY-SA 3.0
Smooth salt formation
Macroscopic · xta11, CC BY-SA 4.0

Mathematical (random) crystal growth
phantom

How can something random have a predictable shape?

Part I

Warm-up in Flatland

Paths on a Grid

Paths may be thought of as tracing a slice of the crystal by a plane parallel to the xy-plane.

Start at (0, 0), end at (a, b)
Move only right (R) or up (U)
How many paths are there?
\(\displaystyle\binom{a+b}{a} = \frac{(a+b)!}{a!\, b!}\) paths

A Random Path

Random path in 10 × 6 grid
As the mesh goes to zero, the random path gets close to a straight line!
Coin flipping: At each step independently, flip a biased coin: R with prob. \(\frac{a}{a+b}\), U with prob. \(\frac{b}{a+b}\). This describes the random boxed path at the local lattice scale.

Fluctuations Around the Diagonal

Central Limit Theorem:
The random path has fluctuations of size \(O(\sqrt{N})\) (where \(N=a+b\) is the number of steps).
At the midpoint:
\(\displaystyle \frac{Y_{\text{mid}} - \mathbb{E}[Y_{\text{mid}}]}{\sqrt{N}} \to \mathcal{N}(0, \sigma^2)\)
The histogram of midpoint heights converges to a Gaussian bell curve.
Brownian bridge:
The full fluctuation process (path minus diagonal) converges to the Brownian bridge — the Brownian motion pinned at both ends.
Universality:
The Brownian motion and Brownian bridge appear as scaling limits in most of probability theory, including random walks, random permutations, queuing theory, Kolmogorov–Smirnov statistics, random matrices, random graphs, financial models, and many other settings.
Midpoint height samples: 0


Path minus diagonal (210 × 120)

Weighted Paths

Area above a path:
Count unit squares above the path
Area = 8
What if some paths are more likely?
Weight by area under the path:
\(\displaystyle\text{Probability} = \frac{q^{\text{area}}}{\sum_{\text{paths}} q^{\text{area}}}\)
Effect of q:
  • \(q = 1\): all paths equally likely (uniform)
  • \(q < 1\): paths with less area are favored
  • \(q > 1\): paths with more area are favored
q-binomial coefficient:
\(\displaystyle \sum_{\text{paths}} q^{\text{area}} = \binom{a+b}{a}_q = \frac{(1-q) (1-q^2)\cdots(1-q^{a+b-1})}{(1-q)\cdots(1-q^a)(1-q)\cdots (1-q^b)}\)
1.00
q = 1: ordinary binomial (35 paths)

Curved Limit Shapes

Limit shape [Vershik 1996]
As \(N \to \infty\) with \(q = e^{-\gamma/N}\), the random path concentrates around a deterministic curve
\(\displaystyle A e^{-\gamma y} + B e^{-\gamma x} = 1\), where
\(\displaystyle A = \tfrac{1-e^{-\gamma}}{1-e^{-\gamma(1+\alpha)}}\), \(\displaystyle B = \tfrac{1-e^{-\gamma \alpha}}{1-e^{-\gamma(1+\alpha)}}\), \(\alpha\) = aspect ratio
Entropy vs gravity:
  • Entropy wants the straight diagonal (straight line has the most discrete paths approximating it)
  • Gravity/energy: pulls the path up for \(q<1\) (less area above is preferred)
  • The compromise is a curve: the limit shape
Variational principle:
The limit shape \(y=f(x)\) maximizes
integral entropy − \(\gamma\) · area, where \(\gamma\) is the normalized area bias, \(q=e^{-\gamma/N}\to1\).
\(\displaystyle \int \sigma(f'(x))\,dx - \gamma \int f(x)\,dx\),
\(\qquad \qquad \qquad \sigma(s)\coloneqq (s+1)\log(s+1) - s\log s\)
3.6 (q = 0.990)
a = 240, b = 180 · Loading...
Fluctuations are still Brownian bridge: Random curve minus limit shape

Summary of Flatland Exploration

Local rules — paths go up or right
Boundary conditions — paths live in a box
Limit shape — from entropy/energy balance;
maximizes \(\displaystyle\int \sigma(f') \, dx - \gamma \int f \, dx\)
Local universality — coin flips with \(p\) = slope of limit shape
Universal global Gaussian fluctuations — Brownian bridge (time-changed) for our 2D situation
a=240, b=180, q=0.99 · Loading...
Local zoom — IID Bernoulli
Path − limit shape ≈ time-changed BB
3.6 (q = 0.990)
a = 180, b = 135 · Loading...
Zoomed in: local segment
Locally, paths don't know about the global weight — still IID Bernoulli!
Why? The limit shape maximizes entropy minus γ × area. This yields IID Bernoulli with position-dependent \(p\) = slope of the limit shape. The weight only affects the global shape, not the local behavior.
-->

Part II

Random Surfaces

From Paths to 3D

surface in 1 × 12 × 9 box = 2D path in 12 × 9 rectangle
Theorem [MacMahon, c. 1900]:
Number of surfaces in \(a \times b \times c\) box is \(\displaystyle\prod_{i=1}^{a}\prod_{j=1}^{b}\prod_{k=1}^{c}\frac{i+j+k-1}{i+j+k-2}\)
For 9 × 12 × 9: 1,340,992,301,315,806,672,824,460,528,012,500
(~1034, about 1000× more than virus particles on Earth)
... and pick one at random from this collection (how? — end of the talk)
More complex boxes:

Limit Shape in 3D

pre-sampled
Concentration: [Cohn–Kenyon–Propp 2000]
As the mesh goes to zero, the random surface concentrates around a deterministic limit shape surface.
Entropy vs boundary conditions leads to curving — even for uniformly random tilings, the limit shape is not flat!
Moreover, it contains flat facets which are parallel to the coordinate planes.

Local Patches

Local patches:
  • Approximate the surface by piecewise-planar patches
  • Each patch is a flat random surface with a given slope \((s_1, s_2)\)
  • For a given slope, the number of approximating 3D surfaces is \(\approx e^{N^2\cdot \sigma(s_1,s_2) }\)
  • Not IID Bernoulli! Unlike 1D coin flips, the 2D patches are determinantal point processes — joint occupation probabilities have determinantal form
Variational principle [Cohn–Kenyon–Propp 2000]:
Limit shape \(z=h(x,y)\) maximizes \[\displaystyle\iint \sigma(\nabla h)\, dxdy\] subject to boundary conditions on \(h(x,y)\)
\(\sigma(s_1,s_2) = \tfrac{1}{\pi}\bigl[\Lambda(\pi s_1) + \Lambda(\pi s_2) + \Lambda(\pi (1\!-\!s_1\!-\!s_2))\bigr]\)
\(\Lambda(\theta)=-\!\int_0^\theta \ln|2\sin t|\,dt\)   (Lobachevsky function)
Other 3D building blocks besides squares of three directions: [Kenyon–Okounkov–Sheffield 2006]
plane z=Ax+By

Lozenge Tilings

3 types of lozenges:
Rhombi of 3 orientations, projections of cube faces onto the \(x+y+z=0\) plane.

Frozen vs Liquid Regions

Hexagon
Sawtooth
General polygon

Height Fluctuations in 3D

Gaussian Free Field
Like the Brownian bridge, but with two-dimensional time. Take two independent random tilings \(h_1, h_2\). Each one:
\(\displaystyle h_1 = \underbrace{\bar{h}}_{\text{limit shape}} + \underbrace{\text{fluctuations}_1}_{??\text{ size}}\)
\(\displaystyle h_2 = \underbrace{\bar{h}}_{\text{limit shape}} + \underbrace{\text{fluctuations}_2}_{??\text{ size}}\)
Subtract: limit shapes cancel!
\(\displaystyle {h_1 - h_2} = {\text{fluct}_1 - \text{fluct}_2} \approx \text{Gaussian}\)
How big are the fluctuations?
  • In 2D (paths): \(\sqrt{N}\)
  • In 3D (surfaces): only \(\sqrt{\log N}\) — surprisingly flat!
GFF fluctuations
Governed by the Gaussian Free Field (GFF) — the analog of a Brownian bridge with 2D "time". Fluctuations grow to infinity at each point, but slowly — we can make sense of a random distribution (generalized function) in the sense of Schwartz.
Universal: same GFF for all simply-connected domains. [Kenyon 2001], [Petrov 2014], [Bufetov–Gorin 2018]

What Else Can We Tile by Lozenges?

Non-simply-connected
Non-planar orientable surface
Non-planar orientable
Non-orientable surface
Non-orientable (Möbius strip)
GFF on Escher-like domain
State of the art [Borot–Gorin–Guionnet 2026]: Universal Gaussian fluctuations for some non-simply-connected, non-planar, or non-orientable domains.
Fluctuations for generic simply-connected domains: still open!

q-volume tilt

q =
q =
q =
q =
q =
\(\displaystyle \text{Probability} = \frac{q^{\text{volume}}}{Z}\)
Same idea as 2D: weight each surface by \(q^{\text{volume}}\)
Entropy vs Energy:
  • \(q < 1\): less volume preferred (gravity pulls down)
  • \(q > 1\): more volume preferred (buoyancy pushes up)
New phenomena in weighted models [Knizel–P. 2025]:
In an interpolation between \(q^{\text{volume}}\) and \(q^{-\text{volume}}\), keeping \(q\) fixed, one sees a collapse of 2D random behavior into a new 1D random interface, the waterfall.

Waterfall eureka: Formulas teach us [Knizel–P. 2025]

Dimensional collapse:
Horizontal lozenges vanish → 2D surface becomes a 1D random sequence of two lozenge types
What does the 1D sequence look like?
Put the physics hat on: write a formal series for the density of lozenges:
\(\displaystyle \rho(t) = \sum_{n \in \frac{1}{2}\mathbb{Z}_{\geq 0}} \frac{(-q^{-n})\mathcal{F}_n(-\frac{t}{2})\mathcal{F}_n(1-\frac{t}{2})}{\|\mathcal{F}_n\|^2_{\ell^2(\mathbb{Z})}}\)
But... this series does not converge!
Diverges "mildly": \(\quad +c\;-c\;+c\;-c\;+c\;-\;\cdots\)
The divergence IS the answer!
\(\rho_0\) := lim (even partial sums) \(\approx 0.5611\)
\(\rho_1\) := lim (odd partial sums) \(\approx 0.4389\)
\(\rho_0+\rho_1 = 1\), both independent of \(t\)
The 1D sequence is two-periodic (!)

Part III

How to Make Pictures

Glauber Dynamics

Elementary move:
Three lozenges around a hexagon can be flipped — rotated 180°. This is the only local move that preserves the tiling.
Markov chain:
  • Pick a random site with a \(1\times 1\times 1\) hexagon on it
  • If it's flippable, flip it
  • Repeat many times
After enough steps, the tiling is approximately uniformly random (close to equilibrium).
How fast?
For the hexagon, Glauber dynamics mixes in \(O(N^{4+\epsilon})\) rotations [Laslier–Toninelli 2015], [Caputo–Martinelli–Toninelli 2012]. Mixing means that the two distributions become statistically indistinguishable.

Exact Sampling: We need to know when to stop waiting

Coupling From The Past (CFTP) [Propp–Wilson 1996]
Surfaces are partially ordered by height function. The Markov chain is monotone, so it suffices to track just the maximal and minimal chains.
  • The stationary chain is somewhere in between the maximal and minimal surfaces
  • Run coupled dynamics — when max and min coalesce, we've caught the stationary chain
  • To ensure exactness, use the backward-doubling protocol: start chains at times \(-1, -2, -4, -8, \ldots\) and check coupling at time \(0\)
  • Output is an exact sample from the stationary distribution
Not approximately random — EXACTLY random.
Compare with calculus: we know \(1/n \to 0\), but for which \(n\) is it equal to zero?
Here we reach the limit exactly — but at a random time.
Hexagon 35 × 35 × 35
QR code to lozenge-draw lpetrov.cc/lozenge-draw (opens in new tab)

When Sampling is Hard

CFTP can be slow:
Exponentially slow in \(N\) — prohibitively long to get an exact sample: the monotone sandwich may take \(\sim e^{cN}\) steps to collapse.
This happens in the waterfall model [Knizel–P. 2025] and also in multi-\(q\)-weighted models.
Why?
The CFTP needs to go through local moves, and some local configurations are exponentially unlikely. For example, in the waterfall, we do not have horizontal pieces in the middle — but to mix, horizontal pieces need to cross this gap.
Alternatives come from algebra in specific cases:
  • Commuting Markov matrices for waterfalls
    [Borodin–Gorin–Rains 2009]
  • Powerful local equivalence relations that transform the lattice and weights in algebraically tractable manner (Yang–Baxter equation / urban renewal / spider move)
Hexagon 35 × 35 × 35

What Algebraic Structures Can Do for You

Commuting transfer matrices [Borodin-Gorin 2008], [Borodin–Gorin–Rains 2009]
The algebraic structure of the weights (highly specific to the hexagon) gives commuting Markov operators that map the desired distribution to the desired distribution on a different hexagon.
Growing the hexagon:
Start from a parallelogram — it has exactly one tiling, so it is already an exact sample.
Then apply the Markov operators one by one, changing the hexagon step by step (increasing one side, decreasing the other).
Each step preserves exactness!
Polynomial time:
Each step costs \(\sim N^2\) (area), and we need \(\sim N\) steps — exact sample in \(O(N^3)\) time.
Compare: Glauber CFTP needed \(\sim e^{cN}\) steps on the previous slide.
Step: 0 / 80

Summary

Limit shapes
Random paths and surfaces concentrate around deterministic shapes as mesh shrinks (\(=\) size grows).
Universality
  • Local statistics everywhere share the same nature, independent of boundary shape.
  • Gaussian fluctuations on the global scale.
Exact sampling
Almost every picture in this talk was an exact random sample, computed live as we went along.
How can something random have a predictable shape?
— Through a careful balance of entropy and energy: Local rules create global structures.
QR code for lozenge-draw
draw any border and lozenge tile! (plus 3D printing exports)
lpetrov.cc/lozenge-draw (opens in new tab)