# SCRIPT: BH-SUNLADDER-16b # Loops around the first n exceptional points of the Sun's dynamo operator: does a loop climb the ladder # (as the near-extremal Kerr cascade does) or is it a product of separate swaps? # Operator: the solar paper's minimal Parker alpha-Omega operator, the same matrix as SUN-OPEN-02 in # suns_v1_5_complete.zip: x in [0, pi], N = 120 interior points, Dirichlet, # dA/dt = alpha0 cos(x) B + A'', dB/dt = A' + B'' (G = eta = 1), second-order central differences. # The leading 16 eigenvalues come from sparse shift-invert at 0 (checked equal to dense eigenvalues to 1e-11). # # Why this decides the question. Every loop starts and ends at the base point alpha0 = 3.0 on the real axis, # below the first birth point, where all tracked eigenvalues are real and labelled 0, 1, 2, ... by size # (0 = the fundamental, the least damped). A loop returns the eigenvalues as a permutation fixed by the # branch points it encloses. A lasso around one branch point gives one transposition (a, b). A loop around # the first n branch points gives the ordered product of the n lassos. The fork, written before the run: # CLIMB if consecutive branch points share a label, e.g. (1,2), (2,3), ...: the product is one long cycle # and a single loop carries an eigenvalue up the ladder, as in Kerr; # SWAPS if the transpositions are disjoint, e.g. (1,2), (3,4), ...: the product is a set of separate swaps # and no eigenvalue moves more than one place. # # Tracking. Steps are adaptive in arc length: a step is accepted only if every tracked eigenvalue moves less # than a quarter of its distance to the nearest other candidate, otherwise the step is halved (down to 1e-6). # This is the safeguard against the coarse-step artefact of the solar paper's v1.3 EP4/EP5 tracking. # Eigenvector overlap || is recorded as a second check; it is not used to assign, because the two # eigenvectors of a pair become parallel at the branch point, where eigenvalue distance stays decisive. # # Kill conditions, written before the run: # K1 control: the real-axis census in (0.5, 13.5] reproduces the solar paper's five birth points # 3.872, 5.280, 7.158, 9.085, 11.032 to 1e-2 (the paper quotes 3 decimals; its loops used radius 0.12); # every real-axis branch point (births AND deaths) is listed, and any not in the paper is reported. # K2 each lasso around one branch point gives exactly one transposition. # K3 the discriminant of the tracked eigenvalues winds once per enclosed branch point: winding = n for a loop # around the first n. If not, the loop caught a branch point off the axis and the permutation is not read. # K4 the loop around the first n gives the same permutation at heights 0.05 and 0.2, equals the ordered # product of the n lassos, and the reversed loop gives the inverse permutation. # Kill (no verdict) if K2, K3 or K4 fails; then report what was found. import numpy as np, scipy.sparse as sp, scipy.sparse.linalg as sla N=120; h=np.pi/(N+1); x=np.arange(1,N+1)*h D2=sp.diags([np.ones(N-1),-2*np.ones(N),np.ones(N-1)],[-1,0,1])/h**2 D1=sp.diags([-np.ones(N-1),np.ones(N-1)],[-1,1])/(2*h) Cc=sp.diags(np.cos(x)) K=16 def spec(a, vec=False): A=sp.bmat([[D2,a*Cc],[D1,D2]]).tocsc().astype(complex) if vec: w,V=sla.eigs(A,k=K+4,sigma=0,which='LM'); o=np.argsort(-w.real)[:K]; return w[o],V[:,o] w=sla.eigs(A,k=K+4,sigma=0,which='LM',return_eigenvectors=False); return w[np.argsort(-w.real)][:K] ok=True def report(name,cond,val=''): global ok print(('PASS ' if cond else 'FAIL ')+name, val, flush=True); ok=ok and cond # ---- dense check of the sparse solver ---- wd=np.linalg.eigvals(sp.bmat([[D2,(3+0.1j)*Cc],[D1,D2]]).toarray()); wd=wd[np.argsort(-wd.real)][:K] print(f' sparse vs dense eigenvalues at alpha0 = 3+0.1i: {np.abs(spec(3+0.1j)-wd).max():.1e}', flush=True) # ---- K1: real-axis census, births and deaths ---- def nreal(a): w=spec(a); return int(np.sum(np.abs(w.imag)<1e-7*np.maximum(1,np.abs(w)))) grid=np.linspace(0.5,13.5,2601); cnt=[nreal(a) for a in grid] bps=[] for i in range(len(grid)-1): if cnt[i+1]!=cnt[i]: lo,hi=grid[i],grid[i+1]; c0=cnt[i] for _ in range(45): m=(lo+hi)/2 if nreal(m)==c0: lo=m else: hi=m bps.append(((lo+hi)/2,'birth' if cnt[i+1]1e-2]}') pts=[b for b,_ in bps] # ---- tracking ---- def track(path_fn, s0, s1, ds0=0.01): s=s0; cur,V=spec(path_fn(s),True); start=cur.copy(); disc=[np.sum(np.log((cur[:,None]-cur[None,:])[np.triu_indices(K,1)]**2))] minov=1.0; ds=ds0 while s0.25*d[1]: okstep=False; break if okstep or step<1e-6: break step/=2 Vn=W[:,asg] ov=np.abs(np.sum(np.conj(V)*Vn,axis=0))/(np.linalg.norm(V,axis=0)*np.linalg.norm(Vn,axis=0)); minov=min(minov,ov.min()) cur=new[asg]; V=Vn; s+=step d=np.sum(np.log((cur[:,None]-cur[None,:])[np.triu_indices(K,1)]**2)) disc.append(disc[-1]+np.angle(np.exp(1j*(d.imag-disc[-1].imag)))*1j+(d.real-disc[-1].real)) ds=min(ds0,step*1.5) perm=[int(np.argmin(np.abs(start-c))) for c in cur] err=max(np.min(np.abs(start-c)) for c in cur) wind=(disc[-1].imag-disc[0].imag)/(2*np.pi) return perm,wind,err,minov base=3.0 def stadium(a1,ht,rev=False): # base -> along the lower edge to a1 -> right cap -> upper edge back to base -> left cap closes at base segs=[] L1=a1-base; Lc=np.pi*ht def P(s): if s<=ht: return base-1j*s # drop from the axis s-=ht if s<=L1: return base+s-1j*ht s-=L1 if s<=Lc: th=-np.pi/2+s/ht; return a1+ht*np.exp(1j*th) s-=Lc if s<=L1: return a1-s+1j*ht s-=L1 return base+1j*(ht-s) # rise back to the axis T=2*ht+2*L1+Lc if rev: return (lambda s: P(T-s)),T return P,T def lasso(k,ht=0.05): # base -> lower edge to just left of point k -> small circle around it -> back the same way p=pts[k]; gaps=[abs(p-q) for j,q in enumerate(pts) if j!=k]; r=min(0.12,0.45*min(gaps)) if gaps else 0.12 a=p-r; L1=a-base; Lc=2*np.pi*r def P(s): if s<=ht: return base-1j*s s-=ht if s<=L1: return base+s-1j*ht s-=L1 if s<=ht: return a-1j*(ht-s) # up to the axis left of the point s-=ht if s<=Lc: th=np.pi+s/r; return p+r*np.exp(1j*th) # around it, counterclockwise s-=Lc if s<=ht: return a-1j*s s-=ht if s<=L1: return a-s-1j*ht s-=L1 return base-1j*(ht-s) return P,4*ht+2*L1+Lc def cyc(perm): seen=set(); out=[] for i in range(len(perm)): if i in seen: continue c=[i]; seen.add(i); j=perm[i] while j!=i: c.append(j); seen.add(j); j=perm[j] if len(c)>1: out.append(c) return out def compose(p,q): # first p then q: element i goes to p[i], then to q[p[i]] return [q[p[i]] for i in range(len(p))] def inverse(p): r=[0]*len(p) for i,j in enumerate(p): r[j]=i return r nmax=min(5,len(pts)) las=[] for k in range(nmax): P,T=lasso(k); perm,w,err,mo=track(P,0,T); las.append(perm) print(f' lasso around point {k+1} (alpha0 = {pts[k]:.4f}, {bps[k][1]}): cycles {cyc(perm)} winding {w:+.3f} return err {err:.1e} min overlap {mo:.3f}', flush=True) report('K2 each lasso is one transposition', all(len(cyc(p))==1 and len(cyc(p)[0])==2 for p in las), '') allok=True; verdicts=[] for n in range(1,nmax+1): a1=(pts[n-1]+pts[n])/2 if n1 sw=all(len(c)==2 for c in cs) and len(cs)==n verdicts.append((n,cs,'CLIMB' if lad else 'SWAPS' if sw else 'OTHER')) print(f' loop around first {n}: h=0.05 cycles {cs} winding {res[0.05][1]:+.3f} | h=0.2 cycles {cyc(res[0.2][0])} winding {res[0.2][1]:+.3f} | reversed {cyc(rv[0])} | product of lassos {cyc(prod)}', flush=True) report('K3+K4 windings, two heights, lasso product and reversal all consistent', allok, '') for n,cs,v in verdicts: print(f' verdict n={n}: {v} {cs}') # K5 (report): is the extra birth point near 12.97 a property of the operator or of the grid? def births_at(Nn,lo=12.0,hi=13.5,M=301): hh=np.pi/(Nn+1); xx=np.arange(1,Nn+1)*hh d2=sp.diags([np.ones(Nn-1),-2*np.ones(Nn),np.ones(Nn-1)],[-1,0,1])/hh**2 d1=sp.diags([-np.ones(Nn-1),np.ones(Nn-1)],[-1,1])/(2*hh); cc=sp.diags(np.cos(xx)) def nr(a): w=sla.eigs(sp.bmat([[d2,a*cc],[d1,d2]]).tocsc().astype(complex),k=K+4,sigma=0,which='LM',return_eigenvectors=False) w=w[np.argsort(-w.real)][:K]; return int(np.sum(np.abs(w.imag)<1e-7*np.maximum(1,np.abs(w)))) g=np.linspace(lo,hi,M); c=[nr(a) for a in g]; out=[] for i in range(M-1): if c[i+1]