A particle dropped into a triangular room travels in straight lines
and bounces off the walls, with angle of incidence equal to angle of
reflection. "Ergodic" means that a typical path, run long enough,
spreads evenly: the fraction of time it spends in any region equals
that region's share of the table's area, and no subregion traps it
forever. That is exactly what happens whenever at least one angle of
the triangle is an irrational multiple of pi.

You only need basic triangle geometry and trigonometry, the idea
of a Fourier series (reviewed when used), and the phrase "almost
every" (defined when first used). Every technical step is proved;
nothing is left as an exercise.

The theorem. The billiard flow in any nondegenerate Euclidean
triangle is ergodic for normalized area x uniform angle measure,
provided at least one angle divided by pi is irrational.

The proof, in five sentences.
  1. Glue the triangle to its mirror image along the boundary to
     get a flat surface M -- the "double" -- with three cone points.
     Straight-line motion on M equals the billiard on Q.
  2. Fourier-expand any bounded invariant f(p, theta) in the angular
     variable. Get coefficient functions f_j on M coupled by a
     first-order differential recurrence.
  3. Smooth f spatially, away from the cone points. An "energy
     estimate" bounds every coefficient gradient by the supremum of
     f alone.
  4. Average f over long segments to get an actually invariant
     approximation g. Two limits (smoothing scale -> 0, segment
     length -> 0) force every coefficient gradient to vanish.
  5. Each f_j is then a parallel section of a flat line bundle.
     Going around a cone of angle 2 alpha multiplies f_j by
     exp(-2 i j alpha). If alpha/pi is irrational, only j = 0
     survives; connectedness of M makes f_0 constant.

Contents (follows the five-sentence proof).
  Part 0 (steps 1-5):    billiards, ergodicity, theorem, history.
  Part 1 (steps 6-12):   the double; X and Y; the reduction to
                         "no nonconstant bounded invariant f".
  Part 2 (steps 13-19):  Fourier modes; the endpoint equation;
                         the low-regularity obstacle.
  Part 3 (steps 20-25):  three estimates (smoothing, energy,
                         bounded-modes).
  Part 4 (steps 26-35):  two-limit rigidity argument.
  Part 5 (steps 36-39):  holonomy; conclusion of the theorem.
  Part 6 (steps 40-46):  consequences, caveats, glossary,
                         compressed proof-at-a-glance.

============================================================
PART 0. THE PROBLEM
============================================================

STEP 1. What is a billiard?
A point particle moves at constant unit speed inside a polygonal
region Q, reflecting specularly at the boundary (angle of
incidence = angle of reflection). Between hits the path is a
straight line. Here Q is the interior of a nondegenerate Euclidean
triangle.

A state is a pair (p, theta) with p in Q, theta in S^1 = [0, 2 pi),
recording both position and direction. The billiard flow Phi_t maps
"state now" to "state t seconds later". Vertex hits are excluded
(STEP 8); after that exclusion, the flow is defined for every
real t.

STEP 2. What does "ergodic" mean?
A subset A of state space is invariant if a trajectory that starts
inside A stays inside A forever, and one that starts outside A
stays outside A. The flow is ergodic if every measurable invariant
set has measure 0 or measure 1 (the whole space).

Equivalent (Birkhoff 1931): for any continuous F on state space,
for almost every starting state, the time average of F along the
trajectory equals the space average of F. Ergodicity is a
statistical statement about typical trajectories, not all.

Ergodicity does NOT mean every orbit is dense; it does NOT mean
the flow mixes. It means exactly: any bounded measurable invariant
function must be constant almost everywhere (STEP 12).

STEP 3. The theorem.
Let Q be a nondegenerate triangle with angles alpha, beta, gamma,
alpha + beta + gamma = pi. Suppose at least one of alpha/pi,
beta/pi, gamma/pi is irrational. Discard the measure-zero set of
initial states whose trajectory hits a vertex in either time
direction. Then Phi_t is ergodic on Q x S^1 for the probability
measure

    dmu(p, theta)  =  dArea(p) dtheta / (2 pi Area(Q)),

where dArea is 2-d Lebesgue area and dtheta is uniform angle on
S^1. Equivalently: if A is measurable and mu(Phi_t^{-1} A
symmetric-difference A) = 0 for every real t, then mu(A) is 0 or 1.

STEP 4. Which triangles does the theorem cover?
The ratios alpha/pi, beta/pi, gamma/pi sum to 1. Rational plus
irrational is irrational; rational plus rational is rational. So
the number of irrational ratios is 0, 2, or 3 -- never 1.

Covered (the theorem applies):
  alpha/pi = 1/4,  beta/pi = 1/e,  gamma/pi = 3/4 - 1/e.
  alpha/pi = 1/sqrt(10),  beta/pi = 1/sqrt(3) - 1/sqrt(5),
      gamma/pi = 1 - 1/sqrt(10) - 1/sqrt(3) + 1/sqrt(5).
Excluded (all three ratios rational):
  equilateral (1/3, 1/3, 1/3); isosceles right (1/4, 1/4, 1/2);
  30-60-90 (1/6, 1/3, 1/2); any rational-multiple-of-pi triangle.

STEP 5. What was known before, and why is this new?
For the (excluded) all-rational case: the flow unfolds to a
straight-line flow on a compact flat surface, KMS 1986 gives
unique ergodicity in almost every direction, and some directions
are genuinely non-ergodic -- so the full flow is not ergodic. For
irrational tables, previous work gave: ergodicity for a dense
G_delta (a Baire category statement, which names no particular
triangle); a faster-than-any-Diophantine approximation condition
(Vorobets 1997) with explicit examples; weak mixing for a dense
G_delta of n-gons (Chaika-Forni 2026). The theorem of STEP 3 is
the first pointwise statement: every irrational triangle, no
genericity, no approximation-speed assumption.

============================================================
PART 1. THE DOUBLE
============================================================

STEP 6. Why can't we just unfold a triangle with an irrational angle?
Unfolding a polygonal billiard: when the trajectory hits a side,
do not reflect; instead continue it straight into a mirrored copy
of the polygon across that side, then again across the next side,
forever. For a rational polygon (all angles rational multiples of
pi), the group generated by reflections is finite, so finitely many
copies suffice and the folded path becomes a straight line on a
compact flat surface. KMS 1986 theory applies.

For a triangle with an irrational angle, the reflection group is
infinite (in fact dense in O(2)). No finite family of reflected
copies contains an unfolded path. The KMS machinery is not
available.

STEP 7. What replaces unfolding?
Exactly TWO copies. Glue each side of one copy to the corresponding
side of the other copy, with opposite orientations. The result M
is a closed, connected flat surface with no edges and no boundary,
and it looks like a piece of the plane everywhere except at three
"cone points" above the three vertices of Q. This is the "double"
of Q. It is compact whatever the angles.

