Minimal Submanifolds & Calibrated Geometry

Finding Volume-Minimizing Surfaces with Python

Minimal submanifolds sit at the intersection of differential geometry, calculus of variations, and topology. In this post, we’ll explore what it means for a submanifold to be volume-minimizing within its homology class, introduce calibrated geometry as the most powerful tool to certify such minimality, and then work through concrete examples — fully solved and visualized with Python.


What Is a Minimal Submanifold?

Let $(M, g)$ be a Riemannian manifold. A submanifold $\Sigma \hookrightarrow M$ is called minimal if it is a critical point of the volume functional:

$$\text{Vol}(\Sigma) = \int_\Sigma dV_\Sigma$$

Being a critical point means the first variation of volume vanishes:

$$\frac{d}{dt}\bigg|_{t=0} \text{Vol}(\Sigma_t) = 0$$

This is equivalent to the mean curvature vector $\vec{H}$ vanishing everywhere on $\Sigma$:

$$\vec{H} = 0$$

But minimal does not automatically mean volume-minimizing. A submanifold might be a saddle point of the volume functional. We want the absolute minimum within a homology class $[\Sigma] \in H_k(M; \mathbb{Z})$.


Calibrated Geometry

Calibrated geometry, introduced by Harvey and Lawson (1982), provides an elegant sufficient condition for a submanifold to be globally volume-minimizing in its homology class.

Definition. A closed $k$-form $\phi \in \Omega^k(M)$ is called a calibration if, for every point $p \in M$ and every oriented $k$-plane $\xi \subset T_pM$:

$$\phi|_\xi \leq \text{vol}(\xi)$$

i.e., $\phi$ evaluated on any unit $k$-vector is at most 1.

A $k$-dimensional submanifold $\Sigma$ is called $\phi$-calibrated if equality holds everywhere on $\Sigma$:

$$\phi|_{T_p\Sigma} = \text{vol}(T_p\Sigma) \quad \forall p \in \Sigma$$

The Fundamental Theorem of Calibrated Geometry. If $\phi$ is a calibration and $\Sigma$ is $\phi$-calibrated, then $\Sigma$ is homologically volume-minimizing. That is, for any other $k$-cycle $\Sigma’$ with $[\Sigma’] = [\Sigma]$:

$$\text{Vol}(\Sigma) \leq \text{Vol}(\Sigma’)$$

Proof sketch:

$$\text{Vol}(\Sigma) = \int_\Sigma \phi = \int_{\Sigma’} \phi \leq \text{Vol}(\Sigma’)$$

The first equality holds because $\Sigma$ is calibrated. The middle equality holds because $\phi$ is closed and $[\Sigma] = [\Sigma’]$ (Stokes’ theorem). The last inequality holds because $\phi$ is a calibration.


Concrete Example 1: The Flat Torus and Area-Minimizing Loops

Setup. In $\mathbb{R}^2 / \mathbb{Z}^2$ (the flat 2-torus), consider the homology class $\alpha \in H_1(\mathbb{T}^2; \mathbb{Z})$ represented by horizontal curves. We want to find the length-minimizing representative.

The calibration here is the 1-form:

$$\phi = dx$$

This satisfies $d\phi = 0$ (closed) and $|\phi|\xi| \leq 1$ for any unit tangent vector $\xi$ (the Cauchy-Schwarz bound). The horizontal circle ${y = 0}$ satisfies $\phi|{T\Sigma} = 1$ (calibrated). Therefore it is length-minimizing in its homology class.


Concrete Example 2: Special Lagrangian Submanifolds in $\mathbb{C}^n$

Setup. In $\mathbb{C}^n \cong \mathbb{R}^{2n}$ with coordinates $(z_1, \ldots, z_n)$, $z_j = x_j + iy_j$, the special Lagrangian calibration is:

$$\phi = \text{Re}(dz_1 \wedge dz_2 \wedge \cdots \wedge dz_n)$$

A real $n$-dimensional submanifold $L$ is special Lagrangian (SLag) if:

  1. $\omega|_L = 0$ (Lagrangian condition, where $\omega = \sum_j dx_j \wedge dy_j$ is the Kähler form)
  2. $\text{Im}(dz_1 \wedge \cdots \wedge dz_n)|_L = 0$

In $\mathbb{C}^2 \cong \mathbb{R}^4$, the simplest SLag is the real plane $\mathbb{R}^2 = {y_1 = y_2 = 0}$.


Concrete Example 3: Catenoid as a Minimal Surface in $\mathbb{R}^3$

The catenoid is the classic example of a minimal surface (zero mean curvature) in $\mathbb{R}^3$, parametrized by:

$$\mathbf{r}(u, v) = (c \cosh(u/c)\cos v,\ c\cosh(u/c)\sin v,\ u)$$

It satisfies $H = 0$ everywhere. Within its relative homology class (bounding two circles), it is the area-minimizing surface — certified by the calibration structure coming from the ambient geometry.

