  Self-gravitating astronomical objects often show a central plateau in the
density profile (core) whose physical origin is hotly debated. Cores are
theoretically expected in N-body systems of maximum entropy, however, they are
not present in the canonical N-body numerical simulations of cold dark matter
(CDM). Our work shows that despite this apparent contradiction between theory
and numerical simulations, they are fully consistent. Simply put, cores are
characteristic of systems in thermodynamic equilibrium, but thermalizing
collisions are purposely suppressed in CDM simulations. When collisions are
allowed, N-body numerical simulations develop cored density profiles, in
perfect agreement with the theoretical expectation. We compare theory and two
types of numerical simulations: (1) when DM particles are self-interacting
(SIDM) with enough cross-section, then the effective two-body relaxation
timescale becomes shorter than the Hubble time resulting in cored DM haloes.
The haloes thus obtained, with masses from dwarf galaxies to galaxy clusters,
collapse to a single shape after normalization, and this shape agrees with the
polytropic density profile theoretically expected. (2) The inner radii in
canonical N-body numerical simulations are always discarded because the use of
finite-mass DM particles artificially increases the two-body collision rate. We
show that the discarded radii develop cores that are larger than the employed
numerical softening and have polytropic shapes independently of halo mass. Our
work suggests that the presence of cores in simulated (or observed) density
profiles can used as evidence for systems in thermodynamic equilibrium.