STEP 8. Cone angles, and why vertex hits do not matter.
At the vertex of Q with angle alpha, M has total angle 2 alpha --
two triangle corners of angle alpha glued back to back. Similarly
2 beta and 2 gamma. Each cone angle lies strictly between 0 and
2 pi. Gauss-Bonnet check (curiosity): the three defects
(2 pi - 2 alpha) + (2 pi - 2 beta) + (2 pi - 2 gamma)
= 6 pi - (2 alpha + 2 beta + 2 gamma) = 6 pi - 2 pi = 4 pi,
matching a sphere's 4 pi of total curvature. So M is
homeomorphic to a sphere with three marked cone points.

Vertex hits. A billiard trajectory in Q hits a vertex iff, for
some finite sequence of sides, an unfolded straight line passes
exactly through a fixed point. For each fixed side sequence this
is one scalar condition on the three parameters of state, so it is
a codimension-1 set. The countable union over side sequences is
null. We discard it. On the complement, the billiard flow on Q
lifts to a smooth unit-speed geodesic flow on SM (STEP 9) that
preserves Liouville measure.

STEP 9. Phase space and invariant measure.
The phase space is SM = M x S^1, a state being (p, theta) with
p on M and theta in S^1 the direction of motion. The flow
Phi_t(p, theta) = (the point t away from p in direction theta,
theta), well-defined a.e. by STEP 8. The invariant probability
measure is dA(p) dtheta, normalized to total mass 1, where dA is
area on M and dtheta is uniform angle. The billiard flow on Q and
the geodesic flow on SM are measure-preservingly conjugate away
from the null vertex set.

STEP 10. Complex derivatives in a flat chart.
Locally, on M in a flat chart with coordinates z = x + i y,
    d     = d/dz    = (1/2)(d/dx - i d/dy),
    dbar  = d/dzbar = (1/2)(d/dx + i d/dy).
Both are complex operators, real on real parts, well-defined up
to the chart's rotation phase (STEP 14).

STEP 11. The vector fields X and Y.
On SM, in a flat chart, define two smooth real vector fields
(horizontal; no d/dtheta component):
    X  =  e^{i theta} d   + e^{-i theta} dbar,
    Y  =  i e^{i theta} d - i e^{-i theta} dbar.
In real coordinates,
    X =  cos(theta) d/dx + sin(theta) d/dy,
    Y = -sin(theta) d/dx + cos(theta) d/dy.
The theta-dependent complex coefficients combine into two real
vector fields on SM. X is the infinitesimal generator of the
geodesic flow Phi_t. Y is X rotated by -90 degrees. In a flat
chart, [X, Y] = 0.

STEP 12. The ergodicity reduction.
Let f be bounded and measurable on SM. By STEP 8, Phi_t preserves
the Liouville measure and is defined a.e. The word "invariant"
means: for every real t, f o Phi_t = f a.e. Formally, invariance
is equivalent to the distributional identity

    X f = 0  on the regular part of M x S^1.

Ergodicity of Phi_t is equivalent to: every bounded measurable
invariant f is constant a.e. (STEP 2 plus Fubini along the flow).

So the whole theorem reduces to: prove that
    f in L^inf(SM) with Xf = 0 in distributions on M x S^1
    is constant almost everywhere.

The proof of that is Parts 2-5. The whole content of the argument
is: f is only in L^inf, so it has no derivatives a priori. The
recurrence Xf = 0 says only that some combination of spatial
derivatives of neighboring Fourier coefficients add up to zero.
We will upgrade this to "each coefficient has vanishing gradient",
by regularizing in space and averaging in time.

============================================================
PART 2. FOURIER MODES
============================================================

STEP 13. Fourier expansion in the angular variable.
The flow never changes theta in a flat chart (that is what "flat"
means geometrically). So for the purpose of the differential
analysis below, theta is just a parameter. But theta lives on a
circle, so every L^2 function on SM can be Fourier-expanded in
theta:
    f(p, theta) = sum_{j in Z} f_j(p) e^{i j theta}.
The coefficient function is
    f_j(p) = (1/2 pi) int_0^{2 pi} f(p, theta) e^{-i j theta} dtheta.
Parseval reads ||f||_{L^2(SM)}^2 = sum_j ||f_j||_{L^2(M)}^2.
For f in L^inf, the sum still converges and |f_j| <= ||f||_inf.

STEP 14. The coefficient f_j is not a function; it is a section.
The angle theta is not globally defined on M. Around a cone point,
the flat frame rotates. Concretely, if a chart change multiplies
z by e^{i phi}, then theta becomes theta + phi, so the coefficient
transforms as
    f_j  ->  e^{-i j phi} f_j.
In other words, each f_j is a section of a flat unitary line bundle
L^{-j}. Under chart changes, the coefficient picks up the unitary
phase e^{-i j phi}. Only f_0 is a genuine function on M.
Consequence: d f_j is a section of L^{-(j+1)} and dbar f_j is a
section of L^{-(j-1)}. This twist is the source of the mode
identities in STEP 15.

STEP 15. Mode identities for X and Y.
Apply X and Y to the Fourier expansion. On u = u_j(p) e^{i j theta}:
    X u = e^{i(j+1)theta} (d u_j) + e^{i(j-1)theta} (dbar u_j),
    Y u = i e^{i(j+1)theta} (d u_j) - i e^{i(j-1)theta} (dbar u_j).
Reading off the n-th Fourier coefficient:
    (Xu)_n  =  d u_{n-1}  +  dbar u_{n+1},
    (Yu)_n  =  i d u_{n-1}  -  i dbar u_{n+1}.
Each of X and Y shifts the index by +/-1. Adjacent coefficients
are coupled through derivatives. The versions above are also valid
distributionally on the regular part of M.

STEP 16. The recurrence from X f = 0.
By STEP 12, Xf = 0. Reading off the n-th Fourier coefficient of
Xf = 0 gives
    d f_{n-1}  +  dbar f_{n+1}  =  0,   for every n in Z.
Equivalently, for every k in Z,
    d f_k  =  -  dbar f_{k+2}.
This single linear recurrence encodes the flow invariance
completely. It is a relation between first derivatives of
coefficients two steps apart.

STEP 17. Why is X f = 0 not enough on its own?
The recurrence d f_k = -dbar f_{k+2} expresses d of one coefficient
in terms of dbar of another. In principle it could be satisfied by
nonconstant coefficients: pick any function, differentiate, and
the recurrence adjusts the neighbor. Formal series solutions
exist whose coefficients never vanish.

We need to prove that, for f in L^inf satisfying Xf = 0, every
coefficient derivative actually vanishes. That is a rigidity
statement that Xf = 0 alone does not give: the coefficients have
no a priori regularity beyond L^2, and the recurrence is a
first-order equation with no obvious boundary conditions.

