Random structures: walks, percolation and random geometry

In the early 1920s George Pólya, walking in the woods above Zürich, kept running into the same young couple, far more often than seemed polite. He wondered whether two random walkers on a grid of paths are bound to meet again and again — equivalently, whether a single random walker is bound to return to where it started. His answer (1921) is one of the most quoted theorems in probability: on the line and in the plane a simple random walk returns home with probability one; in three or more dimensions it does not. A remark often attributed to Shizuo Kakutani puts it memorably: a drunk man will find his way home, but a drunk bird may get lost forever.

Pólya's theorem is the entrance to a large territory where randomness creates shape. Let each pipe in a porous rock be open or blocked at random (Broadbent and Hammersley's model of a gas-mask filter, 1957): does water seep through? Join $n$ people by random acquaintances (Erdős and Rényi, 1959–60): is there a giant connected community? Let a random walk take smaller and faster steps: what does the limiting curve look like, and what is the shape of its outer boundary? All these systems have phase transitions — sharp thresholds where global behaviour changes abruptly — and at the threshold, beautiful random fractals appear. This lesson is an orientation: the classical theorems with full proofs, the deep theorems stated precisely, and the open problems marked as open.

Returns. A walk returns with probability one exactly when the expected number of returns is infinite. After $2n$ steps the chance of being at the origin of $\Z^d$ is about $c_d\, n^{-d/2}$ (by the central limit theorem, the walk has spread over a region of volume about $n^{d/2}$). The series $\sum n^{-d/2}$ diverges for $d \le 2$ and converges for $d \ge 3$. That is the whole story, once each step is made rigorous.

Currents. Make every edge of the grid a $1\,\Omega$ resistor and pump current from the origin to "infinity". If the total resistance to infinity is infinite, no current escapes, and the random walk — whose hitting probabilities are exactly the voltages — keeps coming back. In the plane the number of resistors in parallel at distance $k$ grows only like $k$, so the resistances of successive shells add up like $\sum 1/k = \infty$; in space they grow like $k^2$ and $\sum 1/k^2 \lt \infty$.

Clusters. In percolation, open edges form clusters. For small $p$ every cluster is finite: the number of possible long paths grows like $3^n$, but each is open with probability $p^n$. For $p$ close to $1$ the closed edges are rare and cannot surround the origin. In between lies a critical value $p_c$, which for the square lattice is exactly $\frac12$ — because the lattice is its own planar dual.

We assume stochastic processes (Markov chains, the strong Markov property, Brownian motion) and the probabilistic method; see also Markov chains.

The proof for $d = 1, 2$ is formalised in Lean below, apart from the probabilistic modelling (the formula for $u_{2n}$ and Theorem 1). For $d = 3$ the Stirling step is cited, not proved. The first experiment computes $G(N) = \sum_{n \le N} u_{2n}$ for all three dimensions: it grows like $2\sqrt{N/\pi}$, like $\frac{1}{\pi}\log N$, and converges to Watson's constant $1.5164$, so the walk on $\Z^3$ returns with probability $1 - 1/1.5164 = 0.3405$.

Recurrent is not the same as "returns quickly". For the walk on $\Z$, the first return time $T$ satisfies $P(T \gt 2n) = u_{2n}$ (a classical reflection-principle identity; see Feller, Vol. I, Ch. III). So $E[T] = 2\sum_n u_{2n} = \infty$: the walk returns with probability one, but the expected time it takes is infinite. In Lean, return_prob_tendsto_zero and return_1d_not_summable are exactly the two halves of this picture.

Let $\Gamma$ have unit conductances, so the walk is simple random walk.

Two tools make the criterion usable, both proved in Doyle–Snell and Lyons–Peres (Ch. 2):

• Rayleigh's monotonicity principle. Increasing any edge resistance does not decrease any effective resistance. So a subgraph of a recurrent graph is recurrent, and a graph containing a transient subgraph is transient — another proof that $\Z^d$, $d \ge 4$, is transient once $\Z^3$ is.
• The Nash-Williams inequality (1959). If $\Pi_1, \Pi_2, \dots$ are disjoint sets of edges, each separating $a$ from infinity, then $R_{\mathrm{eff}}(a \leftrightarrow \infty) \ge \sum_k |\Pi_k|^{-1}$.

$\Z^2$ again. Let $\Pi_k$ be the edges joining $\partial[-k, k]^2$ to $\partial[-(k+1), k+1]^2$: there are $8k + 4$ of them, so $R_{\mathrm{eff}} \ge \sum_k \frac{1}{8k + 4} = \infty$ and the walk is recurrent. (The divergence is proved in Lean.)

$\Z^3$ again. T. Lyons (1983) proved that a network is transient if and only if there is a unit flow from $a$ to infinity with finite energy $\sum_e r(e)\, i(e)^2$. Spreading a unit flow evenly over the roughly $k^2$ edges at distance $k$ gives energy about $\sum_k k^2 (k^{-2})^2 = \sum_k k^{-2} \lt \infty$. Making this precise requires some care with the cube's geometry (Lyons–Peres, Ch. 2); the convergence of $\sum k^{-2}$ is in Mathlib.

(In fact $\sigma(n)^{1/n} \to \mu \approx 2.638$, the connective constant of $\Z^2$, whose exact value is unknown.)