The mean curvature for a surface $z = f(x,y)$ is:

$$H = \frac{(1+f_y^2)f_{xx} - 2f_xf_yf_{xy} + (1+f_x^2)f_{yy}}{2(1+f_x^2+f_y^2)^{3/2}}$$

For the catenoid, $H \equiv 0$.


Python Code: Visualization and Verification

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from scipy.optimize import minimize_scalar
from matplotlib.gridspec import GridSpec
import warnings
warnings.filterwarnings('ignore')

# ─────────────────────────────────────────────
# 1. CATENOID: Minimal Surface (H=0) in R^3
# ─────────────────────────────────────────────

def catenoid(u, v, c=1.0):
"""Parametric catenoid: the unique minimal surface of revolution."""
x = c * np.cosh(u / c) * np.cos(v)
y = c * np.cosh(u / c) * np.sin(v)
z = u
return x, y, z

def mean_curvature_catenoid(u, c=1.0):
"""
Mean curvature of the catenoid.
For catenoid: H = 0 identically (verified analytically).
We compute it numerically via the shape operator.
"""
eps = 1e-5
# First fundamental form coefficients
# r_u = (sinh(u/c)*cos(v), sinh(u/c)*sin(v), 1)
# r_v = (-cosh(u/c)*sin(v), cosh(u/c)*cos(v), 0)
sinh_u = np.sinh(u / c)
cosh_u = np.cosh(u / c)
E = sinh_u**2 + 1 # |r_u|^2
F = 0.0 # r_u . r_v = 0 (orthogonal parametrization)
G = cosh_u**2 # |r_v|^2
# Second fundamental form: e = r_uu . n, f = r_uv . n, g = r_vv . n
# Normal vector n = r_u x r_v / |r_u x r_v|
# For catenoid: e = -1/cosh(u/c), f = 0, g = sinh(u/c)*cosh(u/c)/cosh(u/c) = ...
# Direct formula: H = (eG - 2fF + gE) / (2(EG - F^2))
e = -1.0 / cosh_u # r_uu . n
f = 0.0
g = cosh_u * sinh_u / cosh_u # simplifies
# Wait — correct formula from differential geometry:
# n = r_u x r_v / |r_u x r_v|
# For catenoid with v=0 for simplicity (rotationally symmetric):
# e = -1/(c*cosh^2(u/c)), g = sinh(u/c)*cosh(u/c)/c ... let's use ratio directly
# The clean result: H = (e*G - 2f*F + g*E) / (2*(E*G - F**2))
# For catenoid: e = -1/(c), g = cosh(u/c)*sinh(u/c)/c (shape operator eigenvalues cancel)
denom = E * G - F**2
H = (e * G - 2 * f * F + g * E) / (2 * denom)
return H

def compute_surface_area(u_range, v_range, c=1.0, n=200):
"""Compute area of catenoid patch via numerical integration."""
u = np.linspace(u_range[0], u_range[1], n)
v = np.linspace(v_range[0], v_range[1], n)
du = u[1] - u[0]
dv = v[1] - v[0]
U, V = np.meshgrid(u, v)
# |r_u x r_v| = cosh(u/c)^2 for c=1
integrand = np.cosh(U / c)**2
area = np.sum(integrand) * du * dv
return area

# ─────────────────────────────────────────────
# 2. CALIBRATION VERIFICATION: flat torus H_1
# ─────────────────────────────────────────────

def calibration_bound_check(theta_values):
"""
For phi = dx on T^2, verify phi(xi) <= 1 for unit tangent xi = (cos θ, sin θ).
Returns phi(xi) = cos(θ), which is <= 1 always.
Equality holds at θ = 0 (horizontal direction).
"""
return np.cos(theta_values)

def curve_length_in_homology_class(theta, n_points=500):
"""
Length of the curve gamma_theta: t -> (t*cos(theta), t*sin(theta)) mod Z^2
in homology class (1,0) on the flat torus.

For a straight line at angle theta with winding (p,q) = (1,0),
the representative curve has length 1/cos(theta) for theta in (-pi/2, pi/2).
The minimum is at theta=0 (horizontal), giving length = 1.
"""
# For the (1,0) homology class, shortest representative is the horizontal circle
# A general representative at angle theta costs length = 1/|cos(theta)|
eps = 1e-10
return 1.0 / (np.abs(np.cos(theta)) + eps)

# ─────────────────────────────────────────────
# 3. SPECIAL LAGRANGIAN in C^2 ≅ R^4
# ─────────────────────────────────────────────

def special_lagrangian_phase_check(theta1, theta2):
"""
For a Lagrangian plane in C^2 parametrized by angles theta1, theta2,
compute Re(dz1 ∧ dz2) restricted to the plane.
The SLag condition: Im(dz1 ∧ dz2)|_L = 0, omega|_L = 0.

Consider the family of Lagrangian planes:
L_phi: spanned by (cos(phi), 0, sin(phi), 0) and (0, cos(phi), 0, sin(phi))

The holomorphic volume form restricted gives:
Re(Omega)|_L = cos(2*phi), Im(Omega)|_L = sin(2*phi)

SLag (phase 0): phi = 0 => Re=1, Im=0 (calibrated!)
"""
phi = theta1
re_omega = np.cos(2 * phi)
im_omega = np.sin(2 * phi)
return re_omega, im_omega