The argument below supplies the missing rigidity in three
mechanisms, applied in succession:
  (i) smoothing in space, to be able to differentiate,
  (ii) long-time averaging, to produce an intermediate object
       that is exactly flow-invariant,
  (iii) an energy estimate that controls individual coefficient
       gradients by ||f||_inf alone.
These are developed in Parts 3-4.

STEP 18. The tail U and the endpoint B.
Fix an integer m. Let P_j denote projection onto the j-th Fourier
mode. Define
    U :=  sum_{k >= 0} P_{m+2k} f,     (the "tail")
    B :=  e^{i(m-1)theta} dbar f_m.    (the "endpoint")
The series for U converges in L^2(SM) since the f_j are square-
summable. B is a single Fourier mode.

Claim:  XU  =  B,   distributionally.

Proof of the claim. Truncate:
    U_N := sum_{k=0}^{N} P_{m+2k} f.
Then by STEP 15 (with n = m+1, m+3, ..., m+2N-1):
    X U_N  =  e^{i(m-1)theta} dbar f_m
            + sum_{n = m+1, m+3, ..., m+2N-1} (d f_{n-1} + dbar f_{n+1})
            + e^{i(m+2N+1)theta} d f_{m+2N}.
By the recurrence (STEP 16) each bracketed term vanishes, so
    X U_N  =  B  +  R_N,
where R_N := e^{i(m+2N+1)theta} d f_{m+2N}. The term R_N is a
single Fourier mode, its L^2 norm bounded by C ||f||_inf (STEP 25,
bounded-modes lemma). Its angular mode m+2N+1 tends to +infinity
with N. So for any fixed L^2 test field phi,
    ||  <=  C ||f||_inf ||P_{m+2N+1} phi||_2  -->  0
as N -> +infinity. Thus R_N -> 0 weakly in L^2.
Since U_N -> U strongly in L^2, the identity X U_N = B + R_N
passes to the distributional limit X U = B.

STEP 19. Why proving B = 0 finishes the whole proof.
By STEP 16 with k = m-2:
    d f_{m-2}  =  -  dbar f_m.
Now compute (Yf)_{m-1} using STEP 15:
    (Yf)_{m-1}  =  i d f_{m-2}  -  i dbar f_m
                =  i(-dbar f_m) - i dbar f_m
                =  -2 i dbar f_m.
So as single-mode (m-1) fields,
    P_{m-1} Y f  =  -2 i B.
Hence
      =    =  -2 i ||B||_2^2.
If we can prove  = 0 for arbitrary m, then
||B||_2^2 = 0, so B = 0, so dbar f_m = 0 for arbitrary m, so
dbar f_k = 0 for all k. Then the recurrence gives d f_k = 0 for
all k. Each f_j is locally constant on M, and M is connected, so
each f_j is constant. That is the rigidity proposition.

The one target left is therefore: prove  = 0
for every m. This is done in Part 4, after the three estimates
in Part 3.

============================================================
PART 3. THREE ESTIMATES
============================================================

Let r(x) = distance from x in M to the nearest cone point. There
exists a constant K > 0 such that, for every eps in (0, 1), every
x with r(x) > K eps has an embedded flat disk of radius 4 eps.
This is a compactness statement: away from cones the embedded
radius is bounded below; near a cone of angle theta_c the embedded
radius of a disk at distance rho from the apex is c(theta_c) rho
for some c > 0.

Define the safe region S_eps := { r > K eps }.
Also: for each fixed constant L,
    Area( { r <= L eps } ) = O_L( eps^2 ).
This will localize all our error terms to shrinking cone annuli.

STEP 20. The spatial smoothing operator S_eps.
Choose a nonnegative smooth radial bump rho on R^2, supported in
the unit disk, with total mass 1. Set rho_eps(y) = eps^{-2}
rho(y/eps). For x in the safe region, define
    (S_eps u)(x, theta)
        :=  int_{R^2} rho_eps(y) u(x - y, theta) dy,
the convolution taken in an embedded flat chart (direction theta
frozen, only positions averaged). The definition is independent of
the choice of flat chart, by radiality of rho.

So S_eps averages over a disk of radius eps, keeping theta fixed.
It smooths positions, not directions, and it avoids the cone
points because it only acts on the safe region.

STEP 21. Statement of the smoothing lemma.
For all sufficiently small eps > 0:
  (a) Bounds on the safe region:
        ||S_eps u||_{L^2(r > K eps)} <= ||u||_{L^2(SM)},
        ||eps grad S_eps u||_{L^2(r > K eps)} <= C ||u||_{L^2(SM)}.
      Also ||S_eps u||_inf <= ||u||_inf.
  (b) Small scaled derivative: for each fixed u in L^2(SM),
        ||eps Y S_eps u||_{L^2(r > K eps)}  -->  0  as eps -> 0.
  (c) Convergence: on compact subsets of the safe region,
        S_eps u -> u in L^2.
  (d) Support and adjoint: fix L > K+4. If g in L^2(SM) vanishes
      on { r < L eps }, then S_eps g is supported in { r >=
      (L-1) eps } and extends by zero to a compactly supported
      field on M. For u, g in L^2(SM),
          =  -  .
  (e) Commutation with X: if Xg in L^2 distributionally, then
        X (S_eps g)  =  S_eps (Xg),
        ||X S_eps g||_{L^2}  <=  C ||Xg||_{L^2},
      with the same zero-extension convention as in (d).

STEP 22. Proof of the smoothing lemma.
In a flat chart containing the convolution disk, S_eps is ordinary
spatial convolution at fixed angle. Write the kernel
K_eps(x, y) = rho_eps(x - y).

(a) The L^1 mass of K_eps in either variable is 1; that gives
the first bound by Jensen or Young's inequality. The kernel of
eps * (a first derivative) * rho_eps has L^1 mass O(1); that gives
the second bound. The sup bound is immediate since |S_eps u| <=
S_eps |u| <= ||u||_inf.

(b) Key observation: the L^1 mass of (eps * first derivative of
rho_eps) is O(eps). Indeed, rho_eps = eps^{-2} rho(y/eps), so a
first derivative w.r.t. y_i equals eps^{-3} (d rho/dy_i)(y/eps),
which has L^1 mass eps^{-1}. Multiplying by eps gives L^1 mass
O(eps), so
    ||eps Y S_eps u||_{L^2} <= ||eps Y rho_eps||_{L^1} ||u||_2
                             = O(eps) * ||u||_2  -->  0.
This gives (b) directly, for all u in L^2, without invoking
density.

(c) Ordinary approximate-identity argument.

