Add scripts for exact Casimir projectors and recoupling proof

- Implemented `2_exact_casimir_projectors.py` to construct exact SO(3) isotypic projectors using the Casimir operator, replacing Monte-Carlo methods.
- Created `3_apply_exact_projectors.py` to apply the exact projectors to example states, calculating ||A_j||_* estimates with improved precision.
- Developed `4_spherical_basis.py` to build a Condon-Shortley-consistent spherical basis for a single spin-1 leg using ladder operators.
- Introduced `6_six_j_recoupling_proof.py` to provide a complete proof of the cut-recoupling formula for full collective SU(2) symmetry, verifying the relationship between reduced blocks A_j^(1) and A_p^(2).
- Added a README file to document the execution order and purpose of each script in the symmetric states project.

scripts/symmetric_states/4_spherical_basis.py

@@ -0,0 +1,65 @@
+"""
+Step 1: build a Condon-Shortley-consistent spherical basis {|1,-1>,|1,0>,|1,+1>}
+for a single spin-1 leg, starting from the EXACT Cartesian generators
+(J_a)_{bc} = -i eps_{abc} (already verified symbolically in
+exact_casimir_projectors.py), and using the ladder-operator construction
+so we do not have to trust a memorized phase convention.
+"""
+import numpy as np
+
+def J_component(a):
+    eps = np.zeros((3,3,3))
+    eps[0,1,2]=eps[1,2,0]=eps[2,0,1]=1
+    eps[0,2,1]=eps[2,1,0]=eps[1,0,2]=-1
+    M = np.zeros((3,3), dtype=complex)
+    for b in range(3):
+        for c in range(3):
+            M[b,c] = -1j*eps[a,b,c]
+    return M
+
+Jx, Jy, Jz = J_component(0), J_component(1), J_component(2)
+Jp = Jx + 1j*Jy   # raising
+Jm = Jx - 1j*Jy   # lowering
+
+# sanity
+print("[Jx,Jy]-i Jz max err:", np.abs(Jx@Jy-Jy@Jx - 1j*Jz).max())
+print("J^2 (single leg), should be 2*I:")
+print(np.round(Jx@Jx+Jy@Jy+Jz@Jz,6))
+
+# eigenvectors of Jz
+evals, evecs = np.linalg.eigh(Jz)  # Jz Hermitian? check
+print("Jz Hermitian check:", np.abs(Jz - Jz.conj().T).max())
+print("Jz eigenvalues:", np.round(evals,6))
+
+# pick |1,-1> = eigenvector with eigenvalue closest to -1, fix phase: first
+# nonzero component real positive
+idx_m1 = np.argmin(np.abs(evals - (-1)))
+v_m1 = evecs[:, idx_m1]
+# fix global phase
+k = np.argmax(np.abs(v_m1))
+v_m1 = v_m1 * np.exp(-1j*np.angle(v_m1[k]))
+if v_m1[k].real < 0: v_m1 = -v_m1
+print("\n|1,-1> (Cartesian components x,y,z):", np.round(v_m1,4))
+
+# ladder up: |1,0> = Jp|1,-1> / ||...||    (standard CS convention: J+|j,m>=sqrt((j-m)(j+m+1))|j,m+1>, positive real coefficient)
+v0_raw = Jp @ v_m1
+n0 = np.linalg.norm(v0_raw)
+v_0 = v0_raw / n0
+print("|1,0>  raw ladder norm (expect sqrt((1-(-1))*(1+(-1)+1))=sqrt(2)):", n0)
+
+v_p1_raw = Jp @ v_0
+n_p1 = np.linalg.norm(v_p1_raw)
+v_p1 = v_p1_raw / n_p1
+print("|1,+1> raw ladder norm (expect sqrt((1-0)*(1+0+1))=sqrt(2)):", n_p1)
+
+# check orthonormality and Jz eigenvalues
+basis = np.stack([v_m1, v_0, v_p1], axis=1)  # columns
+print("\northonormality check (should be I_3):")
+print(np.round(basis.conj().T @ basis, 6))
+
+for name, v, m in [("|1,-1>", v_m1, -1), ("|1,0>", v_0, 0), ("|1,+1>", v_p1, 1)]:
+    Jzv = Jz @ v
+    print(f"{name}: Jz|.> - {m}|.>  max err = {np.abs(Jzv - m*v).max():.2e}")
+
+np.save("spherical_basis_single_leg.npy", basis)  # columns m=-1,0,+1
+print("\nsaved spherical_basis_single_leg.npy (columns ordered m=-1,0,+1)")