FCMFD: Algebraic Substrate Theory of Gravitational Relaxation

FCMFD: A Unified Algebraic Substrate Theory of Gravitational Relaxation

From Galactic Halos to Binary Mergers — A Constitutive Framework without Dark Matter

Version 1.0 — Public Release
Based on numerical experiments in Google Colab using SPARC rotation curves, matrix network simulations, and LIGO ringdown data.

Abstract

We present the complete mathematical formulation of Finite Coupled Monad Field Dynamics (FCMFD) — a coordinate‑free, non‑abelian matrix network theory in which gravitational phenomena emerge from a finite‑capacity substrate with fractional memory and amplitude saturation. No background spacetime geometry is assumed. The state of the universe is represented by a set of Hermitian operators M_i(t) coupled via an adjacency matrix W_ij. Stress propagates through commutator torque [M_i, M_j]; memory is implemented via a Grünwald‑Letnikov fractional time derivative; and saturation (“phase locking”) is enforced by spectral truncation.

The framework is tested across three distinct physical regimes:

Galactic halos (SPARC database): outer rotation curve residuals are well described by a spatial power‑law ΔV² ∝ R^2γ with median γ = 0.272, statistically indistinguishable from a two‑parameter NFW dark matter halo in 96% of the usable sample.
Single‑core magnetar‑like relaxation: matrix network simulations yield an effective stretched‑exponent β_eff ≈ 0.29, matching the fractional order α = 0.35 used in the memory kernel.
Binary inspiral chirp: time‑varying coupling produces a rising frequency ridge in the lag‑rate spectrogram, with global strain exponent β_eff ≈ 0.21.
Black hole merger ringdown (GW150914): the same pipeline extracts a pure exponential tail (β_eff = 1.00), consistent with the saturation limit of the substrate.

These results are unified by a stress‑dependent transition function β_eff(T) that evolves from the fractional regime (β ≈ 0.3) at low stress to the Markovian limit (β = 1.0) at high stress. The paper provides the complete mathematical core of FCMFD, including all derived parameters and open assumptions.

1. Introduction

General Relativity describes gravity as curvature of a differentiable spacetime manifold. While highly successful, it contains singularities and requires a pre‑existing geometric stage. Finite Coupled Monad Field Dynamics (FCMFD) replaces this geometric picture with a discrete, algebraic substrate of finite capacity. The “Monad Field” is a set of interacting matrix operators whose internal state, memory, and saturation determine all physical interactions. Geometry is not fundamental; it emerges from the stress and coupling patterns of the network.

This white paper documents the complete mathematical framework as implemented in our Colab numerical experiments. All parameters and equations are derived from reproducible simulations using open data (SPARC, LIGO GW150914) and custom matrix network solvers. The goal is to present a self‑contained, falsifiable alternative to dark matter and ΛCDM phenomenology, with explicit caveats and remaining open questions.

2. Core Mathematical Framework

2.1 State Representation

The Monad Field consists of N nodes (N = 4 in our simulations). Each node i at time t is represented by a Hermitian 2×2 matrix M_i(t):

M_i(t) = a_i(t)·I + b_i(t)·σ_x + c_i(t)·σ_y + d_i(t)·σ_z

where I is the identity, σ_x, σ_y, σ_z are Pauli matrices, and a_i, b_i, c_i, d_i ∈ ℝ. The eigenvalues λ₁(t), λ₂(t) of M_i(t) are interpreted as the local “field strength” or stress. In the ground state, all eigenvalues are zero. The maximum allowed eigenvalue is bounded by a saturation ceiling S_max (set to 2.0 in our runs).

Derivation from simulations: The parameters a_i, b_i, c_i, d_i are evolved via the dynamical equations below. Initial conditions for the binary inspiral used:

Node 0: (a,b,c,d) = (1.2, 0.8, 0.6, 0.4)
Node 1: (a,b,c,d) = (1.2, -0.8, -0.6, 0.4)
Nodes 2,3: small off‑diagonal terms to break commutator degeneracy.

2.2 Network Topology (Adjacency Matrix)

Interactions are defined by a symmetric adjacency matrix W_ij. For the binary inspiral, we used a ring topology with background coupling 0.2 between neighbouring nodes, plus a time‑varying binary coupling between nodes 0 and 1:

W₀₁(t) = W₁₀(t) = C(t) · cos²(Φ(t))

where C(t) increases linearly from ~1.5 to ~8.0 over the merger timescale t_merge = 3.5 s, and the phase Φ(t) follows an accelerating chirp:

Φ(t) = 2π [ f_start·t + (f_ramp/3)·t³ ] with f_start = 2.0 Hz, f_ramp = 8.0 Hz/s².

Assumption: The time‑varying adjacency is a phenomenological model of inspiral. A full derivation from first principles (e.g., from a Lagrangian) is still missing.

2.3 Fractional Time Evolution (Memory)

We approximate a Caputo fractional derivative of order α (0 < α ≤ 1) using the Grünwald‑Letnikov definition with alternating signs. The weights w_j are computed recursively:

w₀ = 1, w_j = – w_j‑1·(α – j + 1)/j (j ≥ 1)

The state update for node i at time step n (time t_n = n·Δt) is:

M_i⁽ⁿ⁾ = – Σ_j=1ⁿ w_j · M_i^(n‑j) + (Δt)^α · F_i^(n‑1)

where F_i is the “force” term (commutator torque + damping + driving). For α = 1 the scheme reduces to ordinary Euler integration.

Derived value: All simulations used α = 0.35, which was chosen to match the stretched‑exponent behaviour observed in earlier scalar fractional relaxation studies (Step 4). This value is not derived from first principles but is a fixed parameter that yields β_eff ≈ 0.3.

Placeholder: The rigorous discretisation of the Caputo derivative would require normalisation factors involving Γ(2−α). Our heuristic form is a “GL‑inspired” memory kernel, not a validated fractional PDE solver.

2.4 Dynamical Force Terms

The force term F_i contains three contributions:

Commutator torque (stress propagation):
T_i = Σ_j W_ij · [M_i, M_j] where [M_i, M_j] = M_i M_j – M_j M_i.
This term is coordinate‑free and encodes torque transmission without a spatial Laplacian.
Nonlinear damping (cubic, ensures boundedness):
D_i = –γ · M_i³
with γ = 0.4–0.5 in our runs. This term mimics energy loss and prevents unbounded growth.
External driving J_i(t): for binary inspiral, only nodes 0 and 1 receive a driving term proportional to the chirp phase derivative. In the single‑core relaxation, J = 0.

Thus F_i = T_i + D_i + J_i.

2.5 Saturation (Phase Locking)

After each update, we compute the eigenvalues λ_k of the proposed M_i. If any eigenvalue exceeds S_max in absolute value, we scale the whole matrix so that the maximum eigenvalue becomes exactly S_max:

If max(|λ_k|) > S_max ⇒ λ_k → λ_k · (S_max / max(|λ_k|))

After scaling, the matrix is reassembled: M_i^new = V diag(λ_scaled) V†.

Saturation ceiling: S_max = 2.0 in all simulations. This value is arbitrary; it sets the scale for what “saturated” means. In a future physical theory, S_max would be related to the Planck scale or another fundamental invariant.

Note on clipping: In early versions we used np.clip, which is a hard cut‑off. The spectral scaling method is more physical because it preserves eigenvector directions while limiting the operator norm.

2.6 Observational Mapping (Magnetar Lag Analogue)

To connect the algebraic state to a timing residual (e.g., pulsar timing or gravitational wave strain), we define an effective scalar observable S_eff(t). Three choices were used:

Single‑node: S_eff(t) = λ_max(t) of node 0 (largest eigenvalue).
Global strain: S_eff(t) = (1/N) Σ_i λ_max⁽ⁱ⁾(t).
Commutator torque: S_eff(t) = || [M₀, M₁] ||_F = √( Tr( [M₀,M₁]†[M₀,M₁] ) ).

This scalar is then passed through a phenomenological magnetar‑lag pipeline:

dθ_eff/dt = Ω / (1 + k·S_eff(t)) (Ω = 2π·10 rad/s, k = 5.0)

lag(t) = (θ_eff(t) – Ω t) / (2π)

rate(t) = | d(lag)/dt |

The tail of the rate (t > 0.5 s after the peak) is fitted to a stretched exponential:

rate(t) = A · exp( – (t/τ)^β_eff )

The fitted β_eff is the diagnostic exponent.

Derived values from simulations:

Single saturated core (v6.3): β_eff ≈ 0.292
Binary inspiral – global strain (v6.5): β_eff ≈ 0.206
Binary inspiral – stochastic core (v6.6): β_eff ≈ 0.188
Binary inspiral – commutator norm (v6.8): β_eff ≈ 0.305

2.7 Spatial Power‑Law for SPARC Galaxies (γ)

For galactic rotation curves, we treat radius R as the independent variable. Define baryonic velocity squared:

V_bar² = V_gas² + V_disk² + V_bulge²

Observed squared velocity: V_obs²

Residual excess: ΔV²(R) = V_obs² – V_bar²

We fit the outer region (R > 0.5 R_max, ΔV² > 0) with a power law:

ln(ΔV²) = ln(V₀²) + 2γ · ln(R/R₀)

where γ is the spatial scaling exponent. The fit uses weighted least squares with propagated errors:

σ_ΔV² = √[ (2·V_obs·σ_Vobs)² + (2·σ_baryon·V_bar²)² ]

with σ_baryon = 0.10 (10% systematic uncertainty in baryonic mass‑to‑light ratios).

Derived values from SPARC (92 galaxies):

Median γ = 0.272
Mean γ = 0.361 ± 0.361
Weak correlation with V_max: r = 0.208, p = 0.046

2.8 Model Comparison (BIC)

For each galaxy, we also fit a standard NFW dark matter halo (two parameters: V_c, c) to the same ΔV²(R) data, using the same weighting. The Bayesian Information Criterion (BIC) is computed for both models:

BIC = χ² + k·ln(N_points) (k = 2 for both models)

ΔBIC = BIC_FCMFD – BIC_NFW.

Interpretation:

ΔBIC < –2: FCMFD power‑law preferred
|ΔBIC| ≤ 2: indistinguishable
ΔBIC > 2: NFW preferred

Results: 4 galaxies FCMFD‑preferred, 88 indistinguishable, 0 NFW‑preferred.

3. Cross‑Domain Consistency and the Stress‑Dependent Transition

The recovered exponents form a sequence:

Low stress (galactic halos, single core): β ≈ 0.27–0.29
Moderate stress (binary inspiral): β ≈ 0.21
High stress (LIGO ringdown): β = 1.00

We interpret this as a transition from fractional‑memory behaviour to Markovian exponential decay as the local field stress approaches the saturation ceiling S_max. A phenomenological transition function is proposed:

β_eff(T) = β₀ + (1 – β₀) · [ (T/T_crit)ⁿ / (1 + (T/T_crit)ⁿ) ]

where:

β₀ ≈ 0.29 is the zero‑stress fractional exponent (derived from simulations)
T is a proxy for the local stress (e.g., the trace of the energy‑momentum tensor, or the commutator norm)
T_crit is the critical stress at which saturation dominates
n controls the sharpness of the transition

Placeholder: The exact form of T in terms of the Monad Field operators is not yet defined. In our current work, we simply match the observed β_eff values to the empirical regime. A first‑principles derivation of the transition from the fundamental equations (commutator dynamics + saturation) remains an open problem.

4. Discussion and Limitations

4.1 What the framework currently achieves

A coordinate‑free, background‑independent algebraic formulation of a finite‑capacity substrate.
Reproducible numerical simulations that produce stretched‑exponential relaxation (β ≈ 0.29) and binary chirp spectrograms.
A successful fit of outer rotation curve residuals with a power‑law scaling (γ ≈ 0.27) that is statistically indistinguishable from NFW in 96% of the SPARC sample.
A consistent transition from β ≈ 0.3 to β = 1.0 across three distinct physical regimes.

4.2 Key assumptions and placeholders

Fractional discretisation: The GL scheme used is heuristic; rigorous convergence and stability have not been proven.
Saturation ceiling S_max: The value 2.0 is arbitrary. A physically derived saturation scale (e.g., from Planck units) is missing.
Stress‑dependent transition: The function β_eff(T) is empirical; a derivation from the fundamental equations is absent.
Baryonic systematic uncertainty σ_baryon = 0.10: This is a guess; actual mass‑to‑light uncertainties may be larger.
Outer‑radius threshold (R > 0.5 R_max): This choice influences the γ distribution; sensitivity tests are needed.
No comparison with other halo models (Burkert, pseudo‑isothermal): Only NFW was tested.
Small network size (N=4): Scaling to larger N may alter the exponents.

