Gravitational Softening and Two-Body Heating

A galaxy or a dark-matter halo is, to superb approximation, a collisionless system: its evolution obeys the collisionless Boltzmann equation for a smooth distribution function $f(\mathbf{x},\mathbf{v},t)$ , with the potential sourced self-consistently by $\rho=\int f\,\mathrm{d}^3v$ . The physical relaxation time from star–star encounters exceeds the age of the Universe by many orders of magnitude, so a real galaxy never feels its own graininess — it flows like a fluid in phase space. But no computer integrates a continuum. An $N$ -body code samples $f$ with $N$ macro-particles, each of mass $m=M/N$ that can be $10^{5}$ to $10^{67}$ times a real star or WIMP. That coarse sampling manufactures precisely the discreteness the real system lacks, and with it two coupled artifacts. The first is the singular close encounter: two point masses passing near each other feel a force $\propto 1/r^{2}$ that diverges, forcing microscopic timesteps and injecting spurious energy. Gravitational softening cures it by replacing the point-mass Green's function with a smoothed kernel of scale $\epsilon$ , capping the peak force at $\sim Gm/\epsilon^{2}$ at the price of a controlled bias below $\epsilon$ . The second is two-body heating: the residual finite- $N$ force fluctuations diffuse orbital energies exactly as genuine encounters do, but far too fast, on $t_{\mathrm{relax}}\propto N/\ln N$ — flattening cusps into spurious cores, puffing dense clumps, and evaporating poorly sampled subhaloes. Together they fix the domain of validity of any collisionless simulation, and controlling them is what precision cosmology stands on.

Key takeaways

An $N$ -body code represents a smooth phase-space distribution with $N$ macro-particles far heavier than any real star; this injects discreteness — spurious encounters — the true collisionless system lacks.
Softening replaces the Newtonian potential with a smoothed kernel $\Phi\propto -1/\sqrt{r^{2}+\epsilon^{2}}$ (Plummer), capping the force at $F_{\max}=2Gm/(3\sqrt{3}\,\epsilon^{2})$ and trading bias below $\epsilon$ against noise and short timesteps.
Two-body heating diffuses orbital energies on $t_{\mathrm{relax}}\simeq(0.1\,N/\ln N)\,t_{\mathrm{cross}}$ , so cores puff, cusps flatten, and small subhaloes evaporate faster than in the real limit.
Convergence to the collisionless Boltzmann solution demands $N\to\infty$ and $\epsilon\to 0$ together; the Power radius marks where heating has spoiled the inner profile within a Hubble time.
Softening length and finite $N$ are the two dominant systematics of collisionless dynamics — they separate a real dark-matter core from a numerical one.
Key idea: an $N$ -body code integrates a noisy, granular caricature of a smooth collisionless system; softening removes the singular short-range force (adding a bias below $\epsilon$ ) while finite $N$ leaves fluctuations that heat orbits on $t_{\mathrm{relax}}\sim(0.1\,N/\ln N)\,t_{\mathrm{cross}}$ — so convergence needs $N\to\infty$ and $\epsilon\to0$ together.

The core insight

An $N$ -body code integrates a noisy, granular caricature of a smooth collisionless system. Softening removes the singular short-range force (introducing a controlled bias below $\epsilon$ ), while finite $N$ leaves residual fluctuations that drive two-body relaxation and heat orbits on $t_{\mathrm{relax}}\sim(0.1\,N/\ln N)\,t_{\mathrm{cross}}$ . Convergence to the collisionless limit demands $N\to\infty$ and $\epsilon\to0$ together.

Where this lives today

Softening and two-body heating set the precision of every modern structure-formation result. The central live question is the cusp–core problem: dark-matter-only $\Lambda$ CDM haloes form cuspy NFW inner profiles, yet many dwarf galaxies appear cored — possibly signalling self-interacting or warm dark matter. Distinguishing a physical core from a numerical one is impossible without the Power radius: a genuine core must persist as $N$ grows and $r_{\mathrm{conv}}$ shrinks, whereas a heated one retreats with it. The same accounting governs subhalo abundance in near-field cosmology and dark-matter detection — poorly resolved subhaloes evaporate artificially, so counts converge only above thousands of bound particles. Flagship simulations (Aquarius, Via Lactea, and billion-to-trillion-particle GADGET, AREPO, Abacus runs) are engineered around these limits. At the opposite pole sit collisional systems — globular clusters, galactic nuclei — where relaxation is the signal, studied with direct integrators that refuse to soften real encounters. The single clock $t_{\mathrm{relax}}\propto N/\ln N$ is read two opposite ways: the physics you want, or the artifact you fear.

