  We present a novel perspective on the clustering of dark matter in phase
space by defining the particle phase space average density ($P^2SAD$) as a
two-dimensional extension of the two-point correlation function averaged within
a certain volume in phase space. This statistics is a sensitive measure of
small scale (sub-)structure of dark matter haloes. By analysing the structure
of $P^2SAD$ in Milky-Way-size haloes using the Aquarius simulations, we find it
to be nearly universal at small scales, i.e. small separations in phase space,
where substructures dominate. This remarkable universality occurs across time
and in regions of substantially different ambient densities (by nearly four
orders of magnitude), with typical variations in $P^2SAD$ of a factor of a few.
The maximum variations occur in regions where substructures have been strongly
disrupted. The universality is also preserved across haloes of similar mass but
diverse mass accretion histories and subhalo distributions. The universality is
also broken at large scales, where the smooth dark matter distribution in the
halo dominates. Although at small scales the structure of $P^2SAD$ is roughly
described by a subhalo model, we argue that the simulation data is better
fitted by a family of superellipse contours. This functional shape is inspired
by a model that extends the stable clustering hypothesis into phase space. In a
companion paper, we refine this model and show its advantages as a method to
obtain predictions for non-gravitational signatures of dark matter.