(d) Support: if the center x is more than eps away from the
support of g in the "shrinking" sense, the convolution vanishes by
the support of rho_eps. Triangle inequality for r gives the
claimed support { r >= (L-1) eps }. Combined with the support of
g (away from cones), the extension by zero creates no boundary
derivative.

Adjoint identity: differentiate K_eps(x, y) w.r.t. the first
variable. Since K_eps(x, y) = rho_eps(x - y), we have
Y_x K_eps(x, y) = -Y_y K_eps(x, y) in any flat chart. Integrating
by parts in y (both u and g compactly supported in the relevant
sense):
      =  int int (Y_x K_eps(x, y)) u(y) g(x) dy dx
                    =  -int int K_eps(x, y) u(y) (Y_y g)(y) dy dx
                    =  -.
No derivative of g is required beyond existence a.e. in the
distributional pairing.

(e) Local commutation in flat charts: X S_eps = S_eps X, so
X(S_eps g) = S_eps(Xg) distributionally. The L^2 bound follows
from (a).

STEP 23. Statement of the energy lemma.
Let w in L^2(SM) have compact spatial support in the regular part
of M. Suppose Xw and Yw are in L^2(SM) distributionally. Let
E subset M be measurable, and suppose Xw = 0 on M \ E (i.e. Xw
is supported, as a function on SM, in the base-lift of E). Then
for every integer j:
    ||grad w_j||_2^2  <=  4  ||Xw||_2  ||1_E grad w||_2,
where grad is the flat gradient on M (acting on the base variable;
||grad w_j||_2^2 = sum over j of (||d w_j||^2 + ||dbar w_j||^2)),
and 1_E grad w is the L^2 function equal to grad w on the base
set E and 0 elsewhere.

STEP 24. Proof of the energy lemma.
Two ingredients.

Ingredient 1 (equal Cauchy-Riemann norms). For each coefficient
w_j (compact spatial support),
    ||d w_j||_2^2  =  ||dbar w_j||_2^2  =  (1/4) ||grad w_j||_2^2.
Proof of the equal-norm identity: the difference
    ||d w_j||^2 - ||dbar w_j||^2
is a constant times int_{M} d(w_j) wedge d(conj(w_j)). The
primitive w_j * d(conj(w_j)) is a global 1-form because the
constant unitary transition factors of w_j and conj(w_j) cancel;
its support is compact in the regular part of M. Stokes gives
zero. (Spatially smooth w_j suffice for the computation; ordinary
local Sobolev approximation gives the stated generality.)
Define D_j := ||d w_j||_2^2 = (1/4) ||grad w_j||_2^2.

Ingredient 2 (telescoping). By STEP 15 with n = j+1,
    (Xw)_{j+1}  =  d w_j  +  dbar w_{j+2}.
A short algebraic identity gives
    D_j - D_{j+2}
      = Re <(Xw)_{j+1},  d w_j  -  dbar w_{j+2}>,
where <.,.> is complex-linear in the first entry.

Sum this in k = 0, 1, ..., N-1:
    D_j - D_{j+2N}  =  sum_{k=0}^{N-1} Re <(Xw)_{j+2k+1},
                                d w_{j+2k} - dbar w_{j+2k+2}>.
Since (Xw)_{j+2k+1} is supported in E, we can replace the second
factor by 1_E(d w_{j+2k} - dbar w_{j+2k+2}). Cauchy-Schwarz in the
k-index gives
    D_j  <=  ||Xw||_2 *
        (sum_k  ||1_E (d w_{j+2k} - dbar w_{j+2k+2})||_2^2)^{1/2}
    + D_{j+2N}.
The sum on the right converges because sum_k D_k = (1/4) ||
grad w ||_2^2 < infinity (Parseval), hence D_{j+2N} -> 0. On the
other hand, pointwise in the base variable,
    |d w_{j+2k} - dbar w_{j+2k+2}|^2  <=
        2 (|d w_{j+2k}|^2 + |dbar w_{j+2k+2}|^2),
and summing over k, pointwise on SM, is bounded by
    2 sum_l (|d w_l|^2 + |dbar w_l|^2)  =  |grad w|^2.
Parseval plus Ingredient 1 gives
    sum_k ||1_E (d w_{j+2k} - dbar w_{j+2k+2})||_2^2  <=  ||1_E grad w ||_2^2.
Substituting and using D_j = (1/4) ||grad w_j||_2^2 gives
    ||grad w_j||_2^2  <=  4 ||Xw||_2 ||1_E grad w||_2.
This is the energy lemma.

STEP 25. Bounded-modes lemma (statement + proof).
Statement. Let h in L^inf(SM) with H = ||h||_inf. Suppose
Xh = 0 distributionally on the regular part of M. Then every
Fourier coefficient h_j has a global L^2 gradient, and
    ||grad h_j||_2  <=  C H,   uniformly in j in Z.
The constant C depends only on M and on the smoothing kernel rho.

Proof. Fix a smooth radial cutoff chi_eps on M with
    chi_eps = 0 on { r <= (K+1) eps },
    chi_eps = 1 on { r >= (K+2) eps },
    |grad chi_eps| <= C / eps.
For small eps, the transition regions are disjoint cone annuli.
Let w := chi_eps S_eps h, extended by 0 near the vertices. Then
w has compact spatial support in the regular part of M.
Distributional commutation with X gives
    X w  =  (X chi_eps) S_eps h,
since X(S_eps h) = S_eps(Xh) = 0 on the interior of the safe
region. So Xw is supported in the transition annuli
E_eps := { (K+1) eps < r < (K+2) eps }, whose total area is
O(eps^2). On E_eps, using |S_eps h| <= H and |X chi_eps| <= C/eps:
    |X w|  <=  C H / eps.
Similarly, |grad w| <= C H / eps on E_eps (gradient of chi_eps
bounds, gradient of S_eps h bounded by (a)).

Apply STEP 23 (energy lemma) to w with the base set E := E_eps.
On E_eps we have |X w| <= C H / eps and |grad w| <= C H / eps,
and Area(E_eps) = O(eps^2), so
    ||X w||_2 <= (C H / eps) sqrt(Area(E_eps)) = O(H),
    ||1_{E_eps} grad w||_2 <= (C H / eps) sqrt(Area(E_eps)) = O(H).
By STEP 23:
    ||grad w_j||_2^2  <=  4 ||X w||_2 ||1_{E_eps} grad w||_2
                      <=  C' H^2.
Uniform in eps and j. Hence ||grad w_j||_2 <= C H.

Now let eps -> 0 along any sequence. w -> h in L^2 on compact
subsets of the safe region (STEP 21c), so w_j -> h_j locally in
L^2. The gradients grad w_j are uniformly bounded in L^2, hence
have weakly convergent subsequences in L^2. Testing in regular
charts identifies every subsequential weak limit with grad h_j.
By lower semicontinuity of the L^2 norm under weak limits:
    ||grad h_j||_2  <=  liminf_{eps -> 0} ||grad w_j||_2  <=  C H.