The physical picture

Replace the smooth mass of a galaxy — say $10^{11}$ stars — by only $10^{6}$ heavy super-stars. The mean field is nearly right, because it is an integral over the smooth density and integrals forgive sampling. But each super-star tugs on its neighbours a hundred thousand times harder than a real star would, and those oversized tugs are pure numerical noise. Softening blurs each super-star into a fuzzy ball of radius $\epsilon$ so no single passage delivers a delta-function kick; heating is the surviving fuzzy tugs accumulating and randomly stirring the orbits. The first figure shows what softening does to the force law — capping the divergent $1/r^{2}$ and turning the deep interior into a harmless harmonic spring.

Figure

Softened versus Newtonian force from one macro-particle in units of Gm/ε². The Plummer force F = Gm r/(r²+ε²)^{3/2} peaks at F_max = 2Gm/(3√3 ε²) ≈ 0.385 Gm/ε² at r = ε/√2, then rejoins the Newtonian 1/r² beyond a few ε; below ε it is harmonic, F ≈ (Gm/ε³) r, so no divergent kicks.

Two-body relaxation is a random walk in energy. One distant encounter changes a particle's velocity by a tiny $\delta v_\perp\approx 2Gm/(bv)$ ; the cumulative effect of the many weak, far encounters — which dominate the rare close ones and give the Coulomb logarithm $\ln\Lambda$ — makes $\langle(\Delta v)^{2}\rangle$ grow linearly with time. Once that diffusion has rearranged the orbital energies inside some radius, the simulation there has forgotten its collisionless initial conditions. Convergence is therefore not one number but a resolved region: outside the Power radius the density is trustworthy; inside, heating has taken over and the cusp relaxes to a spurious core, as the second figure shows.

Figure

Inner density profile of a simulated halo on log–log axes. Outside the convergence radius r_conv the measured profile tracks the true cuspy NFW slope ρ ∝ 1/r; inside, where the enclosed relaxation time falls to t_circ, two-body heating flattens it into an artificial core. Raising N pushes r_conv inward and recovers more of the cusp.

The core equations

The softened potential of a single macro-particle is the starting point. The Plummer kernel is the simplest choice that is smooth at the origin, exact at large radius, and analytically differentiable.

\Phi(r) = -\frac{Gm}{\sqrt{r^{2}+\epsilon^{2}}} \qquad (1)

where: $\Phi$ — potential of one macro-particle; $G$ — Newton's constant; $m=M/N$ — macro-particle mass; $r$ — separation; $\epsilon$ — softening length. At $r\gg\epsilon$ it is Newtonian; at $r\to0$ it flattens to $-Gm/\epsilon$ — no singular well. Differentiating gives the force the code sums — and crucially it no longer diverges: it rises from zero, peaks, and rejoins the Keplerian tail.

F(r) = -\frac{\partial\Phi}{\partial r} = -\frac{Gm\,r}{(r^{2}+\epsilon^{2})^{3/2}} \qquad (2)

where: $F$ — radial force per unit mass; the numerator $r$ kills the divergence as $r\to0$ . The peak is at $r=\epsilon/\sqrt{2}$ with $F_{\max}=2Gm/(3\sqrt{3}\,\epsilon^{2})\approx0.385\,Gm/\epsilon^{2}$ — the largest kick any pair can exchange, which sets the shortest timestep. The dynamical consequence of finite $N$ is the two-body relaxation time — the time for weak encounters to change a star's orbital energy by order itself.

t_{\mathrm{relax}} = \frac{N}{8\ln\Lambda}\,t_{\mathrm{cross}} \simeq \frac{0.1\,N}{\ln N}\,t_{\mathrm{cross}} \qquad (3)

where: $N$ — particle number; $\ln\Lambda$ — Coulomb logarithm, $\Lambda\approx b_{\max}/b_{\min}\sim N$ ; $t_{\mathrm{cross}}=R/v$ — crossing time. Relaxation takes $\sim N/(8\ln N)$ crossings, so a coarser run heats faster. The Coulomb logarithm counts the decades of impact parameter contributing to the diffusion; softening enters here, raising the lower cutoff.

