Eddington Inversion Formula

The Eddington inversion formula (or Eddington's formula) is a foundational integral inversion theorem in stellar and galactic dynamics—derived by British astrophysicist Sir Arthur Stanley Eddington in 1916—that reconstructs the unique ergodic (isotropic) six-dimensional phase-space distribution function f(ℰ) of a spherically symmetric system of stars or dark matter particles directly from its three-dimensional spatial density profile ρ(r) and Newtonian gravitational potential Φ(r). Whereas integrating the collisionless Boltzmann equation over velocity space yields the Jeans equations for the second velocity moment σ_r²(r) without ensuring that a non-negative distribution of stellar orbits exists in phase space, Eddington noticed that the integral expressing ρ(r) in terms of f(ℰ) has the mathematical form of an Abel integral equation. Inverting this Abel integral produces an explicit formula for f(ℰ) that tests whether any proposed spherical galaxy or halo model—such as the Plummer model, Hernquist profile, Jaffe profile, Dehnen profile, or Einasto profile—can be supported in collisionless equilibrium by an isotropic velocity distribution (f(ℰ) ≥ 0).

Consider a steady-state, spherically symmetric collisionless system of stars or dark-matter particles moving in a smooth gravitational potential Φ(r) that approaches Φ(∞) = 0 at infinite radius (or Φ(r_t) at a finite tidal boundary r_t, as in the King profile). To work with positive quantities for bound orbits, dynamicists define the relative gravitational potential Ψ(r) and the relative binding energy per unit mass ℰ:

Ψ(r) ≡ −Φ(r) + Φ₀ > 0 and ℰ ≡ −E + Φ₀ = Ψ(r) − ½v²

where Φ₀ = Φ(∞) = 0 for an unbounded system so that bound orbits satisfy ℰ > 0 and unbound escaping particles have ℰ ≤ 0 (where f(ℰ) = 0). By Jeans's theorem for the collisionless Boltzmann equation, a spherically symmetric system has an isotropic velocity dispersion tensor everywhere (σ_r² = σ_θ² = σ_ϕ², or Binney anisotropy parameter β = 0) if and only if its phase-space distribution function depends on the phase-space coordinates (x, v) through the relative energy ℰ alone: f = f(ℰ).

Integrating f(ℰ) over all three velocity components in spherical velocity coordinates (d³v = 4πv² dv) up to the local escape speed v_esc(r) = [2Ψ(r)]^1/2 yields the three-dimensional mass density ρ(r) at radius r:

ρ(r) = ∫₀^v_esc(r) 4πv² f(Ψ(r) − ½v²) dv

Changing the integration variable from speed v to relative energy ℰ = Ψ − ½v² (so that v = [2(Ψ − ℰ)]^1/2 and v dv = −dℰ), and noting that Ψ(r) is a monotonic decreasing function of radius r (dΨ/dr = −GM(r)/r² < 0) so that ρ(r) can be regarded as a single-valued function ρ(Ψ) of the potential Ψ, transforms the density integral into a Riemann–Liouville fractional integral of order 3/2:

ρ(Ψ) = 4π(2)^1/2 ∫₀^Ψ f(ℰ) (Ψ − ℰ)^1/2 dℰ

To solve this integral equation for the unknown function f(ℰ) given a known density–potential relation ρ(Ψ), Eddington (1916) differentiated both sides once with respect to Ψ using the Leibniz rule (noting that the integrand vanishes at the upper limit ℰ = Ψ):

dρ/dΨ = 2π(2)^1/2 ∫₀^Ψ f(ℰ) (Ψ − ℰ)^−1/2 dℰ

This equation has the canonical form of Abel's integral equation, g(Ψ) = ∫₀^Ψ f(ℰ)(Ψ − ℰ)^−α dℰ with exponent α = 1/2 (solved by Niels Henrik Abel in 1823). Multiplying both sides by (ℰ′ − Ψ)^−1/2, integrating with respect to Ψ from 0 to ℰ′, reversing the order of integration via Dirichlet's formula (since ∫_ℰ^ℰ′ (ℰ′ − Ψ)^−1/2(Ψ − ℰ)^−1/2 dΨ = π), and differentiating with respect to ℰ′ yields the Eddington inversion formula:

f(ℰ) = [1 / (8^1/2π²)] d/dℰ ∫₀^ℰ (ℰ − Ψ)^−1/2 (dρ/dΨ) dΨ

When the boundary derivative (dρ/dΨ)_Ψ=0 vanishes at infinite radius (as holds for all finite-mass galaxy models where ρ(r) falls off faster than r⁻¹ at r → ∞), integrating by parts and differentiating under the integral sign yields the second-derivative form of the Eddington inversion formula:

f(ℰ) = [1 / (8^1/2π²)] { ∫₀^ℰ (ℰ − Ψ)^−1/2 (d²ρ/dΨ²) dΨ + ℰ^−1/2(dρ/dΨ)_Ψ=0 }

In numerical calculations where ρ(r) and Ψ(r) are given as functions of radius r rather than as a closed-form composite function ρ(Ψ), changing the integration variable from Ψ to r (where Ψ(r_ℰ) = ℰ defines the apocentric radius r_ℰ of a radial orbit of energy ℰ) avoids inverting Ψ(r) algebraically:

∫₀^ℰ (ℰ − Ψ)^−1/2 (dρ/dΨ) dΨ = −∫_rℰ^∞ [ℰ − Ψ(r)]^−1/2 (dρ/dr) dr

Because f(ℰ) represents a number or mass density of particles in six-dimensional phase space, a valid isotropic model must satisfy the phase-space non-negativity condition f(ℰ) ≥ 0 for all 0 ≤ ℰ ≤ Ψ(0). If ρ(Ψ) is a monotonically increasing, convex function of Ψ (dρ/dΨ ≥ 0 and d²ρ/dΨ² ≥ 0 everywhere), then f(ℰ) ≥ 0 holds globally; furthermore, by Antonov's Second Law (1962), any spherical isotropic stellar system with df(ℰ)/dℰ > 0 (phase-space density increasing with binding energy ℰ, or df(E)/dE < 0) is stable against all radial and non-radial linear perturbation modes. Similarly, for a massless tracer population ν(r) ∝ r^−γ orbiting in the Keplerian potential Ψ(r) = GM_BH/r of a central supermassive black hole, Eddington's formula gives f(ℰ) ∝ ℰ^γ−3/2/Γ(γ − 1/2), proving that an isotropic tracer cusp around a point mass cannot be shallower than γ = 1/2.

In 1979 and 1985, Leonid Osipkov and David Merritt generalized Eddington's inversion to radially anisotropic spherical distribution functions f(Q) that depend on energy ℰ and angular momentum L through Q ≡ ℰ − L²/(2r_a²) (producing outer radial anisotropy β(r) = r²/(r² + r_a²) characteristic of incomplete violent relaxation). Defining the augmented density ρ_Q(r) ≡ (1 + r²/r_a²)ρ(r) as a function of Ψ(r) recovers f(Q) via the same Abel integral transform:

f(Q) = [1 / (8^1/2π²)] d/dQ ∫₀^Q (Q − Ψ)^−1/2 (dρ_Q/dΨ) dΨ

Why is the isotropic distribution function f(ℰ) unique for a given density ρ(r)?

In a spherical gravitational potential Ψ(r) that decreases monotonically with radius, the relationship between ρ(Ψ) and f(ℰ) is an Abel Volterra integral equation of the first kind. Because the Abel integral operator is invertible one-to-one, there is at most one isotropic function f(ℰ) of energy alone that reproduces a given three-dimensional density profile ρ(r) in a given potential Ψ(r).

What does it mean when the Eddington inversion formula gives a negative f(ℰ)?

If evaluating Eddington's integral for a proposed density profile ρ(r) and potential Ψ(r) yields f(ℰ) < 0 over any interval of binding energy ℰ, then no isotropic equilibrium system with that density profile can exist (because a negative number of stars in a phase-space cell is impossible). Such a density profile can only be supported in equilibrium if the stellar orbits are tangentially biased (β < 0).

How is Eddington inversion used in dark-matter direct-detection and N-body simulations?

In N-body simulations of galaxies and star clusters, sampling initial particle speeds from the Eddington distribution f(ℰ) ensures that the initial halo starts in Vlasov–Poisson equilibrium without transient oscillations. In terrestrial dark matter direct-detection experiments, Eddington inversion of the Milky Way's Einasto profile or Burkert profile provides a self-consistent laboratory speed distribution f(v).

1. Eddington, A. S. (1916). "The distribution of stars in globular clusters." Monthly Notices of the Royal Astronomical Society, 76(7), 572–585. https://doi.org/10.1093/mnras/76.7.572
2. Merritt, D. (1985). "Spherical stellar systems with spheroidal velocity distributions." The Astronomical Journal, 90, 1027–1037. https://doi.org/10.1086/113810
3. Cuddeford, P. (1991). "An analytic inversion for anisotropic spherical galaxies." Monthly Notices of the Royal Astronomical Society, 253(3), 414–426. https://doi.org/10.1093/mnras/253.3.414
4. Lacroix, T., Stref, M., & Lavalle, J. (2018). "Anatomy of Eddington-like inversion methods for dark matter searches in the Galaxy." Journal of Cosmology and Astroparticle Physics, 2018(09), 040. https://doi.org/10.1088/1475-7516/2018/09/040
5. Binney, J., & Tremaine, S. (2008). Galactic Dynamics (2nd ed.). Princeton University Press. https://doi.org/10.1515/9781400828722