The Spacetime Resonance Theorem and the Origin of Dark Mass

The Spacetime Resonance Theorem and the Origin of Dark Mass Bryan Hennelly ∗1,2,3 1Department of Electronic Engineering, Maynooth University, Maynooth, Co. Kildare, Ireland 2Department of Computer Science, Maynooth University, Maynooth, Co. Kildare, Ireland 3The Hamilton Institute, Maynooth University, Maynooth, Co. Kildare, Ireland 19 November 2025 Abstract In linearised general relativity the gravitational field is determined by a retarded Green function, but in most applications its “causal tail” is discarded. When this tail is kept and combined with the linear dampedoscillator susceptibility we introduce for baryonic matter (mathematically identical to the response used in acoustical and optical scattering theory), this yields an effective–medium description of the weak–field gravitational response. Coherent, time–varying baryonic motion drives a collective curvature field whose delayed rescattering produces a small but persistent cycle–averaged component. This stored curvature energy appears as an effective dark–mass density and reproduces the observed dark–mass phenomenology. No new fields, particles, or modifications of GR are introduced: all dynamics arise from standard linearised gravity once the usual quasistatic truncation of the retarded kernel is avoided. The resulting causal–resonant framework reproduces the effective mass distributions of spiral galaxies, ellipticals, clusters, and filaments; yields natural “local” and “homogeneous” silence; and remains fully consistent with Solar–System tests and with the CMB. In this picture, dark–mass behaviour emerges as curvature energy stored by the collective resonant response of baryonic structure. ∗Email: [email protected] 1 Contents 1 Introduction 5 2 Spacetime Resonance Theorem 6 2.1 The Fundamental Oscillator–Medium Theorem . . . . . . . . . . . 6 2.2 The Spacetime Resonance Theorem . . . . . . . . . . . . . . . . . 8 2.3 Corollary: Density–Controlled Restoring Frequency (Motivated by Known Gravitational Modes) . . . . . . . . . . . . . . . . . . . 9 2.4 Physical Interpretation and Preview . . . . . . . . . . . . . . . . . 11 3 Dark Mass 12 4 Resonant Field Dynamics 13 4.1 ChannelSeparation.......................... 13 4.2 Spacetime as an Effective Reversible Medium . . . . . . . . . . . 14 4.3 Collective Mode Dynamics and the Masson . . . . . . . . . . . . . 15 4.4 Linear Spatial Perturbations and the Propagation Law . . . . . . 16 5 Gravitational Scattering Theory 18 5.1 Single–Scatterer Response . . . . . . . . . . . . . . . . . . . . . . 19 5.2 Multiple Scattering and the Foldy–Lax Hierarchy . . . . . . . . . 20 5.3 Effective–Medium and Ensemble Averaging . . . . . . . . . . . . . 21 5.4 Transport and Diffusion Limits . . . . . . . . . . . . . . . . . . . 22 5.5 Unified Interpretation . . . . . . . . . . . . . . . . . . . . . . . . . 23 6 Spatiotemporal Feedback and Resonant Coherence 24 6.1 Constituent Entities & States of the Feedback System . . . . . . . 24 6.2 Feedback Channels and Interactions . . . . . . . . . . . . . . . . . 26 6.3 Temporal Coherence and Mutual Susceptibility . . . . . . . . . . 29 6.3.1 Phase Locking and Mutual Susceptibility . . . . . . . . . . 29 6.3.2 Global Frequency Selection in Inhomogeneous Baryonic Systems ............................ 31 6.3.3 Polarisation of Temporal Coherence . . . . . . . . . . . . . 32 7 Spatial Coherence of the Resonant Channel 33 7.1 Kirchhoff–Sommerfeld Representation . . . . . . . . . . . . . . . . 34 7.2 Effective Potential and Energy Storage in the Resonant Medium . 36 7.3 Single Point Source in a Resonant Channel–1 Background . . . . 37 7.4 Interference of Two Point Sources . . . . . . . . . . . . . . . . . . 39 7.5 Temporal Memory and the Steady–State Reduction . . . . . . . . 44 8 Steady–State (Attractor) Solutions 45 8.1 Spiral Galaxies as Coherent m= 1 (Laguerre–Gaussian) Modes . . . . . . . . . . . . . . . . . . . . . 46 2 8.1.1 Observing the Steady–State Solution . . . . . . . . . . . . 46 8.1.2 The First Moment . . . . . . . . . . . . . . . . . . . . . . 49 8.1.3 The Thin–Disk Model: Solving the 2D Helmholtz Equation RevealstheHalo ....................... 51 8.1.4 Extending the Thin–Disk to Three Dimensions . . . . . . . 53 8.1.5 Gain from Incoherent Foldy–Lax Scattering . . . . . . . . 54 8.1.6 Stored Energy, Dark–Mass Equivalence, and EnvironmentalGain ............................ 55 8.1.7 Numerical Results: 3D Thin–Disk Model . . . . . . . . . . 56 8.1.8 Physical Interpretation . . . . . . . . . . . . . . . . . . . . 60 8.2 Ellipticals as Random Resonant Media . . . . . . . . . . . . . . . 64 8.2.1 Observing the Statistical Steady State . . . . . . . . . . . 64 8.2.2 Complex Baryonic Forcing with Memory . . . . . . . . . . 66 8.2.3 Four–dimensional speckle geometry of the gravitational field 68 8.2.4 Ensemble reduction: masson–speckle summation . . . . . . 69 8.2.5 Transport Mean Free Path and Collective Damping . . . . 72 8.2.6 Core Regularization in the Non–Transport Regime . . . . 73 8.2.7 Gradient–Energy Boost and Final Transport–Regularized HaloLaw ........................... 75 8.2.8 Numerical Results . . . . . . . . . . . . . . . . . . . . . . 76 8.2.9 Halo Formation from Coherent Energy Flux . . . . . . . . 79 8.2.10 Physical Interpretation . . . . . . . . . . . . . . . . . . . . 80 8.3 Clusters as Radiative Transport Extensions of Ellipticals . . . . . 83 8.3.1 Observing the Statistical Transport Equilibrium . . . . . . 83 8.3.2 Mathematical Formulation: Overlapping Elliptical TransportEnvelopes ........................ 86 8.3.3 Numerical Results . . . . . . . . . . . . . . . . . . . . . . 91 8.3.4 Physical Interpretation . . . . . . . . . . . . . . . . . . . . 94 9 Local Silence and the Preservation of Classical Gravity 96 10 Prediction: Absence of Gain in Homogeneous Systems 98 10.1 Three Distinct Regimes for Homogeneous Clouds . . . . . . . . . 98 10.2 Theorem: Homogeneous Media Produce No Resonant Gain . . . . 99 10.3 Interpretation and Connection to Perturbation Theory . . . . . . 100 11 Emergence of Multiple Scattering from a Weakly Inhomogeneous Continuous Medium 100 11.1 A Brief Review of Garvitational Multiple Scatter Theory . . . . . 100 11.2 The Resonant Frequency in a Continuous Medium . . . . . . . . . 102 11.3 Weak Random Potential and the Ballistic Regime . . . . . . . . . 103 11.4 Nonlinear Growth Toward the Transport Threshold . . . . . . . . 103 11.5 Emergent Resonant Patches and Effective Discreteness . . . . . . 104 11.6 Linear Limit and Connection to the CMB Power Spectrum . . . . 104 12 Conclusion 106 3 13 Dedication and Acknowledgement 107 A Appendix: Bessel and Hankel Functions Used in the Main Text 108 A.1 Cylindrical Bessel and Hankel Functions . . . . . . . . . . . . . . 108 A.2 Spherical Bessel Functions . . . . . . . . . . . . . . . . . . . . . . 109 A.3 Physical Interpretation . . . . . . . . . . . . . . . . . . . . . . . . 109 B Appendix: Code for Simulating the 3D Thin–Disk Resonant Halo 110 C Appendix: Code for Simulating Elliptical Galaxies - Transport Theory 119 D Appendix: Code for Simulating Cluster Galaxies 129 4 1 Introduction The empirical success of Newtonian and Einsteinian gravity [ 1 , 2 ] on laboratory, Solar–System, and stellar scales is beyond dispute. Yet across galaxies, clusters, and cosmological environments, the observed gravitational field exceeds the response predicted from baryons alone. Rather than postulating a new particulate component, we revisit a well–known but rarely exploited aspect of general relativity: the full retarded Green–function structure of the linearised field. In the weak–field limit, curvature responds through the causal kernel Gret , whose support extends throughout the past light cone. These interior contributions, the “tails”, are genuine features of GR but are normally discarded when the field is approximated by its instantaneous Poisson limit. Restoring the full retarded response exposes a small but physically meaningful effect: curvature reacts to time–varying matter through a delayed, frequency–dependent, finite-memory response that is already present in the linearised theory. A central point—formalised later by the Spacetime Resonance Theorem (Section 2)—is that we do not introduce an oscillator model for the mass distribution. The oscillator structure emerges automatically once two generic ingredients are present: (i) baryonic elements possess intrinsic oscillatory degrees of freedom— well established in gravitational physics through stellar p and g modes, Jeans oscillations, and vertical and breathing modes of disks (Cox [ 3 ]); and (ii) these elements communicate through a causal retarded propagation kernel that transmits their oscillations to neighbouring elements with finite delay. Whenever spatially distributed oscillators interact through such a causal medium—acoustic, electromagnetic, elastic, plasma, optical, or gravitational—linearity and retardation inevitably generate phase lag, delayed mutual forcing, and therefore a collective damped–oscillator susceptibility. Weak–field GR inherits this same universal structure once the Poisson truncation is lifted. Baryonic elements therefore behave as gravitational scatterers: each responds to incident curvature with a linear susceptibility χ⋆ ( ω ) and re–emits a delayed perturbation. When many such elements participate, their retarded interactions produce a gravitational analogue of multiple scattering. The resulting ensemble renormalises the propagator in the standard Dyson sense, yielding an effective medium with a collective restoring frequency Ω, damping rate Γ, and propagation speed ceff . The coherent curvature response, Φ res , is purely dynamical. Because the retarded kernel carries long memory, its oscillatory part does not cancel cycle–to–cycle. A slowly accumulated, reversible curvature energy builds up, defining an effective mass density ρeff ∝ ⟨| Φ res|2⟩ that sources the ordinary Poisson equation. Thus all local tests of gravity remain unchanged, while extended systems (spirals, ellipticals, and clusters) develop substantial effective dark–mass envelopes from stored curvature energy. In this multiple–scattering framework, coherent baryonic motion in rotating disks, orbiting cores, or clustered environments drives constructive rescattering of curvature. Depending on spatial contrast and coherence length, this produces 5 two distinct observational regimes: large–scale coherent modes in spirals, and diffusive effective–medium halos in ellipticals and clusters. Homogeneous or weakly perturbed environments, lacking spatial contrast, exhibit no rescattering and remain resonantly silent, a fact consistent with the early Universe and the CMB. The present work develops this causal–resonant picture from first principles. Beginning with the retarded solution of linearised GR and the linear susceptibility of baryonic matter, we derive a self–consistent effective medium that accounts for spiral galaxies, ellipticals, clusters, and cosmological structure without introducing new fields or non–baryonic particles. The resulting Spacetime Resonance Theorem provides the microscopic foundation for the collective resonant mode and the stored curvature energy that behaves observationally as dark mass. 2 Spacetime Resonance Theorem The central result of this work rests upon a structural feature of all linear oscillator ensembles embedded in any causal medium governed by a retarded Green function. This behaviour is entirely general: it is found in acoustics, electromagnetism, plasma physics, condensed–matter systems, seismology, elastic media, and even in optical systems such as mode–locked lasers, where many cavity modes synchronise through retarded electromagnetic coupling. To make this universality explicit, we begin by formulating a general theorem— independent of gravitational assumptions—which codifies the emergence of collective frequencies, phase locking, renormalised propagation, and the universal damped–oscillator susceptibility arising from causal retardation and spatially distributed oscillatory elements. Weak–field general relativity will then appear as a direct and natural specialisation of this general structure. 2.1 The Fundamental Oscillator–Medium Theorem Theorem (The Fundamental Oscillator–Medium Theorem). Let {qi ( t ) } be an ensemble of linear oscillators (mechanical, electromagnetic, elastic, gravitational, or otherwise) whose mutual interaction is mediated by a medium described by a causal retarded Green function Gret ( x, x′ ). Assume only linearity, causality, and that the oscillatory elements are spatially distributed, so that their sourcing and rescattering of the field occur at distinct locations. Then the following universal properties hold: 1. Collective frequency and phase emergence. Retarded coupling induces mutual phase pulling among the oscillators. If the coupling bandwidth exceeds the natural detuning spread, the system organises into a phase– locked collective mode with a single emergent frequency Ωthat differs from the independent natural frequencies ωi. In strong coupling, Ω≪ωi,(1) 6 reflecting the dominance of long–wavelength collective motion. 2. Universal damped–oscillator susceptibility. The linear response of each oscillatory element to a harmonic component of the field mediated by Gret is universally described by the susceptibility χ(ω) = 1 ω2 0−ω2−iγω,(2) where ω is the driving frequency, ω0 is the natural (undamped) resonance frequency of the isolated oscillator, and γ is its local damping coefficient. This Lorentzian form is the unique analytic structure compatible with linearity and causality and appears identically in mechanical, electromagnetic, plasma, elastic, and gravitational response theories. 3. Multiple scattering and self–energy renormalisation. Because the oscillatory elements are spatially distributed, the field they source is regenerated at distinct locations, producing multiple scattering. The collective effect of these rescattering processes modifies the propagator through the Dyson relation G−1 eff (ω, k) = G−1 0(ω, k)−Σ(ω, k),(3) where: •Geff(ω, k)is the effective (ensemble) propagator, •G0 ( ω, k)is the propagator of the underlying medium in the absence of scattering elements, • Σ( ω, k)is the self–energy, summarising all multiple–scattering and feedback processes, •ωis the temporal frequency, and kis the spatial wavenumber. In the long–wavelength, low–frequency regime the self–energy admits the universal local expansion Σ(ω, k)≃a0+a1(iω)+a2ω2+b2k2,(4) where a0 , a1 , a2 , and b2 are real coefficients determined by the statistics, strength, and distribution of the oscillatory elements. This expansion yields a homogenised second–order evolution equation with collective parameters (Ω , Γ , ceff ), where Ωis the emergent collective frequency, Γthe effective damping rate, and ceff the emergent propagation speed. These statements hold for all oscillator ensembles in all causal media with retarded propagators. They constitute a general structural principle of linear response theory, independent of any particular physical realisation. 7 Theorem 2.1 expresses a pattern that is so widespread across physics that it is rarely stated as a single principle: whenever spatially distributed oscillatory elements communicate through a retarded field, the resulting system acquires collective modes, a renormalised propagator, an effective medium description, and a universal damped–oscillator susceptibility. The emergence of the self–energy Σ is nothing more than the compact representation of multiple scattering between the spatially separated sources. No reference to gravity is required for these conclusions. 2.2 The Spacetime Resonance Theorem Weak–field general relativity fits precisely within the framework of Theorem 2.1. When the Newtonian limit is not imposed and the full retarded gravitational Green function is retained, spacetime curvature plays the role of the causal medium; baryonic mass elements provide the spatially separated oscillatory sources; and gravitational rescattering generates the corresponding self–energy. The gravitational case therefore inherits, without additional assumptions, the entire universal structure: collective frequency, susceptibility, effective propagation, and resonance. Theorem (The Spacetime Resonance Theorem). Let Φ( x, t ) denote the weak–field gravitational potential defined as the linear perturbation of the flat background metric, Φ(x, t) = −4πGZGret(x, x′)ρb(x′, t′)d4x′,(5) where Gret is the full retarded Green function of linearised general relativity. Decompose the baryonic density into a static part and a coherent, time–dependent perturbation, ρb(x, t) = ρb,0(x)+ρb,1(x, t),(6) where ρb,0 sources the ordinary Newtonian potential and ρb,1 drives dynamical curvature. The static background field ΦN(x)generated by ρb,0satisfies ∇2ΦN(x) = 4πG ρb,0(x),(7) while ρb,1 sources only the dynamical channel. By linearity of the weak–field Einstein equation, the total gravitational potential separates exactly as Φ(x, t) = ΦN(x)+Φres(x, t).(8) Because the baryonic distribution is spatially extended, curvature sourced by ρb,1 undergoes gravitational rescattering. The resulting multiple–scattering processes modify the propagator through the Dyson relation G−1 eff (ω, k) = G−1 0(ω, k)−Σ(ω, k),(9) 8 where G0 ( ω, k)is the retarded propagator of the linearised Einstein operator in flat space prior to rescattering. Its inverse has the standard form G−1 0(ω, k) = −ω2+c2k2,(10) up to gauge conventions, and Σ( ω, k)is the self–energy encoding all feedback effects. In the homogenised long–wavelength, low–frequency regime the self–energy adopts the universal expansion Σ(ω, k)≃a0+a1(iω)+a2ω2+b2k2,(11) with real coefficients a0, a1, a2, b2 determined by the baryonic susceptibility and spatial statistics. Then the coherent curvature response Φ res sourced by ρb,1 satisfies the effective– medium evolution equation ∂2 tΦres + 2Γ ∂tΦres + Ω2Φres −c2 eff∇2Φres = 4πG ρb,1(x, t),(12) with identifications Ω2=a0,Γ = a1, c2 eff =c2(1 −b2).(13) Here Ωis the collective resonant frequency of the gravitational medium, Γis the effective damping rate, and ceff is the effective propagation speed for coherent curvature fluctuations. Moreover, in the absence of coherent driving, ρb,1(x, t)→0 =⇒Φres(x, t)→0,(14) so the resonant curvature mode is purely dynamical and produces no static modification to Newton–Poisson gravity. In the static limit Ω → 0, Eq. (12) reduces to Eq. (7), yielding the ordinary Poisson potential sourced by ρb,0alone. Theorem 2.2 is not an additional physical assumption; it is the gravitational realisation of the universal behaviour stated in Theorem 2.1. The full retarded GR propagator supplies the causal medium; baryons act as spatially distributed oscillatory elements; and gravitational rescattering produces the effective self– energy. The resulting resonant curvature field is therefore a natural, inevitable collective mode of weak–field GR in an extended baryonic environment. 2.3 Corollary: Density–Controlled Restoring Frequency (Motivated by Known Gravitational Modes) The Spacetime Resonance Theorem identifies the collective restoring term of the resonant curvature mode with the self–energy coefficient Ω 2 = a0 . We now 9 For clarity, and to distinguish it from the familiar metric gravitational wave of vacuum GR, we refer to this collective excitation as a masson. The term is descriptive rather than ontological: a masson is not a new particle species, but the coherent curvature mode of the effective gravitational propagator generated by the baryonic ensemble itself. Each masson carries phase and reversible energy, propagates with effective speed ceff , and couples to baryonic structure through the coherent susceptibility χ⋆(Ω) introduced in the scattering theory. The reversible energy associated with these modes is the quadratic curvature energy density already defined in Eq. (18) . This energy behaves gravitationally as an effective mass density and later feeds Channel 2 through the dark–mass mechanism developed in Section 6. In the driven system, the coefficient Γ controls linewidth and temporal phase coherence, with the corresponding quality factor Q=Ω 2Γ.(30) High– Q media (Γ ≪ Ω) support long–lived, spatially coherent curvature oscillations; macroscopic signatures of masson dynamics sustained by multiple gravitational scattering. The masson therefore provides the microscopic language for Channel 1: a collective curvature excitation of the effective gravitational medium. Its propagation, linewidth, and energy transport form the basis of the ensemble scattering and coherence theory developed in Section 5. Before constructing that scattering theory, we require the corresponding spatial propagation law for these modes. This follows from linear perturbation of Eq. (24) , expanding the resonant field about a quasi–static background and resolving it into spatial Fourier components. The resulting dispersion relation, group velocity, damping rate, and energy balance are summarised in the next subsection; these constitute the linear inputs for the multiple–scattering framework that follows. 4.4 Linear Spatial Perturbations and the Propagation Law In this section we summarise the spatial perturbation theory of the causal–resonant field. Nothing new is introduced here: this is the standard linearisation used in wave mechanics to extract dispersion, damping, and energy transport. Its purpose is simply to determine how the already–identified collective mode of Channel 1 propagates spatially as a weak gravitational wave within the effective medium established by multiple scattering. For completeness, the few algebraic steps that lead to the relations summarised below follow directly from inserting a small perturbation δ Φ into the effective-medium PDE of Channel 1 and performing the standard Fourier decomposition. 16 Linearisation about a quasi–static background yields the scalar wave–type equation ∇2δΦ + 1 c2 eff −∂2 t−2Γ ∂t+ Ω2δΦ = 4πG c2 eff δρb,(31) where δ Φ is the weak spatial perturbation of the resonant field and δρb its small baryonic source. This perturbative split is purely mathematical and should not be confused with the physical temporal split ρb = ρb,0 + ρb,1 used to define the two gravitational channels; here we examine only the spatial wave content of the established collective mode. A Fourier mode δΦ∝ei(k·x−ωt)satisfies the dispersion relation ω2+ 2iΓω−Ω2=−c2 effk2,(32) from which ω2 r= Ω2+c2 effk2,|δΦ|∝e−Γt,(33) and the phase/group velocities follow as vph =ωr k, vg=c2 effk ωr .(34) Thus Channel 1 supports weakly damped spatial waves with natural frequency Ω and linewidth Γ, exactly as in reversible acoustic or optical media—consistent with the effective–medium analogy established in [4, 5]. The corresponding energy balance is ∂tEΦ+∇·SΦ=1 c2 eff δρb∂tδΦ−2Γ(∂tδΦ)2,(35) where the Γ–term represents the reversible phase–lag attenuation associated with the frequency–dependent self–energy of the effective propagator. Cycle–averaging yields the usual quality factor, Q=Ω 2Γ,(36) which measures temporal phase coherence. Under harmonic driving at frequency ωthe RMS response is |δΦ|RMS =4πG |δρb| p(Ω2−ω2)2+ (2Γω)2,(37) peaking strongly at ω≃Ω with bandwidth 2Γ. For high driving frequencies ω≫Ω, |δΦ|RMS ∼4πG |δρb| ω2→0,(ω≫Ω),(38) 17 demonstrating that the collective mode acts as a strict gravitational low–pass filter. Only frequencies near Ω can excite Channel 1; higher–frequency baryonic fluctuations are strongly suppressed. This mirrors the behaviour of ideal compressible media, where frequency response is controlled by the effective bulk modulus: high–frequency density perturbations cannot drive the medium. Together, Eqs. (31) – (37) show that the resonant field behaves as a weakly damped, energy–storing continuum supporting dispersive waves with finite bandwidth ∆ω= Ω/Q. The associated coherence time, τcoh ∼Q Ω,(39) links the microscopic phase–lag parameter Γ to the macroscopic persistence of curvature oscillations. As demonstrated in Section 6, the effective memory of the system extends to τmem ∼ Genvτcoh , providing the basis for sustained coherent amplification in Channel 1. These linear propagation and damping laws supply the foundation for the single– and multiple–scattering theory developed next in Section 5. 5 Gravitational Scattering Theory Having established that the dynamical (Channel 1) sector of gravity is governed by the effective resonant wave equation derived in the previous section, we now construct the corresponding gravitational scattering theory. Throughout this section we work exclusively with the oscillatory baryonic component ρb,1 , since only temporally varying sources can excite and sustain the collective resonant field Φ res . The static component ρb,0 generates the background Newtonian field but does not contribute to scattering. The goal in what follows is to develop the gravitational analogues of the Lippmann– Schwinger, Foldy–Lax, and Dyson equations using the propagator Gret associated with the effective-medium PDE. This provides the microscopic foundation for the collective behaviour analysed later—including coherent spiral modes and diffusive elliptical/cluster regimes (Section 8). In this framework each baryonic mass acts as a gravitational scatterer characterised by a linear susceptibility χ⋆ (Ω), which encodes its phase-lagged response to an oscillatory potential at frequency Ω. More generally, the gravitational response of a localized mass is a polarizability tensor χij ⋆ (Ω), whose symmetry is set by the local spacetime environment. In isotropic settings the tensor reduces to the scalar χ⋆ (Ω) used throughout this section, whereas in anisotropic geometries the effective medium develops preferred directions of response. This tensorial structure will become essential when comparing planar spiral systems (with lateral polarization of coherent gain) to the isotropic coherence cells of ellipticals, but for the present microscopic theory the scalar form suffices. 18 An incident perturbation δ Φ inc produces a delayed re-emission described by the retarded Green function Gret , and the resulting scattered field propagates to all other masses. Repeated re-scattering among many inclusions generates a hierarchy of coupled responses, precisely analogous to classical multiple-scattering theory in acoustics and electromagnetism [6, 7, 8]. Thus baryonic matter plays a dual role: its oscillatory component ρb,1 acts as the continuous driving source of the resonant field, while the same masses also act as discrete gravitational scatterers that re–emit curvature perturbations through their susceptibility χ⋆ (Ω). The iterative network of scattering and re-scattering links the microscopic delayed response of individual baryons to the emergent, large-scale propagation properties of the collective resonant field. It therefore supplies the necessary bridge between the effective wave operator introduced earlier and the macroscopic gravitational behaviour developed in the remainder of the book. 5.1 Single–Scatterer Response A discrete baryonic mass m⋆ with linear gravitational susceptibility χ⋆ (Ω) responds to an incident oscillatory potential Ainc through the standard Lippmann–Schwinger relation Asc(x)=t⋆(Ω) GΩ(x−x⋆)Ainc(x⋆), t⋆(Ω) = 4πG m⋆χ⋆(Ω),(40) where GΩ is the frequency–domain form of the retarded Green function associated with the effective–medium propagation equation. The quantity t⋆ (Ω) defines the gravitational scattering strength of the inclusion, capturing its phase–lagged re–emission of the incident curvature perturbation. In general, the susceptibility of a localized gravitating mass is a polarizability tensor χij ⋆ (Ω). Local geometric anisotropies (for example, planar kinematics or coherent streaming) can preferentially amplify particular tensor components, thereby polarizing the local gravitational response. In the isotropic approximation used here the tensor reduces to its trace, χ⋆ (Ω), but the full tensor structure will reappear in Section 8.1.8 (spiral galaxies) and Section 8.2.10 (elliptical galaxies). Because the resonant sector is weakly damped (Γ ≪ Ω) and the underlying dynamics are those of linearised GR, the interaction is reversible: no curvature energy is irreversibly absorbed by the scatterer. Instead, the inclusion stores curvature perturbations only transiently, re–emitting them with a phase delay determined by χ⋆(Ω) and the causal structure of GΩ. Equation (40) therefore provides the microscopic building block of the gravitational multiple–scattering theory developed below, linking the individual baryonic response to the collective, ensemble–level behaviour of the resonant field. 19 5.2 Multiple Scattering and the Foldy–Lax Hierarchy In classical wave physics, the interaction of many scatterers with a coherent incident field is described by the Foldy–Lax hierarchy [ 6 , 7 , 8 ], which provides the standard self–consistent formulation of multiple scattering among discrete inclusions. Originally developed for acoustics and later generalised to electromagnetism, elasticity, and photonics, the same framework applies directly to the present gravitational problem once the effective propagator GΩ has been specified by the linearised dynamics of Section 4.4. Throughout this subsection the field amplitude A denotes the oscillatory (resonant) curvature perturbation Φres driven by the baryonic component ρb,1. Each baryonic inclusion re–emits the incident curvature fluctuation with a phase delay determined by its susceptibility χ⋆ (Ω). These delayed re–emissions propagate to all other inclusions via the Green function GΩ , and the resulting network of mutual couplings constitutes a gravitational analogue of standard multiple– scattering theory. When the spacing between inclusions is much smaller than the resonant wavelength ( kΩd≪ 1), their phases remain locked to the incident field and the scattered amplitudes add coherently. In this fully coherent Foldy–Lax limit, Acoh ∝N|χ⋆(Ω)|Ainc,(41) so that the ensemble behaves as a single collective radiator with intensity Icoh ∝N2|χ⋆(Ω)|2|Ainc|2,(42) It is convenient to quantify this as an environmental gain Genv ≡|Atot|2 |Ainc|2,Genv =N2|χ⋆(Ω)|2(coherent limit).(43) For an ensemble of N scatterers, the total field satisfies the discrete Lippmann– Schwinger system A(x)=Ainc(x)+t(Ω) X n GΩ(x−xn)A(xn),(44) with t(Ω) = 4πGm⋆χ⋆(Ω). In operator notation this becomes A= (1 −tGΩ)−1Ainc,(45) whose Neumann expansion, A=Ainc +tGΩAinc + (tGΩ)2Ainc +··· ,(46) sums all orders of mutual gravitational re–scattering. In the limit Γ → 0 (purely real GΩ ), the response is strictly coherent and Eq. (45) reduces to the linear scaling of Eq. (41). 20 If instead the relative phases of different inclusions become effectively random—because of geometric delays, finite memory (Γ > 0), or dynamical phase drift—then scattered amplitudes no longer combine coherently. They add in intensity only: Itot =N|χ⋆(Ω)|2|Ainc|2,(47) yielding the standard incoherent scaling Genv =N|χ⋆(Ω)|2,pGenv =√N|χ⋆(Ω)|.(48) The spiral–disk regime discussed in Section 8.1 lies in the long–wavelength, phase– ordered limit ( kΩd≪ 1), where both spatial phasing and delayed temporal response organise coherently into a global m = 1 pattern. As coherence decreases—due to finite Γ, increasing disorder, or large separations—the system transitions smoothly toward the statistical regime treated in the next subsection, where ensemble averaging of Eq. (44) produces an effective–medium description appropriate to elliptical and cluster environments. This subsection has treated multiple scattering for a fixed arrangement of scatterers, with coherence determined purely by the phase relations among their delayed re–emissions. The following subsection extends the framework to statistically disordered or dynamically evolving ensembles, where averaging over many such configurations yields the emergent effective parameters of the resonant gravitational field. 5.3 Effective–Medium and Ensemble Averaging The Foldy–Lax hierarchy describes the fully coherent interaction of sub-wavelength baryonic scatterers with an incident oscillatory gravitational field. In realistic configurations—particularly in irregular, thick, or dynamically evolving systems—phase coherence among scatterers is degraded as separations approach the resonant wavelength or as differential motions introduce random phase delays. Under such conditions the local fields A (x n ) at different inclusions can no longer be treated as perfectly correlated, and the collective behaviour must be described statistically through an ensemble average over possible configurations. Averaging the discrete Lippmann–Schwinger equation (44) over many realisations of scatterer positions yields the standard effective–medium (Dyson) equation for the coherent mean field: (∇2+k2 eff)⟨A⟩= 0, k2 eff =k2+ Σ(Ω),(49) where Σ(Ω) is the ensemble self–energy. The self–energy encapsulates the average phase delay and finite dwell time associated with multiple gravitational rescatterings among baryonic inclusions. Its real part modifies the dispersion relation, while its imaginary part determines the ensemble attenuation rate: Γeff =ceff Im keff,(50) 21 which is reduced relative to the microscopic damping Γ whenever reversible rescattering retains phase coherence over extended domains. The degree of coherent amplification produced by the environment can be summarised using the environmental gain factor Genv , defined earlier in Eq. (48) as the total intensity gain. In the effective–medium limit this quantity can be conveniently written in terms of the coherence volume Vcell of the emergent Bloch-like mode: Genv ∝1 Vcell =1 ℓ2 ⊥ℓ∥ ,(51) where ( ℓ⊥, ℓ∥ ) denote the transverse and longitudinal coherence lengths of the averaged field. From a microscopic perspective, the local intensity gain in a coherent patch scales as the square of its average coupling efficiency divided by its coherence volume: Genv ∝|X|2 Vcell =|X|2 ℓ2 ⊥ℓ∥ ,(52) where X is the effective (ensemble–averaged) coupling amplitude. Large coherence domains (large ℓ⊥ and ℓ∥ ) thus imply strong environmental gain, while shrinking coherence lengths reduce Genv by increasing the number of independent scattering cells. Systems with strong spatial phase order (e.g. thin, rotating disks) possess large coherence volumes and hence large environmental gain. These satisfy Γeff ≃Γ,(53) recovering the coherent Foldy–Lax limit. Increasing disorder or dynamical randomness decreases the coherence lengths, lowers Genv , and drives the system toward regimes where ensemble damping is dominated by phase randomisation rather than the microscopic Γ. Equation (49) therefore provides the statistical extension of the coherent multiple–scattering framework. Ordered systems correspond to real, phase–aligned self–energies with large environmental gain, whereas disordered ensembles are described by complex, stochastic self–energies whose imaginary parts encode diffusive loss of coherence. This effective–medium formulation supplies the macroscopic link between the microscopic scattering interactions and the global intensity structure analysed later in Section 8.2. 5.4 Transport and Diffusion Limits When the relative phases of individual baryonic scatterers become random and time–dependent, the ensemble field correlation function obeys the standard radiative–transfer equation familiar from acoustic, electromagnetic, and elastic multiple–scattering theory [8]: ˆ s·∇I(x,ˆ s)=−I(x,ˆ s) ℓs +1 ℓsZP(ˆ s,ˆ s′)I(x,ˆ s′)dΩ′+S(x,ˆ s),(54) 22 where I (x ,ˆ s ) is the directional intensity of the oscillatory gravitational field, ℓs is the scattering mean free path, and P is the angular redistribution kernel specifying how gravitational perturbation energy is redirected by baryonic inclusions. Equation (54) is the exact gravitational analogue of the radiative–transfer formulation validated in resonant acoustic experiments by Derode, Tourin, and Fink [9, 4, 5]. In the isotropic limit, Eq. (54) reduces to the usual diffusion equation for the ensemble-averaged intensity: ∇2I(x) = I(x) L2 damp,eff ,(55) where the effective damping (or coherence–decay) length is Ldamp,eff =ceff Γeff .(56) This length characterises the scale over which the coherent component of the oscillatory gravitational field survives repeated random-phase rescattering. In ordered or phase-aligned systems, Ldamp,eff →∞ , recovering the coherent Foldy–Lax regime, whereas in disordered or dynamically fluctuating ensembles it becomes finite, reflecting the diffusive loss of coherence arising solely from phase randomisation in weak-field GR. Thus Eqs. (54) – (56) provide the statistical, transport-level closure of the multiple–scattering framework: they describe how gravitational curvature perturbations propagate once coherent additions among scatterers are no longer maintained. 5.5 Unified Interpretation The coherent (Foldy–Lax) and diffusive (transport) limits developed above arise from the same governing relation, A=1−t(Ω) GΩ−1Ainc,(57) and differ only in the statistical structure of phase correlations among the baryonic scatterers. In the microscopic description the response of each inclusion is encoded in its gravitational susceptibility, which is in general a polarizability tensor χij ⋆ (Ω). In isotropic environments this reduces to the scalar form χ⋆ (Ω) used throughout Section 5, but in anisotropic geometries the tensor develops preferred directions of response that feed directly into the large-scale coherence properties. When relative phases are organised (as in rotating disks), multiple rescattering produces coherent amplitude addition and large environmental gain. This corresponds to a situation in which the tensorial susceptibility is effectively aligned: a dominant in-plane polarisation sustains lateral coherence over large scales, 23 producing the global m = 1 spiral modes described in Section 8.1. Here the medium behaves as a phase-ordered resonator, and the Foldy–Lax limit applies. When phases are randomized (as in irregular, thick, or pressure-supported systems), the tensorial response averages to its isotropic trace and the same operator (57) generates a diffusive intensity field with a finite coherence length. This is the regime analysed in the elliptical and cluster cases, where ensemble averaging produces an isotropic effective medium governed by a Dyson self–energy Σ(Ω) and a transport damping rate Γeff . Thus the unified scattering equation couples seamlessly to both extremes: • coherent reinforcement in spiral galaxies, where geometric and kinematic anisotropies polarize the susceptibility and stabilise large-scale phase locking; • diffusive, transport-dominated behaviour in ellipticals and clusters, where random three-dimensional motions isotropise the tensor response and reduce the coherence length. In both cases the emergent macroscopic behaviour follows directly from how the microscopic polarizability, mediated by gravitational rescattering, organises (or fails to organise) phase correlations across the baryonic ensemble. 6 Spatiotemporal Feedback and Resonant Coherence 6.1 Constituent Entities & States of the Feedback System Spacetime in the weak–field limit does not behave as a purely passive background but as a causal medium whose curvature degrees of freedom respond on two distinct timescales. Channel 1 captures the dynamical, resonant response of the metric to the oscillatory baryonic component ρb,1 ;Channel 2 captures the quasi–static, Newtonian response to the slowly varying mass distribution and to the curvature energy accumulated from past Channel 1 activity. Channel 1 is therefore the dynamical driver of the feedback system: its resonant curvature oscillations exchange energy and phase with the baryons and deposit curvature energy into the long–term state variable ρeff .Channel 2 provides the reactive, quasi–static response of the metric to the slowly varying mass distribution and the accumulated stored curvature energy. The interaction of the two channels forms a coupled spatiotemporal feedback system: Channel 1 continuously modifies the long–term metric state, while Channel 2 sets the background potential within which resonant oscillations occur. The architecture thus consists of two physical entities—baryonic matter and the spacetime response kernel—and two curvature states that evolve on fast (Channel 1) and slow (Channel 2) timescales, summarised in Table 1. 24 Table 1: Constituent entities and states in the two–channel curvature–matter feedback system. Entity / State Description Baryonic mass density ρb=ρb,0+ρb,1 (Entity) The physical source of curvature. ρb,1 is the oscillatory component that drives Channel 1 and acts as both emitter and scatterer of the resonant metric response. ρb,0 is the slowly varying component sourcing the Newtonian Channel 2 field. Spacetime response parameters B= (ceff ,Ω,Γ) (Entity) The parameters describing the causal, frequency–dependent linear response of Channel 1: ceff governs propagation, Ω the restoring term, and Γ the finite memory time. Channel 2 uses the same metric but responds through the instantaneous Poisson relation. effective curvature field Φ = Φres + ΦN (State of the spacetime) Φ res is the resonant Channel 1 field driven by ρb,1 , mediating phase–lagged curvature exchange with the baryons. Φ N is the quasi–static Channel 2 field sourced by ρb,0 together with the accumulated stored curvature energy ρeff. Stored curvature energy (effective ‘dark’ mass) ρeff (Long–term state variable) The slowly evolving curvature energy built up from past Channel 1 activity. It decays only on the timescale Γ −1 and contributes to the Newtonian potential via ∇2 Φ N = 4 πG ( ρb,0 + ρeff ). 25 6.3.3 Polarisation of Temporal Coherence The temporal coherence of the resonant channel is not only a matter of frequency locking; it also possesses a polarisation structure determined by the gravitational susceptibility of the baryonic ensemble. In the general case the linear response of a mass element to an oscillatory curvature perturbation is described by a susceptibility tensor χij(Ω),(63) which maps an incident curvature perturbation Aj into a delayed re–emitted field proportional to χijAj . The effective gravitational polarisation states are given by the eigenvectors of χij(Ω), obtained by diagonalising χijvj (α)=χ(α)vi (α),(64) so that each eigenmode vi (α) behaves as an independent, polarised oscillator with scalar susceptibility χ(α). In a statistically isotropic medium (e.g. pressure–supported ellipticals) the tensor reduces to a multiple of the identity, χij(Ω) = χ(Ω) δij,(65) so all polarisations respond equally; no preferred temporal coherence direction exists. This isotropy leads naturally to three–dimensional coherence domains and the spherical transport envelopes characteristic of ellipticals. In rotating discs, by contrast, the baryonic mass distribution and its dynamical response are strongly anisotropic. Vertical excursions are weak, thin, and largely incoherent, whereas in–plane motions involve substantial coherent mass with long memory. Accordingly the susceptibility tensor decomposes into two unequal sectors, χij = diagχ∥, χ∥, χ⊥,|χ∥|≫|χ⊥|,(66) so the curvature–matter feedback channel amplifies perturbations in the galactic plane far more efficiently than those perpendicular to it. The eigenmodes of Eq. (64) therefore separate into: •adominant in–plane polarisation with large eigenvalue χ∥, and •aweak vertical polarisation with small eigenvalue χ⊥. The resonant feedback loop selects whichever eigenmode yields the greatest coherent reinforcement. Because the environmental gain satisfies Genv,(α)∝ |χ(α)|2, the in–plane mode overwhelmingly dominates in thin discs, forming the planar ( xy ) breathing pattern that drives global m = 1 spiral modes (Section 8.1). Feedback then locks both the baryonic forcing and the effective curvature field 32 into the same polarisation eigenvector: the oscillatory mass distribution drives curvature primarily in the xy plane, and the effective curvature field in turn acts back on the baryons in that same direction, closing the self–consistent loop. In isotropic systems (ellipticals, clusters) no such preferred eigen–direction exists. The susceptibility tensor has three equal eigenvalues, so the resonant mode is unpolarised in the mean: each polarisation contributes equally to the temporal coherence, producing the spatially isotropic speckle cells and large–scale transport envelopes derived in Section 8.2. In mildly anisotropic systems (triaxial ellipticals, S0 galaxies), the eigenvalues of χij split mildly, and the dominant temporal coherence aligns along whichever principal axis maximises |χ(α)|2, giving rise to systematically flattened but non–rotational transport halos. Thus the polarisation of temporal coherence is not an independent ingredient but a direct consequence of the susceptibility tensor of the ensemble. Geometry and kinematics determine χij (Ω); diagonalisation selects the polarisation eigenmodes; and the feedback loop amplifies whichever mode has the largest eigenvalue. Spiral discs, lenticulars, ellipticals, and triaxial systems thereby emerge as distinct polarisation phases of the same curvature–matter resonator. 7 Spatial Coherence of the Resonant Channel The preceding sections established that the resonant sector (Channel 1) acts as a causal, self–locking oscillator whose feedback dynamics drive the system toward a globally synchronized state. Once temporal coherence has been achieved, all baryonic emitters contributing to the oscillatory component ρb,1 acquire a common phase and oscillate at the same resonant frequency Ω = ωres . The purpose of this section is to describe the corresponding spatial coherence of the Channel 1 effective curvature field. Because Channel 1 is defined entirely by the oscillatory source ρb,1 , all coherent structures discussed here belong strictly to the resonant sector. The spatially coherent field Φ res is not directly observable. Its macroscopic imprint enters only through the slowly accumulated, time–averaged stored component ρeff ∝DGΩ∗ρb,1 2Eτmem ,(67) where the average is taken over the long memory time τmem = Γ −1 . Equation (250) expresses the fact that only the stored curvature energy feeds Channel 2; the coherent spatial patterns within Channel 1 remain internally dynamical and observationally silent. The spatial coherence kernel GΩ (x , x ′ ) implicitly defines a spatial polarizability tensor for the resonant channel: it selects which directions of the baryonic fluctuations couple most efficiently into the coherent field. Disk–like coherence domains amplify lateral components of the response, whereas spherical or ellipsoidal coherence cells weight the components almost isotropically. This geometric dependence 33 sets the relative spatial response weights that later enter the susceptibility tensor χij (Ω); its explicit consequences are developed in the spiral and elliptical sections. As proposed in Section 6.3, curvature-mediated frequency pulling drives neighbouring regions toward a common temporal frequency, so a global carrier Ω is already selected. Spatial coherence therefore concerns the geometry of the resonant envelope at this locked frequency, not the mechanism of temporal locking. With temporal coherence established, the resonant field admits the decomposition Φres(x, t) = ReA(x)e−iΩt,(68) where the complex envelope A (x) satisfies a stationary Helmholtz–type equation. Constructive interference among the phase–locked modes produces long–lived, slowly evolving resonant envelopes, the gravitational analogue of standing–wave patterns in optical or acoustic cavities. These envelopes evolve only through slow baryonic advection and the long curvature–memory time of the resonant sector. The next subsection derives the corresponding Helmholtz envelope equation and its exact Kirchhoff–Sommerfeld representation, identifying the spatial coherence scales that govern the internal structure of the resonant field. 7.1 Kirchhoff–Sommerfeld Representation Assuming global synchronization at frequency Ω, all curvature modes carry the common temporal factor e−iΩt . Substituting the ansatz (68) into the linearized field equation (31) and neglecting slow temporal modulations gives the stationary Helmholtz envelope equation ∇2A(x) + Ω2+i2ΓΩ c2 eff A(x) = 4πG c2 eff δρb(x),(69) which determines the spatial profile of the phase–locked resonant field. The exact solution of Eq. (69) is the standard Kirchhoff–Sommerfeld representation: A(x) = G c2 eff Zd3x′eikΩ|x−x′| |x−x′|δρb(x′),(70) with complex wavenumber kΩ=Ω+iΓ ceff , λres =2πceff Ω.(71) The imaginary part of kΩ produces a slow spatial attenuation over the damping length Ldamp =ceff Γ=ceffτmem,(72) which measures the spatial extent of curvature memory. For weak damping (Γ≪Ω) and for r≪Ldamp, one may take kΩ≃2π/λres as real. 34 On galactic and cluster scales the coherence length of the resonant curvature field is of order the resonant wavelength, ℓcoh ∼λres . In spiral galaxies the baryonic disc occupies only a fraction of this coherent envelope, typically adisk ∼ 0 . 1– 0 . 2 λres , but the entire region still belongs to a single phase–locked resonant domain. Within such a domain the curvature phase varies only slowly, and the envelope A (x) behaves as a quasi–static, phase–coherent standing pattern maintained by curvature memory and delayed rescattering, rather than by geometric boundaries as in optical cavities. This definition of the coherence domain is fully consistent with the polarisation–based coherence cells introduced earlier: its size is determined by the curvature wavelength λres , not by the geometric extent of the baryonic material contained within it. In this representation the source term δρb (x) in Eq. (69) should be understood as the coherent Channel 1 driver: the oscillatory component of the baryonic density ρb,1expressed in the global resonant phase frame. Formally, one may write δρb(x) = δρb(x)eiφd(x),(73) where φd (x) encodes the local phase delay of the baryonic motion relative to the curvature carrier e−iΩt . The complex amplitude in Eq. (69) therefore represents the phase–weighted mass density that actually couples to the resonant mode, not a static geometric overdensity. The structure above is directly familiar from scalar Fourier optics 2 (see Goodman [ 11 ]). In a coherent wave system the physical source is represented by acomplex distribution, and the field generated at any point is obtained by a phase–weighted superposition of all emitter contributions. In this language the oscillatory baryonic source δρb (x) = |δρb|eiφd(x) plays the role of a complex aperture distribution, while the curvature kernel GΩ supplies the propagating phase front. The field envelope A (x) is therefore obtained by the standard coherent integral defined in Eq. (70) , which is exactly the Kirchhoff–Sommerfeld representation of the Helmholtz solution. The complex phase of δρb ensures that each baryonic emitter is automatically “rephased” by the propagation kernel before contributing to the field; the integral is simply the coherent superposition of these phase–weighted contributions. This form makes clear that the effective drive of the resonant channel is not determined by the geometric asymmetry of the mass distribution but by its coherent projection onto the curvature phase front. In general, the Green function GΩ(x,x′) = eikΩ|x−x′| |x−x′|(74) 2 Unlike standard free–space optical propagation—which lacks an internal distributed memory—this curvature–mediated regime behaves as a quasi–static, phase–coherent resonator whose geometry and long memory determine the envelope structure. Section 7.2 shows how these coherent envelopes generate stored curvature energy in the Poisson channel. Closer physical analogues exist in acoustical media with delayed multiple scattering; see Section 4.2. 35 acts as the continuous phase reference against which the oscillatory mass distribution is compared. The integral A(x) = G c2 eff Zd3x′GΩ(x,x′)δρb(x′) (75) can thus be interpreted as the phase–matched projection of the complex source δρb onto the resonant effective curvature field. In regions where the baryonic phases φd (x ′ ) are aligned with the kernel, contributions add constructively; where they decorrelate, they cancel statistically. A nearly symmetric configuration may therefore drive the resonant field strongly if its internal phase pattern is coherent, whereas a highly asymmetric distribution may drive it weakly if its phases are random. This viewpoint applies equally to ordered and disordered systems. In a spiral disc, the phase pattern φd (x) may align with a low–order azimuthal structure (treated later via explicit Fourier moments), while in an elliptical galaxy the same formalism describes a random ensemble of local “speckle” phasors whose net drive is governed by their phase–weighted average under Eq. (75) . In all cases it is the coherently projected mass, rather than the raw geometric density, that controls the strength of the Channel 1 field and hence the rate at which stored curvature energy ρeff builds up in Section 7.2. 7.2 Effective Potential and Energy Storage in the Resonant Medium Within the two–channel framework developed above, the physically observable gravitational field is the Channel 2 potential Φ N . This potential is sourced by the time–averaged baryonic density ρb,0 together with the long–term stored curvature energy ρeff, ∇2ΦN= 4πGρb,0+ρeff,(76) and represents the quasi–static geodesic structure of the system. Crucially, Φ N does not contain the instantaneous oscillatory contribution from the resonant channel. The oscillatory effective curvature field Φ res ( t ), driven solely by the coherent baryonic component ρb,1 , averages to zero over a period and is not directly observable. Only its stored curvature energy, accumulated over the memory time τmem = Γ−1, feeds into Channel 2. Because the effective resonant medium possesses finite causal memory, its response to the oscillatory field is phase–lagged. This delay prevents perfect cancellation of curvature stresses over each cycle: part of the resonant curvature energy injected by ρb,1 remains stored in the medium rather than being returned immediately. The result is a slowly evolving curvature bias—the quantity ρeff —that persists on timescales long compared to 2 π/ Ω and sources the observable Channel 2 potential. 36 The local energy density of the Channel 1 resonant field, expressed through the complex spatial envelope A(x) introduced in Eq. (68), is EΦ(x) = 1 8πGhc2 eff|∇A|2+ Ω2|A|2i,(77) representing the reversible partition between curvature–strain gradients and temporal restoring energy. By the equivalence of mass and energy, this stored curvature energy behaves gravitationally as an effective mass density, ρeff(x) = EΦ(x) c2=1 8πGc2hc2 eff|∇A|2+ Ω2|A|2i.(78) This identification ensures that all curvature energy retained in the effective resonant medium enters Channel 2 in precisely the same manner as conventional matter. For statistically incoherent systems—such as elliptical galaxies or cluster environments composed of many locally uncorrelated oscillators—the relevant observable is the ensemble average of Eq. (78), ρeff(x)=1 8πGc2hc2 eff|∇A|2+ Ω2|A|2i,(79) since both intensity and gradient fluctuations contribute to the macroscopic curvature bias. In either the coherent or statistical limit, the stored curvature energy produces a gravitational field identical in form to that attributed in the ΛCDM framework to a dark–matter component, ∇2ΦN= 4πGρb,0+ρeff⇐⇒ ∇2ΦΛCDM = 4πGρb,0+ρDM.(80) Thus ρeff acts observationally as the “missing mass” required for flat rotation curves and weak–lensing signatures, even though its origin lies entirely in the resonant, energy–storing dynamics of Channel 1 rather than in particulate dark matter. The next section examines how this effective resonant medium responds to a single baryonic source, whose localized oscillation constitutes the fundamental building block of all collective resonant structures. 7.3 Single Point Source in a Resonant Channel–1 Background As in the preceding subsections, all results here apply exclusively to Channel 1. The fields derived are purely resonant, oscillatory curvature responses driven by the coherent baryonic component ρb,1 and are therefore not directly observable. Their only physical imprint arises through the stored curvature energy they deposit into the effective resonant medium, which increments the long–term effective density ρeff that later sources Channel 2. 37 For a resonant point mass mb the Channel–1 driver is represented as a Dirac–delta oscillatory source [11]: δρb,1(x)=mbδ(3) x−xp.(81) Inserting this source into the Kirchhoff–Sommerfeld representation (Eq. (70) ) yields the corresponding Channel–1 resonant field δΦp(r): δΦp(x) = Gmb c2 eff eikΩr r, r =|x−xp|.(82) Equation (82) is simply the Green–function evaluation of the Channel–1 Helmholtz operator with a delta–function source. The real, physical Channel–1 oscillation corresponds to the real part of Eq. (82) , giving δΦp(r) = Gmb c2 eff e−r/Ldamp j02πr λres ,(83) where j0 describes spherical spreading and the exponential factor represents attenuation over the curvature–memory length Ldamp . This oscillation remains confined to Channel 1 and has no direct dynamical effect on matter. The quantity of physical interest is the stored curvature energy associated with this oscillation, obtained from the RMS magnitude |δΦp(r)|21/2=Gmb c2 eff e−r/Ldamp r√2j0(2πr/λres).(84) This RMS field is what enters the resonant energy density and ultimately modifies the effective mass density ρeff . In the general case the effective density contains both gradient and curvature contributions, as in Eq. (79) . For the single, spherically symmetric point–source mode considered here we adopt the long–wavelength, slowly varying–envelope approximation ( |∇A|≪kΩ|A| ), so that the gradient term contributes only a small correction to the total stored energy. Retaining only the dominant curvature term ∝Ω2⟨|A|2⟩yields the illustrative expression below. The contribution of this single oscillatory source to the stored curvature energy is therefore ρeff(r) = Ω2 4πGc2|δΦp(r)|2 =G m2 bΩ2 8π c2c4 eff e−2r/Ldamp r2j2 02πr λres . (85) This stored energy is the only mechanism by which a point source modifies observable dynamics: its contribution adds to the total ρeff appearing in the 38 Channel–2 Poisson equation. The instantaneous oscillatory field never contributes directly to ΦN, but the associated ρeff does. Thus the spherical Channel–1 response to a point mass provides the elementary “building block” from which collective resonant fields and the resulting effective density ρeff of extended systems are constructed. In realistic galactic settings the true ρeff arises from the interference of many such point–source contributions to form the global amplitude envelope A (x). The next subsection therefore considers the interference of two oscillatory point sources as the simplest nontrivial example. 7.4 Interference of Two Point Sources We now examine the spatial analogue of temporal coherence. Consider two point masses, m1 and m2 , located at fixed positions x 1 and x 2 within the same resonant domain. Both oscillate coherently at the common frequency Ω, with a fixed relative phase. At an arbitrary observation point xthe distances to the sources are ri=|x−xi|, i = 1,2. Because both sources radiate simultaneously, their Channel–1 oscillatory effective curvature fields overlap throughout the effective resonant medium, producing spatial interference patterns analogous to optical diffraction from two coherent emitters. In the single–source case, interference is absent: the Newtonian background carries no phase and responds only to time–averaged stored energy. With multiple phase–locked emitters, however, the oscillatory components interfere coherently with one another (but never with the quasi–static background). This mutual interference generates spatial modulation of the resonant envelope, and hence of the stored curvature energy. Because the Channel–1 equation is linear in the complex amplitude, the total field is the superposition of the two Helmholtz solutions: δΦ(x) = δΦp1(x)+δΦp2(x),(86) with each contribution given by δΦpi(x) = Gmi c2 eff eikΩri ri , ri=|x−xi|, i = 1,2.(87) Interference arises entirely from the complex phase factors eikΩri; the physically relevant stored energy depends on the RMS intensity |δΦ|2after superposition. Figures 2 and 3 illustrate the curvature interference pattern generated by two equal, coherently oscillating point masses in the same resonant medium (using the same parameters as in the single–source case: ceff = 2 × 10 5 m s −1 , Q = 5, Ω = 7 . 45 × 10 −16 s −1 , and m1,2 = M⊙ ). At a separation of 10 pc (Fig. 2) the spacing is much smaller than the resonant wavelength λres , so the individual fields merge 39 into a single smooth envelope with only slight modulation. At 100 pc 3 (Fig. 3), several wavelengths fit between the sources and distinct interference fringes appear in both amplitude and ρeff . These alternating bright and dark regions correspond to constructive and destructive curvature interference—essentially spacetime diffraction patterns within the effective resonant medium. The observable curvature intensity is the time–averaged magnitude of the total field: |A|2=|δΦp1+δΦp2|2.(88) Substituting Eq. (87) gives the explicit interference pattern |A|2=G2 c4 eff e−2(r1+r2)/Ldamp "m2 1 r2 1 j2 02πr1 λres +m2 2 r2 2 j2 02πr2 λres  +2m1m2 r1r2 j0 2πr1 λres j0 2πr2 λres cos 2π∆r λres #, (89) where ∆ r = r2−r1 is the path–length difference. The cosine term describes mutual interference; the exponential term conveys the damping imposed by the finite curvature–memory length Ldamp. Constructive interference occurs for 2 π ∆ r/λres = 2 nπ , and destructive interference for (2 n + 1) π . In the long–wavelength limit ( ri≪λres ), the field reduces smoothly to the Newtonian form G ( m1 + m2 ) /r . At intermediate scales, however, coherent interference redistributes stored curvature energy, forming stable resonant envelopes that are the building blocks of extended halo–like structures. In general the effective density contains both gradient and curvature contributions, as in Eq. (79) . For the present two–source configuration we again adopt the long–wavelength, slow–envelope approximation so that the gradient term provides only a small correction. Keeping only the dominant curvature contribution yields ρeff(x) = Ω2 4πGc2|A(x)|2,(90) showing directly that interference between oscillatory fields modulates the local effective gravitational density. The cross–term in Eq. (89) thus represents a real redistribution of stored curvature energy within the effective resonant medium; the total energy remains conserved in the absence of damping, but is shifted between regions of constructive and destructive curvature amplitude. The exponential envelope is again controlled by the same damping rate Γ that sets the temporal coherence time. This two–source configuration makes the phase structure of the Kirchhoff– Sommerfeld representation (Section 7.1) fully explicit. The Channel–1 driver takes the form δρb(x) = δρb(x)eiφd(x),(91) 3 Subsequent simulations (see Section 8) indicate that these separations and wavelengths are not realistic for galaxy–scale resonances, but the configuration remains pedagogically useful. 40 Figure 2: Two–source interference (10 pc separation). Amplitude, phase, and effective density ρeff for two coherent point sources. The separation is much smaller than the resonant wavelength, so the two fields merge into a single envelope with gentle modulation. 41 of the surrounding medium be described by the driven, damped Helmholtz envelope equation ∇2+k2 ΩA=4πG c2 eff ρb,(110) where A is the complex Channel–1 curvature–potential amplitude, kΩ = (Ω + iΓ)/ceff with Γ>0, and ceff is the curvature–wave speed. Suppose the baryonic emitters drive the medium at a single temporal frequency Ω. Then Eq. (110) admits nontrivial steady–state solutions whose azimuthal dependence consists of a single Fourier harmonic, A(r, θ, z, t)=ℜne−iΩteiθ ˜ A(r, z)o,(111) for some bounded spatial profile ˜ A ( r, z ). These solutions represent coherent, rotating m= 1 curvature patterns. A self–consistent mode arises when the driving frequency matches the orbital frequency at some radius r0, Ω = ω(r0) = v0 r0 ,(112) so that baryons near r0 experience minimal phase drift relative to the pattern. For flat rotation curves a unique such radius exists for any choice of Ω ∈ (0 , ωmax ), and the weak shear around r0 ensures that a broad annulus contributes coherently to the m= 1 harmonic. Because Γ > 0, the homogeneous part of Eq. (110) is exponentially damped. Thus the driven configuration (111) is a linearly attracting steady state: perturbations decay and the system converges toward the coherent m = 1 pattern. The resulting configuration defines a self–consistent steady state in which the effective dark–mass distribution and the baryonic rotation field remain mutually compatible. Corollary (Resonant Coherence Condition). The corotation relation (112) should be interpreted as a steady–state attractor: for the frequency Ω satisfying Ω = ω ( r0 ), the effective curvature field at each point remains phase coherent with its past, and the driven solution is stable under damping. Consequently the m= 1 Fourier sector dominates the long–time response of the disk. Remark on Vertical structure. For numerical work it is convenient to approximate the vertical dependence of ˜ A ( r, z ) by the lowest vertical eigenfunction Z0 ( z ) of the finite–thickness disk and to assume approximate separability, ˜ A ( r, z ) ≈M1 ( r ) Z0 ( z ). This yields an effective 2D in–plane problem for the radial envelope M1 ( r ); the approximation is useful but not required for the existence of the global, phase–locked m= 1 solution. 48 8.1.2 The First Moment The existence theorem establishes that the global effective curvature field settles into a single, temporally locked m = 1 mode. To determine the strength and radial structure of this mode, we now identify the precise baryonic driver that enters the thin–disk Helmholtz equation of Section 8.1.3. This requires extracting, from the full time–dependent mass distribution, the unique component that survives temporal averaging and couples coherently to the m= 1 eigenpattern. We begin by modelling the baryonic emitters as point masses on circular orbits, each oscillating at the common proper–time frequency Ω, ρb(r, θ, t) = N X n=1 m⋆δ(r−rn)δθ−θn(t)e−i[Ωt+∆n(t)],(113) where the orbital phases evolve according to θn(t) = θn0+ω(rn)t, (114) and differential rotation contributes the drift ∆n(t) = −ω(rn)t. (115) All emitters thus share the same temporal frequency while acquiring radius–dependent phases; perfect temporal locking occurs only at radii satisfying ω(r0) = Ω.(116) To extract the coherent azimuthal structure, expand ρbin a Fourier series, ρb(r, θ, t) = ∞ X m=−∞ ρb,m(r, t)eimθ,(117) with ρb,m(r, t) = 1 2πZ2π 0 ρb(r, θ, t)e−imθ dθ. (118) Using (113) gives ρb,m(r, t) = N X n=1 m⋆δ(r−rn)e−i[Ωt+∆n(t)] e−imθn(t),(119) with total phase Ψm,n(t) = −Ωt−∆n(t)−mθn(t).(120) Substituting the orbital evolution, Ψm,n(t) = −mθn0+t[ω(rn)(1 −m)−Ω] .(121) 49 A mode remains phase–coherent only if ω(rn)(1 −m)−Ω = 0,(122) i.e. ω(rn) = Ω 1−m.(123) Only m= 1 allows a physical solution, ω(r0) = Ω,(124) while all m= 1 dephase under differential rotation. It is convenient to write the source in a form that makes the temporal frequencies explicit, ρb(r, θ, t) = ∞ X m=−∞ ˆρm(r) exp−i[(Ω −ω(r))t−(1 −m)θ],(125) with mode frequency ∆Ωm(r)=Ω−(1 −m)ω(r).(126) Only m = 1 satisfies ∆Ω 1 ( r ) = 0 in a flat–rotation disk, and therefore becomes strictly time–independent, ρb,1(r, θ, t) = ˆρ1(r)eiθ.(127) All other harmonics acquire oscillatory factors e−i∆Ωm(r)t and dephase. In a high– Q medium the effective curvature field responds only to the long–time average, ρb(r, θ, t)T−−−→ T→∞ ˆρ1(r)eiθ.(128) Thus the coherent driver of the effective curvature field is precisely the m = 1 component. On scales large compared with the interstellar spacing, the discrete source is replaced by a smooth surface density Σb(r, θ) with vertical profile Z0(z), ρb(r, θ, z)≃Σb(r, θ)Z0(z).(129) The azimuthal average and first moment are Σ0(r) = 1 2πZ2π 0 Σb(r, θ)dθ, M1(r) =Z2π 0 Σb(r, θ)e−iθ dθ. (130) We parametrize the coherent m= 1 component as Σb,1(r, θ) = ϵ(r) Σ0(r)eiθ,(131) so that M1(r)=2π ϵ(r) Σ0(r), ϵ(r) = M1(r) 2πΣ0(r).(132) 50 This coherent moment M1 ( r ) is the unique long–time driver of the global m = 1 spiral resonance. In Section 8.1.3, the stationary envelope equation is reduced to a 2D Helmholtz problem for the midplane effective curvature field. The source term appearing there is precisely ρb,1 ( r, θ ), or equivalently the radial moment M1 ( r ). The remainder of the analysis therefore proceeds by inserting the coherent first moment into the Helmholtz equation to obtain the radial envelope, the outgoing Green function, and ultimately the full spiral morphology. 8.1.3 The Thin–Disk Model: Solving the 2D Helmholtz Equation Reveals the Halo The coherent first moment M1 ( r ) obtained in the previous subsection is the unique long–time driver of the global m = 1 curvature mode. We now determine the spatial structure of the corresponding field by solving the midplane Helmholtz problem. This constitutes the projection of the full three–dimensional eigenmode established in the Theorem in Section 8.1.1 onto the thin disk. 4 The thin–disk treatment represents the limiting case ℓz/Ldamp → 0, providing an accurate description of the galactic midplane and capturing the dominant lateral structure of the coherent m =1 spacetime mode. The 3D model will be built from this in the section that follows. In this limit ( h≪r ), the stationary envelope equation (69) reduces in the midplane to a two–dimensional Helmholtz form, ∂2 ∂r2+1 r ∂ ∂r +1 r2 ∂2 ∂θ2+k2 ΩA(r, θ) = 4πG c2 eff ρb(r, θ),(133) where ρb(r, θ) denotes the projected baryonic density of the disk. Expanding both the curvature envelope and the baryonic source in azimuthal Fourier modes, A(r, θ) = ∞ X m=−∞ Am(r)eimθ, ρb(r, θ) = ∞ X m=−∞ ρm(r)eimθ,(134) is justified because any square–integrable function on the circle ( θ∈ [0 , 2 π )) admits a complete Fourier representation at each radius r . This decomposition follows from the intrinsic 2π–periodicity in θ. The first azimuthal coefficient, ρ1 ( r ), corresponds to the same phase–weighted driver M1 ( r ) introduced in Eq. (132) . Substituting the expansions (134) into Eq. (133) yields, for each mode m, d2Am dr2+1 r dAm dr +k2 Ω−m2 r2Am=4πG c2 eff ρm(r).(135) 4 In this projection the envelope A ( r, θ ) represents the lateral component of the full mode A(r, θ, z)=A(r, θ)Z0(z), where Z0(z) is the lowest vertical eigenfunction of scale height ℓz. 51 The homogeneous solutions of Eq. (135) are the cylindrical Bessel functions Jm ( kΩr ) and Ym ( kΩr ), or equivalently the outgoing Hankel function H(1) m = Jm + iYm , which represents a radiating, causally damped wave (see Appendix A for definitions and asymptotic forms). The coherent spiral attractor corresponds to the dominant m =1 azimuthal mode. Identifying the stationary envelope as A ( r, θ ) = eiθM1 ( r ) / (Γ + i Ω) in Eq. (133) yields the radial Helmholtz equation for M1(r): d2M1 dr2+1 r dM1 dr +k2 Ω−1 r2M1=4πG(Γ+iΩ) c2 eff ρ1(r).(136) The unique outgoing–wave Green function for the 2D Helmholtz operator, consistent with the 3D normalization (∇2+k2 Ω)G=−4πδ3(r), is G1(r, r′) = iπ J1(kΩr<)H(1) 1(kΩr>). This normalization ensures that the far–field asymptotic behaviour A ( r ) ∝ ( GM1/c2 eff ) eikΩr/r matches the three–dimensional outgoing solution of Eq. (70) . The corresponding Green–function solution of Eq. (136) is A(r, θ) = eiθ 4π2iG c2 eff "J1(kΩr)Z∞ r r′M1(r′)H(1) 1(kΩr′)dr′ +H(1) 1(kΩr)Zr 0 r′M1(r′)J1(kΩr′)dr′#, (137) The first integral represents the interior (source–dominated) region that remains regular at r = 0, while the second describes the exterior response that decays outward and satisfies the outgoing (Sommerfeld) boundary condition. It is convenient to define the complex radial amplitude Ar(r) = 4π2iG c2 eff "J1(kΩr)Z∞ r r′M1(r′)H(1) 1(kΩr′)dr′ +H(1) 1(kΩr)Zr 0 r′M1(r′)J1(kΩr′)dr′#, (138) so that the two–dimensional field can be written compactly as A(r, θ) = eiθ Ar(r).(139) This separation makes explicit that all geometric and resonant information resides in the radial envelope Ar ( r ), while the azimuthal dependence eiθ simply encodes the m =1 rotation. The same definition of Ar ( r ) will be used in Section 8.1.4 to construct the three–dimensional field A ( r, θ, z ) and in Section 8.1.6 to evaluate the stored energy and effective dark–mass distribution. 52 8.1.4 Extending the Thin–Disk to Three Dimensions The two–dimensional Helmholtz formulation derived above captures the dominant radial and azimuthal structure of the resonant halo, but a complete description of the coherent m = 1 mode requires inclusion of the finite vertical extent of the disk. In cylindrical coordinates ( r, θ, z ) the stationary curvature envelope satisfies the full three–dimensional Helmholtz equation, ∂2 ∂r2+1 r ∂ ∂r +1 r2 ∂2 ∂θ2+∂2 ∂z2+k2 ΩA(r, θ, z) = 4πG c2 eff ρb(r, θ, z),(140) where ρb ( r, θ, z ) is the baryonic source confined to the thin midplane and its immediate surroundings. The effective curvature field extends beyond this layer in z , but its vertical profile is limited not by the temporal damping length Ldamp = ceff/ Γ, but by a geometric scale height Lz that reflects the finite thickness of the resonant layer itself. Assuming a separable form A ( r, θ, z ) = Ar ( r, θ ) Z ( z ) leads to the vertical eigenvalue equation d2Z dz2+k2 zZ= 0,(141) in which the vertical wavenumber kz is determined from the local dispersion relation k2 Ω=k2 r+m2 r2+k2 z,so that k2 z=k2 Ω−k2 r−m2 r2.(142) Here kr ( r ) denotes the local radial wavenumber of the midplane solution, which may be estimated in WKB form as k2 r≃ |∂rln Ar|2 . Both kr and kz vary only slowly with r , so the product ansatz A ( r, θ, z ) ≃Ar ( r, θ ) Z0 ( z ; r ) remains accurate in the adiabatic (thin–disk) limit. The lowest–order, even vertical solution consistent with midplane symmetry and exponential decay at large |z|is Z0(z;r) = e−|z|/Lz(r) p2Lz(r),(143) where the vertical scale height follows from the evanescent branch of Eq. (142) as L−1 z(r) = pℜ[−k2 z(r)].(144) The normalization in Eq. (143) ensures R∞ −∞|Z0|2dz = 1, so that R|A|2dV is identical to the midplane two–dimensional integral of |Ar|2 . Typically Lz≫ℓz but Lz≪Ldamp , so the effective curvature field remains vertically stratified rather than filling the full damping length. The same Green–function normalization used in the thin–disk solution - see Eq. (137) - is preserved here, ensuring that the three–dimensional extension retains the correct far–field amplitude and outgoing–wave behavior. 53 The complete three–dimensional eigenmode is therefore A(r, θ, z) = eiθ e−|z|/Lz(r) p2Lz(r)Ar(r),(145) where Ar ( r ) is defined in Eq. (138) . Equation (145) thus defines the properly normalized three–dimensional curvature envelope of the fundamental ( m, n ) = (1 , 0) mode. The prefactor combines the global phase and vertical confinement, while the radial amplitude Ar ( r ) encapsulates the full horizontal response of the resonant disk through the baryonic source distribution M1 ( r ). This vertically self– consistent, azimuthally rotating Laguerre–Gaussian mode therefore represents the fundamental three–dimensional eigenfunction of the earlier Theorem, with Ar ( r ) setting the observable halo morphology and Z0 ( z ; r ) defining a finite, physically determined scale height Lz(r) that follows from the local dispersion relation. 8.1.5 Gain from Incoherent Foldy–Lax Scattering The global m = 1 Laguerre–Gaussian mode is selected coherently at the level of the baryonic source through the first azimuthal moment ρ1 ( r ), which extracts the eiθ rotational harmonic associated with the one–armed spiral. Once launched, however, the subsequent build–up of the curvature field within the galactic medium is not governed by fully coherent interference but by the incoherent multiple–scattering dynamics described in Section 5.2. The long propagation times implied by the reduced wave speed ceff ≪c introduce order–unity phase delays between different scattering paths at the carrier frequency Ω, placing the system firmly in the regime where Foldy–Lax theory predicts intensity addition rather than coherent amplitude summation. Accordingly, the net resonant amplitude obeys the scaling |Aeff| ∝ Q|X|pN⋆ϵ, (146) rather than Q|X|N⋆ , reflecting the randomization of temporal phases from successive rescattering events. Crucially, this temporal incoherence does not destroy the m = 1 spatial structure of the mode. The reason is the rotational symmetry of the Helmholtz Green function that governs curvature propagation in the thin–disk geometry. In polar coordinates the kernel depends only on the separation angle, GΩ(r, θ;r′, θ′) = GΩ(r, r′, θ −θ′),(147) and is therefore invariant under global rotations ( θ, θ′ ) 7→ ( θ + φ, θ′ + φ ). This invariance implies that the propagation operator commutes with the azimuthal generator ∂/∂θ , so the Fourier harmonics eimθ form an orthogonal eigenbasis. Consequently, a field of the form A(0)(r, θ) = A(0) r(r)eiθ,(148) 54 scatters according to A(1)(r, θ) = ZGΩ(r, r′, θ −θ′)A(0) r(r′)eiθ′r′dr′dθ′.(149) Performing the θ′ -integration isolates only the m = 1 harmonic of the kernel, yielding A(1)(r, θ) = eiθ A(1) r(r),(150) with the angular dependence eiθ preserved exactly. Since the kernel contains no m = 1 components, no azimuthal mixing can occur unless the medium itself breaks axisymmetry. Iterating the Foldy–Lax expansion therefore maintains the m = 1 dependence at every scattering order, even though the temporal phases of successive generations of the field become uncorrelated. The resulting physical picture is that the mode ( m = 1) is selected coherently by the stellar disk, but the magnitude of that mode is set by incoherent, phase–randomized multiple scattering in the curvature medium. Rotational symmetry preserves the angular structure, while the long propagation delays enforce the √N⋆ Foldy–Lax scaling encoded in Eq. (146) . Thus the emergent one–armed spiral is a global resonant state whose geometry is fixed by coherent baryonic forcing and symmetry, but whose amplitude is determined by the statistical accumulation of many incoherent scattering events within the galactic environment. 8.1.6 Stored Energy, Dark–Mass Equivalence, and Environmental Gain The observable (quasi–static) curvature bias arises from the stored oscillatory energy of the resonant field. From Section 7.2, the local field energy density and its mass equivalent are EΦ(x) = 1 8πGhc2 eff |∇A|2+ Ω2|A|2i, ρeff(x) = EΦ(x) c2.(151) Substituting the normalized field A ( r, θ, z ) of Eq. (145) and separating variables gives |∇A|2=∂rAr 2|Z0|2+|Ar|2 r2|Z0|2+|Ar|2∂zZ0 2.(152) For the exponential vertical profile Z0 ( z ; r ) = e−|z|/Lz(r)/p2Lz(r) , we have ∂zZ0 2≃ |Z0|2/L2 z. The energy density therefore becomes ρeff(r, θ, z) = |A(r, θ, z)|2 8πGc2hΩ2+c2 eff∂rln Ar 2+1 r2+1 L2 zi.(153) In the adiabatic WKB limit, the local dispersion relation k2 Ω = k2 r + m2/r2 + k2 z (with the vertical term replaced by its evanescent magnitude 1 /L2 z ) renders the 55 bracket approximately constant: Ω2+c2 effk2 r+m2 r2+1 L2 z≃Ω2+c2 effk2 Ω≈2 Ω2,(154) since kΩ= Ω/ceff and Γ≪Ω. Thus, to good approximation, ρeff(r, θ, z)≃2 Ω2 8πGc2|A(r, θ, z)|2.(155) To include the cumulative effect of multiple scattering in a statistically homogeneous ensemble, the field amplitude is scaled by the environmental reinforcement factor, A(r, θ, z)−→ Aeff(r, θ, z) = QpGenv A(r, θ, z),(156) where Q is the temporal quality factor of the curvature resonance, and Genv represents statistical amplification of field amplitude through repeated scattering and re–emission. For a long–wavelength field ( λΩ≫ mean inter–scatterer spacing), individual scatterers respond additively in amplitude, yielding Genv =N⋆|X|2.(157) where N⋆ is the number of baryonic scatterers within the effective coherence volume and |X| ≤ 1 quantifies their mean coupling efficiency to the global oscillation frequency Ω. The factor |X| thus measures how strongly each scatterer participates in the collective resonance. Substituting Eq. (156) into Eq. (155) gives the final expression for the effective energy density: ρeff(r, θ, z) = 2 Ω2 8πGc2Genv Q2|A(r, θ, z)|2,Genv =N⋆|X|2.(158) The result (158) therefore provides the complete description of the statistically amplified dark–mass density in terms of the intrinsic mode amplitude, resonance quality factor, and the effective number of coupled baryonic scatterers. 8.1.7 Numerical Results: 3D Thin–Disk Model The complete numerical results obtained from the three–dimensional thin–disk model of Section 8.1.4 are shown in Figures 4–6. All quantities and algorithms are implemented exactly as specified in the simulation code given in Appendix B. The rotation frequency is fixed by the flat circular speed v0 = 230 km s−1 at a reference radius r0 = 8 . 5 kpc , giving Ω = 8 . 77 × 10 −16 s−1 and a resonant wavelength λres ≃ 46 . 44 kpc . An effective propagation speed ceff = 2 × 10 5m s−1 is adopted, while the quality factor Q = 12 . 6 is chosen such that the resulting damping length Ldamp = ceff/ Γ ≃ 192 . 17 kpc (where Γ = Ω / (2 Q )) reproduces 56 the observed halo extent and the asymptotic flattening of typical spiral–galaxy rotation curves. In this formulation the total resonant amplitude arises directly from the combination of the temporal coherence Q , coupling efficiency X , and geometric projection factor ϵ , giving an overall scaling proportional to Q|X|√N⋆ϵ . Unlike earlier parameterizations in which external gain was inserted by hand, these parameters now enter self–consistently through the integral Helmholtz solution. The resulting field strengths, density contrasts, and rotation curves remain identical in form to those of the original model, confirming the internal normalization of the updated implementation. Figure 4 presents the one–dimensional radial solutions derived from the coherent m = 1 integral formulation. Part (a) shows the field amplitude |A ( r ) | , which rises from zero at the origin to a well–defined maximum near r≃ 6 kpc , before declining exponentially on scales comparable to the damping length. Part (b) displays the logarithmic mass densities: the effective curvature density ρeff (blue) and the baryonic midplane density ρb (orange dashed). The effective component decreases almost linearly in log–space and approaches the background level near r≃ 650 kpc . Part (c) compares the baryonic, resonant, and total rotation curves. The total velocity remains approximately flat beyond r∼ 10 kpc , consistent with observed spiral–galaxy kinematics. Part (d) shows the enclosed–mass ratio Meff/Mb , which increases from ∼ 10 at 7 kpc to ∼ 14 at 100 kpc , ∼ 17 at 150 kpc , and approaches its asymptotic value of ≃ 20 by r≃ 200–250 kpc . These results confirm that the thin–disk model recovers both the outer flattening of the rotation curve and the canonical dark–to–baryon mass ratio inferred for Milky–Way–like galaxies. Figure 5 illustrates the two–dimensional midplane structure of the same m = 1 resonant mode over a 100 × 100 kpc domain. Part (a) shows the normalized amplitude |A| , revealing the expected toroidal morphology peaking at r≃ 6 kpc . Part (b) displays the normalized azimuthal phase arg ( A ), which forms a continuous one–armed spiral with a single 2 π winding, indicating a coherent global m = 1 mode. Part (c) overlays the baryonic (red) and resonant (green) densities, showing their partial spatial displacement within the disk plane. Part (d) shows the total (baryonic + resonant) density on a linear scale; only the bright baryonic core and nearby resonant envelope are visible, as the extended halo falls below the display contrast. Figure 6 shows the corresponding three–dimensional effective–density distribution. Part (a) presents the horizontal (XY) midplane slice, and part (b) the vertical (XZ) cross–section, each over a 200 × 200 kpc domain. The XY projection again displays the toroidal halo structure, whereas the XZ projection reveals a vertically stratified envelope extending to heights of order |z|∼Lz≃ 57 . 65 kpc , comparable to the resonant wavelength λres ≃ 46 . 44 kpc . This scale follows from the evanescent vertical branch of the disk dispersion relation and indicates that the effective curvature field remains confined within a few wavelengths of the midplane rather than filling the entire damping region. The resulting 57 analogue is found in the family of low–order, self–gravitating oscillation modes briefly reviewed in Section ?? and by Cox [ 3 ]. These modes—including the classical radial and low-degree ( ℓ = 1–2) “breathing” oscillations of extended gravitating systems—possess geometries tightly constrained by the underlying mass distribution (e.g. disk-like, spheroidal, or ellipsoidal profiles). Although such stellar and cluster modes are not known to exhibit global phase coherence, their existence demonstrates that large-scale, geometry-aligned gravitational oscillations are dynamically admissible. If analogous motions within a rotating galactic disk retain even partial coherence over many orbital periods, they would naturally produce the predominantly planar, m = 1–polarised forcing required by the curvature resonance developed here. Thus the existence of these geometrically constrained low-order modes provides a plausible physical basis for the baryonic sourcing of the coherent m= 1 curvature resonance in spiral galaxies. Finally, the resonant curvature mode behaves as a damped oscillator with natural frequency Ω and damping rate Γ, so that microscopic dissipation removes effective gravitating density at the rate Ploss = 2Γ ρeff . In the steady state this loss is exactly replenished by coherent orbital driving from the baryons, and the resulting mode amplitude is fixed by the combined action of four physical factors: the quality factor Q = Ω / (2Γ), the single–scatterer coupling strength |X| , the geometric m =1 projection factor ϵ , and the statistical Foldy–Lax gain √N⋆ from the N⋆ baryonic emitters within a coherence volume. Together these give the scaling Aeff ∝Q|X|pN⋆ϵ, (161) expressing that the halo is a self–maintaining resonant structure whose long-term amplitude is determined by the balance between orbital energy injection and dissipation at rate Γ. 8.2 Ellipticals as Random Resonant Media 8.2.1 Observing the Statistical Steady State Elliptical galaxies represent the natural continuation of the resonant framework developed for spiral systems in Section 8.1, but in a regime where global coherence is replaced by statistical transport. Spirals maintain a coherent m =1 curvature mode through organized rotation and phase–locked rescattering; ellipticals, by contrast, lack any global azimuthal order. The random orbital motion of stars continually scrambles relative phases, so that coherence survives only locally while the ensemble as a whole behaves as a diffusive resonant medium. Curvature energy is still stored and re–emitted through the same causal mechanism, but now in a statistically fluctuating rather than phase–locked form. This transition from coherent to diffusive resonance is the same one described in radiative–transport theory for acoustic and elastic wave fields [ 13 , 14 , 9 , 4 , 5 ], where random rescattering replaces global phase locking with ensemble–averaged diffusion. 64 In this limit, the coherent multiple–scattering picture described by the Foldy Lax hierarchy (Section 5.2) must be replaced by its statistical extension—the effective–medium and transport framework developed in Sections 5.3 and 5.4. Each star acts simultaneously as source and scatterer of curvature waves, but its motion introduces stochastic phase drift on the scale of the resonant period. The resulting field is a dense superposition of delayed, phase–drifting phasors whose interference forms a curvature speckle field: a dynamically evolving, granular intensity pattern produced by the mutual scattering of moving masses. This pattern is statistically analogous to optical speckle [ 15 ]. In spiral galaxies the response is dominated by a single, globally coherent m =1 mode with a well–defined complex amplitude, whereas in the elliptical regime no single standing–wave mode dominates: the dynamically relevant quantity is the ensemble–averaged intensity built from the incoherent superposition of many phasors. The continual overlap and randomization of these phasors drive a diffusive redistribution of curvature energy throughout the ensemble, as described by the radiative–transfer and diffusion relations of Section 5.4. The corresponding ensemble parameters (Γ eff, Ldamp,eff, Qeff ) govern the macroscopic propagation and memory of the effective curvature field. Because curvature energy is repeatedly rescattered and reused, the effective damping rate Γ eff becomes much smaller than the microscopic Γ, and both the memory time and damping length ( Qeff, Ldamp,eff ) are strongly enhanced. These extended coherence measures explain the smooth, long–tailed 1/r2intensity profiles characteristic of elliptical and cluster halos. Elliptical galaxies therefore act as statistical curvature resonators: locally coherent, globally diffusive systems in which curvature energy is perpetually recycled through the random rescattering of the very masses that sustain it. Their apparent “dark” halos arise from the ensemble–averaged intensity of this self–generated curvature speckle field—a dynamically sustained, statistical extension of the same causal resonance that governs spirals. Theorem (Statistical Steady–State Theorem for Ellipticals). Consider a pressure–supported stellar system lacking global rotational order, in which the resonant effective curvature field is described by a complex envelope A (x , t )and an associated effective gravitating density ρeff (x , t ). Let ukin (x , t )denote the stellar kinetic–energy density. Then the curvature–feedback loop reaches a statistical steady state when the ensemble–averaged effective density obeys d dt ρeff=G−L= 0,(162) where G and L denote the local generation and loss rates of effective gravitating density. These satisfy G ∝ ukin,L= 2Γeff ρeff,(163) with Γeff the transport (multiple–scattering) damping rate. 65 In the source–free exterior region, where the effective curvature field becomes statistically isotropic and diffusive, the steady–state solution implies an ensemble– averaged density profile ρeff(r)∝exp[−2r/Ldamp,eff ] r2,(164) with transport damping length Ldamp,eff =ceff Γeff ,(165) where ceff is the effective curvature–wave speed. Interpretation. Equation (162) expresses the condition for statistical stationarity: the effective gravitating density generated by random stellar motion balances the density removed by transport damping. Equations (164) – (165) then give the resulting spatial form—a transport–regulated curvature halo with an r−2 intensity tail and an exponential cutoff set by Ldamp,eff. 8.2.2 Complex Baryonic Forcing with Memory A key result of this book is that the curvature–envelope field A (x , t ) is the phononic excitation of the resonant effective curvature field of spacetime, sourced and scattered by baryons. 5 Throughout what follows we therefore interpret A as a continuous field of curvature phasors whose interference carries stress and energy through the system. In regions where these phasors reinforce coherently, the field forms localized packets of compression and rarefaction—phonon-like curvature wavepackets which we refer to as massons (Section 4.3). These are not material quanta but localized solutions of the same linear, causal PDE derived earlier. The governing dynamics are elastic, obeying the same forced–damped Helmholtz equation and transport identities as in acoustic multiple scattering. One may picture a vast ensemble of randomly distributed resonant scatterers, each acting as a kinematic source while simultaneously obeying a common forced–damped oscillator response that defines χ (Ω), and therefore both emitting and resonantly scattering curvature phonons. This is directly analogous to the resonant–rod experiments of Derode et al. [ 4 , 5 ], in which multiple scattering of diffuse ultrasound produces coherent amplification and long–range interference within a random elastic network [16, 17, 18]. 6 5 Recall that A (x , t ) is the slowly varying complex amplitude (envelope) of the resonant effective curvature field, Φ res = A e−iΩt + c.c. . It represents a continuous field of curvature phasors arising from the same forced–damped oscillator equation that governs single-mass susceptibility χ(Ω). 6 In the acoustic case the rods are driven by an external source; in the gravitational system each mass element is a self–driven kinematic emitter, while the same causal forced–damped PDE determines how its emission propagates, scatters, and damps. 66 In such a scattering environment, the finite transport mean free path reduces the spatial coherence of the resonant effective curvature field while simultaneously extending its effective damping length. These two effects (shorter coherence, longer diffusive decay) govern how curvature energy is redistributed through the medium and ultimately define the smooth, anisotropic envelope of the effective density field. As shown below, this transition from coherent to diffusive propagation is precisely what gives rise to the extended, elliptical morphology of the dark–mass halo. In spiral galaxies, the low–entropy m =1 orbital order enforced a nearly uniform geometric phase. The complex source term mneiϕn was therefore a coherent phasor sum: the driving was essentially instantaneous and spatially phase–locked (see Section 8.1). By contrast, elliptical galaxies host a phase–mixed stellar population. The geometric phase ϕn evolves stochastically due to random stellar velocities, destroying long–range coherence. Coherence survives only locally in space (transverse scale ℓ⊥) and in time (memory scale τcoh). The correct source of curvature phasors in an elliptical is therefore the phasor representation of the baryonic forcing term in the same causal PDE: ρb(x, t) = N⋆ X n=1 m⋆eiϕn(t)δ(3) x−xn(t),(166) which drives the forced–damped Helmholtz operator ∇2+k2 ΩA(x, t) = 4πG c2 eff ρb(x, t),(167) where m⋆ is the mass of a representative emitter, x n ( t ) its instantaneous position, and ϕn ( t ) its geometric phase relative to the local field. This is the same driven, damped equation used in the coherent (spiral) case; only the phase statistics differ. The field amplitude A (x , t ) may be written using the causal (single–scatterer) Green function: A(x, t) = 4πG c2 eff N⋆⋆ X n=1 m⋆eiϕn(t)GΩ x−xn(t),(168) The causal Green function GΩ (r) describes the response to a single resonantly driven stellar source. Its amplitude decays in the far field as |GΩ(r)|2=e−2|r|/Ldamp (4π)2|r|2, Ldamp =ceff Γ,(169) representing the causal damping of an individual curvature phasor. Ensemble averaging over many such phasors will later produce the transport–renormalized kernel. 67 In Eq. (168), N⋆⋆ denotes the full population of past emitters within the causal memory window t−τmem < t′≤t, where τmem = Γ −1 is the memory time of the underlying forced–damped oscillator response and τcoh = (2Γ) −1 its shorter phase–decorrelation time. Equation (168) therefore represents a four–dimensional ensemble of phasors, each carrying its causal phase drift across the memory window. Not all of the N⋆⋆ phasors are statistically independent: correlations persist over finite regions of space and time, defining a four–dimensional speckle scale determined by the coherence of the same underlying PDE. The following subsection therefore defines the characteristic size of a four–dimensional speckle cell, within which phasors remain phase–correlated. Determining this coherence volume and timescale allows the number of statistically independent phasors within N⋆⋆ —and the corresponding baryonic mass— to be quantified. These quantities are then used in Section 8.2.4 to rewrite Eq. (168) in terms of the independent phasor ensemble and to compute ensemble averages such as ⟨|A|2⟩ and ⟨ρeff⟩. 8.2.3 Four–dimensional speckle geometry of the gravitational field To extract statistical predictions from Eq. (168) , we must determine over what regions of space and time the gravitational field A (x , t ) retains phase coherence. Because elliptical galaxies lack global orbital order, coherence is neither disk–wide nor steady in time: it is localized and transient. The field therefore forms a granular, four–dimensional speckle pattern—the gravitational analogue of diffuse ultrasound speckle. Transverse coherence. Random stellar motion produces geometric phase drift. When the phase difference between two emitters reaches one radian, their contributions become statistically uncorrelated. This defines the transverse coherence length ℓ⊥=λres 2π,(170) so phasors launched from within a transverse patch of radius ℓ⊥ interfere coherently. Longitudinal coherence. Along a star’s trajectory, orbital motion also scrambles phase in time. During the coherence interval τcoh the typical orbital excursion is ℓ∥=v τcoh =v 2Γ,(171) so emissions separated by more than ℓ∥along an orbit are decorrelated. Space–time coherence cell. The smallest region over which phasors retain a common phase is a transverse disk of radius ℓ⊥ extended over a longitudinal distance ℓ∥. Its spatial volume is Vcell =α ℓ2 ⊥ℓ∥=αλres 2π2v 2Γ,(172) 68 where α is a geometric factor of order unity. The corresponding temporal extent is ∆tcell =τcoh =1 2Γ.(173) This pair ( Vcell, ∆ tcell ) is the gravitational analogue of the isoplanatic patch and coherence time of diffuse acoustics. Population statistics. Let Vgal denote the active baryonic volume. The number of disjoint coherence cells at any moment is N3D =Vgal Vcell ,(174) and each cell contains on average Ncell =n⋆Vcell =N⋆Vcell Vgal ,(175) stars, where n⋆is the mean stellar number density. Effective phasor mass. Within a single cell all stars share, to leading order, a common field phase. Their masses therefore add coherently as one composite phasor of effective mass Mcell =Ncell m⋆.(176) The original stellar sum in Eq. (168) may thus be re-expressed as a sum over statistically independent four–dimensional coherence cells, each weighted by its complex mass Mcell and the appropriate Green kernel. In summary, coherence in ellipticals is local in space and finite in time. The gravitational field is driven not by a global mode but by a four–dimensional mosaic of statistically independent phasor ensembles. Rewriting Eq. (168) in this statistical basis enables the evaluation of ensemble quantities such as ⟨|A|2⟩ and the mean effective density ⟨ρeff ⟩, as developed in the following subsection. 8.2.4 Ensemble reduction: masson–speckle summation We now rewrite the forcing in Eq. (168) in terms of the statistically independent four–dimensional coherence cells derived in Section 8.2.3. Each spatial coherence volume Vcell = α ℓ2 ⊥ℓ∥ defines an elementary unit of the ensemble environment, within which the stellar phasors remain phase–locked and the local effective curvature field evolves coherently over the memory time τcoh ≃ 1 / (2Γ). Cells separated by more than ℓ⊥in space or ∆tcell in time are fully decorrelated. Stars within a common cell therefore act as a single coherent source of effective mass Mcell , while the full galaxy forms a random, time–varying mosaic of such locally coherent “masson” excitations. The effective curvature field can thus be expressed as a superposition of statistically independent phasors: A(x, t) = 4πG c2 eff N3D X c=1 NT X j=1 M(c,j) cell GΩ x−xc,(177) 69 where c indexes spatial cells and j indexes independent temporal realizations within the field’s memory window. Temporal multiplicity: curvature memory. Each spatial cell stores a finite number of statistically distinct phasors within its coherence time. Using the definitions of τcoh and the oscillation period 2 π/ Ω, the number of independent temporal realizations per cell is NT=τcoh ∆tcell = 1 2Γ 2π Ω =Ω 4πΓ=Q 2π,(178) so that each cell contributes NT independent temporal samples of its masson excitation. Intensity expansion. Taking the square modulus of Eq. (177) gives |A(x, t)|2=4πG c2 eff 2 X c,j X c′,j′ M(c,j) cell M(c′,j′) cell ∗GΩx−xcG∗ Ωx−xc′,(179) from which ensemble averages may be computed. Statistical independence of distinct cells and temporal windows gives DM(c,j) cell M(c′,j′) cell ∗E=(|Mcell|2, c =c′, j =j′, 0,otherwise, (180) so all cross terms vanish. The mean field intensity is therefore D|A(x, t)|2E=4πG c2 eff 2 NT N3D X c=1 |Mcell|2GΩ(x−xc) 2.(181) For a Poisson stellar distribution, |Mcell|2=N2 cellm2 ⋆+NcellVar(m⋆),(182) and for Ncell ≫1 this simplifies to |Mcell|2≃M2 cell =N2 cellm2 ⋆.(183) Using N3D =Vgal/Vcell and Ncell =N⋆Vcell/Vgal,we obtain N3D M2 cell =Vgal Vcell N⋆Vcell Vgal 2 m2 ⋆=N2 ⋆m2 ⋆ Vgal Vcell.(184) It is now convenient to collect the ensemble coupling terms into a single dimensionless environmental gain factor Genv . In Section 5.1 this quantity was defined microscopically as the total intensity gain relative to the incident field reducing in the fully coherent Foldy–Lax limit to Genv = N2|χ⋆|2 [Eq. (43) ]. In 70 the effective–medium picture of Section 5.3, the same quantity was shown to scale inversely with the coherence volume Vcell through Genv ∝ 1 /Vcell reflecting the fact that larger coherence cells produce stronger local phase alignment. In the present, fully explicit cell–counting picture we include both the local coupling per cell and the number of cells filling the galaxy. To keep track of the mean stellar coupling efficiency introduced earlier, we write the ensemble environmental gain as Genv := |X|2N3D |Mcell|2 m2 ⋆≃ |X|2N3D N2 cell =|X|2N2 ⋆Vcell Vgal =|X|2N2 ⋆ Vgal α ℓ2 ⊥ℓ∥ =|X|2N2 ⋆ Vgal αλres 2π2v 2Γ,(185) where Vcell = α ℓ2 ⊥ℓ∥ is the spatial coherence volume from Eq. (172) , ℓ⊥ = λres/ (2 π ) is the transverse coherence length (Eq. (170) ), and ℓ∥ = v/ (2Γ) is the longitudinal coherence length (Eq. (171) ). The dimensionless factor |X|2≤ 1 represents the mean stellar coupling efficiency already introduced in the scattering theory: in the coherent Foldy–Lax limit it multiplies the intensity as Genv ∝|X|2N2 , while here it weights the total cell–ensemble gain in exactly the same way. These successive equalities make explicit how the ensemble environmental gain retains the same conceptual structure as in the scattering theory. The first line shows the statistical definition, where the coherent mass content of each cell, |Mcell|2 , is weighted by the mean coupling efficiency |X|2 and multiplied by the number of spatial and temporal realizations. The second line reformulates this in geometric terms, introducing the total stellar population N⋆ and the ratio Vcell/Vgal , which connects the microscopic cell statistics to the macroscopic galactic volume. The third line then expresses the same quantity in explicit physical scales, replacing the coherence lengths ( ℓ⊥, ℓ∥ ) with the resonant wavelength λres and damping time 1 / Γ, thereby revealing how Genv depends on the underlying causal and geometric parameters of the medium. Thus the earlier local relation Genv ∝ |X|2/Vcell for a single coherent patch is naturally extended here to the global ensemble gain: the per–cell response scales as |X|2/Vcell , while the total number of coherently contributing sources scales as N2 ⋆Vcell/Vgal . Their product yields the dimensionless Genv relevant for ellipticals. The replacements Γ → Γ eff and Ldamp →Ldamp,eff will be introduced in the following subsection when transport effects are included. Substituting Eqs. (178), (184), and (185) into Eq. (181) finally gives D|A(x, t)|2E=4πG c2 eff 2Q 2πGenv GΩ(x−xc) 2,(186) 71 which makes explicit how the environmental gain Genv , the temporal memory Q/ (2 π ), and the Green–function geometry combine to determine the observable curvature–field intensity in a statistically mixed stellar system. This result forms the starting point for the transport–enhanced description developed in the next subsection. 8.2.5 Transport Mean Free Path and Collective Damping The ensemble formulation developed above treats the effective curvature field as a superposition of kinematically driven stellar sources, each of which contributes phase-modulated curvature phasors within its local coherence window. In realistic galactic environments, however, each star acts simultaneously as (i) a kinematic source of curvature phonons generated by its orbital motion, and (ii) a compact elastic scatterer governed by the same forced–damped Helmholtz operator that propagates and damps the field. Through continual emission, scattering, and reabsorption, these mass elements generate a diffusive interference network whose statistics mirror those of resonant ultrasound and diffuse field correlations in complex materials [ 13 , 9 , 4 , 5 , 16 , 17 , 18 ]. The governing dynamics are elastic, obeying the same damped Helmholtz equation and transport identities as acoustic phonons in a diffuse medium. Following the formalism of multiple–scattering theory [ 8 , 19 ], the coherent field propagation is not determined by the microscopic damping rate Γ but by an effective transport wavenumber that incorporates the ensemble self–energy Σ(Ω): k2 eff =k2 0+ Σ(Ω), k0=Ω+iΓ ceff .(187) The self–energy Σ(Ω) encodes the cumulative phase delay and dwell time produced by the resonant multiple-scattering network, not by the emitters themselves. Its imaginary part defines the transport damping rate Γeff =ceff Im keff,(188) and the associated effective damping length Ldamp,eff =ceff Γeff .(189) In the diffusive regime, Γeff ≪Γ, Ldamp,eff ≫Ldamp =ceff Γ,(190) expressing the key transport result: resonant multiple scattering shortens spatial coherence but greatly extends the effective damping length. Curvature phonons undergo repeated dwell–time interactions before decorrelating, continuously recycling part of their energy back into the coherent mode. This recycling produces an effective quality factor Qeff =Ω 2 Γeff ≫Q=Ω 2 Γ,(191) 72 so that the observed large Q values in elliptical halos are not microscopic properties of single emitters but collective transport effects of the stellar ensemble. In this regime, all ensemble relations of Section 8.2.4 retain their formal structure under the replacements Γ→Γeff, Q →Qeff, Ldamp →Ldamp,eff,(192) and the environmental gain itself becomes transport–renormalized through the effective longitudinal coherence length ℓ∥,eff =v/(2Γeff ), yielding Genv,eff =|X|2N2 ⋆ Vgal αλres 2π2v 2Γeff .(193) This form preserves the same microscopic structure as Eq. (185) , but with Γ replaced everywhere by its transport value Γ eff to account for the collective damping reduction of the ensemble medium. The angularly averaged Green function of the multiply scattered curvature phonon field defines the causal transport kernel Keff(r) = exp−2r/Ldamp,eff  (4π)2r2,(194) which represents the spherically symmetric energy envelope of the coherent curvature mode. It provides the physical basis for the closed halo law derived in the next subsection: a long–range, slowly decaying effective curvature field whose extended elliptical morphology arises naturally from the diffusive recycling of curvature energy within the transport medium. Substituting Eqs. (193) and (192) into the ensemble intensity (186) , and replacing the microscopic propagator |GΩ|2 by the transport kernel Keff ( r ), yields the compact transport–enhanced expression D|A(r)|2E=4πG c2 eff 2Qeff 2πGenv,eff Keff(r).(195) This form makes the transport renormalization of both temporal and spatial coherence explicit: Genv,eff now embodies the reduced damping rate Γ eff , while Qeff and Keff ( r ) encode the extended memory time and spatial transport length. It provides the physical closure for the ensemble halo law, from which the core regularization and gradient–energy correction will follow. 8.2.6 Core Regularization in the Non–Transport Regime Before introducing the full energy conversion in Section 8.2.7, we first establish the spatial morphology implied by the transport solution for the field intensity ⟨|A ( r ) |2⟩ . At this stage, only the relative scaling of the envelope is needed: the absolute normalization and physical prefactors from the curvature–energy law 73 8.2.10 Physical Interpretation Elliptical galaxies represent the self–sustaining limit of the transport–resonant curvature–phonon framework developed above. In these pressure–supported systems the random three–dimensional stellar motions both generate and scatter the effective curvature field. Each star acts as a moving phase perturbation of the local spacetime curvature, and the ensemble of such perturbations forms a dense, multiply–scattering phononic medium governed by the Dyson self–energy Σ(Ω). Kinetic energy injected by the baryons is continually converted into curvature phonons, which scatter back upon the baryons and help maintain the velocity dispersion that excites them. The galaxy as a whole thus behaves as a gigantic, statistically stationary, self–excited transport resonator in which microscopic generation and ensemble damping remain in dynamic equilibrium. In the diffusive limit, the relevant resonant scale is determined by the density and geometry of the coherent domains—or “4D speckle cells”—that compose the stellar ensemble. Each cell of characteristic size λres encloses a coherent stellar mass Mcell within a coarse–graining kernel of volume Vcell ∼λ3 res . Within such a domain, phase correlations cause curvature amplitudes to add linearly while ensemble–averaged intensities add quadratically, yielding ρeff ∝M2 cell Vcell c2,(211) where ρeff is the local ensemble–averaged effective gravitating density. The resulting profile therefore depends jointly on the number of coherent sources within a speckle cell and the compactness of the cell itself. A central modelling assumption of this framework is that the natural restoring frequency of the resonant curvature mode within any coherence domain scales with the square root of the coarse–grained baryonic density, Ω∝√¯ρb, where ¯ρb denotes the mass density averaged over the local coherence kernel. This assumption mirrors the well–established ω∼√Gρ scaling found in small–amplitude oscillations of self–gravitating systems (stellar pulsations, polytropes, Jeans modes). If this empirical behaviour continues to hold approximately within a finite coherence region, then the density of that region determines the resonant frequency that participates in global mode locking. Under this assumption, the corresponding resonant length scale follows as λres =ceff Ω∝¯ρ−1/2 b,(212) so regions of higher coarse–grained density support shorter resonant wavelengths. In ellipticals, the coherence kernel is effectively isotropic, its volume set by the short mean free path of curvature phonons in a dense, three–dimensional stellar field. The resulting coarse–grained ¯ρb is therefore large, producing small values 80 of λres . In spirals, by contrast, coherence extends anisotropically through the disk: the kernel is oblate and quasi–toroidal, its vertical thickness restricted by the stellar scale height. This larger coherence volume reduces ¯ρb and yields much longer resonant wavelengths. With ceff ≈ 2 × 10 5m s−1 fixed, the observed contrast between systems follows directly from geometry and density. Typical coarse–grained densities of ¯ρb∼ 10 −21 –10 −22 kg m −3 in ellipticals yield λres ∼ 20–30 pc. Disk–averaged densities around 10 −24 kg m −3 in spirals give λres ∼ 30–50 kpc. Cluster–scale cavities with ¯ρb≲ 10 −26 kg m −3 support modes approaching megaparsec scales. Thus the short wavelengths of ellipticals reflect their compact, three–dimensional coherence and high baryonic density, whereas spirals and clusters express progressively larger, lower–density resonant domains. The issue of sustainment is not concerned with the static ensemble profile but is related to the persistence of the microscopic fluctuations that continually regenerate the effective curvature field. In the diffusive–transport regime each scattering event introduces a small delay and phase randomization, preventing the instantaneous field from ever relaxing into a perfectly smooth configuration. The central regions therefore exhibit a shifting speckle pattern of bright and dim curvature zones—localized fluctuations evolving on the coherence timescale τcoh,eff ∼ 1 / (2Γ eff ) (Section 8.2.5). Although minuscule compared with the stellar kinetic energy, these fluctuations form a microscopic phase–randomizing bath that continuously refreshes the stochastic component of the halo. On macroscopic scales this agitation establishes a steady exchange between baryonic kinetic motion and the effective gravitating density carried by curvature modes. In ensemble average the balance takes the form d⟨ρeff⟩ dt =⟨G⟩−⟨L⟩ ≃ 0,(213) where the generation rate is proportional to the local kinetic–energy density, G ∝ ukin , and the loss rate is L = 2Γ effρeff . Stationarity therefore requires that random stellar motion replenish the effective density at precisely the rate at which ensemble damping removes it. This is the same self–consistent energy closure that fixes the transport damping rate Γ eff in the Dyson effective–medium theory. This balance is directly analogous to the classical Jeans equilibrium. In the Jeans picture random stellar velocities generate a pressure gradient balancing self–gravity. Here, the same random motions generate curvature–phonon energy ( G ∝ ukin ), while the resulting effective density feeds back through the transport damping term L = 2Γ effρeff . The steady state is thus the resonant–transport counterpart of Jeans support: kinetic agitation sustains the effective curvature field, and the effective curvature field confines the very motions that regenerate it. 81 Because the speckle bath is continually refreshed, the ensemble–averaged effective density behaves as a causal moving average of the instantaneous field: ⟨ρeff⟩t+∆t= (1 −β)⟨ρeff⟩t+β ρinst(t), β = 1 −e−∆t/τcoh,eff ,(214) with τcoh,eff determining the memory time of the transport process. In stationary regions the average converges to a constant value; in evolving environments it governs how the ensemble field tracks transient curvature fluctuations. Through this continual renewal on the coherence timescale, the system maintains a persistent background of stochastic curvature stresses that both mirror and sustain the random stellar motions themselves—thus closing the feedback loop between baryonic agitation and curvature diffusion. Although the mean effective curvature field is smooth and energetically subdominant, its finite variance implies a nonzero probability of localized curvature intensification. If such peaks persist for several coherence times they may weakly bias the nearby stellar distribution, producing transient clustering or faint substructure. These features arise naturally from the feedback loop itself: the very fluctuations that preserve equilibrium occasionally imprint brief, localized curvature enhancements within an otherwise relaxed halo. The polarisation properties of the curvature–matter resonance in ellipticals follow from the susceptibility tensor χij (Ω) introduced in Section 6.3.3. Although individual four–dimensional speckle cells are locally anisotropic—each having transverse and longitudinal scales ( ℓ⊥, ℓ∥ ) set by the stellar velocity distribution— the ensemble of cells is randomly oriented and statistically well mixed. Because the stellar velocity ellipsoid in pressure–supported ellipticals is nearly isotropic, this random orientation washes out the anisotropy of individual cells in the ensemble average. The linear response therefore reduces to the isotropic form χij(Ω) = χ(Ω) δij,(215) so all tensor eigenmodes share the same susceptibility, χ(x)=χ(y)=χ(z). No preferred polarisation axis exists in this limit. Three–dimensional stellar agitation excites all eigenvectors vi (α) with comparable strength, and rapid phase mixing removes any transient anisotropy. Ellipticals are therefore statistically unpolarised: their curvature–phonon field populates all tensor eigenmodes equally, consistent with the nearly spherical transport envelopes of pressure–supported systems. This framework naturally extends to lenticular, triaxial, and irregular systems. Whenever the stellar velocity distribution departs from isotropy, the ensemble– averaged susceptibility acquires a small but finite eigenvalue splitting, χxx =χyy =χzz,(216) with the principal axes aligned with those of the stellar velocity ellipsoid. The coherence kernel remains three–dimensional but becomes mildly oblate or triaxial, imparting a correspondingly weak polarisation bias to the global effective 82 curvature field. Such systems lie between the isotropic (giant elliptical) and strongly planar (spiral) limits. In S0 galaxies, the residual disk component enhances the in–plane susceptibilities, χxx ≃χyy > χzz,(217) producing an oblate coherence tensor and a flattened transport envelope, though still without the strong phase locking or large eigenvalue contrast characteristic of spiral galaxies. Triaxial ellipticals exhibit three distinct susceptibilities, each aligned with a principal inertial axis, naturally yielding the observed E3–E7 sequence as a continuous family of anisotropic transport–coherence tensors. Thus the projected halo ellipticity traces the eigenvalue hierarchy of χij(Ω). Irregular galaxies represent the partially coherent limit. Low baryonic densities and turbulent, stochastic environments prevent long–lived phase locking, yet short–range curvature correlations arise within transient stellar associations. Each system is therefore a fluctuating ensemble of small coherence patches whose cross–correlations are weak but nonzero, 0<⟨AiA∗ j⟩ ⟨|A|2⟩<1,(218) continually forming, dissolving, and merging without establishing a global mode. This picture agrees with their patchy morphology, irregular kinematics, and diffuse halos: irregulars are turbulence–driven curvature ensembles that occupy the transition between the coherent and fully random extremes of the resonant– speckle hierarchy. 8.3 Clusters as Radiative Transport Extensions of Ellipticals 8.3.1 Observing the Statistical Transport Equilibrium In this section we extend the radiative transport formalism developed for isolated ellipticals and spirals to the cluster scale. Whereas the elliptical case required a transport treatment to capture long-lived diffusive tails, and the spiral case was fully described by coherent Foldy–Lax scattering, the cluster environment involves the superposition of many independent elliptical transport envelopes. These envelopes overlap and partially interfere, producing a collective intensity field that must be modelled statistically. The primary goal of this section is to demonstrate that the resulting effective transport density reproduces, in both form and amplitude, the observed darkmatter distribution of clusters as described by the Navarro–Frenk–White (NFW) profile [21, 22]. The NFW model, derived empirically from cosmological N -body simulations, describes the dark-matter halo as ρNFW(r) = ρs x(1+x)2, x =r rs ,(219) 83 with rs and ρs as characteristic scale radius and density, respectively. Recent extensions include the introduction of a “splashback” radius Rsp , beyond which ρ(r) rapidly steepens [23]. The theoretical basis for the cluster-scale treatment presented in this section, lies in the long-standing radiative transport framework developed for diffuse acoustic and optical propagation. Weaver [ 13 , 14 ] established that acoustic energy in heterogeneous solids obeys a radiative transfer equation for the intensity field, with distinct ballistic, coherent, and diffusive regimes. These works provided the first quantitative formulation of acoustic diffusion and localization, defining the statistical kernel that our own model generalises to astrophysical scales. The subsequent experimental validation was achieved by Derode, Tourin, and Fink [ 9 , 4 , 5 ], who investigated ultrasonic propagation through random ensembles of resonant scatterers. Their experiments confirmed the radiative transport regime predicted by Weaver, demonstrating the separation of coherent, ballistic, and diffuse components, as well as the emergence of reverberation gain in multiply scattered fields. These resonant-rod systems represent direct analogues of our elliptical galaxies: individual mesoscopic resonators embedded in a diffusive host that collectively generate extended energy-density tails. Later developments by Weaver and Lobkis [ 16 , 17 ] established that diffuse fields from independent random ensembles can exhibit weak cross-coherence, recoverable through cross-correlation measurements. This phenomenon—now widely used in seismology and acoustics—formalised the concept of overlapping diffuse transport fields, providing a rigorous physical basis for the small nonlinear coupling parameter ηused in our cluster model. Finally, Larose et al. [ 18 ] unified these concepts across acoustics, optics, and geophysics, showing that the superposition of diffuse fields from independent sources can reproduce effective Green’s functions through reverberant correlations. This provides the closest experimental and theoretical precedent for our treatment of overlapping elliptical transport envelopes within clusters: random, extended resonators (galaxies) embedded in a diffusive medium (intracluster baryons) that generate partially coherent, overlapping transport tails and a measurable ensemble reverberation gain. The ensemble behaviour demonstrated in the acoustic and optical transport experiments of Weaver and Derode provides the direct methodological foundation for the present cluster model. Here we apply the same radiative–transport formalism—previously shown to describe overlapping diffuse fields in random media—to the collective curvature–transport regime of spacetime itself. Each elliptical galaxy sustains its own diffusive envelope, generated by the random stellar motions that define its internal transport equilibrium, and these envelopes coexist within a shared effective curvature field (the common transport continuum of spacetime) that supports their mutual overlap. The resulting ensemble of overlapping transport intensities forms a statistically smooth superposition whose 84 weak cross–coherence produces the modest reverberation gain Grev =1+Ngalη, (220) which, when applied at the cluster scale, reproduces the NFW–like curvature of the effective density profile without invoking any additional physical assumptions. In the sections that follow, we adopt identical methodologies to those used in the acoustic and optical transport studies, applying them directly to the astrophysical transport regime to construct and normalize the cluster-scale effective density field. Each elliptical or lenticular acts as amesoscale driver of the effective curvature field, ρeff,i(r)∝exp(−r/Lmean),(221) where Lmean =ceffτtr (222) is the mean transport distance determined by the effective transport time τtr . Within a cluster, these envelopes overlap and interfere weakly. Their mutual correlation defines a small but finite reverberation parameter η, ⟨AiA∗ j⟩=ηq⟨|Ai|2⟩⟨|Aj|2⟩,(i=j),(223) yielding the statistical amplification Grev =1+Ngalη, (224) where Ngal is the number of contributing resonators. This weak collective gain reproduces the macroscopic curvature excess associated with the observed cluster halo. Theorem (Statistical Transport Theorem for Galaxy Clusters). Let a cluster contain Ngal elliptical or lenticular galaxies, each producing a stationary diffusive envelope ρeff,i ( r ) ∝exp ( −r/Lmean )within a shared effective transport field of speed ceff and damping rate Γ tr . If the ensemble satisfies the statistical balance condition dρeff dt =G −L ≃ 0,(225) then the ensemble–averaged effective density obeys ρeff(r)∝Grev exp[−2r/Lmean] r2,(226) with Lmean =ceff Γtr , Grev =1+Ngalη. (227) The convolution of these overlapping exponential envelopes approaches the observed NFW form ρ(r)∝r−1(1+r/rs)−2for Lmean ∼rs. 85 Interpretation. Equation (244) defines the steady transport equilibrium: the effective curvature energy density, governed by transport, is continually replenished by embedded elliptical and lenticular galaxies. The relevant length scale Lmean depends on the transport properties of the medium, not on the baryonic density or resonant frequency. Ellipticals act as isotropic diffusive emitters; lenticulars, with their oblate morphology, generate aspherical transport tails preferentially elongated toward the cluster centre—a behaviour consistent with X–ray and weak–lensing observations indicating that S0s are radially aligned and concentrated toward cluster cores. Spiral galaxies can contribute negligibly to the cluster effective curvature field because they have no transport tails and are correspondingly observed to populate only the periphery. The same transport equations extend naturally to filaments, which represent the cylindrical continuation of the cluster ensemble. In this geometry, overlapping diffusive envelopes align along a filamentary axis, and the mean gain is reduced by geometric dilution: Gfil ≃Grev Lmean Rfil ,(228) where Rfil is the filament radius. Although weaker in amplitude, these linear transport ensembles preserve the same physical mechanism—curvature energy diffusion and reverberation coupling among extended resonators—and thus naturally explain the emergence and persistence of cosmic filaments. Finally, the effective curvature field that sustains the cluster halo is self–maintaining in the same sense as in smaller systems. The apparent dark mass produced by the collective amplification Grev is also the mechanism that enhances the coupling among member galaxies: the increased curvature energy density focuses transport back toward dense regions, strengthening overlap and maintaining nonlinear equilibrium. The cluster is therefore a self–regulated, radiative–transport system—a statistically sustained effective curvature field whose ensemble dynamics reproduce both the amplitude and morphology of observed dark–matter halos. 8.3.2 Mathematical Formulation: Overlapping Elliptical Transport Envelopes Each elliptical galaxy acts as an effective source of the curvature–phonon field amplitude Ai (r), whose ensemble–averaged intensity ⟨|Ai|2⟩ is given by the single–galaxy transport solution derived in Eq. (195) . At the cluster scale, the total field amplitude is the coherent sum Atot(r) = Ngal X i=1 Ai(r),(229) so that the ensemble–averaged intensity becomes |Atot(r)|2=X i⟨|Ai|2⟩+X i=j⟨AiA∗ j⟩.(230) 86 The cross–correlation term ⟨AiA∗ j⟩ represents the weak mutual coherence between overlapping transport tails. Following the standard radiative–transport treatment of diffuse field superposition [13, 9, 4, 16, 18], we write ⟨AiA∗ j⟩=ηijq⟨|Ai|2⟩⟨|Aj|2⟩,(231) where 0 ≤ηij ≪ 1 measures the fractional coherence of the pair. For statistically independent galaxies, ηij fluctuates randomly about a small mean value η . The total intensity therefore reads |Atot|2≃X i⟨|Ai|2⟩(1+Ngalη),(232) identifying the dimensionless reverberation gain Grev =1+Ngalη. (233) This amplitude–level formulation maintains the same causal hierarchy as the elliptical case: the field intensity arises from the quadratic combination of statistically independent amplitude realizations, and the weak coherence factor η provides a small multiplicative enhancement due to overlapping transport regions. Analogous ensemble treatments appear in radiative–acoustic transport and time–reversal acoustics [ 13 , 9 , 5 , 17 , 18 ], where the reverberant energy of independent diffuse sources produces a measurable intensity gain. The corresponding effective energy–density field for the cluster follows directly from the ensemble intensity law of Eq. (195) . Each elliptical contributes a transport–renormalized amplitude envelope ⟨|Ai|2⟩ ∝ Keff ( | r − r i| ), including the gradient–energy boost factor Ξ ∇ defined in Eq. (204) . The cluster-scale superposition is therefore ρeff(r) = Grev Ξ∇ Ngal X i=1 Keff(|r−ri|),(234) where Grev accounts for the weak reverberant coupling between overlapping transport tails. This discrete form will subsequently be replaced by its continuous analogue through convolution with the galaxy number–density profile ngal ( r ) in the following subsection. This formulation directly parallels the intensity–based radiative transport models developed for multiply scattering acoustic and elastic media. Weaver [ 13 , 14 ] first expressed the diffuse energy field as an ensemble average over incoherent amplitude realizations, while Derode, Tourin, and Fink demonstrated experimentally that weak mutual coherence among independent resonant clusters produces a measurable reverberation gain [ 9 , 4 , 5 ]. Subsequent correlation analyses by Weaver and Lobkis [ 16 , 17 ] and by Larose [ 18 ] showed that such diffuse fields retain small but finite cross–coherence, recoverable through cross–correlation of 87 independent realizations. Our treatment of the cluster intensity follows this same logic: the total field is the incoherent sum of independent elliptical transport envelopes, augmented by a small global coherence factor η representing the ensemble–averaged reverberant coupling between them. In principle, one might attempt to model the total cluster field directly from the ensemble of overlapping galaxy amplitudes, fully accounting for all phase correlations. In practice, however, this approach is both computationally prohibitive and conceptually fragile, as the individual phases are rapidly randomized by multiple scattering and transport effects. We therefore proceed by employing the standard practice of radiative transport theory to model phenomenon of this type: by constructing the effective intensity (or energy–density) field instead, retaining the correct ensemble statistics while averaging over the random phases. This yields a tractable and physically transparent formulation, which naturally leads to the stepwise development outlined below. Continuous Galaxy Distribution and NFW Analogy Observational studies indicate that the radial number–density profiles of cluster galaxies are well described by an NFW–like form, albeit with systematically lower concentrations than the corresponding dark–matter halos [ 24 ]. For analytical convenience and physical realism, we therefore model the spatial distribution of member galaxies as an NFW–type profile smoothly truncated at large radii by a “splashback” [23] term: ngal(r)=n0 Wsp(r) x(1+x)2, Wsp(r) = 1 1 + exp[(r−Rsp)/∆sp],(235) where x = r/rs , and Rsp and ∆ sp represent the splashback radius and its characteristic transition width, respectively. This continuous parameterisation replaces the ∼ 100–500 discrete galaxies of a typical rich cluster with a smooth number–density field appropriate for convolution with the single–galaxy transport kernel. A key interpretative point is that these galaxies are not static particles. Their long diffusive decay times ensure that, over many dynamical periods, the transport envelopes of individual members overlap throughout the cluster volume, forming an effective continuum of resonant sources. The Necessity of 2.5D Convolution A direct three–dimensional convolution of the form ρeff(r) = ZK(|r−r′|)ngal(r′)d3r′(236) does not yield physically consistent results because K ( r ) describes an angle–integrated intensity rather than a Green’s function. This distinction is fundamental: a Green’s function preserves full phase and directional information, whereas an intensity kernel embodies an ensemble average over those angles. The 88 issue is well known in transport theory and has been extensively discussed in the contexts of optical diffusion, ultrasound propagation, and acoustic multiple scattering [8, 4, 5, 13, 19]. The appropriate formalism—long established in both radiative and acoustic transport—is to employ an angularly averaged, reduced–dimensional representation: ρeff(r)=2πZ∞ 0 ngal(r′)r′⟨K(s)⟩µdr′, s2=r2+r′2−2rr′µ, (237) where ⟨K ( s ) ⟩µ denotes the Legendre–weighted average over µ = cos θ . This “2.5D” formulation removes redundant angular degrees of freedom and guarantees mathematical stability while preserving the correct physical form of the transport field. It omits only a single geometric length factor, recoverable through explicit normalization. Such dimensional reductions are standard practice in radiative transfer [ 8 ], ultrasound diffusion [ 4 , 5 , 13 ], and optical multiple scattering [ 19 ]. In all of these contexts, the reduced formalism omits a single geometric length factor, producing an apparent scale ambiguity that is conventionally resolved by introducing a macroscopic normalization parameter, Λ geom , determined through global energy or mass conservation. In our case, because the convolution represents a direct superposition of intensities rather than fields, the appropriate normalization follows naturally from the total enclosed mass of the constituent galaxies. Writing ρeff(r)=Λgeom ρshape(r), we impose Z∞ 0 4πr2ρeff (r)dr =Ngal M(elliptical) DM ,(238) which yields the explicit expression Λgeom =Ngal M(elliptical) DM Z∞ 0 4πr2ρshape(r)dr .(239) This ensures that the integrated effective transport density recovers the combined dark–mass contribution of all elliptical members, anchoring the 2.5D convolution to a physically measurable mass scale without introducing additional free parameters. Reverberation Gain and Partial Coherence The parameter η quantifies weak statistical correlations among overlapping transport tails, representing the residual coherence that remains after ensemble averaging. Such partial coherence is a well–established feature of radiative and acoustic transport systems containing multiple scattering and long–lived resonances. In ultrasonic and optical analogues [ 4 , 5 , 13 , 19 , 17 ], the effect arises from recurrent scattering loops and time–delayed “reverberant” returns that couple neighbouring intensity fields. These processes lead to a small, positive excess in the ensemble–averaged intensity relative to the purely incoherent sum. 89 Because many individual transport tails overlap within a few Lmean, the cluster behaves as a continuous curvature field whose radial symmetry emerges naturally from the convolution of exponentially decaying envelopes. The NFW-like slope of the effective density profile reflects the statistical superposition of these diffusive fields, while its amplitude is set by the reverberation gain Grev. Lenticular (S0) galaxies and transitional systems occupy an intermediate position in the transport hierarchy. Their flattened stellar distributions generate anisotropic transport envelopes, though without the global phase locking of spiral disks. Within clusters, S0s are more centrally concentrated and frequently exhibit weak radial or cluster–axis alignments, reflecting the morphology–density relation and the influence of the cluster potential. These factors naturally bias their transport tails toward the cluster interior, enhancing reverberation in the core without requiring any specific disk–orientation mechanism. Spiral galaxies, by contrast, contribute negligibly to the cluster–scale field since their coherent m= 1 modes do not sustain extended transport tails. The same radiative–transport principles extend to cosmic filaments, which may be regarded as cylindrical continuations of the cluster ensemble. Here, resonators (galaxies, groups) lie along an extended axis, and their diffusive envelopes overlap primarily in one dimension. The reduced symmetry decreases the effective gain but preserves the transport mechanism, producing elongated curvature envelopes aligned with the filament spine. In summary, a galaxy cluster behaves as a reverberant curvature–transport system: a self-averaging ensemble of overlapping diffusive envelopes whose collective behaviour reproduces the NFW form of the apparent dark-matter halo. Curvature energy is sustained statistically rather than coherently, maintained through continual diffusive exchange among the constituent galaxies. Spirals correspond to coherent resonant emitters, ellipticals to self-sustaining diffusive cavities, lenticulars to anisotropic transitional resonators, clusters to radiative–transport continua, and filaments to linear transport ensembles. Together these systems form a continuous hierarchy of curvature–matter coupling, governed throughout by the same causal transport framework. 9 Local Silence and the Preservation of Classical Gravity The two–channel framework developed in Section 4.1 and Section 6.2 immediately explains why all laboratory, Solar–System, and stellar tests of gravity recover the Newton–Einstein law with complete fidelity. Channel 2—the Newtonian Channel—represents the quasi–static stress–response loop and is identical to classical gravity. It determines the slowly varying gravitational potential ∇2ΦN= 4πG (ρb,0+ρeff),(246) 96 where ρb,0 is the coarse–grained baryonic density and ρeff is the cycle–averaged curvature energy stored over cosmological timescales by Channel 1. Since Channel 2 governs all observable accelerations, it reproduces Newtonian and relativistic gravity in every regime accessible to experiment. Channel 1, the Resonant Channel, is entirely different in character. It is driven solely by the coherent component ρb,1 of the baryonic motion, producing an oscillatory curvature perturbation Φ res . The response is given by the causal susceptibility of the effective field, Φres =GΩ∗ρb,1,(247) but this oscillatory field does not enter the directly observable gravitational potential. Instead, its cycle–averaged curvature energy is transferred extremely slowly into the stored component ρeff over the memory time τmem = Γ−1,(248) which exceeds stellar and galactic dynamical times by many orders of magnitude. For all practical and observational purposes, ρeff is therefore constant. The decomposition ρb=ρb,0+ρb,1,Φ = ΦN+ Φres,(249) makes clear that only Φ N governs geodesics. Channel 1 produces rapidly oscillatory curvature, but its only macroscopic imprint is the ultra–slow accumulation of stored curvature energy, ρeff ∝DGΩ∗ρb,1 2Eτmem ,(250) which behaves observationally as a time–independent dark–mass density. This separation of timescales defines the regime of local silence. No laboratory, Solar–System, or stellar experiment can detect Φ res directly, because any change in ρeff unfolds only over the cosmological memory time (248) . Even a dramatic event—such as instantaneously removing every elliptical galaxy from a cluster— would not erase the cluster’s dark envelope: the stored curvature energy would decay only over many billions of years. To any local observer, the gravitational field would appear unchanged, continuing to follow the Newton–Einstein potential ΦN. Hence Channel 1 is not merely dynamically quiet: it is observationally silent. Its influence enters measurable dynamics only through the quasi–static stored component ρeff that feeds Channel 2. This ensures full consistency with all precision tests of Newtonian and relativistic gravity, while permitting resonant curvature to generate effective dark–mass halos only in systems where coherent baryonic motion persists over cosmic times. 97 A useful analogy is supplied by acoustics. The acoustic field oscillates rapidly relative to the slow buildup of time–averaged pressure; what one measures is not the sound field itself but the quasi–static pressure it leaves behind. Similarly, the resonant curvature field Φ res is not observable in isolation. Only its stored cycle– averaged energy density ρeff , acting gravitationally through Channel 2, enters measurable dynamics. This captures precisely why the resonant channel remains locally silent even while shaping the large–scale gravitational environment. 10 Prediction: Absence of Gain in Homogeneous Systems This section establishes a central and falsifiable prediction of the two–channel causal–resonant framework: a perfectly homogeneous baryonic medium cannot generate curvature gain in Channel 1 under any circumstances. Inhomogeneity is required to form the source–rescatter loops, Foldy–Lax coupling, and reverberation processes that define the effective dark density ρeff . A spatially smooth system excites only Channel 2 (the Newtonian Channel) and is therefore resonantly silent. To avoid confusion, we first clarify how the linear perturbation theory of Section 4.4 relates to homogeneous media. The perturbative equation, repeated here for reference, ∇2δΦ + 1 c2 eff −∂2 t−2Γ ∂t+ Ω2δΦ = 4πG c2 eff δρb,(251) presupposes that multiple scattering has already produced an effective medium with renormalised parameters Ω, Γ, and ceff . A homogeneous configuration cannot generate such a medium and therefore cannot satisfy Eq. (251) in isolation. 10.1 Three Distinct Regimes for Homogeneous Clouds A homogeneous cloud may arise in three physically distinct settings: (1) Isolated homogeneous cloud. The cloud possesses a natural frequency ωnat ∝√ρ0 , but generates no scattering, no self–energy, and therefore no renormalised parameters (Ω,Γ, ceff). Thus Channel 1 does not form: Φ(isolated) res = 0.(252) (2) Cloud embedded within an externally generated effective medium but not phase–locked. The surrounding medium possesses (Ω , Γ , ceff ), but the cloud oscillates at its own ωnat and cannot couple coherently through the causal Green function. Hence no Channel 1 field is generated: Φ(embedded,free) res = 0.(253) 98 (3) Cloud embedded within an externally generated effective medium and forced at the global resonance Ω.In this artificial situation the homogeneous cloud may oscillate in temporal phase with the background because the driving field imposes the phase externally. The cloud then satisfies the driven Helmholtz relation ∇2+k2 ΩδΦ = 4πG c2 eff δρb, k2 Ω=Ω2 c2 eff ,(254) but the scattering strength remains σ(x)≡0,(255) because temporal coherence does not create spatial contrast. Thus no Foldy– Lax loops form, no self-energy is generated, and the response remains a strictly single–pass, non-amplifying Helmholtz field. The theorem below applies to all three regimes. 10.2 Theorem: Homogeneous Media Produce No Resonant Gain Theorem (Absence of Resonant Gain in Homogeneous Systems). Let ρb (x) = ρ0 be a perfectly homogeneous baryonic distribution with no spatial contrast. Then, in each of the regimes in Section 10.1, the Channel 1 curvature response exhibits no resonant amplification: 1. Isolated case: No effective medium forms, so Channel 1 does not exist: Φres = 0.(256) 2. Embedded but free-running: Without phase–locking, the cloud does not couple to the resonant field: Φres = 0.(257) 3. Embedded and phase–locked: Even when driven at the global resonance Ω, the response is single–pass and non–amplifying because σ (x) ≡ 0: A=GΩ∗(σA) =⇒A= 0.(258) The reverberation gain vanishes, Grev =1+Nη −→ 1, η = 0,(259) and the stored curvature energy remains at the single–pass level: ρ(hom) eff ∼1 Nρ(inhom) eff ,(260) for systems with the same total mass, where N counts the discrete substructures available in the inhomogeneous case. Therefore a perfectly homogeneous baryonic medium cannot generate measurable curvature amplification in Channel 1 under any physical configuration. 99 10.3 Interpretation and Connection to Perturbation Theory The theorem shows that spatial contrast is essential for Channel 1. Without it, no Foldy–Lax loops form, no self–energy accumulates, and the collective resonant mode does not arise. This explains why the perturbative framework of Section 4.4 cannot be applied to a homogeneous system in isolation: the parameters (Ω , Γ , ceff ) exist only after rescattering has generated an effective medium. In the idealised circumstance of regime (3), where a homogeneous cloud is externally forced to oscillate at the global resonance within an already-existing effective medium, the response reduces to the driven Helmholtz convolution δΦ(x) = ZGΩ(|x−x′|)δρb(x′)d3x′,(261) but its amplitude remains faint, being limited to the single-pass level because no spatial contrast exists to sustain coherent feedback. The causal–resonant model therefore predicts that apparent dark structure correlates with baryonic contrast and temporal coherence, not with total baryonic mass: ρeff ∝(baryonic contrast) ×(coherence time).(262) In cluster collisions such as the Bullet Cluster (1E 0657 − 56), where dense, collisionless galaxies retain their spatial contrast while the diffuse plasma is stripped away, the curvature resonance must remain anchored to the galaxy component rather than to the displaced gas. This behaviour precisely matches observations without invoking non–baryonic matter, and yields a decisive test: if lensing peaks consistently coincide with regions of greatest baryonic contrast, the causal– resonant interpretation gains strong empirical support. 11 Emergence of Multiple Scattering from a Weakly Inhomogeneous Continuous Medium 11.1 A Brief Review of Garvitational Multiple Scatter Theory Before attempting to analyse the early Universe, it is useful to summarise the assumptions that underlie the effective–medium framework developed thus far. Throughout previous sections we have considered a baryonic mass distribution ρb(x) = N X j=1 mjδ(x−xj),(263) 100 consisting of spatially separated stellar or galactic constituents. Within a coherent domain these constituents oscillate at a common curvature frequency Ω2∝G¯ρb,(264) where ¯ρbis the root–mean baryonic density over the coherent volume. Each baryonic element responds to the resonant curvature field as a damped oscillator, ¨ Xj+ 2Γ ˙ Xj+ Ω2Xj=αΦres(xj),(265) which reinjects curvature waves back into the medium. Substituting the induced oscillators into the Green–function representation yields the self–consistent Foldy– Lax equation, Φres(x)=Φinc(x) + ZGΩ(|x−x′|)σ(x′) Φres(x′)d3x′,(266) where the contrast–defined scattering strength is σ(x)∝δρb(x).(267) Applying the Foldy–Lax operator to both sides of Eq. (266) produces the renormalised effective–medium wave equation, ∇2Φres +1 c2 eff −∂2 t−2Γ ∂t+ Ω2Φres =4πG c2 eff ρb,1,(268) where ρb,1 denotes the slowly varying component of the baryonic distribution that drives the coherent mode. All renormalised parameters (Ω , Γ , ceff ) arise self– consistently from repeated rescattering encoded in Eq. (266) . The cycle–averaged stored curvature energy defines the effective density, ρeff(x)∝ ⟨|Φres(x)|2⟩.(269) Crucially, nothing in Eqs. (264) – (269) requires ρb (x) to be discrete. Any baryonic distribution with sufficient spatial contrast—even a continuous density with internal structure— produces a nonzero σ (x) and therefore participates in the multiple–scattering mechanism. This observation is essential for cosmology: a continuous medium can become a resonant, scattering effective medium once it develops spatial inhomogeneity. However, the early Universe does not begin in that regime. The primordial baryon–photon plasma is extremely smooth, with δρb ρb∼10−5,(270) and therefore realises Regime 1 of Section 10.2, namely an isolated, homogeneous system with σ(x)≡0, Grev = 1,Φres = 0.(271) 101 In this state no Foldy–Lax feedback operates, no effective self–energy can form, and the perturbation theory of Section 4.4 does not apply. There is therefore ρeff(x, t) = 0 (early homogeneous era).(272) Only when gravitational instability amplifies density fluctuations to δ∼ O (1) does σ (x) become nonzero and the Foldy–Lax transport mechanism turn on. The remainder of this section traces this transition: how a continuous, nearly homogeneous plasma dynamically evolves into the resonant, multiply–scattering environment required for the effective dark density described in the rest of this work. Although σ = 0 eliminates the multiple–scattering channel and therefore suppresses any renormalised feedback, it does not eliminate the underlying homogeneous curvature resonance. Even in a perfectly uniform baryon–photon plasma the bare oscillation ¨ Φres + Ω 2 0 ( t ) Φ res = 4 πG ¯ρb ( t ) δb remains active as a coherent, spatially homogeneous mode, exactly like the Langmuir oscillation in a uniform electron plasma. In this regime the cycle–averaged energy of the mode need not vanish, despite the absence of a scattering medium. The central assumption of this work is that baryonic matter generically supports gravitational oscillation modes with natural frequency Ω 2∝G¯ρ , mirroring the familiar classical g–mode relation in stellar interiors. At late times such modes are driven by internal heat and pressure imbalances within stars, while in the early Universe they are excited by thermal and pressure fluctuations in the tightly coupled baryon–photon plasma. Thus, although the physical driver varies between epochs, the same gravitational restoring resonance is present throughout: baryonic matter at any density supports an oscillatory gravitational response, and any process capable of producing coherent compressions—whether stellar heat at late times or photon–baryon pressure at early times—can excite the same underlying curvature mode. 11.2 The Resonant Frequency in a Continuous Medium Before an effective medium forms, the local (bare) curvature resonance is determined solely by the coarse–grained baryon density: Ω2(x, t)∝G¯ρb(x, t),(273) In the early plasma, the coarse–grained density takes the form ¯ρb(x, t) = ¯ρb(t) [1 + δ(x, t)],|δ| ≪ 1,(274) so that Ω(x, t) = Ω0(t)1 + δ(x, t) 2+O(δ2).(275) 102 Even tiny density perturbations therefore imprint a weak spatial modulation on the local curvature resonance. However, as the homogeneity theorem states, such fluctuations do not produce a nonzero scattering strength: σ(x)∝δρb(x)⇒σ≡0 (δ∼10−5).(276) Thus the early medium lies strictly in the non-scattering regime of Regime 1. 11.3 Weak Random Potential and the Ballistic Regime Formally inserting (275) into the resonant operator yields a Helmholtz equation of the form ∇2+k2 0+V(x, t)Φres =S, k0=Ω0(t) ceff ,(277) with weak random potential V(x, t)=2k0 δΩ(x, t) ceff +O(δ2)∝δ(x, t).(278) For a continuum with RMS amplitude σV and correlation length ℓc , random– medium theory gives the single–scattering mean free path ℓ−1 s(t)∝k4 0(t)σ2 V(t)ℓ3 c(t).(279) In the primordial plasma one has σV∼δ∼10−5, giving k0ℓs≫1,(280) so curvature propagation is entirely ballistic. There is no rescattering and hence no possibility of forming an effective medium. This resolves the apparent tension with Section 4.4: since no effective medium exists, the perturbative theory there does not apply at early times. 11.4 Nonlinear Growth Toward the Transport Threshold As gravitational instability amplifies density fluctuations, δ→ O(1), σV→ O(1), ℓc→galactic scales,(281) the combination in Eq. (279) grows by many orders of magnitude. There exists a cosmological epoch t∗defined implicitly by k0(t∗)ℓs(t∗)∼1,(282) which is the analogue of the Ioffe–Regel criterion for waves in random media. At t∗ curvature–phonon propagation transitions from the ballistic to the multiple– scattering regime. Beyond this point, the density field is no longer a weak perturbation, and Channel 1 becomes dynamically activated. 103 11.5 Emergent Resonant Patches and Effective Discreteness When δ∼ 1, the medium fragments into nonlinear overdensities. Within each patch, Ω(x)≃Ωpatch ∝√ρpatch,(283) so curvature energy becomes locally trapped. The cycle–averaged effective density grows as ρeff(x)∝|Φres|2,(284) leading to positive feedback: curvature trapping enhances mass inflow via Channel 2, sharpening contrast until each region behaves as a discrete, mesoscale resonant scatterer. These “resonant patches” constitute the progenitors of ellipticals and the basic units required for the Foldy–Lax and transport formalisms of Sections 5–8.3.4. This mechanism closely parallels modulational instability in plasmas, where spatial variations in ωp ( x ) ∝pne(x) produce a refractive index landscape that traps Langmuir waves. Here, fluctuations in Ω(x) trap curvature–phonons, giving rise to long–lived resonant envelopes that function as dark–matter analogues within the present theory. The predicted cosmological progression is therefore homogeneous plasma (Regime 1) −→ weak random medium (ballistic) −→ nonlinear resonant patches −→ multiple–scattering resonators (elliptical progenitors) −→ coherent disks (spirals) −→ transport ensembles (clusters, filaments). The onset of multiple scattering is not assumed but rather emerges dynamically when the mean free path falls to the resonant scale, enabling the same effective– medium physics that governs the present-day Universe. 11.6 Linear Limit and Connection to the CMB Power Spectrum The developments of Sections 5–8.3.4 demonstrate that multiple scattering, renormalised propagation, and ρeff arise only once the contrast σ (x) becomes nonzero. In the primordial plasma, however, Eqs. (270) – (272) imply that the Universe lies strictly in Regime 1, so that σ(x) = 0,Φres = 0,(Ω,Γ, ceff) undefined,(285) and no effective medium exists. Despite the absence of a scattering medium, the underlying bare oscillator— the curvature–phonon degree of freedom with frequency Ω 0 ( t ) from Eq. (273) — remains a well-defined linear response field. In Regime 1 all renormalised parameters revert to their bare values: ceff →c, Γ→0,Ω→Ω0(t),(286) 104 and the effective–medium PDE (268) consistently reduces to the ballistic Helmholtz equation ∇2+1 c2−∂2 t+ Ω2 0(t)Φres(x, t) = 4πG c2¯ρb(t)δb(x, t),(287) which is simply the linear curvature response to the baryon overdensity. Fourier transforming Eq. (287) gives ¨ Φres(k, t) + Ω2 0(t) Φres(k, t)=4πG ¯ρb(t)δb(k, t),(288) showing that Φ res isaforced, linear oscillator sourced by the baryonic overdensity. Importantly, because σ = 0, there is no self-energy, no renormalisation, and no feedback loop: Eq. (288) adds no nonlinear gravitational effect. It is important to emphasise that Eq. (288) does not introduce an additional gravitational potential beyond the usual Newtonian metric perturbations (Φ N, Ψ N ). Rather, Φ res is a linear curvature response field driven by the same baryon fluctuations that also generate the Newtonian potentials. There is no independent source term and no modification of the Einstein equations in Regime 1. Consequently, the combination ΦN(k, t)−→ ΦN(k, t) + ϵ(t) Φres(k, t) (289) should be viewed merely as a bookkeeping device for how curvature–phonon oscillations produce a tiny correction to the effective driving of the photon Boltzmann hierarchy. It does not represent new gravitational physics or an additional source term. The dimensionless prefactor ϵ(t)≡Ω2 0(t) k2c2(290) quantifies the relative size of this correction. It is possible to perform a basic parametric estimate at recombination. Using ¯ρb(z∼1100) ∼10−21 kg m−3gives Ω0(t∗)∼pG¯ρb∼10−17 s−1.(291) The comoving acoustic scale corresponds to k∼0.01–0.1 Mpc−1, i.e. k c ∼10−13–10−12 s−1.(292) Thus ϵ(t∗) = Ω2 0 k2c2∼10−8–10−10,(293) many orders of magnitude below unity. Therefore: ϵΦres ≪ΦN,(recombination era),(294) 105 41 # ================================================================ 42 43 import numpy as np 44 import mpmath as mp 45 import matplotlib.pyplot as plt 46 47 # ---------- Global plot style ---------- 48 plt.rcParams.update({ 49 ’font.size’: 16, 50 ’axes.titlesize’: 18, 51 ’axes.labelsize’: 16, 52 ’xtick.labelsize’: 14, 53 ’ytick.labelsize’: 14, 54 ’legend.fontsize’: 14 55 }) 56 57 # ---------- Physical constants ---------- 58 G = 6.67430e-11 59 Msun = 1.98847e30 60 pc = 3.085677581491367e16 61 kpc = 1.0e3 * pc 62 c_light = 2.99792458e8 63 64 # ================================================================ 65 # USER-TUNABLE MODEL PARAMETERS 66 # ------------------------------------------------ 67 # These are the only physically meaningful dials: 68 # Omega (fixed by v0, r0), Q, X, and eps. 69 # ================================================================ 70 v0 = 230e3 # flat rotation speed [m/s] 71 r0 = 8.5 * kpc # reference radius [m] 72 Omega = v0 / r0 # angular frequency [s^-1] 73 N_star = 6.0e10 # number of baryonic emitters (stars) 74 75 c_eff = 2.0e5 # effective wave speed [m/s] 76 Q = 12.6 # quality factor (temporal coherence) 77 X = 1 # mean coupling efficiency (|X| 1) 78 G_env_root = np.sqrt(X*X*N_star) # amplitude-level environmental gain 79 Gamma = Omega / (2 * Q) # damping rate 80 eps = 0.3 # geometric m=1 projection factor 81 # ================================================================ 82 83 # ---------- Derived scales ---------- 84 k_complex = (Omega + 1j * Gamma) / c_eff 85 lambda_res = 2 * np.pi * c_eff / Omega 86 L_damp = c_eff / Gamma 87 L_z = 0.3 * L_damp # vertical confinement scale << L_damp 88 89 print(f"Omega= {Omega}, L_damp = {L_damp/kpc:.2f} kpc, _res = ,→{lambda_res/kpc:.2f} kpc, L_z = {L_z/kpc:.2f} kpc") 90 91 # ---------- Exponential baryonic disk ---------- 92 M_disk_total = N_star * Msun 93 R_d = 3.0 * kpc 112 94 h_z = 0.30 * kpc 95 Sigma0 = M_disk_total / (2.0 * np.pi * R_d**2) 96 97 def Sigma_r(r): return Sigma0 * np.exp(-r / R_d) 98 def rho_b_midplane(r): return Sigma_r(r) / (2.0 * h_z) 99 100 # ---------- Radial grid ---------- 101 r_max = 1000.0 * kpc 102 N = 1800 103 r = np.linspace(1.0e-3 * kpc, r_max, N) 104 105 # ---------- Special functions ---------- 106 def J1(z): return mp.besselj(1, z) 107 def H1(z): return mp.hankel1(1, z) 108 109 J1_vals = np.array([complex(J1(k_complex * ri)) for ri in r]) 110 H1_vals = np.array([complex(H1(k_complex * ri)) for ri in r]) 111 112 # ---------- Source moment ---------- 113 M1_vol = 2.0 * np.pi * eps * rho_b_midplane(r) 114 115 # ---------- Complex trapezoidal integration ---------- 116 def cumtrapz_complex(y, x): 117 out = np.zeros_like(y, dtype=complex) 118 for iin range(1, len(x)): 119 out[i] = out[i-1] + 0.5 * (y[i] + y[i-1]) * (x[i] - x[i-1]) 120 return out 121 122 integrand_J = r * M1_vol * J1_vals 123 integrand_H = r * M1_vol * H1_vals 124 125 I_J_cum = cumtrapz_complex(integrand_J, r) 126 I_H_cum = cumtrapz_complex(integrand_H, r) 127 I_H_total = I_H_cum[-1] 128 129 # ---------- Envelope field A(r) ---------- 130 prefac = 4.0 * np.pi**2 * 1j * G / (c_eff**2) 131 A_r = prefac * (J1_vals * (I_H_total - I_H_cum) + H1_vals * I_J_cum) 132 133 # ---------- Apply coherence gain ---------- 134 # Environmental amplification (Q (X N_star)) 135 gain_field = Q * G_env_root 136 A_r *= gain_field 137 138 # ---------- Effective and baryonic densities ---------- 139 # Field energy density |A|, already depth-normalized. 140 rho_eff_mid = (2*Omega**2 / (8.0 * np.pi * G * c_light**2)) * np.abs(A_r)**2 141 rho_b = rho_b_midplane(r) 142 143 # ---------- Vertically averaged surface densities ---------- 144 # The vertical column depth is represented by L_damp (no extra L_z ,→normalization) 145 Sigma_eff = L_damp * rho_eff_mid 146 113 147 # ---------- Enclosed masses ---------- 148 M_b, M_eff = np.zeros_like(r), np.zeros_like(r) 149 for iin range(1, len(r)): 150 dM_b = 2 * np.pi * 0.5 * (r[i]*rho_b[i] + r[i-1]*rho_b[i-1]) * ,→(r[i]-r[i-1]) * (2*h_z) 151 dM_eff = 2 * np.pi * 0.5 * (r[i]*Sigma_eff[i] + r[i-1]*Sigma_eff[i-1]) * ,→(r[i]-r[i-1]) 152 M_b[i] = M_b[i-1] + dM_b 153 M_eff[i] = M_eff[i-1] + dM_eff 154 155 print(f"Check: integrated baryonic mass = {M_b[-1]/Msun:.2e} Msun (target ,→{M_disk_total/Msun:.2e})") 156 157 158 # ================================================================ 159 # FIGURE 1 Radial Profiles and Rotation Curves 160 # ================================================================ 161 162 r_kpc = r / kpc 163 fig, axs = plt.subplots(2, 2, figsize=(11, 9)) 164 plt.subplots_adjust(wspace=0.35, hspace=0.35) 165 166 plt.rcParams.update({ 167 ’font.size’: 21, 168 ’axes.titlesize’: 23, 169 ’axes.labelsize’: 21, 170 ’xtick.labelsize’: 19, 171 ’ytick.labelsize’: 19, 172 ’legend.fontsize’: 19 173 }) 174 175 # (a) Amplitude |A| 176 axs[0,0].plot(r_kpc, np.abs(A_r), color=’black’) 177 axs[0,0].set_title("(a) Amplitude $|A|$") 178 axs[0,0].set_xlabel("r [kpc]") 179 axs[0,0].set_ylabel(r"$|A(r)|$") 180 181 # (b) Mass density 182 Msun_pc3 = (pc**3) / Msun 183 axs[0,1].plot(r_kpc, rho_eff_mid*Msun_pc3, label=r"$\rho_{\rm eff}$", ,→color=’blue’) 184 axs[0,1].plot(r_kpc, rho_b*Msun_pc3, label=r"$\rho_{\rm b}$", ,→color=’orange’, linestyle=’--’) 185 axs[0,1].set_yscale(’log’) 186 axs[0,1].set_ylim(1e-8, 1e2) 187 axs[0,1].set_title("(b) Mass density") 188 axs[0,1].set_xlabel("r [kpc]") 189 axs[0,1].set_ylabel(r"Density [$M_\odot\,{\rm pc}^{-3}$]") 190 axs[0,1].legend() 191 192 # (c) Rotation curves 193 limit = r_kpc <= 50 194 axs[1,0].plot(r_kpc[limit], np.sqrt(G*M_b[limit]/r[limit])/1e3, ’--’, ,→label=r"$v_{\rm b}$", color=’orange’) 114 195 axs[1,0].plot(r_kpc[limit], np.sqrt(G*M_eff[limit]/r[limit])/1e3, ,→label=r"$v_{\rm eff}$", color=’blue’) 196 axs[1,0].plot(r_kpc[limit], ,→np.sqrt(G*(M_b[limit]+M_eff[limit])/r[limit])/1e3, label=r"$v_{\rm ,→tot}$", color=’black’) 197 axs[1,0].set_title("(c) Rotation curves (50 kpc)") 198 axs[1,0].set_xlabel("r [kpc]") 199 axs[1,0].set_ylabel("Velocity [km s$^{-1}$]") 200 axs[1,0].legend(loc=’center right’, frameon=True, facecolor=’white’, ,→framealpha=0.3, edgecolor=’none’, fontsize=16) 201 202 # (d) Mass ratio 203 axs[1,1].plot(r_kpc, M_eff / np.maximum(M_b, 1e-30), color=’purple’) 204 axs[1,1].set_title("(d) Enclosed mass ratio") 205 axs[1,1].set_xlabel("r [kpc]") 206 axs[1,1].set_ylabel(r"$M_{\rm eff}/M_{\rm b}$") 207 208 fig.suptitle("3D Helmholtz Halo Model (m = 1 Mode)", 209 fontsize=26, fontweight=’bold’, y=0.98) 210 plt.tight_layout(rect=[0, 0, 1, 0.96]) 211 plt.savefig("C:/Users/BryanH/Downloads/fig1_spiral_profiles.png", dpi=300, ,→facecolor=’white’) 212 plt.show() 213 214 # ================================================================ 215 # FIGURE 2 2D Midplane Maps (Amplitude, Phase, and Densities) 216 # ================================================================ 217 218 Nmap = 400 219 x = np.linspace(-50, 50, Nmap) * kpc 220 y = np.linspace(-50, 50, Nmap) * kpc 221 X, Y = np.meshgrid(x, y) 222 R = np.sqrt(X**2 + Y**2) 223 TH = np.arctan2(Y, X) 224 225 A_radial = np.interp(R, r, A_r) 226 rho_b_2D = np.interp(R, r, rho_b) 227 # --- corrected factor of 2 for consistency with rho_eff_mid --- 228 rho_eff_2D = (2 * Omega**2 / (8.0 * np.pi * G * c_light**2)) * ,→np.abs(A_radial)**2 229 A_2D = A_radial * np.exp(1j * TH) 230 231 amp_map = np.abs(A_2D) 232 phase_map = (np.angle(A_2D) + np.pi) / (2.0 * np.pi) 233 rho_tot_2D = rho_b_2D + rho_eff_2D 234 235 amp_norm = amp_map / np.max(amp_map) 236 phase_norm = phase_map 237 rho_b_norm = rho_b_2D / np.max(rho_b_2D) 238 rho_eff_norm = rho_eff_2D / np.max(rho_eff_2D) 239 rho_tot_norm = rho_tot_2D / np.max(rho_tot_2D) 240 241 fig, axs = plt.subplots(2, 2, figsize=(10, 10)) 242 115 243 axs[0,0].imshow(amp_norm, extent=[-50, 50, -50, 50], origin=’lower’, ,→cmap=’viridis’) 244 axs[0,0].set_title("(a) |A| amplitude\n(normalized)", fontsize=20) 245 axs[0,0].axis(’off’) 246 247 axs[0,1].imshow(phase_norm, extent=[-50, 50, -50, 50], origin=’lower’, ,→cmap=’gray’) 248 axs[0,1].set_title("(b) Phase arg(A)\n(normalized)", fontsize=20) 249 axs[0,1].axis(’off’) 250 251 rgb = np.zeros((*rho_eff_norm.shape, 3)) 252 rgb[..., 0] = rho_b_norm 253 rgb[..., 1] = rho_eff_norm 254 axs[1,0].imshow(rgb, extent=[-50, 50, -50, 50], origin=’lower’, vmin=0, ,→vmax=1) 255 axs[1,0].set_title("(c) Baryonic (red)\n and resonant (green)", ,→fontsize=20) 256 axs[1,0].axis(’off’) 257 258 axs[1,1].imshow(rho_tot_norm, extent=[-50, 50, -50, 50], origin=’lower’, ,→cmap=’gray’) 259 axs[1,1].set_title("(d) Combined density\n (normalized)", fontsize=20) 260 axs[1,1].axis(’off’) 261 262 fig.suptitle("Spiral Helmholtz Solutions\nNormalized Results (100 kpc 100 ,→kpc)", 263 fontsize=22, fontweight=’bold’, y=0.93) 264 plt.tight_layout(rect=[0, 0, 1, 0.95]) 265 plt.subplots_adjust(wspace=-0.3, hspace=0.25) 266 plt.savefig("C:/Users/BryanH/Downloads/fig2_spiral_maps.png", dpi=300, ,→facecolor=’white’) 267 plt.show() 268 269 # ================================================================ 270 # FIGURE 3 Dark Mass Density Slices (XY and XZ) 271 # ================================================================ 272 # The vertical XZ slice shows exponential decay exp(-2|z|/L_z), 273 # where L_z represents the geometric confinement scale. The vertical 274 # contrast therefore illustrates the halos stratification typically 275 # a few 10 kpc thick and is not rescaled by any normalization factor. 276 # ================================================================ 277 278 Nxy = 600 279 WIDTH = 200 280 extent_xy = [-WIDTH/2, WIDTH/2, -WIDTH/2, WIDTH/2] 281 extent_xz = [-WIDTH/2, WIDTH/2, -WIDTH/2, WIDTH/2] 282 283 x = np.linspace(-WIDTH/2, WIDTH/2, Nxy) * kpc 284 y = np.linspace(-WIDTH/2, WIDTH/2, Nxy) * kpc 285 z = np.linspace(-WIDTH/2, WIDTH/2, Nxy) * kpc 286 287 # --- Horizontal (XY) slice at z = 0 --- 288 X, Y = np.meshgrid(x, y) 289 R_xy = np.sqrt(X**2 + Y**2) 116 290 rho_eff_xy = np.interp(R_xy, np.real(r), np.real(rho_eff_mid)) 291 rho_eff_xy_norm = rho_eff_xy / np.max(rho_eff_xy) 292 293 # --- Vertical (XZ) slice using true L_z(r) profile --- 294 Xz, Z = np.meshgrid(x, z) 295 R_xz = np.abs(Xz) 296 297 # --- Ensure L_z and r are real 1D arrays --- 298 r_real = np.atleast_1d(np.real(r)) 299 rho_mid_real = np.atleast_1d(np.real(rho_eff_mid)) 300 301 # --- Case 1: L_z is an array --- 302 if np.ndim(L_z) > 0 and len(np.atleast_1d(L_z)) > 1: 303 Lz_real = np.atleast_1d(np.real(L_z)) 304 rho_mid_map = np.interp(R_xz.ravel(), r_real, ,→rho_mid_real).reshape(R_xz.shape) 305 Lz_map = np.interp(R_xz.ravel(), r_real, Lz_real).reshape(R_xz.shape) 306 else: 307 # --- Case 2: L_z is scalar --- 308 Lz_scalar = float(np.real(L_z)) if np.ndim(L_z) == 0 else ,→float(np.real(L_z[0])) 309 rho_mid_map = np.interp(R_xz.ravel(), r_real, ,→rho_mid_real).reshape(R_xz.shape) 310 Lz_map = np.full_like(R_xz, Lz_scalar) 311 312 # --- Apply vertical exponential envelope exp(-2|z|/L_z(r)) --- 313 rho_eff_xz = rho_mid_map * np.exp(-2.0 * np.abs(Z) / Lz_map) 314 rho_eff_xz_norm = rho_eff_xz / np.max(rho_eff_xz) 315 316 # --- Plot results --- 317 fig, axs = plt.subplots(1, 2, figsize=(12, 6), constrained_layout=False) 318 319 # (a) Horizontal slice XY 320 axs[0].imshow(rho_eff_xy_norm, extent=extent_xy, origin=’lower’, ,→cmap=’gray’, vmin=0, vmax=1) 321 axs[0].set_title("(a) Dark mass density\nHorizontal slice XY", fontsize=20) 322 axs[0].axis(’off’) 323 axs[0].set_aspect(’equal’, adjustable=’box’) 324 325 # (b) Vertical slice XZ (true or scalar L_z) 326 axs[1].imshow(rho_eff_xz_norm, extent=extent_xz, origin=’lower’, ,→cmap=’gray’, vmin=0, vmax=1) 327 axs[1].set_title("(b) Dark mass density\nVertical slice XZ", fontsize=20) 328 axs[1].axis(’off’) 329 axs[1].set_aspect(’equal’, adjustable=’box’) 330 331 fig.suptitle("Spiral 3D Halo Normalized |A| (200 kpc 200 kpc)", 332 fontsize=22, fontweight=’bold’, y=0.97) 333 fig.subplots_adjust(left=0.04, right=0.96, bottom=0.05, top=0.80, ,→wspace=0.10) 334 335 plt.savefig("C:/Users/BryanH/Downloads/fig3_darkmass_xy_xz.png", 336 dpi=300, facecolor=’white’) 337 plt.show() 117 Listing 1: Python code for the full 3D Helmholtz resonant–halo simulation. 118 C Appendix: Code for Simulating Elliptical Galaxies - Transport Theory This appendix provides the exact numerical implementation used to generate the ensemble–averaged resonant–halo predictions shown in Section 8.2.8. The simulation realizes the closed transport–law curvature–phonon halo derived analytically in Section 8.2.4 onwards, including the multiplicative gradient–energy contribution of Section 8.2.7. The stellar mass distribution is modeled as a Hernquist sphere with total stellar population N⋆ , individual mass m⋆ , and scale length a = Re/ 1 . 8153, giving the effective baryonic volume Vgal = 2πa3.(C.1) The elliptical transport closure is governed by four control parameters: •the transport damping length Ldamp,eff =ceff Γeff ,(C.2) fixed directly from the observed halo extent; •the resonant curvature wavelength λres,(C.3) which sets both the carrier frequency Ω = 2πceff λres (C.4) and the effective quality factor Qeff =Ω 2Γeff ; (C.5) • the speckle–shape factor α∼ 1, which controls the mean geometric coherence of the stellar sources; and • the mean coupling efficiency |X|2∼ 1, which sets the relative weighting of individual sources in the coherent ensemble field. Here Γ eff is the transport–damping rate associated with multiple stellar rescattering, and ceff is the curvature–phonon propagation speed; both are defined exactly as in the main text. 119 The ensemble–averaged effective gravitating density implemented in the code follows the closed transport–halo law, ρeff(r)=4πG c2 eff 2Qeff 2πΞ∇Genv,eff      r, r < Rmatch, e−2r/Ldamp,eff (4π)2r2, r ≥Rmatch. (C.6) where Genv,eff is the ensemble environmental gain, Qeff = Ω / (2Γ eff ) is the effective quality factor, and Ξ ∇ is the transport–enhanced gradient–energy boost all defined in Section 8.2.7. For numerical regularity, the small–radius linear core ⟨ρeff(r)⟩∝ris enforced for r < Rmatch (see Section 8.2.6). For numerical consistency and physical regularity, a small–radius linear core, ⟨ρeff(r)⟩∝r(r < Rmatch),(C.7) is enforced as in Section 8.2.6. All symbols in Eqs. (C.1) – (C.7) retain the exact same meaning and notation as in the main text; no new variables are introduced in the numerical implementation. The code uses α=|X|2=1 for the baseline models. The code below: •constructs ρeff(r), enclosed Meff(r), and circular velocities; •compares directly to the stellar mass distribution; and • produces both midplane composite maps and full–scale XY/XZ slices of the resonant dark halo. 1 2# ===================================================================== 3# Elliptical Galaxy Transport-Law Curvature-Phonon Speckle Halo 4# Written by: Bryan Hennelly 28 October 2025 5# ===================================================================== 6# 7# This script computes the ensemble-averaged effective halo density 8# generated by curvature-phonon transport in a densely populated, 9# phase-mixed stellar medium (elliptical galaxy). 10 # 11 # Physics summary: 12 # 13 # Stars are both sources and scatterers of curvature phonons. 14 # Multiple resonant scattering increases the coherent path length 15 # far beyond single-pass damping: _eff << . 16 # The coherent transport tail decays as exp(-2r/L_damp_eff)/r. 17 # The ensemble-stored curvature energy manifests as an 18 # effective gravitating density _eff(r). 19 # 20 # There are FOUR principal control parameters: 120 21 # 22 # (1) L_damp_eff transport damping length of the coherent tail 23 # (DATA-DRIVEN from observed halo extent) 24 # 25 # (2) _res structural resonant wavelength of curvature 26 # (PHYSICAL MICRO-INPUT controlling amplitude) 27 # 28 # (3) speckle-shape factor (~1), controls phase-averaged 29 # geometric coherence of stellar sources 30 # 31 # (4) |X| mean stellar coupling efficiency (~1), sets 32 # relative weighting of each source in the coherent ,→field 33 # 34 # From _res we obtain: 35 # 36 # = 2 c_eff / _res (carrier frequency) 37 # Q_eff = / (2 _eff) (effective quality factor) 38 # 39 # where _eff = c_eff / L_damp_eff is implied by transport theory, 40 # and c_eff is the curvature-phonon propagation speed we take to be 2e5 m/s 41 # 42 # FINAL transport-enhanced halo law: 43 # 44 # _eff(r) = 45 # _nabla * 46 # [ G v /(2 c c_eff) ] * 47 # L_damp_eff * Q_eff * 48 # (N_* m_* / V_gal) * |X| * 49 # { r , r<R_match ; exp[-2r/L_damp_eff]/r , rR_match } 50 # 51 # with: 52 # _eff = c_eff / L_damp_eff 53 # Q_eff = / (2 _eff) 54 # _nabla = 1 + (2 L_damp_eff / _res) 55 # G_env,eff = |X| (N_*/V_gal) (_res/2) (v / 2_eff) 56 # 57 # Interpretation: 58 # Tail SHAPE is dictated solely by L_damp_eff (observed). 59 # Tail AMPLITUDE is dictated solely by _res (microphysics). 60 # and |X| enter only as weak, O(1) normalization factors. 61 # 62 # Parameter-economical closure: 63 # One observational input halo size fixed 64 # One physical parameter halo mass fixed 65 # 66 # A small-r linear core is enforced by regularity of coherent flux. 67 # 68 # Outputs: 69 # _eff(r), circular velocities, enclosed masses 70 # 2D composite baryonic vs dark mass maps 71 # large XY/XZ slices of the full halo field 72 # 73 # All variable names follow the books notation for 1:1 correspondence 121 396 img_xy = log_norm2(rho_xy) 397 img_xz = log_norm2(rho_xz) 398 399 fig, axs = plt.subplots(1, 2, figsize=(12, 6)) 400 401 # (a) Horizontal slice XY 402 axs[0].imshow(img_xy, extent=[-extent, extent, -extent, extent], 403 origin=’lower’, cmap=’gray’, vmin=0, vmax=1) 404 axs[0].set_title("(a) Dark mass density\n Horizontal slice XY", 405 fontsize=20) 406 axs[0].axis(’off’) 407 axs[0].set_aspect(’equal’, adjustable=’box’) 408 409 # (b) Vertical slice XZ 410 axs[1].imshow(img_xz, extent=[-extent, extent, -extent, extent], 411 origin=’lower’, cmap=’gray’, vmin=0, vmax=1) 412 axs[1].set_title("(b) Dark mass density\n Vertical slice XZ", 413 fontsize=20) 414 axs[1].axis(’off’) 415 axs[1].set_aspect(’equal’, adjustable=’box’) 416 417 fig.suptitle("Elliptical Transport Halo Log-scaled Density Slices (1 Mpc)", 418 fontsize=22, fontweight=’bold’, y=0.999) 419 420 fig.subplots_adjust(left=0.02, right=0.98, 421 bottom=0.03, top=0.82, wspace=0.10) 422 423 plt.savefig("C:/Users/BryanH/Downloads/elliptical_fig3_slices_1Mpc.png", 424 dpi=300, facecolor=’white’) 425 plt.show() Listing 2: Python code for the ensemble-averaged resonant halo of an elliptical galaxy. — 128 D Appendix: Code for Simulating Cluster Galaxies This appendix provides the complete Python implementation used to generate Figures 10–12. The code follows directly from the radiative–transport formalism developed in Section 8.3.1, implementing the overlapping elliptical–halo convolution in reduced “2.5D” form for numerical stability. Each galaxy is represented by an exponential transport kernel K(s;L) = e−2s/L (4π)2s2,(D.1) ensemble–averaged over a lognormal distribution of damping lengths L with logarithmic dispersion σln L . The angularly averaged 2.5D convolution is evaluated as ρshape(r)=2πZ∞ 0 ngal(r′)r′⟨K(s)⟩µdr′, s2=r2+r′2−2rr′µ, (D.2) using Gauss–Legendre quadrature for the µ–integration. The galaxy number–density profile adopts an NFW-like form with a smooth splashback taper, ngal(r)=n0 Wsp(r) x(1+x)2, Wsp(r) = 1 1 + exp(r−Rsp)/∆sp,(D.3) ensuring a finite outer extent without discontinuities. The reduced convolution result ρshape(r) is normalized by the geometric factor Λgeom =Ngal M(elliptical) DM Z∞ 0 4πr2ρshape(r)dr ,(D.4) and scaled by the nonlinear reverberation gain Grev =1+Ngalη, (D.5) to yield the final effective density ρeff(r)=Grev Λgeom Sbase ρshape(r).(D.6) For the fiducial 1015 M⊙cluster, the adopted parameters are Lmean = 636 kpc, σln L= 2.5, Ngal = 300, η = 8.3×10−3, corresponding to Grev ≃ 3 . 5. The code produces three diagnostic figures: (i) cluster density and enclosed mass, (ii) circular velocity profile, and (iii) 2D distributions of ρeff and vc, reproducing the results presented in Figures 10–12. 129 1 2# ===================================================================== 3# Galaxy Cluster Static Transport-Overlap vs. NFW Benchmark 4# 5# Purpose: 6# Compute the cluster-scale effective transport density _eff(r) 7# arising from the collective overlap of mesoscopic transport halos 8# surrounding individual galaxies, and compare its emergent profile 9# with the canonical NFW dark-matter distribution. 10 # 11 # Model Overview: 12 # 1. K_tail(s) defines the dimensionless single-galaxy transport kernel, 13 # representing the exponential intensity tail of an elliptical halo. 14 # 15 # 2. The elliptical-halo microphysics determine the absolute scaling 16 # through Sbar_base, which sets the normalization of _eff(r) 17 # once the effective damping length L_damp_eff is chosen. 18 # Here, we adopt the same L_damp_eff as used in the single-elliptical 19 # simulation for internal consistency. 20 # 21 # 3. The spatial convolution is performed in a reduced (~2.5D) form 22 # for numerical stability. One geometric length factor is omitted 23 # and later restored by _geom, derived via global mass conservation. 24 # 25 # 4. A nonlinear reverberation gain factor, 26 # G_rev = 1 + N_gal * , 27 # accounts for weak cross-coherence between overlapping transport 28 # fields, where is the fractional coupling coefficient. 29 # 30 # Key Control Parameters: 31 # _lnL Logarithmic dispersion of transport lengths (heterogeneity): 32 # broadens the effective overlap and governs the curvature 33 # of _eff(r). Crucial for matching the NFW slope. 34 # 35 # Reverberation coupling coefficient: 36 # controls the degree of collective amplification among 37 # overlapping galaxy halos. Values in the range 1010 38 # yield physically plausible coherence gains. 39 # 40 # Once the galaxy microphysics and damping scale are fixed, 41 # _lnL and become the two primary tuning parameters determining 42 # the cluster-scale density and velocity profiles. 43 # 44 # Final Expression: 45 # _eff(r) = G_rev * _geom * Sbar_base * _shape(r) 46 # ===================================================================== 47 48 import numpy as np 49 import matplotlib.pyplot as plt 50 51 # --------------------------------------------------------------------- 52 # Fundamental constants and units 53 # --------------------------------------------------------------------- 130 54 G = 6.67430e-11 # m^3 kg^-1 s^-2 55 c = 2.99792458e8 # m/s 56 Msun= 1.98847e30 # kg 57 pc = 3.085677581491367e16 # m 58 kpc = 1.0e3 * pc 59 Mpc = 1.0e6 * pc 60 pi = np.pi 61 62 # --------------------------------------------------------------------- 63 # Cluster-scale parameters 64 # --------------------------------------------------------------------- 65 L_bar = 636.0 * kpc # main transport length [m] 66 sigma_lnL = 2.5 # lognormal scatter in L 67 R_cl = 2.0 * Mpc # cluster radius [m] 68 N_gal = 300 # number of elliptical galaxies 69 70 # Structural (cluster-scale) parameters 71 R_sp = 2.5 * Mpc # splashback radius [m] 72 Delta_sp = 0.5 * Mpc # splashback transition width [m] 73 r_s = 1.0 * Mpc # satellite NFW scale radius [m] 74 n0_sat = 2.0e-63 # satellite number-density scale [m^-3] 75 76 # --------------------------------------------------------------------- 77 # Elliptical-halo microphysics Sbar_base 78 # --------------------------------------------------------------------- 79 L_damp_eff = L_bar 80 c_eff_E = 2.0e5 81 Gamma_eff = c_eff_E / L_damp_eff 82 lambda_res = 27.0 * pc 83 Omega_E = 2.0 * pi * c_eff_E / lambda_res 84 Q_eff = Omega_E / (2.0 * Gamma_eff) 85 Xi_nabla = 1.0 + (2.0 * pi * L_damp_eff / lambda_res)**2 86 87 M_tot_E = 2.0e11 * Msun # total mass of one elliptical 88 m_star_E = Msun 89 N_star_E = M_tot_E / m_star_E 90 Re_E = 5.0 * kpc 91 a_h_E = Re_E / 1.8153 92 V_gal_E = 2.0 * pi * a_h_E**3 93 v_E = 200.0e3 94 alpha_E = 1.0 95 96 K_base = (G * alpha_E * v_E) / (2.0 * c**2 * c_eff_E**3) \ 97 * L_damp_eff * Q_eff * (N_star_E**2 * m_star_E**2 / V_gal_E) 98 Sbar_base = Xi_nabla * K_base 99 100 print("=== Single Elliptical Kernel Parameters ===") 101 print(f"L_bar = {L_bar/kpc:.1f} kpc, Q_eff = {Q_eff:.3e}, Xi_nabla = ,→{Xi_nabla:.3e}") 102 print(f"Sbar_base = {Sbar_base:.3e} [kg/m^3]\n") 103 104 # --------------------------------------------------------------------- 105 # Dimensionless transport kernel (shape only) 106 # --------------------------------------------------------------------- 131 107 def K_tail(s, L): 108 """Dimensionless transport kernel.""" 109 s_eps = 1.0 * pc 110 s_clp = np.maximum(s, s_eps) 111 return np.exp(-2.0 * s_clp / L) / (((4.0 * pi)**2) * s_clp**2) 112 113 def K_tail_ensemble(s, L_mean, sigma_ln): 114 """Lognormal-averaged kernel ensemble.""" 115 if sigma_ln <= 0.0: 116 return K_tail(s, L_mean) 117 xi = np.array([-2, -1, 0, 1, 2]) * sigma_ln 118 w = np.array([1, 4, 10, 4, 1], dtype=float) 119 w /= w.sum() 120 out = np.zeros_like(s) 121 for wi, Li in zip(w, L_mean * np.exp(xi)): 122 out += wi * K_tail(s, Li) 123 return out 124 125 Kfun = lambda s: K_tail_ensemble(s, L_bar, sigma_lnL) 126 127 # --------------------------------------------------------------------- 128 # Galaxy number-density profile (satellite population) 129 # --------------------------------------------------------------------- 130 def W_splash(r): 131 """Smooth splashback taper.""" 132 return 1.0 / (1.0 + np.exp((r - R_sp) / Delta_sp)) 133 134 def n_sat_NFW(r): 135 """Satellite population ~ 1/[x(1+x)^2].""" 136 x = np.maximum(r / r_s, 1e-9) 137 return W_splash(r) * n0_sat / (x * (1.0 + x)**2) 138 139 # --------------------------------------------------------------------- 140 # Reduced spherical convolution (~2.5D form) 141 # --------------------------------------------------------------------- 142 def angular_average_K(r, rp, Kfun): 143 """Angularly averaged kernel integral K(s) d.""" 144 mu, w = np.polynomial.legendre.leggauss(48) 145 r, rp = np.atleast_1d(r).reshape(-1,1), np.atleast_1d(rp).reshape(1,-1) 146 mu3, w3 = mu[:,None,None], w[:,None,None] 147 s = np.sqrt(r**2 + rp**2 - 2.0 * r * rp * mu3) 148 return np.sum(w3 * Kfun(s), axis=0) 149 150 def rho_from_population(n_pop): 151 """Compute reduced convolution _shape(r).""" 152 I = angular_average_K(r, r, Kfun) 153 integrand = n_pop * r 154 return 2.0 * pi * (I @ (integrand * dr)) 155 156 # --------------------------------------------------------------------- 157 # Radial grid and convolution 158 # --------------------------------------------------------------------- 159 r = np.geomspace(0.2 * kpc, 3.0 * Mpc, 480) 160 dr = np.gradient(r) 132 161 rho_eff_shape = rho_from_population(n_sat_NFW(r)) 162 163 # --------------------------------------------------------------------- 164 # (1) Geometric normalization _geom by mass conservation 165 # --------------------------------------------------------------------- 166 r_single = np.geomspace(0.1 * kpc, 100.0 * kpc, 400) 167 rho_single = Sbar_base * K_tail(r_single, L_bar) 168 M_DM_ellipt = 4.0 * np.pi * np.trapz(rho_single * r_single**2, r_single) 169 170 sum_M_DM_ellipticals = N_gal * M_DM_ellipt 171 M_cluster_shape = 4.0 * np.pi * np.trapz(rho_eff_shape * r**2, r) 172 Lambda_geom = sum_M_DM_ellipticals / (Sbar_base * M_cluster_shape) 173 174 print(f"Geometric normalization _geom = {Lambda_geom:.3e}") 175 176 # --------------------------------------------------------------------- 177 # (2) Reverberation gain factor G_rev = 1 + N_gal * 178 # --------------------------------------------------------------------- 179 eta = 0.00825 # fractional coherence (~1e-31e-2 typical) 180 G_rev = 1.0 + N_gal * eta 181 print(f"Reverberation gain G_rev = {G_rev:.2f}\n") 182 183 # --------------------------------------------------------------------- 184 # Final physical density 185 # --------------------------------------------------------------------- 186 rho_eff = G_rev * Lambda_geom * Sbar_base * rho_eff_shape 187 188 # --------------------------------------------------------------------- 189 # Enclosed mass and circular velocity 190 # --------------------------------------------------------------------- 191 M_eff = 4.0 * pi * np.cumsum(rho_eff * r**2 * dr) 192 v_c = np.sqrt(G * np.maximum(M_eff, 0.0) / r) / 1.0e3 # [km/s] 193 194 # --------------------------------------------------------------------- 195 # NFW benchmark 196 # --------------------------------------------------------------------- 197 H0 = 70.0 * 1000.0 / (1.0e6 * pc) 198 rho_crit = 3.0 * H0**2 / (8.0 * pi * G) 199 M200 = 1.0e15 * Msun 200 c200 = 4.0 201 R200 = (3.0 * M200 / (4.0 * pi * 200.0 * rho_crit))**(1/3) 202 r_s_DM = R200 / c200 203 rho_s = (200.0/3.0) * rho_crit * (c200**3) / (np.log(1+c200) - ,→c200/(1+c200)) 204 205 def rho_nfw(ri): 206 x = np.maximum(ri / r_s_DM, 1e-12) 207 return rho_s / (x * (1.0 + x)**2) 208 209 rho_nfw_vals = rho_nfw(r) 210 v_nfw = np.sqrt(G * (4.0*pi*np.cumsum(rho_nfw_vals*r**2*dr)) / r) / 1.0e3 211 212 # --------------------------------------------------------------------- 213 # Diagnostics 133 214 # --------------------------------------------------------------------- 215 print("=== Diagnostics ===") 216 print(f"L_bar = {L_bar/kpc:.1f} kpc") 217 print(f"sigma_lnL = {sigma_lnL:.2f} (lognormal width)") 218 print(f" (eta) = {eta:.4f} (reverberation coupling strength)") 219 print(f"G_rev = {G_rev:.2f} (1 + N_gal , nonlinear reverberation ,→gain)") 220 print(f"_geom = {Lambda_geom:.3e} (mass-conserving normalization)") 221 print(f"max _eff = {rho_eff.max():.3e} kg/m at r = ,→{r[np.argmax(rho_eff)]/kpc:.2f} kpc") 222 print(f"max v_c = {v_c.max():.2f} km/s at r = ,→{r[np.argmax(v_c)]/kpc:.2f} kpc\n") 223 224 225 r_Mpc = r / Mpc 226 # ================================================================ 227 # FIGURE 1 Density and Enclosed Mass (Linear + Log, raised titles) 228 # ================================================================ 229 230 from matplotlib.ticker import ScalarFormatter 231 formatter = ScalarFormatter(useMathText=True) 232 formatter.set_powerlimits((-2, 4)) 233 234 r_Mpc = r / Mpc 235 rho_eff_msun_pc3 = rho_eff * (pc**3 / Msun) 236 rho_nfw_msun_pc3 = rho_nfw_vals * (pc**3 / Msun) 237 238 # ---------- Figure setup ---------- 239 fig1, axs1 = plt.subplots(2, 2, figsize=(13, 11)) 240 plt.subplots_adjust(wspace=0.35, hspace=0.35) 241 242 # ---------- Font sizes ---------- 243 title_fs = 26 244 label_fs = 22 245 tick_fs = 18 246 legend_fs = 20 247 suptitle_fs = 32 248 title_pad = 28 # raised titles for all panels 249 250 # (a) Density (linear) 251 axs1[0,0].plot(r_Mpc, rho_eff_msun_pc3, color=’navy’, lw=2, ,→label=r"$\rho_{\rm eff}$") 252 axs1[0,0].plot(r_Mpc, rho_nfw_msun_pc3, ’--’, color=’darkorange’, lw=2, ,→label="NFW") 253 axs1[0,0].set_title("(a) Density (linear)", fontsize=title_fs, pad=title_pad) 254 axs1[0,0].set_xlabel("r [Mpc]", fontsize=label_fs) 255 axs1[0,0].set_ylabel(r"Density [$M_\odot\,{\rm pc^{-3}}$]", ,→fontsize=label_fs) 256 axs1[0,0].legend(fontsize=legend_fs) 257 axs1[0,0].tick_params(axis=’both’, which=’major’, labelsize=tick_fs) 258 axs1[0,0].xaxis.set_major_formatter(formatter) 259 axs1[0,0].yaxis.set_major_formatter(formatter) 260 261 # (b) Density (loglog) 134 262 axs1[0,1].loglog(r_Mpc, rho_eff_msun_pc3, color=’navy’, lw=2) 263 axs1[0,1].loglog(r_Mpc, rho_nfw_msun_pc3, ’--’, color=’darkorange’, lw=2) 264 axs1[0,1].set_title("(b) Density (loglog)", fontsize=title_fs, pad=title_pad) 265 axs1[0,1].set_xlabel("r [Mpc]", fontsize=label_fs) 266 axs1[0,1].set_ylabel(r"Density [$M_\odot\,{\rm pc^{-3}}$]", ,→fontsize=label_fs) 267 axs1[0,1].tick_params(axis=’both’, which=’major’, labelsize=tick_fs) 268 269 # (c) Enclosed mass (linear) 270 M_nfw = 4.0 * np.pi * np.cumsum(rho_nfw_vals * r**2 * dr) 271 axs1[1,0].plot(r_Mpc, M_eff/Msun, color=’crimson’, lw=2, label=r"$M_{\rm ,→eff}$") 272 axs1[1,0].plot(r_Mpc, M_nfw/Msun, ’--’, color=’orange’, lw=2, label="NFW") 273 axs1[1,0].set_title("(c) Enclosed mass (linear)", fontsize=title_fs, ,→pad=title_pad) 274 axs1[1,0].set_xlabel("r [Mpc]", fontsize=label_fs) 275 axs1[1,0].set_ylabel(r"$M(<r)\,[M_\odot]$", fontsize=label_fs) 276 axs1[1,0].legend(fontsize=legend_fs) 277 axs1[1,0].tick_params(axis=’both’, which=’major’, labelsize=tick_fs) 278 axs1[1,0].xaxis.set_major_formatter(formatter) 279 axs1[1,0].yaxis.set_major_formatter(formatter) 280 axs1[1,0].yaxis.get_offset_text().set_fontsize(tick_fs) 281 282 # (d) Enclosed mass (loglog) 283 axs1[1,1].loglog(r_Mpc, M_eff/Msun, color=’crimson’, lw=2) 284 axs1[1,1].loglog(r_Mpc, M_nfw/Msun, ’--’, color=’orange’, lw=2) 285 axs1[1,1].set_title("(d) Enclosed mass (loglog)", fontsize=title_fs, ,→pad=title_pad) 286 axs1[1,1].set_xlabel("r [Mpc]", fontsize=label_fs) 287 axs1[1,1].set_ylabel(r"$M(<r)\,[M_\odot]$", fontsize=label_fs) 288 axs1[1,1].tick_params(axis=’both’, which=’major’, labelsize=tick_fs) 289 290 # ---------- Main title ---------- 291 fig1.suptitle("Cluster Density and Enclosed Mass", 292 fontsize=suptitle_fs, fontweight=’bold’, y=0.98) 293 294 # ---------- Save & display ---------- 295 plt.tight_layout(rect=[0, 0, 1, 0.96]) 296 plt.savefig("C:/Users/BryanH/Downloads/fig1_density_mass.png", 297 dpi=300, facecolor=’white’) 298 plt.show() 299 300 # ================================================================ 301 # FIGURE 2 Circular Velocity Profile (Linear Only, matching Figure 1 font ,→sizes) 302 # ================================================================ 303 304 fig2, ax2 = plt.subplots(figsize=(14, 6)) 305 306 # Match Figure 1 font sizes 307 title_fs = 26 308 label_fs = 22 309 tick_fs = 18 310 legend_fs = 20 135 311 suptitle_fs = 32 312 313 ax2.plot(r_Mpc, v_c, color=’crimson’, lw=2, label="Transport overlap") 314 ax2.plot(r_Mpc, v_nfw, ’--’, color=’orange’, lw=2, label="NFW") 315 316 ax2.set_title("(a) Circular velocity (linear)", fontsize=title_fs, pad=20) 317 ax2.set_xlabel("r [Mpc]", fontsize=label_fs) 318 ax2.set_ylabel(r"$v_c$[km s$^{-1}$]", fontsize=label_fs) 319 ax2.legend(fontsize=legend_fs, loc=’best’) 320 ax2.tick_params(axis=’both’, which=’major’, labelsize=tick_fs) 321 ax2.xaxis.set_major_formatter(formatter) 322 ax2.yaxis.set_major_formatter(formatter) 323 324 fig2.suptitle("Circular Velocity Profile Transport vs NFW", 325 fontsize=suptitle_fs, fontweight=’bold’, y=0.97) 326 327 plt.tight_layout(rect=[0.02, 0, 0.98, 0.95]) 328 plt.savefig("C:/Users/BryanH/Downloads/fig2_velocity.png", 329 dpi=300, facecolor=’white’) 330 plt.show() 331 332 333 # ================================================================ 334 # FIGURE 3 2D Distributions over r_max (final with extra top space) 335 # ================================================================ 336 337 # ---------- Construct 2D grid ---------- 338 n = 700 339 r_max_Mpc = r.max() / Mpc 340 extent = r_max_Mpc 341 x2 = np.linspace(-extent, extent, n) 342 y2 = np.linspace(-extent, extent, n) 343 X, Y = np.meshgrid(x2, y2) 344 R2 = np.sqrt(X**2 + Y**2) * Mpc # convert to meters 345 346 # ---------- Interpolate 1D radial profiles ---------- 347 rho_2D = np.interp(R2, r, rho_eff, right=np.nan) 348 v_2D = np.interp(R2, r, v_c, right=np.nan) # [km/s] 349 350 # ---------- Normalize effective density ---------- 351 rho_log = np.log10(rho_2D + 1e-40) 352 rho_norm = (rho_log - np.nanmin(rho_log)) / (np.nanmax(rho_log) - ,→np.nanmin(rho_log)) 353 354 # ---------- Figure setup ---------- 355 fig3, axs3 = plt.subplots(1, 2, figsize=(10, 6)) 356 title_fs = 22 357 358 # (a) Effective density 359 im0 = axs3[0].imshow( 360 rho_norm, extent=[-r_max_Mpc, r_max_Mpc, -r_max_Mpc, r_max_Mpc], 361 origin=’lower’, cmap=’inferno’, vmin=0, vmax=1 362 ) 363 axs3[0].set_title(r"(a) $\log_{10}\rho_{\mathrm{eff}}$" "\n(normalized)", 136 364 fontsize=title_fs, pad=10) 365 axs3[0].axis(’off’) 366 fig3.colorbar(im0, ax=axs3[0], fraction=0.046, pad=0.04) 367 368 # (b) Velocity field (linear) 369 im1 = axs3[1].imshow( 370 v_2D, extent=[-r_max_Mpc, r_max_Mpc, -r_max_Mpc, r_max_Mpc], 371 origin=’lower’, cmap=’viridis’ 372 ) 373 axs3[1].set_title(r"(b) Circular velocity $v_c$" "\n(linear)", 374 fontsize=title_fs, pad=10) 375 axs3[1].axis(’off’) 376 cbar1 = fig3.colorbar(im1, ax=axs3[1], fraction=0.046, pad=0.04) 377 cbar1.set_label(r"$v_c$[km s$^{-1}$]", fontsize=14) 378 379 # ---------- Main title ---------- 380 fig3.suptitle( 381 rf"2D Distributions over {r_max_Mpc:.1f} Mpc", 382 fontsize=28, fontweight=’bold’, y=0.995 383 ) 384 385 # ---------- Layout & spacing ---------- 386 plt.subplots_adjust(left=0.03, right=0.97, top=0.86, bottom=0.05, 387 wspace=0.20, hspace=0.25) 388 389 # ---------- Save & display ---------- 390 plt.savefig("C:/Users/BryanH/Downloads/fig3_2D_distributions.png", 391 dpi=300, facecolor=’white’, bbox_inches=’tight’) 392 plt.show() Listing 3: Python code for the quarter–wave spherical cluster resonator. References [1] Albert Einstein. Die grundlagen der allgemeinen. Relativitatsteorie, Annale der Physic, 49:769, 1916. [2] Robert M Wald. General relativity. University of Chicago press, 2024. [3] John P. Cox. Theory of Stellar Pulsation. Princeton Series in Astrophysics. Princeton University Press, 1980. [4] Arnaud Derode, Arnaud Tourin, and Mathias Fink. Random multiple scattering of ultrasound. i. coherent and ballistic waves. Physical Review E, 64(3):036605, 2001. [5] Arnaud Derode, Arnaud Tourin, and Mathias Fink. Random multiple scattering of ultrasound. ii. is time reversal a self-averaging process? Physical Review E, 64(3):036606, 2001. 137