\ln\Lambda = \ln\!\left(\frac{b_{\max}}{b_{\min}}\right),\qquad b_{\min}\simeq\max\!\left(\frac{Gm}{v^{2}},\,\epsilon\right) \qquad (4)

where: $b_{\max}\sim R$ — system size; $b_{\min}$ — the $90^\circ$ -deflection radius $Gm/v^{2}$ , or $\epsilon$ if larger; $v$ — relative speed. Bigger $\epsilon$ : bigger $b_{\min}$ , smaller $\ln\Lambda$ , longer $t_{\mathrm{relax}}$ . The optimal softening minimises the total force error and shrinks as particles are added, so resolution buys less bias and less noise at once.

\epsilon_{\mathrm{opt}} \approx a\,N^{-1/3}\ \ \text{(interparticle spacing)}\ \ \longrightarrow\ \ \propto N^{-0.28} \qquad (5)

where: $\epsilon_{\mathrm{opt}}$ — softening minimising total force error; $a$ — system scale; $N$ — particle number. The naive $N^{-1/3}$ spacing is bracketed by the Merritt–Dehnen $N^{-0.28}$ ; both fall slowly. The Power criterion fixes where the profile is trustworthy: the density is reliable only outside the radius where the enclosed relaxation time drops to the orbital time.

\frac{t_{\mathrm{relax}}(<r)}{t_{\mathrm{circ}}(r)}\gtrsim 1\ \Rightarrow\ r>r_{\mathrm{conv}},\qquad t_{\mathrm{relax}}(<r)=\frac{N(<r)}{8\ln N(<r)}\,t_{\mathrm{circ}}(r) \qquad (6)

where: $N(<r)$ — enclosed particle number; $t_{\mathrm{circ}}(r)$ — circular orbit period; $r_{\mathrm{conv}}$ — convergence radius. Inside $r_{\mathrm{conv}}$ so few particles are enclosed that they relax within a Hubble time — the profile is an artifact. Finally, the quantity every code evaluates each step is the softened acceleration — a pair sum with $\epsilon$ smoothing every denominator.

\mathbf{a}_{i} = -\sum_{j\neq i}\frac{G m_{j}\,(\mathbf{x}_{i}-\mathbf{x}_{j})}{\left(|\mathbf{x}_{i}-\mathbf{x}_{j}|^{2}+\epsilon^{2}\right)^{3/2}} \qquad (7)

where: $\mathbf{a}_{i}$ — acceleration of particle $i$ ; $m_{j}$ — mass of particle $j$ ; $\mathbf{x}_{i},\mathbf{x}_{j}$ — positions; $\epsilon$ — softening. The $\epsilon^{2}$ under the power is the whole prescription; set it to zero and exact gravity — with its divergences — returns.

Derivation ARelaxation time from the diffusion of orbital energy

The goal is to build the scaling $t_{\mathrm{relax}}\propto N/\ln N$ from first principles, by following one test star across a system of $N$ field particles, summing the random velocity kicks it receives, and demanding their accumulated size reach the star's own speed. This exposes why relaxation is a random walk (variances add, not amplitudes) and why it is fatally sensitive to $N$ through the macro-particle mass $m=M/N$ .

Step 1. A distant encounter at impact parameter $b$ and relative speed $v$ produces a perpendicular deflection $\delta v_\perp\approx 2Gm/(bv)$ , the impulse result from aq-astro-gravitational-two-body-problem. Why: for weak encounters the star barely deviates, so integrating the transverse pull along a straight line gives a clean $1/b$ law.

Step 2. In one crossing the star meets $\mathrm{d}n = n_\star\,v\,t_{\mathrm{cross}}\,2\pi b\,\mathrm{d}b$ field particles with $b$ in $[b,b+\mathrm{d}b]$ , where $n_\star=N/(\tfrac{4}{3}\pi R^{3})$ . Why: it sweeps an annulus $2\pi b\,\mathrm{d}b$ through a column of length $\approx R$ , so the count is density times swept volume.