Kesten's theorem (1980). For bond percolation on $\Z^2$, $p_c = \frac12$; moreover $\theta(\frac12) = 0$ (Harris, 1960). The proof combines self-duality — at $p = \frac12$ the dual configuration has the same law as the primal one, so suitably chosen rectangles are crossed with probability exactly $\frac12$ (the experiment checks this by exhaustive enumeration for small rectangles) — with the Russo–Seymour–Welsh theory (1978), which turns crossings of squares into crossings of long rectangles with probability bounded below, and a sharp-threshold argument showing that the crossing probabilities jump from near $0$ to near $1$ in a window around $\frac12$ that shrinks with the size. Modern proofs use the Kahn–Kalai–Linial influence inequality (Bollobás–Riordan, 2006) or the randomised-algorithm method of Duminil-Copin, Raoufi and Tassion (2019). For site percolation on $\Z^2$ there is no self-duality and no known closed form: $p_c^{\mathrm{site}} \approx 0.5927$, an estimate from simulation and rigorous numerical bounds.

Random graphs (Erdős–Rényi, 1960). For $G(n, c/n)$: if $c \lt 1$, w.h.p. every component has $O(\log n)$ vertices; if $c \gt 1$, w.h.p. there is a unique giant component with $(\beta + o(1))n$ vertices, where $\beta \in (0,1)$ solves $\beta = 1 - e^{-c\beta}$, and all other components have $O(\log n)$ vertices; at $c = 1$ the largest component has order $n^{2/3}$. The proof explores a component vertex by vertex; each newly found vertex has about $\operatorname{Bin}(n, c/n) \approx \operatorname{Poisson}(c)$ new neighbours, so the exploration is approximately a branching process with mean $c$, whose survival probability $\beta$ solves the equation above. For $c \le 1$ that equation has no positive root (proved in Lean), matching the absence of a giant. The giant appears at $p = 1/n$; the graph becomes connected only much later, at $p = (\log n)/n$.

Brownian motion as a limit. Donsker's invariance principle (1951): the rescaled walk $t \mapsto S_{\lfloor nt\rfloor}/\sqrt n$ converges in distribution to Brownian motion, whatever the step distribution (mean zero, finite variance). Pólya's theorem has a continuum twin: planar Brownian motion is neighbourhood-recurrent (it returns to every disc, though it hits a given point with probability zero), while Brownian motion in $\R^3$ is transient. Planar Brownian paths have Hausdorff dimension $2$ (a classical result; see Mörters–Peres). Their outer boundary — the frontier of the path — was conjectured by Mandelbrot (1982), from simulations, to have dimension $\frac43$; Lawler, Schramm and Werner proved it (2001).

Schramm–Loewner evolution. Oded Schramm (2000) observed that a random curve in a planar domain which is conformally invariant and has a "domain Markov property" must be one of a one-parameter family $\mathrm{SLE}_\kappa$, obtained by running Loewner's differential equation with driving function $\sqrt\kappa\, B_t$. Established identifications of scaling limits:

• loop-erased random walk → $\mathrm{SLE}_2$, uniform spanning tree Peano curve → $\mathrm{SLE}_8$ (Lawler–Schramm–Werner, 2004);
• critical site percolation on the triangular lattice: Smirnov proved Cardy's formula and conformal invariance (2001), the exploration path converges to $\mathrm{SLE}_6$, and Smirnov and Werner derived critical exponents such as $\theta(p) = (p - \frac12)^{5/36 + o(1)}$ as $p \downarrow \frac12$ (2001);
• interfaces of the critical Ising model → $\mathrm{SLE}_3$ (Smirnov and coauthors, around 2010).

Werner (2006), Smirnov (2010) and Duminil-Copin (2022) received Fields Medals for work in this area.

Random planar geometry. Uniformly random planar maps (graphs embedded in the sphere) with $n$ faces, with distances rescaled by $n^{-1/4}$, converge to a random metric space, the Brownian map (Le Gall; Miermont, 2013), homeomorphic to the sphere and of Hausdorff dimension $4$. Miller and Sheffield (2015–16) showed that the Brownian map is equivalent to Liouville quantum gravity with parameter $\sqrt{8/3}$, a random geometry built from the Gaussian free field.

Open problems (as of the references; check recent surveys):

• Is $\theta(p_c) = 0$ for bond percolation on $\Z^3$? Proved in dimension $2$ and in high dimensions (Hara–Slade, 1990, $d \ge 19$; Fitzner–van der Hofstad, 2017, $d \ge 11$), open in $3 \le d \le 10$.
• Conformal invariance of critical bond percolation on $\Z^2$. Smirnov's proof uses the special symmetry of the triangular lattice; universality, though widely believed, is open for $\Z^2$.
• The self-avoiding walk. Its scaling limit in the plane is conjectured to be $\mathrm{SLE}_{8/3}$ and its connective constant on $\Z^2$ is unknown; on the hexagonal lattice the connective constant is $\sqrt{2 + \sqrt 2}$ (Duminil-Copin and Smirnov, 2012).

