The Quaternion Screw:
Galaxy Rotation, Spiral Arms, and the Origin of Dark Matter
Galaxy Rotation, Spiral Arms, and the Origin of Dark Matter
Independent Researcher April 2026
We show that the non-commutativity of Hamilton’s quaternion algebra, applied to gravitational interactions in the framework of the exponential infall cosmology developed in Paper 1 of this series, generates angular momentum from initially non-rotating matter distributions. The composition of gravitational influences from multiple sources is a quaternion product, not a vector sum. The Baker–Campbell–Hausdorff commutator of non-collinear gravitational quaternions produces a net torque perpendicular to the plane of interaction—a screw force intrinsic to the algebra. This screw force provides a natural explanation for three phenomena currently attributed to dark matter: (i) flat galaxy rotation curves, (ii) the spontaneous formation of spiral arms, and (iii) the MOND acceleration scale a0. We derive a0 = cH/(2e) = 1.253 × 10–10 m/s2, matching the observed value 1.2 × 10–10 m/s2 to 4.4%, with zero free parameters. Every quantity in the formula—
1Introduction: The Missing Mass
In 1933, Fritz Zwicky measured the velocities of galaxies in the Coma Cluster and found that they were moving too fast [1]. The visible mass of the cluster—all the stars, all the gas—could not generate enough gravity to hold the galaxies in their orbits. He concluded that the cluster must contain unseen dunkle Materie: dark matter, outweighing the visible matter by a factor of at least ten. Forty years later, Vera Rubin and Kent Ford measured rotation curves of individual spiral galaxies [2]. The results were equally disturbing. In Newtonian gravity, the orbital velocity of a star at distance r from the galactic center should fall as v ∝ 1/√r once it is outside the bulk of the mass—just as planets farther from the Sun orbit more slowly. Instead, Rubin found that v(r) stays approximately constant out to the farthest measurable radius. The rotation curves are flat. The standard resolution is the same as Zwicky’s: postulate a vast halo of invisible matter surrounding each galaxy, carefully arranged so that the enclosed mass grows linearly with radius, M(r) ∝ r, producing v = constant. This dark matter halo must contain roughly five times as much mass as all visible matter [3]. Despite decades of searching—underground detectors, particle colliders, satellite experiments—no dark matter particle has ever been directly observed. In 1983, Mordehai Milgrom proposed an alternative: Modified Newtonian Dynamics, or MOND [4]. He showed that all observed rotation curves could be explained if Newton’s second law is modified below a critical acceleration a0 ≈ 1.2 × 10–10 m/s2. MOND works remarkably well—but the acceleration scale a0 has no known origin. It is a number that falls from the sky. We propose that it falls from the algebra.
2The Quaternion Composition of Gravity
In Cartesian mechanics, the gravitational forces from multiple sources add as vectors: F = F1 + F2 + … . Vector addition is commutative: F1 + F2 = F2 + F1. The order does not matter. In the quaternion framework of this series, spacetime is a quaternion Q = ct + xi + yj + zk, and the gravitational influence of a mass on a test particle is a quaternion operator. The combined effect of two gravitational sources is not a sum but a product: exp(H1) · exp(H2). Quaternion multiplication is not commutative: ij = k but ji = –k. The order matters. The Baker–Campbell–Hausdorff formula quantifies the difference. For two quaternion operators Ha and Hb:
exp(H_a) · exp(H_b) = exp(H_a + H_b + ½[H_a, H_b] + …)
The commutator [Ha, Hb] = HaHb – HbHa is computed from the quaternion product formula. For two gravitational quaternions with imaginary (spatial) parts u and v:
[H_a, H_b] = 2(u × v)
This is a vector perpendicular to both gravitational pulls. It is a torque. It generates angular momentum in a direction orthogonal to the plane defined by the two sources and the test particle. In Cartesian mechanics, this term does not exist—it is identically zero because vector addition commutes. In quaternion mechanics, it is as fundamental as the forces themselves.
3Angular Momentum from Nothing
Consider a uniform distribution of matter—a gas cloud, a primordial disk—with zero net angular momentum. In Cartesian gravity, if the system starts with L = 0, it stays at L = 0 forever. Angular momentum is strictly conserved, and symmetric initial conditions produce symmetric collapse: every particle falls radially inward, and the cloud implodes to a point. In quaternion gravity, the same initial conditions produce a different outcome. Each particle is subject to gravitational pulls from all other particles. These pulls come from different directions in quaternion space. For each pair of sources (j, k) acting on particle i, the commutator [Hj, Hk] contributes a torque perpendicular to the plane of the three bodies. The total torque is the sum over all pairs:
τ_total = ∑_{j For a perfectly uniform, perfectly symmetric distribution, this sum vanishes by symmetry. But any slight perturbation—a density fluctuation, an off-center clump, the kind of irregularity that the real universe is full of—breaks the cancellation and produces a net torque. Angular momentum is generated from zero. The cloud begins to rotate. This is not a violation of any conservation law. Total angular momentum in quaternion space (the full four-dimensional quantity) is conserved. What we observe as angular momentum in three dimensions is the projection of a four-dimensional quantity. The “missing” angular momentum lives in the time component of the quaternion—the fourth axis that Cartesian mechanics does not have. The commutator torque has a natural physical interpretation: it is a screw. A particle falling radially inward under the combined quaternion influence of multiple off-axis sources acquires a tangential velocity. The infall spirals. The analogy is precise. A screw converts linear motion into rotational motion through a helix—a curve that advances along one axis while rotating around it. The quaternion commutator does exactly this: it converts radial gravitational infall (linear, along the spatial quaternion axes) into tangential motion (rotational, perpendicular to the radial direction). The pitch of the screw—the ratio of advance per turn—is set by the strength of the gravitational coupling relative to the commutator correction. The screw force is not a new fundamental interaction. It is a geometric consequence of gravity being quaternion-valued rather than vector-valued. It costs nothing—no new particles, no new fields, no new parameters. It is latent in Hamilton’s algebra, waiting for 180 years to be applied to galactic dynamics. At what acceleration does the screw force become significant? The Newtonian gravitational acceleration falls as 1/r2. The screw force, arising from the commutator, is a second-order correction—it scales as the product of two gravitational influences. At small radii (strong gravity), Newton dominates overwhelmingly and the screw is negligible. At large radii (weak gravity), the screw becomes comparable to and eventually exceeds the Newtonian contribution. The transition occurs at a critical acceleration a0. We can determine this scale from the constants already established in this series. The exponential infall cosmology of Paper 1 [5] established that the metric scale factor is: a(d) = e^(Hd/ Three constants appear in this formula. The speed of light c is the conversion factor between space and time. The Hubble constant H is the curvature parameter—derived, not assumed, in Paper 1 from the requirement of constant Gaussian curvature. And the base of the exponential is Euler’s number e, because constant curvature means exponential: d(ex)/dx = ex is the definition. The quaternion structure adds one more number: 2. This is (2e) Numerically, with The observed MOND value is a0 = 1.2 × 10–10 m/s2 [4]. The agreement is 4.4%. Table 1. The four ingredients of the MOND acceleration. None is adjustable. Honesty requires one caveat. Dimensional analysis fixes cH as the scale of In Newtonian gravity, the centripetal acceleration for a circular orbit is a = v2/r = GM(r)/r2. When a > a0 (inner galaxy, strong field), the Newtonian expression holds and v ∝ 1/√r. When a < a0 (outer galaxy, weak field), the screw force dominates and the effective acceleration becomes: a_eff = √(a_Newton · This is exactly the MOND interpolation. Setting aeff = v2/r and aNewton = GM/r2: v⁴ = GM · The velocity depends on M and a0 but not on r. The rotation curve is flat. This is the Tully–Fisher relation [6]—the observed correlation between galaxy luminosity (proportional to mass) and rotation velocity—derived here from first principles. For the Milky Way (M ≈ 6 × 1010 M☉, vobs ≈ 220 km/s): v = (GM · Within 5% of the observed value, using only the visible mass and the derived acceleration scale. No dark matter halo required. The screw force explains the rotation. But galaxies are not featureless rotating disks—they have spiral arms. Where do the arms come from? The answer is local gravity. Once the screw sets the disk rotating, neighboring particles attract each other gravitationally. Particles that happen to be slightly closer together fall toward each other, forming clumps. These clumps, now denser than their surroundings, attract still more material, growing into elongated streams. The differential rotation of the disk—inner regions orbiting faster than outer regions—stretches these streams into arcs. An initially radial overdensity is sheared into a trailing spiral. The arms are not rigid structures; they are density waves, continuously forming, stretching, and dissolving as matter flows through them. One classical objection must be met head-on before anything else: the winding problem. If arms were material structures — the same stars riding the same arm — differential rotation would wrap them into tight coils within two or three revolutions, and a mature disk has made dozens (in an eternal galaxy, the problem is infinitely acute). Arms therefore cannot be objects. They must be patterns, continuously regenerated, through which matter flows. Any account of spiral structure that does not say where the pattern is regenerated has not explained it. The modern account supplies the machinery, and it is worth stating because it favors this framework. Lin and Shu proposed quasi-stationary density waves [10]; the picture that survived simulation is Toomre’s swing amplification [11]: in a differentially rotating, self-gravitating disk near marginal stability (Toomre Q ≈ 1–2), any leading disturbance is sheared into a trailing one and amplified by factors of tens on the way. The arms this produces are transient and recurrent — forming, winding, dissolving, reforming — exactly as long-run The quaternion commutator has a sign. The product ij = k defines a right-handed coordinate system; ji = –k defines the opposite. Left multiplication by exp( Let us trace the logical chain from Hamilton to galaxy rotation curves: Step 1 (Hamilton, 1843): Quaternions exist. Four-dimensional numbers with the multiplication rule ij = k, jk = i, ki = j. The multiplication is non-commutative [7]. Step 2 (Paper 1, 2026): Spacetime is a quaternion. The gravitational metric contracts exponentially: a(d) = eHd/ Galaxy rotation curves are the most famous evidence for dark matter, but not the only evidence. We briefly address the other pillars: Galaxy cluster dynamics (Zwicky). The velocity dispersion of galaxies in clusters exceeds what visible mass can bind. In our framework, the same screw force that flattens rotation curves also provides additional effective binding in clusters. The MOND acceleration a0 applies at the cluster scale just as it does at the galaxy scale—it is a universal constant derived from Table 2. Predictions of the quaternion screw model. Dark matter was invented to explain why galaxies rotate too fast. The explanation we offer is simpler: they rotate because quaternion multiplication does not commute. The non-commutativity of Hamilton’s quaternions, applied to the composition of gravitational influences, produces a screw force that generates angular momentum from initially non-rotating matter. The critical acceleration below which this force dominates is a0 = cH/(2e) = 1.253 × 10–10 m/s2, matching the empirical MOND value to 4.4%. Spiral arms form from gravitational clumping in the screw-driven rotating disk. The handedness of the screw is fixed by the algebra, offering a geometric origin for parity violation. In the completed series the screw also has a named geometric home: torsion — the antisymmetric part of the connection in the Einstein–Cartan extension of relativity, where spacetime twist is sourced by spin. Carrying the screw from the force law into the connection is part of [1] F. Zwicky, “Die Rotverschiebung von extragalaktischen Nebeln,” Helvetica Physica Acta, vol. 6, pp. 110–127, 1933. [2] V. C. Rubin and W. K. Ford Jr., “Rotation of the Andromeda Nebula from a Spectroscopic Survey of Emission Regions,” The Astrophysical Journal, vol. 159, pp. 379–403, 1970. [3] S. Navas et al. (Particle Data Group), “Review of Particle Physics,” Physical Review D, vol. 110, 030001, 2024. [4] M. Milgrom, “A modification of the Newtonian dynamics as a possible alternative to the hidden mass hypothesis,” The Astrophysical Journal, vol. 270, pp. 365–370, 1983. [5] M. Scholl, “Cosmological Redshift as Gravitational Metric Contraction: An Exponential Infall Model with Quaternion Geometry,” unpublished manuscript, 2026. [6] R. B. Tully and J. R. Fisher, “A new method of determining distances to galaxies,” Astronomy and Astrophysics, vol. 54, no. 3, pp. 661–673, 1977. [7] W. R. Hamilton, “On quaternions; or on a new system of imaginaries in algebra,” Philosophical Magazine, vol. 25, no. 3, pp. 489–495, 1844. [8] M. Scholl, “The Quantum Leap as Quaternion Transfiguration: E = hf as Geometric Identity,” unpublished manuscript, 2026. [9] G. W. Angus, B. Famaey, and H. S. Zhao, “Can MOND take a bullet? Analytical comparisons of three versions of MOND beyond spherical symmetry,” Monthly Notices of the Royal Astronomical Society, vol. 371, no. 1, pp. 138–146, 2006. [10] C. C. Lin and F. H. Shu, “On the Spiral Structure of Disk Galaxies,” The Astrophysical Journal, vol. 140, pp. 646–655, 1964. [11] A. Toomre, “What amplifies the spirals?,” in The Structure and Evolution of Normal Galaxies, S. M. Fall and D. Lynden-Bell, eds., Cambridge University Press, pp. 111–136, 1981. [12] J. A. Sellwood and R. G. Carlberg, “Spiral instabilities provoked by accretion and star formation,” The Astrophysical Journal, vol. 282, pp. 61–74, 1984. [13] J. P. Ostriker and P. J. E. Peebles, “A Numerical Study of the Stability of Flattened Galaxies: or, can Cold Galaxies Survive?,” The Astrophysical Journal, vol. 186, pp. 467–480, 1973. [14] R. Brada and M. Milgrom, “The stability of disc galaxies in the modified dynamics,” The Astrophysical Journal, vol. 519, pp. 590–598, 1999. [15] K. Land et al., “Galaxy Zoo: the large-scale spin statistics of spiral galaxies in the Sloan Digital Sky Survey,” Monthly Notices of the Royal Astronomical Society, vol. 388, pp. 1686–1692, 2008. [16] M. J. Longo, “Detection of a dipole in the handedness of spiral galaxies with redshifts z ~ 0.04,” Physics Letters B, vol. 699, pp. 224–229, 2011. The program that produced Figures 1 and 2 is supplied alongside this paper as QuaternionScrew_simulation.py, and it is the program, not a sketch of it. It is a particle-mesh code: forty thousand particles deposited by cloud-in-cell onto a 512×512 grid, the potential solved in Fourier space with the thin-disk kernel Φ_k = −2πGΣ_k/k and a Gaussian softening of 0.25 length units, the acceleration interpolated back to the particles. Every parameter is declared in one disclosed block at the head of the file — The two stages are distinct experiments and are worth separating. Stage A starts from a uniform disk at rest — every velocity exactly zero, so the total angular momentum L_z is exactly zero, not small — and applies the screw as a rotation of the velocity itself: a_screw = ω(v_y, −v_x) with ω = s·√(a_N·4The Screw Force
5The MOND Acceleration Scale
6Flat Rotation Curves
7Why Galaxies Have Arms
8Parity Violation and the Weak Force
9The Chain of Derivation
10What Dark Matter Was Supposed to Explain
11Summary of Predictions
12Conclusion
References
Appendix A. N-Body Simulation Code
QuaternionScrew_simulation.py — QuaternionScrew_simulation
[recorded run truncated — this script is checkpointed and continues past the recording window]
# -*- coding: utf-8 -*-
"""
Quaternion-screw disk galaxy simulation — particle-mesh (FFT) version.
N = 40,000 particles on a 512x512 grid, thin-disk 3D gravity kernel Phi_k = -2 pi G Sigma_k / k.
Stage A: rotation from exactly zero angular momentum (velocity-coupled screw).
Stage B: arm formation in the rotationally supported screw disk (radial MOND enhancement, no halo).
Checkpointing: resumes from state.npz; run repeatedly until 'DONE'.
All parameters disclosed below.
"""
import numpy as np, time, os, sys
# ---------------- parameters (all disclosed) ----------------
G=1.0; Mtot=1.0; a0=0.02
N=40000; NG=512; BOX=60.0 # grid, box size (disk occupies centre)
soft=0.25 # softening (k-space Gaussian), in length units
dt=0.10
STEPS_A=600 # stage A: spin-up from L=0
STEPS_B=3000 # stage B: arm formation
Mb=0.35; Md=0.65; bp=1.2 # fixed bulge (Plummer), disk mass
s_screw=0.8 # stage-A screw frequency scaling
Q_disp=0.22 # stage-B init dispersion fraction of vc
vflat=(G*Mtot*a0)**0.25
TIME_BUDGET=33.0 # seconds per invocation
cell=BOX/NG
kx=2*np.pi*np.fft.fftfreq(NG,d=cell); ky=kx.copy()
KX,KY=np.meshgrid(kx,ky,indexing='ij'); K=np.hypot(KX,KY); K[0,0]=1.0
KERN=(-2*np.pi*G/K)*np.exp(-K*soft) # thin-disk kernel + softening
KERN[0,0]=0.0
def pm_accel(x,y,m):
gx=(x+BOX/2)/cell; gy=(y+BOX/2)/cell
i0=np.floor(gx).astype(np.int32)%NG; j0=np.floor(gy).astype(np.int32)%NG
fx=gx-np.floor(gx); fy=gy-np.floor(gy)
i1=(i0+1)%NG; j1=(j0+1)%NG
rho=np.zeros((NG,NG))
np.add.at(rho,(i0,j0),m*(1-fx)*(1-fy)); np.add.at(rho,(i1,j0),m*fx*(1-fy))
np.add.at(rho,(i0,j1),m*(1-fx)*fy); np.add.at(rho,(i1,j1),m*fx*fy)
rho/=cell*cell
phik=np.fft.fft2(rho)*KERN
axg=np.real(np.fft.ifft2(-1j*KX*phik)); ayg=np.real(np.fft.ifft2(-1j*KY*phik))
ax=(axg[i0,j0]*(1-fx)*(1-fy)+axg[i1,j0]*fx*(1-fy)+axg[i0,j1]*(1-fx)*fy+axg[i1,j1]*fx*fy)
ay=(ayg[i0,j0]*(1-fx)*(1-fy)+ayg[i1,j0]*fx*(1-fy)+ayg[i0,j1]*(1-fx)*fy+ayg[i1,j1]*fx*fy)
return ax,ay
def bulge(x,y):
rb2=x*x+y*y+bp*bp
f=-G*Mb/rb2**1.5
return f*x,f*y
def full_accel(x,y,m,mond=True):
ax,ay=pm_accel(x,y,m)
bx,by=bulge(x,y)
ax+=bx; ay+=by
if mond:
aN=np.hypot(ax,ay)+1e-12
boost=0.5*(1+np.sqrt(1+4*a0/aN))
ax*=boost; ay*=boost
return ax,ay
ST='state.npz'
if os.path.exists(ST):
S=dict(np.load(ST,allow_pickle=True))
phase=str(S['phase']); step=int(S['step'])
x,y,vx,vy=S['x'],S['y'],S['vx'],S['vy']
LzA=list(S['LzA']); ttA=list(S['ttA'])
snapsteps=list(S['snapsteps'])
else:
np.random.seed(7)
th=np.random.uniform(0,2*np.pi,N); r0=9.0*np.sqrt(np.random.uniform(0.03,1,N))
x=r0*np.cos(th); y=r0*np.sin(th)
vx=np.zeros(N); vy=np.zeros(N)
phase='A'; step=0; LzA=[]; ttA=[]; snapsteps=[]
m=np.full(N,Md/N)
t0=time.time()
while True:
if time.time()-t0>TIME_BUDGET:
np.savez(ST,phase=phase,step=step,x=x,y=y,vx=vx,vy=vy,LzA=LzA,ttA=ttA,snapsteps=snapsteps)
print(f'CHECKPOINT phase={phase} step={step}'); sys.exit(0)
if phase=='A':
if step>=STEPS_A:
# transition: re-init stage B in rotational support (keep the A result recorded)
np.random.seed(11)
th=np.random.uniform(0,2*np.pi,N); r0=9.0*np.sqrt(np.random.uniform(0.03,1,N))
x=r0*np.cos(th); y=r0*np.sin(th)
ax,ay=full_accel(x,y,m)
ri=np.hypot(x,y)+1e-9; aR=np.hypot(ax,ay); vc=np.sqrt(aR*ri)
tx=-y/ri; ty=x/ri
rng=np.random.default_rng(3)
vx=vc*tx+rng.normal(0,1,N)*Q_disp*vc; vy=vc*ty+rng.normal(0,1,N)*Q_disp*vc
phase='B'; step=0
np.savez('stageB_t0.npz',x=x,y=y)
continue
ax,ay=full_accel(x,y,m,mond=False)
aN=np.hypot(ax,ay)+1e-12
w=s_screw*np.sqrt(aN*a0)/vflat
asx=+w*vy; asy=-w*vx
vx+=(ax+asx)*dt; vy+=(ay+asy)*dt
x+=vx*dt; y+=vy*dt
if step%5==0: LzA.append(float((m*(x*vy-y*vx)).sum())); ttA.append(step*dt)
step+=1
else:
if step>=STEPS_B:
np.savez('final.npz',x=x,y=y,vx=vx,vy=vy,LzA=LzA,ttA=ttA)
print('DONE'); sys.exit(0)
ax,ay=full_accel(x,y,m)
vx+=ax*dt; vy+=ay*dt
x+=vx*dt; y+=vy*dt
if step in (0,1000,2000,2999):
np.savez(f'snapB_{step}.npz',x=x,y=y)
step+=1
Appendix B. The Continuous Run, Open
QuaternionScrew_simulation_continuous.py — QuaternionScrew_simulation_continuous
[recorded run truncated — this script is checkpointed and continues past the recording window]
# -*- coding: utf-8 -*-
"""
CONTINUOUS quaternion-screw run: motionless uniform cloud (L = 0) -> rotation -> arms, one story.
Particle-mesh, N = 40,000, checkpointed (run repeatedly until DONE).
STATUS AND FINDINGS (July 2026 session; documented for further work):
The two-stage demonstration in the paper is solid. This continuous mode is an open
numerical problem. Four variants were run; each taught something:
1. Screw as raw tangential push -> runaway spin-up, disk ejected. The screw must saturate.
2. Screw as velocity rotation (all of v) -> after spin-up it scrambles circular orbits
(Stage A's own late-time |L_z| decline shows this).
3. Screw on radial velocity, both signs -> post-bounce OUTFLOW is counter-torqued; ejected
particles at large lever arm poison the L_z budget.
4. Screw on INFALL only (this file) -> inner disk spins up to ~the (G M a0)^(1/4) prediction,
but the cold-cloud collapse is violent (core bounce
ejects most mass). Violent relaxation, not the screw,
is now the obstacle.
NEXT: tame the collapse, not the screw. Two routes: (a) gradual mass assembly - grow the
cloud by accretion instead of dropping it whole (add particles over time); (b) a gas
pressure term (e.g. local-density repulsion) so the collapse is subsonic. Either should
let the infall-only screw (line ~77) convert calmly. Parameters s_screw ~ 0.8-1.2,
damp ~ 0.9999 were the best-behaved region.
"""
import numpy as np, time, os, sys
G=1.0; Mtot=1.0; a0=0.02
N=40000; NG=512; BOX=60.0; soft=0.25; dt=0.10
STEPS=4500
Mb=0.35; Md=0.65; bp=1.2
s_screw=1.2; damp=0.9999
vflat=(G*Mtot*a0)**0.25
TIME_BUDGET=33.0
SNAPS=(0,600,1200,2000,3000,4500)
cell=BOX/NG
kx=2*np.pi*np.fft.fftfreq(NG,d=cell)
KX,KY=np.meshgrid(kx,kx,indexing='ij'); K=np.hypot(KX,KY); K[0,0]=1.0
KERN=(-2*np.pi*G/K)*np.exp(-K*soft); KERN[0,0]=0.0
def accel(x,y,m):
gx=(x+BOX/2)/cell; gy=(y+BOX/2)/cell
i0=np.floor(gx).astype(np.int32)%NG; j0=np.floor(gy).astype(np.int32)%NG
fx=gx-np.floor(gx); fy=gy-np.floor(gy); i1=(i0+1)%NG; j1=(j0+1)%NG
rho=np.zeros((NG,NG))
np.add.at(rho,(i0,j0),m*(1-fx)*(1-fy)); np.add.at(rho,(i1,j0),m*fx*(1-fy))
np.add.at(rho,(i0,j1),m*(1-fx)*fy); np.add.at(rho,(i1,j1),m*fx*fy)
rho/=cell*cell
phik=np.fft.fft2(rho)*KERN
axg=np.real(np.fft.ifft2(-1j*KX*phik)); ayg=np.real(np.fft.ifft2(-1j*KY*phik))
ax=(axg[i0,j0]*(1-fx)*(1-fy)+axg[i1,j0]*fx*(1-fy)+axg[i0,j1]*(1-fx)*fy+axg[i1,j1]*fx*fy)
ay=(ayg[i0,j0]*(1-fx)*(1-fy)+ayg[i1,j0]*fx*(1-fy)+ayg[i0,j1]*(1-fx)*fy+ayg[i1,j1]*fx*fy)
rb2=x*x+y*y+bp*bp; f=-G*Mb/rb2**1.5
ax+=f*x; ay+=f*y
aN=np.hypot(ax,ay)+1e-12
boost=0.5*(1+np.sqrt(1+4*a0/aN))
return ax*boost,ay*boost,aN
def diag(x,y,step,Ms=(1,2,3,4,5,6)):
r=np.hypot(x,y)
lo,hi=np.percentile(r,[20,85])
sel=(r>lo)&(r<hi)
th=np.arctan2(y,x)[sel]; lr=np.log(r[sel])
A={M:abs(np.exp(1j*M*th).sum())/max(sel.sum(),1) for M in Ms}
ps=np.linspace(-30,30,601)
amp=np.array([abs(np.exp(1j*(2*th+p*lr)).sum()) for p in ps])
pb=ps[amp.argmax()]
import math
return [step*dt]+[A[M] for M in Ms]+[math.degrees(math.atan2(2,abs(pb)))]
ST='screw_continuous_state.npz' # unique per script: two runs must never share a checkpoint
if os.path.exists(ST):
S=dict(np.load(ST)); step=int(S['step'])
x,y,vx,vy=S['x'],S['y'],S['vx'],S['vy']
Lz=list(S['Lz']); tt=list(S['tt']); D=[list(row) for row in S['D']]
else:
np.random.seed(7)
th=np.random.uniform(0,2*np.pi,N); r0=9.0*np.sqrt(np.random.uniform(0.03,1,N))
x=r0*np.cos(th); y=r0*np.sin(th)
vx=np.zeros(N); vy=np.zeros(N)
step=0; Lz=[]; tt=[]; D=[]
m=np.full(N,Md/N); t0=time.time()
while True:
if time.time()-t0>TIME_BUDGET:
np.savez(ST,step=step,x=x,y=y,vx=vx,vy=vy,Lz=Lz,tt=tt,D=np.array(D))
print(f'CHECKPOINT step={step}'); sys.exit(0)
if step>STEPS:
np.savez('screw_continuous_final.npz',x=x,y=y,vx=vx,vy=vy,Lz=Lz,tt=tt,D=np.array(D))
print('DONE'); sys.exit(0)
ax,ay,aN=accel(x,y,m)
w=s_screw*np.sqrt(aN*a0)/vflat
ri=np.hypot(x,y)+1e-9
rx=x/ri; ry=y/ri
vr=vx*rx+vy*ry # radial velocity component
tx=-ry; ty=rx
vin=np.minimum(vr,0.0); asx=-w*vin*tx; asy=-w*vin*ty # deflect infall into rotation; dormant on circular orbits
vx=(vx+(ax+asx)*dt)*damp; vy=(vy+(ay+asy)*dt)*damp
x+=vx*dt; y+=vy*dt
if step%25==0: Lz.append(float((m*(x*vy-y*vx)).sum())); tt.append(step*dt)
if step%250==0: D.append(diag(x,y,step))
if step in SNAPS: np.savez(f'screw_continuous_snap_{step}.npz',x=x,y=y)
step+=1
Symbols & Terms