Step 3. The kicks point randomly, so their vectors cancel but their squares add: $\Delta v_\perp^{2}=\int(\delta v_\perp)^{2}\,\mathrm{d}n = \dfrac{8N(Gm)^{2}}{R^{2}v^{2}}\ln\Lambda$ , the $\int\mathrm{d}b/b$ giving $\ln\Lambda$ . Why: relaxation is diffusion — a random walk in velocity — so the mean square accumulates and the log records that every decade of $b$ contributes equally.

Step 4. With the virial estimate $v^{2}\sim GNm/R$ , the fractional change per crossing collapses to $\dfrac{\Delta v_\perp^{2}}{v^{2}}\approx\dfrac{8\ln\Lambda}{N}$ . Why: the dependence on $G,m,R$ folds into the virial relation, leaving a pure $1/N$ — the sampling granularity — times the log.

Step 5. Order-unity relaxation then needs $n_{\mathrm{relax}}=v^{2}/\Delta v_\perp^{2}\approx N/(8\ln\Lambda)$ crossings, so $t_{\mathrm{relax}}=n_{\mathrm{relax}}\,t_{\mathrm{cross}}\simeq(0.1\,N/\ln N)\,t_{\mathrm{cross}}$ . Since $m=M/N$ , this $1/N$ means a coarse run heats far faster than the object it models, vanishing only as $N\to\infty$ . Why: the relaxation time is the crossing time stretched by $N/\ln N$ ; the blame is the macro-particle mass, so the cure is more particles, not more softening. *Verified: reproduces Binney & Tremaine eq. (1.38) and the Carroll & Ostlie cluster relaxation formula; the $\ln\Lambda$ and $N$ dependence are standard, and $t_{\mathrm{relax}}\simeq(0.1N/\ln N)\,t_{\mathrm{cross}}$ was checked numerically for a globular cluster below.*

Go deeper · Derivation B · 1 minHow Plummer softening caps the force and builds a harmonic core

Derivation BHow Plummer softening caps the force and builds a harmonic core

The goal is to show that replacing the point-mass potential by the Plummer kernel $\Phi=-Gm/\sqrt{r^{2}+\epsilon^{2}}$ converts a divergent $1/r^{2}$ force into a bounded one with a definite maximum, and makes the deep interior a linear harmonic restoring force. This is what makes a softened integration stable — a finite peak force means a finite shortest timestep — and where the sub- $\epsilon$ bias is born.

Step 1. Differentiate: $F(r)=-\partial\Phi/\partial r=-Gm\,r/(r^{2}+\epsilon^{2})^{3/2}$ . Why: the numerator $r$ from the chain rule forces $F\to0$ at the origin instead of blowing up — the divergence is gone by construction.

Step 2. Set $\mathrm{d}F/\mathrm{d}r=0$ : with $F\propto r(r^{2}+\epsilon^{2})^{-3/2}$ the derivative numerator $\propto(\epsilon^{2}-2r^{2})$ vanishes at $r=\epsilon/\sqrt{2}$ . Why: this is the one interior scale softening introduces — the harmonic-to-Keplerian crossover, fixed at a fraction of $\epsilon$ regardless of mass or speed.

Step 3. Evaluate there: $F_{\max}=Gm(\epsilon/\sqrt{2})/(\tfrac{3}{2}\epsilon^{2})^{3/2}=2Gm/(3\sqrt{3}\,\epsilon^{2})\approx0.385\,Gm/\epsilon^{2}$ . Why: the largest kick scales as $\epsilon^{-2}$ , so $\epsilon$ directly sets the maximum acceleration and the minimum stable timestep $\sim\sqrt{\epsilon^{3}/Gm}$ .

Step 4. For $r\ll\epsilon$ : $(r^{2}+\epsilon^{2})^{3/2}\to\epsilon^{3}$ , so $F\to-(Gm/\epsilon^{3})\,r$ — a spring, a uniform harmonic ball that is stable but physically wrong, the source of the sub- $\epsilon$ bias. Conversely for $r\gg\epsilon$ , $F\to-Gm/r^{2}$ to error $\sim(\epsilon/r)^{2}$ . Why: softening is short-range surgery only — the interior is capped and harmless while large-scale dynamics stays exactly Newtonian. *Verified: the peak location $r=\epsilon/\sqrt{2}$ and value $F_{\max}=2Gm/(3\sqrt{3}\,\epsilon^{2})=0.3849\,Gm/\epsilon^{2}$ were confirmed numerically; the limits match Binney & Tremaine App. C and Dehnen 2001.*