This holds for every fixed j, independently of eps. Also, it
gives ||grad h_j||_2 < infinity for every j in Z, so h_j has a
global L^2 gradient. QED.

============================================================
PART 4. RIGIDITY: THE TWO-LIMIT ARGUMENT
============================================================

Fix f in L^inf(SM) with Xf = 0 distributionally. Let H := ||f||_inf.
We keep m in Z fixed (STEP 18) and prove  = 0.
STEP 19 then closes the rigidity proposition.

The proof runs two limits, deeply nested.
  Limit 1 (eps -> 0, s fixed): kill the smoothing error, producing
      an exactly invariant g^(s), with ||g^(s) - f||_2 <= C sqrt(s).
  Limit 2 (s -> 0): upgrade "exactly invariant" to "constant".
The first limit uses the energy lemma. The second uses the
bounded-modes lemma on the exact invariant g^(s). The two limits
have different eps/s scalings and cannot be combined.

STEP 26. Long-segment average g_eps.
Fix s in (0, 1), to be sent to 0 later. Put T := s / eps. Fix
L > K+4 (STEP 21d). Define the good set:
    G_eps := { z in SM : r(Phi_t z) > L eps for all 0 <= t <= T }.
Let
    a_eps := f * 1_{G_eps},
    g_eps := (1/T) int_0^T a_eps o Phi_{-t} dt.
g_eps is a long backward time-average of f, only on trajectories
that stay clear of the cones. Note: g_eps vanishes on { r < L eps }
-- if g_eps(z) != 0, then some Phi_{-t} z is in G_eps, hence
r(z) = r(Phi_t Phi_{-t} z) > L eps.

STEP 27. Quantitative bounds on g_eps.
  (a) ||g_eps||_inf <= H.
  (b) ||g_eps - f||_2 <= C H sqrt(s + eps^2).
  (c) X g_eps = (a_eps - a_eps o Phi_{-T}) / T weakly, so
      ||X g_eps||_2 <= 2 H / T = 2 H eps / s.

Proof of (b). Flow invariance: f o Phi_{-t} = f a.e. So
    g_eps - f = (1/T) int_0^T (a_eps - f) o Phi_{-t} dt.
The integrand a_eps - f is supported on SM minus G_eps and bounded
by 2H in absolute value. Minkowski in L^2:
    ||g_eps - f||_2 <= (1/T) int_0^T ||(a_eps - f) o Phi_{-t}||_2 dt
                    <= 2H * sqrt(nu(SM minus G_eps)),
where nu is the normalized Liouville measure (Phi_t preserves nu).
Now estimate nu(SM minus G_eps). A state is in SM minus G_eps iff
its trajectory leaves { r > L eps } at some t in [0, T]. Since r
is 1-Lipschitz along the unit-speed trajectory, this is detected
by sampling at times 0, eps, 2 eps, ..., T+eps -- at most
T/eps + 2 samples, each contributing at most
nu({ r <= (L+1) eps }) = O(eps^2) to the preimage measure. Hence
    nu(SM minus G_eps) <= (T/eps + 2) O(eps^2) = O(T eps + eps^2)
                      = O(s + eps^2).
Substituting gives (b).

Proof of (c). Differentiate g_eps under the flow. For each t the
function t |-> a_eps o Phi_{-t} is differentiable a.e., so
    d/dt (a_eps o Phi_{-t}) = - X(a_eps o Phi_{-t}).
Integrating the average in t, and using the fundamental theorem of
calculus, gives the weak identity (c). It does not differentiate
the step function 1_{G_eps}. Both |a_eps| and |a_eps o Phi_{-T}|
are bounded by H, so
    ||X g_eps||_2 <= (H + H) / T = 2H/T.

STEP 28. The smoothed equation.
On the safe region, X and Y commute and S_eps commutes with X, so
    X (Y S_eps U) = Y (X S_eps U) = Y S_eps (XU) = Y S_eps B,
where U and B are the endpoint objects from STEP 18. This is an
exact identity on the safe region. Both sides are spatially smooth
away from the cone points.

STEP 29. The segment identity.
Pair STEP 28 against g_eps and integrate by parts. First:
     = 
                       = -
                       = -(1/T) .
The last pairing expands (using the definition (c) of X g_eps) to
    = (1/T) int_{G_eps} f(z) ( (Y S_eps U)(Phi_T z)
                              - (Y S_eps U)(z) ) dnu(z).
(Change variables z -> Phi_T z in the second term, using flow
invariance of nu and Phi_{-T} Phi_T = identity; and f o Phi_{-T}
= f on G_eps.)

Exception handling. For each fixed eps, cover the safe region by
finitely many flat charts big enough for the S_eps convolution.
On each chart, for a.e. angular slice, the smoothed fields have
smooth spatial representatives and the differential identity of
STEP 28 holds there. During a chart visit a regular orbit stays at
constant chart angle, so an orbit hits an exceptional slice only if
it stays there over an open time interval, which contains a rational
time. For each chart and each rational time, the set of initial
states whose orbit visits the exceptional slice at that time has
measure zero by measure preservation. Discard this countable union
of null sets; the calculus above now runs along almost every good
segment. This bookkeeping does not affect any of the quantitative
bounds below.

STEP 30. The RHS of the segment identity is small.
Bound the integrand pointwise by H times the absolute difference
of two evaluations of Y S_eps U. Change variables on the second
evaluation using measure preservation of Phi_T:
    |RHS| <= (2H/T) ||Y S_eps U||_{L^2(r > K eps)}.
The base points z and Phi_T z both sit in the safe region on the
support of the integrand (z in G_eps implies r(Phi_s z) > L eps for
all s in [0, T]). With T = s / eps,
    |RHS| <= (2H/s) * eps * ||Y S_eps U||_{L^2(r > K eps)}.
By STEP 21(b), eps * ||Y S_eps U||_{L^2(r > K eps)} -> 0 as eps -> 0,
at fixed s. Hence RHS = o(1).

STEP 31. Adjoint:  = o(1).
Set w_eps := S_eps g_eps. By STEP 21(d), w_eps is supported in
{r >= (L-1) eps}, compactly supported in M. Applying the adjoint
identity (STEP 21d):
     = -.
The right side is exactly minus the segment identity from STEP 29,
which equals minus RHS. So
    || = o(1)  as eps -> 0, s fixed.

STEP 32. Gradient bounds on w_eps.
By STEP 27(a) and STEP 21(a):
    ||w_eps||_inf <= H,
    ||grad w_eps||_2 <= C H / eps.
By STEP 27(c) and STEP 21(e):
    ||X w_eps||_2 <= C ||X g_eps||_2 <= C H eps / s.
