"""
Landau minimization with winding suppression: three-generation window.
We minimize
F({s_n}) = sum_n [ (1/2) r_n s_n + (u/4) s_n^2 ] + v * sum_{n<m} s_n s_m,
with s_n = |a_n|^2 >= 0.
Key modeling choice (winding suppression):
r_n(c_n) = alpha + c_n * (n-1)^2,
so increasing c_n penalizes higher winding (n).
For suitable parameters (alpha < 0, small v > 0), there exists an interval
in c_n where exactly three modes condense (s_n > 0 for n=1,2,3; others 0).
This script:
• solves the KKT system analytically for every active set (no local minima),
• scans c_n, plots r_n(c_n), the optimized s_n(c_n), and the condensed-mode count,
• shades the interval where exactly three modes are condensed,
• prints Hessian eigenvalues (convexity) and summary of the 3-mode window.
Dependencies: numpy, matplotlib (SciPy optional; not required).
"""
import itertools
import numpy as np
import matplotlib.pyplot as plt
# -------------------------
# Model parameters (edit here if desired)
# -------------------------
N = 5 # number of modes (n = 1..N)
u = 1.0 # quartic self-coupling (>0)
v = 0.05 # inter-mode coupling (small, >=0). Keep v < u/2 for convexity.
alpha = -1.0 # base offset for r_n after curvature/back-reaction (alpha < 0)
# c_n scan (choose a range that crosses the 3-mode window)
CN_MIN, CN_MAX, CN_NUM = 0.05, 0.35, 200
# Numerical tolerances
TOL_NEG = 1e-9
TOL_ACT = 1e-9
# -------------------------
# Helpers: r_n(c_n), energy, analytic active-set solver
# -------------------------
modes = np.arange(1, N+1)
def r_of_cn(cn: float) -> np.ndarray:
# Winding suppression: higher n receive larger (positive) offset as c_n grows
return alpha + cn * (modes - 1)**2
def energy(s: np.ndarray, r: np.ndarray) -> float:
# F = 1/2 r·s + u/4 sum s^2 + v * sum_{i<j} s_i s_j
diag = 0.5 * np.dot(r, s) + 0.25 * u * np.dot(s, s)
cross = v * np.sum([s[i] * s[j] for i in range(N) for j in range(i+1, N)])
return float(diag + cross)
def solve_active_set(r: np.ndarray):
"""
Enumerate all active sets A ⊆ {0..N-1}, enforce KKT conditions analytically:
For i in A:
0 = ∂F/∂s_i = 1/2 r_i + 1/2 u s_i + v sum_{j in A, j≠i} s_j
For k not in A:
∂F/∂s_k |_{s_k=0} = 1/2 r_k + v sum_{j in A} s_j >= 0
and s_i >= 0.
With s_i = |a_i|^2, define S_A ≡ sum_{j in A} s_j.
Solve: s_i = -(r_i + 2 v S_A) / (u - 2 v), with
S_A = - sum_{i in A} r_i / (u - 2 v + 2 v |A|).
Returns (s*, A*, F*).
"""
if u - 2*v <= 0:
raise ValueError("Convexity requires u - 2 v > 0.")
best_F = np.inf
best_s = np.zeros(N)
best_A = tuple()
for mask in range(1 << N):
A = tuple(i for i in range(N) if (mask >> i) & 1)
if not A:
# All-zero candidate valid if gradients at 0 are nonnegative: r_i >= 0 for all i
if np.all(r >= -TOL_NEG):
return np.zeros(N), A, 0.0
continue
denom = (u - 2*v)
S_A = -np.sum(r[list(A)]) / (denom + 2*v*len(A))
s = np.zeros(N)
ok = True
# Active components
for i in A:
s_i = -(r[i] + 2*v*S_A) / denom
if s_i < -TOL_NEG:
ok = False
break
s[i] = max(0.0, s_i)
if not ok:
continue
# Inactive KKT: gradient at zero must be >= 0
grad_ok = True
for k in range(N):
if k in A:
continue
if 0.5*r[k] + v*S_A < -TOL_NEG:
grad_ok = False
break
if not grad_ok:
continue
F = energy(s, r)
if F < best_F - 1e-12:
best_F, best_s, best_A = F, s, A
return best_s, best_A, best_F
# -------------------------
# Scan c_n and collect results
# -------------------------
cns = np.linspace(CN_MIN, CN_MAX, CN_NUM)
R = np.zeros((CN_NUM, N))
S = np.zeros((CN_NUM, N))
Aks = [()] * CN_NUM
for k, cn in enumerate(cns):
r = r_of_cn(cn)
s, A, F = solve_active_set(r)
R[k, :] = r
S[k, :] = s
Aks[k] = A
counts = (S > TOL_ACT).sum(axis=1)
idx3 = np.where(counts == 3)[0]
cn_window = (cns[idx3[0]], cns[idx3[-1]]) if idx3.size > 0 else (None, None)
# -------------------------
# Convexity / stability: Hessian eigenvalues
# H_ii = u/2; H_ij = v (i ≠ j). Eigenvalues:
# λ1 = u/2 + (N-1) v (multiplicity 1) and λ⊥ = u/2 - v (multiplicity N-1).
# -------------------------
lam1 = 0.5*u + (N-1)*v
lamper = 0.5*u - v
print("=== Hessian eigenvalues (global convexity) ===")
print(f"λ_parallel = {lam1:.6f}, λ_perp = {lamper:.6f} -> positive definite? {lam1>0 and lamper>0}")
# Report the 3-mode window and a representative point
print("\n=== Three-generation window (exactly three condensed modes) ===")
if idx3.size > 0:
print(f"c_n in [{cn_window[0]:.6f}, {cn_window[1]:.6f}] (length ≈ {cn_window[1]-cn_window[0]:.6f})")
mid = idx3[len(idx3)//2]
print(f"Example at c_n = {cns[mid]:.6f}:")
print(" r_n:", np.round(R[mid], 6))
print(" |a_n|^2:", np.round(S[mid], 6))
print(" active modes (1-indexed):", [i+1 for i in Aks[mid]])
else:
print("No 3-mode window found with current parameters. Increase c_n range or adjust alpha.")
# -------------------------
# Plotting
# -------------------------
fig, axes = plt.subplots(3, 1, figsize=(10, 10), constrained_layout=True)
# 1) r_n vs c_n
for i in range(N):
axes[0].plot(cns, R[:, i], label=f"r_{i+1}")
axes[0].axhline(0.0, linestyle="--")
axes[0].set_xlabel("Winding suppression c_n")
axes[0].set_ylabel("Quadratic coefficients r_n")
axes[0].set_title("Quadratic Coefficients vs c_n")
axes[0].legend()
# 2) VEVs |a_n|^2 vs c_n
for i in range(N):
axes[1].plot(cns, S[:, i], label=f"|a_{i+1}|^2")
axes[1].set_xlabel("Winding suppression c_n")
axes[1].set_ylabel("Optimized |a_n|^2")
axes[1].set_title("Vacuum Expectation Values")
axes[1].legend()
# 3) Count of condensed modes
axes[2].plot(cns, counts, marker="o")
axes[2].axhline(3, linestyle="--", label="Three generations")
axes[2].set_ylim(0, N + 1)
axes[2].set_xlabel("Winding suppression c_n")
axes[2].set_ylabel("Number of condensed modes")
axes[2].set_title("Mode Count")
axes[2].legend()
# Shade the 3-mode window, if present
if idx3.size > 0:
for ax in axes:
ax.axvspan(cn_window[0], cn_window[1], alpha=0.12)
plt.show()
=== Hessian eigenvalues (global convexity) ===
λ_parallel = 0.700000, λ_perp = 0.450000 -> positive definite? True
=== Three-generation window (exactly three condensed modes) ===
c_n in [0.087688, 0.208291] (length ≈ 0.120603)
Example at c_n = 0.147990:
r_n: [-1. -0.85201 -0.40804 0.33191 1.367839]
|a_n|^2: [0.901847 0.737414 0.244114 0. 0. ]
active modes (1-indexed): [1, 2, 3]