The file src/32-dark/deep_mond.py of the validation package of the paper 32-dark, rendered with an anchor per line; the file itself is beside this page.
1# Construction 2 (LABELLED continuum confrontation): the two-chart Gauss law.2# Conserved flux per shell = enclosed baryonic source (Gauss, exact, both charts).3# position chart: registered force = flux -> g_eff = g_bar (Newtonian)4# spectral chart: registered force = Born amplitude -> g_eff = sqrt(g_bar*a0) (deep)5# crossover at the floor a0 = c H0 / 2pi (the outward scale horizon). The standard6# rotation-angle weighting gives the interpolation g_eff = g_bar/(1-exp(-sqrt(g_bar/a0))).7import numpy as np8c=2.998e8; G=6.674e-11; Msun=1.989e30; Mpc=3.0857e22; kpc=Mpc/10009H0=67.4*1000/Mpc; a0=c*H0/(2*np.pi) # H0 = 67.4 km/s/Mpc, the rate the octant lemma entails (14-entr); a0 = 1.042e-1010print(f"a0 = c H0/2pi = {a0:.3e} m/s^2 (RAR fit 1.20e-10)")1112def g_eff(g_bar): # interpolation (deep slope 1/2, Newtonian limit)13 return g_bar/(1-np.exp(-np.sqrt(g_bar/a0)))1415# RAR deep slope16gb=np.logspace(-13,-8,200); go=g_eff(gb)17deep=gb<1e-1218slope=np.polyfit(np.log(gb[deep]),np.log(go[deep]),1)[0]19print(f"RAR deep-regime slope d ln g_obs/d ln g_bar = {slope:.3f} (target 0.5)")2021# BTFR: flat-velocity vs mass, point-mass deep limit v^4 = G M a022M=np.logspace(8,12,40)*Msun23vflat=(G*M*a0)**0.2524bt=np.polyfit(np.log(M),np.log(vflat),1)[0]25print(f"BTFR slope d ln v_flat/d ln M = {bt:.3f} (target 0.25 -> v^4 ~ M)")26print(f" M=5e10 Msun -> v_flat = {(G*5e10*Msun*a0)**0.25/1000:.0f} km/s")2728# Rotation curve for an exponential disk, Sigma0 with scale length Rd29def disk_curve(Mdisk, Rd):30 r=np.linspace(0.1,30,300)*kpc31 x=r/(2*Rd)32 from scipy.special import i0,i1,k0,k133 # exact exponential-disk Newtonian v_bar^234 vb2=(G*Mdisk/Rd)*x**2*(i0(x)*k0(x)-i1(x)*k1(x))*235 gbar=vb2/r36 gobs=g_eff(gbar)37 vobs=np.sqrt(gobs*r)38 return r/kpc, np.sqrt(vb2)/1000, vobs/100039try:40 import scipy41 r,vb,vo=disk_curve(6e10*Msun,3*kpc)42 print(f"disk: baryonic v peaks {vb.max():.0f} km/s -> observed flattens at {vo[-50:].mean():.0f} km/s")43 have_disk=True44except Exception as e:45 have_disk=False; print("scipy unavailable, skipping disk curve panel")4647import matplotlib; matplotlib.use('Agg'); import matplotlib.pyplot as plt48n=3 if have_disk else 249fig,ax=plt.subplots(1,n,figsize=(4.2*n,3.6))50ax[0].loglog(gb,go,'b',label='two-chart interpolation')51ax[0].loglog(gb,gb,'k--',lw=0.8,label='Newtonian (slope 1)')52ax[0].loglog(gb,np.sqrt(gb*a0),'g-',lw=1,label='deep $\\sqrt{g_{bar}a_0}$ (slope 1/2)')53ax[0].axvline(a0,color='gray',ls=':',lw=0.8); ax[0].set_xlabel('$g_{bar}$'); ax[0].set_ylabel('$g_{obs}$')54ax[0].set_title('Radial acceleration relation'); ax[0].legend(fontsize=6)55ax[1].loglog(M/Msun,vflat/1000,'r'); ax[1].set_xlabel('$M_{bar}/M_\\odot$'); ax[1].set_ylabel('$v_{flat}$ [km/s]')56ax[1].set_title(f'Baryonic Tully-Fisher (slope {bt:.2f})')57if have_disk:58 ax[2].plot(r,vb,'k--',label='baryonic'); ax[2].plot(r,vo,'b',label='observed (two-chart)')59 ax[2].set_xlabel('r [kpc]'); ax[2].set_ylabel('v [km/s]'); ax[2].set_title('Rotation curve (exp. disk)'); ax[2].legend(fontsize=7)60fig.tight_layout(); fig.savefig('cons2_rar_btfr.pdf'); print("figure saved")61