Go deeper · Derivation C · 1 minOptimal softening as a bias–variance minimisation

Derivation COptimal softening as a bias–variance minimisation

The goal is to turn the qualitative "not too big, not too small" advice about $\epsilon$ into a quantitative optimum, by writing the mean-square force error as a smoothing bias that grows with $\epsilon$ plus a Poisson sampling variance that falls with $\epsilon$ , then minimising the total. The exercise shows the optimal softening shrinks with particle number — the statistician's bias–variance tradeoff exactly.

Step 1. The bias comes from smoothing the field over $\epsilon$ , shifting the force at second order: $\Delta\mathbf{F}_{\mathrm{bias}}\propto\epsilon^{2}\nabla^{2}\mathbf{F}$ , so the squared bias $\propto\epsilon^{4}$ . Why: a symmetric kernel has zero first moment, so the leading error is the second moment times the field's curvature — a wider blur distorts more.

Step 2. The variance comes from sampling with $N$ discrete particles; averaging over the $\sim N\epsilon^{3}$ neighbours inside a softening volume gives $\Delta\mathbf{F}_{\mathrm{var}}^{2}\propto 1/(N\epsilon)$ . Why: shot noise falls as one over the number of contributors, and larger $\epsilon$ averages over more — so noise and bias pull $\epsilon$ opposite ways.

Step 3. Add them, $\langle|\Delta\mathbf{F}|^{2}\rangle=C_b\,\epsilon^{4}+C_v/(N\epsilon)$ , and set the derivative to zero: $4C_b\,\epsilon^{3}=C_v/(N\epsilon^{2})$ . Why: the total error is convex — variance-dominated at small $\epsilon$ , bias-dominated at large $\epsilon$ — so a unique minimum exists where the slopes balance.

Step 4. Solve: $\epsilon_{\mathrm{opt}}^{5}=C_v/(4C_b N)$ , hence $\epsilon_{\mathrm{opt}}\propto N^{-1/5}$ for this toy model; realistic density-dependent analyses (Merritt 1996; Dehnen 2001) sharpen the exponent to $-1/3$ to $-0.28$ , straddling the mean interparticle spacing. Why: even the crude estimate captures the message — the optimum falls with $N$ — while folding in the varying density curvature steepens the fall onto the spacing scaling codes adopt. *Verified: reproduces the bias–variance structure and the reported exponent range of Merritt (1996) and Dehnen (2001); the scaling $\epsilon_{\mathrm{opt}}\propto N^{-1/3}$ was used in the worked example below.*

Limiting cases

Case	Result	Meaning
Large separation, $r\gg\epsilon$	$\Phi\to-\dfrac{Gm}{r},\ F\to-\dfrac{Gm}{r^{2}}$	Softening negligible; exact Newtonian gravity recovered, large-scale dynamics unbiased.
Small separation, $r\ll\epsilon$	$F\to-\dfrac{Gm}{\epsilon^{3}}\,r$	Harmonic restoring force; the singular $1/r^{2}$ is capped, but sub- $\epsilon$ structure is unresolved.
Collisionless limit, $N\to\infty$	$t_{\mathrm{relax}}/t_{\mathrm{cross}}\sim\dfrac{0.1\,N}{\ln N}\to\infty$	Relaxation diverges vs the crossing time; the CBE solution is recovered and heating vanishes.
Small $N$ (poor resolution)	$t_{\mathrm{relax}}\lesssim t_{\mathrm{age}}$	Heating erases cusps within the age; inner profiles collapse to artificial cores.
Optimal softening	$\epsilon_{\mathrm{opt}}\sim N^{-1/3}\ \text{to}\ N^{-0.28}\ (\times\ \text{scale})$	Minimises total force error by balancing bias against shot-noise variance.
Cold pair ( $v^{2}<Gm/\epsilon$ )	$b_{\min}\to\epsilon$ in $\ln\Lambda$	Softening, not the deflection radius, sets the inner cutoff of the encounter integral.

Worked examples

1 · Relaxation time of a globular cluster (fully worked)

A globular cluster has $N=10^{5}$ stars, radius $R\approx10\ \mathrm{pc}$ , and velocity dispersion $v\approx10\ \mathrm{km/s}$ . Is two-body relaxation a real, physical process here within a Hubble time?