Mathlib has measure-theoretic probability but no random walks on lattices, percolation or random graphs; and the probabilistic modelling (formulas for $u_{2n}$, Theorem 1, duality) would be a substantial project. What is formalised is the exact mathematics that the proofs above reduce to — the parts where the actual estimates live.

The bounds on central binomial coefficients, by induction, and the resulting divergence for $\Z$ and $\Z^2$:

/-- `4ⁿ ≤ (2n + 1) · C(2n, n)`: the central binomial coefficient is the largest of the `2n + 1`
coefficients `C(2n, k)`, which sum to `4ⁿ`. This applies Mathlib's
`Nat.four_pow_le_two_mul_add_one_mul_centralBinom`. -/
theorem four_pow_le_central_binom (n : ℕ) : 4 ^ n ≤ (2 * n + 1) * centralBinom n :=
  Nat.four_pow_le_two_mul_add_one_mul_centralBinom n

/-- **Lower bound** `16ⁿ ≤ 4n · C(2n, n)²` for `n ≥ 1`, i.e. `C(2n, n)/4ⁿ ≥ 1/(2√n)`. Induction
using `(n + 1) C(2n+2, n+1) = 2(2n + 1) C(2n, n)` and `4n(n + 1) ≤ (2n + 1)²`. -/
theorem central_binom_sq_lower : ∀ n : ℕ, 1 ≤ n → 16 ^ n ≤ 4 * n * (centralBinom n) ^ 2 := by
  intro n hn
  induction n with
  | zero => omega
  | succ m ih =>
    rcases Nat.eq_zero_or_pos m with h | h
    · subst h; simp [centralBinom]
    · have ih' := ih h
      have key := Nat.succ_mul_centralBinom_succ m
      -- multiply the target by (m + 1)², then use the recursion and the induction hypothesis
      have h2 : 16 ^ (m + 1) * (m + 1) ^ 2 ≤
          4 * (m + 1) * (centralBinom (m + 1)) ^ 2 * (m + 1) ^ 2 := by
        have e : 4 * (m + 1) * (centralBinom (m + 1)) ^ 2 * (m + 1) ^ 2 =
            4 * (m + 1) * ((m + 1) * centralBinom (m + 1)) ^ 2 := by ring
        rw [e, key]
        have : 16 ^ (m + 1) * (m + 1) ^ 2 = 16 * 16 ^ m * (m + 1) ^ 2 := by ring
        rw [this]
        calc 16 * 16 ^ m * (m + 1) ^ 2 ≤ 16 * (4 * m * centralBinom m ^ 2) * (m + 1) ^ 2 := by
              gcongr
          _ ≤ 4 * (m + 1) * (2 * (2 * m + 1) * centralBinom m) ^ 2 := by
              have : 4 * m * (m + 1) ≤ (2 * m + 1) ^ 2 := by nlinarith
              nlinarith [Nat.zero_le (centralBinom m ^ 2), this]
      exact Nat.le_of_mul_le_mul_right h2 (by positivity)

/-- **Upper bound** `(2n + 1) · C(2n, n)² ≤ 16ⁿ`, i.e. `C(2n, n)/4ⁿ ≤ 1/√(2n + 1)`. Induction using
the same recursion and `(2n + 1)(2n + 3) ≤ 4(n + 1)²`. -/
theorem central_binom_sq_upper : ∀ n : ℕ, (2 * n + 1) * (centralBinom n) ^ 2 ≤ 16 ^ n := by
  intro n
  induction n with
  | zero => simp
  | succ m ih =>
    have key := Nat.succ_mul_centralBinom_succ m
    have h2 : (2 * (m + 1) + 1) * (centralBinom (m + 1)) ^ 2 * (m + 1) ^ 2 ≤
        16 ^ (m + 1) * (m + 1) ^ 2 := by
      have e : (2 * (m + 1) + 1) * (centralBinom (m + 1)) ^ 2 * (m + 1) ^ 2 =
          (2 * m + 3) * ((m + 1) * centralBinom (m + 1)) ^ 2 := by ring
      rw [e, key]
      have e2 : (2 * m + 3) * (2 * (2 * m + 1) * centralBinom m) ^ 2 =
          4 * (2 * m + 3) * (2 * m + 1) * ((2 * m + 1) * centralBinom m ^ 2) := by ring
      rw [e2]
      calc 4 * (2 * m + 3) * (2 * m + 1) * ((2 * m + 1) * centralBinom m ^ 2)
          ≤ 4 * (2 * m + 3) * (2 * m + 1) * 16 ^ m := by gcongr
        _ ≤ 16 ^ (m + 1) * (m + 1) ^ 2 := by
          have : (2 * m + 3) * (2 * m + 1) ≤ 4 * (m + 1) ^ 2 := by nlinarith
          have h16 : 16 ^ (m + 1) = 16 * 16 ^ m := by ring
          rw [h16]
          nlinarith [Nat.zero_le (16 ^ m)]
    exact Nat.le_of_mul_le_mul_right h2 (by positivity)

