Date of last update at top of file.
The partition function
We start with the N-particle partition function: $$ Z_N = \int_{-\infty}^\infty dp_0 ... dp_N \int_0^\infty dx_0 \int_0^{x_0} dx_1 ... dx_N \thinspace \delta(E - H(x,p)) $$
with the Hamiltonian $$ H(x, p) = Mgx_0 + \frac{1}{2M}p_0^2 + \sum_i mgx_i + \sum_i \frac{1}{2m}p_i^2 $$
To make our lives easier we'll want to nondimensionalize this. Start by taking an \( E \) out of the delta function, remembering
$$ \delta\left(E - H(x,q)\right) = \frac{1}{E} \delta\left( \frac{Mg}{E}x_0 + \frac{1}{2ME}p_0^2 + \sum_i \frac{mg}{E} x_i + \sum_i \frac{1}{2mE}p_i^2\right) $$
Then change to dimensionless variables:
$$ \begin{align*} p_0'^2 &= \frac{1}{2ME}p_0^2 \implies dp_0 = (2ME)^\frac{1}{2} dp_0' \\ p_i'^2 &= \frac{1}{2mE}p_i^2 \implies dp_i = (2mE)^\frac{1}{2} dp_i'\\ \ x_0' &= \frac{Mg}{E}x_0 \implies dx_0 = \frac{E}{Mg} x_0' \\ x_i' &= \frac{mg}{E}x_0 \implies dx_i = \frac{E}{mg} dx_i' \end{align*} $$
This leaves:
$$ Z_N = \frac{1}{E}(2mE)^\frac{N}{2}(2ME)^\frac{1}{2} \left( \frac{E}{Mg} \right) \left( \frac{E}{mg} \right)^N \int_{-\infty}^\infty dp_0' ... dp_N' \int_0^\infty dx_0' \int_0^{\frac{m}{M}x_0'} dx_1' ... dx_N' \thinspace \delta\left(1 - x_0' - \sum_i x_i - p_0'^2 - \sum p_i'^2\right) $$
Define the total dimensionless momentum:
$$ p^2 = p_0'^2 + \sum p_i'^2 $$
To simplify again:
$$ Z_N = \frac{1}{E}(2mE)^\frac{N}{2}(2ME)^\frac{1}{2} \left( \frac{E}{Mg} \right) \left( \frac{E}{mg} \right)^N \int_{-\infty}^\infty dp_0' ... dp_N' \int_0^\infty dx_0' \int_0^{\frac{m}{M}x_0'} dx_1' ... dx_N' \thinspace \delta(1 - x_0' - \sum_i x_i - p^2) $$
Give the prefactor a name:
$$ K_N = \frac{1}{E}(2ME)^\frac{1}{2} \left( \frac{E}{Mg} \right) \left( \frac{E}{mg} \right)^N (2mE)^\frac{N}{2} $$
Rewrite:
$$ Z_N = K_N \int_{-\infty}^\infty dp_0' ... dp_N' \int_0^\infty dx_0' \int_0^{\frac{m}{M}x_0'} dx_1' ... dx_N' \thinspace \delta \left( 1 - x_0' - \sum_i x_i - p^2\right) $$
Since the \( p \) integral is spherically symmetric,
$$ \int_{-\infty}^\infty dp_0' ... dp_N' = S_N \int_0^\infty dp \thinspace p^N $$
Therefore:
$$ \boxed{Z_N = K_N S_N \int_0^\infty dp \thinspace p^N \int_0^\infty dx_0' \int_0^{\frac{m}{M}x_0'} dx_1' ... dx_N' \thinspace \delta \left( 1 - x_0' - \sum_i x_i - p^2\right)} $$
This is the first result.
Computing Z
We'll do this a few ways. First let's do it starting with the momentum. We need to evaluate
$$ \int_0^\infty dp \thinspace p^N \delta(1 - x_0' - \sum_i x_i - p^2) $$
Define:
$$ c^2 = 1 - x_0' - \sum_i x_i' $$
With this substitution,
$$ \int_0^\infty dp \thinspace p^N \delta(p^2 - c^2) $$
Let me show how this can be done. Do a \( u \) substitution on \( p \),
$$ \delta(p^2 - c^2) = \delta(u - c^2) $$
Then
$$ p = \pm \sqrt{u} \implies dp = \pm \frac{du}{2\sqrt{u}} $$
This is a doubly-branched integral, but we take only the positive branch, since \( p \) is always positive. Therefore \( dp = du/(2\sqrt{u})\), and the integration limits stay the same:
$$ \int_0^\infty dp \thinspace p^N \delta(p^2 - c^2) = \frac{1}{2} \int_0^\infty du \frac{1}{\sqrt{u}}u^\frac{N}{2} \delta(u - c^2) = \frac{1}{2} \int_{-\infty}^\infty du \thinspace \theta(u) u^\frac{N-1}{2} \delta(u - c^2) = \frac{1}{2} c^{N-1} \theta(c^2) $$
or
$$ \frac{1}{2} \left(1 - x_0' - \sum_i x_i'\right)^\frac{N-1}{2} \theta\left(1 - x_0' - \sum_i x_i'\right) $$
This leaves, for \( Z_N \),
$$ Z_N = \left( \frac{1}{2} \right) K_N S_N \int_0^\infty dx_0' \int_0^{\frac{m}{M}x_0'} dx_1' ... dx_N' \thinspace \left(1 - x_0' - \sum_i x_i' \right)^\frac{N-1}{2} \theta\left(1 - x_0' - \sum_i x_i'\right) $$
We can't do the \( x_0 \) integral until we do the others. Let's start with \( x_N' \).
$$ Z_N = \left( \frac{1}{2} \right) K_N S_N \int_0^\infty dx_0' \int_0^{\frac{m}{M}x_0'} dx_1' ... dx_N' \left(1 - x_0' - x_1' - ... - x_N' \thinspace \right)^\frac{N-1}{2} \theta\left(1 - x_0' - x_1' - ... - x_N'\right) $$
Rewrite this as
$$ Z_N = \left( \frac{1}{2} \right) K_N S_N \int_0^\infty dx_0' \int_0^{\frac{m}{M}x_0'} dx_1' ... dx_{N-1} \int_0^{\frac{m}{M}x_0'} dx_N' \thinspace \left(q - x_N' \right)^\frac{N-1}{2} \theta\left(q - x_N'\right) $$
where
$$ q = 1 - x_0' - x_1' - ... - x_{N-1}'. $$
Let's do this first integral.
$$ \int_0^{\frac{m}{M}x_0'} dx_N' \thinspace \left(q - x_N' \right)^\frac{N-1}{2} \theta\left(q - x_N'\right) $$
Change variables \( u = q - x_N'\), so that \( du = - dx_N'\), $$ \int_{{q-\frac{m}{M}x_0'}}^{{q}} du \thinspace u^\frac{N-1}{2} \theta(u) = \left( \frac{2}{N+1} \right) u^\frac{N+1}{2} \theta(u)\Big|_{{q-\frac{m}{M}x_0'}}^{q} = \frac{2}{N+1} \left( q^\frac{N+1}{2}\theta({q}) - \left(q - \frac{m}{M}x_0'\right) ^\frac{N+1}{2}\theta\left({q - \frac{m}{M}x_0'}\right) \right) $$
Resubbing in \( q \),
$$ \frac{2}{N+1} \left[ \left(1 - x_0' - x_1' - ... - x_{N-1}'\right)^\frac{N+1}{2}\theta \left( 1 - x_0' - x_1' - ... - x_{N-1}'\right) - \left(1 - x_0' - x_1' - ... - x_{N-1}' - \frac{m}{M}x_0'\right) ^\frac{N+1}{2}\theta\left(1 - x_0' - x_1' - ... - x_{N-1}' - \frac{m}{M}x_0'\ \right) \right] $$
Leaves us with this mess:
$$ Z_N = \frac{1}{2} K_N S_N \left( \frac{2}{N+1}\right) \int_0^\infty dx_0' \int_0^{\frac{m}{M}x_0'} dx_1' ... dx_{N-1}' \left( (1 - x_0' - x_1' - ... - x_{N-1}')^\frac{N+1}{2} \theta \left( 1 - x_0' - x_1' - ... - x_{N-1}'\right) - \left(1 - x_0' - x_1' - ... - x_{N-1}' - \frac{m}{M}x_0'\right)^\frac{N+1}{2}\theta \left( 1 - x_0' - x_1' - ... - x_{N-1} - \frac{m}{M}x_0'\right)\right) $$
Let's go to \( x_{N-1} \) and see where that gets us.
So the integral we need to do is: $$ \int_0^{\frac{m}{M}x_0'} dx_{N-1}' \left( (1 - x_0' - x_1' - ... - x_{N-1}')^\frac{N+1}{2} \theta(...) - \left(1 - x_0' - x_1' - ... - x_{N-1}' - \frac{m}{M}x_0'\right)^\frac{N+1}{2} \theta(...) \right) $$
If you do the same thing as above, you get:
$$ Z_N = \frac{1}{2} K_N S_N \left( \frac{2}{N+1}\right) \left( \frac{2}{N+3 }\right) \int_0^\infty dx_0' \int_0^{\frac{m}{M}x_0'} dx_1' ... dx_{N-2}' \left[ (1 - x_0' - ... - x_{N-2}')^\frac{N+3}{2}\theta(...) - 2 \left( 1 - x_0' - ... - x_{N-2} - \frac{m}{M}x_0'\right)^\frac{N+3}{2}\theta(...) + \left( 1 - x_0' - ... - x_{N-2}- 2\frac{m}{M}x_0'\right)^\frac{N+3}{2}\theta(...)\right] $$
Where the thing inside the theta is always the argument of what it's multiplying. (This pattern will continue.)
If you think about how we got to these integrals, you'll see that they form Pascal's triangle. If we repeat what we just did \( N \) times, we'll get:
$$ Z_N = \frac{1}{2} K_N S_N \left[ \prod_{k=0}^{N-1} \left( \frac{2}{N + 2k + 1}\right) \right] \int_0^{\infty} dx_0' \thinspace \sum_{k=0}^{N} \binom{N}{k}(-1)^k \left( 1 - x_0' - k \frac{m}{M}x_0'\right)^{\frac{3N-1}{2}} \theta(...) $$
Let's do the integral,
$$ \int_0^\frac{1}{\lambda} dx_0' \thinspace \left( 1 - x_0' - k \frac{m}{M}x_0'\right)^{\frac{3N-1}{2}} $$
Change variables \( u = 1 - \lambda x_0' \),
$$ \int_0^\frac{1}{\lambda} dx_0' \left( 1 - \lambda x_0'\right)^{\frac{3N-1}{2}} = \frac{1}{\lambda} \left( \frac{2}{3N + 1} \right) = \left( \frac{2}{3N +1} \right) \left( \frac{1}{1+ \frac{m}{M}k} \right) $$
Putting this back in \( Z \),
$$ Z_N = \frac{1}{2} K_N S_N \left( \frac{2}{3N +1} \right) \left[ \prod_{k=0}^{N-1} \left( \frac{2}{N + 2k + 1}\right) \right] \sum_{k=0}^{N} \binom{N}{k}(-1)^k \left( \frac{1}{1+ \frac{m}{M}k} \right) $$
A result that we'll need (and that I got from WolframAlpha) is:
$$ \sum_{k=0}^{N} \binom{N}{k}(-1)^k \left( \frac{1}{1+ \frac{m}{M}k} \right) = \frac{ \Gamma\left(N+1\right) \Gamma\left(\frac{M}{m} + 1\right)} {\Gamma\left(N + \frac{M}{m} + 1\right) } $$
Another one is $$ \prod_{k=0}^{N} \left( \frac{2}{N + 2k + 1}\right) = \frac{\Gamma((N+1)/2)}{\Gamma(3(N+1)/2)} $$
So, after simplifying:
$$ \boxed{Z_N = \frac{1}{2} K_N S_N \frac{ \Gamma\left(N+1\right) \Gamma\left(\frac{M}{m} + 1\right)} {\Gamma\left(N + \frac{M}{m} + 1\right) } \frac{\Gamma((N+1)/2)}{\Gamma(3(N+1)/2)}} $$
Checking our work by going the other way
Let's start by integrating over the particle positions this time. We'll refer to the original non-dimensionalized partition function that we derived earlier.
$$ Z_N = K_N S_N \int_0^\infty dp \thinspace p^N \int_0^\infty dx_0' \int_0^{\frac{m}{M}x_0'} dx_1' ... dx_N' \thinspace \delta(1 - x_0' - \sum_i x_i - p^2) $$
Let's evaluate the \( x \) integrals one at a time. For the first, one, \( x_1' \), we need to do
$$ \int_0^{\frac{m}{M}x_0'}dx_1'\thinspace\delta(q_1-x_1') $$
where \( q_1 = 1 - x_0' - p^2 - x_2' - ... - x_N') \). With a \( u \) substitution this becomes
$$ \int_{q_1-\frac{m}{M}x_0'}^{q_1}du \thinspace\delta(u) = \theta(u)\Big|_{q_1-\frac{m}{M}x_0'}^{q_1} = \theta\left(q_1\right) - \theta\left(q_1-\frac{m}{M}x_0'\right) $$
For completeness, let's do \(x_2' \). Let \( q_2 - x_2' = q_1 \). We repeat the process:
$$ \int_0^{\frac{m}{M}x_0'}dx_2'\thinspace \left[ \theta\left(q_2 - x_2' \right) - \theta\left(q_2 - x_2 ' -\frac{m}{M}x_0'\right) \right] $$
Again by \( u \) substitution and with the identity
$$ \int_a^b du \thinspace u^n \theta(u) = \frac{u^{n+1}}{n+1} \theta(u) \Big|_a^b $$
(in this case, \(n = 1\) but in the later cases we'll need the general form), we have
$$ \int_0^{\frac{m}{M}x_0'}dx_2'\thinspace \left( \theta\left(q_2 - x_2' \right) - \theta\left(q_2 - x_2 ' -\frac{m}{M}x_0'\right) \right) = q_2 \theta(q_2) - 2 \left(q_2 - \frac{m}{M}x_0' \right)\theta\left(q_2 - \frac{m}{M}x_0' \right) +\left(q_2 - 2\frac{m}{M}x_0' \right)\theta\left(q_2 - 2\frac{m}{M}x_0' \right) $$
We see Pascal again! The general formula for all \( N \) integrals is then
$$ \int_0^{\frac{m}{M}x_0'} dx_1' ... dx_N' \thinspace \delta(1 - x_0' - \sum_i x_i - p^2) = \frac{1}{(N-1)!}\sum_{k=0}^{N}\binom{N}{k}(-1)^k\left(1-p^2 - x_0' - \frac{m}{M}k x_0'\right)^{N-1} \theta\left(1-p^2 - x_0' - \frac{m}{M}k x_0'\right) $$
We have a choice for how we proceeed and in the end it doesn't matter too much which integral we do first, so I'll do the \( x_0' \) one. Let's simplify our lives a little:
$$ \int_0^\infty dx_0' \left(1-p^2 - x_0' - \frac{m}{M}k x_0'\right)^{N-1} \theta\left(1-p^2 - x_0' - \frac{m}{M}k x_0'\right) = \int_a^\infty dx \left(a - bx\right)^{N-1} \theta\left(a - bx\right) $$
Doing the integral:
$$ \begin{align*} \int_a^\infty dx \left(a - bx\right)^{N-1} \theta\left(a - bx\right) &= \frac{1}{b}\int_{-\infty}^{a}du \thinspace u^{N-1} \theta(u) \ &= \frac{1}{b}\frac{u^N}{N} \theta(u) \Big|_{-\infty}^{a} \ &= \frac{1}{b} \frac{a^N}{N}\theta(a) \ &= \frac{1}{1+\frac{m}{M}k} \frac{(1-p^2)^N}{N}\theta(1-p^2) \end{align*} $$
Plugging this back in,
$$ Z_N = K_N S_N \int_0^\infty dp \thinspace p^N (1-p^2)^N \theta(1-p^2) \left[ \frac{1}{N!} \sum_{k=0}^N\binom{N}{k}(-1)^k\frac{1}{1 + \frac{m}{M}k} \right] $$
Again, we can simplify the sum with
$$ Z_N = K_N S_N \frac{ \Gamma\left(\frac{M}{m} + 1\right)} {\Gamma\left(N + \frac{M}{m} + 1\right) } \int_0^\infty dp \thinspace p^N (1-p^2)^N \theta(1-p^2) $$
The thing on the right is a secret beta function,
$$ \int_0^1 dp \thinspace p^N (1-p^2)^N = \frac{1}{2} \int_0^1 dt \thinspace t^\frac{N-1}{2} (1-t)^N = \frac{1}{2} \frac{\Gamma((N+1)/2) \Gamma(N+1)}{\Gamma(3(N+1)/2)} $$
Since \(\Gamma(N+1) = N!\), we get a cancellation, and
$$ \boxed{Z_N = \frac{1}{2} K_N S_N \frac{\Gamma\left((N+1)/2\right)}{\Gamma(3(N+1)/2)} \frac{ \Gamma\left(N+1\right) \Gamma\left(\frac{M}{m} + 1\right)} {\Gamma\left(N + \frac{M}{m} + 1\right) } } $$
Position of piston
The marginal distribution is everything except for \(x_0 \).
$$ Z_N \rho(x_0) = \frac{1}{2} K_N S_N \left[ \prod_{k=0}^{N-1} \left( \frac{2}{N + 2k + 1}\right) \right] \left(\frac{Mg}{E}\right) \sum_{k=0}^{N} \binom{N}{k}(-1)^k \left( 1 - \frac{Mg}{E}x_0 - k \frac{m}{M} \frac{Mg}{E} x_0\right)^{\frac{3N-1}{2}} $$
Canceling out a lot of stuff on the left,
$$ \left( \frac{2}{3N + 1} \right) \frac{ \Gamma\left(N+1\right) \Gamma\left(\frac{M}{m} + 1\right)} {\Gamma\left(N + \frac{M}{m} + 1\right) } \rho(x_0) = \left(\frac{Mg}{E}\right) \sum_{k=0}^{N} \binom{N}{k}(-1)^k \left( 1 - \frac{Mg}{E}x_0 - k \frac{m}{M} \frac{Mg}{E} x_0\right)^{\frac{3N-1}{2}} $$
We identify
$$ \rho(x_0) = \left( \frac{3N + 1}{2} \right) \frac{\Gamma\left(N + \frac{M}{m} + 1\right)}{\Gamma\left(N+1\right) \Gamma\left(\frac{M}{m} + 1\right)} \left(\frac{Mg}{E}\right) \sum_{k=0}^{N} \binom{N}{k}(-1)^k \left( 1 - \frac{Mg}{E}x_0 - k \frac{m}{M} \frac{Mg}{E} x_0\right)^{\frac{3N-1}{2}} $$
Now let's try to calculate:
$$ \langle x_0 \rangle = \int dx_0 \thinspace x_0 \rho(x_0) $$
We'll need this integral
$$ \int_0^\frac{1}{a} dx \left[ \left( 1 - \left( \frac{Mg + kmg}{E}\right)x\right)^\frac{3N -1}{2} x \right] = \left( \frac{2}{3N + 1}\right)\left( \frac{2}{3N + 3}\right)\left( \frac{E}{Mg}\right)^2\left( \frac{1}{1 + \frac{m}{M}k}\right)^2 $$
Sticking this back in,
$$ \langle x_0 \rangle = \left( \frac{3N + 1}{2} \right) \frac{\Gamma\left(N + \frac{M}{m} + 1\right)}{\Gamma\left(N+1\right) \Gamma\left(\frac{M}{m} + 1\right)} \left(\frac{Mg}{E}\right) \sum_{k=0}^{N} \binom{N}{k}(-1)^k \left( \frac{2}{3N + 1}\right)\left( \frac{2}{3N + 3}\right)\left( \frac{E}{Mg}\right)^2\left( \frac{1}{1 + \frac{m}{M}k}\right)^2 $$
and simplifying,
$$ \langle x_0 \rangle = \left( \frac{2}{3N + 3}\right)\frac{\Gamma\left(N + \frac{M}{m} + 1\right)}{\Gamma\left(N+1\right) \Gamma\left(\frac{M}{m} + 1\right)} \left(\frac{E}{Mg}\right) \sum_{k=0}^{N} \binom{N}{k}(-1)^k \left( \frac{1}{1 + \frac{m}{M}k}\right)^2 $$
Do the sum,
$$ \sum_{k=0}^{N} \binom{N}{k}(-1)^k \left( \frac{1}{1 + \frac{m}{M}k}\right)^2 = \left(\frac{M}{m}\right)\frac{\Gamma(\frac{M}{m} + 1)\Gamma{(N+1)}}{\Gamma(N+\frac{M}{m} + 1)}\left( \psi^{(0)}\left(N + \frac{M}{m} + 1\right) - \psi^{(0)}\left(\frac{M}{m}\right)\right) $$
and plug this back in. All the gammas cancel:
$$ \langle x_0 \rangle = \frac{2}{3}\left( \frac{E}{Mg}\right) \left( \frac{1}{N+1}\right)\left(\frac{M}{m}\right) \left[ \psi\left(N + \frac{M}{m} + 1\right) - \psi\left(\frac{M}{m}\right)\right] $$
Here, \(\psi \) is the digamma function.
$$ \psi(z) = \frac{d}{dz} \ln \Gamma(z). $$
This turns out to be equivalent to
$$ \langle x_0 \rangle =\frac{2}{3}\left( \frac{E}{Mg}\right) \left( \frac{1}{N+1}\right) \left(\frac{M}{m}\right)\sum_{k=1}^N\frac{1}{1+\frac{m}{M}k} $$
Checking this result by going the other way
Let's start, this time, with
$$ Z_N = K_N S_N \frac{1}{(N-1)!} \sum_{k=0}^{N}\binom{N}{k}(-1)^k\int_0^\infty dp \thinspace p^N \int_0^\infty dx_0' \thinspace \left(1-p^2 - x_0' - \frac{m}{M}k x_0'\right)^{N-1} \theta\left(1-p^2 - x_0' - \frac{m}{M}k x_0'\right) $$
We've already computed this entirely. What we want is the merginal distribution \( \rho(x_0) \). Since \( x_0' = Mgx_0/E \), we can just replace it back in the equation above.
$$ Z_N \rho(x_0) = K_N S_N \frac{1}{(N-1)!} \sum_{k=0}^{N}\binom{N}{k}(-1)^k\int_0^\infty dp \thinspace p^N \left( \frac{Mg}{E} \right) \thinspace \left(1-p^2 - \frac{Mg}{E}x_0 - \frac{mgk}{E}x_0 \right)^{N-1} \theta\left(1-p^2 - \frac{Mg}{E}x_0 - \frac{mgk}{E}x_0\right) $$
Notice the lack of an integral over \( x_0 \). When you integrate, this will sum to one. Let's go ahead and do this integral. We'll make a few substitutions for good QoL. We'll call
$$ a^2 = 1 - \frac{Mg}{E}x_0 - \frac{mgk}{E}x_0 $$
so that
$$ Z_N \rho(x_0) = K_N S_N \left( \frac{Mg}{E} \right) \frac{1}{(N-1)!} \sum_{k=0}^{N}\binom{N}{k}(-1)^k\int_0^\infty dp \thinspace p^N \thinspace \left(a^2-p^2 \right)^{N-1} \theta\left(a^2-p^2\right) $$
So we need to evaluate the integral. Looking at it, you already know it's gonna be a beta function. Let's change variables to \( p = a t \) so that
$$ a^{3N-1} \int_0^1 dt \thinspace t^N \thinspace \left(1-t^2 \right)^{N-1} \theta(a^2 - a^2t^2)= \frac{1}{2} \left(1 - \frac{Mg}{E}x_0 - \frac{mgk}{E}x_0 \right)^\frac{3N-1}{2} \frac{\Gamma(N)\Gamma((N+1)/2)}{\Gamma((3N+1)/2)} \theta(...) $$
Putting this all together (and rembering the gamma cancels the factorial),
$$ Z_N \rho(x_0) = \frac{1}{2} K_N S_N \left( \frac{Mg}{E} \right) \frac{\Gamma((N+1)/2)}{\Gamma((3N+1)/2)} \sum_{k=0}^{N}\binom{N}{k}(-1)^k \left(1 - \frac{Mg}{E}x_0 - \frac{mgk}{E}x_0 \right)^\frac{3N-1}{2} \theta\left(1 - \frac{Mg}{E}x_0 - \frac{mgk}{E}x_0 \right) $$
It's a bit busy, but let's remember that we still have to cancel with what's on the left hand side,
$$ Z_N = \frac{1}{2} K_N S_N \frac{\Gamma\left((N+1)/2\right)}{\Gamma(3(N+1)/2)} \frac{ \Gamma\left(N+1\right) \Gamma\left(\frac{M}{m} + 1\right)} {\Gamma\left(N + \frac{M}{m} + 1\right) } $$
I must admit, some of those cancellations look awfully convenient.
$$ \boxed{\rho(x_0) = \left( \frac{Mg}{E} \right) \left( \frac{3N+1}{2}\right) \frac{\Gamma\left(N + \frac{M}{m} + 1\right)}{\Gamma\left(N+1\right) \Gamma\left(\frac{M}{m} + 1\right)}\sum_{k=0}^{N}\binom{N}{k}(-1)^k \left(1 - \frac{Mg}{E}x_0 - \frac{mgk}{E}x_0 \right)^\frac{3N-1}{2} \theta\left(1 - \frac{Mg}{E}x_0 - \frac{mgk}{E}x_0 \right) } $$
That looks surprisingly right, and it integrates to one (I'm pretty sure, but ya never know!)
The rest of the calculation is the same, so we can say with some confidence
$$ \boxed{\langle x_0 \rangle = \frac{2}{3}\left( \frac{E}{Mg}\right) \left( \frac{1}{N+1}\right)\left(\frac{M}{m}\right) \left[ \psi\left(N + \frac{M}{m} + 1\right) - \psi\left(\frac{M}{m}\right)\right]} $$
Momentum distribution function
Oy gevalt. Ok.
Let's start with
$$ Z \rho(p_N) = \int_{-\infty}^\infty dp_0 ... dp_{N-1} \int_0^\infty dx_0 \int_0^{x_0} dx_1 ... dx_N \thinspace \delta(E - H(p,q)) $$
We can go ahead and do the integrals over all the positions just like before. $$ Z_N \rho(p_N) = K_N \frac{1}{(2mE)^\frac{1}{2}} S_{N-1} \frac{ \Gamma\left(\frac{M}{m} + 1\right)} {\Gamma\left(N + \frac{M}{m} + 1\right) } \int_0^\infty dp \thinspace p^{N-1} \left(1- \frac{p_N^2}{2mE} - p^2\right)^{N} \theta\left(1- \frac{p_N^2}{2mE} - p^2\right) $$
Let's do a substitution on the integral again,
$$ \begin{align*} \int_0^\infty dp \thinspace p^{N-1} \left(1- \frac{p_N^2}{2mE} - p^2\right)^{N} \theta\left(1- \frac{p_N^2}{2mE} - p^2\right) &= \int_0^\infty dp \thinspace p^{N-1} \left(a^2 - p^2\right)^{N} \theta\left(a^2 - p^2\right) \\ &= a^{3N} \int_0^\infty dt \thinspace t^{N-1}(1-t^2)^N\theta(a^2(1-t^2)) \\ &= \frac{1}{2} \frac{\Gamma(N/2) \Gamma(N + 1)}{\Gamma((3N + 2)/2)} \left( 1- \frac{p_N^2}{2mE} \right)^\frac{3N}{2} \theta\left( 1- \frac{p_N^2}{2mE} \right) \end{align*} $$
That gives:
$$ Z_N \rho(p_N) = \frac{1}{2} K_N \frac{1}{(2mE)^\frac{1}{2}} S_{N-1} \frac{ \Gamma\left(\frac{M}{m} + 1\right)} {\Gamma\left(N + \frac{M}{m} + 1\right) } \frac{\Gamma(N/2) \Gamma(N + 1)}{\Gamma((3N + 2)/2)} \left( 1- \frac{p_N^2}{2mE} \right)^\frac{3N}{2} \theta\left( 1- \frac{p_N^2}{2mE} \right) $$
Dividing through (this took me a while! like 45 minutes to check everything),
$$ \boxed{\rho(p) = \frac{1}{\sqrt{2 \pi m E}} \frac{\Gamma(3N/2+3/2)}{\Gamma(3N/2 + 1)} \left( 1- \frac{p^2}{2mE} \right)^\frac{3N}{2} \theta\left( 1- \frac{p^2}{2mE} \right)} $$
This integrates to one.
Piston total microstate weight
Considering the zero gravity case, with plans to generalize later.
$$ \begin{align*} W(p_0, x_0) &= \frac{1}{E} \left( 2 m E\right)^\frac{N}{2} x_0^N \int_{-\infty}^\infty dp_1' ... dp_N' \thinspace \thinspace \delta\left( 1 - \frac{Mg}{E} x_0 - \frac{1}{2ME} p_0^2 - p^2 \right) \\ &= \frac{1}{E} \left( 2 m E\right)^\frac{N}{2} x_0^N S_{N-1}\int_{0}^\infty dp \thinspace p^{N-1} \thinspace \thinspace \thinspace \delta\left( 1 - \frac{Mg}{E} x_0 - \frac{1}{2ME} p_0^2 - p^2 \right) \\ &= \frac{1}{E} \left( 2 m E\right)^\frac{N}{2} x_0^N S_{N-1} \frac{1}{2} \int_{0}^\infty du \thinspace u^\frac{N-2}{2} \thinspace \thinspace \thinspace \delta\left( c^2 - u \right) \\ &= \frac{1}{2E} \left( 2 m E\right)^\frac{N}{2} x_0^N S_{N-1} \int_{-\infty}^\infty du \thinspace u^\frac{N-2}{2} \thinspace \thinspace \thinspace \delta\left( c^2 - u \right) \theta(u) \\ &= \frac{1}{2E} \left( 2 m E\right)^\frac{N}{2} x_0^N S_{N-1} c^{N-1} \theta(c^2) \\ &= \frac{1}{2E} \left( 2 m E\right)^\frac{N}{2} x_0^N S_{N-1} \left(1 - \frac{Mg}{E} x_0 - \frac{1}{2ME} p_0^2 \right)^\frac{N-1}{2} \theta\left( 1 - \frac{Mg}{E} x_0 - \frac{1}{2ME} p_0^2 \right) \\ % &= \frac{S_{N-1}}{2} \left( 2 m E\right)^\frac{N}{2} x_0^N \left(1 - \frac{Mg}{E} x_0 - \frac{1}{2ME} p_0^2 \right)^\frac{N-1}{2} \theta\left( E - Mg x_0 - \frac{1}{2M} p_0^2 \right) \\ \end{align*} $$
In terms of dimensionless variables, and up to an overall prefactor:
$$ W(\tilde{x}_0, \tilde{p}_0) = \lambda \tilde{x}_0^N \left(1 - \tilde{x}_0 - \tilde{p}_0^2 \right)^\frac{N-1}{2} \theta\left( 1 - \tilde{x}_0 - \tilde{p}_0^2 \right) $$
The relative entropy between two microstates \(A\) and \(B \) is then
$$ \ln \left(\frac{W(\tilde{x}{0,A}, \tilde{p}{0,A})}{W(\tilde{x}{0,B}, \tilde{p}{0,B})} \right) = N \ln\left( \frac{\tilde{x}{0,A}}{\tilde{x}{0,B} } \right) + \frac{N-1}{2}\ln \left(\frac{1 - \tilde{x}{0,A} - \tilde{p}{0,A}^2 }{1 - \tilde{x}{0,B} - \tilde{p}{0,B}^2}\right) $$
In the most likely microstate, \( x_{0,B} = 2/3 \) and \( p_{0,B}\) = 0, so relative to that some microstate with variables \( \tilde{x}_0 \) and \( \tilde{p}_0 \) has entropy
$$ \ln \left(\frac{W(\tilde{x}{0}, \tilde{p}{0})}{W_{max}} \right) = N \ln\left( \frac{3\tilde{x}{0}}{2 } \right) + \frac{N-1}{2}\ln \left(3(1 - \tilde{x}{0} - \tilde{p}_{0}^2) \right) $$
Particle height
This is an interesting one, because I don't know what to expect.
$$ Z \rho(x_N) = \int_{-\infty}^\infty dp_0 ... dp_{N} \int_0^\infty dx_0 \int_0^{x_0} dx_1 ... dx_{N-1} \thinspace \delta(E - H(p,q)) $$
Let's do the position integrals first as before, stopping short of going all the way to \( N \):
$$ Z_N \rho(p_N) = K_N \left( \frac{mg}{E} \right) S_{N} \frac{ \Gamma\left(1 + \frac{M}{m} \right)} {\Gamma\left(N + \frac{M}{m} \right) } \int_0^\infty dp \thinspace p^{N} \left(1- \frac{mg}{E}x_N - p^2\right)^{N-1} \theta\left(1- \frac{mg}{E}x_N - p^2\right) $$
Do a substitution as before:
$$ \begin{align*} \int_0^\infty dp \thinspace p^{N} \left(1- \frac{mg}{E}x_N - p^2\right)^{N-1} \theta\left(1- \frac{mg}{E}x_N - p^2\right) &= \int_0^\infty dp \thinspace p^{N} \left(a^2 - p^2\right)^{N-1} \theta\left(a^2 - p^2\right) \\ &= a^{3N-1}\int_0^\infty dt \thinspace t^N(1-t^2)^{N-1}\theta(a^2(1-t^2)) \\ &= \frac{1}{2} \frac{\Gamma(N) \Gamma((N + 1)/2)}{\Gamma((3N + 1)/2)} \left( 1- \frac{mg}{E}x_N \right)^\frac{3N-1}{2} \theta\left( 1- \frac{mg}{E}x_N \right) \\ \end{align*} $$
Plugging this stuff back in,
$$ Z \rho(x_N) = \frac{1}{2} K_N S_N \left( \frac{mg}{E} \right) \frac{ \Gamma\left(1 + \frac{M}{m} \right)} {\Gamma\left(N + \frac{M}{m} \right) } \frac{\Gamma(N) \Gamma((N + 1)/2)}{\Gamma((3N + 1)/2)} \left( 1- \frac{mg}{E}x_N \right)^\frac{3N-1}{2} \theta\left( 1- \frac{mg}{E}x_N \right) $$
Substituting for \( Z_N \),
$$ \rho(x_N) = \left(\frac{mg}{E}\right) \left( \frac{3N+1}{2} \right) \left( 1 + \frac{M}{Nm} \right)\left( 1- \frac{mg}{E}x_N \right)^\frac{3N-1}{2} \theta\left( 1- \frac{mg}{E}x_N \right) $$
Does this integrate to one? It doesn't. Hmph. Without that middle factor everything would work out, but I can't figure out how to make it go away.
Arbitrary potential, no gravity
Let's solve the system where we don't know the potential of the piston, but we do know that the gas is ideal and not subject to that potential. For convenience, we'll only remove the dimensions from the momenta:
$$ Z_N = \frac{1}{E}(2mE)^\frac{N}{2}(2ME)^\frac{1}{2} \int_{-\infty}^\infty dp_0' ... dp_N' \int_0^\infty dx_0 \int_0^{x_0} dx_1 ... dx_N \thinspace \delta\left(1 - p^2 - \frac{U(x_0)}{E} \right) $$
We can do the \( x_i \) integrals right away,
$$ Z_N = \frac{1}{E}(2mE)^\frac{N}{2}(2ME)^\frac{1}{2} \int_{-\infty}^\infty dp_0' ... dp_N' \int_0^\infty dx_0 \thinspace x_0^N \thinspace \delta\left(1 - p^2 - \frac{U(x_0)}{E} \right) $$
Now swap the order of integration,
$$ \begin{align*} Z_N &= \frac{1}{E}(2mE)^\frac{N}{2}(2ME)^\frac{1}{2}S_N \int_0^\infty dx_0 \thinspace x_0^N \int_{0}^\infty dp \thinspace p^N \delta\left(1 - \frac{U(x_0)}{E} - p^2 \right) \\ &= T_N S_N \int_0^\infty dx_0 \thinspace x_0^N \int_{0}^\infty dp \thinspace p^N \delta\left(1 - \frac{U(x_0)}{E} - p^2 \right)\\ &= T_N S_N \int_0^\infty dx_0 \thinspace x_0^N \int_{0}^\infty dp \thinspace p^N \delta\left(a^2 - p^2 \right) \\ &= \frac{1}{2} T_N S_N \int_0^\infty dx_0 \thinspace x_0^N \int_{0}^\infty du \thinspace u^\frac{N-1}{2} \delta\left(a^2 - u \right) \\ &= \frac{1}{2} T_N S_N \int_0^\infty dx_0 \thinspace x_0^N \int_{-\infty}^\infty du \thinspace u^\frac{N-1}{2} \delta\left(a^2 - u \right)\theta(u) \\ &= \frac{1}{2} T_N S_N \int_0^\infty dx_0 \thinspace x_0^N \left( 1 - \frac{U(x_0)}{E} \right)^\frac{N-1}{2} \theta\left( 1 - \frac{U(x_0)}{E} \right) \\ \end{align*} $$
For one, we can quickly identify
$$ \rho(x_0) = \frac{x_0^N \left( 1 - \frac{U(x_0)}{E} \right)^\frac{N-1}{2} \theta\left( 1 - \frac{U(x_0)}{E} \right)}{\int_0^\infty dx_0 \thinspace x_0^N \left( 1 - \frac{U(x_0)}{E} \right)^\frac{N-1}{2} \theta\left( 1 - \frac{U(x_0)}{E} \right)} $$
What about the mean height? One strategy might be to consider the following: in the thermodynamic limit, the distribution
$$ \rho(x_o) = x_0^N \left( 1 - \frac{U(x_0)}{E} \right)^\frac{N-1}{2} $$
has an extremely sharp peak around its mean value, kind of like a delta function. That's going to be where it's most likely to find the piston. So let's see how we can estimate that. One simple way is to write this as a pure exponential
$$ \begin{align*} x_0^N \left( 1 - \frac{U(x_0)}{E} \right)^\frac{N-1}{2} &= \exp\left( \left( \frac{N-1}{2} \right) \ln (1 - U(x_0)/E) + N \ln(x_0) \right) \\ &\approx \exp(N M(x_0)) \end{align*} $$
where
$$ M(x_0) = \frac{1}{2} \ln(1-U(x_0)/E) + N \ln(x_0) $$
As the exponential is monotone, it means when \( M \) is increasing, so is $\rho$, and same for when it's decreasing. So we can estimate the peak of the distribution by checking when the derivative of \( M \) is equal to zero,
$$ \frac{d M }{dx_0} = - \frac{1}{2} \frac{U'(x_0)/E}{1-U/E}+ \frac{1}{x} = 0 $$
Or
$$ \overline{x}_0 = 2 \left( \frac{1-U/E}{U'(\overline{x}_0 )/E} \right) $$
where the overbar represents the mean value. Let's see how this works for the potential \( U(x) = Mgx \),
$$ \begin{align*} \overline{x}_0 &= 2 \left( \frac{E}{Mg} \right) \left(1-\frac{Mg \overline{x}_0 }{E} \right) \\ &= 2 \left( \frac{E}{Mg} \right) - 2 \overline{x}_0 \\ \end{align*} $$
or
$$ \overline{x}_0 = \frac{2}{3}\frac{E}{Mg} $$
We can also do it with a harmonic potential, \( U(x) = \frac{1}{2} k x^2 \). Then
$$ \begin{align*} \overline{x}_0 &= 2 \left( \frac{1-U/E}{U'(\overline{x}_0 )/E} \right) \\ &= 2 \left( \frac{E}{k\overline{x}_0 } \right) \left(1-\frac{1}{2E}k\overline{x}_0 ^2 \right) \\ \end{align*} $$
which implies
$$ k\overline{x}_0 ^2 = 2E - k\overline{x}_0 ^2 \implies \overline{x}_0 = \sqrt{E/k} $$
In fact this works for any potential of the form \( U(x) = \frac{1}{n} \lambda x^n \),
$$ \begin{align*} \overline{x}_0 &= 2 \left( \frac{1-U/E}{U'(\overline{x}_0 )/E} \right) \\ &= 2 \left( \frac{E}{\lambda \overline{x}_0^{n-1} } \right) \left(1-\frac{1}{nE}\lambda x^n \right) \\ \end{align*} $$
which implies
$$ \lambda \overline{x}_0^n = 2E - \frac{2}{n} \lambda \overline{x}_0^n \implies \overline{x}_0 = \left( \frac{2E}{\lambda} \frac{n}{n + 2} \right)^\frac{1}{n} $$
for the mean height of the piston for any potential, as long as there are enough particles in the box.
Nondimensionalized Hamiltonian
Per Josh, we should try writing the hamiltonian in terms of nondimensional variables. $$ \begin{align*} \frac{H(x,p)}{E}&= \frac{Mg}{E}x + \sum_i \frac{mg'}{E} x_i + \frac{p_0^2}{2ME} + \sum_i \frac{p_i^2}{2mE} \\ &= \frac{x}{L} + \sum_i \frac{\mu}{L}x_i + \frac{p_0^2}{P^2} + \sum_i\frac{1}{\mu} \frac{p_i^2}{P^2} \end{align*} $$
where \( \mu = m/M \), \( \gamma = g'/g \), \(L = Mg/E \), and \( P^2 = 2ME \). Then we define the dimensionless variables with tildes.
$$ \begin{align*} \tilde{x}_0 &= \frac{Mg}{E} x_0 && = x_0/L \\ \tilde{x}_i &= \frac{mg}{E} x_i && = \mu x_i/L \\ \tilde{p}_0 &= \frac{p_0}{2ME} &&= p_0/P \\ \tilde{p}_i &= \frac{p_i}{2ME} &&= p_i/P \\ \end{align*} $$
This leaves the transformed hamiltonian as
$$ H(\tilde{x}, \tilde{p}) = E \tilde{x} + E \mu \sum_i \tilde{x}_i + E\tilde{p}_0^2 + \frac{E}{\mu} \sum_i \tilde{p}_i^2 $$
We need to find the equations of motion for this system. My strategy was to start from the hamiltonian equations that I knew to be correct, then just add the constants in where necessary. The steps are below:
$$ \begin{align*} \frac{\partial H}{\partial x_0} &= - \dot{p}_0 \implies \frac{1}{L} \frac{\partial H}{\partial \tilde{x}_0} = -P \dot{\tilde{p}_0} && \implies \frac{\partial H}{\partial \tilde{x}_0} = - LP \dot{\tilde{p}}_0 \implies E = - LP\dot{\tilde{p}}_0 \\ \frac{\partial H}{\partial p_0} &= \dot{x}_0 \implies \frac{1}{P} \frac{\partial H}{\partial \tilde{p}_0} = L \dot{\tilde{x}_0} && \implies \frac{\partial H}{\partial \tilde{p}_0} = LP \dot{\tilde{x}}_0 \implies 2E \tilde{p}_0 = LP\dot{\tilde{x}}_0 \\ \frac{\partial H}{\partial x_i} &= - \dot{p}_i \implies \frac{\mu}{L} \frac{\partial H}{\partial \tilde{x}_i} = -P \dot{\tilde{p}_i} && \implies \frac{\partial H}{\partial \tilde{x}_i} = - \frac{LP}{\mu} \dot{\tilde{p}}_i \implies E \mu = - \frac{LP}{\mu}\dot{\tilde{p}}_0 \\ \frac{\partial H}{\partial p_i} &= \dot{x}_i \implies \frac{1}{P} \frac{\partial H}{\partial \tilde{p}_i} = L \dot{\tilde{x}_i} && \implies \frac{\partial H}{\partial \tilde{p}_i} = LP \dot{\tilde{x}}_i \implies \frac{2E}{\mu} \tilde{p}_i = LP\dot{\tilde{x}}_i \\ \end{align*} $$
$$ \begin{align*} \dot{\tilde{p}_0} &= - \frac{E}{LP} &&= - \sqrt{\frac{g}{2L}} && \implies \tilde{p}_0(t) = \tilde{p}_0(0) - \sqrt{\frac{g}{2L}}t \\ \dot{\tilde{x}_0} &= \frac{2E}{LP} \tilde{p}_0 &&= \sqrt{\frac{2g}{L}}\tilde{p}_0 && \implies \tilde{x}_0(t) = \tilde{x}_0(0) + \sqrt{\frac{2g}{L}} \tilde{p}_0(0) t - \frac{g}{L} t^2 \\ \dot{\tilde{p}_i} &= 0 && && \implies \tilde{p}_i(t) = \text{const}\\ \dot{\tilde{x}_i} &= \frac{2E}{\mu LP} \tilde{p}_i &&= \frac{1}{\mu} \sqrt{\frac{2g}{L}}\tilde{p}_i && \implies \tilde{x}_i(t) = \tilde{x}_i(0) + \frac{1}{\mu}\sqrt{\frac{2g}{L}} \tilde{p}_i(0) t \\ \end{align*} $$
The elastic collision formulas, for a piston and particle, are:
$$ \begin{align*} p_0 + p_i &= p_0' + p_i' &&\implies \tilde{p_0} + \tilde{p_i} = \tilde{p_0}' + \tilde{p_i}' \\ \frac{p_0^2}{2M} + \frac{p_i^2}{2m} &= \frac{p_0'^2}{2M} + \frac{p_i'^2}{2m} &&\implies \tilde{p_0}^2 + \frac{1}{\mu}\tilde{p}_i^2 = \tilde{p_0}'^2 + \frac{1}{\mu}\tilde{p}_i'^2 \end{align*} $$
Now we need to solve these. This is the least fun part of my day.
$$ \begin{align*} \tilde{p}_0' &= \frac{2}{\mu + 1} \tilde{p}_i - \frac{\mu - 1}{\mu + 1}\tilde{p}_0 \\ \tilde{p}_i' &= \frac{\mu - 1}{\mu + 1} \tilde{p}_i + \frac{2 \mu }{\mu + 1}\tilde{p}_0 \\ \end{align*} $$
So this will allow us to rewrite the simulation just in terms of a few parameters. Let's think about how to do this. First let's remember that we want all length scales in the problem to be dimensionless. That means we'll measure positions in terms of the natural length scale \( L \).
Let's see how long it takes to hit the ground using the new variables:
$$ \tilde{x}_i(t_g) = 0 \implies t_g = - \sqrt{\frac{L}{2g}} \frac{\mu \tilde{x}}{\tilde{p}} $$
A kind of natural time scale appears here:
$$ T = \sqrt{\frac{L}{ 2 g}} $$
so that we can rewrite that time to ground as
$$ \tau_g = \frac{t}{T} =- \mu \tilde{x} / \tilde{p} $$
We can also recompute the time it takes to hit the piston:
$$ \tilde{x}_i(t_p) = \tilde{x}_0(t_p) \implies \tilde{x}_i + \frac{1}{\mu} \sqrt{\frac{2g}{L}}\tilde{p}_i t_p = \tilde{x}_0 + \sqrt{\frac{2g}{L}}\tilde{p}_0t_p-\frac{g}{2L}t_p^2 $$
or more simply, in terms of the natural time,
$$ \tilde{x}_i + \frac{1}{\mu} \tilde{p}_i \tau_p = \tilde{x}_0 + \tilde{p}_0 \tau_p -\frac{1}{4} \tau_p^2 $$
Solving for \( \tau_p = t_p / T \),
$$ \tau_p = 2 \left( \tilde{p}_0 - \frac{1}{\mu} \tilde{p}_i \right) + 2 \sqrt{ \left( \tilde{p}_0 - \frac{1}{\mu} \tilde{p}_i \right)^2 + (\tilde{x}_0 - \tilde{x}_i)} $$
Note that in the general case, that is, where the particle has its own \( g' \), the collision times are
$$ \begin{align*} \tau_p &= \frac{2 \left( \tilde{p}_0 - \frac{1}{\mu} \tilde{p}_i \right) + 2 \sqrt{ \left( \tilde{p}_0 - \frac{1}{\mu} \tilde{p}_i \right)^2 + (1 - \gamma) (\tilde{x}_0 - \tilde{x}_i)}}{1 - \gamma} \\ \tau_g &= \frac{2}{\mu \gamma} \tilde{p}_i + \frac{2}{\gamma} \sqrt{\left( \frac{1}{\mu} \tilde{p}_i \right)^2 + \gamma \tilde{x}} \end{align*} $$
The equations of motion in terms of \( \tau = t / T \),
$$ \begin{align*} \tilde{x}_0(\tau) &= \tilde{x}_0(0) + \tilde{p}_0(0) \tau - \frac{1}{4} \tau ^2 \\ \tilde{x}_i(\tau) &= \tilde{x}_i(0) + \frac{1}{\mu}\tilde{p}_i(0) \tau - \frac{\gamma}{4} \tau^2 \\ \tilde{p}_0(\tau) &= \tilde{p}_0(0) - \frac{1}{2} \tau \\ \tilde{p}_i(\tau) &= \tilde{p}_0(0) - \frac{\gamma \mu }{2} \tau \end{align*} $$
So now what we have is an entire kinematic system that only depends on \( \mu \), really. I guess if you have a \( \gamma \) then it also depends on that.
Let's look at the momentum distribution, which is independent of \( \gamma \) incidentally:
$$ \rho(p) = \rho(p) = \frac{1}{\sqrt{2 \pi m E}} \frac{\Gamma(3N/2 + 3/2)}{\Gamma(3N/2 + 1)} \left( 1- \frac{p^2}{2mE} \right)^\frac{3N}{2} \theta\left( 1- \frac{p^2}{2mE} \right) $$
This becomes, for dimensionless \( p \),
$$ \boxed{\rho(\tilde{p}) = \frac{1}{\sqrt{\pi \mu}} \frac{\Gamma(3N/2 + 3/2)}{\Gamma(3N/2 + 1)} \left( 1- \frac{1}{\mu}\tilde{p}^2 \right)^\frac{3N}{2} \theta\left( 1- \frac{1}{\mu}\tilde{p}^2 \right)} $$
Again, only depends on \(\mu\). Pretty cool.
Originally I had this equation wrong: I thought that the gamma functions differed by 1, and not 1/2, and wrote down a prefactor of \((3N+1)/2 \) (expecting something the same as for the piston.) But then when I integrated this numerically I got all the wrong answers. That caused me to go back and double check against Josh's notes, and when he didn't simplify, I realized my mistake.
The bad part about this is for simulation purposes, it's pretty impossible to compute (exactly) the ratio of these gammas, because each \( \Gamma( 100 ) \) is basically infinity. But a pretty good approximation for large \( x \) is that the ratio is about
$$ \frac{\Gamma(x + 1/2)}{\Gamma(x)} \approx \sqrt{x} $$ and that's what I use in the sim.
Graveyard
What follows are all unfinished notes.
Coupling
Let's integrate out the degrees of freedom of the piston from the partition function in the zero gravity case.
$$ \begin{align*} Z_N &= \int_{-\infty}^\infty dp_0 ... dp_N \int_0^\infty dx_0 \thinspace x_0^N \delta\left(E - Mgx_0 - \frac{p_0^2}{2M} + \sum_i \frac{p_i^2}{2m}\right) \\ &= \frac{1}{E} \int_{-\infty}^\infty dp_0 ... dp_N \int_0^\infty dx_0 \thinspace x_0^N \delta\left(1 - \frac{Mg}{E} x_0 - \frac{p_0^2}{2ME} + \sum_i \frac{p_i^2}{2mE}\right) \\ &= \frac{1}{E} \left(\frac{E}{mg}\right)^N \left(2ME\right)^\frac{1}{2}\left(2mE\right)^\frac{N}{2} \int_{-\infty}^\infty dp'_0 ... dp'_N \int_0^\infty dx'_0 \thinspace x_0'^N \delta\left(1 - x_0' - p_0'^2 + \sum_i p_i'^2\right) \\ &= K_N \int dp_0'...dp_N' \int dx_0' \thinspace x_0'^N \delta( x_0' - 1 + p^2) \theta(x_0') \\ &= K_N \int dp_0'...dp_N'\thinspace (1 - p^2)^N \theta(1-p^2) \\ &= K_N \int dp_1'...dp_N' \int dp_0' \thinspace \left( 1 - \sum_i p_i'^2 - p_0'^2 \right)^N \theta\left((1- \sum_i p_i'^2) - p_0'^2 \right) \\ &= K_N \int dp_1'...dp_N' \int dp_0' \thinspace \left( a^2 - p_0'^2 \right)^N \theta\left(a^2 - p_0'^2 \right) \\ &= K_N \int dp_1'...dp_N' a^{2N} \int dp_0' \thinspace \left(1 - (p_0'/a)^2 \right)^N \theta\left(a^2 - p_0'^2 \right) \\ \end{align*} $$
This last integral is just
$$ \int dp_0' \thinspace \left(1 - (p_0'/a)^2 \right)^N \theta\left(a^2 - p_0'^2 \right) = a \thinspace \theta(a^2) \int_{-1}^{1} dt \thinspace \left(1 - t^2 \right)^N = a \theta(a^2) \frac{\sqrt{\pi}\Gamma(N+1)}{\Gamma(N + 3/2)} $$
For a total partition function:
$$ Z_N = K_N\frac{\sqrt{\pi}\Gamma(N+1)}{\Gamma(N + 3/2)} \int dp_1'...dp_N' \left( 1 - \sum_i p_i'^2 \right)^{N+ \frac{1}{2}} \theta\left(1 - \sum_i p_i'^2 \right) $$
Grad-H effect
$$ H(x, p) = Mgx_0 + \frac{1}{2M}p_0^2 + \sum_i \frac{1}{2m}p_i^2 $$
and
$$ |v| = | \nabla H | = \left| \left(\frac{\partial H}{\partial x}\right)^2 + \left(\frac{\partial H}{\partial p}\right)^2 \right|^\frac{1}{2} $$
Therefore
$$ |v| = \left| Mg + p_0/M + \sum_i p_i/m \right|^\frac{1}{2} $$
If $$ \langle p_i \rangle = \frac{1}{N} \sum_i p_i $$
Then $$ |v| = \left| Mg + p_0/M + N \langle p_i \rangle /m \right|^\frac{1}{2} $$
So the velocity through phase space depends on \( G \), the average momentum of the particles and the momentum of the piston.
Lagrangian formalism
No gravity case:
$$ L = \frac{1}{2} M \dot{x_0}^2 + \sum_i \frac{1}{2} m \dot{x_i}^2 - Mgx_0 + \sum_i \gamma_i \theta(x_i) + \sum_i \lambda_i \theta (x_0 - x_i) $$
$$ \frac{d}{dt} \frac{\partial L}{\partial \dot{x_i}} - \frac{\partial L }{\partial x_i} = m \ddot{x_i} - \gamma_i \delta(x_i) + \lambda_i \delta(x_0 - x_i) = 0 $$
$$ \frac{d}{dt} \frac{\partial L}{\partial \dot{x_0}} - \frac{\partial L }{\partial x_0} = M \ddot{x_0} - \sum_i \lambda_i \delta(x_0 - x_i) = 0 $$
Kinetics
I want to think a little bit about how you might approach the problem of equilibrium from a kinetics point of view.
In perfect mechanical equilibrium the piston is exactly fixed. But per unit time \( dt \) it feels an impulse from both the particles and gravity. Let's see if we can see how these balance. The impulse from gravity is just \( M g dt \), but for the particles it's a bit trickier to calculate.
If the piston is at perfect mechanical equilibrium, then it basically functions as a wall, so the impulse from any particle is just equal to twice the momentum of that particle, assuming it hits the piston. So we just want to know how many particles hit the piston in time \( dt \) and what their momenta are.
The particles lie somewhere between \( 0 \) and \( x_0 \), where \( x_0 \) is the current height of the piston. Let's divide the particles into ones which are going forwards (i.e. torwards the piston, \( v > 0 \) ) and backwards (\( v < 0 \)). The particles going forwards, at position \( x \), are \( x_0 - x \) units of distance away from the piston. The ones going backwards are effectively \(2x + (x_0 - x) = x_0 + x \) away from the piston.
In a time \( dt \), which of these particles will hit the piston? This is actually a really tricky problem. If the particle velocity is totally unbounded, then in a given interval of time, some particles will hit once, some twice, and some arbitrarily many times. But we're working in the case that there's an upper bound for momentum, which simplifies things a little. We can take \( dt \) small enough so that in that time interval, the only particles which hit the piston are the ones in the immediate neighborhood of the piston and traveling towards it. That is, we can consider the impulses as a sum of the impulses,
$$ \int_{x_0 - v_{max} dt}^{x_0}dx \thinspace i(x) $$
But the impulse \( i(x) \) is just twice the momentum of the particles arriving from there. At \( x \) within \( dx \) of \(x_0 \) we can sum over all the particles, up to the maximum velocity.
$$ i(x_0)dx = 2 dx \int_0^{v_{max}} \rho(v) dv = dx $$
But if we go \( dx \) further away, we need to exclude the particles that no longer collide with the piston because their velocity is too low:
$$ i(x_0 - dx)dx = 2 dx \int_{\frac{dx}{dt}}^{v_{max}} \rho(v) dv = dx $$
In a given instant, the energy change of the piston is balanced by the energy change of the particles.
$$ \frac{N}{2} m \langle v \rangle^2 + Mg x_0 = E $$
Then
Entropy flow
$$ S = \int_\Omega dp dx \thinspace \rho(x, p) \ln\rho(x, p) $$
$$ \frac{dS}{dt} = \int_\Omega dp dx \thinspace \left[ \dot{\rho}\ln \rho + \dot{\rho}\right] $$
Arbitrary potential for both [unfinished]
We want to approximate the partition function for an arbitrary potential for the particles. Let's assume they're all in the same potential \( u(x) \).
$$ Z_N = T_N S_N \int_0^\infty dp \thinspace p^N \int_0^\infty dx_0 \int_0^{x_0} dx_1 ... dx_N \thinspace \delta \left( 1 - \frac{U(x_0)}{E} - \sum_i \frac{u(x_i)}{E} - p^2 \right) $$
Let's do the momentum integral as before,
$$ \begin{align*} Z_N &= T_N S_N \int_0^\infty dx_0 \int_0^{x_0} dx_1 ... dx_N \thinspace \int_0^\infty dp \thinspace p^N \delta \left( 1 - \frac{U(x_0)}{E} - \sum_i \frac{u(x_i)}{E} - p^2 \right) \\ &= T_N S_N \int_0^\infty dx_0 \int_0^{x_0} dx_1 ... dx_N \thinspace \int_0^\infty dp \thinspace p^N \delta \left(a^2 - p^2 \right) \\ &= \frac{1}{2} T_N S_N \int_0^\infty dx_0 \int_0^{x_0} dx_1 ... dx_N \thinspace \left( 1 - \frac{U}{E} - \sum_i \frac{u(x_i)}{E} \right)^\frac{N-1}{2} \theta\left( 1 - \frac{U}{E} - \sum_i \frac{u(x_i)}{E} \right) \\ \end{align*} $$
Let's do the same trick for large \( N \),
$$ \left( 1 - \frac{U(x_0)}{E} - \sum_i \frac{u(x_i)}{E} \right)^\frac{N-1}{2} \approx \exp \left( \frac{N}{2} \ln \left( 1 - \frac{U(x_0)}{E} - \sum_i \frac{u(x_i)}{E} \right) \right) $$
At this point we invoke Wikipedia. Yes, I've fooled you, I am not original, this is a trick from the 1700s. But it was hot then and it's hot now. Let's integrate over \( x_N \),
$$ \int_0^{x_0} dx_N \exp \left( \frac{N}{2} \ln \left( 1 - \frac{U(x_0)}{E} - \sum_i \frac{u(x_i)}{E} \right) \right) \approx \sqrt{\frac{2\pi}{N||}} $$
Other relations
$$ \frac{dS}{dE} = \frac{d}{dE} \ln(Z) = \frac{d}{dE} \ln \left(E^\frac{3N+1}{2} \right) = \frac{3N+1}{2E} = \beta $$
Question: how does time to thermalization depend on N?
Question:
Is it possible, for functions \( U(x) \) and \( v(y) \), to write
$$ \int_0^\infty dy \exp\left( - \frac{N}{2} \ln \left( 1 - U(x) - v(y) \right) \right) = \exp\left( - \frac{N}{2} \ln \left( 1 - U_\textrm{eff}(x)\right) \right) $$
and what is the value of the "effective" potential \( U_\textrm{eff}(x) \)?