Partition function for weakly interacting gases I'm studying the Susskind lecture on statistical mechanics. A potential energy of pairwise interactions has been defined:
$$\sum_{n>m}U(|x_n -x_m|)$$
We want to calculate the partition function, and we start by evaluating the potential energy term by considering the interaction of a single pair, through the integral of potential energy between two particles over all possible positions:
$$\int dx_1^3 dx_2^3 U(|x_1-x_2|)$$
Now,there are 2 steps,

*

*Assuming that $U$ goes to zero at large distances, hold particle $1$ fixed and integrate over all possible positions of particle $2$:
$\int dx^3 U(|x|)=U_0$

*Integrating over the position of particle $1$ gives a factor of volume, $V$.

Conclusion:
$$\int dx_1^3 dx_2^3 U(|x_1-x_2|)=VU_0$$
I cannot understand steps 1 and 2. It looks pure magic to me. Is there a formal way to explain how the integral is broken up in two parts, why we can consider the position fixed, and why the results are equivalent?
 A: For simplicity, let us first assume that the potential $U$ has finite range: $U(|x|)=0$ whenever $|x|>R$. Then, for any $x_1$ located at a distance at least $R$ from the boundary of the vessel enclosing your system,
$$
\int_V\mathrm{d}x_2\, U(|x_1-x_2|) = \int_{\mathbb{R}^3}\mathrm{d}y\, U(|y|) = U_0,
$$
where the first identity follows from the change of variables $y=x_2-x_1$, (which has Jacobian $1$). Then, splitting the integral over $x_1$ into two parts according to whether $x_1$ is at a distance at least $R$ from the boundary or not, we obtain
$$
\int_V\mathrm{d}x_1\int_V\mathrm{d}x_2\, U(|x_1-x_2|) = U_0 V + O(R\partial V),
$$
where the error term is of the order of the area $\partial V$ of the vessel.
The situation is entirely similar if $U$ does not have finite range (but is, of course, integrable). First, we can fix some $\epsilon>0$ and observe that there exists $R(\epsilon)$ such that
$$
\Bigl| \int_{|y|\leq R(\epsilon)}\mathrm{d}y\, U(|y|) - U_0 \Bigr| \leq \epsilon.
$$
From this, we get
$$
\int_V\mathrm{d}x_1\int_V\mathrm{d}x_2\, U(|x_1-x_2|) = U_0 V + O(\epsilon V) + O(R(\epsilon)\partial V).
$$
Letting $V\to\infty$ (say, along a sequence of cubes or other reasonable shapes) and then letting $\epsilon\to 0$, we obtain
$$
\lim_{V\to\infty} \frac1V \int_V\mathrm{d}x_1\int_V
\mathrm{d}x_2\, U(|x_1-x_2|) = U_0,
$$
which is a precise way of stating the identity.
