Phase space volume and relativity Much of statistical mechanics is derived from Liouville's theorem, which can be stated as "the phase space volume occupied by an ensemble of isolated systems is conserved over time." (I'm mostly interested in classical systems for the moment.) 
It's clear that special relativity doesn't change this, since relativity just adds a different set of invariances to the Hamiltonian. So under special relativity, the phase space volume of an ensemble must remain constant over time, as long as we're consistent in using a given inertial reference frame.  However, this doesn't tell us how the phase space volume would change if we boosted to a different reference frame. So I'd like to know the following: 


*

*Can anything useful be said about how the phase space volume occupied by an ensemble of systems changes under a Lorentz boost? I guess this would entail taking a different time-slice through the system, as explained below. I would suspect that phase space volume doesn't behave very nicely under Lorentz transforms, since phase space is defined in terms of 3-position and 3-momentum.

*If phase space volume doesn't behave nicely under Lorentz boosts, is there a generalisation (e.g. in terms of 4-position and 4-momentum) that does? And if so, where can I go to read more about it?

*Optional bonus questions: what effect do general relativity and quantum mechanics have on all of the above?
Additional clarification 
Below are some diagrams which should make it clearer what I'm asking for. (a) shows a space-time diagram of a system consisting of several classical particles. We take a simultaneous time-slice though the system (red line) and note the 3-position and 3-momentum of each particle. This gives us a point in a phase space (diagram (b)), which I've drawn as 2 dimensional for convenience. An ensemble of similar systems (meaning systems with the same dynamics but different initial conditions) can be thought of as occupying a region of this phase space. Liouville's theorem tells us that if we do the same thing at different values of $t$ (for the same ensemble of systems), the shape of this region may change but its volume will be the same.
(c) shows a Lorentz-transformed version of the same system as (a). We can take a simultaneous time-slice through the system in the new $(x',t')$ reference frame, but it will not be parallel to the one in the $(x, t)$ frame. The particles' 3-positions and 3-momenta will also be different. We can plot the system's position in the new phase space generated by doing the same procedure in the new reference frame (d). We can also plot the region of phase space occupied by the ensemble. Doing the same thing at different values of $t'$ will produce different regions with the same volume as each other. However, my question is about whether volume of the region plotted in (d) must equal the volume of the region in (b).

