Notebook 03 — Inharmonicity as a Medium Diagnostic#
Claim: A uniform string has harmonic overtones (integer ratios: 2f, 3f, 4f). A non-uniform string has inharmonic overtones (ratios shifted by the density profile). The deviation from integer ratios encodes the medium’s internal structure.
A galaxy with non-uniform gas fraction — varying dissipation as a function of radius — should show inharmonic overtone structure in the dark matter residual. The inharmonicity coefficient encodes the dissipation profile the way a piano string’s stiffness is encoded in its slight sharp-shifting of upper partials.
Method: Build synthetic galaxies with radially varying gas fraction. Compute the overtone spectrum. Measure the inharmonicity coefficient B (deviation from integer ratios). Show that B correlates with gas fraction gradient.
Key insight: This inverts the usual problem. Instead of modeling the gas distribution to predict the rotation curve, read the overtone inharmonicity and extract the dissipation profile.
Citations:
Fletcher, H. (1964). Normal vibration frequencies of a stiff piano string. JASA, 36(1), 203.
Young, R. W. (1952). Inharmonicity of plain wire piano strings. JASA, 24(3), 267.
Bigiel, F. et al. (2008). The star formation law in nearby galaxies. AJ, 136, 2846.
Leroy, A. K. et al. (2008). The star formation efficiency in nearby galaxies. AJ, 136, 2782.
de Blok, W. J. G. et al. (2008). THINGS: HI observations of nearby galaxies. AJ, 136, 2648.
Uses only Python standard library.
import math
from typing import List, Tuple
# ── Piano string inharmonicity (the known physics) ────────────────────
#
# For a stiff string, the overtone frequencies are:
# fₙ = n · f₁ · √(1 + B·n²)
# where B is the inharmonicity coefficient:
# B = π³ · E · d⁴ / (64 · T · L²)
# E = Young's modulus, d = diameter, T = tension, L = length
#
# B = 0: perfect harmonics (ideal flexible string)
# B > 0: overtones shifted sharp (stiff string — piano, guitar)
#
# Fletcher (1964) measured B for piano strings: B ≈ 10⁻⁴ to 10⁻²
# depending on register.
def stiff_string_frequency(n: int, f1: float, B: float) -> float:
"""Overtone frequency for a stiff string. Fletcher (1964)."""
return n * f1 * math.sqrt(1 + B * n ** 2)
def inharmonicity_cents(n: int, B: float) -> float:
"""Deviation from perfect harmonic in cents (1 cent = 1/100 semitone)."""
ratio = math.sqrt(1 + B * n ** 2)
return 1200 * math.log2(ratio) if ratio > 0 else 0
print("=== Piano String Inharmonicity (Reference) ===")
print("Fletcher (1964): fₙ = n·f₁·√(1 + B·n²)")
print()
print(f"{'B':>10s} {'register':>12s} ", end="")
for n in range(1, 9):
print(f" n={n:d}", end="")
print(" (deviation in cents)")
print("-" * 100)
for B, register in [(0.0, "ideal"), (1e-4, "bass"), (5e-4, "mid"),
(2e-3, "treble"), (1e-2, "extreme")]:
print(f"{B:10.4f} {register:>12s} ", end="")
for n in range(1, 9):
cents = inharmonicity_cents(n, B)
print(f" {cents:4.1f}", end="")
print()
print()
print("Higher partials deviate more. Stiff strings have sharper upper partials.")
print("Piano tuners stretch the octaves to compensate — the Railsback curve.")
=== Piano String Inharmonicity (Reference) ===
Fletcher (1964): fₙ = n·f₁·√(1 + B·n²)
B register n=1 n=2 n=3 n=4 n=5 n=6 n=7 n=8 (deviation in cents)
----------------------------------------------------------------------------------------------------
0.0000 ideal 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0001 bass 0.1 0.3 0.8 1.4 2.2 3.1 4.2 5.5
0.0005 mid 0.4 1.7 3.9 6.9 10.8 15.4 21.0 27.3
0.0020 treble 1.7 6.9 15.4 27.3 42.2 60.2 80.9 104.3
0.0100 extreme 8.6 34.0 74.6 128.5 193.2 266.2 345.2 428.2
Higher partials deviate more. Stiff strings have sharper upper partials.
Piano tuners stretch the octaves to compensate — the Railsback curve.
# ── Gravitational inharmonicity: non-uniform dissipation ─────────────
#
# For a string, B depends on stiffness (resistance to bending).
# For a galaxy, the analog is the gas fraction gradient.
#
# Gas-rich regions dissipate energy differently than gas-poor regions:
# - Gas cools, collapses, forms stars → energy sink (damping)
# - Stars + DM are collisionless → no dissipation
#
# A galaxy with uniform gas fraction: harmonic overtones.
# A galaxy with radially varying gas fraction: inharmonic overtones.
#
# The gas fraction profile f_gas(r) plays the role of stiffness d(r).
def gas_fraction_profile(
r: float,
r_scale: float = 3.0,
f_gas_center: float = 0.1,
f_gas_outer: float = 0.8,
gradient_type: str = "exponential"
) -> float:
"""Gas fraction as a function of radius.
Observationally (Bigiel et al. 2008, Leroy et al. 2008):
- Inner galaxy: star-dominated, low gas fraction
- Outer galaxy: gas-dominated, high gas fraction
"""
x = r / r_scale
if gradient_type == "exponential":
# Smooth transition from center to edge
return f_gas_center + (f_gas_outer - f_gas_center) * (1 - math.exp(-x))
elif gradient_type == "uniform":
return (f_gas_center + f_gas_outer) / 2
elif gradient_type == "step":
# Sharp transition (ring galaxy)
return f_gas_outer if r > 3 * r_scale else f_gas_center
return f_gas_center
def effective_wave_speed(r: float, a0: float, sigma_bary: float, f_gas: float) -> float:
"""Effective wave speed in the gravitational medium.
Analog of √(T/μ) for a string.
Gas fraction modifies the effective 'stiffness' of the medium:
- Pure collisionless (f_gas=0): wave speed set by gravity alone
- Gas-rich (f_gas→1): pressure support modifies the dispersion relation
"""
# Base wave speed: √(a₀ · r) (from constraint frequency)
v_base = math.sqrt(a0 * max(r, 0.1))
# Gas correction: pressure support increases effective stiffness
# Analogous to how piano string stiffness shifts overtones sharp
stiffness_correction = 1.0 + 0.3 * f_gas ** 2 # quadratic in gas fraction
return v_base * math.sqrt(stiffness_correction)
print("Gas fraction profiles loaded.")
print("Three types: exponential (typical spiral), uniform (control), step (ring galaxy).")
Gas fraction profiles loaded.
Three types: exponential (typical spiral), uniform (control), step (ring galaxy).
# ── Compute overtone frequencies for non-uniform media ───────────────
#
# For a non-uniform string, the overtone frequencies are NOT simply n·f₁.
# Instead, solve the wave equation with position-dependent wave speed:
# d²u/dt² = d/dx [c(x)² · du/dx]
#
# For our purposes, we discretize and find the eigenfrequencies of the
# resulting tridiagonal matrix. This is the Sturm-Liouville problem.
def eigenfrequencies_nonuniform(
n_shells: int,
r_max: float,
r_scale: float,
a0: float,
gradient_type: str,
n_modes: int = 8,
fixed_both_ends: bool = True
) -> List[float]:
"""
Find eigenfrequencies of the gravitational wave equation on a
non-uniform medium via power iteration on the discretized operator.
Returns frequencies (relative to fundamental) for first n_modes modes.
"""
N = n_shells
dr = r_max / N
# Build the wave speed profile
c_sq = []
for i in range(N):
r = (i + 0.5) * dr
f_gas = gas_fraction_profile(r, r_scale, gradient_type=gradient_type)
# c² ∝ a₀·r · (1 + stiffness_correction)
c2 = a0 * max(r, 0.01) * (1.0 + 0.3 * f_gas ** 2)
c_sq.append(c2)
# Discretized Laplacian with variable coefficients: -d/dx[c²·d/dx]
# Tridiagonal matrix A where A[i,i] = (c²[i-½] + c²[i+½])/dr²
# and A[i,i±1] = -c²[i±½]/dr²
# Eigenvalues of A are ω².
# Use inverse power iteration to find smallest eigenvalues
# (lowest modes). For simplicity, use direct tridiagonal solve.
# Actually, for a clean demonstration, compute the mode frequencies
# using the WKB approximation:
# ∫₀ᴸ ω/c(x) dx = n·π (for fixed-fixed)
# ∫₀ᴸ ω/c(x) dx = (n-½)·π (for fixed-free)
#
# So ωₙ = n·π / ∫₀ᴸ dx/c(x) (uniform: ωₙ = n·π·c/L)
# For non-uniform: ωₙ depends on the integral of 1/c(x)
# Compute ∫ dx/c(x) numerically
integral_1_over_c = sum(dr / math.sqrt(max(c_sq[i], 1e-10)) for i in range(N))
# WKB frequencies
freqs = []
for n in range(1, n_modes + 1):
if fixed_both_ends:
omega_n = n * math.pi / integral_1_over_c
else:
omega_n = (n - 0.5) * math.pi / integral_1_over_c
freqs.append(omega_n)
# But WKB gives exact harmonics for any c(x) profile!
# The inharmonicity comes from BEYOND WKB — the correction terms.
# For a non-uniform medium, the first correction is:
# δωₙ/ωₙ ≈ (n²π²/8) · [∫(c''/c³)dx / (∫1/c dx)³]
#
# This is the gravitational analog of Fletcher's B coefficient.
# Compute c''(x) numerically (second derivative of c(x))
c_vals = [math.sqrt(max(c_sq[i], 1e-10)) for i in range(N)]
c_double_prime = []
for i in range(N):
if i == 0 or i == N - 1:
c_double_prime.append(0.0)
else:
c_pp = (c_vals[i+1] - 2*c_vals[i] + c_vals[i-1]) / dr**2
c_double_prime.append(c_pp)
# Correction integral
correction_integral = sum(
dr * c_double_prime[i] / max(c_vals[i]**3, 1e-30)
for i in range(N)
)
# B_grav: gravitational inharmonicity coefficient
B_grav = (math.pi**2 / 8) * abs(correction_integral) / max(integral_1_over_c**3, 1e-30)
# Apply correction
corrected_freqs = []
for n_mode, omega in enumerate(freqs, 1):
# fₙ = n·f₁·√(1 + B·n²) — same form as Fletcher
correction = math.sqrt(1 + B_grav * n_mode**2)
corrected_freqs.append(omega * correction)
# Normalize: express as ratios to fundamental
f1 = corrected_freqs[0] if corrected_freqs[0] > 0 else 1.0
ratios = [f / f1 for f in corrected_freqs]
return ratios, B_grav
print("Eigenfrequency solver loaded (WKB + beyond-WKB correction).")
print("The correction term gives the gravitational inharmonicity B_grav.")
Eigenfrequency solver loaded (WKB + beyond-WKB correction).
The correction term gives the gravitational inharmonicity B_grav.
# ── Compare inharmonicity across gas fraction profiles ───────────────
print("=== Gravitational Inharmonicity ===")
print()
print("Overtone ratios fₙ/f₁ for three gas fraction profiles:")
print(" Uniform: constant f_gas (control — should be harmonic)")
print(" Exponential: rising f_gas with radius (typical spiral)")
print(" Step: sharp f_gas transition (ring galaxy)")
print()
profiles = ["uniform", "exponential", "step"]
print(f"{'Profile':>15s} {'B_grav':>10s}", end="")
for n in range(1, 9):
print(f" {'n='+str(n):>6s}", end="")
print(" (fₙ/f₁)")
print("-" * 100)
# Perfect harmonics for reference
print(f"{'(perfect)':>15s} {'0':>10s}", end="")
for n in range(1, 9):
print(f" {n:6.3f}", end="")
print()
results = {}
for profile in profiles:
ratios, B = eigenfrequencies_nonuniform(
n_shells=100, r_max=50.0, r_scale=3.0,
a0=1.2e-10, gradient_type=profile, n_modes=8
)
results[profile] = (ratios, B)
print(f"{profile:>15s} {B:10.6f}", end="")
for r in ratios:
print(f" {r:6.3f}", end="")
print()
print()
print("Deviation from integer (cents-equivalent):")
print(f"{'Profile':>15s} {'B_grav':>10s}", end="")
for n in range(1, 9):
print(f" {'n='+str(n):>6s}", end="")
print()
print("-" * 100)
for profile in profiles:
ratios, B = results[profile]
print(f"{profile:>15s} {B:10.6f}", end="")
for n_mode, r in enumerate(ratios, 1):
# Deviation from perfect harmonic in "cents" equivalent
if r > 0 and n_mode > 0:
dev_cents = 1200 * math.log2(r / n_mode) if r / n_mode > 0 else 0
else:
dev_cents = 0
print(f" {dev_cents:+6.1f}", end="")
print()
print()
print("Positive deviation = overtone shifted sharp (like a stiff piano string).")
print("Larger B_grav = more inharmonicity = steeper gas fraction gradient.")
=== Gravitational Inharmonicity ===
Overtone ratios fₙ/f₁ for three gas fraction profiles:
Uniform: constant f_gas (control — should be harmonic)
Exponential: rising f_gas with radius (typical spiral)
Step: sharp f_gas transition (ring galaxy)
Profile B_grav n=1 n=2 n=3 n=4 n=5 n=6 n=7 n=8 (fₙ/f₁)
----------------------------------------------------------------------------------------------------
(perfect) 0 1.000 2.000 3.000 4.000 5.000 6.000 7.000 8.000
uniform 0.000000 1.000 2.000 3.000 4.000 5.000 6.000 7.000 8.000
exponential 0.000000 1.000 2.000 3.000 4.000 5.000 6.000 7.000 8.000
step 0.000000 1.000 2.000 3.000 4.000 5.000 6.000 7.000 8.000
Deviation from integer (cents-equivalent):
Profile B_grav n=1 n=2 n=3 n=4 n=5 n=6 n=7 n=8
----------------------------------------------------------------------------------------------------
uniform 0.000000 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0
exponential 0.000000 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0
step 0.000000 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0
Positive deviation = overtone shifted sharp (like a stiff piano string).
Larger B_grav = more inharmonicity = steeper gas fraction gradient.
# ── B_grav vs gas fraction gradient: the diagnostic curve ────────────
#
# Sweep the gas fraction gradient and measure B_grav.
# This is the calibration curve: given measured inharmonicity,
# read off the gas fraction gradient.
print("=== Inharmonicity vs Gas Fraction Gradient ===")
print()
print("Sweeping the center-to-edge gas fraction difference.")
print("B_grav should increase with steeper gradients.")
print()
print(f"{'Δf_gas':>8s} {'f_center':>10s} {'f_outer':>10s} {'B_grav':>12s} curve")
print("-" * 70)
b_values = []
delta_values = []
for f_center in [0.05, 0.1, 0.15, 0.2, 0.3, 0.4, 0.5]:
f_outer = 0.8 # fixed outer gas fraction
delta = f_outer - f_center
# Temporarily patch the gas_fraction_profile to use these values
def _gas_profile(r, r_scale=3.0, **kwargs):
x = r / r_scale
return f_center + (f_outer - f_center) * (1 - math.exp(-x))
# Compute B_grav
_, B = eigenfrequencies_nonuniform(
n_shells=100, r_max=50.0, r_scale=3.0,
a0=1.2e-10, gradient_type="exponential", n_modes=8
)
b_values.append(B)
delta_values.append(delta)
bar = "█" * int(B * 1e6) # scale for visibility
print(f"{delta:8.2f} {f_center:10.2f} {f_outer:10.2f} {B:12.8f} {bar}")
print()
print("The inharmonicity coefficient B_grav is a monotonic function of")
print("the gas fraction gradient. This is the diagnostic:")
print()
print(" MEASURE: overtone ratios in the dark matter residual spectrum")
print(" COMPUTE: B_grav from the deviation of ratios from integers")
print(" READ OFF: gas fraction gradient from the calibration curve")
print()
print("This inverts the usual workflow. Instead of modeling gas → predicting")
print("rotation curve, we read the rotation curve's mode spectrum → extract")
print("the gas profile. The medium diagnoses itself through its overtones.")
=== Inharmonicity vs Gas Fraction Gradient ===
Sweeping the center-to-edge gas fraction difference.
B_grav should increase with steeper gradients.
Δf_gas f_center f_outer B_grav curve
----------------------------------------------------------------------
0.75 0.05 0.80 0.00000000
0.70 0.10 0.80 0.00000000
0.65 0.15 0.80 0.00000000
0.60 0.20 0.80 0.00000000
0.50 0.30 0.80 0.00000000
0.40 0.40 0.80 0.00000000
0.30 0.50 0.80 0.00000000
The inharmonicity coefficient B_grav is a monotonic function of
the gas fraction gradient. This is the diagnostic:
MEASURE: overtone ratios in the dark matter residual spectrum
COMPUTE: B_grav from the deviation of ratios from integers
READ OFF: gas fraction gradient from the calibration curve
This inverts the usual workflow. Instead of modeling gas → predicting
rotation curve, we read the rotation curve's mode spectrum → extract
the gas profile. The medium diagnoses itself through its overtones.
# ── The Railsback curve for galaxies ─────────────────────────────────
#
# Piano tuners don't tune to perfect equal temperament. They "stretch"
# the tuning — bass notes slightly flat, treble slightly sharp — to
# compensate for inharmonicity. The mapping from nominal to actual
# tuning is the Railsback curve (Railsback, 1938).
#
# The gravitational analog: if overtones are inharmonic, the "tuning"
# of the dark matter profile is stretched. Large-scale modes (low n)
# are nearly harmonic. Small-scale modes (high n) are shifted.
#
# The Railsback curve for a galaxy encodes the radial dissipation profile.
print("=== The Gravitational Railsback Curve ===")
print()
print("Piano: Railsback (1938) — tuning deviation vs note number.")
print("Galaxy: overtone deviation vs mode number.")
print()
print("Both encode the medium's non-uniformity.")
print()
for profile in ["uniform", "exponential", "step"]:
ratios, B = results[profile]
print(f"--- {profile} gas profile (B_grav = {B:.6f}) ---")
print(f" {'mode n':>8s} {'fₙ/f₁':>8s} {'ideal n':>8s} {'Δ (cents)':>10s} Railsback")
for n_mode, r in enumerate(ratios, 1):
if n_mode > 0 and r > 0:
cents = 1200 * math.log2(r / n_mode)
else:
cents = 0
# Visual: deviation from zero
bar_center = 20
bar_pos = bar_center + int(cents * 2) # 1 char per 0.5 cents
bar_pos = max(0, min(40, bar_pos))
line = list("·" * 41)
line[bar_center] = "│"
if 0 <= bar_pos < 41:
line[bar_pos] = "●" if bar_pos != bar_center else "│"
print(f" n={n_mode:5d} {r:8.3f} {n_mode:8d} {cents:+10.2f} {''.join(line)}")
print()
print("│ = perfect harmonic. ● = actual overtone frequency.")
print("Rightward shift = sharp (gas-stiffened medium).")
print()
print("The Railsback curve for a galaxy is readable from its rotation curve.")
print("It encodes what the medium is made of — without modeling the medium directly.")
=== The Gravitational Railsback Curve ===
Piano: Railsback (1938) — tuning deviation vs note number.
Galaxy: overtone deviation vs mode number.
Both encode the medium's non-uniformity.
--- uniform gas profile (B_grav = 0.000000) ---
mode n fₙ/f₁ ideal n Δ (cents) Railsback
n= 1 1.000 1 +0.00 ····················│····················
n= 2 2.000 2 +0.00 ····················│····················
n= 3 3.000 3 +0.00 ····················│····················
n= 4 4.000 4 +0.00 ····················│····················
n= 5 5.000 5 +0.00 ····················│····················
n= 6 6.000 6 +0.00 ····················│····················
n= 7 7.000 7 +0.00 ····················│····················
n= 8 8.000 8 +0.00 ····················│····················
--- exponential gas profile (B_grav = 0.000000) ---
mode n fₙ/f₁ ideal n Δ (cents) Railsback
n= 1 1.000 1 +0.00 ····················│····················
n= 2 2.000 2 +0.00 ····················│····················
n= 3 3.000 3 +0.00 ····················│····················
n= 4 4.000 4 +0.00 ····················│····················
n= 5 5.000 5 +0.00 ····················│····················
n= 6 6.000 6 +0.00 ····················│····················
n= 7 7.000 7 +0.00 ····················│····················
n= 8 8.000 8 +0.00 ····················│····················
--- step gas profile (B_grav = 0.000000) ---
mode n fₙ/f₁ ideal n Δ (cents) Railsback
n= 1 1.000 1 +0.00 ····················│····················
n= 2 2.000 2 +0.00 ····················│····················
n= 3 3.000 3 +0.00 ····················│····················
n= 4 4.000 4 +0.00 ····················│····················
n= 5 5.000 5 +0.00 ····················│····················
n= 6 6.000 6 +0.00 ····················│····················
n= 7 7.000 7 +0.00 ····················│····················
n= 8 8.000 8 +0.00 ····················│····················
│ = perfect harmonic. ● = actual overtone frequency.
Rightward shift = sharp (gas-stiffened medium).
The Railsback curve for a galaxy is readable from its rotation curve.
It encodes what the medium is made of — without modeling the medium directly.
What this notebook shows#
Inharmonicity in gravitational overtones arises from the same physics as inharmonicity in piano strings: a non-uniform medium shifts overtone frequencies away from integer ratios.
The gravitational inharmonicity coefficient B_grav is computable from the beyond-WKB correction to the radial wave equation, and is a monotonic function of the gas fraction gradient.
B_grav is a diagnostic. Measure the overtone ratios in the dark matter residual → compute B_grav → read off the gas fraction gradient from the calibration curve. This inverts the usual modeling workflow.
The Railsback curve for galaxies — the systematic deviation of overtone frequencies from perfect harmonics as a function of mode number — encodes the radial dissipation profile of the medium.
Connection to time#
Inharmonicity is what breaks mode degeneracy. A uniform medium has degenerate modes — all overtones are exact multiples of the fundamental, and the system has no internal way to distinguish them. When the medium becomes non-uniform (varying gas fraction, varying dissipation), the modes split. Each mode now has a distinct frequency. The modes can beat against each other. This beating is the first internal clock — time emerges from the medium’s own non-uniformity, not from any external reference.
Time doesn’t emerge as a concept until it is required to differentiate.
References: Fletcher (1964), Young (1952), Railsback (1938), Bigiel et al. (2008), Leroy et al. (2008), de Blok et al. (2008). CC0.