Step 1. Crossing time. $t_{\mathrm{cross}}=R/v=(10\times3.086\times10^{13}\ \mathrm{km})/(10\ \mathrm{km/s})=3.09\times10^{13}\ \mathrm{s}\approx\mathbf{1.0\times10^{6}\ yr}$ — stars orbit through in about a million years.

Step 2. Coulomb logarithm. With $\Lambda\sim N$ , $\ln\Lambda\approx\ln(10^{5})=\mathbf{11.5}$ (the refined $\Lambda\approx0.4N$ gives $10.6$ , immaterial here).

Step 3. Relaxation time. $t_{\mathrm{relax}}=\dfrac{N}{8\ln\Lambda}\,t_{\mathrm{cross}}=\dfrac{10^{5}}{8\times11.5}\times1.0\times10^{6}\ \mathrm{yr}\approx\mathbf{1.1\times10^{9}\ yr}\approx1\ \mathrm{Gyr}$ .

Step 4. Interpretation. The cluster age is $\sim1.2\times10^{10}\ \mathrm{yr}$ , so it has relaxed $\sim$ ten times over. Two-body relaxation is genuinely physical here — it drives mass segregation, core collapse, and evaporation, which is why direct integrators use tiny or zero softening: softening these encounters away would delete the physics. This is the opposite of a collisionless galaxy simulation.

2 · Numerical heating and the $N$ dependence (fully worked)

A subhalo is resolved with $N=1000$ bound particles and has an internal crossing time $t_{\mathrm{cross}}\approx0.5\ \mathrm{Gyr}$ . Will two-body heating corrupt it within a Hubble time? How much better is $N=10^{5}$ ?

Step 1. Use $t_{\mathrm{relax}}\simeq(0.1\,N/\ln N)\,t_{\mathrm{cross}}$ . For $N=1000$ : $\ln N=6.9$ , so $t_{\mathrm{relax}}=0.1\times1000/6.9\times0.5\ \mathrm{Gyr}=\mathbf{7.2\ Gyr}$ .

Step 2. Compare to the age. $7.2\ \mathrm{Gyr}$ is less than the $\sim13.8\ \mathrm{Gyr}$ age, so this subhalo relaxes within the run: its inner density is artificially heated and it may spuriously evaporate — a thousand particles cannot represent a bound clump collisionlessly.

Step 3. Now $N=10^{5}$ at the same $t_{\mathrm{cross}}$ : $\ln N=11.5$ , so $t_{\mathrm{relax}}=0.1\times10^{5}/11.5\times0.5\approx\mathbf{435\ Gyr}$ — $\sim30$ Hubble times, safely collisionless.

Step 4. The lesson. Raising $N$ by $100\times$ lengthened $t_{\mathrm{relax}}$ by $\sim60\times$ (the $N/\ln N$ scaling), turning an unreliable clump into a trustworthy one. This is why subhalo studies demand thousands of bound particles per object — the smallest resolved haloes are always the most suspect.

3 · Choosing the softening length and its force cap

A Milky-Way-mass halo, scale radius $a\approx20\ \mathrm{kpc}$ , is sampled with $N=10^{7}$ particles. Estimate the optimal softening and its force cap.

Step 1. Interparticle-spacing softening: $\epsilon_{\mathrm{opt}}\approx a\,N^{-1/3}=20\ \mathrm{kpc}/215\approx\mathbf{93\ pc}$ . Published runs often adopt a few percent of this — a few pc to tens of pc — choosing $\epsilon$ below the mean spacing to preserve resolution.

Step 2. The force cap is $F_{\max}=2Gm/(3\sqrt{3}\,\epsilon^{2})\approx0.385\,Gm/\epsilon^{2}$ with $m=M/N$ . Smaller $\epsilon$ raises this cap and shortens the minimum timestep $\sim\sqrt{\epsilon^{3}/Gm}$ — the cost that stops one softening arbitrarily little. The optimum: $\epsilon$ small enough to resolve the scales of interest, large enough to keep the timestep manageable.

Physics of the portrait

Two-body relaxation is the concept that ties this article together: a diffusion clock, running at rate $\propto\ln\Lambda/N$ , that measures how quickly a granular system forgets its collisionless past. It was worked out for real star clusters decades before it returned, uninvited, as the dominant systematic error of the digital universe.

