Consider an isolated system of $N$ identical nonrelativistic classical particles that move within a $d$-dimensional spatial manifold according to a given Hamiltonian function.

Liouville's theorem states that "the phase-space distribution function is constant along the trajectories of the system." But I'm not entirely clear on how "the phase-space distribution function" should be interpreted, and whether "the phase space" under consideration is the full "microscopic" $2Nd$-dimensional phase space or the "macroscopic" $2d$-dimensional single-particle phase space.

There seem to be two conceptually quite different formulations of the theorem:

In version 1, the phase space under consideration is the $2d$-dimensional single-particle phase space. In this case, the phase space distribution function $\rho({\bf q}, {\bf p})$ has a very simple interpretation. It's simply a number density: $\rho({\bf q}, {\bf p})\ d^d{\bf q}\ d^d{\bf p}$ equals the number of particles in the infinitesimal real-space volume $d^d {\bf q}$ at position ${\bf q}$ whose canonical momenta lie within the infinitesimal volume $d^d {\bf p}$ in momentum space at momentum ${\bf p}$. So $\int d^d{\bf q}\ d^d{\bf p}\ \rho({\bf q},{\bf p}) = N$. We can also trivially divide $\rho({\bf q}, {\bf p})$ by $N$ to get a normalized probability density function that gives the probability (density) that a particle selected uniformly at random has position ${\bf q}$ and momentum ${\bf p}$.

In version 2, the phase space under consideration is the full $2Nd$-dimensional phase space. In this case, the interpretation of the phase space distribution function $\rho({\bf q}, {\bf p})$ is more abstract: it corresponds to an abstract statistical uncertainty over microstates that does not correspond to any kind of number density or count of microstates (since the individual points in the full $2Nd$-dimensional phase space themselves are the microstates).

To pick a few random examples: this article seems to be taking the first perspective (since it talks about the "density of system points" - although the article is inconsistent in its use of the variable $n$, so it isn't entirely clear). This article, this one, and this one explicitly say that the distribution function is a number density (and one of them clearly implies that it is defined on the single-particle phase space). This article identifies the distribution function with a count of microstates, so it seems to be taking the first perspective as well.

On the other hand, this article seems to be taking the second perspective, because it says that the phase space is "$6N$-dimensional" (although it doesn't define $N$).

Each formulation seems to have its pros and cons. The first formulation is concrete and conceptually straightforward to me, since we're just talking about the density of particles evolving over time. But if "the Hamiltonian" is just defined on the single-particle phase space, then it isn't clear to me how it can incorporate interactions between the particles. So it seems to only apply to noninteracting ideal gases, which is a pretty serious limitation.

The second formulation can handle arbitrary interacting many-body Hamiltonians, but I have a hard time conceptualizing exactly what it means. It deals with the time-evolution of an abstract "phase-space distribution function" that (I think) cannot represent any kind of number density over microstates (since we're already operating at the "microstate" level). So what exactly is the theorem saying?

My problem is that some explanations seem to conflate these two pictures, by first describing the phase-space distribution function as representing some kind of number density function over macrostates, but then when actually doing the math, suddenly switching to a $2Nd$-dimensional "microscopic" phase space that seems to already be working at the microstate level. At this microstate level, shouldn't the state of the system just be single deterministic point in the "microscopic" phase space?

At some level, I guess the answer is just "The same mathematical derivation works in both cases, so either interpretation is valid." But this answer isn't fully satisfying to me. Is either version the "better" or "more useful" interpretation of Liouville's theorem? The first version seems extremely limited, since it can't handle interparticle interactions, but the second version seems conceptually confusing and less useful in practice.