Apply the energy lemma (STEP 23) to w_eps with E = M (full space;
X w_eps has the same spatial support as w_eps). For every j in Z:
    ||grad (w_eps)_j||_2^2 <= 4 ||X w_eps||_2 ||grad w_eps||_2
                           <= 4 (C H eps / s)(C H / eps)
                           = C H^2 / s.
So
    ||grad (w_eps)_j||_2 <= C H / sqrt(s),
uniform in eps and j.

STEP 33. The first limit: eps -> 0, s fixed.
Fix s > 0. Along a sequence eps -> 0, w_eps has an L^2 bound and
individual coefficient gradients uniformly bounded in L^2 (STEP
32). By weak compactness, a subsequence converges weakly in L^2 to
a limit g^(s). The following hold.
  (a) ||g^(s)||_inf <= H. The closed convex set { |u| <= H } of
      L^2 is weakly closed.
  (b) X g^(s) = 0 distributionally. From STEP 21(e),
      X w_eps = S_eps X g_eps, and ||X g_eps||_2 <= 2H eps / s -> 0,
      so X w_eps -> 0 in L^2; the weak limit satisfies X g^(s) = 0.
  (c) ||g^(s) - f||_2 <= C H sqrt(s). By STEP 27(b), g_eps is
      close to f, and by STEP 21(c) plus symmetry of convolution,
      w_eps is close to g_eps. Their common weak L^2 limit g^(s)
      therefore satisfies ||g^(s) - f||_2 <= liminf ||g_eps - f||_2
      <= C H sqrt(s), using weak lower semicontinuity.
  (d) The two coefficient gradients that matter. B lives in mode
      m-1; only modes m-2 and m appear in the pairing .
      On a subsequence, d(w_eps)_{m-2} and dbar(w_eps)_m converge
      weakly in L^2 (uniform bound from STEP 32). Testing against
      smooth compactly supported fields identifies the limits with
      d(g^(s))_{m-2} and dbar(g^(s))_m respectively.
  (e) Since B is fixed in L^2, the weak limit of 
      is . Combined with STEP 31:
          = 0,  for every s > 0.

STEP 34. The second limit: s -> 0.
We cannot use the bound ||grad (w_eps)_j|| <= C H / sqrt(s) here
because s -> 0. Instead, apply the bounded-modes lemma (STEP 25)
directly to the exact invariant g^(s):
  - g^(s) in L^inf, ||g^(s)||_inf <= H;
  - X g^(s) = 0 distributionally (STEP 33b).
STEP 25 gives
    ||grad (g^(s))_j||_2 <= C H,  uniformly in s and j.
No 1/sqrt(s) blowup. In parallel:
    g^(s) -> f strongly in L^2  as s -> 0  (STEP 33c).

Choose a sequence s_n -> 0 along which the two relevant coefficient
gradients d(g^(s_n))_{m-2} and dbar(g^(s_n))_m have weak L^2
limits. The uniform bound says these limits exist as L^2 elements.
Because g^(s_n) -> f strongly in L^2, testing in regular charts
identifies d(g^(s_n))_{m-2} -> d f_{m-2};
dbar(g^(s_n))_m -> dbar f_m. So d f_{m-2} and dbar f_m exist as
L^2 elements and equal the distributional gradients.

Passing to the limit in  = 0 (Step 34e):
     = 0.

STEP 35. Close the argument.
By STEP 19, P_{m-1} Y f = -2 i B, so
    0 =  =  = -2 i ||B||_2^2.
Hence ||B||_2 = 0, so B = 0, so dbar f_m = 0. Since m in Z was
arbitrary, dbar f_k = 0 for every k. The recurrence (STEP 16) then
gives d f_k = -dbar f_{k+2} = 0 for every k. Each Fourier
coefficient f_k is locally constant on M; since M is connected,
f_k is constant on M. This is the rigidity proposition.

Why two limits. The first limit uses the energy lemma to get
||grad (w_eps)_j|| <= C H / sqrt(s), uniform in eps at fixed s.
The second limit avoids this bound (it blows up as s -> 0) by
relying on the exact invariance of g^(s), which activates the
bounded-modes lemma and gives the uniform-in-s bound C H. Reverse
the order and nothing works: without spatial smoothing, the first
limit (eps -> 0) is not justified; without exact invariance, the
second limit (s -> 0) has no uniform-in-s gradient bound. Both
limits are essential, and they must run in this order.

============================================================
PART 5. HOLONOMY AND CONCLUSION
============================================================

STEP 36. Each f_j is a parallel section of its line bundle.
By STEP 35, d f_j = dbar f_j = 0 as distributions on M. In a
regular chart, a function with both d and dbar zero is locally
constant. On the whole surface M with the transition rules from
STEP 14, d f_j = dbar f_j = 0 makes f_j a covariantly constant
("parallel") section of the flat line bundle L^{-j}: locally it is
a constant, and these local constants glue together through the
transition phases e^{-ij phi}.

"Parallel" is not "globally constant" once j != 0. The difference
is exactly a holonomy phase, measured by going around a cone point.

STEP 37. Holonomy around the alpha-cone.
Fix the vertex of Q with angle alpha (the alpha-cone). The cone on
M there has total angle 2 alpha. Cut out a small disk around the
cone point and unroll it flat: it becomes a sector with total angle
2 alpha. Going once around, the flat frame rotates by 2 alpha.
Concretely, going around the cone point, the angular coordinate on
a chart changes by 2 alpha. The Fourier coefficient f_j, a section
of L^{-j}, therefore transforms (STEP 14) by
    f_j  ->  e^{-2 i j alpha} f_j.
For f_j to be a single-valued section on M, the phase must be one:
    e^{-2 i j alpha} = 1.

STEP 38. Irrationality kills every nonzero mode.
Assume (WLOG after relabeling) alpha/pi is irrational. The condition
e^{-2 i j alpha} = 1 is equivalent to 2 j alpha = 2 pi k for some
integer k, i.e. alpha/pi = k/j. Since alpha/pi is irrational, no
nonzero integer j works. Hence f_j = 0 for every j != 0.

The zeroth coefficient f_0 is an ordinary function on M (STEP 14)
with d f_0 = dbar f_0 = 0. Since M is connected, f_0 is a constant
c on M. Therefore f = f_0 = c almost everywhere.

STEP 39. The theorem, complete.
Let A be a measurable flow-invariant set of states. Let f = 1_A.
Then f is bounded measurable and Xf = 0 distributionally. By
STEPS 36-38, f is constant a.e., hence f = 0 a.e. or f = 1 a.e.
Equivalently, mu(A) is 0 or 1. Ergodicity holds for every triangle
with at least one irrational angle/pi. The proof of STEP 3 is
complete.

============================================================
PART 6. CONSEQUENCES AND REMARKS
============================================================