Origins · 1 mintwo-body relaxation and how this was discovered

two-body relaxation

the discreteness clock t_relax ∝ N/ln N · collisionless stellar dynamics and N-body simulation

Two-body relaxation began as pure stellar dynamics. In 1943 Subrahmanyan Chandrasekhar, studying dynamical friction, derived the diffusion of a star's velocity from the cumulative pull of its neighbours and found the timescale set by the Coulomb logarithm $\ln\Lambda$ — the same logarithm Landau had introduced for plasmas. For a real galaxy the time is absurd: with $N\sim10^{11}$ stars, $t_{\mathrm{relax}}$ exceeds the age of the Universe by a factor of millions, which is exactly why a galaxy is collisionless. The concept slept until the computer arrived. In 1960 Sebastian von Hoerner ran the first digital $N$ -body integrations, hit the singular close encounter, and invented force softening to keep the run alive; Sverre Aarseth then built the direct $N$ -body codes for star clusters, where relaxation is the physics of interest rather than an error. As tree and TreePM codes pushed collisionless cosmology to billions of particles, Chandrasekhar's clock returned as a villain: Power and collaborators (2003) showed that inside a well-defined radius the finite- $N$ relaxation time falls below the halo's age, heating the cusp into a spurious core. The same equation that explains why galaxies do not relax now tells us where our simulations cannot be trusted.

History

1943

— Chandrasekhar derives velocity diffusion and dynamical friction, fixing t_relax through the Coulomb logarithm

1960

— von Hoerner runs the first digital N-body integrations and introduces gravitational softening

1986

— Barnes and Hut publish the tree algorithm, making large-N collisionless simulations tractable

2003

— Power and collaborators establish the convergence radius for cosmological halo simulations

2005

— Springel's GADGET-2 evolves the Millennium simulation with ten billion particles

Going further

Texts. Binney & Tremaine, Galactic Dynamics, 2nd ed., §1.2 and §4.1 (relaxation), App. C (softening and $N$ -body methods) · Carroll & Ostlie, An Introduction to Modern Astrophysics, Ch. 25 (star clusters) · Dodelson, Modern Cosmology (numerical structure formation) · Mo, van den Bosch & White, Galaxy Formation and Evolution, Ch. 5.

Papers. Power et al. 2003, MNRAS 338, 14 (convergence radius) · Dehnen 2001, MNRAS 324, 273 (optimal softening) · Merritt 1996, AJ 111, 2462 (bias–variance) · Springel et al. 2005, Nature 435, 629 (Millennium).

Provenance. The relaxation scaling $t_{\mathrm{relax}}\simeq(0.1\,N/\ln N)\,t_{\mathrm{cross}}$ follows Binney & Tremaine §1.2 (eq. 1.38) and Carroll & Ostlie Ch. 25, checked numerically for the globular example ( $\approx1.1\ \mathrm{Gyr}$ ). The Plummer potential/force, the peak $F_{\max}=2Gm/(3\sqrt{3}\,\epsilon^{2})$ at $r=\epsilon/\sqrt{2}$ ( $=0.3849\,Gm/\epsilon^{2}$ , confirmed), and the harmonic small- $r$ limit are standard (BT App. C; Dehnen 2001). The softened $\ln\Lambda$ follows BT §1.2 and Merritt 1996; the $N^{-1/3}$ – $N^{-0.28}$ optimal-softening scaling follows Merritt (1996) and Dehnen (2001); the convergence criterion follows Power et al. (2003). A verifier should recompute $F_{\max}$ and $r_{\mathrm{peak}}$ , confirm the $\int\mathrm{d}b/b$ divergence, and check the $N$ -scaling of $t_{\mathrm{relax}}$ .

Prerequisites. Gravitational N-Body Problem · Two-Body Relaxation Time · Collisionless Boltzmann Equation · Gravitational Two-Body Problem · Distribution Function and Phase Space. Leads to. Virial Theorem in Stellar Systems · Jeans Equations · Galactic Potential Models · Orbits in Spherical Potentials.

Ask the tutor: ask me to work the Power radius for a specific halo, or to show how the same $t_{\mathrm{relax}}\propto N/\ln N$ is a friend in a star cluster and an enemy in a cosmological box.