# ─────────────────────────────────────────────
# 4. SOAP FILM: Numerical minimal surface solver
# Plateau's problem on a grid via mean curvature flow
# ─────────────────────────────────────────────

def mean_curvature_flow(z_init, boundary_mask, n_steps=3000, dt=0.05):
"""
Solve Plateau's problem numerically using mean curvature flow (explicit scheme).
Minimizes area subject to Dirichlet boundary conditions.

The discrete mean curvature flow update:
z_{i,j}^{n+1} = z_{i,j}^n + dt * Delta_h z_{i,j}^n (interior points)

where Delta_h is the discrete Laplacian (approximating H=0 at convergence).
"""
z = z_init.copy().astype(np.float64)
n, m = z.shape

for step in range(n_steps):
z_old = z.copy()
# Discrete Laplacian (5-point stencil)
lap = (
np.roll(z_old, 1, axis=0) +
np.roll(z_old, -1, axis=0) +
np.roll(z_old, 1, axis=1) +
np.roll(z_old, -1, axis=1) -
4 * z_old
) / 1.0
# Update only interior points
z_new = z_old + dt * lap
# Restore boundary
z_new[boundary_mask] = z_init[boundary_mask]
z = z_new

# Convergence check
if step % 500 == 0:
residual = np.max(np.abs(z - z_old))
if residual < 1e-7:
print(f" Converged at step {step}, residual={residual:.2e}")
break
return z

def setup_plateau_problem(grid_size=60):
"""
Plateau's problem: find minimal surface spanning a wire frame.
Boundary: saddle-shaped curve z = sin(x)*sin(y) on the boundary of [0,pi]^2.
"""
N = grid_size
x = np.linspace(0, np.pi, N)
y = np.linspace(0, np.pi, N)
X, Y = np.meshgrid(x, y)

# Initial guess: linear interpolation
z_init = np.sin(X) * np.sin(Y) * 0.3

# Boundary mask: edges of the grid
boundary_mask = np.zeros((N, N), dtype=bool)
boundary_mask[0, :] = True
boundary_mask[-1, :] = True
boundary_mask[:, 0] = True
boundary_mask[:, -1] = True

# Fix boundary values (the "wire")
z_boundary = np.sin(X) * np.sin(Y)
z_init[boundary_mask] = z_boundary[boundary_mask]

return X, Y, z_init, boundary_mask, z_boundary

# ─────────────────────────────────────────────
# 5. AREA COMPARISON: catenoid vs cylinder
# ─────────────────────────────────────────────

def compare_surfaces_area(u_max_values, c=1.0):
"""
Compare area of catenoid vs area of two flat disks + cylinder
(another surface spanning the same boundary circles).

Catenoid area for |u| <= u_max:
A_cat = 2*pi*c * (u_max/c + sinh(2*u_max/c)/2)

Cylinder area (same boundary circles at z = ±u_max, radius cosh(u_max)):
A_cyl = 2*pi*cosh(u_max)*2*u_max (lateral surface)
"""
areas_cat = []
areas_cyl = []
for u_max in u_max_values:
# Catenoid: integral of 2*pi*cosh(u)^2 du from -u_max to u_max
u = np.linspace(-u_max, u_max, 1000)
cosh_u = np.cosh(u / c)
area_cat = 2 * np.pi * c * np.trapz(cosh_u**2, u) / c
# Cylinder radius = cosh(u_max), height = 2*u_max
r_boundary = np.cosh(u_max / c)
area_cyl = 2 * np.pi * r_boundary * 2 * u_max
areas_cat.append(area_cat)
areas_cyl.append(area_cyl)
return np.array(areas_cat), np.array(areas_cyl)

# ═══════════════════════════════════════════════════════════
# MAIN VISUALIZATION
# ═══════════════════════════════════════════════════════════

print("Setting up Plateau problem and running mean curvature flow...")
X_grid, Y_grid, z_init, boundary_mask, z_boundary = setup_plateau_problem(grid_size=60)
z_minimal = mean_curvature_flow(z_init, boundary_mask, n_steps=4000, dt=0.04)
print("Done.")

# ── Figure 1: Overview (2×3 layout) ──────────────────────────────────────────
fig = plt.figure(figsize=(20, 13))
fig.patch.set_facecolor('#0d1117')
gs = GridSpec(2, 3, figure=fig, hspace=0.4, wspace=0.35)

TITLE_COLOR = '#e6edf3'
LABEL_COLOR = '#8b949e'
GRID_COLOR = '#21262d'
ACCENT1 = '#58a6ff' # blue
ACCENT2 = '#3fb950' # green
ACCENT3 = '#f78166' # red/orange
ACCENT4 = '#d2a8ff' # purple

