The correct phase space $X$ of a classical system describing an arbitrary number of identical particles is
$$
X = \coprod_{n \in \mathbb{N}} \mathbb{R}^{6n},
$$
where $\coprod$ denotes the disjoint union. An element of $X$ has the form $(n,x)$ where $n\in \mathbb{N}$ is the number of particles and $x\in \mathbb{R}^{6n}$ is the vector of momenta/positions of the $n$ particles.
Now there is not a single probability density function. Instead there is a function $\rho_n : \mathbb{R}^{6n} \to [0,\infty)$ for every $n \in \mathbb{N}$
such that the probability $P(n,\Gamma)$ to observe the system consisting of $n$ particles with their momenta and positions in the phase space volume $\Gamma \subset \mathbb{R}^{6n}$ is
$$
P(n,\Gamma) = \int_\Gamma \rho_n (x) dx^{6n}.
$$
Note that the normalisation condition here is, that
$$
1 = \sum_{n=0}^\infty \int_{\mathbb{R}^{6n}} \rho_n (x) dx^{6n},
$$
i.e. the probability of observing any number of particles with any position and momenta is 1.
I am not sure about the "ergodic hypothesis" in this context, maybe you can elaborate what exactly the "ergodic hypothesis" means to you.
Elaborating on the ideas discussed in the comments (in particular HTNW's comment):
In the case of multiple kind (species) of particles the phase space can be constructed as follows:
Let $S$ be the (finite) set of species. Define the state space $X$ of a classical system describing an arbitrary number of particles from different species (from $S$) as
$$
X = \prod_{s \in S} \coprod_{n \in \mathbb{N}} \mathbb{R}^{6n}.
$$
So an element of $X$ (a state of the system) is a tuple of the form $(n_s,x_s)_{s \in S}$, where $n_s$ is the number of particles of species $s$ in the state and $x_s$ is the vector storing positions and momenta of these $n_s$ particles.
A probability measure can be constructed on this space by considering the product of the probability measures defined above for each species.
