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.

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.

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,
    |<R_N, 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
    <B, P_{m-1} Y f>  =  <B, -2 i B>  =  -2 i ||B||_2^2.
If we can prove <B, P_{m-1} Y f> = 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 <B, P_{m-1} Y f> = 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),
        <Y S_eps u, g>  =  -  <u, Y S_eps g>.
  (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):
    <Y S_eps u, g>  =  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
                    =  -<u, Y S_eps g>.
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 minus 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 <B, P_{m-1} Y f> = 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:
    <Y S_eps B, g_eps> = <X(Y S_eps U), g_eps>
                       = -<Y S_eps U, X g_eps>
                       = -(1/T) <Y S_eps U, a_eps - a_eps o Phi_{-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: <B, Y w_eps> = 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):
    <B, Y w_eps> = -<Y S_eps B, g_eps>.
The right side is exactly minus the segment identity from STEP 29,
which equals minus RHS. So
    |<B, Y w_eps>| = 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 <B, Y ...>.
      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 <B, Y w_eps>
      is <B, P_{m-1} Y g^(s)>. Combined with STEP 31:
         <B, P_{m-1} Y g^(s)> = 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 <B, P_{m-1} Y g^(s_n)> = 0 (Step 34e):
    <B, P_{m-1} Y f> = 0.

STEP 35. Close the argument.
By STEP 19, P_{m-1} Y f = -2 i B, so
    0 = <B, P_{m-1} Y f> = <B, -2 i B> = -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.

============================================================