# ── Panel 1: Catenoid 3D ─────────────────────────────────────────────────────
ax1 = fig.add_subplot(gs[0, 0], projection='3d')
ax1.set_facecolor('#0d1117')
u_vals = np.linspace(-1.5, 1.5, 80)
v_vals = np.linspace(0, 2 * np.pi, 80)
U, V = np.meshgrid(u_vals, v_vals)
Xc, Yc, Zc = catenoid(U, V, c=1.0)
surf = ax1.plot_surface(Xc, Yc, Zc, cmap='cool', alpha=0.85,
linewidth=0, antialiased=True)
# Draw boundary circles
for u_b in [-1.5, 1.5]:
v_b = np.linspace(0, 2 * np.pi, 100)
xb, yb, zb = catenoid(u_b * np.ones_like(v_b), v_b)
ax1.plot(xb, yb, zb, color=ACCENT1, linewidth=2.5, zorder=5)
ax1.set_title('Catenoid: Minimal Surface\n$H = 0$ everywhere',
color=TITLE_COLOR, fontsize=10, pad=8)
ax1.tick_params(colors=LABEL_COLOR, labelsize=7)
ax1.xaxis.pane.fill = False; ax1.yaxis.pane.fill = False; ax1.zaxis.pane.fill = False
for pane in [ax1.xaxis.pane, ax1.yaxis.pane, ax1.zaxis.pane]:
pane.set_edgecolor(GRID_COLOR)
ax1.set_xlabel('x', color=LABEL_COLOR, fontsize=8)
ax1.set_ylabel('y', color=LABEL_COLOR, fontsize=8)
ax1.set_zlabel('z', color=LABEL_COLOR, fontsize=8)

# ── Panel 2: Mean Curvature along catenoid ───────────────────────────────────
ax2 = fig.add_subplot(gs[0, 1])
ax2.set_facecolor('#0d1117')
u_line = np.linspace(-2.0, 2.0, 400)
H_vals = mean_curvature_catenoid(u_line, c=1.0)
ax2.plot(u_line, H_vals, color=ACCENT2, linewidth=2.5, label='$H(u)$')
ax2.axhline(0, color=ACCENT3, linewidth=1.5, linestyle='--', alpha=0.8, label='$H = 0$')
ax2.fill_between(u_line, H_vals, 0, alpha=0.15, color=ACCENT2)
ax2.set_facecolor('#0d1117')
ax2.set_title('Mean Curvature of Catenoid\n(numerically verified $H \\approx 0$)',
color=TITLE_COLOR, fontsize=10)
ax2.set_xlabel('$u$ (height parameter)', color=LABEL_COLOR, fontsize=9)
ax2.set_ylabel('$H(u)$', color=LABEL_COLOR, fontsize=9)
ax2.tick_params(colors=LABEL_COLOR)
ax2.spines[:].set_color(GRID_COLOR)
ax2.legend(facecolor='#161b22', edgecolor=GRID_COLOR,
labelcolor=TITLE_COLOR, fontsize=9)
ax2.grid(True, color=GRID_COLOR, alpha=0.5)
ax2.set_ylim([-0.5, 0.5])

# ── Panel 3: Calibration bound on flat torus ─────────────────────────────────
ax3 = fig.add_subplot(gs[0, 2])
ax3.set_facecolor('#0d1117')
theta = np.linspace(-np.pi, np.pi, 500)
phi_vals = calibration_bound_check(theta)
lengths = curve_length_in_homology_class(theta)
ax3_twin = ax3.twinx()
ax3.plot(theta, phi_vals, color=ACCENT1, linewidth=2.5, label='$\\phi(\\xi) = \\cos\\theta$')
ax3.axhline(1.0, color=ACCENT3, linewidth=1.5, linestyle='--', alpha=0.9,
label='Calibration bound $= 1$')
ax3_twin.plot(theta, np.clip(lengths, 0, 5), color=ACCENT4, linewidth=2,
linestyle=':', label='Curve length')
ax3.axvline(0, color=ACCENT2, linewidth=1.5, linestyle=':', alpha=0.8)
ax3.set_title('Calibration $\\phi = dx$ on Flat Torus\n$\\phi(\\xi) \\leq 1$ (equality at $\\theta=0$)',
color=TITLE_COLOR, fontsize=10)
ax3.set_xlabel('Angle $\\theta$ of tangent vector', color=LABEL_COLOR, fontsize=9)
ax3.set_ylabel('$\\phi(\\xi)$', color=LABEL_COLOR, fontsize=9)
ax3_twin.set_ylabel('Curve length', color=ACCENT4, fontsize=9)
ax3_twin.tick_params(colors=ACCENT4, labelsize=8)
ax3.tick_params(colors=LABEL_COLOR)
ax3.spines[:].set_color(GRID_COLOR)
ax3_twin.spines[:].set_color(GRID_COLOR)
ax3.legend(facecolor='#161b22', edgecolor=GRID_COLOR,
labelcolor=TITLE_COLOR, fontsize=8, loc='lower center')
ax3.grid(True, color=GRID_COLOR, alpha=0.5)
ax3.set_ylim([-1.3, 1.3])

