The Spacetime Resonance Theorem and the Origin of Dark Mass
Abstract
This work applies a well-known idea from acoustics and optics—multiple scattering in a medium—to gravity itself. In linearised general relativity the gravitational field is determined by a retarded Green function, but in most applications its “causal tail’’ (the part inside the light cone) is discarded. When this tail is kept and combined with the linear damped-oscillator susceptibility we introduce for baryonic matter (mathematically identical to the response used in acoustical and optical scattering theory) spacetime behaves as an effective transport medium for weak curvature waves. From this simple starting point we derive a single theorem—the Spacetime Resonance Theorem—showing how coherent, time-varying baryonic motion generates a small resonant curvature field. Multiple gravitational scattering causes part of this oscillatory field to accumulate as a slow, cycle-averaged component. This stored curvature energy behaves exactly like an effective dark-mass density. The resulting framework naturally reproduces the observed dark-mass profiles of spiral galaxies, ellipticals, clusters, and filaments; predicts “local silence’’ and “homogeneous silence’’ (matching the Bullet Cluster and the CMB); and remains fully consistent with standard GR. No new particles or modifications of gravity are introduced: dark-mass behaviour emerges from the resonant, multiple-scattering response of baryonic structure within ordinary weak-field GR.
Full text
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 We show that the full retarded Green–function structure of linearised general relativity, when combined with the ordinary damped–oscillator susceptibility of baryonic matter, 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 Constitution of the Effective Gravitational Medium . . . . . . . . 6 2.2 Relation to Known Gravitational Oscillation Modes and Precedents 10 3 Resonant Field Dynamics 11 3.1 ChannelSeparation.......................... 11 3.2 Spacetime as an Effective Reversible Medium . . . . . . . . . . . 12 3.3 Collective Mode Dynamics and the Masson . . . . . . . . . . . . . 13 3.4 Linear Spatial Perturbations and the Propagation Law . . . . . . 14 4 Gravitational Scattering Theory 16 4.1 Single–Scatterer Response . . . . . . . . . . . . . . . . . . . . . . 17 4.2 Multiple Scattering and the Foldy–Lax Hierarchy . . . . . . . . . 17 4.3 Effective–Medium and Ensemble Averaging . . . . . . . . . . . . . 19 4.4 Transport and Diffusion Limits . . . . . . . . . . . . . . . . . . . 20 4.5 Unified Interpretation . . . . . . . . . . . . . . . . . . . . . . . . . 21 5 Spatiotemporal Feedback and Resonant Coherence 22 5.1 Constituent Entities & States of the Feedback System . . . . . . . 22 5.2 Feedback Channels and Interactions . . . . . . . . . . . . . . . . . 24 5.3 Temporal Coherence and Mutual Susceptibility . . . . . . . . . . 27 5.3.1 Phase Locking and Mutual Susceptibility . . . . . . . . . . 27 5.3.2 Global Frequency Selection in Inhomogeneous Discs . . . . 29 5.3.3 Polarisation of Temporal Coherence . . . . . . . . . . . . . 30 6 Spatial Coherence of the Resonant Channel 31 6.1 Kirchhoff–Sommerfeld Representation . . . . . . . . . . . . . . . . 32 6.2 Effective Potential and Energy Storage in the Resonant Medium . 34 6.3 Single Point Source in a Resonant Channel–1 Background . . . . 35 6.4 Interference of Two Point Sources . . . . . . . . . . . . . . . . . . 37 6.5 Temporal Memory and the Steady–State Reduction . . . . . . . . 42 7 Steady–State (Attractor) Solutions 43 7.1 Spiral Galaxies as Coherent m= 1 (Laguerre–Gaussian) Modes . . . . . . . . . . . . . . . . . . . . . 44 7.1.1 Observing the Steady–State Solution . . . . . . . . . . . . 44 7.1.2 The First Moment . . . . . . . . . . . . . . . . . . . . . . 47 7.1.3 The Thin–Disk Model: Solving the 2D Helmholtz Equation RevealstheHalo ....................... 49 7.1.4 Extending the Thin–Disk to Three Dimensions . . . . . . . 51 7.1.5 Gain from Incoherent Foldy–Lax Scattering . . . . . . . . 52 2
7.1.6 Stored Energy, Dark–Mass Equivalence, and EnvironmentalGain ............................ 53 7.1.7 Numerical Results: 3D Thin–Disk Model . . . . . . . . . . 54 7.1.8 Physical Interpretation . . . . . . . . . . . . . . . . . . . . 58 7.2 Ellipticals as Random Resonant Media . . . . . . . . . . . . . . . 62 7.2.1 Observing the Statistical Steady State . . . . . . . . . . . 62 7.2.2 Complex Baryonic Forcing with Memory . . . . . . . . . . 64 7.2.3 Four–dimensional speckle geometry of the gravitational field 66 7.2.4 Ensemble reduction: masson–speckle summation . . . . . . 67 7.2.5 Transport Mean Free Path and Collective Damping . . . . 70 7.2.6 Core Regularization in the Non–Transport Regime . . . . 71 7.2.7 Gradient–Energy Boost and Final Transport–Regularized HaloLaw ........................... 73 7.2.8 Numerical Results . . . . . . . . . . . . . . . . . . . . . . 74 7.2.9 Halo Formation from Coherent Energy Flux . . . . . . . . 77 7.2.10 Physical Interpretation . . . . . . . . . . . . . . . . . . . . 78 7.3 Clusters as Radiative Transport Extensions of Ellipticals . . . . . 81 7.3.1 Observing the Statistical Transport Equilibrium . . . . . . 81 7.3.2 Mathematical Formulation: Overlapping Elliptical TransportEnvelopes ........................ 84 7.3.3 Numerical Results . . . . . . . . . . . . . . . . . . . . . . 89 7.3.4 Physical Interpretation . . . . . . . . . . . . . . . . . . . . 92 8 Local Silence and the Preservation of Classical Gravity 94 9 Prediction: Absence of Gain in Homogeneous Systems 96 9.1 Three Distinct Regimes for Homogeneous Clouds . . . . . . . . . 96 9.2 Theorem: Homogeneous Media Produce No Resonant Gain . . . . 97 9.3 Interpretation and Connection to Perturbation Theory . . . . . . 98 10 Emergence of Multiple Scattering from a Weakly Inhomogeneous Continuous Medium 98 10.1 A Brief Review of Garvitational Multiple Scatter Theory . . . . . 98 10.2 The Resonant Frequency in a Continuous Medium . . . . . . . . . 100 10.3 Weak Random Potential and the Ballistic Regime . . . . . . . . . 100 10.4 Nonlinear Growth Toward the Transport Threshold . . . . . . . . 101 10.5 Emergent Resonant Patches and Effective Discreteness . . . . . . 101 10.6 Linear Limit and Connection to the CMB Power Spectrum . . . . 102 11 Conclusion 104 12 Dedication and Acknowledgement 104 A Appendix: Effective–Medium Derivation 105 3
B Appendix: Bessel and Hankel Functions Used in the Main Text 107 B.1 Cylindrical Bessel and Hankel Functions . . . . . . . . . . . . . . 107 B.2 Spherical Bessel Functions . . . . . . . . . . . . . . . . . . . . . . 108 B.3 Physical Interpretation . . . . . . . . . . . . . . . . . . . . . . . . 108 C Appendix: Code for Simulating the 3D Thin–Disk Resonant Halo 109 D Appendix: Code for Simulating Elliptical Galaxies - Transport Theory 118 E Appendix: Code for Simulating Cluster Galaxies 128 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 finite–memory, dispersive kernel already present in linearised GR once no quasistatic truncation is applied. The central idea of this work is to describe how baryonic matter interacts with this causal kernel. Each gravitating element responds to periodic curvature perturbations as a damped gravitational oscillator with a frequency–dependent linear susceptibility χ⋆ ( ω ). When many such elements participate, their susceptibility–controlled delayed re-emission of curvature perturbations produces a gravitational analogue of classical multiple scattering. Thus the effective response of an extended baryonic system does not arise from modifying GR, but from the interaction between (i) the retarded Green function of linearised GR and (ii) the damped-oscillator susceptibility of ordinary baryonic matter. From these two ingredients an effective medium emerges. Repeated rescattering of weak curvature perturbations renormalises the propagator in the standard Dyson sense, producing a collective field with a characteristic restoring frequency Ω, damping rate Γ, and effective propagation speed ceff . Motivated by well–known gravitational oscillators—stellar pulsation modes, Jeans modes, disk breathing modes—it is natural to take Ω 2∝G¯ρb , where ¯ρb is a suitable coarse–grained baryonic density. This identification is not imposed by the field equation; it is a physically well–motivated correspondence with established gravitational phenomena. The resulting collective curvature response, Φ res , is oscillatory and dynamically sustained. Because the retarded kernel carries finite memory, this oscillatory field does not cancel perfectly from cycle to cycle. Its cycle–averaged curvature energy accumulates extremely slowly, defining an effective mass density ρeff ∝ ⟨| Φ res|2⟩ . This stored component contributes to the usual Poisson equation while the oscillatory part remains locally silent. Thus the theory preserves all laboratory, Solar–System, and stellar tests of gravity while allowing extended systems—spirals, ellipticals, and clusters—to acquire substantial effective dark–mass envelopes. Structured baryonic systems also act as gravitational scatterers. Coherent motion 5
in rotating disks, orbiting cores, and clustered assemblies generates constructive rescattering of curvature perturbations. Depending on spatial contrast and coherence scale, this yields two observationally distinct regimes: large–scale coherent modes in spirals, and diffusive effective–medium halos in ellipticals and clusters. In contrast, homogeneous or weakly perturbed environments possess no spatial contrast and therefore no rescattering: they remain resonantly silent. This homogeneous silence explains why the early Universe and the CMB follow standard cosmology: no contrast means no resonant channel. 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 capable of describing spiral galaxies, ellipticals, clusters, and cosmological structure without new fields or non–baryonic particles. The resulting Spacetime Resonance Theorem formalises the emergence of the renormalised resonant mode and provides the microscopic foundation for the stored curvature energy that behaves observationally as dark mass. 2 Spacetime Resonance Theorem 2.1 Constitution of the Effective Gravitational Medium The analysis in this paper remains entirely within the weak–field limit of general relativity and introduces no modifications of the Einstein equations, no new fields, and no physical properties assigned to the vacuum. The only refinement beyond the standard Newtonian approximation is the decision to retain the full retarded Green–function structure of the linearised gravitational field, rather than truncating it to the instantaneous Poisson form. In extended, inhomogeneous baryonic systems this retarded kernel produces nontrivial feedback: coherent mass motions both source and rescatter weak curvature perturbations. When many baryonic elements participate, repeated rescattering generates a collective response analogous to effective–medium behaviour in acoustics or electromagnetic wave transport. In this sense the matter distribution, not the vacuum, plays the role of the medium. All dispersion, memory, attenuation, and resonant behaviour arise solely from the baryonic ensemble’s gravitational susceptibility. The effective–medium viewpoint used throughout this paper is therefore entirely conventional: it derives from applying standard multiple scattering, homogenisation, and Green–function methods to the weak–field Einstein propagator. The following structural principles summarise this interpretation: 1. Weak–field GR with nonlocal propagators. The linearised potential is a retarded convolution of the baryonic source with the general–relativistic Green function. This nonlocality introduces finite response time and allows 6
multiple scattering to accumulate, but it does not endow spacetime with material properties. 2. Baryonic susceptibility under periodic forcing. A time–periodic component of the baryonic density, ρb,1 (x , t ), acts as the coherent source of the curvature field. The baryonic ensemble responds through a frequency– dependent gravitational susceptibility χ⋆ ( ω ), which determines how efficiently curvature perturbations are regenerated and rescattered and therefore provides the microscopic origin of resonant behaviour. 3. Lossless gravitational rescattering by inhomogeneities. Spatial variations in the baryonic density act as weak gravitational scatterers. Collectively they produce interference, transport, and partial coherence of curvature perturbations, directly analogous to acoustic or electromagnetic multiple scattering. 4. Self–energy and effective propagation. The ensemble of scatterers modifies the propagator through a Dyson–type relation, G−1 eff =G−1 0−Σ(ω, k), where Σ encodes collective rescattering. A long–wavelength, low–frequency expansion of Σ yields a local, second–order PDE for the collective curvature response with coefficients (Ω,Γ, ceff ). 5. Resonant self–organisation. A central working hypothesis of this paper is that the effective macroscopic coefficient a0 in the self–energy expansion yields a restoring frequency scaling as Ω2∝G¯ρb, where ¯ρb is the mean baryonic density in a coherence domain. This mirrors the known √Gρ scaling of stellar and galactic oscillation modes and leads naturally to long–lived, phase–locked collective oscillations. 6. Effective curvature field and stored energy. The coherent curvature perturbation Φ res generated by multiple scattering defines a renormalised weak–field potential—the effective curvature field. Its cycle–averaged intensity stores gravitational energy, EΦ∝ |Φres|2, and this stored energy appears in the Poisson channel as an effective mass density, ρeff ∝ ⟨|Φres|2⟩. No new gravitational degree of freedom is introduced: the effective curvature field is the macroscopic, homogenised response of standard weak–field GR in the presence of baryonic susceptibility and rescattering, and the associated energy is entirely a property of that field. 7
These principles depict galaxies, clusters, and other extended baryonic systems as self–consistent gravitational resonators whose macroscopic behaviour follows from standard weak–field GR combined with multiple scattering and homogenisation. The following theorem formalises this effective–medium resonance and derives the local resonant evolution equation for the collective curvature response. Theorem (Effective–Medium Resonance of the Weak— Field Gravitational Potential). Consider the linearised weak–field gravitational potential Φ(x) = −4πGZGret(x, x′)ρb(x′)d4x′,(1) where Gret is the full GR retarded Green function. Decompose the matter density into a static and coherent component, ρb(x, t) = ρb,0(x)+ρb,1(x, t).(2) Assume the baryons form a statistically inhomogeneous ensemble of weak gravitational scatterers whose periodic internal rearrangement is described by a linear susceptibility χ⋆ ( ω ). Multiple scattering renormalises the propagator according to the Dyson relation G−1 eff (ω, k) = G−1 0(ω, k)−Σ(ω, k),(3) where Σis the ensemble self–energy. In the long–wavelength, low–frequency limit, Σadmits the local expansion Σ(ω, k)≃a0+a1(iω)+a2ω2+b2k2,(4) with coefficients determined by χ⋆(ω)and the spatial statistics of the scatterers. Transforming Eq. (3) back to configuration space using (4) yields the homogenised evolution equation for the coherent oscillatory field, ∂2 tΦres + 2Γ ∂tΦres + Ω2Φres −c2 eff∇2Φres = 4πG ρb,1(x, t),(5) with identifications Ω2=a0,Γ = a1, c2 eff =c2(1 −b2).(6) The resonant field is purely dynamical: when the coherent driving vanishes, ρb,1(x, t)→0 =⇒Φres(x, t)→0,(7) and the remaining gravitational potential is the usual Poisson response to ρb,0. Thus, standard weak–field GR, when applied to an inhomogeneous baryonic ensemble with gravitational susceptibility and multiple scattering, supports a macroscopic resonant mode with parameters (Ω , Γ , ceff )determined entirely by the matter distribution. No modifications of general relativity or additional gravitational fields are required. 8
Condensed derivation (full details in Appendix A). We outline the steps connecting the weak–field retarded potential to the local resonant equation (5). The resonant component of the potential is obtained by restricting the retarded solution (1) to the coherent source: Φres(x, t) := −4πGZGret(x, x′)ρb,1(x′, t′)d4x′,(C1) which vanishes whenever ρb,1= 0, consistent with (7). Fourier transforming gives the linear response relation ˜ Φres(ω, k) = −4πG ˜ Geff(ω, k) ˜ρb,1(ω, k),(C2) where the effective propagator obeys the Dyson equation ˜ G−1 eff =˜ G−1 0−Σ(ω, k).(C3) In the long–wavelength, low–frequency regime the self–energy takes the local form Σ(ω, k)≃a0+a1(iω)+a2ω2+b2k2,(C4) with coefficients set by the microscopic susceptibility χ⋆ ( ω ) and the scatterer statistics. Substituting (C4) into (C3) yields the standard resonant structure ˜ G−1 eff ∝ −ω2−2iΓω+ Ω2+c2 effk2,(C5) with Ω2=a0,Γ = a1, c2 eff =c2(1 −b2),(C6) as in (6). Inserting (C5) into (C2) and transforming back to configuration space gives the homogenised resonant equation (∂2 t+ 2Γ ∂t+ Ω2−c2 eff∇2)Φres = 4πG ρb,1(x, t),(C7) identical to Eq. (5). Since (C2) implies ˜ Φres → 0 whenever ˜ρb,1→ 0, the resonant field has no static limit, again matching (7) . The parameters (Ω , Γ , ceff ) arise solely from baryonic susceptibility and multiple scattering; no modification of GR is involved. 9
compressible media, where frequency response is controlled by the effective bulk modulus: high–frequency density perturbations cannot drive the medium. Together, Eqs. (13) – (19) 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 Ω,(21) links the microscopic phase–lag parameter Γ to the macroscopic persistence of curvature oscillations. As demonstrated in Section 5, 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 4. 4 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 7). 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. 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 16
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 paper. 4.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⋆χ⋆(Ω),(22) 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 7.1.8 (spiral galaxies) and Section 7.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 (22) 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. 4.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 17
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 3.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,(23) so that the ensemble behaves as a single collective radiator with intensity Icoh ∝N2|χ⋆(Ω)|2|Ainc|2,(24) It is convenient to quantify this as an environmental gain Genv ≡|Atot|2 |Ainc|2,Genv =N2|χ⋆(Ω)|2(coherent limit).(25) 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),(26) with t(Ω) = 4πGm⋆χ⋆(Ω). In operator notation this becomes A= (1 −tGΩ)−1Ainc,(27) whose Neumann expansion, A=Ainc +tGΩAinc + (tGΩ)2Ainc +··· ,(28) sums all orders of mutual gravitational re–scattering. In the limit Γ → 0 (purely real GΩ ), the response is strictly coherent and Eq. (27) reduces to the linear scaling of Eq. (23). 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,(29) 18
yielding the standard incoherent scaling Genv =N|χ⋆(Ω)|2,pGenv =√N|χ⋆(Ω)|.(30) The spiral–disk regime discussed in Section 7.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. (26) 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. 4.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 (26) 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+ Σ(Ω),(31) 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,(32) 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. (30) 19
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 ⊥ℓ∥ ,(33) 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 ⊥ℓ∥ ,(34) 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 ≃Γ,(35) 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 (31) 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 7.2. 4.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),(36) 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 20
specifying how gravitational perturbation energy is redirected by baryonic inclusions. Equation (36) 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. (36) reduces to the usual diffusion equation for the ensemble-averaged intensity: ∇2I(x) = I(x) L2 damp,eff ,(37) where the effective damping (or coherence–decay) length is Ldamp,eff =ceff Γeff .(38) 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. (36) – (38) 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. 4.5 Unified Interpretation The coherent (Foldy–Lax) and diffusive (transport) limits developed above arise from the same governing relation, A=1−t(Ω) GΩ−1Ainc,(39) 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 4, 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, producing the global m = 1 spiral modes described in Section 7.1. Here the medium behaves as a phase-ordered resonator, and the Foldy–Lax limit applies. 21
When phases are randomized (as in irregular, thick, or pressure-supported systems), the tensorial response averages to its isotropic trace and the same operator (39) 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. 5 Spatiotemporal Feedback and Resonant Coherence 5.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. 22
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 ). 23
5.2 Feedback Channels and Interactions The curvature–matter system evolves through two coupled feedback channels, both arising from the same weak–field gravitational dynamics. Channel 1 is the resonant, wave–based response driven by the oscillatory baryonic density ρb,1 and governed by the causal, memory–bearing propagator derived in Section 6. Channel 2 is the Newtonian, effectively instantaneous response generated by the slowly varying baryonic density ρb,0 together with the accumulated effective (dark) mass ρeff . Their interaction forms the dual feedback architecture shown in Fig. 1. (1) Resonant Channel (Channel 1). The oscillatory component ρb,1 excites the resonant curvature field Φ res through the causal Green function GΩ of the damped resonant curvature equation. As developed in Section 6, GΩ encodes propagation, phase delay, and finite causal memory at the effective speed ceff . Because many baryons act as coherent emitters and scatterers, the resonant field is statistically amplified by the environmental gain Genv, giving Φres ∝pGenv GΩ∗ρb,1,(40) with ∗ denoting causal convolution. The interaction is reversible: resonant curvature perturbations exchange energy and phase with the baryons that generate them. Over the memory time τmem = Γ −1 , this oscillatory exchange accumulates into an effective (dark) mass density, ρeff ∝ ⟨|Φres|2⟩t,(41) which then enters the Newtonian response of Channel 2. (2) Newtonian Channel (Channel 2). The background baryonic density ρb,0 together with the stored effective mass ρeff source the Newtonian potential Φ N via the Poisson equation. This potential propagates effectively instantaneously at speed c within the weak–field regime and establishes the local geodesic structure. The resulting velocity field vredistributes ρb,0 and modulates the oscillatory component ρb,1 , thereby resetting the pattern of resonant sources that feed Channel 1. This constitutes the stress–advection loop: the accumulated effective mass stabilizes the quasi–static potential, while baryonic motion continually reorganizes the drivers of the resonant field. Long–term coherence and equilibrium. Channel 1 continually deposits curvature energy into ρeff , while Channel 2 continually redistributes matter and reshapes the sources of resonant excitation. Together these processes drive the system toward a self–organized spatiotemporal equilibrium in which ρb,1 , Φ res , ρeff , and Φ N remain mutually consistent across their respective timescales. Figure 1 illustrates these causal links: red arrows denote the Newtonian stress–advection response; black and blue arrows denote resonant excitation and slow memory accumulation. Together these links form the dual–channel architecture shown in Fig. 1. The figure compresses the resonant feed and slow accumulation into a single arrow, 24
Channel 1 Resonant channel feeding dark mass Channel 2 Newtonian channel feeding continuity ρb,1Φres ρeff GΩ rescat./prop. (∝ Genv) ρb,0, ρeff ΦNvρb,0, ρb,1 ∇2−∇ Continuity Advection long-term energy storage memory ∝(τmem) Figure 1: Two-channel gravitational feedback network. Channel 1 (top, black →blue): the oscillatory baryonic component ρb,1 drives the resonant gravitational field Φ res through the causal propagator GΩ . Multiple scattering produces an intensity gain Genv , and the resonant field deposits accumulated curvature into the slow, long-lived component ρeff over the memory time τmem = Γ −1 .Channel 2 (bottom, red): the slowly varying density ( ρb,0, ρeff ) sources the Newtonian potential Φ N , which induces bulk velocities and advects the baryonic distribution. Blue dashed arrows denote slow memory feedback; red arrows denote instantaneous Newtonian response and mass transport. Together the two channels form a closed spatiotemporal feedback loop linking baryonic motion, resonant gravitational response, and long-term curvature accumulation. 25
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 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 5.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) = ReA(x)e−iΩt,(53) 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. Throughout each coherent domain we therefore adopt the locking approximation Ω≃ωres,(54) neglecting slow amplitude variations except near steep density gradients or strong baryonic asymmetry. 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. 6.1 Kirchhoff–Sommerfeld Representation Assuming global synchronization at frequency Ω, all curvature modes carry the common temporal factor e−iΩt . Substituting the ansatz (53) into the linearized field equation (13) 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),(55) which determines the spatial profile of the phase–locked resonant field. The exact solution of Eq. (55) is the standard Kirchhoff–Sommerfeld representation: A(x) = G c2 eff Zd3x′eikΩ|x−x′| |x−x′|δρb(x′),(56) 32
with complex wavenumber kΩ=Ω+iΓ ceff , λres =2πceff Ω.(57) The imaginary part of kΩ produces a slow spatial attenuation over the damping length Ldamp =ceff Γ=ceffτmem,(58) which measures the spatial extent of curvature memory. For weak damping (Γ≪Ω) and for r≪Ldamp, one may take kΩ≃2π/λres as real. On galactic and cluster scales the characteristic size a of a coherent structure spans only a modest fraction of the resonant wavelength, a/λres ∼ 0 . 1–0 . 2. Thus the usual far–field diffraction hierarchy does not apply: the effective curvature field never samples many 2 π cycles across a structure. Instead each galactic disk, elliptical ensemble, or local speckle cell acts as a phase–locked resonant domain: a slowly varying envelope maintained by baryonic motion, curvature memory, and delayed feedback. All regions within a coherence cell share a common temporal phase and a spatial correlation scale set by λres and Ldamp. Unlike optical propagation, this regime behaves as a quasi–static, phase–coherent resonator whose geometry and memory determine the spatial profile of the field. Section 6.2 shows how these coherent envelopes translate, through long–time averaging, into the stored curvature energy that feeds the observable Newtonian channel. In this representation the source term δρb (x) in Eq. (55) 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),(59) 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. (55) 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 (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. (56) , 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. 33
This form makes clear that the effective drive of the resonant channel is determined not 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′|(60) 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′) (61) 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. (61) . 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 6.2. 6.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,(62) 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 34
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. The local energy density of the Channel 1 resonant field, expressed through the complex spatial envelope A(x) introduced in Eq. (53), is EΦ(x) = 1 8πGhc2 eff|∇A|2+ Ω2|A|2i,(63) 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.(64) 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. (64), ρeff(x)=1 8πGc2hc2 eff|∇A|2+ Ω2|A|2i,(65) 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.(66) 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. 6.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 35
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. 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.(67) Inserting this source into the Kirchhoff–Sommerfeld representation (Eq. (56) ) yields the corresponding Channel–1 resonant field δΦp(r): δΦp(x) = Gmb c2 eff eikΩr r, r =|x−xp|.(68) Equation (68) 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. (68) , giving δΦp(r) = Gmb c2 eff e−r/Ldamp j02πr λres ,(69) 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)|21/2=Gmb c2 eff e−r/Ldamp r√2j0(2πr/λres).(70) 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. (65) . 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 02πr λres . (71) 36
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 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. 6.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),(72) with each contribution given by δΦpi(x) = Gmi c2 eff eikΩri ri , ri=|x−xi|, i = 1,2.(73) 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, Ω = 37
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 into a single smooth envelope with only slight modulation. At 100 pc 2 (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.(74) Substituting Eq. (73) gives the explicit interference pattern |A|2=G2 c4 eff e−2(r1+r2)/Ldamp "m2 1 r2 1 j2 02πr1 λres +m2 2 r2 2 j2 02πr2 λres +2m1m2 r1r2 j0 2πr1 λres j0 2πr2 λres cos 2π∆r λres #, (75) 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. (65) . 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,(76) showing directly that interference between oscillatory fields modulates the local effective gravitational density. The cross–term in Eq. (75) 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. 2 Subsequent simulations (see Section 7) indicate that these separations and wavelengths are not realistic for galaxy–scale resonances, but the configuration remains pedagogically useful. 38
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. 39
Figure 3: Two–source interference (100 pc separation). Same setup as Fig. 2, but with the sources several wavelengths apart, producing strong constructive and destructive interference fringes. 40
This two–source configuration makes the phase structure of the Kirchhoff– Sommerfeld representation (Section 6.1) fully explicit. The Channel–1 driver takes the form δρb(x) = δρb(x)eiφd(x),(77) where φd (x) encodes the local phase delay of the baryonic motion relative to the global carrier e−iΩt . The propagation kernel GΩ (x , x ′ ) rephases each emitter before summation, so that A (x) is the phase–matched projection of the source onto the resonant field. For two isolated emitters this becomes transparent. Assign complex driver amplitudes δρ1eiϕ1, δρ2eiϕ2,(78) at x1and x2. Each generates a complex effective curvature field δΦp1(x)eiϕ1, δΦp2(x)eiϕ2,(79) so that the total Channel–1 envelope is A(x)=δΦp1(x)eiϕ1+δΦp2(x)eiϕ2.(80) The relative phase ∆ϕ=ϕ1−ϕ2(81) determines whether the superposition is constructive or destructive. If the sources maintain a fixed configuration or relative motion, ∆ ϕ remains constant in the global resonant frame, producing a stationary interference pattern and a time–independent composite envelope. In this discrete setting the role of the complex amplitudes is fully visible: δρ1,2eiϕ1,2 describes the phase–weighted baryonic drive, the kernel applies the additional propagation phase eikΩri , and Eq. (80) is the resulting phasor sum. This is precisely the two–element analogue of the coherent integral A(x) = G c2 eff Zd3x′GΩ(x,x′)δρb(x′) (82) introduced in Eq. (56). If many emitters satisfy ∆ϕ(x)≈const,(83) their contributions add coherently and the envelope amplitude grows approximately with the number of participants. If instead their phases are random, ei∆ϕ= 0,(84) the cross–terms cancel statistically, leaving only self–intensity contributions and producing an incoherent, speckle–like background. This phase–locking behaviour is the discrete prototype of the coherent projection mechanism governing extended baryonic systems. It forms the conceptual bridge to the collective resonant modes developed in the following sections, where the interference of many emitters shapes the large–scale resonant envelope and the emergent effective density ρeff. 41
A mode remains phase–coherent only if ω(rn)(1 −m)−Ω = 0,(108) i.e. ω(rn) = Ω 1−m.(109) Only m= 1 allows a physical solution, ω(r0) = Ω,(110) 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)θ],(111) with mode frequency ∆Ωm(r)=Ω−(1 −m)ω(r).(112) 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θ.(113) 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θ.(114) 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).(115) 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θ. (116) We parametrize the coherent m= 1 component as Σb,1(r, θ) = ϵ(r) Σ0(r)eiθ,(117) so that M1(r)=2π ϵ(r) Σ0(r), ϵ(r) = M1(r) 2πΣ0(r).(118) 48
This coherent moment M1 ( r ) is the unique long–time driver of the global m = 1 spiral resonance. In Section 7.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. 7.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 7.1.1 onto the thin disk. 3 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 (55) 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, θ),(119) 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θ,(120) 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. (118) . Substituting the expansions (120) into Eq. (119) yields, for each mode m, d2Am dr2+1 r dAm dr +k2 Ω−m2 r2Am=4πG c2 eff ρm(r).(121) 3 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. 49
The homogeneous solutions of Eq. (121) 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 B 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. (119) yields the radial Helmholtz equation for M1(r): d2M1 dr2+1 r dM1 dr +k2 Ω−1 r2M1=4πG(Γ+iΩ) c2 eff ρ1(r).(122) 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. (56) . The corresponding Green–function solution of Eq. (122) 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′#, (123) 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′#, (124) so that the two–dimensional field can be written compactly as A(r, θ) = eiθ Ar(r).(125) 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 7.1.4 to construct the three–dimensional field A ( r, θ, z ) and in Section 7.1.6 to evaluate the stored energy and effective dark–mass distribution. 50
7.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),(126) 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,(127) 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.(128) 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),(129) where the vertical scale height follows from the evanescent branch of Eq. (128) as L−1 z(r) = pℜ[−k2 z(r)].(130) The normalization in Eq. (129) 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. (123) - is preserved here, ensuring that the three–dimensional extension retains the correct far–field amplitude and outgoing–wave behavior. 51
The complete three–dimensional eigenmode is therefore A(r, θ, z) = eiθ e−|z|/Lz(r) p2Lz(r)Ar(r),(131) where Ar ( r ) is defined in Eq. (124) . Equation (131) 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. 7.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 4.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⋆ϵ, (132) 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′, θ −θ′),(133) 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θ,(134) 52
scatters according to A(1)(r, θ) = ZGΩ(r, r′, θ −θ′)A(0) r(r′)eiθ′r′dr′dθ′.(135) Performing the θ′ -integration isolates only the m = 1 harmonic of the kernel, yielding A(1)(r, θ) = eiθ A(1) r(r),(136) 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. (132) . 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. 7.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 6.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.(137) Substituting the normalized field A ( r, θ, z ) of Eq. (131) and separating variables gives |∇A|2=∂rAr 2|Z0|2+|Ar|2 r2|Z0|2+|Ar|2∂zZ0 2.(138) 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 zi.(139) 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 53
bracket approximately constant: Ω2+c2 effk2 r+m2 r2+1 L2 z≃Ω2+c2 effk2 Ω≈2 Ω2,(140) since kΩ= Ω/ceff and Γ≪Ω. Thus, to good approximation, ρeff(r, θ, z)≃2 Ω2 8πGc2|A(r, θ, z)|2.(141) 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),(142) 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.(143) 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. (142) into Eq. (141) 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.(144) The result (144) 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. 7.1.7 Numerical Results: 3D Thin–Disk Model The complete numerical results obtained from the three–dimensional thin–disk model of Section 7.1.4 are shown in Figures 4–6. All quantities and algorithms are implemented exactly as specified in the simulation code given in Appendix C. 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 54
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 55
Figure 4: Radial profiles and rotation curves. (a) Field amplitude |A ( r ) | ; (b) effective and baryonic mass densities (logarithmic); (c) rotation curves for baryonic, resonant, and total components; (d) enclosed–mass ratio Meff/Mb . Results from the 3D thin–disk Helmholtz model (Section 7.1.4); numerical parameters as defined in Appendix C. 56
Figure 5: Midplane amplitude and density maps. (a) Normalized amplitude |A| showing the toroidal resonant structure; (b) normalized phase arg ( A ) exhibiting the spiral m = 1 winding; (c) combined baryonic (red) and resonant (green) densities; (d) total normalized density on a linear scale. Spatial domain: 100 ×100 kpc. 57
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,(150) with transport damping length Ldamp,eff =ceff Γeff ,(151) where ceff is the effective curvature–wave speed. Interpretation. Equation (148) expresses the condition for statistical stationarity: the effective gravitating density generated by random stellar motion balances the density removed by transport damping. Equations (150) – (151) 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. 7.2.2 Complex Baryonic Forcing with Memory A key result of this paper 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. 4 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 3.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]. 5 4 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 χ(Ω). 5 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. 64
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 7.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),(152) which drives the forced–damped Helmholtz operator ∇2+k2 ΩA(x, t) = 4πG c2 eff ρb(x, t),(153) 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),(154) 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 Γ,(155) representing the causal damping of an individual curvature phasor. Ensemble averaging over many such phasors will later produce the transport–renormalized kernel. 65
In Eq. (154), 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 (154) 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 7.2.4 to rewrite Eq. (154) in terms of the independent phasor ensemble and to compute ensemble averages such as ⟨|A|2⟩ and ⟨ρeff⟩. 7.2.3 Four–dimensional speckle geometry of the gravitational field To extract statistical predictions from Eq. (154) , 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π,(156) 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Γ,(157) 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π2v 2Γ,(158) 66
where α is a geometric factor of order unity. The corresponding temporal extent is ∆tcell =τcoh =1 2Γ.(159) 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 ,(160) and each cell contains on average Ncell =n⋆Vcell =N⋆Vcell Vgal ,(161) 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⋆.(162) The original stellar sum in Eq. (154) 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. (154) 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. 7.2.4 Ensemble reduction: masson–speckle summation We now rewrite the forcing in Eq. (154) in terms of the statistically independent four–dimensional coherence cells derived in Section 7.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,(163) 67
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π,(164) so that each cell contributes NT independent temporal samples of its masson excitation. Intensity expansion. Taking the square modulus of Eq. (163) 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−xcG∗ Ωx−xc′,(165) 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, (166) 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|2GΩ(x−xc) 2.(167) For a Poisson stellar distribution, |Mcell|2=N2 cellm2 ⋆+NcellVar(m⋆),(168) and for Ncell ≫1 this simplifies to |Mcell|2≃M2 cell =N2 cellm2 ⋆.(169) 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.(170) It is now convenient to collect the ensemble coupling terms into a single dimensionless environmental gain factor Genv . In Section 4.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. (25) ]. In 68
the effective–medium picture of Section 4.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π2v 2Γ,(171) where Vcell = α ℓ2 ⊥ℓ∥ is the spatial coherence volume from Eq. (158) , ℓ⊥ = λres/ (2 π ) is the transverse coherence length (Eq. (156) ), and ℓ∥ = v/ (2Γ) is the longitudinal coherence length (Eq. (157) ). 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. (164), (170), and (171) into Eq. (167) finally gives D|A(x, t)|2E=4πG c2 eff 2Q 2πGenv GΩ(x−xc) 2,(172) 69
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. 7.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 .(173) 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,(174) and the associated effective damping length Ldamp,eff =ceff Γeff .(175) In the diffusive regime, Γeff ≪Γ, Ldamp,eff ≫Ldamp =ceff Γ,(176) 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 Γ,(177) 70
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 7.2.4 retain their formal structure under the replacements Γ→Γeff, Q →Qeff, Ldamp →Ldamp,eff,(178) 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π2v 2Γeff .(179) This form preserves the same microscopic structure as Eq. (171) , 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,(180) 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. (179) and (178) into the ensemble intensity (172) , 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).(181) 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. 7.2.6 Core Regularization in the Non–Transport Regime Before introducing the full energy conversion in Section 7.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 71
will be applied in the next subsection. Accordingly, the discussion below treats ρeff ( r ) as a quantity proportional to ⟨|A ( r ) |2⟩ , with the understanding that the proportionality constant will be specified later. The transport theory described in the previous section predicts a far–field profile proportional to exp ( − 2 r/Ldamp,eff ) /r2 . This expression is exact for the ensemble–averaged coherent mode of a multiply scattered curvature–phonon field in the diffusive regime; that is, at radii beyond several coherence lengths where angular isotropy and statistical decorrelation hold [ 19 , 8 ]. However, the same form cannot apply arbitrarily close to the centre. At small radii: (i) stellar scatterers are closely packed, (ii) direct (ballistic) propagation dominates over diffusion, and (iii) the local speckle field becomes dominated by the coherent mean–field curvature rather than by ensemble statistics. In this inner region, the diffusion approximation implicit in the ensemble–transport kernel Keff ( r ) ∝ 1 /r2 breaks down. As r→ 0, Eq. (181) would therefore diverge unphysically, signalling the transition from the statistical (Dyson) regime back to the coherent Foldy–Lax limit. Curvature flux must remain finite at the centre of any bound system. For a spherically symmetric, stationary envelope, the net curvature energy flux crossing a sphere of radius rmust vanish smoothly as r→0; in particular, there can be no singular point source or sink at the origin. Writing the flux schematically as F(r)∝4πr2ρeff(r),(182) regularity requires that F(r) decrease at least as fast as r3, F(r)∝r3,(r→0),(183) which implies the small–rscaling ρeff(r)∝r, (r≪ℓ∥),(184) where ℓ∥ = v/ (2Γ eff ) is the longitudinal coherence length of the ensemble–averaged field. This linear–core behaviour is therefore a natural, nonsingular choice: it regularizes the 1 /r2 divergence of the diffusive kernel while ensuring finite central flux and a smooth approach to the coherent (inner) regime. We choose a matching radius Rmatch of order the outer stellar scale radius, where the field has undergone many scatterings and the transport description becomes valid: Rmatch ∼ O(few Re)∼5 kpc for typical ellipticals.(185) The regularized envelope then takes the generic form ρeff(r)∝ r, r < Rmatch, exp[−2r/Ldamp,eff ]/r2, r ≥Rmatch. (186) Both the density and its effective flux 4 πr2ρeff are continuous at r = Rmatch by construction. 72
The inner linear law is thus not an ad hoc empirical patch but a regularized limit of the transport solution: it removes the unphysical central divergence of the diffusive kernel, keeps the central flux finite, and matches smoothly onto the exponential–over–r2transport tail. It predicts: •cored density profiles in ellipticals and clusters, •a flat circular–velocity contribution at small radii, and •a smooth transition to the exponential–over–r2transport tail. The core radius is fixed by the coherence scale ℓ∥ = v/ (2Γ eff ), while the outer halo is controlled by the transport damping length Ldamp,eff . Together, these two scales define a two–zone structure—a coherent core joined continuously to a diffusive halo—that remains fully consistent with the transport hierarchy. In what follows, we introduce the curvature–energy conversion and gradient–energy correction that complete the physical normalization of this envelope, yielding the final quantitative halo law. 7.2.7 Gradient–Energy Boost and Final Transport–Regularized Halo Law Having established the spatial morphology of the coherent curvature envelope, we now complete its physical normalization by restoring the curvature–energy relation and including the gradient–energy contribution. This step connects the ensemble field intensity ⟨|A ( r ) |2⟩ to the observable effective density ρeff ( r ) and yields the final closed halo law. The local curvature–energy density is ρeff(x) = 1 8πGc2Ω2|A(x)|2+c2 eff |∇A(x)|2,(187) where the first term represents the potential energy of the global curvature oscillation and the second its elastic curvature (gradient) energy. Defining a correlation length ℓcby ℓ−2 c:= |∇A|2 |A|2,(188) Averaging Eq. (187) and using the definition ℓ−2 c=⟨|∇A|2⟩/⟨|A|2⟩gives ρeff(r)=Ω2 8πGc2D|A(r)|2E1 + c2 eff Ω2ℓ2 c≡Ξ∇ Ω2 8πGc2D|A(r)|2E,(189) where the dimensionless gradient–energy boost is Ξ∇= 1 + 2π Ldamp,eff λres 2 ,(190) 73
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 ,(200) 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 5.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,(201) 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,(202) 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 80
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,(203) 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,(204) 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. 7.3 Clusters as Radiative Transport Extensions of Ellipticals 7.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 ,(205) 81
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 82
weak cross–coherence produces the modest reverberation gain Grev =1+Ngalη, (206) 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),(207) where Lmean =ceffτtr (208) 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),(209) yielding the statistical amplification Grev =1+Ngalη, (210) 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,(211) then the ensemble–averaged effective density obeys ρeff(r)∝Grev exp[−2r/Lmean] r2,(212) with Lmean =ceff Γtr , Grev =1+Ngalη. (213) The convolution of these overlapping exponential envelopes approaches the observed NFW form ρ(r)∝r−1(1+r/rs)−2for Lmean ∼rs. 83
Interpretation. Equation (230) 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 ,(214) 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. 7.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. (181) . At the cluster scale, the total field amplitude is the coherent sum Atot(r) = Ngal X i=1 Ai(r),(215) so that the ensemble–averaged intensity becomes |Atot(r)|2=X i⟨|Ai|2⟩+X i=j⟨AiA∗ j⟩.(216) 84
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⟩,(217) 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η),(218) identifying the dimensionless reverberation gain Grev =1+Ngalη. (219) 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. (181) . Each elliptical contributes a transport–renormalized amplitude envelope ⟨|Ai|2⟩ ∝ Keff ( | r − r i| ), including the gradient–energy boost factor Ξ ∇ defined in Eq. (190) . The cluster-scale superposition is therefore ρeff(r) = Grev Ξ∇ Ngal X i=1 Keff(|r−ri|),(220) 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 85
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],(221) 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′(222) 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 86
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′µ, (223) 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 ,(224) which yields the explicit expression Λgeom =Ngal M(elliptical) DM Z∞ 0 4πr2ρshape(r)dr .(225) 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. 87
We model this effect phenomenologically through a nonlinear amplification, ρeff,rev(r)=Grev Λgeom ρshape(r), Grev =1+Ngalη, (226) where ρshape ( r ) is the normalized 2.5D convolution result and Λ geom enforces the global mass normalization. The reverberation gain Grev thus captures the small but finite cross–coherence among the overlapping transport fields of neighbouring galaxies. Laboratory measurements in acoustics and optics typically find η∼ 10 −3 –10 −2 , corresponding to total amplification factors Grev ≃3–10 for systems with many independent scatterers. This range coincides remarkably well with the amplification required for clusters: the ratio between the summed dark–mass contribution of all ellipticals and the total mass inferred from NFW fits is of precisely the same order. In this sense, the “reverberation gain” formalism provides a natural, physically motivated bridge between the mesoscopic transport behaviour of individual galaxies and the macroscopic gravitational potential of the cluster as a whole. Lognormal Distribution of Transport Lengths In many studies of radiative, optical, and acoustic transport through disordered or heterogeneous media, it is common practice to represent variability in local scattering or damping lengths using a lognormal ensemble rather than a single mean value. This reflects the multiplicative nature of scattering paths and the intrinsic heterogeneity of local environments, which naturally produce a lognormal distribution of effective transport lengths. Such formulations are well established in wave-diffusion and multiple-scattering analyses across acoustics [ 4 , 13], optics [19, 15], and mesoscopic transport theory [8]. Following this convention, we extend the kernel to an ensemble-averaged form, Kens(r) = ZK(r;L)P(L;Lmean, σln L)dL, (227) where P ( L ) is a lognormal probability density with logarithmic dispersion σln L . This procedure reproduces the empirically observed heavy tails in intensity transport and regularizes the large-scale overlap behaviour without introducing any ad hoc smoothing. Typical acoustic or optical diffusion studies employ 0 . 3 ≲σln L≲ 1 . 0 for moderately heterogeneous media, and occasionally σln L∼ 2 in strongly scattering or porous systems where transport paths vary over several orders of magnitude [ 19 , 8 ]. In the present context, a value of σln L = 2 . 5 yields a physically consistent and observationally well-matched profile: it preserves the correct NFW-like slope at large radii while capturing the required convex curvature near the core. The long-tail outliers in this distribution carry direct physical meaning—they represent rare instances in which a neighbouring galaxy’s transport tail penetrates a baryonic region, re-initiating a secondary radiative-transport process with an extended effective damping length. 88
While laboratory studies typically infer 0 . 5 ≲σln L≲ 2 . 0 for single-phase diffusive media [ 4 , 13 ], strongly resonant or dual-channel systems exhibit effective spreads approaching σln L∼ 2–3 [ 8 ], fully consistent with the broad heterogeneity expected in the cluster-scale double-transport regime considered here. Effective truncation of the ensemble. In practice, radiative-transport models rarely integrate the lognormal distribution to infinity. Rather, they approximate the ensemble average by sampling over a finite range of L/Lmean , typically within a few standard deviations of the logarithmic mean. This suppresses the statistical weight of extreme outliers without altering the intrinsic shape of K(r;L). Such effective truncation is standard practice in acoustic and optical diffusion analyses, where the physical coherence length of the medium imposes a natural limit on valid transport paths [ 4 , 13 , 8 ]. Our five-point sampling across ± 2 σln L follows this convention: it retains the full lognormal curvature while strongly suppressing contributions from exceedingly long paths ( L≫Lmean ). This preserves the physical broadness of the ensemble and the fidelity of the kernel’s form, exactly as in established radiative and ultrasonic transport treatments. 7.3.3 Numerical Results The numerical implementation follows directly from the formalism developed above. All calculations employ the reduced 2.5D convolution to evaluate the transport–overlap density ρeff ( r ), followed by application of the geometric normalization Λ geom to enforce total mass conservation, and the reverberation gain Grev = 1 + Ngalη to incorporate partial coherence between overlapping tails. The complete algorithm and parameter definitions are given in Appendix E, together with the Python code used to generate the figures. For a fiducial 10 15 M⊙ cluster, we adopt Ngal = 300 ellipticals, Lmean = 636 kpc , σln L = 2 . 5, and η = 8 . 3 × 10 −3 , corresponding to Grev ≃ 3 . 5. The resulting one-dimensional transport profiles are presented in Figure 10, which compares the model’s effective density ρeff ( r ) and enclosed mass Meff ( r ) with a standard NFW benchmark. All quantities were computed on a logarithmic grid spanning 0 . 2 kpc ≤r≤ 3 Mpc , using the 2.5D convolution for ρshape ( r ) and normalized by Λ geom as defined in Eq. (225) to preserve total enclosed mass. The reference NFW halo corresponds to M200 = 1015 M⊙and c200 = 4. Panels (a) and (b) of Figure 10 show the effective density in linear and logarithmic scales, respectively. The transport–overlap profile reproduces the NFW slope and amplitude across three orders of magnitude in radius, maintaining close agreement from ∼ 30 kpc to 3 Mpc . The slightly increased concavity of ρeff ( r ) at small radii ( r≲ 50 kpc ) possibly arises from the omission of the central brightest cluster galaxy (BCG), whose extended stellar envelope would naturally flatten the inner slope. The shallow dip near r→ 0 originates from the discrete angular integration and is purely numerical. Panels (c) and (d) display the corresponding enclosed-mass profiles. The total mass derived from the transport model follows the NFW prediction almost 89
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. 9 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 3.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,(237) 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. (237) in isolation. 9.1 Three Distinct Regimes for Homogeneous Clouds A homogeneous cloud may arise in three physically distinct settings: [label=(1)] 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.(238) 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.(239) 96
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 ,(240) but the scattering strength remains σ(x)≡0,(241) 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. 9.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 9.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.(242) 2. Embedded but free-running: Without phase–locking, the cloud does not couple to the resonant field: Φres = 0.(243) 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.(244) The reverberation gain vanishes, Grev =1+Nη −→ 1, η = 0,(245) and the stored curvature energy remains at the single–pass level: ρ(hom) eff ∼1 Nρ(inhom) eff ,(246) 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. 97
9.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 3.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′,(247) 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).(248) 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. 10 Emergence of Multiple Scattering from a Weakly Inhomogeneous Continuous Medium 10.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),(249) 98
consisting of spatially separated stellar or galactic constituents. Within a coherent domain these constituents oscillate at a common curvature frequency Ω2∝G¯ρb,(250) 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),(251) 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′,(252) where the contrast–defined scattering strength is σ(x)∝δρb(x).(253) Applying the Foldy–Lax operator to both sides of Eq. (252) produces the renormalised effective–medium wave equation, ∇2Φres +1 c2 eff −∂2 t−2Γ ∂t+ Ω2Φres =4πG c2 eff ρb,1,(254) 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. (252) . The cycle–averaged stored curvature energy defines the effective density, ρeff(x)∝ ⟨|Φres(x)|2⟩.(255) Crucially, nothing in Eqs. (250) – (255) 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,(256) and therefore realises Regime 1 of Section 9.2, namely an isolated, homogeneous system with σ(x)≡0, Grev = 1,Φres = 0.(257) 99
In this state no Foldy–Lax feedback operates, no effective self–energy can form, and the perturbation theory of Section 3.4 does not apply. There is therefore ρeff(x, t) = 0 (early homogeneous era).(258) 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. 10.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),(259) In the early plasma, the coarse–grained density takes the form ¯ρb(x, t) = ¯ρb(t) [1 + δ(x, t)],|δ| ≪ 1,(260) so that Ω(x, t) = Ω0(t)1 + δ(x, t) 2+O(δ2).(261) 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).(262) Thus the early medium lies strictly in the non-scattering regime of Regime 1. 10.3 Weak Random Potential and the Ballistic Regime Formally inserting (261) into the resonant operator yields a Helmholtz equation of the form ∇2+k2 0+V(x, t)Φres =S, k0=Ω0(t) ceff ,(263) with weak random potential V(x, t)=2k0 δΩ(x, t) ceff +O(δ2)∝δ(x, t).(264) 100
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).(265) In the primordial plasma one has σV∼δ∼10−5, giving k0ℓs≫1,(266) 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 3.4: since no effective medium exists, the perturbative theory there does not apply at early times. 10.4 Nonlinear Growth Toward the Transport Threshold As gravitational instability amplifies density fluctuations, δ→ O(1), σV→ O(1), ℓc→galactic scales,(267) the combination in Eq. (265) grows by many orders of magnitude. There exists a cosmological epoch t∗defined implicitly by k0(t∗)ℓs(t∗)∼1,(268) 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. 10.5 Emergent Resonant Patches and Effective Discreteness When δ∼ 1, the medium fragments into nonlinear overdensities. Within each patch, Ω(x)≃Ωpatch ∝√ρpatch,(269) so curvature energy becomes locally trapped. The cycle–averaged effective density grows as ρeff(x)∝|Φres|2,(270) 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 4–7.3.4. This mechanism closely parallels modulational instability in plasmas, where spatial variations in ωp ( x ) ∝pne(x) produce a refractive index landscape that 101
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. 10.6 Linear Limit and Connection to the CMB Power Spectrum The developments of Sections 4–7.3.4 demonstrate that multiple scattering, renormalised propagation, and ρeff arise only once the contrast σ (x) becomes nonzero. In the primordial plasma, however, Eqs. (256) – (258) imply that the Universe lies strictly in Regime 1, so that σ(x) = 0,Φres = 0,(Ω,Γ, ceff) undefined,(271) 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. (259) — remains a well-defined linear response field. In Regime 1 all renormalised parameters revert to their bare values: ceff →c, Γ→0,Ω→Ω0(t),(272) and the effective–medium PDE (254) 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),(273) which is simply the linear curvature response to the baryon overdensity. Fourier transforming Eq. (273) gives ¨ Φres(k, t) + Ω2 0(t) Φres(k, t)=4πG ¯ρb(t)δb(k, t),(274) 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. (274) adds no nonlinear gravitational effect. 102
It is important to emphasise that Eq. (274) 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) (275) 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(276) 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.(277) The comoving acoustic scale corresponds to k∼0.01–0.1 Mpc−1, i.e. k c ∼10−13–10−12 s−1.(278) Thus ϵ(t∗) = Ω2 0 k2c2∼10−8–10−10,(279) many orders of magnitude below unity. Therefore: ϵΦres ≪ΦN,(recombination era),(280) so the curvature–phonon contribution is far too small to affect the acoustic oscillations or the resulting CMB temperature and polarisation spectra. It is automatically consistent with Planck-precision constraints and cannot lead to double counting of metric perturbations. In summary, the Regime 1 limit of the model maps seamlessly onto the standard linear CMB hierarchy. The curvature–phonon field acts only as a tiny, strictly linear response driven by the same baryonic fluctuations already present in Φ N , and its contribution is O (10 −8 ) relative to the dominant gravitational potential. Only once δ→ O (1) does σ (x) become nonzero and the renormalised effective medium of Sections 4–7.3.4 begin to form. 103
11 Conclusion This work has shown that a small but physically essential feature of linearised general relativity—the full retarded Green–function response— combined with the ordinary damped–oscillator susceptibility of baryonic matter, yields a self–consistent and observationally successful extension of weak–field gravitational dynamics. No new fields, particles, or modifications of GR are introduced. Once the usual quasistatic truncation of the retarded kernel is avoided, the delayed interaction between baryons and curvature generates a collective oscillatory field whose cycle–averaged energy behaves as an effective dark–mass density. The resulting causal–resonant framework produces a renormalised response characterised by three effective parameters (Ω , Γ , ceff ) determined by baryonic susceptibility and multiple scattering. A single evolution equation for the collective curvature field reproduces the observed effective mass distributions of spiral galaxies, elliptical galaxies, clusters, and filaments. Coherent motion in structured systems drives constructive rescattering and large–scale organised modes, while homogeneous or weakly perturbed environments provide no spatial contrast and therefore remain resonantly silent. This dichotomy accounts naturally for flat rotation curves, exponential elliptical envelopes, cluster–scale halos, filamentary coherence, and the absence of resonant effects in the early Universe and in the CMB. The underlying mechanism is unified and minimal: time–varying baryonic structure sources curvature oscillations through its susceptibility; multiple rescattering shapes their large–scale behaviour; and the finite memory of the retarded kernel allows a slow accumulation of cycle–averaged curvature energy. This stored component enters the Poisson channel as an effective mass density while the oscillatory field remains locally silent, ensuring full consistency with laboratory, Solar–System, and stellar tests of gravity. In summary, the causal–resonant framework reveals that the observed dark–mass phenomena need not signal new gravitational laws or unseen particles. They emerge naturally from the interplay between the retarded response of general relativity and the ordinary dynamical behaviour of baryonic matter. The large–scale gravitational structure of the Universe—its halos, clusters, and filaments—can thus be understood as a manifestation of the persistent memory and collective resonant dynamics already present within weak–field GR. 12 Dedication and Acknowledgement I would like to dedicate this paper to the late Prof. John T. Sheridan, my former PhD supervisor and friend, who introduced me to wave optics and is sadly missed. I acknowledge with gratitude the work of Prof. Joseph Goodman, whose seminal texts Fourier Optics and Statistical Optics have remained page-turners for me for 104
more than twenty years. His contributions underpin much of my understanding of wave theory and coherence. I am indebted to the work of Sheng, Derode, Roux, Fink, Weaver, van Tiggelen, and others in wave transport and multiple-scattering theory—an area I was not very familiar with when I began this project, yet essential for deriving the behaviours of elliptical galaxies and clusters. I sincerely hope I have not misinterpreted their theories. I gratefully acknowledge the support of Research Ireland / Science Foundation Ireland and the European Union, whose funding over the past twenty years enabled my research in optics and allowed me to deepen my understanding of physics. Last and most importantly, I thank my wife Karen and my daughters for their unwavering support as I wrestled with this theory over these past several months. A Appendix: Effective–Medium Derivation For completeness we sketch the standard weak–scattering effective–medium reduction leading to the homogenised resonant equation used in the main text. No modification of general relativity is introduced; all nonlocality arises from the ordinary retarded response of the linearised gravitational field and from the rescattering of weak perturbations by baryonic inhomogeneities. Retarded solution and decomposition of the source In the weak–field limit the gravitational potential satisfies Φ(x) = −4πGZGret(x−x′)ρb(x′)d4x′,(A.1) where Gret is the standard retarded Green function in a slowly varying background. Decomposing the density into a slowly varying part and a coherent oscillatory part, ρb(x, t) = ρb,0(x)+ρb,1(x, t),(A.2) we isolate the oscillatory response Φres driven by ρb,1. Multiple scattering and Dyson renormalisation Treating the baryons as weak gravitational scatterers with linear frequency–dependent susceptibility χ⋆ ( ω ), the ensemble renormalises the propagator through the usual Dyson series, G−1 eff (ω, k) = G−1 0(ω, k)−Σ(ω, k),(A.3) where G0 is the free weak–field propagator and Σ is the self–energy describing coherent rescattering by the inhomogeneous baryonic distribution. No new dynamical fields appear: Σ is determined fully by the statistics and susceptibility of the scatterers. 105
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 112
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’) 113
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 114
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) 115
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() 116
Listing 1: Python code for the full 3D Helmholtz resonant–halo simulation. 117
D 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 7.2.8. The simulation realizes the closed transport–law curvature–phonon halo derived analytically in Section 7.2.4 onwards, including the multiplicative gradient–energy contribution of Section 7.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.(D.1) The elliptical transport closure is governed by four control parameters: •the transport damping length Ldamp,eff =ceff Γeff ,(D.2) fixed directly from the observed halo extent; •the resonant curvature wavelength λres,(D.3) which sets both the carrier frequency Ω = 2πceff λres (D.4) and the effective quality factor Qeff =Ω 2Γeff ; (D.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. 118
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. (D.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 7.2.7. For numerical regularity, the small–radius linear core ⟨ρeff(r)⟩∝ris enforced for r < Rmatch (see Section 7.2.6). For numerical consistency and physical regularity, a small–radius linear core, ⟨ρeff(r)⟩∝r(r < Rmatch),(D.7) is enforced as in Section 7.2.6. All symbols in Eqs. (D.1) – (D.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: 119
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 papers notation for 1:1 correspondence 120
74 # between analytic expressions and numerical quantities. 75 # ===================================================================== 76 77 78 79 # ===================================================================== 80 # Elliptical Galaxy Transport-Law Statistical Resonant Halo 81 # With Multiple-Scattering Effective Damping + Gradient Boost 82 # ===================================================================== 83 84 import numpy as np 85 import matplotlib.pyplot as plt 86 87 # ---------- Constants ---------- 88 G = 6.67430e-11 89 c = 2.99792458e8 90 Msun = 1.98847e30 91 pc = 3.085677581491367e16 92 kpc = 1.0e3 * pc 93 Mpc = 1.0e6 * pc 94 pi = np.pi 95 96 # ============================================================ 97 # TRANSPORT MEDIUM (effective) PARAMETERS [DATA-DRIVEN] 98 # ============================================================ 99 100 L_damp_eff = 636.0 * kpc # [m] SET by observed halo e-folding 101 c_eff = 2.0e5 # [m/s] TUNABLE within plausible range 102 Gamma_eff = c_eff / L_damp_eff # [1/s] 103 tau_coh_eff = 1/(2 * Gamma_eff) # [s] 104 105 print("\n--- Transport Medium (effective) ---") 106 print(f"L_damp_eff = {L_damp_eff/kpc:.1f} kpc [SET]") 107 print(f"c_eff = {c_eff:.3e} m/s [TUNE]") 108 print(f"Gamma_eff = {Gamma_eff:.3e} s^-1 [IMPLIED]") 109 print(f"tau_coh_eff = {tau_coh_eff/3.15576e13:.2f} Myr") 110 111 # ============================================================ 112 # RESONANT MICROPHYSICS 113 # ============================================================ 114 115 lambda_res = 27.0 * pc # [m] 116 Omega = 2*np.pi*c_eff/lambda_res # [1/s] 117 Q_eff = Omega/(2*Gamma_eff) # dimensionless 118 Xi_nabla = 1.0 + (2.0*pi*L_damp_eff / lambda_res)**2 # gradient-energy ,→boost 119 120 print("\n--- Resonant Microphysics ---") 121 print(f"lambda_res = {lambda_res/pc:.2f} pc") 122 print(f"Omega = {Omega:.3e} s^-1") 123 print(f"Q_eff = {Q_eff:.3e}") 124 print(f"Xi_nabla = {Xi_nabla:.3e}") 125 126 # ============================================================ 121
E 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 7.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,(E.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′µ, (E.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,(E.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 ,(E.4) and scaled by the nonlinear reverberation gain Grev =1+Ngalη, (E.5) to yield the final effective density ρeff(r)=Grev Λgeom Sbase ρshape(r).(E.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. 128
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 # --------------------------------------------------------------------- 129
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 # --------------------------------------------------------------------- 130
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) 131
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 132
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) 133
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 134
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)", 135
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. 136
[6] Leslie L Foldy. The multiple scattering of waves. i. general theory of isotropic scattering by randomly distributed scatterers. Physical review, 67(3-4):107, 1945. [7] Melvin Lax. Multiple scattering of waves. ii. the effective field in dense systems. Physical Review, 85(4):621, 1952. [8] Ping Sheng and Bart van Tiggelen. Introduction to wave scattering, localization and mesoscopic phenomena., 2007. [9] Arnaud Derode, Philippe Roux, and Mathias Fink. Robust acoustic time reversal with high-order multiple scattering. Physical review letters, 75(23):4206, 1995. [10] Anthony E Siegman. Lasers. University science books, 1986. [11] Joseph W Goodman. Introduction to Fourier optics. Roberts and Company publishers, 2005. [12] Alison M Yao and Miles J Padgett. Orbital angular momentum: origins, behavior and applications. Advances in optics and photonics, 3(2):161–204, 2011. [13] Richard L Weaver. Diffusivity of ultrasound in polycrystals. Journal of the Mechanics and Physics of Solids, 38(1):55–86, 1990. [14] Richard L Weaver. Anderson localization of ultrasound. Wave motion, 12(2):129–142, 1990. [15] Joseph W Goodman. Statistical optics. John Wiley & Sons, 2015. [16] Richard L Weaver and Oleg I Lobkis. Ultrasonics without a source: Thermal fluctuation correlations at mhz frequencies. Physical Review Letters, 87(13):134301, 2001. [17] Richard L Weaver and Oleg I Lobkis. Diffuse fields in open systems and the emergence of the green’s function (l). The Journal of the Acoustical Society of America, 116(5):2731–2734, 2004. [18] Eric Larose, Ludovic Margerin, Arnaud Derode, Bart van Tiggelen, Michel Campillo, Nikolai Shapiro, Anne Paul, Laurent Stehly, and Mickael Tanter. Correlation of random wavefields: An interdisciplinary review. Geophysics, 71(4):SI11–SI21, 2006. [19] MCW van van Rossum and Th M Nieuwenhuizen. Multiple scattering of classical waves: microscopy, mesoscopy, and diffusion. Reviews of Modern Physics, 71(1):313, 1999. 137