Progress
In this document (J. Goodman, Topics in High-Energy Astrophysics, 2012, p.12-13; unfortunately Goodman gives no further reference) there is a proof that infinitesimal phase space volumes are Lorentz invariant. It looks legit, but the author assumes that every particle's $x$ position is within the same small interval $[x, x+dx]$, which means that he doesn't have to take account of the fact that you take a different time slice when you change the reference frame. Additionally, I'm using an integrated version of Liouville's theorem, in which the ensemble has a finite rather than infinitesimal phase space volume, and it isn't immediately clear to me whether this makes a difference. So this seems to suggest that the phase space volumes sketched in figures (b) and (d) above will be equal, but I'm still not sure and would like to know where I can find the full proof, if it exists.
Further Progress
I'm awarding a 100 point bounty to Qmechanic for his helpful answer, which expresses Goodman's argument (discussed above) in more formal language. However, I don't think the question has been answered yet. This is because the argument doesn't just assume that the system occupies an infinitesimal volume of phase space. On its own this would be fine, since you can just integrate over phase space to get the finite phase-space volume version. The problem is that this argument also assumes the system occupies only an infinitesimal volume of real space. In general, a system that occupies an infinitesimal volume of phase space can occupy a finite area of actual space (think of a system composed of multiple, spatially separated, classical particles, as illustrated above) and so this argument seems to cover only a restricted subset of the cases I'm interested in. Goodman's argument is fine if you're considering non-interacting particles or fluids in equilibrium (which makes sense given its origins in astrophysics) but I'm also interested in multi-particle systems that may or may not be in thermal equilibrium.  It is my strong intuition that the argument can be extended to deal with all cases, and I'll happily award an additional 100 point bounty to anyone who can show how to do so.
I think one route to an answer might be to note that, in special relativity, classical particles cannot interact unless they actually collide. This probably means that in the times in between collisions they can be thought of as occupying their own individual phase spaces after all, and perhaps with the aid of the proof that's been presented and some careful accounting about what happens during collisions, it will all work out OK. But of course then it has to be extended to deal with (classical) fields, which could be tricky. One gets the feeling there ought to be a nice simple way to derive it from Lorentz invariance of the Lagrangian, without worrying about what type of system we're dealing with. (I'll award a 200 point bounty to a satisfactory answer that takes this latter approach!)
 A: EDIT: Upon thinking about it a bit, i am no longer sure this proof is actually correct. Nonetheless I'll leave this answer up for bookkeeping until I figure it out.

Here is a rather simple proof. We first make some assumptions:

*

*The momentum $\mathbf p$ and the Hamiltonian $H(\mathbf x,\mathbf p,t)$ behave as a four-vector together under a Lorentz transformation, i.e. as $$ p^\mu=(H,\mathbf p). $$

*Since invariance under rotations and translations is obvious, it is sufficient to consider only boosts in one spatial direction. Hence I will work with a two dimensional spacetime with coordinates $(t,x)$ ($c=1$)

A Lorentz transformation then has the form $$ \left(\begin{matrix} t^\prime \\ x^\prime\end{matrix}\right)=\left(\begin{matrix}\cosh(\phi) & -\sinh(\phi) \\ -\sinh(\phi) & \cosh(\phi)\end{matrix}\right)\left(\begin{matrix} t \\ x\end{matrix}\right), $$ and the same is true with $t$ and $x$ replaced by $H$ and $p$ respectively. In particular, since the Hamiltonian is a function, this means $$ H^\prime(x^\prime,p^\prime,t^\prime)=\cosh(\phi)H(x,p,t)-\sinh(\phi)p \\ p^\prime=\cosh(\phi)p-\sinh(\phi)H(x,p,t) \\ t^\prime=\cosh(\phi)t-\sinh(\phi)x \\ x^\prime=\cosh(\phi)x-\sinh(\phi)t. $$
We need to prove that Lorentz transformations are canonical. To that end we use the odd-dimensional formulation of Hamiltonian dynamics. Phase space has the coordinates $(x,p,t)$, and the Poincaré-Cartan form is $$ \theta=pdx-H(x,p,t)dt. $$ We use the extended definition of canonical transformations, whereby the transformation $$ x^\prime=x^\prime(x,p,t) \\ p^\prime=p^\prime(x,p,t) \\ t^\prime=t^\prime(x,p,t) $$ is canonical if $$ \theta^\prime=p^\prime dx^\prime-H^\prime(x^\prime,p^\prime,t^\prime)dt^\prime=\theta+dW $$for some function $W$.
Inserting the Lorentz transformation explicitly gives $$ \theta^\prime=\left(\cosh(\phi)p-\sinh(\phi)H\right)\left(\cosh(\phi)dx-\sinh(\phi)dt\right) \\ -\left( \cosh(\phi)H-\sinh(\phi)p \right)\left(\cosh(\phi)dt-\sinh(\phi)dx\right) \\ =pdx-Hdt=\theta, $$ i.e. the Poincaré-Cartan form is invariant, thus Lorentz transformations are canonical.
We are essentially done, but to show the invariance of volume more explicitly, the symplectic potential (canonical $1$-form) for a hypersurface $\sigma:\Sigma\rightarrow P$ ($P$ is the odd-dimensional phase space) is $$ \theta_\Sigma=\sigma^\ast\theta, $$ and the symplectic form is $$ \omega_\Sigma=d\theta_\Sigma. $$
Since the equal-time phase space is two dimensional, the symplectic form is the Liouville volume form.
The equal-time phase space for the unprimed observer is given by $t=c$, while the equal-time phase space for the primed observer is $t^\prime=c$. We obtain for the symplectic forms $$ \theta_{t=c}=pdx,\quad \omega_{t=c}=d\theta_{t=c}=dp\wedge dx, $$ and $$ \theta_{t^\prime=c}=p^\prime dx^\prime,\quad \omega_{t^\prime=c}=d\theta_{t^\prime=c}=dp^\prime\wedge dx^\prime, $$ which are of the same form, therefore the phase space volume is Lorentz invariant.
A: Phase space is Lorentz invariant. You can prove this by writing an integral over $d^4$ x $d^4$ p  and doing a Lorentz transformation, but there's a nice short proof in Padmanabhan's book "Gravitation", p.26 : 
For an observer moving with four-velocity $u_i$, the proper three-volume element
is given by $d^3V = u_0d^3x$ which is a scalar invariant. To prove this, note that the
quantity $d^4V = dx \ dy \ dz \ dt$ is a scalar. Multiplying this by $1 = u_0\frac{d\tau}{dt}$ and
noting that $d\tau$ is invariant, we conclude that $d^3V = u_0d^3x$ is an invariant.
