  We present a method for computing the 6-dimensional coarse-grained
phase-space density $f(x,v)$ in an N-body system, and derive its distribution
function $v(f)$. The method is based on Delaunay tessellation, where $v(f)$ is
obtained with an effective fixed smoothing window over a wide $f$ range. The
errors are estimated, and $v(f)$ is found to be insensitive to the sampling
resolution or the simulation technique. We find that in gravitationally relaxed
haloes built by hierarchical clustering, $v(f)$ is well approximated by a
robust power law, $v(f) \propto f^{-2.5 \pm 0.05}$, over more than 4 decades in
$f$, from its virial level to the numerical resolution limit. This is tested to
be valid in the $\Lambda$CDM cosmology for haloes with masses
$10^9-10^{15}\msun$, indicating insensitivity to the slope of the initial
fluctuation power spectrum. By mapping the phase-space density in position
space, we find that the high-$f$ end of $v(f)$ is dominated by the ``cold''
subhaloes rather than the parent-halo central region and its global spherical
profile. The value of $f$ in subhaloes near the virial radius is typically
$>100$ times higher than its value at the halo centre, and it decreases
gradually from outside in toward its value at the halo centre. This seems to
reflect phase mixing due to mergers and tidal effects involving puffing up and
heating. The phase-space density can thus provide a sensitive tool for studying
the evolution of subhaloes during the hierarchical buildup of haloes. It
remains to be understood why the evolved substructure adds up to the actual
universal power law of $v(f) \propto f^{-2.5}$. It seems that this behaviour
results from the hierarchical clustering process and is not a general result of
violent relaxation.