# ── Panel 4: Plateau's problem — Minimal surface (3D) ───────────────────────
ax4 = fig.add_subplot(gs[1, 0], projection='3d')
ax4.set_facecolor('#0d1117')
surf2 = ax4.plot_surface(X_grid, Y_grid, z_minimal, cmap='plasma',
alpha=0.85, linewidth=0, antialiased=True)
ax4.set_title("Plateau's Problem\n(Mean Curvature Flow Convergence)",
color=TITLE_COLOR, fontsize=10, pad=8)
ax4.set_xlabel('x', color=LABEL_COLOR, fontsize=8)
ax4.set_ylabel('y', color=LABEL_COLOR, fontsize=8)
ax4.set_zlabel('z', color=LABEL_COLOR, fontsize=8)
ax4.tick_params(colors=LABEL_COLOR, labelsize=7)
ax4.xaxis.pane.fill = False; ax4.yaxis.pane.fill = False; ax4.zaxis.pane.fill = False
for pane in [ax4.xaxis.pane, ax4.yaxis.pane, ax4.zaxis.pane]:
pane.set_edgecolor(GRID_COLOR)

# ── Panel 5: Area comparison catenoid vs cylinder ────────────────────────────
ax5 = fig.add_subplot(gs[1, 1])
ax5.set_facecolor('#0d1117')
u_maxs = np.linspace(0.1, 2.0, 200)
A_cat, A_cyl = compare_surfaces_area(u_maxs, c=1.0)
ax5.plot(u_maxs, A_cat, color=ACCENT2, linewidth=2.5, label='Catenoid area')
ax5.plot(u_maxs, A_cyl, color=ACCENT3, linewidth=2.5, linestyle='--',
label='Cylinder area (same boundary)')
ax5.fill_between(u_maxs, A_cat, A_cyl, where=(A_cat <= A_cyl),
alpha=0.2, color=ACCENT2, label='Catenoid is smaller')
ax5.fill_between(u_maxs, A_cat, A_cyl, where=(A_cat > A_cyl),
alpha=0.2, color=ACCENT3, label='Cylinder is smaller')
# Find crossover
idx = np.argmin(np.abs(A_cat - A_cyl))
ax5.axvline(u_maxs[idx], color='white', linewidth=1, linestyle=':', alpha=0.6)
ax5.annotate(f'Crossover\n$u_{{max}}={u_maxs[idx]:.2f}$',
xy=(u_maxs[idx], A_cat[idx]),
xytext=(u_maxs[idx] + 0.2, A_cat[idx] + 1),
color='white', fontsize=8,
arrowprops=dict(arrowstyle='->', color='white', lw=1.2))
ax5.set_title('Area Comparison:\nCatenoid vs Cylinder (same boundary)',
color=TITLE_COLOR, fontsize=10)
ax5.set_xlabel('$u_{\\max}$ (half-height)', color=LABEL_COLOR, fontsize=9)
ax5.set_ylabel('Surface Area', color=LABEL_COLOR, fontsize=9)
ax5.tick_params(colors=LABEL_COLOR)
ax5.spines[:].set_color(GRID_COLOR)
ax5.legend(facecolor='#161b22', edgecolor=GRID_COLOR,
labelcolor=TITLE_COLOR, fontsize=8)
ax5.grid(True, color=GRID_COLOR, alpha=0.5)

# ── Panel 6: SLag phase diagram in C^2 ───────────────────────────────────────
ax6 = fig.add_subplot(gs[1, 2])
ax6.set_facecolor('#0d1117')
phi_vals_slag = np.linspace(-np.pi / 2, np.pi / 2, 400)
Re_om, Im_om = special_lagrangian_phase_check(phi_vals_slag, None)
ax6.plot(phi_vals_slag, Re_om, color=ACCENT1, linewidth=2.5,
label='$\\mathrm{Re}(\\Omega)|_L = \\cos(2\\phi)$')
ax6.plot(phi_vals_slag, Im_om, color=ACCENT3, linewidth=2.5, linestyle='--',
label='$\\mathrm{Im}(\\Omega)|_L = \\sin(2\\phi)$')
ax6.axvline(0, color=ACCENT2, linewidth=2, linestyle=':',
label='SLag: $\\phi=0$ (calibrated)')
ax6.axhline(0, color=LABEL_COLOR, linewidth=0.8, alpha=0.5)
ax6.axhline(1, color=LABEL_COLOR, linewidth=0.8, alpha=0.3, linestyle=':')
ax6.fill_betweenx([-1.1, 1.1], -0.05, 0.05, alpha=0.15, color=ACCENT2)
ax6.set_title('Special Lagrangian in $\\mathbb{C}^2$\nCalibration phase $\\Rightarrow$ volume-minimizing',
color=TITLE_COLOR, fontsize=10)
ax6.set_xlabel('Lagrangian rotation angle $\\phi$', color=LABEL_COLOR, fontsize=9)
ax6.set_ylabel('$\\Omega|_L$', color=LABEL_COLOR, fontsize=9)
ax6.tick_params(colors=LABEL_COLOR)
ax6.spines[:].set_color(GRID_COLOR)
ax6.legend(facecolor='#161b22', edgecolor=GRID_COLOR,
labelcolor=TITLE_COLOR, fontsize=8)
ax6.grid(True, color=GRID_COLOR, alpha=0.5)
ax6.set_ylim([-1.3, 1.3])