STEP 40. Quantum ergodicity: interior.
Fix Q as in STEP 3, with Dirichlet or Neumann boundary condition.
Let (u_j) be any orthonormal eigenbasis of L^2(Q) for the
corresponding Laplacian, ordered by frequency:
    -Delta u_j = lambda_j^2 u_j,  0 <= lambda_1 <= lambda_2 <= ...
counted with multiplicity. Then there is a density-one subsequence
S_int of positive integers such that, for every classical order-zero
pseudodifferential operator A with compactly supported Schwartz
kernel in Q x Q,
    lim_{j -> infinity, j in S_int} _{L^2(Q)}
        = int_{Q x S^1} sigma_0(A) dmu.
Here "sigma_0(A)" is the principal symbol restricted to the unit
cosphere bundle, and "density one" means #(S_int intersect
{1,...,N}) / N -> 1 as N -> infinity.

Interpretation: eigenfunctions equidistribute in phase space
along a density-one subsequence. Not every j, but a typical one.
The result is due to Zelditch-Zworski 1996 (Dirichlet form from
the control theory approach), extended to the boundary in Burq 2003.

STEP 41. Quantum ergodicity: boundary.
Let Gamma be the union of the three open sides of Q. For
lambda_j > 0, set h_j := 1/lambda_j and define
    u_j^b := h_j * d_v u_j on Gamma,  (Dirichlet)
    u_j^b := u_j on Gamma             (Neumann),
where v is the outward unit normal. Let B*Gamma := {(s, eta) in
T*Gamma : |eta| < 1}, using arc length s and signed tangential
covector eta. Then there is a density-one subsequence S_B of
positive integers such that, for every a in C_c^infty(T*Gamma) and
a properly supported order-zero semiclassical quantization
A_h = Op_h(a),
    lim_{j -> infinity, j in S_B} _{L^2(Gamma)}
        = (2 / (pi Area(Q))) int_{B*Gamma} a(s, eta) w_B(eta)
                                    ds deta,
with weights
    w_D(eta) = sqrt(1 - eta^2),   w_N(eta) = (1 - eta^2)^{-1/2}.
The Cauchy data are NOT normalized to have unit L^2(Gamma) norm.

Reference: Hassell-Zelditch 2004.

STEP 42. Nodal domains grow.
Assume the eigenbasis is real-valued. Then there is a density-one
subsequence S_nod such that the number of nodal domains of u_j
(connected components of Q minus { u_j = 0 }) tends to infinity as
lambda_j -> infinity with j in S_nod.

Interpretation: the number of sign-change lobes grows without bound
along a typical subsequence. Yau's conjecture predicts a rate; our
theorem only says unbounded.

Reference: Hezari 2016.

STEP 43. What the theorem does NOT say.
- Not unique ergodicity. The theorem says the Liouville measure mu
  is ergodic. It does NOT say mu is the only invariant probability
  measure (rational polygons, e.g., have many).
- Not mixing. Mixing (decay of correlations) is strictly stronger.
  For a specific irrational triangle, mixing is open. Chaika-Forni
  2026 proved weak mixing for a dense G_delta of n-gons.
- Not about every orbit. Individual orbits may be periodic, or
  dense, or anything; ergodicity is statistical.
- No quantitative rate. Time averages converge, but slowly or
  quickly is not specified.
- Density-one subsequences are not canonical. The specific sets
  S_int, S_B, S_nod depend on Q, on the boundary condition, and on
  the eigenbasis.
- The nodal count has no rate. Only "tends to infinity", not a
  growth rate in lambda_j.
- Not incompatible with scarring. Individual eigenfunctions might
  concentrate on special sets; the theorem only rules out
  positive-density scarring.

STEP 44. Notation glossary.
Geometry:
  alpha, beta, gamma: triangle angles, alpha + beta + gamma = pi.
  pi: 3.14159...
  S^1: unit circle, identified with [0, 2 pi).
  Q: interior of the triangle.
  M: the double of Q (compact flat surface with three cone points).
  r(x): distance from x in M to the nearest cone point.
  K: constant making embedded eps-disks available when r > K eps.
  eps: smoothing scale, tends to 0.
  rho, rho_eps: radial smoothing bump.
  S_eps: spatial smoothing operator (STEP 20).
Phase space:
  SM: M x S^1, the unit tangent bundle (= phase space).
  (p, theta): state (position in M, direction in S^1).
  Phi_t: geodesic (billiard) flow.
  nu: normalized Liouville measure dA dtheta on SM.
  X, Y: horizontal vector fields on SM (flow direction, and
        its 90-degree rotation). [X, Y] = 0 in a flat chart.
  d, dbar: complex derivatives (1/2)(d/dx -+ i d/dy).
  grad: flat horizontal gradient on M.
  H: ||f||_inf.
Fourier machinery:
  f_j: j-th angular Fourier coefficient of f on SM.
  P_j: projection onto the j-th Fourier mode e^{ij theta}.
  U: sum_{k >= 0} P_{m+2k} f (the tail).
  B: e^{i(m-1)theta} dbar f_m (the endpoint).
  L^{-j}: flat unitary line bundle whose sections include f_j.
  m: fixed parity integer (the mode we are testing).
Rigidity machinery:
  s: small parameter, averaging length, tends to 0.
  T: averaging time, T = s / eps.
  L: constant with L > K+4 (segments must clear cones by L eps).
  G_eps: good set { r(Phi_t z) > L eps for all t in [0, T] }.
  a_eps: f * 1_{G_eps}.
  g_eps: (1/T) int_0^T a_eps o Phi_{-t} dt.
  w_eps: S_eps g_eps.
  g^(s): weak L^2 limit of w_eps as eps -> 0, at fixed s.
Spectral machinery:
  lambda_j: eigenvalue of -Delta on Q, listed with multiplicity.
  u_j: eigenfunction.
  Gamma: union of the three open sides of Q (regular boundary).
  B*Gamma: cosphere bundle {(s, eta): |eta| < 1}.
  w_D, w_N: boundary weights for Dirichlet and Neumann.

STEP 45. Review questions for a high-school reader.
1.  What is a billiard?
    A point particle at unit speed reflecting specularly inside a
    polygon; between hits it moves in a straight line.
2.  What does ergodic mean?
    Every measurable invariant set has measure 0 or 1. Equivalently,
    every bounded invariant function is constant a.e.
3.  What is the double of a triangle?
    Two copies of the triangle glued along corresponding sides,
    giving a compact flat surface with three cone points of cone
    angles 2 alpha, 2 beta, 2 gamma.
4.  Why do we need the double?
    Straight-line motion on the double is equivalent to billiard
    motion on Q (specular reflection = gluing across the seam).
5.  Which triangles does the theorem cover?
    Everything except the all-rational ones. That means: 2 or 3
    irrational angle/pi ratios; never exactly 1 (they sum to 1).
