Is number density uniform ($n = N/V$) in equilibrium for Boltzmann distribution? It is often said in the physics papers or textbooks that in equilibrium if there is no external potential then number density is uniform $n = N/V$, where $N$ is number of particles, $V$ is volume. I would like to see if this follows from statistical mechanics (or maybe contradicts statistical mechanics).
Let us assume that we are given $N$ particles under Hamiltonian that has only pairwise interaction.
$$ H = \sum \limits_{i=1}^{N} \frac{\vec{p}_i^2}{2m} + \sum \limits_{i=1}^{N} \sum \limits_{j=i+1}^N \Phi(|\vec{x}_i - \vec{x}_j|)$$
Assume that distribution in the phase space is given by Boltzmann distribution.
$$ f(\vec{x}_1, ..., \vec{x}_N, \vec{p}_1, ..., \vec{p}_N) = \begin{cases} \frac{1}{Z} e^{-\beta H(\vec{x}_1, ..., \vec{x}_N, \vec{p}_1, ..., \vec{p}_N)} \quad \textrm{ if all } \vec{x}_i \in V  \\ 0 
 \textrm{,} \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \textrm{else.}  \end{cases} $$
Define number density as the following quantity.
$$ n(\vec{x}) = N \int \delta(\vec{x}-\vec{x}_1) f(\vec{x}_1, ..., \vec{x}_N, \vec{p}_1, ..., \vec{p}_N) 
\prod \limits_{i=1}^N \mathrm{d}\vec{x}_i \mathrm{d}\vec{p}_i$$
Because of properties of Boltzmann distribution, it is possible to integrate out all momentum coordinates. In that case, we have the following simplification, where $Q_N$ is defined in the following equation.
$$ n(\vec{x}) = \frac{N}{Q_N} \int \limits_{V^N} \delta(\vec{x} - \vec{x}_1) e^{-\beta \sum \limits_{i=1}^{N} \sum \limits_{j=i+1}^N \Phi(|\vec{x}_i - \vec{x}_j|)} \prod \limits_{i=1}^N \mathrm{d}\vec{x}_i$$
$$ Q_N = \int \limits_{V^N} e^{-\beta \sum \limits_{i=1}^{N} \sum \limits_{j=i+1}^N \Phi(|\vec{x}_i - \vec{x}_j|)} \prod \limits_{i=1}^N \mathrm{d}\vec{x}_i$$
$\textbf{Question}$: Is it possible using these definitions and assumptions to show that in this Boltzmann equilibrium model number density is uniform (or maybe that it is false that it is uniform)?
 A: (I am assuming that there is no external potential, apart from the box confining the particles; see Alexander's answer for a discussion of the case when this is not true.)
The density is uniform as long as the free energy is strictly convex in the density. Namely assume that the overall density $\rho = N/V$ belongs to an interval $I$ such that
$$
f_\beta(\lambda \rho_1+(1-\lambda)\rho_2)< \lambda
f_\beta(\rho_1)+(1-\lambda) f_\beta(\rho_2)
$$
for all $0<\lambda<1$ and all $\rho_1<\rho_2$ in $I$. Then the density is uniform over the whole system.

Although providing a full proof of the above statement would be too much for this site (a simple version of this claim is proved (for the lattice gas) in Theorem 4.19 in this book if you are interested), it is easy to explain heuristically why strict convexity should be relevant.
Consider a system in a (nice) box of volume $V$ with density $\rho=N/V$. Consider a (nice) sub-volume $\Delta$ of $V$ with volume $\alpha V$ for some $0 < \alpha < 1$.
Denoting $Q_{W;\beta,\rho'}$ the canonical partition function in a (nice) box of volume $W$ with density $\rho'$ (at inverse temperature $\beta$), we have
$$
Q_{W;\beta,\rho'} = e^{-\beta f_\beta(\rho') W + O(\partial W)},
$$
where the term $O(\partial W)$ takes into account boundary effects and is irrelevant in the thermodynamic limit.
Now assume that the density inside the sub-volume $\Delta$ is given by $\rho_1 < \rho$. Since the overall density in $V$ must be $\rho$, we conclude that the density inside the sub-box $V\setminus \Delta$ has to be given by $\rho_2>\rho$, where $\rho_2$ satisfies
$$
\rho_1 \alpha V + \rho_2 (1-\alpha) V = \rho V .
\tag{$\star$}
$$
On the one hand, as pointed out above, the full partition function is (up to negligible boundary terms)
$$
Q_{V;\beta,\rho} = e^{-\beta f_\beta(\rho) V}.
$$
On the other hand, the contribution to the partition function coming from configurations with density $\rho_1$ in $\Delta$ (and thus $\rho_2$ in $V\setminus\Delta$) is given (up to negligible boundary terms) by
$$
Q_{ \Delta;\beta,\rho_1}Q_{V\setminus\Delta;\beta,\rho_2}
= e^{-\beta f_\beta(\rho_1) \alpha V - \beta f_\beta(\rho_2) (1-\alpha) V}
= e^{-\beta \bigl( \alpha f_\beta(\rho_1) + (1-\alpha) f_\beta(\rho_2) \bigr) V}.
$$
If the free energy is strictly convex on the interval $(\rho_1,\rho_2)$, we have
$$
\alpha f_\beta(\rho_1) + (1-\alpha) f_\beta(\rho_2)
>
f_\beta\bigl( \alpha\rho_1 + (1-\alpha)\rho_2 \bigr)
=
f_\beta(\rho),
$$
where we have used ($\star$). This shows that
$$
Q_{ \Delta;\beta,\rho_1}Q_{V\setminus\Delta;\beta,\rho_2}
\leq e^{-c V} Q_{V;\beta,\rho}
$$
for some constant $c>0$. This shows that the probability of observing a density $\rho_1<\rho$ in $\Delta$ is smaller than $e^{-cV}$. In particular, in the thermodynamic limit, typical configurations must be spatially homogeneous.

Note that this is not true anymore when strict convexity fails. This is the case when $T<T_c$ and the density $\rho$ belongs to the interval $(\rho_g,\rho_l)$ where $\rho_g$ and $\rho_l$ are the densities of the gas and liquid phases respectively (although convex, the free energy is affine between $\rho_g$ and $\rho_l$ and thus not strictly convex). Then there will be phase separation: part of the system will be in the liquid phase, the rest in the gas phase. Spatial homegeneity is then lost.
This is discussed informally (with relevant bibliographical references) in Section 4.12.1 of the same book. Here's a picture of a typical configuration for the lattice gas (taken from the book):

A: In equilibrium state of non-interacting particles, number density is uniform only for box potential with some boundary conditions, e.g.
$$V(x)=\cases{0,x\in V\\ \infty,\text{othersize}}$$
For general potential the distribution is not uniform. Consider non-interacting particles in harmonic potential, that distribute as Gaussian. This is exactly $n(x) \sim \exp(-\beta V(x))$ distribution.
Pairwise interactions modify the $n(x) \sim \exp(-\beta V(x))$ distribution, however as a general case it is not uniform. In your equation is means the $n(x)$ is independent of $x$, which requires some engineering of the interaction.