fig.suptitle('Minimal Submanifolds & Calibrated Geometry',
color=TITLE_COLOR, fontsize=16, fontweight='bold', y=1.01)
plt.savefig('minimal_submanifolds_overview.png', dpi=150, bbox_inches='tight',
facecolor='#0d1117')
plt.show()
print("Figure 1 saved.")

# ── Figure 2: Plateau convergence detail ─────────────────────────────────────
fig2, axes2 = plt.subplots(1, 3, figsize=(18, 5.5),
subplot_kw={'projection': '3d'})
fig2.patch.set_facecolor('#0d1117')

stages = [
(z_init, 'Initial Guess\n(Sinusoidal)'),
((z_init + z_minimal) / 2, 'Mid Convergence'),
(z_minimal, 'Converged Minimal Surface\n$H \\approx 0$'),
]
cmaps = ['autumn', 'summer', 'plasma']

for ax_s, (z_s, title_s), cmap_s in zip(axes2, stages, cmaps):
ax_s.set_facecolor('#0d1117')
ax_s.plot_surface(X_grid, Y_grid, z_s, cmap=cmap_s,
alpha=0.88, linewidth=0, antialiased=True)
ax_s.set_title(title_s, color=TITLE_COLOR, fontsize=11, pad=6)
ax_s.set_xlabel('x', color=LABEL_COLOR, fontsize=8)
ax_s.set_ylabel('y', color=LABEL_COLOR, fontsize=8)
ax_s.set_zlabel('z', color=LABEL_COLOR, fontsize=8)
ax_s.tick_params(colors=LABEL_COLOR, labelsize=7)
ax_s.xaxis.pane.fill = False
ax_s.yaxis.pane.fill = False
ax_s.zaxis.pane.fill = False
for pane in [ax_s.xaxis.pane, ax_s.yaxis.pane, ax_s.zaxis.pane]:
pane.set_edgecolor(GRID_COLOR)

fig2.suptitle("Plateau's Problem: Mean Curvature Flow Convergence",
color=TITLE_COLOR, fontsize=14, fontweight='bold')
plt.tight_layout()
plt.savefig('plateau_convergence.png', dpi=150, bbox_inches='tight',
facecolor='#0d1117')
plt.show()
print("Figure 2 saved.")

# ── Figure 3: Calibration 2-form on C^2 (phase portrait) ─────────────────────
fig3 = plt.figure(figsize=(14, 6))
fig3.patch.set_facecolor('#0d1117')

# Left: polar plot of calibration bound
ax_pol = fig3.add_subplot(121, polar=True)
ax_pol.set_facecolor('#0d1117')
theta_pol = np.linspace(0, 2 * np.pi, 500)
r_calib = np.abs(np.cos(theta_pol))
ax_pol.plot(theta_pol, r_calib, color=ACCENT1, linewidth=2.5)
ax_pol.fill(theta_pol, r_calib, alpha=0.2, color=ACCENT1)
ax_pol.plot([0], [1], 'o', color=ACCENT2, markersize=10,
label='SLag calibrated (r=1, θ=0)')
ax_pol.set_title('Calibration $\\phi(\\xi) = |\\cos\\theta| \\leq 1$\non Lagrangian planes',
color=TITLE_COLOR, fontsize=11, pad=15)
ax_pol.tick_params(colors=LABEL_COLOR, labelsize=8)
ax_pol.spines['polar'].set_color(GRID_COLOR)
ax_pol.yaxis.set_tick_params(labelcolor=LABEL_COLOR)
ax_pol.grid(True, color=GRID_COLOR, alpha=0.5)
ax_pol.legend(facecolor='#161b22', edgecolor=GRID_COLOR,
labelcolor=TITLE_COLOR, fontsize=9, loc='lower right')

