I am trying to solve the following problem from an old qualifying exam:

"Ideal gas in gravitational potential"

Consider an ideal gas of N indistinguishable molecules of mass m in a cylindrical volume $V=Ah$ with base area A and height h.

Calculate the canonical partition function of the ideal gas including the effect of gravity. (Hint: You may find useful the integral $\int_0^\infty t^2 e^{-t^2}dt=\sqrt{\pi/4}$).

My work so far:

Since the partition function of a total system is the product of the partition function of the subsystems, i.e. if there are N subsystems, we'd have

$$ Z_{total} = Z_1 Z_2 Z_3 ... Z_N = \prod_i^N Z_i $$

Moreover, if (as is the case in this instance) the subsystems are indistinguishable, we can (after correct Boltzmann counting) reduce this to,

$$ Z_{total} = \frac{1}{N!} (Z_1)^N $$

where $Z_1$ is the partition function of 1 molecule of an ideal gas subject to gravity.

Question #1:

The solution accompanying this problem makes the following adjustment to my above expression for $Z_{total}$. Notably, they use Stirling's formula to write out N!. However, I'm either making a silly mistake, or the solution is wrong.

Using Stirling's formula

$$ N! \propto (\frac{N}{e})^N $$

I would think that

$$ Z_{total} = \frac{1}{N!} (Z_1)^N = \bigg(\frac{Z_1 e}{N}\bigg)^N $$

but the solution says

$$ Z_{total} = \bigg(\frac{Z_1 N}{e}\bigg)^N $$

Question #2:

$Z_1$, the partition function for one of the gas molecules subject to gravity is

$$ Z_1 = \frac{1}{h^3} \int d^3p d^3q e^{-\beta H} $$

where $h$ is Plank's constant, $\beta = 1/T$, and $H$ is the Hamiltionian for a single molecule accounting for both the particles momentum and gravity ($H=\frac{p^2}{2m} + mgy$, with y being the height of the particle).

The first step that the solution takes in evaluating this integral is the following

$$ \begin{align} Z_1 &= \frac{1}{h^3} \int d^3p d^3q e^{-\beta H} \\ &= \frac{1}{(2\pi \hbar)^3} \int e^{-\beta \frac{p^2}{2m}} d^3p \int e^{-\beta m g y} d^3q\\ &=\frac{1}{(2\pi \hbar)^3} \int e^{-\beta \frac{p^2}{2m}} d^3p \cdot A \int_0^h e^{-\beta m g y}dy \end{align} $$

where recall that A is the area of the cylinder and h in the height.

Now the next step is where I'm thoroughly confused,

$$ \begin{align} Z_1 &=\frac{1}{(2\pi \hbar)^3} \int e^{-\beta \frac{p^2}{2m}} d^3p \cdot A \int_0^h e^{-\beta m g y} dy \\ &=\frac{4 \pi}{(2\pi \hbar)^3} \int_0^\infty p^2 e^{-\beta \frac{p^2}{2m}} dp \cdot A \int_0^h e^{-\beta m g y} dy \end{align} $$

My only guess is switching to polar and changing the limits of integration?

Can someone clarify these two questions for me?

Cheers