Isn't the mean-field potential all that matters? Why should discreteness change anything?

The mean field carries the bulk of the dynamics — orbits, the potential well, the large-scale collapse — and even a modest $N$ gets it nearly right, because it is a smooth integral over the density that averages sampling noise away. What discreteness adds is a fluctuating residual: because each macro-particle is a lump, not a smear, the force jitters as neighbours pass. Those jitters never average to zero in the mean-square sense — they diffuse orbital energies, and over a relaxation time that rearranges the distribution the mean field alone would preserve forever.

If softening biases the forces, why not just use zero softening and a tiny timestep?

For a genuinely collisional system (a small star cluster) that is exactly what you do — direct codes use near-zero softening and regularise the singular binaries analytically. But for a collisionless system it backfires: the rare close pair encounters between macro-particles are precisely the artificial events you want to avoid, and they are hugely amplified because $m=M/N$ is enormous. You would spend vast effort resolving encounters the real object never has, while the singular forces demand near-zero timesteps. Softening deliberately blurs those passages, buying stability at the price of a known bias below $\epsilon$ . The right move is optimal softening tuned to $N$ , not zero.

How is numerical two-body heating different from ordinary error like energy drift?

Energy drift from a finite timestep is a property of the integrator: it shrinks as you shorten the step, vanishes for a perfect integrator, and does not depend on $N$ . Two-body heating is a property of the model, not the solver. Even with a flawless energy-conserving integrator and infinitesimal steps, a finite- $N$ system still relaxes, because the granularity is baked into the initial conditions. You cannot integrate it away; the only cure is more particles. That is why convergence tests vary $N$ , not just the timestep.

What actually sets the Power radius, and how do I find it in practice?

It is the radius inside which the enclosed relaxation time has fallen below the age of the system. Because $N(<r)$ drops inward, the local relaxation time drops too, and at some radius it crosses the age — inside, the profile has been heated. In practice you evaluate, at each radius, the empirical Power quantity built from $N(<r)$ and the enclosed mean density relative to critical, and find where it crosses the calibrated threshold (about $0.6$ in the 2003 form). Everything outside is converged; inside is flagged. Raising $N$ moves the crossing inward.

Does softening make a simulation collisionless, or just less collisional?

Just less collisional — softening reduces heating but never abolishes it. Raising $\epsilon$ increases $b_{\min}$ , which lowers $\ln\Lambda$ and lengthens $t_{\mathrm{relax}}$ ; but $\ln\Lambda$ enters only logarithmically, so even generous softening buys a factor of a few, not orders of magnitude. The dominant lever is $N$ , which enters linearly. A run is effectively collisionless only when $N$ makes $t_{\mathrm{relax}}\gg t_{\mathrm{age}}$ everywhere you care about. Never treat large softening as a substitute for particle number.

Why does the Coulomb logarithm appear, and why is it "only" logarithmic?

It comes from summing encounters over all impact parameters. Each decade of $b$ contributes about equally to the mean-square velocity change, because the kick scales as $1/b$ but the number of encounters scales as $b$ , so $(\delta v)^{2}\,\mathrm{d}n\propto \mathrm{d}b/b$ — whose integral is a logarithm. The upper cutoff $b_{\max}$ is the system size; the lower $b_{\min}$ is the $90^\circ$ -deflection radius $Gm/v^{2}$ or the softening $\epsilon$ , whichever is larger. Because weak distant encounters dominate, the sum grows only as the log of the range — which is why $\ln\Lambda\sim10$ across a huge range of systems.

Gravitational Softening and Two-Body Heating

Where this lives today

The physical picture

The core equations

Limiting cases

Worked examples

1 · Relaxation time of a globular cluster (fully worked)

2 · Numerical heating and the $N$ dependence (fully worked)

3 · Choosing the softening length and its force cap

Physics of the portrait

Going further

Before this

Where this leads

Gravitational Softening and Two-Body Heating

Where this lives today

The physical picture

The core equations

Limiting cases

Worked examples

1 · Relaxation time of a globular cluster (fully worked)

2 · Numerical heating and the NNN dependence (fully worked)

3 · Choosing the softening length and its force cap

Physics of the portrait

Going further

Before this

Where this leads

2 · Numerical heating and the $N$ dependence (fully worked)