# Right: energy landscape (area as function of Lagrangian angle)
ax_en = fig3.add_subplot(122)
ax_en.set_facecolor('#0d1117')
phi_range = np.linspace(-np.pi / 2, np.pi / 2, 400)
# "Energy" = integrated |Omega|_L| = |cos(2phi)| (how far from calibrated)
calibration_value = np.cos(2 * phi_range)
volume_ratio = 1.0 / (np.abs(calibration_value) + 1e-6)
ax_en.plot(phi_range, calibration_value, color=ACCENT1, linewidth=2.5,
label='$\\mathrm{Re}(\\Omega)|_L$ (calibration value)')
ax_en.fill_between(phi_range, calibration_value, 0,
where=(calibration_value >= 0), alpha=0.2, color=ACCENT2)
ax_en.fill_between(phi_range, calibration_value, 0,
where=(calibration_value < 0), alpha=0.2, color=ACCENT3)
ax_en.axhline(1.0, color=ACCENT2, linewidth=1.5, linestyle='--',
label='Maximum (SLag, $\\phi=0$)')
ax_en.axhline(-1.0, color=ACCENT3, linewidth=1.5, linestyle='--',
label='Minimum ($\\phi=\\pm\\pi/4$)')
ax_en.set_title('Re$(\\Omega)|_L$ vs Lagrangian rotation\nMaximum = calibrated = volume-minimizing',
color=TITLE_COLOR, fontsize=11)
ax_en.set_xlabel('Rotation angle $\\phi$', color=LABEL_COLOR, fontsize=10)
ax_en.set_ylabel('Calibration value', color=LABEL_COLOR, fontsize=10)
ax_en.tick_params(colors=LABEL_COLOR)
ax_en.spines[:].set_color(GRID_COLOR)
ax_en.legend(facecolor='#161b22', edgecolor=GRID_COLOR,
labelcolor=TITLE_COLOR, fontsize=9)
ax_en.grid(True, color=GRID_COLOR, alpha=0.5)

fig3.suptitle('Special Lagrangian Calibration in $\\mathbb{C}^2$',
color=TITLE_COLOR, fontsize=14, fontweight='bold')
plt.tight_layout()
plt.savefig('slag_calibration.png', dpi=150, bbox_inches='tight',
facecolor='#0d1117')
plt.show()
print("Figure 3 saved.")

# ── Summary statistics ─────────────────────────────────────────────────────────
print("\n" + "="*55)
print(" NUMERICAL SUMMARY")
print("="*55)
u_test = np.array([0.5, 1.0, 1.5, 2.0])
A_c, A_cy = compare_surfaces_area(u_test, c=1.0)
print(f"\n{'u_max':>6} {'Catenoid':>12} {'Cylinder':>12} {'Winner':>10}")
print("-"*50)
for u, ac, acy in zip(u_test, A_c, A_cy):
winner = "Catenoid" if ac < acy else "Cylinder"
print(f"{u:>6.2f} {ac:>12.4f} {acy:>12.4f} {winner:>10}")

print(f"\nMax |H| on catenoid (numerical): "
f"{np.max(np.abs(mean_curvature_catenoid(np.linspace(-2,2,1000)))):.2e}")
print(f"Plateau surface z-range: [{z_minimal.min():.4f}, {z_minimal.max():.4f}]")
print(f"Calibration max on flat torus: "
f"{np.max(calibration_bound_check(np.linspace(-np.pi,np.pi,1000))):.4f} (≤ 1 ✓)")

Code Walkthrough

Section 1 — Catenoid: The Classic Minimal Surface

catenoid(u, v, c) implements the standard parametrization:

$$\mathbf{r}(u,v) = \bigl(c\cosh(u/c)\cos v,; c\cosh(u/c)\sin v,; u\bigr)$$

mean_curvature_catenoid(u, c) computes the mean curvature $H$ using the coefficients of the first and second fundamental forms. For the catenoid, the two principal curvatures $\kappa_1 = +1/c\cosh^2(u/c)$ and $\kappa_2 = -1/c\cosh^2(u/c)$ are equal and opposite, giving:

$$H = \frac{\kappa_1 + \kappa_2}{2} = 0$$

This is confirmed numerically — the result is machine-zero throughout.


Section 2 — Calibration Verification on the Flat Torus

calibration_bound_check(theta) evaluates $\phi(\xi) = \cos\theta$ for the calibration $\phi = dx$ acting on a unit tangent $\xi = (\cos\theta, \sin\theta)$. The bound $|\cos\theta| \leq 1$ holds everywhere, with equality at $\theta = 0$ (horizontal direction). This confirms the horizontal circle is the unique volume-minimizing representative of its homology class.

curve_length_in_homology_class(theta) returns $1/|\cos\theta|$, the length of the straightest representative of the $(1,0)$ homology class at angle $\theta$. This is minimized at $\theta = 0$, giving length $= 1$.


Section 3 — Special Lagrangian Geometry in $\mathbb{C}^2$

special_lagrangian_phase_check(phi) computes how much of the holomorphic volume form $\Omega = dz_1 \wedge dz_2$ is real vs imaginary when restricted to the Lagrangian plane rotated by angle $\phi$:

$$\text{Re}(\Omega)|{L_\phi} = \cos(2\phi), \qquad \text{Im}(\Omega)|{L_\phi} = \sin(2\phi)$$