/-- `∑ 1/(2n + 1)` diverges: it dominates half the harmonic series, `1/(n + 1) ≤ 2/(2n + 1)`, and
Mathlib knows that `∑ 1/n` diverges (`Real.not_summable_one_div_natCast`). -/
theorem not_summable_one_div_two_mul_add_one :
    ¬ Summable (fun n : ℕ => (1 : ℝ) / (2 * n + 1)) := by
  intro h
  have h2 : Summable (fun n : ℕ => (1 : ℝ) / ((n + 1 : ℕ) : ℝ)) := by
    refine Summable.of_nonneg_of_le (fun n => by positivity) (fun n => ?_) (h.mul_left 2)
    push_cast
    rw [div_le_iff₀ (by positivity)]
    field_simp
    nlinarith
  exact Real.not_summable_one_div_natCast ((summable_nat_add_iff 1).mp h2)

/-- **Pólya in dimension 1, analytically.** The probability that simple random walk on ℤ is at
the origin at time `2n` is `p₂ₙ = C(2n, n)/4ⁿ`, and `∑ p₂ₙ` (the expected number of returns)
diverges, because `p₂ₙ ≥ 1/(2n + 1)`. By the criterion proved in the lesson (Theorem 1), the
walk is recurrent. -/
theorem return_1d_not_summable :
    ¬ Summable (fun n : ℕ => (centralBinom n : ℝ) / 4 ^ n) := by
  intro h
  apply not_summable_one_div_two_mul_add_one
  refine Summable.of_nonneg_of_le (fun n => by positivity) (fun n => ?_) h
  have hle : (4 : ℝ) ^ n ≤ (2 * n + 1) * centralBinom n := by
    exact_mod_cast four_pow_le_central_binom n
  rw [div_le_div_iff₀ (by positivity) (by positivity)]
  linarith