4.3 Falsifiability

FCMFD makes the following testable predictions:

Intermediate‑mass mergers (e.g., GW170817) should yield β_eff in the range 0.3 < β_eff < 1.0, depending on the system’s mass and compactness.
Low‑surface‑brightness (LSB) galaxies should exhibit γ ≈ 0.27 with no radial gradient, whereas massive spirals may show a decreasing γ from centre to edge.
High‑precision gravitational wave spectrograms should reveal deviations from pure exponential decay at high frequencies, characteristic of the fractional memory tail.

If future observations contradict any of these, the framework would be falsified.

5. Conclusions

We have presented the complete mathematical core of Finite Coupled Monad Field Dynamics (FCMFD), a coordinate‑free algebraic substrate theory of gravitational relaxation. The framework replaces spacetime geometry with a network of Hermitian matrices, fractional memory, and amplitude saturation. Numerical simulations reproduce stretched‑exponential exponents (β ≈ 0.29), binary chirp spectrograms, and a stress‑dependent transition to exponential decay (β = 1.0) at high energies.

Analyses of SPARC rotation curves show that a simple power‑law scaling (γ ≈ 0.27) describes outer‑halo residuals as well as a two‑parameter NFW dark matter profile in 96% of the usable sample, with a weak positive correlation with galaxy mass. The numerical coincidence between the temporal β (≈0.29) and the spatial median γ (≈0.27) suggests a constitutive constant of the substrate, though the large scatter precludes a claim of universality.

All equations and parameters are explicitly documented, and all results are reproducible from the Colab scripts provided in the Appendix (not included here). The framework is falsifiable and open to further development, particularly in deriving the stress‑dependent transition from first principles and in extending to larger network sizes and continuum limits.

Acknowledgments

This work synthesises collaborative numerical experiments conducted in Google Colab, using public data from the SPARC database (Lelli et al. 2016) and LIGO GW150914 (GWOSC). The authors thank the open‑source Python ecosystem (NumPy, SciPy, Matplotlib) for enabling the simulations.

References

Lelli, F., McGaugh, S. S., & Schombert, J. M. (2016). SPARC: Mass Models for 175 Disk Galaxies with Spitzer Photometry and Accurate Rotation Curves. Astronomical Journal, 152, 157.
LIGO Scientific Collaboration & Virgo Collaboration (2016). Observation of Gravitational Waves from a Binary Black Hole Merger. Physical Review Letters, 116, 061102.
Oldham, K. B., & Spanier, J. (1974). The Fractional Calculus. Academic Press.
Grünwald, A. (1967). Über begrenzte Derivationen und deren Anwendung. Z. Angew. Math. Phys., 18, 609–620.

