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!)