Geodesic Monism II
Here is the second post on Geodesic Monism, a completely geometrical theory of reality.
As the first post I asked Claude, Fable, to compute what I was lazy to finish for years/decades,
and it comes up with its own interpretation as well, so the tone below is pessimistic,
where I am very optimistic about it! I like it's being doubtful anyway so let's publish it!
I also gave it to itself for audit and it only could complain about and compare with trashy theories,
so I am happy with the achievements here!
I think the core idea would get zero out of 10 in the Sabina Hossenfelder's bs meter,
but Fable didn't do only the calculations as I mentioned, so let it be!
For instance, I have no idea why Calabi–Yau is in this post.
Even human slop needs your filter, right?
If I get good feedback I will explain a little bit of Quantum Mechanics with this theory, sooner than later! I am specifically proud of it being the first theory to show you exactly what the higher dimensions looks like and where to find them, check the galaxy bar section. Let's read Fable's writings.
In the Geodesic Monism post1 we put two principles on the table, found the hydrogen solution, and watched it face the galactic and cosmological data. This post is the promised deepening of its electromagnetic sector: we run the Kaluza8 program on the theory all the way through, with a completely arbitrary transverse metric, the way Kaluza himself would have demanded. Three things fall out. The Maxwell equations, derived instead of read off. The spin back-reaction coefficients \( 1/12 \) and \( 1/6 \), which the machine discovered numerically last time, now derived by hand in two lines. And one theorem we did not see coming: the theory forces the internal space to be Ricci-flat, which upgrades our gauge-group daydreams from numerology to something with a spine, and hands us, uninvited, the same class of internal spaces that string theory compactifies on. And past the Kaluza program, four more hunts report in below: a five dimensional sibling for the trapped sector and the group-theory whisper it carries, the empirical law that the galactic coupling \( \beta \) turns out to obey, what fourth-order electrodynamics does to light itself, and a first quantitative confrontation between the galactic bar and the collapse boundary of the extra dimensions.
We also pay two debts. Last time we compared the rotation curves only against Kepler, and a reader rightly called that a weak baseline; this time the hydrogen faces MOND and the dark matter halos on the same data. And three numbers from the previous post were wrong or inconsistent; they are fixed below, in the open, with the corrected code.
The two principles, restated
Principle 1: Everything moves on geodesics. There is no matter, there is only spacetime. What we call matter is a pattern in the metric field, and what we call dynamics is geodesic motion in that pattern.
Principle 2: The action is the square of the Ricci curvature,
\[ S = \int R_{ab}R^{ab} \sqrt{-g} \mathrm{d}^n x \]
whose equations of motion \( E_{ab} = 0 \) are
\[ E_{ab} = \Box R_{ab} + \tfrac{1}{2} g_{ab} \Box R - \nabla_a \nabla_b R + 2 R_{acbd} R^{cd} - \tfrac{1}{2} g_{ab} R_{cd}R^{cd} = 0 \]
Every term contains at least one Ricci tensor, so every vacuum solution of General Relativity is inherited for free, and every solution that is not Ricci-flat is new physics.
Is this action new? The honesty ledger
Short answer: the action is not new; the reading of it is. Quadratic curvature actions go back a century, to Weyl and Bach, and the general family \( \alpha R^2 + \beta R_{ab}R^{ab} \) is the quadratic gravity that Stelle proved renormalizable in 19772 and solved classically in 19783. Our action is the pure \( \beta \) corner of that family. Even our beloved solution mechanism has ancestors: pp-waves and their Kundt cousins are known exact solutions of generic quadratic gravity, studied by Málek and Pravda4, by Pravda, Pravdová, Podolský and Švarc5, and the vanishing of all curvature invariants on them is the VSI property classified by Pravda, Pravdová, Coley and Milson6. Critical gravity7 lives one shelf over. Anyone who tells you a two-derivative-squared action is virgin territory is selling something.
What we could not find anywhere, and what this program actually is: taking the pure Ricci-squared theory as the entire ontology, with no matter fields at all, not even as sources; reading the non-Ricci-flat solutions as the matter content of the world; and running that reading against real astronomical data. The hydrogen solution family with its trapped extra dimensions, the null Kaluza program of this post, and the theorems below are, to the best of our literature hunting, new. If a reader finds them elsewhere, the honesty section at the bottom is where the correction will live.
Einstein's equation, read backwards
Einstein wrote
\[ G_{ab} = T_{ab} \]
with geometry on the left and a postulated inventory of stuff on the right: dust here, fields there, an equation of state everywhere, each with its own continuity equation bolted on. In Geodesic Monism the same equation is read backwards. Take any solution of \( E_{ab} = 0 \), compute its Einstein tensor, and define
\[ T^{\mathrm{eff}}_{ab} \equiv G_{ab} = R_{ab} - \tfrac12 g_{ab} R \]
The right hand side is built from the metric and nothing else: a completely geometrical equation for \( T_{ab} \). Two things come for free that Einstein had to demand. Conservation, \( \nabla^a T^{\mathrm{eff}}_{ab} = 0 \), holds identically by the contracted Bianchi identity, so the effective matter never needs a continuity equation, it cannot help being conserved. And the type of matter is not postulated but computed: move the terms of \( E_{ab} = 0 \) across the equals sign and the stress tensor is found to obey its own wave equation,
\[ \Box T_{ab} = \nabla_a \nabla_b R - 2 R_{acbd} R^{cd} + \tfrac{1}{2} g_{ab} R_{cd}R^{cd} \]
Matter is not sourced here; it propagates, a wave of curvature whose self-interaction is curvature squared. (Bookkeeping for the sticklers: the exact rearrangement of \( E_{ab} = 0 \) carries an extra \( -g_{ab}\Box R \), and the trace of the field equations makes \( \Box R \) itself curvature squared, \( \Box R = \tfrac{n-4}{n} R_{cd}R^{cd} \), so in four dimensions the equation is exactly as displayed, while in \( n \) dimensions the last coefficient reads \( \tfrac{8-n}{2n} \).) For the whole solution class of this post the computation collapses further (Theorem 1 below, with \( k = \mathrm{d}u \) the null direction and hatted operators living on the transverse space): the effective matter is a null dust with energy density \( -\hat\Delta H + \tfrac14 F_{kl}F^{kl} \) streaming along \( k \), carrying momentum flux \( W_i = \tfrac12 \hat\nabla^k F_{ik} \). And read through the Kaluza dictionary \( A_i = g_{ui} \), that momentum flux is the electric current, \( W_i = -\tfrac12 J_i \): the effective matter current and the effective charge current are the same component of the same tensor. In this theory the old suspicion that charge always rides on matter is not a coincidence of our universe, it is an identity of the bookkeeping.
The Kaluza method, all the way through
Last time we hard-coded the transverse space flat and extracted electromagnetism from a linearized kernel. The honest Kaluza method keeps the base metric fully dynamical, so this time the ansatz is
\[ ds^2 = 2 \mathrm{d}t \mathrm{d}u + (1+2H) \mathrm{d}u^2 + 2A_i \mathrm{d}x^i \mathrm{d}u + h_{ij}(x) \mathrm{d}x^i \mathrm{d}x^j \]
with \( H \), \( A_i \) and the transverse metric \( h_{ij} \) arbitrary functions of the transverse coordinates. One structural remark before the theorem: in this null version of Kaluza the role of the four dimensional spacetime block is played by a Riemannian metric \( h_{ij} \), because the Lorentzian signature is entirely consumed by the null pairing \( 2 \mathrm{d}t \mathrm{d}u \). That signature fact will grow teeth in a moment.
Two lemmas from the previous post survive the generalization untouched, because their proofs never used flatness: the entire \( t \) row of the metric is still the constant \( \delta^u_b \), so \( g^{uu} = g^{ui} = 0 \) and \( {\Gamma^u}_{ab} = 0 \) exactly (Proposition 3 of the previous post), and therefore \( k = \mathrm{d}u \) is covariantly constant. Chasing the Christoffel symbols, every one with a lower \( t \) index vanishes, the transverse ones are pure, \( {\Gamma^k}_{ij} = {\hat\Gamma^k}{}_{ij}[h] \), and the cross ones carry the field strength, \( {\Gamma^k}_{ui} = \tfrac12 F_i{}^k \), indices moved by \( h \).
That cross Christoffel deserves its own sentence, because dropped into the geodesic equation it is the Lorentz force:
\[ \ddot x^k = \dot u F^k{}_i \dot x^i + \partial^k H \dot u^2 + \cdots \]
with the conserved drift \( \dot u \) playing exactly the role of the electric charge (charge-to-mass ratio, strictly): a particle with no hidden-velocity memory feels no magnetic force, a particle with much of it curls hard, and electromagnetism's signature force is revealed as ordinary free fall reading the cross term of the metric. Out of this bookkeeping comes the main result.
Theorem 1 (null Kaluza reduction). For the ansatz above, the field equations \( E_{ab} = 0 \) hold if and only if
(i) \( \mathcal{R}_{ij}[h] = 0 \): the transverse space is exactly Ricci-flat;
(ii) \( \hat\Box \hat\nabla^k F_{ki} = 0 \): biharmonic Maxwell on \( (M_\perp, h) \), with \( \hat\Box = h^{kl}\hat\nabla_k \hat\nabla_l \) the rough Laplacian;
(iii) \( \hat\Delta\big( \hat\Delta H - \tfrac14 F^2 \big) = \tfrac12 J_k J^k + F^{ik}\hat\nabla_i J_k \), with \( J_i = \hat\nabla^k F_{ki} \) the Maxwell-reading current.
Proof. By hand for the structure, by machine for every component identity, on a deliberately curved non-Ricci-flat base with non-solution profiles so that nothing vanishes by accident, and again on Euclidean Schwarzschild as a curved Ricci-flat base: Appendices A and B. \( \blacksquare \)
The skeleton of the proof deserves daylight, because each vertebra is a small surprise.
The transverse block is pure and closed. The lowered Riemann tensor of the full metric restricted to transverse indices is exactly the Riemann tensor of \( h \): the \( F_{ik}F_{jl} \) contamination that ordinary spacelike Kaluza–Klein8 feeds into the base curvature is killed by the null fiber, because every such term must contract through \( g^{uu} = 0 \). Consequently \( R_{ij} = \mathcal{R}_{ij}[h] \), \( R = \mathcal{R}[h] \), and the transverse components of the field equations close on themselves:
\[ E_{ij} = \hat\Box\mathcal{R}_{ij} + \tfrac12 h_{ij}\hat\Delta\mathcal{R} - \hat\nabla_i\hat\nabla_j\mathcal{R} + 2\hat R_{ikjl}\mathcal{R}^{kl} - \tfrac12 h_{ij}\mathcal{R}_{kl}\mathcal{R}^{kl} = \mathcal{E}_{ij}[h] \]
which is the same Ricci-squared theory, one dimension down, with \( H \) and \( A_i \) completely absent. The theory dimensionally reduces into itself, and the electromagnetic sector cannot back-react on the base at all. In Einstein–Maxwell theory the \( F^2 \) stress-energy curves the base; here the coupling is strictly one-way, geometry dresses electromagnetism but effective charge cannot gravitate transversally, only feed \( g_{uu} \). This is exactly the consistency condition that the slogan "there is no charge, there is geometry" was quietly assuming.
A constraint appears, and it has Riemannian teeth. Because \( g_{tu} = 1 \), the \( (t,u) \) component of the field equations is no longer trivial:
\[ E_{tu} = \tfrac12\big( \hat\Delta\mathcal{R} - \mathcal{R}_{kl}\mathcal{R}^{kl} \big) \]
Combine with the trace of the transverse block, \( h^{ij}\mathcal{E}_{ij} = \tfrac d2 \hat\Delta\mathcal{R} + (2 - \tfrac d2)\mathcal{R}_{kl}\mathcal{R}^{kl} = 0 \), and \( E_{tu} = 0 \) becomes \( \mathcal{R}_{kl}\mathcal{R}^{kl} = 0 \). And now the signature remark bites: \( h \) is Riemannian, so \( \mathcal{R}_{kl}\mathcal{R}^{kl} \) is a sum of squares, and it vanishes only if \( \mathcal{R}_{ij} = 0 \) at every point. The "completely arbitrary" base metric of the Kaluza method is eaten by the field equations: the internal space must be exactly Ricci-flat. Flat was never a simplification in the previous post; it was one member of the only allowed family. On a Ricci-flat base all polynomial invariants die again, \( E_{ab} = \Box R_{ab} \) exactly, and clauses (ii) and (iii) are what \( \Box R_{ab} = 0 \) says in the Kaluza dictionary. One bookkeeping note so nobody accuses the split of being chart luck: under a gauge shift \( t \rightarrow t + \lambda \) the components mix as \( E_{ui} \rightarrow E_{ui} - \partial_i\lambda E_{tu} \), so the clean separation of the clauses is gauge-covariant precisely because \( E_{ta} = 0 \) and \( E_{tu} = 0 \) on-shell.
Maxwell is the decaying sector, exactly. Clause (ii) is fourth-order electrodynamics whose second-order factor is the Maxwell operator, with gauge invariance a theorem (Proposition 2 of the previous post) rather than an axiom. Read it in Maxwell's language: the current \( J_i \) does not have to vanish, it has to be harmonic. On \( \mathbb{R}^3 \), or on any base where the rough Laplacian has no decaying kernel, a harmonic current that decays at infinity is zero, so every globally decaying solution of the theory is an exact vacuum Maxwell solution. Maxwell is not approximately inside this theory; Maxwell is the decaying sector, and the extra modes are all growing modes. Which is precisely the structure of the gravitational sector, where Newton \( c_0/r \) is the decaying mode and the dark-sector physics \( c_2 r, c_3 r^2 \) are the growing biharmonic partners. One mechanism, both forces: fourth order equals second order plus growing modes.
The previous post's linearized kernel is now explained line by line. In the axial magnetic sector \( a_\phi = \sin^2\theta P(r) \) the Maxwell operator is \( M[P] = P'' - 2P/r^2 \), and clause (ii) reduces to \( M^2[P] = 0 \). Check the five monomials the machine was fed: \( M[1/r] = 0 \) and \( M[r^2] = 0 \), the Maxwell dipole and uniform field, free. \( M[r] = -2/r \) and \( M[-2/r] = 0 \), a genuine fourth-order partner, free. \( M[1] = -2/r^2 \) with \( M[-2/r^2] \neq 0 \), killed; \( M[r^3] = 4r \) with \( M[4r] \neq 0 \), killed. Free: \( b_0, b_2, b_3 \). Killed: \( b_1, b_4 \). That is, symbol for symbol, the kernel the machine printed last time, now a two-line calculation. (That the exact nonlinear hunts killed \( b_2 \) while the linear theory allows it stays interesting: growing partners can be obstructed nonlinearly, a fact worth its own hunt.)
The spin coefficients, derived. The previous post celebrated that the machine discovered the back-reaction \( H_J = \frac{J^2}{12 r^4} + \frac{J^2}{6}\frac{P_2(\cos\theta)}{r^4} \) over two finite fields. Here is the same result by hand. The spin term \( g_{u\phi} = J\sin^2\theta/r \) is, in the Kaluza dictionary, exactly the vector potential of a magnetic dipole of moment \( J \), whose energy density is \( \tfrac14 F^2 = \frac{J^2(1+3\cos^2\theta)}{2r^6} \). The dipole is an exact Maxwell solution, \( J_i = 0 \), so clause (iii) collapses to \( \hat\Delta(\hat\Delta H - \tfrac14 F^2) = 0 \), and the decaying solution is the Einstein-level magnetostatic identity
\[ \hat\Delta H_J = \tfrac14 F_{kl}F^{kl} \]
Apply \( \Delta (r^{-4}) = 12 r^{-6} \) and \( \Delta(P_2 r^{-4}) = 6 P_2 r^{-6} \), and the left side is \( J^2(1 + P_2)/r^6 = J^2(1+3\cos^2\theta)/(2r^6) \). The coefficients \( 1/12 \) and \( 1/6 \) are the magnetostatic self-energy of the dipole, nothing more. The machine's discovery is now an identity a student can check on paper, and the sentence from the previous post, "like the self-energy of a magnetic dipole", loses its "like".
One correction to our own working notes, in public because that is the deal. In deriving clause (iii) by hand we first obtained \( -\tfrac12 J^2 \) where the machine insists on \( +\tfrac12 J^2 \); the machine is right, and the fitting run that caught it is in Appendix B. The fit also returned a one-parameter family of equivalent forms, which is not ambiguity but the Bianchi identity at work: \( \hat\Delta F_{jk} = \hat\nabla_j J_k - \hat\nabla_k J_j \) plus curvature terms makes \( F\cdot\hat\Delta F \) and \( F\cdot\hat\nabla J \) proportional on-shell of \( \mathrm{d}F = 0 \), so the same equation wears two coats.
A free enlargement of the solution space. Any Ricci-flat base now serves: the verifier below literally exhibits hydrogen-type profiles living on a Euclidean Schwarzschild throat, and Taub–NUT, Eguchi–Hanson, every gravitational instanton, and every Calabi–Yau come along for free.
Eight trapped dimensions: adding five \( z \)'s
The previous post trapped three extra dimensions and got an \( SO(3) \). The obvious next hunt: trap five more, \( z_1,\dots,z_5 \), and watch what the equations do with eight. The result is now cheap to state rigorously, because Theorem 1 does the heavy lifting that used to require a bespoke high-dimensional curvature run.
Theorem 2 (the \( z \)-extended hydrogen). On the thirteen dimensional spacetime with coordinates \( (t, r, \theta, \phi, y_1..y_3, z_1..z_5, u) \) and metric
\[ ds^2 = 2\mathrm{d}t\mathrm{d}u + (1+2H)\mathrm{d}u^2 + \frac{2J\sin^2\theta}{r}\mathrm{d}\phi \mathrm{d}u + \mathrm{d}r^2 + r^2\mathrm{d}\Omega_2^2 + \mathrm{d}\vec y^{ 2} + \mathrm{d}\vec z^{ 2} \]
\[ H = \frac{c_0}{r} + c_1 + c_2 r + c_3 r^2 + \frac{J^2}{12r^4} + \frac{J^2 P_2(\cos\theta)}{6 r^4} - \Big(\gamma_0 + \frac{\gamma_1}{r}\Big)|\vec y|^2 - \gamma_2\big(r|\vec y|^2 - r^3\big) - \Big(\delta_0 + \frac{\delta_1}{r}\Big)|\vec z|^2 - \delta_2\Big(r|\vec z|^2 - \tfrac{5}{3}r^3\Big) - \sum_{a=1}^{5} f_a(r) z_a Q_a(\vec y) \]
where the \( Q_a \) are the five traceless quadratics of \( \vec y \) and each \( f_a \) is any biharmonic radial profile \( \epsilon^a_0 + \epsilon^a_1/r + \epsilon^a_2 r + \epsilon^a_3 r^2 \), the field equations \( E_{ab} = 0 \) hold identically for all parameter values.
Proof. Theorem 1 with flat transverse \( \mathbb{R}^{11} \). Clause (i) is trivial. Clause (ii) holds because the spin term is the gyraton9 cross term of the magnetic dipole, whose current vanishes. Clause (iii) splits by linearity into the dipole identity \( \hat\Delta H_J = \tfrac14 F^2 \), proven by hand in the Kaluza section, and the biharmonicity of every remaining term of \( H \), verified by machine with all parameters symbolic and with negative controls: Appendix D. \( \blacksquare \)
Notice what happened methodologically: the previous post needed an eight dimensional curvature computation to certify the hydrogen; Theorem 1 certifies the thirteen dimensional version from a handful of scalar Laplacians. The reduction theorem is now the workhorse, and the honest cost of a new trapped sector has dropped from a machine hunt to a lemma.
The anatomy carries dimension fingerprints. The trap identity generalizes to \( \Delta^2\big(g(r)|\vec z|^2\big) = |\vec z|^2\Delta^2 g + 4n \Delta g \) for \( n \) trapped dimensions, so the cross-coefficient is \( 12 \) for the \( y \)'s and \( 20 \) for the \( z \)'s, and the trap must be harmonic, \( \delta_0 + \delta_1/r \), in every dimension. The coupled mode obeys a clean law: for \( r|\vec z|^2 - \lambda r^3 \), one line of Laplacians gives \( \Delta^2 = (8n - 24\lambda)/r \), so
\[ \lambda_n = \frac{n}{3} \]
which is \( 1 \) for three dimensions, recovering the previous post's mode, and \( 5/3 \) for five; the machine confirms both and, as a negative control, correctly rejects each coefficient in the other's slot. Finally the surprise family: the coupling modes \( f_a(r) z_a Q_a(\vec y) \) are exact because \( z_a Q_a \) is harmonic in the eight trapped coordinates, so the radial profile only needs to be biharmonic, growing modes included. These couplings tie the \( \vec z \)'s to the quadrupole of the \( \vec y \)'s, and they will matter in the symmetry section.
Proposition 1 (two clocks, one rate). In the runaway regime, every trapped coordinate obeys the same Cauchy–Euler equation18, \( \ddot y = -(4\gamma_1/c_2)\tau^{-2} y \) and \( \ddot z = -(4\delta_1/c_2)\tau^{-2} z \), so both sectors are forgotten with the same exponent: envelopes \( \tau^{1/2} \), velocities \( \tau^{-1/2} \), coordinate rates \( dy/dt \sim dz/dt \sim r^{-5/4} \), independent of the number of trapped dimensions and of the trap strengths. Only the log-frequencies differ:
\[ \mu_y = \sqrt{\frac{4\gamma_1}{c_2} - \frac14}, \qquad \mu_z = \sqrt{\frac{4\delta_1}{c_2} - \frac14} \]
The exact geodesic integration (Appendix D) confirms it over three decades of proper time while \( r \) grows by five: both envelopes \( \tau^{1/2}|\dot y| \) and \( \tau^{1/2}|\dot z| \) stay bounded with no trend, the ratio \( |\dot z|/|\dot y| \) oscillates around one with no hierarchy, and the zeros of each sector space themselves by their own universal ratio:
y1 zero-crossing ratios (predict e^{pi/mu}=3.0910): [3.1316 3.0877 3.088 3.0899]
z1 zero-crossing ratios (predict e^{pi/mu}=2.1079): [2.106 2.1065 2.1072 2.1075]
The verdict is now a computed fact rather than a suspicion: a faster-fading sector cannot be built from more flat trapped dimensions. Flat traps produce a hierarchy of frequencies, never of exponents, which upgrades the previous post's strong-force flag from a worry to a theorem-adjacent statement, and leaves curved internal geometry, Ricci-flat by Theorem 1, as the only remaining road. The observable fingerprint of the \( 3+5 \) split is delicious though: two incommensurate log-clocks ticking in one and the same geodesic, a two-frequency discrete scale invariance.
\( SU(3) \times SU(2) \times U(1) \), honestly
The previous post ended its symmetry story at \( SU(2)\times U(1) \) and flagged \( SU(3) \) as needing "curved internal geometry". Time to audit the whole claim at the standard of rigor the rest of this program demands, and then to report what Theorem 1 adds.
What survives scrutiny: the three \( y \) dimensions carry an exact \( SO(3) \) isometry, built in by the common trap. The \( u \) direction carries the gauge shifts of the previous post's Proposition 2. So the honest symmetry statement is an isometry group \( SO(3) \times \mathbb{R} \) whose Lie algebra is \( \mathfrak{su}(2) \oplus \mathfrak{u}(1) \).
What does not survive: three separate overreaches. First, promoting \( SO(3) \) to its double cover \( SU(2) \) is only physics when spinorial degrees of freedom exist to feel the difference, and this theory has none; metric perturbations live in tensor representations, which cannot tell the two groups apart. Second, our \( u \) is noncompact, on purpose, so the gauge transformations \( A \rightarrow A + \mathrm{d}\lambda \) form \( \mathbb{R} \), not a circle; without compactness there is no honest \( U(1) \) and no charge quantization. Third, and worst, calling the product "the electroweak group" was a category error even on its own terms: the electroweak \( U(1) \) is hypercharge, while the photon's \( U(1) \) is the Weinberg-angle mixture; having identified our gauge direction with electromagnetism via gauge-invariance-as-coordinate-freedom, we were not entitled to relabel it hypercharge one section later. Consider that paragraph of the previous post corrected.
And one genuinely new piece of structure from the \( z \) hunt, stated at the same standard.
Proposition 2 (the \( \mathbf{3}\oplus\mathbf{5} \) skeleton). The coupling modes \( f_a(r) z_a Q_a(\vec y) \) of Theorem 2 break the generic isometry \( SO(3)_y \times SO(5)_z \) of the trapped sector down to the diagonal \( SO(3) \) acting on \( \vec y \) in the spin-1 and on \( \vec z \) in the spin-2 representation. Under this \( SO(3) \), the eight trapped dimensions transform as \( \mathbf{3}\oplus\mathbf{5} \), which is precisely the branching of the adjoint representation of \( SU(3) \) under its principal \( SO(3) \) subgroup.
Proof. The coupling term is invariant only under simultaneous rotations \( \vec y \rightarrow R\vec y \), \( z_a \rightarrow D^{(2)}(R)_{ab} z_b \), while the traps are invariant under anything orthogonal, so the residual continuous symmetry is this diagonal \( SO(3) \). For the branching: \( \mathbf{3}\otimes\bar{\mathbf{3}} = \mathbf{8}\oplus\mathbf{1} \), and restricted to the principal \( SO(3) \) this is \( (\mathrm{spin} 1)^{\otimes 2} = \mathrm{spin} 0 \oplus \mathrm{spin} 1 \oplus \mathrm{spin} 2 \), hence \( \mathbf{8} \rightarrow \mathbf{3}\oplus\mathbf{5} \). \( \blacksquare \)
So there exists an exact solution whose trapped sector is one adjoint of \( \mathfrak{su}(3) \)'s worth of dimensions, organized by a single \( SO(3) \), with the organization enforced by an explicit mode of the metric rather than by our wishes. Kept at the honesty level of this section: the actual symmetry is \( SO(3) \), the decomposition \( \mathbf{3}\oplus\mathbf{5} \) is also nothing more exotic than \( \ell=1 \oplus \ell=2 \), and no eight gauge bosons appear anywhere. One notch above numerology, several notches below dynamics.
Now the new part on the curved side. Theorem 1 says curved internal geometry is allowed, but only Ricci-flat. The hydrogen's transverse space is six dimensional, and the compact Ricci-flat six-manifolds with the richest structure are exactly those with holonomy11 group \( SU(3) \): the Calabi–Yau threefolds12, the same internal spaces that string theory reached by demanding unbroken supersymmetry10. Our field equations reach the same class by a completely different argument, a constraint in the null sector, with no supersymmetry anywhere in sight. So the corrected symmetry skeleton reads: \( \mathfrak{su}(2) \oplus \mathfrak{u}(1) \) as the isometry algebra of the flat trapped sector, and \( SU(3) \) available as the holonomy group of the allowed curved internal sector. Full honesty, in bold so it cannot be skimmed past: this is a structurally motivated skeleton, not a derivation of the standard model. There are no gauge bosons of the \( SO(3) \), no chiral fermions, no hypercharge, and holonomy is not a gauge symmetry. The improvement over the previous post is that the \( SU(3) \) slot is now backed by a theorem about which internal geometries the theory permits, instead of by wishful counting. The gap between skeleton and dynamics is the honest size of the remaining mountain.
What the fourth order does to light
We now own the electromagnetic field equations, so we can finally ask the theory what light does instead of assuming it. Four propositions follow; one of them killed our own favorite mechanism on arrival, and another handed us a better one.
Proposition 3 (charge is conserved, exactly, always). For any antisymmetric \( F_{\mu\nu} \) on any spacetime, \( \nabla_\mu \nabla_\nu F^{\nu\mu} = 0 \); hence the geometric current \( J^\mu = \nabla_\nu F^{\nu\mu} \) of this theory is identically conserved.
Proof. Split \( \nabla_\mu\nabla_\nu \) into its symmetric and antisymmetric parts; the symmetric part dies against the antisymmetry of \( F \), and the commutator produces one Ricci contraction against \( F \), zero because Ricci is symmetric, and one Riemann contraction, zero by the first Bianchi identity. Machine check for completely arbitrary \( A_\mu \): Appendix E. \( \blacksquare \)
The proposition closes one tempting door in every metric theory at once: what the theory abandons is everything classical about charge: \( J \) need not vanish in vacuum, only be harmonic, by clause (ii); it is not attached to worldlines, not localized, not quantized; the charge inside a sphere is a boundary condition that can grow with the sphere. Charge without charges, but conserved charge.
Proposition 4 (light inherits the redshift). In the eikonal limit \( A_\mu = a_\mu e^{iS/\varepsilon} \), the leading symbol of clause (ii) is \( \big( g^{ab} \partial_a S \partial_b S \big)^2 = 0 \), a doubled null characteristic. Phases ride null geodesics, so the frequency ratio between comoving observers is fixed by the Killing charges: exactly the endpoint redshift formula of the previous post, now derived for genuine electromagnetic waves rather than assumed for test photons.
Proposition 5 (partner waves break the inverse square law). The spherical field
\[ \Phi_p = \frac{(t+r) G(t-r)}{4r} \]
satisfies \( \Box \Phi_p = G'(t-r)/r \), an ordinary outgoing Maxwell wave, and \( \Box^2 \Phi_p = 0 \); its amplitude at fixed retarded time tends to \( G/2 \), with no \( 1/r \) falloff at all. (Machine check: Appendix E.)
The doubled characteristic of Proposition 4 is a Jordan block, and these are its off-diagonal modes: waves whose flux does not dilute. If astrophysical sources excite partner modes, luminosity bookkeeping changes while frequency does not, so the theory predicts violations of the distance duality \( d_L = d(1+z) \), the Etherington relation19. That moves the supernova side of the dark energy argument, not the redshift side, and it is falsifiable against existing distance-duality constraints. The fully propagating \( t \)-dependent analysis stays the announced sequel; these two propositions are its fixed boundary posts.
Proposition 6 (the back-reaction hierarchy). Clause (iii) makes the magnetic modes gravitate, and the growth of the sourced potential tracks how far the mode sits from the decaying Maxwell sector:
- dipole \( b_0 \): \( J = 0 \) and \( \hat\Delta H = \tfrac14 F^2 \) exactly, the bounded Einstein-level back-reaction, which is the spin sector of the hydrogen;
- uniform field \( b_3 \): \( J = 0 \) and \( F \) constant, right hand side zero, no back-reaction at all;
- fourth-order partner \( b_2 \): the current terms are exactly monopole-free, \( \tfrac12 J^2 + F\cdot\hat\nabla J = -8 b_2^2 P_2(\cos\theta)/r^4 \), but the Maxwell-energy term \( \tfrac14\hat\Delta F^2 \) keeps a monopole, and the full back-reaction is
\[ H = b_2^2 \Big( \ln r - \tfrac12 P_2(\cos\theta) \Big) \]
a logarithmically growing potential;
- a uniform harmonic current \( J_0 \): right hand side \( J_0^2 \), giving \( H = J_0^2 r^4/120 \), quartic growth.
Proof. By machine, all four cases: Appendix E. \( \blacksquare \)
So the constructive statement, sharper than any non-conservation mechanism could have been: nondecaying geometric currents and partner fields do feed the growing potentials that redshift light, through the \( J^2 \) and \( F^2 \) sourcing of clause (iii), and with a clean hierarchy: dipole bounded, partner logarithmic, uniform current quartic. The Hubble modes \( c_2 r \) and \( c_3 r^2 \) themselves remain source-free vacuum modes; the electromagnetic sector adds distinguishable textures on top, a logarithmic term in the Hubble relation being the most testable of them.
Dark matter: facing the stronger rivals
The previous post beat Kepler in 143 of 143 galaxies and a reader rightly pointed out that beating a one-parameter model with a two-parameter model is a low bar. So this time the hydrogen prediction \( v^2 = \alpha/r + \beta r \) faces the real competition on the same SPARC1314 data, same outer-region selection, same published error bars: MOND in its radial-acceleration-relation form15 with one free parameter (the disk mass-to-light ratio), and the two standard dark matter halos, NFW16 and Burkert17, each with two free halo parameters on top of the fixed-M/L baryonic curves. First, the audit: the original five-model pipeline reproduces every number of the previous post to the digit (code and logs in Appendix C), which is its own small comfort. For the consolidated table below, every row then runs through one unified code path, with exact closed-form nonnegative least squares replacing the earlier least-squares-then-clip; the unification moves the hydrogen baseline from 0.70 to 0.67 on the outer region and from 11.90 to 11.98 on full curves, disclosed here so that no row gains from solver luck. Five new rows join the original five: the constant-mode hydrogen \( \alpha/r + c_1 + \beta r \), the per-galaxy radius shift \( \alpha/(r+\delta r) + \beta(r+\delta r) \) of the bar section below, our own disciplined single-global-coupling variants (\( \delta r = \kappa r_1 \) and the softened radius \( \sqrt{r^2 + (\kappa r_1)^2} \), not part of the bar prescription, charged \( n-2 \) per galaxy with the one global \( \kappa \) amortized across the zoo), and the free softened-radius map. Code and the full per-galaxy table: Appendix F.
The row names, so the table reads without guessing. hydrogen(2p) is \( v^2 = \alpha/r + \beta r \), two parameters per galaxy. Kepler(1p) keeps only \( \alpha/r \). MOND/RAR(1p) is the radial acceleration relation with the disk mass-to-light ratio as its one free parameter. NFW(2p+bar) and Burkert(2p+bar) are the halo profiles, two free halo parameters each on top of the fixed mass-to-light baryonic curves, which is what the "+bar" flags: those rows consume the photometry. hydro+c1(3p) adds the constant mode, \( v^2 = \alpha/r + c_1 + \beta r \). hydro+dr(3p) is the radius shift \( v^2 = \alpha/(r+\delta r) + \beta(r+\delta r) \) with \( \delta r \) free per galaxy. hydro+dr=kr1(2p+g) is the same shift disciplined: \( \delta r = \kappa r_1 \), with \( r_1 = \sqrt{\alpha/\beta} \) taken from the outer fit and one global \( \kappa \) shared by all 143 galaxies, hence two per-galaxy parameters plus one global, which is what "(2p+g)" means. hydro+soft(3p) softens the radius instead, \( r \rightarrow \sqrt{r^2 + R_b^2} \) with \( R_b \) free per galaxy, and hydro+soft=kr1(2p+g) disciplines it the same way, \( R_b = \kappa r_1 \).
== outer region == galaxies: 143 (global kappa: shift=1.8, soft=1.3)
model median chi2/dof <1 <2 beats hydrogen
hydrogen(2p) 0.67 59% 77%
Kepler(1p) 23.85 0% 1% 0/143 (0%)
MOND/RAR(1p) 1.94 35% 52% 33/143 (23%)
NFW(2p+bar) 0.34 73% 85% 104/143 (73%)
Burkert(2p+bar) 0.28 76% 84% 99/143 (69%)
hydro+c1(3p) 0.37 76% 87% 120/143 (84%)
hydro+dr(3p) 0.32 76% 87% 122/143 (85%)
hydro+dr=kr1(2p+g) 0.40 69% 82% 98/143 (69%)
hydro+soft(3p) 0.40 73% 82% 90/143 (63%)
hydro+soft=kr1(2p+g) 0.47 67% 80% 105/143 (73%)
== full region == galaxies: 143 (global kappa: shift=0.5, soft=0.5)
model median chi2/dof <1 <2 beats hydrogen
hydrogen(2p) 11.98 7% 15%
Kepler(1p) 313.29 0% 0% 0/143 (0%)
MOND/RAR(1p) 3.51 15% 31% 114/143 (80%)
NFW(2p+bar) 1.85 36% 54% 133/143 (93%)
Burkert(2p+bar) 1.02 50% 64% 136/143 (95%)
hydro+c1(3p) 4.09 18% 30% 125/143 (87%)
hydro+dr(3p) 3.93 13% 29% 125/143 (87%)
hydro+dr=kr1(2p+g) 7.27 7% 18% 88/143 (62%)
hydro+soft(3p) 4.30 11% 22% 84/143 (59%)
hydro+soft=kr1(2p+g) 7.27 8% 17% 86/143 (60%)
The headline of the new rows: on the outer region, the third biharmonic slot brings the hydrogen to halo parity without any baryonic input. The shift at 0.32 and the constant mode at 0.37 straddle NFW's 0.34 and sit beside Burkert's 0.28, while consuming only \( (r, v, \sigma_v) \) where the halos consume the full photometric mass models; each beats the two-parameter hydrogen in about 85 percent of the zoo. On full curves the third slot buys a factor of three, from 11.98 to 3.93, overtaking MOND, while the halos keep the inner regions, exactly as the superposition caveat predicts. That the shift row and the constant row score as near-twins in both regions is not an accident but a theorem, Proposition 8 of the bar section: they are the same signal in two costumes, and only an external referee can split them.
Read it slowly, both ways. In the outer region, where the dark matter mystery lives, the two-parameter hydrogen sits comfortably inside the error bars (median 0.67) and beats the one-parameter MOND fit on the median, but the halo models, given two parameters and the baryonic photometry, fit better still in about three quarters of the galaxies. On full curves the single hydrogen fails, exactly as the previous post's own honesty section predicted: a real galaxy is a superposition of sources and the single-source fit is only the far field. Three caveats that cut in different directions, all on the table. The halo fits consume the baryonic data which the hydrogen does not even look at, so the comparison is not a clean parameter count in either direction. Halo parameters sprawl over decades as galaxy mass varies while \( \beta \) stays inside 0.64 dex, but the right response to that is not to declare \( \beta \) a constant of nature; it is to measure what \( \beta \) actually depends on, which we do next. And the models diverge where data runs out: beyond the last measured point the hydrogen curve keeps rising as \( \sqrt{\beta r} \) while every halo bends down toward Kepler, so the falsifiable discriminant of the previous post stands, sharpened: it is now the only regime where the hydrogen is distinguishable from a Burkert halo at these error bars.
The law of \( \beta \). The fits determine each galaxy's \( \beta \) to a median precision of 8 percent, while the population spreads by 0.27 dex, so \( \beta \) is a genuine per-galaxy property, well measured and genuinely varying, not a universal constant. Regressing it against the observables available in the mass-models table itself (code in Appendix C):
log beta vs log(Vout^2/Rmax) slope +0.76 r = +0.85 scatter 0.14 dex
best 2-var fit: log beta = 1.67 log Vout - 1.07 log Rmax + 0.70 (0.11 dex)
pure dimensional law beta = Vout^2/Rmax: offset -0.05 dex, scatter 0.16 dex
The free exponents land on the dimensional law, \( \beta \approx 0.9 V_{\mathrm{out}}^2 / R_{\mathrm{max}} \) to 0.16 dex, which says in words: \( \beta \) is the observed centripetal acceleration at the last measured point of each galaxy. The near-universality of the previous post dissolves into the narrow dynamic range of \( g_{\mathrm{obs}}(R_{\mathrm{max}}) \) across the SPARC zoo, which is the radial acceleration relation15 wearing our notation. Two honest caveats, cutting in opposite directions. For perfectly flat outer curves this relation is partly forced by the fit geometry itself, so the real information lives in the 0.16 dex of residuals where the outer slopes hide. And it cuts against a purely environmental reading of \( \beta \) in the superposition picture: the ambient growth mode around each galaxy tracks that galaxy's own \( V^2/R \), which pushes toward the sources of \( \beta \) being the galaxies themselves, exactly the kind of self-sourcing that Proposition 6 exhibits in the electromagnetic sector and that the mass sector should now be hunted for.
And one promised correction of our own numbers, because the deal is that wounds get shown: the previous post claimed the curve \( \alpha/r + \beta r \) is flat "within about seven percent over a full decade of radius"; the true numbers are seven percent over a factor of 2.9 in radius, and 32 percent over a full decade. The flatness-from-geometry point survives, the arithmetic did not. The cosmological corrections have moved to the dark energy section, where they belong.
The galaxy bar
Nearly every disk galaxy, watched long enough, shows a bar, and every section above kept tripping over the same inner region where the single hydrogen fails. This section takes one audacious reading of that coincidence seriously enough to give it a real test: the bar marks the region where the collapse from the higher dimensional metric to effective 3+1 is incomplete. Light crossing that region is refracted by the not-yet-collapsed geometry, so the radii we assign to matter inside it are not measurements at all, while every radius outside it carries one common offset, \( \delta r = r_{\mathrm{actual}} - r_{\mathrm{measured}} \), accumulated across the region between the source and the center we measure from. The gas streaming, shocks and dust lanes that real bars display are, on this reading, matter negotiating a physical dimensional interface, not evidence against an optical component. Before the data get a vote, two propositions pin down what the theory itself has to say.
Proposition 7 (the intrinsic radius). The two parameter far field owns exactly one length, \( r_1 = \sqrt{\alpha/\beta} \), and it is distinguished three times over: \( H(r_1) = 0 \) exactly, so \( g_{uu} = 1 \) and the null-pair block of the metric is Minkowski-normalized at that radius and only there; the curve \( v^2 = \alpha/r + \beta r \) attains its minimum there; and it is the fixed point of the symmetry \( r \rightarrow r_1^2/r \).
Proof. With \( \alpha = c_0\dot u^2 \) and \( \beta = -c_2\dot u^2 \), the potential of the galactic band is \( H(r) = (\alpha/r - \beta r)/\dot u^2 \), which vanishes iff \( r^2 = \alpha/\beta \). The other two claims are one derivative and one substitution. \( \blacksquare \)
So if any radius deserves to anchor a transition region, it is \( r_1 \): inside it \( H > 0 \) and the \( g_{uu} \) block is inflated, outside it \( |H| \) stays small across the flat decade. But the honest fits below keep \( r_1 \) out of the machinery entirely and let it reappear, or not, in the correlations.
Proposition 8 (a radius shift is the constant mode, at first order). For \( |\delta r| \ll r \),
\[ \frac{\alpha}{r+\delta r} + \beta (r+\delta r) = \frac{\alpha}{r} + \beta \delta r + \beta r - \frac{\alpha \delta r}{r^2} + O(\delta r^2) \]
so a common shift of all measured radii is, to first order, exactly the constant biharmonic mode \( c_1 = \beta \delta r \) that the two parameter fit has been discarding, plus a locked \( -\alpha \delta r/r^2 \) correction carrying the same \( \delta r \).
Proof. Taylor expansion of both terms. \( \blacksquare \)
This proposition is a warning label, and the data honor it: fitting the shift and fitting a free constant on the same outer windows gives medians 0.32 against 0.37, statistically indistinguishable in 78 percent of the galaxies, with the measured constants and shifts obeying the predicted identity in sign and size, median \( c_1 = +5486\ (\mathrm{km/s})^2 \) against median \( \beta \delta r = +3093\ (\mathrm{km/s})^2 \), the factor 1.8 gap owned by the locked correction and higher orders (code and logs: Appendix G). Whatever the fits below achieve is therefore one measured improvement with two competing readings, a refraction shift or a vacuum mode, and no amount of rotation-curve fitting can split them. Only an external referee can, and this section ends by naming one.
The prescription, exactly as the hypothesis demands. Per galaxy: all three of \( \alpha, \beta, \delta r \) free; the points inside the bar region excluded as non-measurements; and the bar radius itself defined by the smallest-cut criterion: \( R_{\mathrm{bar}} \) is the smallest radius such that the shifted hydrogen fits every retained point \( r > R_{\mathrm{bar}} \) at \( \chi^2/\mathrm{dof} \le 1.5 \). The smallness demand is the discipline, since excluding data always helps; the model is forced to keep every point it possibly can. Nothing tied to \( r_1 \), no universal coupling, no baryonic input anywhere. The verdict (code: Appendix G):
no cut needed (R_bar = 0): 27 galaxies (19%)
points excluded: median 14% (16-84%: 0% - 50%)
R_bar (cut>0): median 2.66 kpc (16-84%: 1.08 - 9.84) median R_bar/Rmax = 0.19
chi2/dof on retained: median 1.11 on ALL points incl. cut: median 6.82
delta_r: 82% positive median +10.9 kpc
log R_bar vs log Rmax slope +1.07 r = +0.78 log R_bar vs log r1 slope +0.33 r = +0.54
robustness (threshold 1.0): 13% no-cut, median 22% excluded, median R_bar 3.13 kpc
Read with care, because the encouraging parts and the definitional parts sit close together. The model needs remarkably little forgiveness: a fifth of the galaxies require no bar region at all, the median galaxy surrenders one point in seven, and what remains fits at halo-class quality, 1.11, still using zero baryonic information. The retained-region quality is partly by construction, since the threshold defines retention; the honest content is the smallness of the exclusions, not the goodness of the remainder. And then the result nothing in the fit was told to produce: the required exclusion radius lands squarely in the observed bar range, median 2.66 kpc at \( 0.19 R_{\mathrm{max}} \), against real S4G bar semi-lengths of typically one to five kiloparsecs2021, and it scales with galaxy size at slope 1.07, which sounds like a sampling artifact until one recalls that measured bars also scale with disk size. An exclusion radius derived purely from rotation-curve misfit reproduces the two most basic facts of bar phenomenology. The sign of the shift is the one the mechanism predicts, positive in 82 percent of the zoo, and its magnitude is the standing problem: a median of eleven kiloparsecs implies true galaxies nearly twice their measured size outside the bar, a claim that collides with every independent size and distance measure until an actual imaging computation through the metric, null geodesics traced across the trap region, says otherwise. That computation, not another fit family, is the theoretical debt of this section.
The referee. The prescription outputs a per-galaxy, zero-free-constant prediction: \( R_{\mathrm{bar}} \), in kiloparsecs, for each of the 143 galaxies, to be regressed directly against measured bar semi-lengths from the S4G morphology catalogue20, its deprojected companion21, and the classic independent set22. One cross-match exists so far: ESO116-G012, fitted \( R_{\mathrm{bar}} = 3.9 \) kpc against a measured \( 24.0'' \times 13.0\ \mathrm{Mpc} = 1.5 \) kpc, off by a factor 2.6 on the catalogue's lowest-quality flag, an edge-on galaxy with no deprojection: zero statistical weight, pipeline proven, verdict pending the NGC and UGC pages. Two sharper observables stand behind the catalogue test. Measured bars rotate, with pattern speeds obtained by the Tremaine-Weinberg method23; a transition region anchored at a static radius does not, unless it is carried by the theory's own exact \( m=2 \) growth modes, \( (x^2 - y^2) \) and \( r^2(x^2 - y^2) \), verified biharmonic in Appendix G, which produce genuinely barred, genuinely dynamical potentials and could give the interface its rotation. The optical reading and the dynamical reading may yet be one object; the pattern speed is where they must either merge or part.
Facing the rivals of the bar. The standard theory of bars is a success story, and we compare against it with respect. Bars grow spontaneously by a global instability of self-gravitating disks30, so violently that Ostriker and Peebles proposed massive halos partly to suppress them31; they then evolve secularly, trading angular momentum with the halo3334, and dynamical friction against a dense halo slows them down32. That last link is the standard picture's open bruise: measured pattern speeds come out fast, corotation just beyond the bar end, in most galaxies, which by the friction argument demands low central halo densities, in tension with the cuspy halos the same paradigm predicts. The fast-bar problem. Note what this theory does to it: no halo, no friction, so fast bars are not a tension here, they are the default. One comparative point, scored cheaply. Against that, the instability picture holds two cards we currently cannot match: it explains why bars rotate at all, since an instability pattern is born with a pattern speed, while our interface must borrow rotation from the \( m=2 \) modes; and it naturally accommodates the observed decline of the bar fraction toward high redshift35 through instability growth times, while a collapse interface should exist whenever the trap does, at every epoch, a discriminant that presently leans against the interface reading. And one card is uniquely ours: no other theory of bars extracts the bar radius from the rotation curve alone. The prescription above recovered per-galaxy bar radii, in kiloparsecs, in the observed range and with the observed size scaling, from misfit data that knows nothing of photometric bars. If the catalogue regression confirms the match, that is a prediction the instability picture never made; if it refutes it, the interface reading dies by exactly the sword it chose.
Dark energy
The previous post's dark energy argument stood on two legs: the exact endpoint redshift of comoving observers, and the static-space luminosity distance \( d_L = d(1+z) \), which together gave \( q_0 = -1 \) exactly for \( c_3 = 0 \) and the observed \( q_0 \approx -0.55 \) at \( c_3/c_2^2 \approx 0.225 \)1. This post upgrades the first leg, corrects two numbers, and arms the second leg with falsifiers it did not have.
The upgrade: the redshift formula is no longer an assumption about test photons. By Proposition 4, genuine electromagnetic waves of the derived fourth-order theory ride the doubled null characteristic, so their frequencies redshift by the Killing charges, which is the previous post's formula, now a consequence of the field equations, with one previously unstated condition made explicit: the exact endpoint form holds for photons with \( p_u = 0 \), a choice the null constraint permits.
The corrections, because wounds get shown: matching the Hubble slope with \( \sigma = 1 \) gives \( c_2^{\mathrm{cosmo}} = H_0/c = 2.3\times 10^{-7}\ \mathrm{kpc}^{-1} \), not the previous post's \( 3.5\times 10^{-7} \), which silently corresponded to \( \sigma^2 = 1/2 \). The corrected ratio between the cosmological coupling and the median galactic \( |c_2\dot u^2| = 7.8\times 10^{-9}\ \mathrm{kpc}^{-1} \) is a factor of 30, not 45, bridged by the unmeasured drifts; the unification coincidence gets slightly better under correction.
And the new falsifiers, both born in the light section. First, the partner waves of Proposition 5 do not dilute as \( 1/r \), so if astrophysical sources excite them, the theory predicts violations of the Etherington distance duality19: the supernova side of the dark energy argument moves while the redshift side stays put, and existing distance-duality constraints bound the excitation directly. Second, the back-reaction hierarchy of Proposition 6 makes the electromagnetic growing modes gravitate with distinguishable textures, the \( b_2 \) partner sourcing a logarithmically growing potential \( H = b_2^2(\ln r - \tfrac12 P_2) \) on top of the source-free Hubble modes \( c_2 r \) and \( c_3 r^2 \): a logarithmic term in the Hubble relation is the most testable trace this theory could leave in cosmological data, and it is not a term the standard model of cosmology contains. The fully propagating, time-dependent analysis remains the announced sequel; these are its fixed goalposts.
Facing the rivals of dark energy. The dark matter section earned its keep by standing next to MOND and the halos; this section owes the same courtesy to the rivals of dark energy. Two propositions first, because the oldest static-redshift proposal died of exactly the diseases they rule out.
Proposition 9 (metric redshift dilates and does not blur). In the stationary metric, the redshift factor between two comoving observers is a ratio of Killing-energy projections, linear in the photon momentum and therefore identical for every frequency: the redshift is achromatic and involves no scattering. And because translation in \( t \) is an isometry, two signal features emitted a proper time \( \Delta\tau_s \) apart arrive \( \Delta\tau_o = (1+z)\Delta\tau_s \) apart: every waveform, a supernova light curve included, stretches by exactly \( 1+z \).
Proof. The frequency \( \omega = -g(U,p) \) is linear in \( p \) and the Killing charges \( p_t, p_u \) are conserved, so \( \omega_s/\omega_o \) depends only on the two observers. For the dilation: shift the whole emission history by \( \delta t \); stationarity maps solutions to solutions and endpoints to endpoints, so the arrival history shifts by the same \( \delta t \), and converting coordinate separations to proper time on each comoving worldline supplies the factor \( \omega_s/\omega_o = 1+z \). \( \blacksquare \)
Proposition 10 (static photometry: \( d_L \) derived, and the Tolman exponent). Photon energies redshift by \( 1+z \) and arrival rates dilate by \( 1+z \) (Proposition 9), while areas and angles in the static space are undistorted, so the flux from a source at distance \( d \) is \( L/[4\pi d^2(1+z)^2] \): the luminosity distance is \( d_L = d(1+z) \), exactly the law the previous post assumed1, now derived. And with \( d_A = d \), surface brightness scales as \( (1+z)^{-2} \), against \( (1+z)^{-4} \) in any expanding spacetime.
Proof. Two factors of \( 1+z \) from Proposition 9, none from geometry; the expanding exponent is the Etherington relation \( d_L = (1+z)^2 d_A \)19. \( \blacksquare \)
The scorecard, in the spirit of the dark matter table:
LCDM quintessence tired light this theory
redshift mechanism expansion expansion photon energy metric (Killing
loss en route charges)
SN time dilation yes yes no (falsified) yes (Prop 9)
image blurring none none scattering blur none (Prop 9)
distance duality exact exact violated violable (Prop 5)
Tolman SB exponent 4 4 ~1 2 + partner modes
CMB blackbody natural natural unexplained unaddressed
q0 today -0.55 (fit) w-dependent n/a 2 c3/c2^2 - 1 (fit)
new ingredients Lambda scalar field scattering law c2, c3 vacuum modes
Tired light24 was the historical static-redshift proposal, and the checklist above is the record of its execution: no time dilation, though dilation is measured cleanly in supernova light curves25 and spectra26; scattering blur that is not observed; no blackbody. This theory is not tired light, and the distinction is structural, not rhetorical: the redshift here is metric, so Proposition 9 passes the dilation and blurring tests the same way general relativity does. Worth saying loudly, because every static-universe proposal gets filed under Zwicky by default, and this one belongs in the gravitational-redshift drawer instead. Against \( \Lambda \)CDM2829 the honest scoreboard reads: we match the two headline numbers, the Hubble slope by construction and \( q_0 \) with the one ratio \( c_3/c_2^2 \); we add two falsifiers that \( \Lambda \)CDM forbids, the Etherington violation of Proposition 5 and the logarithmic Hubble texture of Proposition 6; and we lose, today, on two fronts stated without disguise. The Tolman surface-brightness test27 favors the expanding exponent once luminosity evolution is modeled, and our baseline exponent is 2. And the cosmic microwave background, blackbody and acoustic peaks both, is the largest dataset in cosmology and this theory has not addressed it at all. The one lever we own is Proposition 5 itself: partner-mode conversion moves amplitude between the decaying and non-decaying channels, so it moves the effective Tolman exponent, but in which direction depends on the sign of the conversion, which only the propagating sequel can compute. A dial exists; whether it turns the right way is not yet ours to claim.
Open wounds
Carried forward and new, so nobody accuses us of hiding them. The inner galaxies still need the superposition fit, and the halos will be the benchmark to beat there, not Kepler. The sign flip of \( c_2 \) between the local and global story stands. The electric and propagating sector still needs the \( u \)-dependent analysis; the lemmas \( {\Gamma^u}_{ab} = 0 \) and \( \nabla k = 0 \) survive \( u \)-dependence, so Theorem 1's machinery is ready for it, and that is the announced sequel. The \( SU(3) \) slot is holonomy, not dynamics, and the \( \mathbf{3}\oplus\mathbf{5} \) skeleton of Proposition 2 is \( SO(3) \) group theory, not \( SU(3) \) dynamics. The law of \( \beta \) is measured, not derived: the theory now owes the data a demonstration that superposed hydrogens produce \( \beta \approx g_{\mathrm{obs}}(R_{\mathrm{max}}) \) with the observed 0.16 dex of scatter, and until that fit exists the law is a constraint, not a success. The back-reaction hierarchy of Proposition 6 needs matching to cosmology: which geometric currents are actually excited, and whether the logarithmic texture survives in the Hubble diagram, is an open hunt with an open falsifier. And clause (i) of Theorem 1 leans on the Riemannian signature of the transverse block, so it is a statement about this null-pair ansatz, not about every conceivable fibration. And the bar prescription owes three debts of its own: the fitted \( \delta r \) implies galaxies nearly twice their measured size, indefensible until null geodesics traced through the trap region either produce it or kill it; the smallest-cut criterion makes the retained-region fit quality partly definitional, so only the smallness of the exclusions counts as evidence; and the external referee has so far seen a single lowest-quality data point, which is to say nothing at all. Until the catalogue votes, Proposition 8 keeps the dark matter table's best new rows officially ambiguous between a shift and a vacuum mode.The dark energy comparison adds two wounds larger than any before them: the Tolman exponent of the static reading is 2 against measurements favoring the expanding 4 once evolution is modeled, with the partner-mode dial of Proposition 5 pointing in an as-yet-uncomputed direction; and the cosmic microwave background is the largest dataset in cosmology and this theory has not addressed it at all. And the bar comparison concedes the instability picture two cards, the birth of pattern speeds and the declining bar fraction at high redshift, that the interface reading has not earned.
Conclusion
We ran the Kaluza method on the Ricci-squared theory without training wheels, an arbitrary base metric and all, and the theory answered with more structure than we asked for. The base cannot be arbitrary: it is forced Ricci-flat, which decouples it, protects the "charge is geometry" reading, vindicates the flat ansatz of the previous post as the simplest member of the only allowed family, and selects Calabi–Yau geometry for the curved internal option. Electromagnetism came out derived: gauge invariance as coordinate freedom, Maxwell as the exact decaying sector, the extra modes as growing partners in perfect analogy with the dark sector, the linearized kernel of the previous post explained term by term, and the machine-discovered spin coefficients reduced to a paper-and-pencil identity about dipole self-energy. The effective stress tensor became a formula instead of a metaphor, conserved by Bianchi, with matter flux and electric current revealed as the same object. The trapped sector grew to eight dimensions under one theorem instead of one machine hunt, forgetting its dimensions at one rate on two incommensurate log-clocks, and carrying an \( \mathfrak{su}(3) \)-adjoint's worth of directions filed, honestly, under suggestive. Light itself joined the theory: charge conservation held as an identity and killed our favorite redshift mechanism, the endpoint redshift formula got re-derived for genuine waves, the partner modes broke the inverse square law instead of the frequency, and the growing electromagnetic modes turned out to gravitate in a clean hierarchy, from the dipole's bounded self-energy through the partner's logarithm to the uniform current's quartic. And the data section grew up three times: the hydrogen now stands honestly among MOND and the halos, winning some regimes, losing others; \( \beta \) confessed its law, the observed acceleration at each galaxy's edge, a measured per-galaxy property that the theory must now derive; and the bar hypothesis received exactly the test it asked for, free shift, explicit cut, smallest-cut discipline, and answered by excluding a median of one point in seven while predicting per-galaxy bar radii, in kiloparsecs, that land in the observed range and now await the catalogue referee. The equations touched the sky last time; this time they shook hands with a century of literature and survived the introduction.
Appendix A: the SageMath verifier
The same fraction-field trick as the previous post's engine: everything in our charts is a rational function, so all tensor algebra runs over exact fraction fields of polynomial rings where every element is automatically canonical and the zero test is == 0, no simplify heuristics anywhere. The script verifies every component identity behind Theorem 1, twice: on a generic curved, deliberately non-Ricci-flat base with non-solution profiles (so nothing vanishes by accident), and on Euclidean Schwarzschild as a curved Ricci-flat base. Its output, ten True lines, is shown after the listing. (One subtlety the first version got wrong, caught by a reader's stack trace: the derivative helper for cyclic coordinates must return the field's own zero, x*0, not a bare integer zero, or later derivatives choke on an element that has fallen out of the fraction field.) The logically identical checks were also executed in an independent CAS, SymPy, with the logs shown in Appendix B.
def geometry(FF, gens, g):
"""Christoffels, Riemann, Ricci, R, Ric^2 and E_ab over a fraction field.
gens[k] is the polynomial generator of coordinate k, or None for cyclic coordinates."""
dim = len(g); D = range(dim); zero = FF(0)
def d(f, k):
return f.derivative(gens[k]) if gens[k] is not None else zero
gm = matrix(FF, dim, dim, [g[i][j] for i in D for j in D]).inverse()
gi = [[gm[i,j] for j in D] for i in D]
Gam = [[[zero]*dim for _ in D] for _ in D]
for b in D:
for c in range(b, dim):
col = [d(g[e][c],b)+d(g[e][b],c)-d(g[b][c],e) for e in D]
for a in D:
v = sum(gi[a][e]*col[e] for e in D)/2
Gam[a][b][c]=v; Gam[a][c][b]=v
Riem = [[[[zero]*dim for _ in D] for _ in D] for _ in D]
for a in D:
for b in D:
for c in D:
for e in range(c+1, dim):
v = d(Gam[a][e][b],c)-d(Gam[a][c][b],e)+sum(
Gam[a][c][f]*Gam[f][e][b]-Gam[a][e][f]*Gam[f][c][b] for f in D)
Riem[a][b][c][e]=v; Riem[a][b][e][c]=-v
Ric = [[sum(Riem[a][b][a][e] for a in D) for e in D] for b in D]
Rs = sum(gi[b][e]*Ric[b][e] for b in D for e in D)
Cd = [[[ d(Ric[a][b],c)-sum(Gam[f][c][a]*Ric[f][b]+Gam[f][c][b]*Ric[a][f]
for f in D) for b in D] for a in D] for c in D]
BoxRic = [[zero]*dim for _ in D]
for a in D:
for b in range(a, dim):
v = zero
for c in D:
for e in D:
if gi[c][e]==zero: continue
Dt = d(Cd[c][a][b],e)-sum(Gam[f][e][c]*Cd[f][a][b]
+Gam[f][e][a]*Cd[c][f][b]+Gam[f][e][b]*Cd[c][a][f] for f in D)
v += gi[c][e]*Dt
BoxRic[a][b]=v; BoxRic[b][a]=v
dRs = [d(Rs,a) for a in D]
Hess = [[d(dRs[b],a)-sum(Gam[f][a][b]*dRs[f] for f in D) for b in D] for a in D]
BoxRs= sum(gi[a][b]*Hess[a][b] for a in D for b in D)
RicU = [[sum(gi[c][e]*gi[dd][f]*Ric[e][f] for e in D for f in D) for dd in D] for c in D]
Rl = [[[[sum(g[a][f]*Riem[f][b][c][e] for f in D) for e in D] for c in D] for b in D] for a in D]
RR = [[sum(Rl[a][c][b][e]*RicU[c][e] for c in D for e in D) for b in D] for a in D]
RicSq = sum(Ric[a][b]*RicU[a][b] for a in D for b in D)
E = [[BoxRic[i][j]+(BoxRs/2)*g[i][j]-Hess[i][j]+2*RR[i][j]-(RicSq/2)*g[i][j]
for j in D] for i in D]
return {'gi':gi, 'Gam':Gam, 'Ric':Ric, 'Rs':Rs, 'RicSq':RicSq, 'E':E}
def hat_lap(f, gens, gi, Gam, dim):
def d(x, k): return x.derivative(gens[k]) if gens[k] is not None else x*0
return sum(gi[i][j]*(d(d(f,j),i) - sum(Gam[k][i][j]*d(f,k) for k in range(dim)))
for i in range(dim) for j in range(dim))
def rough_box(V, gens, gi, Gam, dim):
def d(x, k): return x.derivative(gens[k]) if gens[k] is not None else x*0
dV = [[d(V[i],l) - sum(Gam[m][l][i]*V[m] for m in range(dim))
for i in range(dim)] for l in range(dim)]
return [sum(gi[k][l]*(d(dV[l][i],k)
- sum(Gam[m][k][l]*dV[m][i] + Gam[m][k][i]*dV[l][m] for m in range(dim)))
for k in range(dim) for l in range(dim)) for i in range(dim)]
# ---- Test 1: generic curved base (NOT Ricci-flat), non-solution profiles, 5D ----
PR = PolynomialRing(QQ, ['x','y','z']); FF = PR.fraction_field()
x, y, z = [FF(v) for v in PR.gens()]; X = list(PR.gens())
H = x*y - z^2 + x*z
A = [y*z + x^2, x*z - y^2, x*y + z^2]
h = [[FF(1),FF(0),FF(0)],[FF(0),1+x^2,FF(0)],[FF(0),FF(0),1+y^2]]
g5 = [[FF(0)]*5 for _ in range(5)]
g5[0][1]=FF(1); g5[1][0]=FF(1); g5[1][1]=1+2*H
for i in range(3):
g5[1][2+i]=A[i]; g5[2+i][1]=A[i]
for j in range(3): g5[2+i][2+j]=h[i][j]
full = geometry(FF, [None,None]+X, g5)
base = geometry(FF, X, h)
gi3, Gam3 = base['gi'], base['Gam']
d3 = lambda f,i: f.derivative(X[i])
E, Ric = full['E'], full['Ric']
print("E_tt and E_ti == 0:", E[0][0]==0 and all(E[0][2+i]==0 for i in range(3)))
print("E_tu == (hatLap R - Ric^2)/2 of the base:",
E[0][1] == (hat_lap(base['Rs'], X, gi3, Gam3, 3) - base['RicSq'])/2)
print("E_ij == E_ij[h], the self-reduction:",
all(E[2+i][2+j]==base['E'][i][j] for i in range(3) for j in range(3)))
print("R_ij == Ricci[h], no F or H leakage:",
all(Ric[2+i][2+j]==base['Ric'][i][j] for i in range(3) for j in range(3)))
Fm = [[d3(A[j],i)-d3(A[i],j) for j in range(3)] for i in range(3)]
cdF = [[[d3(Fm[k][i],l)-sum(Gam3[m][l][k]*Fm[m][i]+Gam3[m][l][i]*Fm[k][m]
for m in range(3)) for i in range(3)] for k in range(3)] for l in range(3)]
W = [sum(gi3[k][l]*cdF[l][i][k] for k in range(3) for l in range(3))/2 for i in range(3)]
F2 = sum(gi3[i][a]*gi3[j][b]*Fm[i][j]*Fm[a][b]
for i in range(3) for j in range(3) for a in range(3) for b in range(3))
print("R_ui == (1/2) hat-div F:", all(Ric[1][2+i]==W[i] for i in range(3)))
print("R_uu == -hatLap H + F^2/4:",
Ric[1][1] == -hat_lap(H, X, gi3, Gam3, 3) + F2/4)
# ---- Test 2: Euclidean Schwarzschild base (curved, Ricci-flat), 6D ----
PR2 = PolynomialRing(QQ, ['r','x','M']); F6 = PR2.fraction_field()
r, xc, M = [F6(v) for v in PR2.gens()]
fS = 1 - 2*M/r
h4 = [[fS,F6(0),F6(0),F6(0)],[F6(0),1/fS,F6(0),F6(0)],
[F6(0),F6(0),r^2/(1-xc^2),F6(0)],[F6(0),F6(0),F6(0),r^2*(1-xc^2)]]
gens4 = [None, PR2.gen(0), PR2.gen(1), None] # (w, r, x, phi)
aP = r^3 + r; HP = r^2 + 1/r # deliberately NOT solutions
A4 = [aP, F6(0), F6(0), F6(0)]
g6 = [[F6(0)]*6 for _ in range(6)]
g6[0][1]=F6(1); g6[1][0]=F6(1); g6[1][1]=1+2*HP
for i in range(4):
g6[1][2+i]=A4[i]; g6[2+i][1]=A4[i]
for j in range(4): g6[2+i][2+j]=h4[i][j]
full6 = geometry(F6, [None,None]+gens4, g6)
base4 = geometry(F6, gens4, h4)
gi4, Gam4 = base4['gi'], base4['Gam']
def d4(f,k): return f.derivative(gens4[k]) if gens4[k] is not None else f*0
print("base Ricci-flat:", all(base4['Ric'][i][j]==0 for i in range(4) for j in range(4)))
E6, R6 = full6['E'], full6['Ric']
print("E_tu == 0 and E_ij == 0, base decoupled and solving itself:",
E6[0][1]==0 and all(E6[2+i][2+j]==0 for i in range(4) for j in range(4)))
Fm4 = [[d4(A4[j],i)-d4(A4[i],j) for j in range(4)] for i in range(4)]
cdF4 = [[[d4(Fm4[k][i],l)-sum(Gam4[m][l][k]*Fm4[m][i]+Gam4[m][l][i]*Fm4[k][m]
for m in range(4)) for i in range(4)] for k in range(4)] for l in range(4)]
W4 = [sum(gi4[k][l]*cdF4[l][i][k] for k in range(4) for l in range(4))/2 for i in range(4)]
J4 = [-2*W4[i] for i in range(4)]
F24 = sum(gi4[i][a]*gi4[j][b]*Fm4[i][j]*Fm4[a][b]
for i in range(4) for j in range(4) for a in range(4) for b in range(4))
lap = lambda f: hat_lap(f, gens4, gi4, Gam4, 4)
BW = rough_box(W4, gens4, gi4, Gam4, 4)
cdJ = [[d4(J4[k],i)-sum(Gam4[m][i][k]*J4[m] for m in range(4)) for k in range(4)]
for i in range(4)]
FdJ = sum(gi4[i][a]*gi4[k][b]*Fm4[a][b]*cdJ[i][k]
for i in range(4) for k in range(4) for a in range(4) for b in range(4))
J2 = sum(gi4[i][j]*J4[i]*J4[j] for i in range(4) for j in range(4))
print("clause (ii): E_ui == roughBox W_i:",
all(E6[1][2+i]==BW[i] for i in range(4)))
print("clause (iii): E_uu == hatLap(-hatLap H + F^2/4) + J^2/2 + F.hatDel J:",
E6[1][1] == lap(-lap(HP) + F24/4) + J2/2 + FdJ)
Its output, executed in Sage 10:
E_tt and E_ti == 0: True
E_tu == (hatLap R - Ric^2)/2 of the base: True
E_ij == E_ij[h], the self-reduction: True
R_ij == Ricci[h], no F or H leakage: True
R_ui == (1/2) hat-div F: True
R_uu == -hatLap H + F^2/4: True
base Ricci-flat: True
E_tu == 0 and E_ij == 0, base decoupled and solving itself: True
clause (ii): E_ui == roughBox W_i: True
clause (iii): E_uu == hatLap(-hatLap H + F^2/4) + J^2/2 + F.hatDel J: True
Appendix B: the SymPy cross-check
An independent CAS, a different codebase, a different canonicalization strategy, and this time fully self-contained: the engine (iszero, curv, E_tensor) is included in the script itself, so nothing needs to be copied from the previous post. Save it with a .py extension and run python3 appendix_b.py (SymPy required). The extension is the whole trick: the Sage preparser only fires on .sage files, so even sage appendix_b.py executes it verbatim, while a .sage copy gets every literal rewritten and the zero tests break.
# appendix_b.py -- run with plain CPython: python3 appendix_b.py
# (requires sympy; do NOT pass this file through the Sage preparser)
import sympy as sp, time
# ---- The post's own E_tensor machinery (Appendix A, SymPy cross-check) ----
def iszero(e):
n, _ = sp.fraction(sp.cancel(sp.expand(e)))
return sp.expand(n) == 0
def curv(coords, g):
n = len(coords)
gm = sp.Matrix(g); gi = gm.inv()
d = lambda f,i: sp.diff(f, coords[i])
Gam = [[[sp.cancel(sum(gi[a,e]*(d(gm[e,c],b)+d(gm[e,b],c)-d(gm[b,c],e))
for e in range(n))/2) for c in range(n)] for b in range(n)] for a in range(n)]
Riem = [[[[sp.cancel(d(Gam[a][dd][b],c)-d(Gam[a][c][b],dd)
+ sum(Gam[a][c][e]*Gam[e][dd][b]-Gam[a][dd][e]*Gam[e][c][b] for e in range(n)))
for dd in range(n)] for c in range(n)] for b in range(n)] for a in range(n)]
Ric = [[sp.cancel(sum(Riem[a][b][a][dd] for a in range(n))) for dd in range(n)] for b in range(n)]
Rs = sp.cancel(sum(gi[b,dd]*Ric[b][dd] for b in range(n) for dd in range(n)))
return gm, gi, Gam, Riem, Ric, Rs
def E_tensor(coords, g):
n = len(coords)
gm, gi, Gam, Riem, Ric, Rs = curv(coords, g)
d = lambda f,i: sp.diff(f, coords[i])
C = [[[sp.cancel(d(Ric[a][b],c) - sum(Gam[e][c][a]*Ric[e][b]+Gam[e][c][b]*Ric[a][e]
for e in range(n))) for b in range(n)] for a in range(n)] for c in range(n)]
BoxRic = [[sp.cancel(sum(gi[c,dd]*(d(C[c][a][b],dd)
- sum(Gam[e][dd][c]*C[e][a][b]+Gam[e][dd][a]*C[c][e][b]+Gam[e][dd][b]*C[c][a][e]
for e in range(n)))
for c in range(n) for dd in range(n))) for b in range(n)] for a in range(n)]
dRs = [d(Rs,a) for a in range(n)]
Hess = [[sp.cancel(d(dRs[b],a)-sum(Gam[e][a][b]*dRs[e] for e in range(n)))
for b in range(n)] for a in range(n)]
BoxRs = sp.cancel(sum(gi[a,b]*Hess[a][b] for a in range(n) for b in range(n)))
RicU = [[sp.cancel(sum(gi[c,e]*gi[dd,f]*Ric[e][f] for e in range(n) for f in range(n)))
for dd in range(n)] for c in range(n)]
Rlow = [[[[sp.cancel(sum(gm[a,f]*Riem[f][b][c][e] for f in range(n)))
for e in range(n)] for c in range(n)] for b in range(n)] for a in range(n)]
RR = [[sp.cancel(sum(Rlow[a][c][b][e]*RicU[c][e] for c in range(n) for e in range(n)))
for b in range(n)] for a in range(n)]
RicSq = sp.cancel(sum(Ric[a][b]*RicU[a][b] for a in range(n) for b in range(n)))
return [[sp.cancel(BoxRic[a][b] + BoxRs/2*gm[a,b] - Hess[a][b]
+ 2*RR[a][b] - RicSq/2*gm[a,b]) for b in range(n)] for a in range(n)], Ric, Rs, RicSq, RR
# ============ Test 1: generic curved, NON-Ricci-flat base, 5D ============
t,u,x,y,z = sp.symbols('t u x y z'); X = [x,y,z]
# generic CURVED transverse metric (not Ricci-flat) + non-solution profiles
h = sp.diag(1, 1 + x**2, 1 + y**2)
H = x*y - z**2 + x*z
A = [y*z + x**2, x*z - y**2, x*y + z**2]
g = sp.zeros(5,5); g[0,1] = 1; g[1,0] = 1; g[1,1] = 1 + 2*H
for i in range(3):
g[1,2+i] = A[i]; g[2+i,1] = A[i]
for j in range(3): g[2+i,2+j] = h[i,j]
t0 = time.time()
E, Ric, Rs, RicSq, RR = E_tensor([t,u,x,y,z], g)
print("5D E_tensor done in %.0fs" % (time.time()-t0))
# transverse-only machinery: literally the same engine one dimension down
hE, hRic, hRs, hRicSq, hRR = E_tensor([x,y,z], h)
hgm, hgi, hGam, hRiem, _, _ = curv([x,y,z], h)
D = lambda f,i: sp.diff(f, X[i])
def hatLap(f): # Laplace-Beltrami via Hessian trace (all-rational, no sqrt)
Hs = [[D(D(f,j),i) - sum(hGam[k][i][j]*D(f,k) for k in range(3))
for j in range(3)] for i in range(3)]
return sp.cancel(sum(hgi[i,j]*Hs[i][j] for i in range(3) for j in range(3)))
# --- structural checks on the full 5D E ---
print("E_tt == 0:", iszero(E[0][0]))
print("E_ti == 0:", all(iszero(E[0][2+i]) for i in range(3)))
print("E_tu == (LapR - Ric^2)/2 of h:", iszero(E[0][1] - (hatLap(hRs) - hRicSq)/2))
print("E_ij == transverse self-reduction E_ij[h]:",
all(iszero(E[2+i][2+j] - hE[i][j]) for i in range(3) for j in range(3)))
print("R_ij == Ricci[h] (no F,H leakage):",
all(iszero(Ric[2+i][2+j] - hRic[i][j]) for i in range(3) for j in range(3)))
# --- Ricci closed forms with covariant operators ---
F = [[D(A[j],i) - D(A[i],j) for j in range(3)] for i in range(3)]
covdF = [[[sp.cancel(D(F[k][i],l) - sum(hGam[m][l][k]*F[m][i] + hGam[m][l][i]*F[k][m]
for m in range(3))) for i in range(3)] for k in range(3)] for l in range(3)]
Wc = [sp.cancel(sum(hgi[k,l]*covdF[l][i][k] for k in range(3) for l in range(3))/2)
for i in range(3)] # (1/2) nabla^k F_{ik}
F2 = sp.cancel(sum(hgi[i,a]*hgi[j,b]*F[i][j]*F[a][b]
for i in range(3) for j in range(3) for a in range(3) for b in range(3)))
print("R_ui == (1/2) hat-div F:", all(iszero(Ric[1][2+i] - Wc[i]) for i in range(3)))
print("R_uu == -hatLap H + F^2/4:", iszero(Ric[1][1] - (-hatLap(H) + F2/4)))
# ============ Test 2: Euclidean Schwarzschild base, 6D ============
t,u,w,r,x,ph,M = sp.symbols('t u w r x phi M')
XT = [w, r, x, ph] # transverse coordinates
f = 1 - 2*M/r
# Euclidean Schwarzschild: curved, Ricci-flat, all-rational in the x=cos(theta) chart
h = sp.diag(f, 1/f, r**2/(1-x**2), r**2*(1-x**2))
# deliberately NON-solution profiles (axisymmetric, so the run stays light)
a = r**3 + r # A = a(r) dw -> F_{rw} = a'
Hp = r**2 + 1/r # H(r)
A = [a, 0, 0, 0]
g = sp.zeros(6,6); g[0,1] = 1; g[1,0] = 1; g[1,1] = 1 + 2*Hp
for i in range(4):
g[1,2+i] = A[i]; g[2+i,1] = A[i]
for j in range(4): g[2+i,2+j] = h[i,j]
t0 = time.time()
E, Ric, Rs, RicSq, RR = E_tensor([t,u,w,r,x,ph], g)
print("6D E_tensor done in %.0fs" % (time.time()-t0))
hgm, hgi, hGam, hRiem, hRic, hRs = curv(XT, h)
print("transverse space Ricci-flat:", all(iszero(hRic[i][j]) for i in range(4) for j in range(4)))
D = lambda F_,i: sp.diff(F_, XT[i])
def hatLap(fn):
Hs = [[D(D(fn,j),i) - sum(hGam[k][i][j]*D(fn,k) for k in range(4))
for j in range(4)] for i in range(4)]
return sp.cancel(sum(hgi[i,j]*Hs[i][j] for i in range(4) for j in range(4)))
F = [[D(A[j],i) - D(A[i],j) for j in range(4)] for i in range(4)]
covdF = [[[sp.cancel(D(F[k][i],l) - sum(hGam[m][l][k]*F[m][i] + hGam[m][l][i]*F[k][m]
for m in range(4))) for i in range(4)] for k in range(4)] for l in range(4)]
Wc = [sp.cancel(sum(hgi[k,l]*covdF[l][i][k] for k in range(4) for l in range(4))/2)
for i in range(4)] # W_i = (1/2) nabla^k F_{ik}
Jc = [sp.cancel(-2*Wc[i]) for i in range(4)] # J_i = nabla^k F_{ki}
F2 = sp.cancel(sum(hgi[i,a]*hgi[j,b]*F[i][j]*F[a][b]
for i in range(4) for j in range(4) for a in range(4) for b in range(4)))
# rough Laplacian on a transverse covector
def roughBox(V):
dV = [[sp.cancel(D(V[i],l) - sum(hGam[m][l][i]*V[m] for m in range(4)))
for i in range(4)] for l in range(4)]
ddV = [[[sp.cancel(D(dV[l][i],k) - sum(hGam[m][k][l]*dV[m][i] + hGam[m][k][i]*dV[l][m]
for m in range(4))) for i in range(4)] for l in range(4)] for k in range(4)]
return [sp.cancel(sum(hgi[k,l]*ddV[k][l][i] for k in range(4) for l in range(4)))
for i in range(4)]
BW = roughBox(Wc)
covdJ = [[sp.cancel(D(Jc[k],i) - sum(hGam[m][i][k]*Jc[m] for m in range(4)))
for k in range(4)] for i in range(4)]
FdJ = sp.cancel(sum(hgi[i,a]*hgi[k,b]*F[a][b]*covdJ[i][k]
for i in range(4) for k in range(4) for a in range(4) for b in range(4)))
J2 = sp.cancel(sum(hgi[i,j]*Jc[i]*Jc[j] for i in range(4) for j in range(4)))
print("R_ui == (1/2) hat-div F:", all(iszero(Ric[1][2+i] - Wc[i]) for i in range(4)))
print("R_uu == -hatLap H + F^2/4:", iszero(Ric[1][1] - (-hatLap(Hp) + F2/4)))
print("E_tu == 0 and E_ij == 0 (base decoupled & solves itself):",
iszero(E[0][1]) and all(iszero(E[2+i][2+j]) for i in range(4) for j in range(4)))
print("E_ui == roughBox W_i (curved biharmonic Maxwell):",
all(iszero(E[1][2+i] - BW[i]) for i in range(4)))
print("E_uu == hatLap(-hatLap H + F^2/4) + J^2/2 + F.hatDel J:",
iszero(E[1][1] - (hatLap(sp.cancel(-hatLap(Hp) + F2/4)) + J2/2 + FdJ)))
Its output on our machine, verbatim:
5D E_tensor done in 17s
E_tt == 0: True
E_ti == 0: True
E_tu == (LapR - Ric^2)/2 of h: True
E_ij == transverse self-reduction E_ij[h]: True
R_ij == Ricci[h] (no F,H leakage): True
R_ui == (1/2) hat-div F: True
R_uu == -hatLap H + F^2/4: True
6D E_tensor done in 1s
transverse space Ricci-flat: True
R_ui == (1/2) hat-div F: True
R_uu == -hatLap H + F^2/4: True
E_tu == 0 and E_ij == 0 (base decoupled & solves itself): True
E_ui == roughBox W_i (curved biharmonic Maxwell): True
E_uu == hatLap(-hatLap H + F^2/4) + J^2/2 + F.hatDel J: True
Note the timings: seventeen seconds on the curved non-Ricci-flat base, one second on Euclidean Schwarzschild, because on a Ricci-flat base the invariants collapse and the expression swell collapses with them. The theorem is visible in the profiler. And one recorded output kept as a trophy: the fitting run that caught the sign of clause (iii) during this work returned
residual = a1*J^2 + a2*F.dJ + a3*|dF|^2 + a4*F.LapF with: [{a1: 1/2, a2: 1 - 2*a4, a3: 0}]
the one-parameter family being the Bianchi identity in disguise, as explained in the main text.
Appendix C: the rivals, and the law of \( \beta \)
The comparison runs on the same MassModels_Lelli2016c.mrt table as the previous post (mirrored at the sparc-rotation-curves repository), with the identical outer-region selection and error weighting. MOND is fitted in its radial acceleration relation form with \( g_\dagger = 1.2\times 10^{-10}\ \mathrm{m/s^2} \) and one free disk mass-to-light ratio (bulge fixed at 0.7); NFW and Burkert each add two free halo parameters on top of the baryonic curves at the standard \( \Upsilon_d = 0.5 \).
import numpy as np
from scipy.optimize import least_squares, minimize_scalar
FN = "MassModels_Lelli2016c.mrt"
gal = {}
for line in open(FN, errors='ignore'):
p = line.split()
if len(p) < 8: continue
try: vals = [float(v) for v in p[1:8]]
except ValueError: continue
gal.setdefault(p[0], []).append(vals)
KPC = 3.0857e19
GDAG = 1.2e-10 / (1e6/KPC) # 1.2e-10 m/s^2 in (km/s)^2/kpc
UB = 0.7
def chi2(vm, vo, eo, k): return np.sum(((vm-vo)/eo)**2)/max(len(vo)-k, 1)
def fit_ours(r, vo, eo):
w = 1.0/(2*vo*eo)
c, *_ = np.linalg.lstsq(np.column_stack([1/r, r])*w[:,None], vo**2*w, rcond=None)
c = np.clip(c, 0, None)
return chi2(np.sqrt(np.clip(c[0]/r + c[1]*r, 0, None)), vo, eo, 2), c[1]
def fit_kepler(r, vo, eo):
w = 1.0/(2*vo*eo)
c, *_ = np.linalg.lstsq((1/r)[:,None]*w[:,None], vo**2*w, rcond=None)
return chi2(np.sqrt(np.clip(max(c[0],0)/r, 0, None)), vo, eo, 1)
def vbar2(ud, vg, vd, vb): return vg*np.abs(vg) + ud*vd**2 + UB*vb**2
def fit_mond(r, vo, eo, vg, vd, vb):
def c2(ud):
gbar = np.clip(vbar2(ud, vg, vd, vb)/r, 1e-8, None)
gobs = gbar/(1.0 - np.exp(-np.sqrt(gbar/GDAG)))
return chi2(np.sqrt(gobs*r), vo, eo, 1)
return minimize_scalar(c2, bounds=(0.05, 5.0), method='bounded').fun
def vnfw2(r, V200, c):
x = r/(V200/0.7)
return V200**2*(np.log(1+c*x)-c*x/(1+c*x))/(x*(np.log(1+c)-c/(1+c)))
def fit_halo(r, vo, eo, vg, vd, vb, vh2, starts, lo, hi):
vb2 = vbar2(0.5, vg, vd, vb)
def resid(q):
v2 = np.clip(vb2 + vh2(r, 10**q[0], 10**q[1]), 0, None)
return (np.sqrt(v2) - vo)/eo
best = np.inf
for q0 in starts:
try:
s = least_squares(resid, q0, bounds=(lo, hi))
best = min(best, np.sum(s.fun**2)/max(len(r)-2, 1))
except Exception: pass
return best
def vburk2(r, rho0, r0):
G = 4.301e-6; x = r/r0
M = 2*np.pi*rho0*r0**3*(0.5*np.log(1+x**2)+np.log(1+x)-np.arctan(x))
return G*M/r
def run(region):
rows = []
for name, data in gal.items():
a = np.array(data)
r,v,ev,vg,vd,vb = a[:,1],a[:,2],a[:,3],a[:,4],a[:,5],a[:,6]
ok = (r>0)&(v>0)&(ev>0)
r,v,ev,vg,vd,vb = r[ok],v[ok],ev[ok],vg[ok],vd[ok],vb[ok]
if len(r) < 8: continue
if region == "outer":
m = r >= 0.4*r.max()
if m.sum() < 5: m = np.argsort(np.argsort(r)) >= len(r)-5
r,v,ev,vg,vd,vb = r[m],v[m],ev[m],vg[m],vd[m],vb[m]
cO, beta = fit_ours(r, v, ev)
rows.append((name, cO, fit_kepler(r,v,ev), fit_mond(r,v,ev,vg,vd,vb),
fit_halo(r,v,ev,vg,vd,vb, vnfw2, ([1.9,.9],[2.3,.5],[1.5,1.2]),
[1.0,0.0],[3.0,1.7]),
fit_halo(r,v,ev,vg,vd,vb, vburk2, ([7.5,.5],[6.5,1.2],[8.5,0.]),
[4.0,-1.0],[10.5,2.5]), beta))
R = np.array([[t[1],t[2],t[3],t[4],t[5]] for t in rows])
names = ["hydrogen(2p)","Kepler(1p)","MOND/RAR(1p)","NFW(2p+bar)","Burkert(2p+bar)"]
print(f"== {region} region == galaxies: {len(rows)}")
for j, nm in enumerate(names):
b = "" if j==0 else "%d/%d (%.0f%%)"%((R[:,j]<R[:,0]).sum(),len(R),
100*(R[:,j]<R[:,0]).mean())
print("%-17s median %6.2f <1: %3.0f%% <2: %3.0f%% beats hydrogen: %s"
% (nm, np.median(R[:,j]), 100*(R[:,j]<1).mean(), 100*(R[:,j]<2).mean(), b))
return rows
rows = run("outer"); print(); run("full")
beta = np.array([t[6] for t in rows if t[6] > 1e-6])
q = np.percentile(beta, [10,50,90])
print("beta 10/50/90: %.0f/%.0f/%.0f (km/s)^2/kpc, scatter %.2f dex" %
(q[0],q[1],q[2],np.log10(q[2]/q[0])))
print("|c2 udot^2| median = %.2e /kpc" % (q[1]/2.998e5**2))
Its output on our machine is the table quoted in the main text, plus
beta 10/50/90: 357/697/1562 (km/s)^2/kpc, scatter 0.64 dex
|c2 udot^2| median = 7.75e-09 /kpc
identical to the previous post's values, which is the audit of our own pipeline passing before the rivals get their turn.
And the regression behind the law of \( \beta \), run on the same table, with the per-galaxy uncertainties coming from the weighted least-squares covariance and the galaxy observables built from the mass models themselves at the standard mass-to-light ratios:
import numpy as np
FN = "MassModels_Lelli2016c.mrt"
gal = {}
for line in open(FN, errors='ignore'):
p = line.split()
if len(p) < 8: continue
try: vals = [float(x) for x in p[1:8]]
except ValueError: continue
gal.setdefault(p[0], []).append(vals)
G = 4.301e-6 # kpc (km/s)^2 / Msun
rows = []
for name, data in gal.items():
a = np.array(data)
r,v,ev,vg,vd,vb = a[:,1],a[:,2],a[:,3],a[:,4],a[:,5],a[:,6]
ok = (r>0)&(v>0)&(ev>0)
r,v,ev,vg,vd,vb = r[ok],v[ok],ev[ok],vg[ok],vd[ok],vb[ok]
if len(r) < 8: continue
m = r >= 0.4*r.max()
if m.sum() < 5: m = np.argsort(np.argsort(r)) >= len(r)-5
ro,vo,eo = r[m],v[m],ev[m]
w = 1.0/(2*vo*eo)
A = np.column_stack([1/ro, ro]); Aw = A*w[:,None]
c, *_ = np.linalg.lstsq(Aw, vo**2*w, rcond=None)
cov = np.linalg.inv(Aw.T @ Aw)
beta, ebeta = c[1], np.sqrt(cov[1,1])
if beta <= 0: continue
# galaxy properties from the same table (M/L = 0.5 disk, 0.7 bulge)
Rmax = r.max()
Vout = np.average(vo, weights=1/eo**2) # outer velocity
vbar2 = vg*np.abs(vg) + 0.5*vd**2 + 0.7*vb**2
Mbar = np.max(np.clip(vbar2, 0, None)*r)/G # enclosed baryonic mass
rows.append((name, beta, ebeta, Rmax, Vout, Mbar))
name = [x[0] for x in rows]
B = np.log10([x[1] for x in rows]); eB = np.array([x[2]/x[1]/np.log(10) for x in rows])
lR = np.log10([x[3] for x in rows]); lV = np.log10([x[4] for x in rows])
lM = np.log10([x[5] for x in rows])
lVR = np.log10(np.array([x[4] for x in rows])**2/np.array([x[3] for x in rows]))
def reg(xn, x):
p = np.polyfit(x, B, 1)
res = B - np.polyval(p, x)
r = np.corrcoef(x, B)[0,1]
print("log beta vs %-14s slope %+6.2f r = %+5.2f residual scatter %.2f dex"
% (xn, p[0], r, np.std(res)))
return np.std(res)
print("galaxies with beta > 0:", len(rows))
print("raw scatter of log beta: %.2f dex (std), %.2f dex (10-90)" %
(np.std(B), np.percentile(B,90)-np.percentile(B,10)))
print("median relative uncertainty on beta from the fit: %.0f%%" %
(100*np.median(eB)*np.log(10)))
print()
reg("log Rmax", lR)
reg("log Vout", lV)
reg("log Mbar", lM)
s = reg("log(Vout^2/Rmax)", lVR)
# two-variable fit: beta ~ Vout^a Rmax^b
X = np.column_stack([lV, lR, np.ones(len(B))])
coef, *_ = np.linalg.lstsq(X, B, rcond=None)
res = B - X@coef
print("\nbest 2-var fit: log beta = %.2f log Vout %+.2f log Rmax %+.2f (scatter %.2f dex)"
% (coef[0], coef[1], coef[2], np.std(res)))
print("compare: pure prediction beta = Vout^2/Rmax gives offset %.2f, scatter %.2f dex"
% (np.mean(B - lVR), np.std(B - lVR)))
Its output on our machine:
galaxies with beta > 0: 143
raw scatter of log beta: 0.27 dex (std), 0.64 dex (10-90)
median relative uncertainty on beta from the fit: 8%
log beta vs log Rmax slope -0.06 r = -0.08 residual scatter 0.27 dex
log beta vs log Vout slope +0.40 r = +0.41 residual scatter 0.25 dex
log beta vs log Mbar slope +0.08 r = +0.26 residual scatter 0.26 dex
log beta vs log(Vout^2/Rmax) slope +0.76 r = +0.85 residual scatter 0.14 dex
best 2-var fit: log beta = 1.67 log Vout -1.07 log Rmax +0.70 (scatter 0.11 dex)
compare: pure prediction beta = Vout^2/Rmax gives offset -0.05, scatter 0.16 dex
Appendix D: the \( z \)-extension verifier and the two clocks
Theorem 2 rests on Theorem 1 plus scalar biharmonicity, so its verifier is a page of Laplacians instead of a thirteen dimensional curvature run; Proposition 1 rests on the exact geodesic integration. Both live in one script, with the negative controls that make the zeros earned. Save it with a .py extension and run python3 appendix_d.py (SymPy, NumPy, SciPy). The extension matters: the Sage preparser only fires on .sage files, so sage appendix_d.py also runs it verbatim, while a .sage copy gets its integer literals rewritten and every zero test turns false. And read the expected output carefully: the two negative control lines are supposed to print False, that is their job; every other check must print True.
# appendix_d.py -- SAVE WITH A .py EXTENSION and run: python3 appendix_d.py
# (or: sage appendix_d.py -- Sage only preparses .sage files, so .py runs verbatim)
# the two 'negative control' lines are SUPPOSED to print False; all else must be True
import sympy as sp
# ---------- A. biharmonic verification of the z-extended hydrogen ----------
r = sp.symbols('r', positive=True)
ys = sp.symbols('y1:4'); zs = sp.symbols('z1:6')
c0,c1,c2,c3,g0,g1,g2,d0,d1,d2,e0 = sp.symbols('c0 c1 c2 c3 g0 g1 g2 d0 d1 d2 e0')
def lap(f): # radial 3D Laplacian in r + flat Laplacians in y (3) and z (5)
out = sp.diff(f, r, r) + sp.diff(f, r)*2/r
out += sum(sp.diff(f, v, 2) for v in ys) + sum(sp.diff(f, v, 2) for v in zs)
return sp.together(out)
Y = sum(v**2 for v in ys); Z = sum(v**2 for v in zs)
Q1 = ys[0]*ys[1] # two of the five l=2 harmonics of y
Q2 = ys[0]**2 - ys[1]**2
H = ( c0/r + c1 + c2*r + c3*r**2
- (g0 + g1/r)*Y - g2*(r*Y - r**3) # 3-dim trap sector (previous post)
- (d0 + d1/r)*Z - d2*(r*Z - sp.Rational(5,3)*r**3) ) # 5-dim trap sector (new)
for f_, nm in [(1/r,'1/r'), (sp.Integer(1),'1'), (r,'r'), (r**2,'r^2')]:
H += e0*f_*(zs[0]*Q1 + zs[3]*Q2) # y-z coupling riding every radial mode
print("coupling mode f=%s kept, Lap^2 H == 0:" % nm,
sp.simplify(lap(lap(H))) == 0)
H -= e0*f_*(zs[0]*Q1 + zs[3]*Q2)
# negative controls: the previous post's lambda=1 fails for 5 dims, 5/3 fails for 3
print("negative control r*Z - r^3:", sp.simplify(lap(lap(r*Z - r**3))) == 0)
print("negative control r*Y - 5/3 r^3:",
sp.simplify(lap(lap(r*Y - sp.Rational(5,3)*r**3))) == 0)
# the identity behind the trap coefficients: Lap^2(g(r)|z|^2) = |z|^2 Lap^2 g + 4n Lap g
gg = sp.Function('g')(r)
lhs = sp.simplify(lap(lap(gg*Z)) - (Z*lap(lap(gg)) + 20*lap(gg)))
print("Lap^2(g Z) == Z Lap^2 g + 20 Lap g (4n, n=5):", lhs == 0)
# ---------- B. forgetting-rate numerics: y (3 dims) vs z (5 dims) ----------
import numpy as np
from scipy.integrate import solve_ivp
ud, c2n, g1n, d1n = 1.0, 0.2, 0.4, 0.9
mu_y = np.sqrt(4*g1n/c2n - 0.25)
mu_z = np.sqrt(4*d1n/c2n - 0.25)
def rhs(tau, s):
R, Rd = s[0], s[1]; y = s[2:5]; yd = s[5:8]; z = s[8:13]; zd = s[13:18]
dHdR = c2n + (g1n*np.dot(y,y) + d1n*np.dot(z,z))/R**2
return [Rd, dHdR*ud**2, *yd, *(-2*(g1n/R)*ud**2*y), *zd, *(-2*(d1n/R)*ud**2*z)]
s0 = [5.0, 0.1, 0.3,0.2,-0.25, 0,0,0, 0.2,-0.15,0.1,0.25,-0.2, 0,0,0,0,0]
sol = solve_ivp(rhs, [1.0, 4000.0], s0, rtol=1e-11, atol=1e-13, dense_output=True)
print("\n tau R t^.5|yd| t^.5|zd| |zd|/|yd|")
for tau in np.geomspace(3, 4000, 8):
st = sol.sol(tau)
nyd = np.linalg.norm(st[5:8]); nzd = np.linalg.norm(st[13:18])
print(f"{tau:8.1f} {st[0]:10.3e} {np.sqrt(tau)*nyd:9.4f} {np.sqrt(tau)*nzd:9.4f} {nzd/nyd:9.3f}")
tt = np.linspace(3, 4000, 400000)
for idx, mu, lab in [(2, mu_y, 'y1'), (8, mu_z, 'z1')]:
q = sol.sol(tt)[idx]
zc = tt[np.where(np.diff(np.sign(q)) != 0)[0]]
print("%s zero-crossing ratios (predict e^{pi/mu}=%.4f): %s"
% (lab, np.exp(np.pi/mu), np.round(zc[-5:][1:]/zc[-5:][:-1], 4)))
Its output on our machine:
coupling mode f=1/r kept, Lap^2 H == 0: True
coupling mode f=1 kept, Lap^2 H == 0: True
coupling mode f=r kept, Lap^2 H == 0: True
coupling mode f=r^2 kept, Lap^2 H == 0: True
negative control r*Z - r^3: False
negative control r*Y - 5/3 r^3: False
Lap^2(g Z) == Z Lap^2 g + 20 Lap g (4n, n=5): True
tau R t^.5|yd| t^.5|zd| |zd|/|yd|
3.0 5.615e+00 0.2095 0.3921 1.872
8.4 1.130e+01 0.2737 0.3317 1.212
23.4 5.812e+01 0.4499 0.6492 1.443
65.5 4.296e+02 0.4784 0.4050 0.847
183.1 3.346e+03 0.4605 0.3524 0.765
511.9 2.618e+04 0.4069 0.6651 1.634
1431.0 2.047e+05 0.3218 0.1399 0.435
4000.0 1.600e+06 0.2115 0.5608 2.651
y1 zero-crossing ratios (predict e^{pi/mu}=3.0910): [3.1316 3.0877 3.088 3.0899]
z1 zero-crossing ratios (predict e^{pi/mu}=2.1079): [2.106 2.1065 2.1072 2.1075]
The envelope columns bounded with no trend while \( r \) grows by five decades, the ratio column refusing to develop a hierarchy, and each sector's zeros converging onto its own universal ratio: one rate, two clocks, measured.
Appendix E: the light verifier
Every claim of Propositions 3, 5 and 6 in one script. Save it with a .py extension and run python3 appendix_e.py; the Sage preparser only fires on .sage files, so sage appendix_e.py also works, while a .sage copy lets Sage integers seize the sympy expressions by left-multiplication and ends in the sympify traceback a reader kindly sent us. The listing below is additionally hardened, sympy objects kept on the left of every literal, but the extension is the real fix: the conservation identity for completely arbitrary \( A_\mu \), the partner wave and its refusal to decay, and the full back-reaction hierarchy of the magnetic modes through clause (iii), each mode checked with its exact sourced potential:
# appendix_e.py -- SAVE WITH A .py EXTENSION and run: python3 appendix_e.py
# (or: sage appendix_e.py -- Sage only preparses .sage files, so .py runs verbatim)
import sympy as sp
x, y, z = sp.symbols('x y z'); X = [x, y, z]
rr = sp.sqrt(x**2 + y**2 + z**2)
P2 = (z**2*3 - rr**2)/(rr**2*2)
lap = lambda f: sum(sp.diff(f, v, 2) for v in X)
def em(A):
F = [[sp.diff(A[j], X[i]) - sp.diff(A[i], X[j]) for j in range(3)] for i in range(3)]
J = [sp.simplify(sum(sp.diff(F[k][i], X[k]) for k in range(3))) for i in range(3)]
F2 = sum(F[i][j]**2 for i in range(3) for j in range(3))
rhs = sp.Rational(1,4)*lap(F2) + sp.Rational(1,2)*sum(Ji**2 for Ji in J) \
+ sum(F[i][k]*sp.diff(J[k], X[i]) for i in range(3) for k in range(3))
return F, J, sp.simplify(F2), sp.simplify(rhs)
# ---- Proposition (conservation identity): div J == 0 for ANY A ----
Agen = [sp.Function(f'A{i}')(x, y, z) for i in (1, 2, 3)]
Fg = [[sp.diff(Agen[j], X[i]) - sp.diff(Agen[i], X[j]) for j in range(3)] for i in range(3)]
Jg = [sum(sp.diff(Fg[k][i], X[k]) for k in range(3)) for i in range(3)]
print("div J == 0 identically, arbitrary A:",
sp.simplify(sum(sp.diff(Jg[i], X[i]) for i in range(3))) == 0)
# ---- Proposition (partner waves) ----
t, r = sp.symbols('t r', positive=True); G = sp.Function('G')
box = lambda f: sp.diff(f, t, t) - sp.diff(f, r, r) - sp.diff(f, r)*2/r
Phi = (t + r)*G(t - r)/(r*4)
chi = sp.simplify(box(Phi))
print("box(partner) == outgoing Maxwell wave G'/r:",
sp.simplify(chi - sp.diff(G(t-r), t)/r) == 0)
print("box^2(partner) == 0:", sp.simplify(box(chi)) == 0)
print("asymptotic amplitude (no 1/r):",
sp.limit(Phi.subs(t, sp.Symbol('v') + r), r, sp.oo))
# ---- Proposition (back-reaction hierarchy), clause (iii): Lap^2 H = RHS ----
b0, b2, b3, J0 = sp.symbols('b0 b2 b3 J0', positive=True)
# (a) dipole b0: A = b0 (-y, x, 0)/r^3 ; J = 0 and Lap H = F^2/4 (Einstein level)
A = [-b0*y/rr**3, b0*x/rr**3, sp.Integer(0)]
F, J, F2, rhs = em(A)
Hd = b0**2/(rr**4*12) + b0**2*P2/(rr**4*6)
print("\n[dipole] J == 0:", all(sp.simplify(Ji) == 0 for Ji in J),
"| Lap H = F^2/4:", sp.simplify(lap(Hd) - F2/4) == 0)
# (b) uniform field b3: A = b3 (-y, x, 0) ; J = 0, F const -> RHS = 0
A = [-b3*y, b3*x, sp.Integer(0)]
F, J, F2, rhs = em(A)
print("[uniform B] J == 0:", all(sp.simplify(Ji) == 0 for Ji in J),
"| RHS == 0:", sp.simplify(rhs) == 0)
# (c) partner mode b2: A = b2 (-y, x, 0)/r
A = [-b2*y/rr, b2*x/rr, sp.Integer(0)]
F, J, F2, rhs = em(A)
print("[b2 mode] F^2 == 4 b2^2 (1+P2)/r^2:",
sp.simplify(F2 - b2**2*(P2 + 1)*4/rr**2) == 0)
Scur = sp.simplify(rhs - sp.Rational(1,4)*lap(F2))
print("[b2 mode] current terms == -8 b2^2 P2 / r^4 (monopole-free):",
sp.simplify(Scur + b2**2*P2*8/rr**4) == 0)
print("[b2 mode] full RHS == b2^2 (2 - 12 P2)/r^4:",
sp.simplify(rhs - b2**2*(sp.Integer(2) - P2*12)/rr**4) == 0)
Hb2 = b2**2*(sp.log(rr) - P2/2)
print("[b2 mode] H = b2^2 (ln r - P2/2) solves Lap^2 H = RHS:",
sp.simplify(lap(lap(Hb2)) - rhs) == 0)
# (d) uniform current J0: F_yz = -J0 y/2, F_zx = J0 x/2 (B of a straight current)
F = [[sp.Integer(0), sp.Integer(0), -J0*x/2],
[sp.Integer(0), sp.Integer(0), -J0*y/2],
[J0*x/2, J0*y/2, sp.Integer(0)]]
print("[uniform J] Bianchi dF == 0:",
sp.simplify(sp.diff(F[1][2],X[0]) - sp.diff(F[0][2],X[1]) + sp.diff(F[0][1],X[2])) == 0)
J = [sp.simplify(sum(sp.diff(F[k][i], X[k]) for k in range(3))) for i in range(3)]
print("[uniform J] current == (0,0,-J0), harmonic:", J == [0, 0, -J0])
F2 = sum(F[i][j]**2 for i in range(3) for j in range(3))
rhs = sp.Rational(1,4)*lap(F2) + sp.Rational(1,2)*sum(Ji**2 for Ji in J) \
+ sum(F[i][k]*sp.diff(J[k], X[i]) for i in range(3) for k in range(3))
print("[uniform J] full RHS == J0^2:", sp.simplify(rhs - J0**2) == 0)
print("[uniform J] H = J0^2 r^4 / 120 solves Lap^2 H = RHS:",
sp.simplify(lap(lap(J0**2*rr**4/sp.Integer(120))) - rhs) == 0)
Its output on our machine:
div J == 0 identically, arbitrary A: True
box(partner) == outgoing Maxwell wave G'/r: True
box^2(partner) == 0: True
asymptotic amplitude (no 1/r): G(v)/2
[dipole] J == 0: True | Lap H = F^2/4: True
[uniform B] J == 0: True | RHS == 0: True
[b2 mode] F^2 == 4 b2^2 (1+P2)/r^2: True
[b2 mode] current terms == -8 b2^2 P2 / r^4 (monopole-free): True
[b2 mode] full RHS == b2^2 (2 - 12 P2)/r^4: True
[b2 mode] H = b2^2 (ln r - P2/2) solves Lap^2 H = RHS: True
[uniform J] Bianchi dF == 0: True
[uniform J] current == (0,0,-J0), harmonic: True
[uniform J] full RHS == J0^2: True
[uniform J] H = J0^2 r^4 / 120 solves Lap^2 H = RHS: True
Fourteen True lines, including the one that corrected our own working notes twice: the current terms of the \( b_2 \) mode are monopole-free, and it is the Maxwell-energy term that carries the logarithm.
Appendix F: the consolidated model table
One code path for every row of the dark matter table: the five original models and the five extensions, both regions, with exact closed-form nonnegative least squares for all linear-in-parameters models so that no row is helped or harmed by solver heuristics, MOND and the halos exactly as in Appendix C, and the single global \( \kappa \) of the universal rows optimized on the median across the zoo. The script also writes the per-galaxy table, 2,860 rows: every galaxy, region and model with its \( \chi^2/\mathrm{dof} \) and fitted parameter.
# appendix_f.py -- run: python3 appendix_f.py (NumPy + SciPy; writes model_comparison_per_galaxy.csv)
import numpy as np
from scipy.optimize import least_squares, minimize_scalar
# ---------------- data ----------------
FN = "MassModels_Lelli2016c.mrt"
gal = {}
for line in open(FN, errors='ignore'):
p = line.split()
if len(p) < 8: continue
try: vals = [float(x) for x in p[1:8]]
except ValueError: continue
gal.setdefault(p[0], []).append(vals)
KPC = 3.0857e19
GDAG = 1.2e-10 / (1e6/KPC) # (km/s)^2/kpc
UB, UD = 0.7, 0.5
Gn = 4.301e-6 # kpc (km/s)^2 / Msun
# ---------------- exact tiny NNLS helpers ----------------
def nnls2(A, y):
"""min ||Ax-y||, x>=0, A has 2 cols: interior or boundary, exact."""
best = (np.inf, np.zeros(2))
x, *_ = np.linalg.lstsq(A, y, rcond=None)
cands = []
if np.all(x >= 0): cands.append(x)
for j in (0, 1): # one variable active at 0
xj = np.zeros(2)
c = A[:, 1-j] @ y / (A[:, 1-j] @ A[:, 1-j])
xj[1-j] = max(c, 0.0)
cands.append(xj)
for xc in cands:
r = A @ xc - y; s = r @ r
if s < best[0]: best = (s, xc)
return best[1]
def fit_lin(cols, v, ev, nonneg):
"""weighted fit of v^2 on given columns; nonneg = tuple of indices with x>=0."""
w = 1.0/(2*v*ev)
A = np.column_stack(cols)*w[:,None]; y = v**2*w
x, *_ = np.linalg.lstsq(A, y, rcond=None)
if all(x[i] >= 0 for i in nonneg): return x
# active-set enumeration over the nonneg indices (<=2 of them here)
best = (np.inf, x*0)
import itertools
for act in itertools.chain.from_iterable(
itertools.combinations(nonneg, k) for k in range(len(nonneg)+1)):
keep = [j for j in range(A.shape[1]) if j not in act]
if keep:
xs, *_ = np.linalg.lstsq(A[:, keep], y, rcond=None)
else:
xs = np.array([])
xf = np.zeros(A.shape[1]); xf[keep] = xs
if any(xf[i] < -1e-12 for i in nonneg): continue
xf[list(act)] = 0.0
r = A @ xf - y; s = r @ r
if s < best[0]: best = (s, xf)
return best[1]
def chi(vm, v, ev, k):
return np.sum(((vm - v)/ev)**2)/max(len(v)-k, 1)
# ---------------- the models ----------------
def m_hyd(r, v, ev, X):
x = nnls2(np.column_stack([1/r, r])/(2*v*ev)[:,None], v**2/(2*v*ev))
vm = np.sqrt(np.clip(x[0]/r + x[1]*r, 0, None))
return chi(vm, v, ev, 2), dict(al=x[0], be=x[1])
def m_kep(r, v, ev, X):
w = 1.0/(2*v*ev)
a = max(((1/r)*w) @ (v**2*w) / (((1/r)*w) @ ((1/r)*w)), 0.0)
return chi(np.sqrt(np.clip(a/r, 0, None)), v, ev, 1), {}
def m_mond(r, v, ev, X):
vg, vd, vb = X['vg'], X['vd'], X['vb']
def c2(ud):
gb = np.clip((vg*np.abs(vg) + ud*vd**2 + UB*vb**2)/r, 1e-8, None)
go = gb/(1.0 - np.exp(-np.sqrt(gb/GDAG)))
return chi(np.sqrt(go*r), v, ev, 1)
return minimize_scalar(c2, bounds=(0.05, 5.0), method='bounded').fun, {}
def _halo(r, v, ev, X, vh2, starts, lo, hi):
vb2 = X['vg']*np.abs(X['vg']) + UD*X['vd']**2 + UB*X['vb']**2
def resid(q):
vv = np.clip(vb2 + vh2(r, 10**q[0], 10**q[1]), 0, None)
return (np.sqrt(vv) - v)/ev
best = np.inf
for q0 in starts:
try:
s = least_squares(resid, q0, bounds=(lo, hi))
best = min(best, np.sum(s.fun**2)/max(len(r)-2, 1))
except Exception: pass
return best, {}
def vnfw2(r, V200, c):
x = r/(V200/0.7)
return V200**2*(np.log(1+c*x)-c*x/(1+c*x))/(x*(np.log(1+c)-c/(1+c)))
def vburk2(r, rho0, r0):
x = r/r0
M = 2*np.pi*rho0*r0**3*(0.5*np.log(1+x**2)+np.log(1+x)-np.arctan(x))
return Gn*M/r
def m_nfw(r, v, ev, X):
return _halo(r, v, ev, X, vnfw2, ([1.9,.9],[2.3,.5],[1.5,1.2]), [1.0,0.0], [3.0,1.7])
def m_burk(r, v, ev, X):
return _halo(r, v, ev, X, vburk2, ([7.5,.5],[6.5,1.2],[8.5,0.]), [4.0,-1.0], [10.5,2.5])
def m_c1(r, v, ev, X):
x = fit_lin([1/r, np.ones_like(r), r], v, ev, nonneg=(0, 2))
vm = np.sqrt(np.clip(x[0]/r + x[1] + x[2]*r, 0, None))
return chi(vm, v, ev, 3), dict(c1=x[1])
def m_shift(r, v, ev, X, dr_fixed=None):
grid = [dr_fixed] if dr_fixed is not None else \
np.linspace(-0.9*r.min(), 3*max(X['r1'], 1.0), 280)
best = (np.inf, {})
k = 2 if dr_fixed is not None else 3
for dr in grid:
g = r + dr
if g.min() < 0.05: continue
x = nnls2(np.column_stack([1/g, g])/(2*v*ev)[:,None], v**2/(2*v*ev))
vm = np.sqrt(np.clip(x[0]/g + x[1]*g, 0, None))
c = chi(vm, v, ev, k)
if c < best[0]: best = (c, dict(dr=dr))
return best
def m_soft(r, v, ev, X, Rb_fixed=None):
grid = [Rb_fixed] if Rb_fixed is not None else np.geomspace(0.05, r.max(), 30)
best = (np.inf, {})
k = 2 if Rb_fixed is not None else 3
for Rb in grid:
g = np.sqrt(r**2 + Rb**2)
x = nnls2(np.column_stack([1/g, g])/(2*v*ev)[:,None], v**2/(2*v*ev))
vm = np.sqrt(np.clip(x[0]/g + x[1]*g, 0, None))
c = chi(vm, v, ev, k)
if c < best[0]: best = (c, dict(Rb=Rb))
return best
# ---------------- assemble per-galaxy data ----------------
G = {}
for name, rows in gal.items():
a = np.array(rows)
r, v, ev, vg, vd, vb = a[:,1], a[:,2], a[:,3], a[:,4], a[:,5], a[:,6]
ok = (r>0)&(v>0)&(ev>0)
r, v, ev, vg, vd, vb = r[ok], v[ok], ev[ok], vg[ok], vd[ok], vb[ok]
if len(r) < 8: continue
i = np.argsort(r)
r, v, ev, vg, vd, vb = r[i], v[i], ev[i], vg[i], vd[i], vb[i]
m = r >= 0.4*r.max()
if m.sum() < 5: m = np.argsort(np.argsort(r)) >= len(r)-5
# r1 from the outer 2-param fit (an outer-field property, used by both regions)
x = nnls2(np.column_stack([1/r[m], r[m]])/(2*v[m]*ev[m])[:,None], v[m]**2/(2*v[m]*ev[m]))
al, be = np.clip(x, 1e-9, None)
G[name] = dict(r=r, v=v, ev=ev, vg=vg, vd=vd, vb=vb, outer=m,
r1=np.sqrt(al/be), al=al, be=be)
def region(d, which):
m = d['outer'] if which == 'outer' else np.ones(len(d['r']), bool)
X = dict(vg=d['vg'][m], vd=d['vd'][m], vb=d['vb'][m], r1=d['r1'])
return d['r'][m], d['v'][m], d['ev'][m], X
# ---------------- universal-kappa optimization (one global parameter) ----------------
def best_kappa(model, which):
best = (np.inf, None)
for kap in np.arange(0.2, 4.01, 0.1):
chis = []
for d in G.values():
r, v, ev, X = region(d, which)
if model == 'shift':
c, _ = m_shift(r, v, ev, X, dr_fixed=kap*d['r1'])
else:
c, _ = m_soft(r, v, ev, X, Rb_fixed=kap*d['r1'])
chis.append(c)
med = np.median(chis)
if med < best[0]: best = (med, kap)
return best[1]
kaps = {}
for which in ('outer', 'full'):
for model in ('shift', 'soft'):
kaps[(model, which)] = best_kappa(model, which)
# ---------------- run everything ----------------
MODELS = [
("hydrogen(2p)", lambda r,v,e,X: m_hyd(r,v,e,X)),
("Kepler(1p)", lambda r,v,e,X: m_kep(r,v,e,X)),
("MOND/RAR(1p)", lambda r,v,e,X: m_mond(r,v,e,X)),
("NFW(2p+bar)", lambda r,v,e,X: m_nfw(r,v,e,X)),
("Burkert(2p+bar)", lambda r,v,e,X: m_burk(r,v,e,X)),
("hydro+c1(3p)", lambda r,v,e,X: m_c1(r,v,e,X)),
("hydro+dr(3p)", lambda r,v,e,X: m_shift(r,v,e,X)),
("hydro+dr=kr1(2p+g)",None), # filled per region with global kappa
("hydro+soft(3p)", lambda r,v,e,X: m_soft(r,v,e,X)),
("hydro+soft=kr1(2p+g)", None),
]
csv = ["galaxy,region,model,chi2dof,param"]
def run(which):
R = {nm: [] for nm, _ in MODELS}
P = {}
for name, d in G.items():
r, v, ev, X = region(d, which)
for nm, fn in MODELS:
if nm == "hydro+dr=kr1(2p+g)":
c, pr = m_shift(r, v, ev, X, dr_fixed=kaps[('shift', which)]*d['r1'])
elif nm == "hydro+soft=kr1(2p+g)":
c, pr = m_soft(r, v, ev, X, Rb_fixed=kaps[('soft', which)]*d['r1'])
else:
c, pr = fn(r, v, ev, X)
R[nm].append(c)
pv = next(iter(pr.values())) if pr else ""
csv.append("%s,%s,%s,%.4f,%s" % (name, which, nm, c,
("%.4g" % pv) if pv != "" else ""))
base = np.array(R["hydrogen(2p)"])
print("== %s region == galaxies: %d (global kappa: shift=%.1f, soft=%.1f)"
% (which, len(base), kaps[('shift', which)], kaps[('soft', which)]))
print("%-21s%16s%7s%7s%17s" % ("model", "median chi2/dof", "<1", "<2", "beats hydrogen"))
for nm, _ in MODELS:
A = np.array(R[nm])
b = "" if nm == "hydrogen(2p)" else "%d/%d (%.0f%%)" % (
(A < base).sum(), len(A), 100*(A < base).mean())
print("%-21s%16.2f%6.0f%%%6.0f%%%17s"
% (nm, np.median(A), 100*(A<1).mean(), 100*(A<2).mean(), b))
print()
run('outer'); run('full')
open('model_comparison_per_galaxy.csv','w').write("\n".join(csv)+"\n")
print("per-galaxy table written: model_comparison_per_galaxy.csv (%d rows)" % (len(csv)-1))
Its output on our machine, verbatim:
== outer region == galaxies: 143 (global kappa: shift=1.8, soft=1.3)
model median chi2/dof <1 <2 beats hydrogen
hydrogen(2p) 0.67 59% 77%
Kepler(1p) 23.85 0% 1% 0/143 (0%)
MOND/RAR(1p) 1.94 35% 52% 33/143 (23%)
NFW(2p+bar) 0.34 73% 85% 104/143 (73%)
Burkert(2p+bar) 0.28 76% 84% 99/143 (69%)
hydro+c1(3p) 0.37 76% 87% 120/143 (84%)
hydro+dr(3p) 0.32 76% 87% 122/143 (85%)
hydro+dr=kr1(2p+g) 0.40 69% 82% 98/143 (69%)
hydro+soft(3p) 0.40 73% 82% 90/143 (63%)
hydro+soft=kr1(2p+g) 0.47 67% 80% 105/143 (73%)
== full region == galaxies: 143 (global kappa: shift=0.5, soft=0.5)
model median chi2/dof <1 <2 beats hydrogen
hydrogen(2p) 11.98 7% 15%
Kepler(1p) 313.29 0% 0% 0/143 (0%)
MOND/RAR(1p) 3.51 15% 31% 114/143 (80%)
NFW(2p+bar) 1.85 36% 54% 133/143 (93%)
Burkert(2p+bar) 1.02 50% 64% 136/143 (95%)
hydro+c1(3p) 4.09 18% 30% 125/143 (87%)
hydro+dr(3p) 3.93 13% 29% 125/143 (87%)
hydro+dr=kr1(2p+g) 7.27 7% 18% 88/143 (62%)
hydro+soft(3p) 4.30 11% 22% 84/143 (59%)
hydro+soft=kr1(2p+g) 7.27 8% 17% 86/143 (60%)
per-galaxy table written: model_comparison_per_galaxy.csv (2860 rows)
Appendix G: the bar prescription and its verifiers
Three pieces of code carry the bar section. The prescription itself, exactly as the hypothesis demands: per galaxy, all of \( \alpha, \beta, \delta r \) free with a data-driven grid, the bar region excluded, and the exclusion radius forced to be as small as the data allow. The transparency column, the model evaluated on every point including the excluded ones, is printed on principle.
# appendix_g_barcut.py -- run: python3 appendix_g_barcut.py (writes bar_cut_fits.csv)
import numpy as np
FN = "MassModels_Lelli2016c.mrt"
gal = {}
for line in open(FN, errors='ignore'):
p = line.split()
if len(p) < 8: continue
try: vals = [float(x) for x in p[1:8]]
except ValueError: continue
gal.setdefault(p[0], []).append(vals)
def nnls2(A, y):
best = (np.inf, np.zeros(2))
x, *_ = np.linalg.lstsq(A, y, rcond=None)
cands = [x] if np.all(x >= 0) else []
for j in (0, 1):
xj = np.zeros(2)
c = A[:,1-j] @ y / (A[:,1-j] @ A[:,1-j])
xj[1-j] = max(c, 0.0)
cands.append(xj)
for xc in cands:
s = np.sum((A @ xc - y)**2)
if s < best[0]: best = (s, xc)
return best[1]
def fit_shift(r, v, ev):
"""free (alpha, beta, delta_r); delta_r grid is data-driven only."""
w = 1.0/(2*v*ev); y = v**2*w
best = (np.inf, 0.0, 0.0, 0.0)
for dr in np.linspace(-0.9*r.min(), 2.0*r.max(), 240):
g = r + dr
if g.min() < 0.05: continue
x = nnls2(np.column_stack([1/g, g])*w[:,None], y)
vm = np.sqrt(np.clip(x[0]/g + x[1]*g, 0, None))
ch = np.sum(((vm - v)/ev)**2)/max(len(r)-3, 1)
if ch < best[0]: best = (ch, dr, x[0], x[1])
return best
def bar_cut_fit(r, v, ev, thresh):
"""smallest cut R_bar such that the shifted fit on r > R_bar has chi2/dof <= thresh.
Returns (R_bar, n_cut, chi_retained, dr, al, be); R_bar=0 means no cut needed."""
n = len(r)
for k in range(0, n - 4): # keep at least 5 points
ch, dr, al, be = fit_shift(r[k:], v[k:], ev[k:])
if ch <= thresh:
Rb = 0.0 if k == 0 else 0.5*(r[k-1] + r[k]) # boundary between cut & kept
return Rb, k, ch, dr, al, be
ch, dr, al, be = fit_shift(r[n-5:], v[n-5:], ev[n-5:])
return 0.5*(r[n-6]+r[n-5]), n-5, ch, dr, al, be
rows = []
for name, dat in gal.items():
a = np.array(dat)
D, r, v, ev = a[0,0], a[:,1], a[:,2], a[:,3]
ok = (r>0)&(v>0)&(ev>0)
r, v, ev = r[ok], v[ok], ev[ok]
if len(r) < 8: continue
i = np.argsort(r); r, v, ev = r[i], v[i], ev[i]
Rb, k, ch, dr, al, be = bar_cut_fit(r, v, ev, thresh=1.5)
# transparency: the final model evaluated on ALL points, cut ones included
g = np.maximum(r + dr, 0.05)
vm = np.sqrt(np.clip(al/g + be*g, 0, None))
ch_full = np.sum(((vm - v)/ev)**2)/max(len(r)-3, 1)
r1 = np.sqrt(al/be) if (al > 0 and be > 0) else np.nan
rows.append(dict(name=name, Rb=Rb, kcut=k, n=len(r), frac=k/len(r),
chi=ch, chif=ch_full, dr=dr, al=al, be=be, r1=r1,
Rmax=r.max()))
Rb = np.array([x['Rb'] for x in rows])
frac = np.array([x['frac'] for x in rows])
chi = np.array([x['chi'] for x in rows])
chif = np.array([x['chif'] for x in rows])
dr = np.array([x['dr'] for x in rows])
r1 = np.array([x['r1'] for x in rows])
Rmax = np.array([x['Rmax'] for x in rows])
print("== bar-cut + free delta_r prescription (thresh chi2/dof <= 1.5), n = %d ==" % len(rows))
print("no cut needed (R_bar = 0): %d galaxies (%.0f%%)"
% ((Rb == 0).sum(), 100*(Rb == 0).mean()))
print("points excluded: median %.0f%% 16-84%%: %.0f%% - %.0f%%"
% (100*np.median(frac), *[100*q for q in np.percentile(frac, [16,84])]))
has = Rb > 0
print("R_bar (cut>0): median %.2f kpc 16-84%%: %.2f - %.2f median R_bar/Rmax = %.2f"
% (np.median(Rb[has]), *np.percentile(Rb[has], [16,84]),
np.median(Rb[has]/Rmax[has])))
print("chi2/dof on retained: median %.2f on ALL points (cut included): median %.2f"
% (np.median(chi), np.median(chif)))
print("delta_r: %.0f%% positive median %+.2f kpc" % (100*np.mean(dr>0), np.median(dr)))
# a-posteriori correlations only (nothing here entered the fits)
ok2 = has & np.isfinite(r1) & (dr > 0.05)
for lab, x in [("r1", r1), ("Rmax", Rmax), ("delta_r", dr)]:
m = has & np.isfinite(x) & (x > 0.05)
lx, ly = np.log10(x[m]), np.log10(Rb[m])
p = np.polyfit(lx, ly, 1)
print("a-posteriori: log R_bar vs log %-7s slope %+5.2f r = %+5.2f (n=%d)"
% (lab, p[0], np.corrcoef(lx, ly)[0,1], m.sum()))
print("median R_bar/r1 = %.2f median R_bar/delta_r = %.2f"
% (np.median((Rb/r1)[ok2]), np.median((Rb/dr)[ok2])))
# robustness: stricter threshold
rows1 = []
for name, dat in gal.items():
a = np.array(dat)
r, v, ev = a[:,1], a[:,2], a[:,3]
ok = (r>0)&(v>0)&(ev>0); r, v, ev = r[ok], v[ok], ev[ok]
if len(r) < 8: continue
i = np.argsort(r); r, v, ev = r[i], v[i], ev[i]
Rb1, k1, *_ = bar_cut_fit(r, v, ev, thresh=1.0)
rows1.append((Rb1, k1/len(r)))
Rb1 = np.array([x[0] for x in rows1]); f1 = np.array([x[1] for x in rows1])
print("\nrobustness (thresh 1.0): no-cut %.0f%%, median excluded %.0f%%, median R_bar %.2f kpc"
% (100*np.mean(Rb1==0), 100*np.median(f1), np.median(Rb1[Rb1>0])))
# the one external point so far, under the proper prescription
for x in rows:
if x['name'] == 'ESO116-G012':
print("\nESO116-G012: fitted R_bar = %.2f kpc, delta_r = %+.2f kpc, "
"chi2 retained %.2f, excluded %d/%d pts | S4G bar = 1.51 kpc (quality 3, sky)"
% (x['Rb'], x['dr'], x['chi'], x['kcut'], x['n']))
hdr = "galaxy,n_points,n_cut,R_bar_kpc,delta_r_kpc,alpha,beta,r1_kpc,chi2dof_retained,chi2dof_all"
out = [hdr] + ["%s,%d,%d,%.3f,%.3f,%.5g,%.5g,%.3f,%.3f,%.3f"
% (x['name'], x['n'], x['kcut'], x['Rb'], x['dr'], x['al'], x['be'],
x['r1'], x['chi'], x['chif']) for x in rows]
open('bar_cut_fits.csv','w').write("\n".join(out)+"\n")
print("\nwrote bar_cut_fits.csv (%d galaxies)" % len(rows))
Its output on our machine:
== bar-cut + free delta_r prescription (thresh chi2/dof <= 1.5), n = 143 ==
no cut needed (R_bar = 0): 27 galaxies (19%)
points excluded: median 14% 16-84%: 0% - 50%
R_bar (cut>0): median 2.66 kpc 16-84%: 1.08 - 9.84 median R_bar/Rmax = 0.19
chi2/dof on retained: median 1.11 on ALL points (cut included): median 6.82
delta_r: 82% positive median +10.89 kpc
a-posteriori: log R_bar vs log r1 slope +0.33 r = +0.54 (n=62)
a-posteriori: log R_bar vs log Rmax slope +1.07 r = +0.78 (n=116)
a-posteriori: log R_bar vs log delta_r slope +0.45 r = +0.57 (n=103)
median R_bar/r1 = 0.10 median R_bar/delta_r = 0.14
robustness (thresh 1.0): no-cut 13%, median excluded 22%, median R_bar 3.13 kpc
ESO116-G012: fitted R_bar = 3.91 kpc, delta_r = +19.72 kpc, chi2 retained 0.22, excluded 7/15 pts | S4G bar = 1.51 kpc (quality 3, sky)
wrote bar_cut_fits.csv (143 galaxies)
The degeneracy check behind Proposition 8's warning label, shift against free constant on identical outer windows with identical solvers:
import numpy as np
FN = "MassModels_Lelli2016c.mrt"
gal = {}
for line in open(FN, errors='ignore'):
p = line.split()
if len(p) < 8: continue
try: vals = [float(x) for x in p[1:8]]
except ValueError: continue
gal.setdefault(p[0], []).append(vals)
def nnls2(A, y):
best = (np.inf, np.zeros(2))
x, *_ = np.linalg.lstsq(A, y, rcond=None)
cands = [x] if np.all(x >= 0) else []
for j in (0, 1):
xj = np.zeros(2)
c = A[:,1-j] @ y / (A[:,1-j] @ A[:,1-j])
xj[1-j] = max(c, 0.0); cands.append(xj)
for xc in cands:
s = np.sum((A @ xc - y)**2)
if s < best[0]: best = (s, xc)
return best[1]
def fit_c1(r, v, ev):
"""v^2 = a/r + c1 + b r with a,b >= 0, c1 free: active-set enumeration."""
import itertools
w = 1.0/(2*v*ev)
A = np.column_stack([1/r, np.ones_like(r), r])*w[:,None]; y = v**2*w
best = (np.inf, np.zeros(3))
for act in itertools.chain.from_iterable(
itertools.combinations((0, 2), k) for k in range(3)):
keep = [j for j in range(3) if j not in act]
xs, *_ = np.linalg.lstsq(A[:, keep], y, rcond=None)
xf = np.zeros(3); xf[keep] = xs
if xf[0] < -1e-12 or xf[2] < -1e-12: continue
s = np.sum((A @ xf - y)**2)
if s < best[0]: best = (s, xf)
x = best[1]
vm = np.sqrt(np.clip(x[0]/r + x[1] + x[2]*r, 0, None))
return np.sum(((vm-v)/ev)**2)/max(len(r)-3, 1), x[1]
cs, cc, drs, bds = [], [], [], []
for name, rows in gal.items():
a = np.array(rows)
r, v, ev = a[:,1], a[:,2], a[:,3]
ok = (r>0)&(v>0)&(ev>0); r, v, ev = r[ok], v[ok], ev[ok]
if len(r) < 8: continue
i = np.argsort(r); r, v, ev = r[i], v[i], ev[i]
m = r >= 0.4*r.max()
if m.sum() < 5: m = np.argsort(np.argsort(r)) >= len(r)-5
ro, vo, eo = r[m], v[m], ev[m]
w = 1.0/(2*vo*eo); y = vo**2*w
best = (np.inf, 0, 0, 0)
for dr in np.linspace(-0.9*ro.min(), 2.0*ro.max(), 240):
g = ro + dr
if g.min() < 0.05: continue
x = nnls2(np.column_stack([1/g, g])*w[:,None], y)
vm = np.sqrt(np.clip(x[0]/g + x[1]*g, 0, None))
ch = np.sum(((vm-vo)/eo)**2)/max(len(ro)-3, 1)
if ch < best[0]: best = (ch, dr, x[0], x[1])
chS, dr, alS, beS = best
chC, c1 = fit_c1(ro, vo, eo)
cs.append(chS); cc.append(chC); drs.append(dr); bds.append(beS*dr)
cc[-1] = chC
cs, cc, drs, bds = map(np.array, (cs, cc, drs, bds))
print("outer window, n = %d galaxies" % len(cs))
print("3-param shift (alpha,beta,dr): median chi2/dof %.2f" % np.median(cs))
print("3-param const (alpha,c1,beta): median chi2/dof %.2f" % np.median(cc))
print("indistinguishable (|diff|<=0.05): %.0f%% shift wins: %.0f%% const wins: %.0f%%"
% (100*np.mean(np.abs(cs-cc)<=0.05), 100*np.mean(cs<cc-0.05), 100*np.mean(cc<cs-0.05)))
pos = drs > 0.05
print("Proposition-8 identity, measured: median c1 of the constant fits vs median beta*dr"
" of the shift fits:")
# recompute c1 medians from the constant fits stored above
c1s = []
for name, rows in gal.items():
a = np.array(rows)
r, v, ev = a[:,1], a[:,2], a[:,3]
ok = (r>0)&(v>0)&(ev>0); r, v, ev = r[ok], v[ok], ev[ok]
if len(r) < 8: continue
i = np.argsort(r); r, v, ev = r[i], v[i], ev[i]
m = r >= 0.4*r.max()
if m.sum() < 5: m = np.argsort(np.argsort(r)) >= len(r)-5
c1s.append(fit_c1(r[m], v[m], ev[m])[1])
print(" median c1 = %+.0f (km/s)^2 median beta*dr = %+.0f (km/s)^2"
% (np.median(c1s), np.median(bds[pos])))
outer window, n = 143 galaxies
3-param shift (alpha,beta,dr): median chi2/dof 0.32
3-param const (alpha,c1,beta): median chi2/dof 0.37
indistinguishable (|diff|<=0.05): 78% shift wins: 15% const wins: 6%
Proposition-8 identity, measured: median c1 of the constant fits vs median beta*dr of the shift fits:
median c1 = +5486 (km/s)^2 median beta*dr = +3093 (km/s)^2
And the two-line existence proof of the dynamical alternative, the exact \( m=2 \) bar-shaped growth modes of the theory:
import sympy as sp
x, y, z = sp.symbols('x y z')
lap = lambda f: sum(sp.diff(f, u_, 2) for u_ in (x, y, z))
S22 = x**2 - y**2
print(sp.simplify(lap(lap(S22))) == 0,
sp.simplify(lap(lap((x**2 + y**2 + z**2)*S22))) == 0)
True True
References
2. K. S. Stelle, Renormalization of higher-derivative quantum gravity, Phys. Rev. D 16, 953 (1977)
3. K. S. Stelle, Classical gravity with higher derivatives, Gen. Rel. Grav. 9, 353 (1978)
4. T. Málek, V. Pravda, Type III and N solutions to quadratic gravity, Phys. Rev. D 84, 024047 (2011)
5. V. Pravda, A. Pravdová, J. Podolský, R. Švarc, Exact solutions to quadratic gravity, Phys. Rev. D 95, 084025 (2017)
6. V. Pravda, A. Pravdová, A. Coley, R. Milson, All spacetimes with vanishing curvature invariants, Class. Quant. Grav. 19, 6213 (2002)
7. H. Lü, C. N. Pope, Critical Gravity in Four Dimensions, Phys. Rev. Lett. 106, 181302 (2011)
9. V. P. Frolov, D. V. Fursaev, Gravitational field of a spinning radiation beam-pulse
10. P. Candelas, G. Horowitz, A. Strominger, E. Witten, Vacuum configurations for superstrings, Nucl. Phys. B 258, 46 (1985)
11. Holonomy group and the Berger classification
13. SPARC database: Spitzer Photometry and Accurate Rotation Curves
14. F. Lelli, S. S. McGaugh, J. M. Schombert, SPARC: Mass Models for 175 Disk Galaxies, Astron. J. 152, 157 (2016)
15. S. S. McGaugh, F. Lelli, J. M. Schombert, The Radial Acceleration Relation in Rotationally Supported Galaxies, Phys. Rev. Lett. 117, 201101 (2016)
16. J. F. Navarro, C. S. Frenk, S. D. M. White, The Structure of Cold Dark Matter Halos, Astrophys. J. 462, 563 (1996)
17. A. Burkert, The Structure of Dark Matter Halos in Dwarf Galaxies, Astrophys. J. 447, L25 (1995)
19. Etherington's distance-duality relation
20. M. Herrera-Endoqui, S. Díaz-García, E. Laurikainen, H. Salo, Catalogue of the morphological features in the S4G, A&A 582, A86 (2015); VizieR J/A+A/582/A86
21. S. Díaz-García, H. Salo, E. Laurikainen, M. Herrera-Endoqui, Characterization of galactic bars from 3.6 μm S4G imaging, A&A 587, A160 (2016); VizieR J/A+A/587/A160
22. P. Erwin, How large are the bars in barred galaxies?, MNRAS 364, 283 (2005)
23. M. Tremaine, S. Weinberg, A kinematic method for measuring the pattern speed of barred galaxies, Astrophys. J. 282, L5 (1984)
24. F. Zwicky, On the Red Shift of Spectral Lines through Interstellar Space, PNAS 15, 773 (1929)
25. G. Goldhaber et al., Timescale Stretch Parameterization of Type Ia Supernova B-band Light Curves, Astrophys. J. 558, 359 (2001)
26. S. Blondin et al., Time Dilation in Type Ia Supernova Spectra at High Redshift, Astrophys. J. 682, 724 (2008)
27. L. M. Lubin, A. Sandage, The Tolman Surface Brightness Test for the Reality of the Expansion. IV, Astron. J. 122, 1084 (2001)
28. A. G. Riess et al., Observational Evidence from Supernovae for an Accelerating Universe, Astron. J. 116, 1009 (1998)
29. S. Perlmutter et al., Measurements of Omega and Lambda from 42 High-Redshift Supernovae, Astrophys. J. 517, 565 (1999)
30. F. Hohl, Numerical Experiments with a Disk of Stars, Astrophys. J. 168, 343 (1971)
31. J. P. Ostriker, P. J. E. Peebles, A Numerical Study of the Stability of Flattened Galaxies, Astrophys. J. 186, 467 (1973)
32. V. P. Debattista, J. A. Sellwood, Constraints from Dynamical Friction on the Dark Matter Content of Barred Galaxies, Astrophys. J. 543, 704 (2000)
33. E. Athanassoula, What determines the strength and the slowdown rate of bars?, MNRAS 341, 1179 (2003)
34. J. A. Sellwood, Secular Evolution in Disk Galaxies, Rev. Mod. Phys. 86, 1 (2014)
35. K. Sheth et al., Evolution of the Bar Fraction in COSMOS: Quantifying the Assembly of the Hubble Sequence, Astrophys. J. 675, 1141 (2008)
Cite
If you found this work useful, please consider citing:
@misc{hadilq2026MonismII,
author = {{Hadi Lashkari Ghouchani}},
note = {Published electronically at \url{https://hadilq.com/posts/geodesic-monism-ii/}},
gitlab = {Gitlab source at \href{https://gitlab.com/hadilq/hadilq.gitlab.io/-/tree/main/content/posts/2026-07-28-geodesic-monism-ii}},
title = {Geodesic Monism and its hydrogen solution},
year={2026},
}