# ============================================================================== # MONAD-FIELD v7.6 – FINAL SPARC NULL-MODEL SELECTION ENGINE # Fixed‑width parser + NFW comparison + BIC + safe plotting # ============================================================================== import numpy as np import matplotlib.pyplot as plt from scipy.optimize import curve_fit from scipy.stats import pearsonr import os # ------------------------------------------------------------------------------ # 1. Proven fixed‑width SPARC parser (Gemini version – works on your file) # ------------------------------------------------------------------------------ def parse_sparc_mrt_fixed_width(filename): galaxies = [] with open(filename, 'r') as f: for line in f: if len(line) < 60 or line.startswith('#') or line.startswith('Bytes'): continue try: name = line[0:11].strip() if not name or name.replace('.', '', 1).isdigit(): continue r = float(line[19:25]) # kpc v_obs = float(line[26:32]) # km/s e_v = float(line[33:38]) # km/s v_gas = float(line[39:45]) # km/s v_disk = float(line[46:52]) # km/s v_bul = float(line[53:59]) # km/s # find or create galaxy gal = next((g for g in galaxies if g['name'] == name), None) if gal is None: gal = {'name': name, 'R': [], 'Vobs': [], 'Verr': [], 'Vgas': [], 'Vdisk': [], 'Vbul': []} galaxies.append(gal) gal['R'].append(r) gal['Vobs'].append(v_obs) gal['Verr'].append(e_v) gal['Vgas'].append(v_gas) gal['Vdisk'].append(v_disk) gal['Vbul'].append(v_bul) except (ValueError, IndexError): continue # convert to arrays and keep only galaxies with at least 6 points result = [] for gal in galaxies: if len(gal['R']) >= 6: for key in ['R', 'Vobs', 'Verr', 'Vgas', 'Vdisk', 'Vbul']: gal[key] = np.array(gal[key]) result.append(gal) print(f"Loaded {len(result)} galaxies using fixed‑width parser.") return result def generate_synthetic_catalog(): """Fallback synthetic catalog if file missing.""" np.random.seed(42) mock = {} for i in range(106): name = f"Gal_{i:03d}" N = np.random.randint(12, 35) R = np.linspace(0.5, np.random.uniform(8, 25), N) Vb = 130 * (1 - np.exp(-R/4)) Vx = 40 * (R/1.0)**0.2919 Vobs = np.sqrt(Vb**2 + Vx**2) Verr = np.random.uniform(1.5, 6, N) Vobs_noisy = np.maximum(Vobs + np.random.normal(0, Verr), 10) mock[name] = { 'R': R, 'Vobs': Vobs_noisy, 'Verr': Verr, 'Vgas': Vb*0.4, 'Vdisk': Vb*0.6, 'Vbul': np.zeros(N) } print("Using synthetic catalog (106 galaxies).") return mock # ------------------------------------------------------------------------------ # 2. Model functions # ------------------------------------------------------------------------------ def fcmfd_power_law(R, V0, gamma, R0=1.0): return (V0**2) * (R/R0)**(2*gamma) def nfw_profile(R, Vc, c): x = R / np.max(R) x = np.maximum(x, 1e-9) num = np.log(1 + c*x) - (c*x)/(1 + c*x) den = np.log(1 + c) - c/(1 + c) return (Vc**2) * (num / den) / x # ------------------------------------------------------------------------------ # 3. Single galaxy evaluation # ------------------------------------------------------------------------------ def evaluate_galaxy(gal, outer_frac=0.5, sigma_baryon=0.10): R = gal['R'] Vobs = gal['Vobs'] Verr = gal['Verr'] Vgas = gal['Vgas'] Vdisk = gal['Vdisk'] Vbul = gal['Vbul'] Vmax = np.max(Vobs) # baryonic quadrature Vbar2 = Vgas**2 + Vdisk**2 + Vbul**2 Vobs2 = Vobs**2 deltaV2 = Vobs2 - Vbar2 # outer mask mask = (R > outer_frac * np.max(R)) & (deltaV2 > 0) & np.isfinite(deltaV2) if np.sum(mask) < 6: return None R_fit = R[mask] delta = deltaV2[mask] Vobs_fit = Vobs[mask] Verr_fit = Verr[mask] Vbar_fit = np.sqrt(Vbar2[mask]) # error propagation sigma_delta = np.sqrt((2*Vobs_fit*Verr_fit)**2 + (2*sigma_baryon*Vbar_fit**2)**2) safe_delta = np.maximum(delta, 1e-6) sigma_ln = sigma_delta / safe_delta weights = 1.0 / (sigma_ln**2 + 1e-12) # ---- FCMFD power law (log-linear) ---- logR = np.log(R_fit) log_delta = np.log(safe_delta) try: slope, intercept = np.polyfit(logR, log_delta, 1, w=np.sqrt(weights)) gamma = slope / 2.0 V0 = np.sqrt(np.exp(intercept)) y_pred_pl = fcmfd_power_law(R_fit, V0, gamma) chi2_pl = np.sum(((delta - y_pred_pl)/sigma_delta)**2) bic_pl = chi2_pl + 2.0 * np.log(len(R_fit)) except: return None # ---- NFW fit ---- try: p0 = [np.max(Vobs_fit), 12.0] bounds = ([0.0, 0.1], [600.0, 150.0]) popt, _ = curve_fit(nfw_profile, R_fit, delta, sigma=sigma_delta, absolute_sigma=True, p0=p0, bounds=bounds, maxfev=5000) Vc, c = popt y_pred_nfw = nfw_profile(R_fit, Vc, c) chi2_nfw = np.sum(((delta - y_pred_nfw)/sigma_delta)**2) bic_nfw = chi2_nfw + 2.0 * np.log(len(R_fit)) except: return None return { 'gamma': gamma, 'Vmax': Vmax, 'bic_pl': bic_pl, 'bic_nfw': bic_nfw, 'chi2_pl': chi2_pl, 'chi2_nfw': chi2_nfw, 'n_points': len(R_fit) } # ------------------------------------------------------------------------------ # 4. Main pipeline # ------------------------------------------------------------------------------ def main(): filename = "MassModels_Lelli2016c.mrt" if os.path.exists(filename) and os.path.getsize(filename) > 1000: galaxies = parse_sparc_mrt_fixed_width(filename) else: print(f"File {filename} not found or too small. Using synthetic data.") galaxies_data = generate_synthetic_catalog() galaxies = [{'name': k, **v} for k, v in galaxies_data.items()] results = [] for gal in galaxies: res = evaluate_galaxy(gal) if res is not None: results.append(res) if not results: print("No successful fits. Exiting.") return gammas = np.array([r['gamma'] for r in results]) Vmaxs = np.array([r['Vmax'] for r in results]) deltas = np.array([r['bic_pl'] - r['bic_nfw'] for r in results]) # statistics r_val, p_val = pearsonr(Vmaxs, gammas) fcmfd_pref = np.sum(deltas < -2.0) nfw_pref = np.sum(deltas > 2.0) tie = np.sum(np.abs(deltas) <= 2.0) print("\n" + "="*80) print("FINAL CONSTITUTIVE SELECTION REPORT") print("="*80) print(f"Successful fits: {len(gammas)}") print(f"Gamma: median = {np.median(gammas):.4f} mean = {np.mean(gammas):.4f} +/- {np.std(gammas):.4f}") print(f"Gamma range: {np.min(gammas):.4f} – {np.max(gammas):.4f}") print(f"Correlation gamma vs Vmax: r = {r_val:.4f} (p = {p_val:.4f})") print(f"Model preference (ΔBIC = BIC_FCMFD - BIC_NFW):") print(f" FCMFD preferred (ΔBIC < -2): {fcmfd_pref}") print(f" Indistinguishable (|ΔBIC| <= 2): {tie}") print(f" NFW preferred (ΔBIC > 2): {nfw_pref}") print("="*80) # Plot (safe labels – no LaTeX crashes) fig, axes = plt.subplots(1, 2, figsize=(14, 5)) # scatter: gamma vs Vmax colored by ΔBIC sc = axes[0].scatter(Vmaxs, gammas, c=deltas, cmap='coolwarm', edgecolor='black', alpha=0.8, s=50, vmin=-10, vmax=10) axes[0].axhline(0.292, color='gray', linestyle='--', label='FCMFD simulation (beta=0.292)') axes[0].set_xlabel('Vmax (km/s) - mass proxy') axes[0].set_ylabel('Gamma (spatial scaling exponent)') axes[0].set_title('Galaxy scaling vs mass') axes[0].legend() axes[0].grid(alpha=0.3) cbar = fig.colorbar(sc, ax=axes[0]) cbar.set_label('Delta BIC = BIC_FCMFD - BIC_NFW') # histogram of ΔBIC axes[1].hist(deltas, bins=20, color='purple', alpha=0.7, edgecolor='black', range=(-15,15)) axes[1].axvline(-2, color='crimson', linestyle='--', label='FCMFD favoured (ΔBIC < -2)') axes[1].axvline(2, color='royalblue', linestyle='--', label='NFW favoured (ΔBIC > 2)') axes[1].set_xlabel('Delta BIC') axes[1].set_ylabel('Number of galaxies') axes[1].set_title('Model selection distribution') axes[1].legend() axes[1].grid(alpha=0.3) plt.tight_layout() plt.show() if __name__ == "__main__": main()

Conspiranon

FCMFD: From Galactic Halos to Binary Mergers — A Constitutive Framework without Dark Matter