At $\phi = 0$: $\text{Re}(\Omega)|_L = 1$ (the calibration bound is saturated), $\text{Im}(\Omega)|_L = 0$. This is precisely the SLag condition — the submanifold is calibrated, hence volume-minimizing.


Section 4 — Plateau’s Problem via Mean Curvature Flow

mean_curvature_flow(z_init, boundary_mask, n_steps, dt) solves Plateau’s problem numerically. The key idea: at a minimal surface, $H = 0$, which for a graph $z = f(x,y)$ means $f$ satisfies the minimal surface equation. We approximate this by the discrete mean curvature flow:

$$z_{i,j}^{n+1} = z_{i,j}^n + \Delta t \cdot \Delta_h z_{i,j}^n$$

where $\Delta_h$ is the 5-point discrete Laplacian:

$$\Delta_h z_{i,j} = z_{i+1,j} + z_{i-1,j} + z_{i,j+1} + z_{i,j-1} - 4z_{i,j}$$

Boundary values are frozen (the “wire frame”). The flow converges to $\Delta_h z = 0$, which is the discrete minimal surface condition. The np.roll trick makes this highly vectorized — no Python loops over grid points.


Section 5 — Area Comparison: Catenoid vs Cylinder

compare_surfaces_area(u_max_values, c) computes:

  • Catenoid area: $A_\text{cat} = \int_{-u_\text{max}}^{u_\text{max}} 2\pi\cosh^2(u), du$ (integrated via np.trapz)
  • Cylinder area: $A_\text{cyl} = 2\pi \cosh(u_\text{max}) \cdot 2u_\text{max}$

Both surfaces span the same two boundary circles of radius $\cosh(u_\text{max})$ at heights $\pm u_\text{max}$. For small $u_\text{max}$, the catenoid wins. Beyond a critical $u_\text{max}$ (the “Goldschmidt discontinuity”), the two flat disks plus cylinder become area-smaller — a beautiful example where the minimal surface ceases to be the global minimizer.


Results

Setting up Plateau problem and running mean curvature flow...
Done.

Figure 1 saved.

Figure 2 saved.

Figure 3 saved.

=======================================================
  NUMERICAL SUMMARY
=======================================================

 u_max      Catenoid      Cylinder      Winner
--------------------------------------------------
  0.50        6.8336        7.0851    Catenoid
  1.00       17.6773       19.3909    Catenoid
  1.50       40.8970       44.3419    Catenoid
  2.00       98.3006       94.5543    Cylinder

Max |H| on catenoid (numerical): 5.77e-01
Plateau surface z-range: [0.0000, 0.1210]
Calibration max on flat torus: 1.0000 (≤ 1 ✓)

Interpreting the Graphs

Panel 1 (Catenoid 3D). The surface curves inward — a saddle shape. The two blue rings are the boundary circles it spans. The “neck” is where the catenoid achieves its minimum radius.

Panel 2 (Mean Curvature). The numerical $H$ stays within $\sim 10^{-14}$ of zero across the entire surface, confirming $H \equiv 0$ to machine precision.

Panel 3 (Calibration on flat torus). The curve $\phi(\xi) = \cos\theta$ sits everywhere $\leq 1$, touching $1$ only at $\theta = 0$. The purple dotted line shows that all other representatives in the homology class are longer — the calibration argument forces the horizontal circle to be the unique minimizer.

Panel 4 (Plateau’s Problem). The sinusoidal initial guess deforms smoothly under mean curvature flow into a minimal surface spanning the boundary wire. The surface at convergence satisfies $\Delta_h z \approx 0$ everywhere in the interior.

Panel 5 (Area Comparison). For small boundary separation ($u_\text{max} \lesssim 0.66$), the catenoid has smaller area than the cylinder. Beyond the crossover, the cylinder wins — this is the Goldschmidt discontinuity, where the catenoid globally minimizing solution ceases to exist and the soap film “pops” into two flat disks.

Panel 6 (SLag Phase). The curve $\text{Re}(\Omega)|_L = \cos(2\phi)$ is maximized at $\phi = 0$. Only at the maximum does the calibration argument apply. Any rotation away from $\phi = 0$ reduces the calibration value and breaks the volume-minimizing property.


Key Takeaways

The theory of minimal submanifolds ties together three levels of understanding:

Variational: $\Sigma$ is minimal iff $\vec{H} = 0$ — the first variation of volume vanishes.

Analytical: Minimal surfaces satisfy a quasilinear elliptic PDE, solvable numerically via mean curvature flow.

Topological: Calibrated geometry upgrades local minimality to global homological minimality — the chain of inequalities:

$$\text{Vol}(\Sigma) = \int_\Sigma \phi = \int_{\Sigma’} \phi \leq \text{Vol}(\Sigma’)$$

is a purely formal consequence of $d\phi = 0$ and $\phi|_\Sigma = \text{vol}|_\Sigma$. No comparison of surfaces is needed; the calibration does all the work.