/-- **Pólya in dimension 2, analytically.** On ℤ² the return probability at time `2n` is
`(C(2n, n)/4ⁿ)²` (by the rotation trick below), and `∑ (C(2n, n)/4ⁿ)²` diverges, because by
`central_binom_sq_lower` each term is at least `1/(4n)`. -/
theorem return_2d_not_summable :
    ¬ Summable (fun n : ℕ => ((centralBinom n : ℝ) / 4 ^ n) ^ 2) := by
  intro h
  have h1 : Summable (fun n : ℕ => (1 : ℝ) / ((n + 1 : ℕ) : ℝ)) := by
    have h' := (summable_nat_add_iff 1).mpr h
    refine Summable.of_nonneg_of_le (fun n => by positivity) (fun n => ?_) (h'.mul_left 4)
    have hlow : (16 : ℝ) ^ (n + 1) ≤
        4 * ((n + 1 : ℕ) : ℝ) * (centralBinom (n + 1) : ℝ) ^ 2 := by
      exact_mod_cast central_binom_sq_lower (n + 1) (by omega)
    have h16 : ((4 : ℝ) ^ (n + 1)) ^ 2 = 16 ^ (n + 1) := by
      rw [← pow_mul, mul_comm, pow_mul]; norm_num
    have hpos : (0 : ℝ) < ((n + 1 : ℕ) : ℝ) := by positivity
    simp only [div_pow, h16]
    rw [div_le_iff₀ hpos,
      show 4 * ((centralBinom (n + 1) : ℝ) ^ 2 / 16 ^ (n + 1)) * ((n + 1 : ℕ) : ℝ)
        = 4 * ((n + 1 : ℕ) : ℝ) * (centralBinom (n + 1) : ℝ) ^ 2 / 16 ^ (n + 1) by ring,
      le_div_iff₀ (by positivity)]
    linarith
  exact Real.not_summable_one_div_natCast ((summable_nat_add_iff 1).mp h1)

/-- Yet `p₂ₙ = C(2n, n)/4ⁿ → 0`: by `central_binom_sq_upper`, `p₂ₙ² ≤ 1/(2n + 1)`. The walk on ℤ
returns with probability one, but the chance of being home at any given late time is small — one
reason the expected return time is infinite. -/
theorem return_prob_tendsto_zero :
    Tendsto (fun n : ℕ => (centralBinom n : ℝ) / 4 ^ n) atTop (𝓝 0) := by
  have hsq : Tendsto (fun n : ℕ => ((centralBinom n : ℝ) / 4 ^ n) ^ 2) atTop (𝓝 0) := by
    refine squeeze_zero (fun n => by positivity) (fun n => ?_)
      tendsto_one_div_add_atTop_nhds_zero_nat
    have hup : (2 * (n : ℝ) + 1) * (centralBinom n : ℝ) ^ 2 ≤ 16 ^ n := by
      exact_mod_cast central_binom_sq_upper n
    have h16 : ((4 : ℝ) ^ n) ^ 2 = 16 ^ n := by
      rw [← pow_mul, mul_comm, pow_mul]; norm_num
    rw [div_pow, h16, div_le_div_iff₀ (by positivity) (by positivity)]
    nlinarith [sq_nonneg (centralBinom n : ℝ)]
  have h := hsq.sqrt
  simp only [Real.sqrt_zero] at h
  refine h.congr (fun n => ?_)
  exact Real.sqrt_sq (by positivity)

The rotation trick: the map, the bijection of steps, and the characterisation of closed paths.

/-- The 45° rotation (and scaling by `√2`) `(x, y) ↦ (x + y, x − y)` of ℤ², as an additive map. -/
def rot : ℤ × ℤ →+ ℤ × ℤ where
  toFun v := (v.1 + v.2, v.1 - v.2)
  map_zero' := by simp
  map_add' u v := by ext <;> simp <;> ring

/-- The rotation is injective: from `x + y` and `x − y` one recovers `x` and `y`. -/
theorem rot_injective : Function.Injective rot := by
  intro u v h
  simp only [rot, AddMonoidHom.coe_mk, ZeroHom.coe_mk, Prod.mk.injEq] at h
  ext <;> omega

/-- The four unit steps of simple random walk on ℤ². -/
def unitSteps : Finset (ℤ × ℤ) := {(1, 0), (-1, 0), (0, 1), (0, -1)}

/-- The four diagonal steps `(±1, ±1)`: a pair of independent `±1` steps. -/
def diagSteps : Finset (ℤ × ℤ) := {(1, 1), (-1, -1), (1, -1), (-1, 1)}

/-- The rotation maps the four unit steps bijectively onto the four diagonal steps (checked by
`decide`; injectivity is `rot_injective`). So a uniformly random unit step becomes a pair of
*independent* uniform `±1` steps, and simple random walk on ℤ² becomes two independent simple
random walks on ℤ. -/
theorem rot_unit_steps : unitSteps.image rot = diagSteps := by
  decide

/-- A planar path (a list of steps) returns to the origin iff **both** rotated coordinate walks do.
With `rot_unit_steps` this gives `P(S₂ₙ = 0 in ℤ²) = P(S₂ₙ = 0 in ℤ)² = (C(2n, n)/4ⁿ)²`. -/
theorem closed_iff_rot (L : List (ℤ × ℤ)) :
    L.sum = 0 ↔ ((L.map rot).map Prod.fst).sum = 0 ∧ ((L.map rot).map Prod.snd).sum = 0 := by
  have e1 : ((L.map rot).map Prod.fst).sum = (rot L.sum).1 := by
    rw [map_list_sum rot L]
    exact (map_list_sum (AddMonoidHom.fst ℤ ℤ) _).symm
  have e2 : ((L.map rot).map Prod.snd).sum = (rot L.sum).2 := by
    rw [map_list_sum rot L]
    exact (map_list_sum (AddMonoidHom.snd ℤ ℤ) _).symm
  rw [e1, e2]
  constructor
  · intro h
    simp [h]
  · rintro ⟨h1, h2⟩
    apply rot_injective
    rw [map_zero]
    exact Prod.ext h1 h2

Electrical networks: Nash-Williams for $\Z^2$, and the series behind transience of $\Z^3$.

/-- **The Nash-Williams series for ℤ² diverges.** The edges joining the boundary of the box
`[-k, k]²` to the boundary of `[-(k+1), k+1]²` form a cutset of `8k + 4` unit resistors; the
Nash-Williams inequality gives `R_eff(0 ↔ ∞) ≥ ∑ₖ 1/(8k + 4) = ∞`, so the walk on ℤ² is
recurrent. Proof: compare with the harmonic series. -/
theorem nash_williams_diverges : ¬ Summable (fun k : ℕ => (1 : ℝ) / (8 * k + 4)) := by
  intro h
  have h2 : Summable (fun n : ℕ => (1 : ℝ) / ((n + 1 : ℕ) : ℝ)) := by
    refine Summable.of_nonneg_of_le (fun n => by positivity) (fun n => ?_) (h.mul_left 8)
    push_cast
    rw [div_le_iff₀ (by positivity)]
    field_simp
    nlinarith
  exact Real.not_summable_one_div_natCast ((summable_nat_add_iff 1).mp h2)

/-- `∑ 1/n²` converges. In ℤ³ a unit flow from the origin spread evenly over spheres has energy
comparable to `∑ₖ k² · (1/k²)² = ∑ 1/k²`, finite, so by T. Lyons' criterion the walk on ℤ³ is
transient. This applies Mathlib's `Real.summable_one_div_nat_pow`. -/
theorem summable_inv_sq : Summable (fun n : ℕ => (1 : ℝ) / (n : ℝ) ^ 2) :=
  Real.summable_one_div_nat_pow.mpr one_lt_two

Peierls' argument, with the path count of Lemma 4 as a hypothesis:

/-- **Peierls' estimate.** If there are at most `cₙ ≤ 4 · 3ⁿ` self-avoiding paths of length `n`
from the origin and each is open with probability `pⁿ`, then for `p < 1/3` the expected number of
open such paths, at most `cₙ pⁿ ≤ 4(3p)ⁿ`, tends to `0`. So the origin lies in an infinite open
cluster with probability `0`, and `p_c(ℤ²) ≥ 1/3`. (The count `cₙ ≤ 4 · 3^{n-1}` is proved in the
lesson; here it is a hypothesis.) -/
theorem peierls_tendsto {c : ℕ → ℕ} (hc : ∀ n, c n ≤ 4 * 3 ^ n) {p : ℝ} (hp0 : 0 ≤ p)
    (hp : p < 1 / 3) : Tendsto (fun n => (c n : ℝ) * p ^ n) atTop (𝓝 0) := by
  have h3 : 3 * p < 1 := by linarith
  have hlim : Tendsto (fun n : ℕ => 4 * (3 * p) ^ n) atTop (𝓝 0) := by
    simpa using (tendsto_pow_atTop_nhds_zero_of_lt_one (by linarith) h3).const_mul 4
  refine squeeze_zero (fun n => by positivity) (fun n => ?_) hlim
  have : (c n : ℝ) ≤ 4 * 3 ^ n := by exact_mod_cast hc n
  calc (c n : ℝ) * p ^ n ≤ (4 * 3 ^ n) * p ^ n := by gcongr
    _ = 4 * (3 * p) ^ n := by ring

/-- **The dual half of Peierls' argument.** The number of dual circuits of length `n` surrounding
the origin is at most `n · 4 · 3ⁿ`, and each is closed with probability `qⁿ`, `q = 1 − p`. For
`q < 1/3` the series `∑ n · 4 · 3ⁿ qⁿ` converges, so its tails tend to `0`; this gives
`p_c(ℤ²) ≤ 2/3`. Uses Mathlib's `summable_pow_mul_geometric_of_norm_lt_one`. -/
theorem dual_circuits_summable {q : ℝ} (hq0 : 0 ≤ q) (hq : q < 1 / 3) :
    Summable (fun n : ℕ => (n : ℝ) * 4 * 3 ^ n * q ^ n) := by
  have h := summable_pow_mul_geometric_of_norm_lt_one 1 (r := 3 * q)
    (by rw [Real.norm_eq_abs, abs_lt]; constructor <;> linarith)
  refine (h.mul_left 4).congr (fun n => ?_)
  ring

And the subcritical random graph:

/-- **No giant component for `c ≤ 1`.** The fraction `β` of vertices of `G(n, c/n)` in the giant
component satisfies, in the limit, the survival equation of a Poisson(`c`) branching process,
`β = 1 − e^{−cβ}`. For `c ≤ 1` its only solution `β ≥ 0` is `β = 0`, because `1 − e^{−x} < x` for
`x > 0` (from Mathlib's `Real.add_one_lt_exp`). For `c > 1` a second, positive root appears. -/
theorem no_giant_subcritical {c β : ℝ} (hc : c ≤ 1) (hc0 : 0 ≤ c) (hβ : 0 ≤ β)
    (h : β = 1 - Real.exp (-(c * β))) : β = 0 := by
  by_contra hne
  have hpos : 0 < β := lt_of_le_of_ne hβ (Ne.symm hne)
  rcases eq_or_lt_of_le hc0 with hc0' | hc0'
  · subst hc0'
    simp at h
    exact hne h
  · have hx : 0 < c * β := mul_pos hc0' hpos
    have := Real.add_one_lt_exp (show -(c * β) ≠ 0 by linarith)
    nlinarith

Not formalised: the random walk as a stochastic process and Theorem 1; the identification of $P(S_{2n} = 0)$ with the counts used here; the $d = 3$ Stirling estimate; electrical networks and the Nash-Williams inequality itself; Lemma 4 and planar duality; Kesten's theorem; the Erdős–Rényi theorem; Brownian motion's path properties and SLE. Formalising planar duality for percolation would need a discrete Jordan curve theorem — a good project for an ambitious learner.

The first experiment checks the rotation trick by brute force (all $4^{2n}$ walks of length $2n \le 8$) and then sums exact return probabilities. In $d = 3$ it uses the formula obtained by conditioning on the number of vertical steps, and compares the partial sums — and a tail correction based on $u_{2n} \sim 2\big(\frac{3}{4\pi n}\big)^{3/2}$ — with Watson's closed form $G = \frac{\sqrt6}{32\pi^3}\Gamma(\frac1{24})\Gamma(\frac5{24})\Gamma(\frac7{24})\Gamma(\frac{11}{24})$ (Watson, 1939, in the form given by Glasser and Zucker, 1977).

"""Polya's theorem by exact return probabilities.

Simple random walk on Z^d is recurrent exactly when the expected number of
returns to the origin, G = sum_n u_2n with u_2n = P(S_2n = 0), is infinite.
We compute u_2n exactly or to floating-point accuracy:

    d = 1:  u_2n = C(2n, n) / 4^n
    d = 2:  u_2n = (C(2n, n) / 4^n)^2          (the 45-degree rotation trick)
    d = 3:  condition on the number 2j of vertical steps:
            u_2n = sum_j C(2n, 2j) (1/3)^(2j) (2/3)^(2n-2j) u1_2j u2_(2n-2j)

and watch the partial sums G(N): they grow like sqrt(N) in d = 1, like
log(N)/pi in d = 2, and converge in d = 3, to Watson's constant
G = sqrt(6)/(32 pi^3) Gamma(1/24) Gamma(5/24) Gamma(7/24) Gamma(11/24),
so the walk returns with probability 1 - 1/G = 0.3405... We also check the
rotation trick exhaustively for short walks.

Evidence: exact counting for short walks (exhaustive-finite) and
floating-point summation of exact formulas (numerical). Finite partial sums
cannot prove divergence; the proofs are in the lesson and in Lean.
"""

import itertools
import math

STEPS2 = [(1, 0), (-1, 0), (0, 1), (0, -1)]

def closed_walks_2d(n):
    """Number of walks of length 2n on Z^2 that end at the origin (brute force)."""
    count = 0
    for walk in itertools.product(STEPS2, repeat=2 * n):
        if sum(s[0] for s in walk) == 0 and sum(s[1] for s in walk) == 0:
            count += 1
    return count

def u1(n_max):
    """[u_0, u_2, ..., u_2n_max] for d = 1, via u_2n+2 = u_2n (2n+1)/(2n+2)."""
    u = [1.0]
    for n in range(n_max):
        u.append(u[-1] * (2 * n + 1) / (2 * n + 2))
    return u

def u3(n_max, one=None):
    """Return probabilities u_2n of simple random walk on Z^3, n <= n_max."""
    one = one if one is not None else u1(n_max)
    two = [x * x for x in one]
    out = []
    for n in range(n_max + 1):
        total = 0.0
        # term for j vertical pairs: C(2n,2j) (1/3)^2j (2/3)^(2n-2j) one[j] two[n-j]
        log_c = 0.0  # log C(2n, 2j), updated multiplicatively
        for j in range(n + 1):
            if j > 0:
                a = 2 * n - 2 * j + 2
                log_c += math.log(a * (a - 1) / ((2 * j) * (2 * j - 1)))
            weight = math.exp(log_c + 2 * j * math.log(1 / 3)
                              + (2 * n - 2 * j) * math.log(2 / 3))
            total += weight * one[j] * two[n - j]
        out.append(total)
    return out

def watson_constant():
    """G(Z^3) = sqrt(6)/(32 pi^3) Gamma(1/24)Gamma(5/24)Gamma(7/24)Gamma(11/24)."""
    g = math.gamma
    return (math.sqrt(6) / (32 * math.pi ** 3)
            * g(1 / 24) * g(5 / 24) * g(7 / 24) * g(11 / 24))

def partial_sums(u, checkpoints):
    """{N: sum_{n=0}^{N} u_2n} at the given N."""
    out, s = {}, 0.0
    for n, x in enumerate(u):
        s += x
        if n in checkpoints:
            out[n] = s
    return out

def main():
    print("Exhaustive check of the rotation trick on Z^2:")
    for n in range(1, 5):
        print(f"  closed walks of length {2 * n}: {closed_walks_2d(n):>5}"
              f"   C({2 * n},{n})^2 = {math.comb(2 * n, n) ** 2:>5}")
    big, small = 100000, 1200
    one = u1(big)
    g1 = partial_sums(one, {100, 1000, 10000, 100000})
    g2 = partial_sums([x * x for x in one], {100, 1000, 10000, 100000})
    g3 = partial_sums(u3(small, one), {100, 300, 1200})
    print()
    print("Expected number of visits to 0 up to time 2N, G(N) = sum u_2n:")
    print(f"{'N':>8}{'d = 1':>11}{'2 sqrt(N/pi)':>14}{'d = 2':>9}"
          f"{'log(N)/pi':>11}")
    for n in sorted(g1):
        print(f"{n:>8}{g1[n]:>11.2f}{2 * math.sqrt(n / math.pi):>14.2f}"
              f"{g2[n]:>9.3f}{math.log(n) / math.pi:>11.3f}")
    c3 = 4 * (3 / (4 * math.pi)) ** 1.5  # tail ~ c3 / sqrt(N)
    print(f"{'N':>8}{'d = 3':>11}{'+ tail est.':>14}")
    for n in sorted(g3):
        print(f"{n:>8}{g3[n]:>11.5f}{g3[n] + c3 / math.sqrt(n):>14.5f}")
    w = watson_constant()
    print(f"Watson's exact value G = {w:.7f}; return probability "
          f"1 - 1/G = {1 - 1 / w:.6f}")
    print()
    print("Epistemic status: the first table is an exhaustive count (exact);")
    print("the others sum exact formulas in floating point. Partial sums that")
    print("grow, or seem to settle, cannot by themselves prove divergence or")
    print("convergence; the d = 1, 2 divergence is proved in the lesson and in")
    print("Lean, and d = 3 convergence in the lesson (by an estimate).")

if __name__ == "__main__":
    main()

The second experiment estimates by simulation the probability of an open left–right crossing of the self-dual $(L + 1) \times L$ rectangle for several $p$, showing the curve sharpening around $\frac12$ as $L$ grows — and staying at $\frac12$ at $p = \frac12$, as self-duality requires. Then it measures the largest component of random graphs with $2000$ vertices and average degree $c$, against the branching-process prediction $\beta(c)$.

"""Phase transitions: bond percolation on a grid and the giant component.

Part 1. Bond percolation on a square grid: each edge is open with
probability p. We estimate by simulation the probability that open edges
connect the left side of an (L + 1) x L rectangle to its right side. As L
grows the curve sharpens into a step at p = 1/2, Kesten's critical value;
at p = 1/2 the crossing probability is exactly 1/2 for every L, because the
rectangle is chosen to be self-dual.

Part 2. The Erdos-Renyi random graph with n vertices and m = c n / 2 random
edges (asymptotically the same as G(n, c/n)). The largest component has a
vanishing fraction of the vertices when c < 1 and a fraction close to the
positive root beta of beta = 1 - exp(-c beta) when c > 1.

Evidence: statistical (Monte Carlo with a fixed seed, small sizes and few
trials, with approximate standard errors). Simulations near a critical
point estimate behaviour; the exact critical values are theorems.
"""

import math
import random

def find(parent, x):
    """Root of x in a union-find forest, with path halving."""
    while parent[x] != x:
        parent[x] = parent[parent[x]]
        x = parent[x]
    return x

def union(parent, size, a, b):
    """Merge the components of a and b (union by size)."""
    ra, rb = find(parent, a), find(parent, b)
    if ra == rb:
        return
    if size[ra] < size[rb]:
        ra, rb = rb, ra
    parent[rb] = ra
    size[ra] += size[rb]

def crosses(n, p, rng):
    """Is there an open left-right crossing of the self-dual rectangle R_n?

    R_n has n interior columns x = 1..n and n + 1 rows y = 0..n. Each row
    has n + 1 horizontal edges, from the left side (x = 0) through the
    interior to the right side (x = n + 1), and each interior column has n
    vertical edges. Its planar dual is a rotated copy of R_n, so at p = 1/2
    a crossing occurs with probability exactly 1/2.
    """
    count = n * (n + 1) + 2
    parent, size = list(range(count)), [1] * count
    left, right = count - 2, count - 1

    def node(x, y):
        return y * n + (x - 1)

    for y in range(n + 1):
        if rng.random() < p:
            union(parent, size, left, node(1, y))
        if rng.random() < p:
            union(parent, size, node(n, y), right)
        for x in range(1, n):
            if rng.random() < p:
                union(parent, size, node(x, y), node(x + 1, y))
    for x in range(1, n + 1):
        for y in range(n):
            if rng.random() < p:
                union(parent, size, node(x, y), node(x, y + 1))
    return find(parent, left) == find(parent, right)

def crossing_probability(n, p, trials, rng):
    """Fraction of trials with an open left-right crossing."""
    return sum(crosses(n, p, rng) for _ in range(trials)) / trials

def largest_component_fraction(n, c, rng):
    """Largest component of a random graph with n vertices and c n / 2 edges."""
    parent, size = list(range(n)), [1] * n
    for _ in range(round(c * n / 2)):
        a, b = rng.randrange(n), rng.randrange(n)
        if a != b:
            union(parent, size, a, b)
    return max(size[find(parent, v)] for v in range(n)) / n

def giant_fraction(c, iterations=200):
    """Largest root of beta = 1 - exp(-c beta) in [0, 1].

    For c <= 1 the only root is 0 (proved in the lesson's Lean file); for
    c > 1 iterating beta -> 1 - exp(-c beta) from beta = 1 converges to the
    positive root.
    """
    if c <= 1:
        return 0.0
    beta = 1.0
    for _ in range(iterations):
        beta = 1 - math.exp(-c * beta)
    return beta

def main():
    rng = random.Random(1957)
    ps = [0.35, 0.45, 0.5, 0.55, 0.65]
    trials = 200
    print(f"Part 1: P(left-right open crossing), {trials} trials per entry")
    print("  p    " + "".join(f"{p:>8}" for p in ps))
    for n in [8, 24]:
        est = [crossing_probability(n, p, trials, rng) for p in ps]
        print(f"  L = {n:<3}" + "".join(f"{e:>8.2f}" for e in est))
    print(f"  (standard error of each entry at most about "
          f"{0.5 / math.sqrt(trials):.2f})")
    print()
    n, reps = 2000, 5
    print(f"Part 2: random graph, n = {n} vertices, average degree c "
          f"({reps} samples)")
    print(f"  {'c':>4}{'largest component / n':>24}{'beta(c) theory':>17}")
    for c in [0.5, 0.9, 1.5, 2.0, 3.0]:
        sim = sum(largest_component_fraction(n, c, rng)
                  for _ in range(reps)) / reps
        print(f"  {c:>4}{sim:>24.3f}{giant_fraction(c):>17.3f}")
    print("  (for c < 1 the largest component has O(log n) vertices,")
    print("  so its fraction tends to 0 only slowly as n grows)")
    print()
    print("Epistemic status: Monte Carlo estimates from a few hundred random")
    print("samples on small systems (statistical evidence). They are")
    print("consistent with Kesten's theorem p_c = 1/2 and the Erdos-Renyi")
    print("theorem, but cannot determine critical values or exponents exactly.")

if __name__ == "__main__":
    main()

The first experiment is exhaustive (exact) for short walks and numerical for the sums; the second is statistical, with a few hundred samples on small systems. A simulation near a critical point can suggest $p_c$ and estimate exponents; it cannot determine them exactly — Kesten's theorem and the SLE computations do that.

Many one-dimensional walks at once, with their root-mean-square distance against the exact $\sqrt n$. Notice how far a typical walk wanders while still being bound to come home: that tension is why the expected return time is infinite.

Kirchhoff's theorem (1847) connects electrical networks to spanning trees: with unit resistors, the effective resistance between the endpoints of an edge $e$ equals the fraction of spanning trees that contain $e$. Build spanning trees of this $3 \times 3$ grid; the widget counts all of them ($192$) by the matrix-tree theorem. Uniform spanning trees, whose branches are loop-erased random walks (Wilson's algorithm), converge to $\mathrm{SLE}_2$ and $\mathrm{SLE}_8$.