6.  What is X, geometrically?
    The vector field generating the geodesic flow. In a flat
    chart, X = cos(theta) d/dx + sin(theta) d/dy.
7.  What is Y, geometrically?
    The vector field obtained by rotating X by 90 degrees.
    In a flat chart, Y = -sin(theta) d/dx + cos(theta) d/dy.
    In a flat chart, [X, Y] = 0.
8.  What is the mode identity?
    (Xu)_n = d u_{n-1} + dbar u_{n+1};  (Yu)_n = i d u_{n-1}
    - i dbar u_{n+1}. X and Y shift the angular index by +/-1.
9.  Why does Xf = 0 imply d f_k = -dbar f_{k+2}?
    Read off the (k+1)-th Fourier coefficient of Xf = 0.
10. Why is Xf = 0 not enough?
    f is only in L^inf, so coefficient derivatives have no a priori
    regularity beyond L^2. Formal series satisfy the recurrence
    without forcing any derivative to vanish.
11. What is the tail U?
    U = sum_{k >= 0} P_{m+2k} f. Infinite half of the arithmetic
    progression of modes with parity m.
12. What is the endpoint B, and why does B = 0 matter?
    B = e^{i(m-1)theta} dbar f_m. Then XU = B, and B = 0 forces
    dbar f_m = 0, which (for arbitrary m) triggers the whole
    rigidity conclusion via STEP 19.
13. What does S_eps do?
    Averages a field over a disk of radius eps in position, keeping
    direction theta frozen. Only acts on the safe region r > K eps.
14. What does the energy lemma say?
    ||grad w_j||^2 <= 4 ||Xw|| * ||1_E grad w||. It bounds each
    coefficient's gradient by the size of Xw on a "cutoff set" E.
15. What does the bounded-modes lemma say?
    If Xh = 0 and ||h||_inf = H, then every ||grad h_j||_2 <= C H,
    uniformly in j. No L^2-type blowup.
16. Why are there two limits?
    First limit (eps -> 0, s fixed) uses the energy lemma; gives
    ||grad (w_eps)_j|| <= C H / sqrt(s) uniform in eps. Second
    limit (s -> 0) uses the bounded-modes lemma on the exact
    invariant g^(s); gives ||grad (g^(s))_j|| <= C H uniform in s.
    The two limits have different eps/s scalings and cannot be
    combined.
17. Where does the irrationality hypothesis enter?
    STEP 37-38. Going around the alpha-cone changes theta by 2
    alpha, so a parallel section f_j of L^{-j} must satisfy
    e^{-2 i j alpha} = 1. If alpha/pi is irrational, no nonzero
    j works.
18. Does the theorem prove mixing? No. Weaker (and open for a
    specific irrational triangle).
19. Does the theorem say every orbit is dense? No; only statistical
    statements about typical starting states.
20. Does the theorem prove unique ergodicity? No. Only that the
    Liouville measure mu is ergodic (cannot be decomposed into two
    invariant pieces of positive measure). Not that it is the only
    invariant probability measure.

STEP 46. Compressed proof at a glance (for review).
Let Q be a nondegenerate triangle with at least one angle/pi
irrational, and let M be the double (STEP 7). On SM = M x S^1, let
X and Y be the horizontal vector fields (STEP 11): [X, Y] = 0.
Let f in L^inf(SM) with Xf = 0 distributionally. Goal: prove f is
constant a.e.

Fourier-expand f in theta (STEP 13): f = sum_j f_j(p) e^{ij theta},
and (Xu)_n = d u_{n-1} + dbar u_{n+1} (STEP 15). Then Xf = 0 gives
d f_{m-1} + dbar f_{m+1} = 0 for every m (STEP 16).
Fix m. Set U = sum_{k >= 0} P_{m+2k} f, B = e^{i(m-1)theta} dbar f_m.
Finite tails + weak limit give X U = B (STEP 18). And
P_{m-1} Y f = -2 i B (STEP 19). So if we prove  = 0,
then ||B||^2 = 0, hence B = 0 and dbar f_m = 0.

Show  = 0 by two limits.
Setup: fix s > 0, T = s/eps, L > K+4.
  G_eps := { r(Phi_t z) > L eps, 0 <= t <= T },
  a_eps := f 1_{G_eps},
  g_eps := (1/T) int_0^T a_eps o Phi_{-t} dt.
Bounds (STEP 27):
  ||g_eps||_inf <= H,
  ||g_eps - f||_2 <= C H sqrt(s + eps^2),
  ||X g_eps||_2 <= 2 H / T = 2 H eps / s,
  g_eps supported off { r < L eps }.
Smooth: w_eps := S_eps g_eps (STEP 20). Adjoint identity
 = - (STEP 21d). Small-scale derivative
eps ||Y S_eps U||_{L^2(r > K eps)} -> 0 (STEP 21b).
Integrate X(Y S_eps U) = Y S_eps B along good segments to get
   = o(1)  (STEPS 28-30),
hence (adjoint)  = o(1)  (STEP 31).
Energy lemma with E = M:
  ||grad (w_eps)_j||^2 <= 4 ||X w_eps|| ||grad w_eps||
                       <= C H^2 / s  (STEP 32).
So ||grad (w_eps)_j|| <= C H / sqrt(s), uniform in eps at fixed s.

First limit: eps -> 0 at fixed s. Subsequence w_eps -> g^(s) weakly
in L^2. Then
  X g^(s) = 0 distributionally,
  ||g^(s)||_inf <= H,
  ||g^(s) - f||_2 <= C H sqrt(s),
   = 0,  for every s > 0.  (STEP 33)

Second limit: s -> 0. Bounded-modes lemma on g^(s):
  ||grad (g^(s))_j|| <= C H, uniform in s and j.
And g^(s) -> f strongly in L^2 (from ||g^(s) - f||_2 <= C H sqrt(s)).
Taking a subsequence s_n -> 0 gives  = 0,
because the weak limits of the relevant coefficient gradients exist
(as L^2 elements) and equal d f_{m-2} and dbar f_m.  (STEP 34)

Close:
  0 =  =  = -2 i ||B||^2,
so B = 0, so dbar f_m = 0 (for arbitrary m), so d f_k = 0 for all k.
Every f_j is a parallel section of L^{-j} (STEP 36). Around the
alpha-cone: single-valuedness of f_j requires e^{-2 i j alpha} = 1.
If alpha/pi irrational, this forces f_j = 0 for every j != 0.
And f_0 is constant on the connected M. Hence f is constant.
Ergodic. (STEP 38-39)

For the spectral consequences (STEPS 40-42), just apply the
classical theorems: ZZ 1996 (interior QE), HZ 2004 (boundary QE),
Hezari 2016 (nodal count).