[proofplan]
We test the minimizer against compactly supported smooth perturbations with zero total mass, so that the probability constraint is preserved. For sufficiently small perturbation parameter, positivity and finite energy remain valid because the perturbation is supported on a compact subset of $\Omega$. Minimality then forces the first derivative of the energy along every such perturbation to vanish. Computing this first variation gives zero pairing against every compactly supported zero-mean test function, and a zero-mean fundamental lemma on connected domains implies that the continuous first variation is constant.
[/proofplan]
[step:Choose compactly supported mass-preserving perturbations]
Let $\sigma\in C_c^\infty(\Omega;\mathbb R)$ satisfy
\begin{align*}
\int_\Omega \sigma(x)\,d\mathcal L^n(x)=0.
\end{align*}
Let $K:=\operatorname{supp}\sigma$, which is a compact subset of $\Omega$. Since $\rho_\infty\in C^\infty(\Omega;(0,\infty))$, the continuous function $\rho_\infty|_K$ has a positive minimum
\begin{align*}
m_K:=\min_{x\in K}\rho_\infty(x)>0.
\end{align*}
Define
\begin{align*}
M_\sigma:=\sup_{x\in K}|\sigma(x)|.
\end{align*}
If $M_\sigma=0$, then $\sigma=0$ and the perturbation is harmless. If $M_\sigma>0$, choose
\begin{align*}
\varepsilon_0:=\frac{m_K}{2M_\sigma}.
\end{align*}
For every $\varepsilon\in(-\varepsilon_0,\varepsilon_0)$, define
\begin{align*}
\rho_\varepsilon:\Omega&\to(0,\infty)
\end{align*}
\begin{align*}
x&\mapsto \rho_\infty(x)+\varepsilon\sigma(x).
\end{align*}
Then $\rho_\varepsilon\in C^\infty(\Omega;(0,\infty))$. Indeed, on $K$ we have
\begin{align*}
\rho_\varepsilon(x)\ge \rho_\infty(x)-|\varepsilon|\,|\sigma(x)|\ge m_K-\varepsilon_0M_\sigma=\frac{m_K}{2}>0,
\end{align*}
and on $\Omega\setminus K$ we have $\rho_\varepsilon=\rho_\infty$. Moreover,
\begin{align*}
\int_\Omega \rho_\varepsilon(x)\,d\mathcal L^n(x)
=
\int_\Omega \rho_\infty(x)\,d\mathcal L^n(x)
+
\varepsilon\int_\Omega \sigma(x)\,d\mathcal L^n(x)
=
1.
\end{align*}
Thus $\rho_\varepsilon$ is a smooth positive probability density for all sufficiently small $\varepsilon$.
[/step]
[step:Verify that compact perturbations have finite energy]
We next check that $\mathcal E[\rho_\varepsilon]$ is finite for sufficiently small $\varepsilon$.
For the internal energy, $\rho_\varepsilon=\rho_\infty$ on $\Omega\setminus K$. Define
\begin{align*}
a_K:=\frac{m_K}{2}
\end{align*}
and
\begin{align*}
b_K:=\max_{x\in K}\rho_\infty(x)+\varepsilon_0M_\sigma.
\end{align*}
On $K$, the functions $\rho_\varepsilon$ take values in the compact interval $[a_K,b_K]\subset(0,\infty)$ for all $|\varepsilon|<\varepsilon_0$, because $\rho_\infty$ and $\sigma$ are continuous on $K$. Since $U\in C^1((0,\infty);\mathbb R)$, $U$ is bounded on $[a_K,b_K]$. As $\mathcal L^n(K)<\infty$, the difference
\begin{align*}
\int_K U(\rho_\varepsilon(x))\,d\mathcal L^n(x)
-
\int_K U(\rho_\infty(x))\,d\mathcal L^n(x)
\end{align*}
is finite. Since the internal energy of $\rho_\infty$ is finite as part of $\mathcal E[\rho_\infty]<\infty$, the internal energy of $\rho_\varepsilon$ is finite.
For the potential energy, $\rho_\varepsilon=\rho_\infty+\varepsilon\sigma$ gives
\begin{align*}
\int_\Omega V(x)\rho_\varepsilon(x)\,d\mathcal L^n(x)
=
\int_\Omega V(x)\rho_\infty(x)\,d\mathcal L^n(x)
+
\varepsilon\int_K V(x)\sigma(x)\,d\mathcal L^n(x).
\end{align*}
The first term is finite because $\mathcal E[\rho_\infty]<\infty$, and the second term is finite because $V$ and $\sigma$ are continuous on the compact set $K$.
For the interaction energy, expand the quadratic expression:
\begin{align*}
\frac{1}{2}\int_\Omega\int_\Omega W(x-y)\rho_\varepsilon(x)\rho_\varepsilon(y)\,d\mathcal L^n(y)\,d\mathcal L^n(x)
\end{align*}
equals the interaction energy of $\rho_\infty$, plus the two first-order cross terms, plus
\begin{align*}
\frac{\varepsilon^2}{2}\int_K\int_K W(x-y)\sigma(x)\sigma(y)\,d\mathcal L^n(y)\,d\mathcal L^n(x).
\end{align*}
The last integral is finite because $W$ is continuous on the compact set $K-K:=\{x-y:x,y\in K\}$ and $\sigma$ is bounded on $K$.
The stated local absolute-integrability hypothesis gives
\begin{align*}
\int_K\int_\Omega |W(x-y)|\rho_\infty(y)\,d\mathcal L^n(y)\,d\mathcal L^n(x)<\infty.
\end{align*}
Since $\sigma$ is bounded on $K$, Tonelli's theorem applies to the absolute value of the first cross term and shows that
\begin{align*}
\int_K\int_\Omega |W(x-y)|\,|\sigma(x)|\rho_\infty(y)\,d\mathcal L^n(y)\,d\mathcal L^n(x)<\infty.
\end{align*}
Thus Fubini's theorem identifies that cross term with
\begin{align*}
\frac{\varepsilon}{2}\int_K \sigma(x)\Phi_\infty(x)\,d\mathcal L^n(x).
\end{align*}
For the second cross term, use the evenness $W(x-y)=W(y-x)$ and the same absolute-integrability estimate with the compact variable in $K$ to obtain
\begin{align*}
\frac{\varepsilon}{2}\int_K \sigma(y)\Phi_\infty(y)\,d\mathcal L^n(y).
\end{align*}
Both integrals are finite because $\Phi_\infty$ and $\sigma$ are continuous on $K$. Therefore $\mathcal E[\rho_\varepsilon]<\infty$ for every sufficiently small $\varepsilon$.
[guided]
The only delicate point in using variations is admissibility: the theorem minimizes over smooth positive probability densities with finite energy, so each perturbation must remain inside that class.
We start with a test perturbation $\sigma\in C_c^\infty(\Omega;\mathbb R)$ satisfying
\begin{align*}
\int_\Omega \sigma(x)\,d\mathcal L^n(x)=0.
\end{align*}
Its compact support $K:=\operatorname{supp}\sigma$ is contained in $\Omega$. Since $\rho_\infty$ is smooth and strictly positive, it has a positive minimum on $K$:
\begin{align*}
m_K:=\min_{x\in K}\rho_\infty(x)>0.
\end{align*}
For sufficiently small $|\varepsilon|$, the function
\begin{align*}
\rho_\varepsilon:\Omega&\to(0,\infty)
\end{align*}
\begin{align*}
x&\mapsto \rho_\infty(x)+\varepsilon\sigma(x)
\end{align*}
remains positive. Indeed, if $M_\sigma:=\sup_{x\in K}|\sigma(x)|$ and $M_\sigma>0$, then for $|\varepsilon|<m_K/(2M_\sigma)$,
\begin{align*}
\rho_\varepsilon(x)\ge m_K-|\varepsilon|M_\sigma>\frac{m_K}{2}
\end{align*}
for every $x\in K$, while outside $K$ the perturbation vanishes and $\rho_\varepsilon=\rho_\infty$. The zero-mean condition preserves the mass:
\begin{align*}
\int_\Omega \rho_\varepsilon(x)\,d\mathcal L^n(x)
=
\int_\Omega \rho_\infty(x)\,d\mathcal L^n(x)
+
\varepsilon\int_\Omega \sigma(x)\,d\mathcal L^n(x)
=
1.
\end{align*}
It remains to check finite energy. This is where compact support matters. The internal energy changes only on $K$. Define
\begin{align*}
a_K:=\frac{m_K}{2}
\end{align*}
and
\begin{align*}
b_K:=\max_{x\in K}\rho_\infty(x)+\varepsilon_0M_\sigma.
\end{align*}
On $K$, the values of $\rho_\varepsilon$ remain inside the compact interval $[a_K,b_K]\subset(0,\infty)$, and $U$ is bounded on that interval because $U\in C^1((0,\infty);\mathbb R)$. Therefore
\begin{align*}
\int_K U(\rho_\varepsilon(x))\,d\mathcal L^n(x)
\end{align*}
is finite, and outside $K$ the internal energy is unchanged.
For the potential term, linearity gives
\begin{align*}
\int_\Omega V(x)\rho_\varepsilon(x)\,d\mathcal L^n(x)
=
\int_\Omega V(x)\rho_\infty(x)\,d\mathcal L^n(x)
+
\varepsilon\int_K V(x)\sigma(x)\,d\mathcal L^n(x).
\end{align*}
The first integral is finite by the finite-energy hypothesis on $\rho_\infty$, and the second is finite because $V$ and $\sigma$ are continuous on the compact set $K$.
For the interaction term, expanding the product gives
\begin{align*}
\rho_\varepsilon(x)\rho_\varepsilon(y)
=
\rho_\infty(x)\rho_\infty(y)
+
\varepsilon\sigma(x)\rho_\infty(y)
+
\varepsilon\rho_\infty(x)\sigma(y)
+
\varepsilon^2\sigma(x)\sigma(y).
\end{align*}
The zeroth-order part is finite by the finite-energy hypothesis. The quadratic perturbation part is supported on $K\times K$, and
\begin{align*}
\int_K\int_K W(x-y)\sigma(x)\sigma(y)\,d\mathcal L^n(y)\,d\mathcal L^n(x)
\end{align*}
is finite because $W$ is continuous on the compact difference set $K-K$ and $\sigma$ is bounded.
For the cross terms, the local absolute-integrability hypothesis gives
\begin{align*}
\int_K\int_\Omega |W(x-y)|\rho_\infty(y)\,d\mathcal L^n(y)\,d\mathcal L^n(x)<\infty.
\end{align*}
Because $\sigma$ is bounded on $K$, Tonelli's theorem gives absolute integrability of
\begin{align*}
(x,y)\mapsto W(x-y)\sigma(x)\rho_\infty(y)
\end{align*}
on $K\times\Omega$. Fubini's theorem then permits the order of integration used in the definition
\begin{align*}
\Phi_\infty(x)=\int_\Omega W(x-y)\rho_\infty(y)\,d\mathcal L^n(y)
\end{align*}
and gives
\begin{align*}
\int_K \sigma(x)\Phi_\infty(x)\,d\mathcal L^n(x).
\end{align*}
The other cross term is handled in the same way after using the evenness $W(x-y)=W(y-x)$ to place the compactly supported variable in the first slot of the same estimate. Since $\Phi_\infty$ and $\sigma$ are continuous on the compact set $K$, the resulting compact integral is finite. Hence every sufficiently small perturbation $\rho_\varepsilon$ is an admissible competitor.
[/guided]
[/step]
[step:Use minimality to force the first variation to vanish]
For the fixed test function $\sigma$, define
\begin{align*}
e_\sigma:(-\varepsilon_0,\varepsilon_0)&\to\mathbb R
\end{align*}
\begin{align*}
\varepsilon&\mapsto \mathcal E[\rho_\varepsilon].
\end{align*}
By the previous step, $e_\sigma$ is finite on a neighbourhood of $0$. Since $\rho_\infty$ minimizes $\mathcal E$ among admissible densities and each $\rho_\varepsilon$ is admissible for sufficiently small $\varepsilon$, the real-valued function $e_\sigma$ has a minimum at $\varepsilon=0$. Therefore, once we compute the derivative at $0$, it must satisfy
\begin{align*}
e_\sigma'(0)=0.
\end{align*}
[/step]
[step:Compute the first variation of the energy]
Define the first-variation function
\begin{align*}
F:\Omega&\to\mathbb R
\end{align*}
\begin{align*}
x&\mapsto U'(\rho_\infty(x))+V(x)+\Phi_\infty(x).
\end{align*}
The function $F$ is continuous because $U'\in C((0,\infty);\mathbb R)$, $\rho_\infty\in C^\infty(\Omega;(0,\infty))$, $V\in C(\Omega;\mathbb R)$, and $\Phi_\infty\in C(\Omega;\mathbb R)$.
For the internal energy, the difference quotient is supported in $K$:
\begin{align*}
\frac{1}{\varepsilon}
\left[
\int_\Omega U(\rho_\varepsilon(x))\,d\mathcal L^n(x)
-
\int_\Omega U(\rho_\infty(x))\,d\mathcal L^n(x)
\right]
=
\int_K
\frac{U(\rho_\infty(x)+\varepsilon\sigma(x))-U(\rho_\infty(x))}{\varepsilon}
\,d\mathcal L^n(x).
\end{align*}
Define
\begin{align*}
a_K:=\frac{m_K}{2}
\end{align*}
and
\begin{align*}
b_K:=\max_{x\in K}\rho_\infty(x)+\varepsilon_0M_\sigma.
\end{align*}
For $|\varepsilon|<\varepsilon_0$, the arguments of $U$ lie in the compact interval $[a_K,b_K]\subset(0,\infty)$. By the mean value theorem and boundedness of $U'$ on $[a_K,b_K]$, the integrand is dominated by
\begin{align*}
\left(\sup_{r\in[a_K,b_K]}|U'(r)|\right)|\sigma(x)|.
\end{align*}
This dominating function is integrable on $K$ with respect to $\mathcal L^n$. Passing to the limit by dominated convergence gives
\begin{align*}
\frac{d}{d\varepsilon}\bigg|_{\varepsilon=0}
\int_\Omega U(\rho_\varepsilon(x))\,d\mathcal L^n(x)
=
\int_K U'(\rho_\infty(x))\sigma(x)\,d\mathcal L^n(x).
\end{align*}
For the potential energy, linearity gives
\begin{align*}
\frac{d}{d\varepsilon}\bigg|_{\varepsilon=0}
\int_\Omega V(x)\rho_\varepsilon(x)\,d\mathcal L^n(x)
=
\int_K V(x)\sigma(x)\,d\mathcal L^n(x).
\end{align*}
For the interaction energy, the expansion from the previous step is justified by absolute convergence of the original interaction integral, the compact $K\times K$ quadratic term, and the local absolute-integrability estimate for the two cross terms. The quadratic perturbation term is multiplied by $\varepsilon^2$, so its derivative at $0$ is zero. The two first-order cross terms are identified by Fubini's theorem and the evenness of $W$, as above, and their sum is
\begin{align*}
\int_K \Phi_\infty(x)\sigma(x)\,d\mathcal L^n(x).
\end{align*}
Therefore
\begin{align*}
\frac{d}{d\varepsilon}\bigg|_{\varepsilon=0}
\frac{1}{2}
\int_\Omega\int_\Omega W(x-y)\rho_\varepsilon(x)\rho_\varepsilon(y)\,d\mathcal L^n(y)\,d\mathcal L^n(x)
=
\int_K \Phi_\infty(x)\sigma(x)\,d\mathcal L^n(x).
\end{align*}
Combining the three derivative computations yields
\begin{align*}
e_\sigma'(0)
=
\int_K
\left[
U'(\rho_\infty(x))+V(x)+\Phi_\infty(x)
\right]\sigma(x)\,d\mathcal L^n(x).
\end{align*}
Since $\sigma=0$ on $\Omega\setminus K$, this is equivalently
\begin{align*}
e_\sigma'(0)
=
\int_\Omega F(x)\sigma(x)\,d\mathcal L^n(x).
\end{align*}
The minimality condition from the previous step therefore gives
\begin{align*}
\int_\Omega F(x)\sigma(x)\,d\mathcal L^n(x)=0
\end{align*}
for every $\sigma\in C_c^\infty(\Omega;\mathbb R)$ satisfying
\begin{align*}
\int_\Omega \sigma(x)\,d\mathcal L^n(x)=0.
\end{align*}
[/step]
[step:Prove the zero-mean fundamental lemma on a connected domain]
[claim:Continuous functions annihilating zero-mean tests are constant on connected domains]
Let $\Omega\subset\mathbb R^n$ be connected and open. If $F\in C(\Omega;\mathbb R)$ satisfies
\begin{align*}
\int_\Omega F(x)\sigma(x)\,d\mathcal L^n(x)=0
\end{align*}
for every $\sigma\in C_c^\infty(\Omega;\mathbb R)$ with
\begin{align*}
\int_\Omega \sigma(x)\,d\mathcal L^n(x)=0,
\end{align*}
then there exists $c\in\mathbb R$ such that $F(x)=c$ for every $x\in\Omega$.
[/claim]
[proof]
Choose a nonnegative function $\eta\in C_c^\infty(\Omega;\mathbb R)$ such that
\begin{align*}
\int_\Omega \eta(x)\,d\mathcal L^n(x)=1.
\end{align*}
Such a function exists because $\Omega$ is nonempty and open. Define
\begin{align*}
c:=\int_\Omega F(x)\eta(x)\,d\mathcal L^n(x).
\end{align*}
For any $\varphi\in C_c^\infty(\Omega;\mathbb R)$, define
\begin{align*}
\sigma_\varphi:\Omega&\to\mathbb R
\end{align*}
\begin{align*}
x&\mapsto \varphi(x)-\left(\int_\Omega \varphi(y)\,d\mathcal L^n(y)\right)\eta(x).
\end{align*}
Then $\sigma_\varphi\in C_c^\infty(\Omega;\mathbb R)$ and
\begin{align*}
\int_\Omega \sigma_\varphi(x)\,d\mathcal L^n(x)
=
\int_\Omega \varphi(x)\,d\mathcal L^n(x)
-
\left(\int_\Omega \varphi(y)\,d\mathcal L^n(y)\right)
\int_\Omega \eta(x)\,d\mathcal L^n(x)
=
0.
\end{align*}
By hypothesis,
\begin{align*}
0
=
\int_\Omega F(x)\sigma_\varphi(x)\,d\mathcal L^n(x)
=
\int_\Omega F(x)\varphi(x)\,d\mathcal L^n(x)
-
c\int_\Omega \varphi(x)\,d\mathcal L^n(x).
\end{align*}
Hence
\begin{align*}
\int_\Omega (F(x)-c)\varphi(x)\,d\mathcal L^n(x)=0
\end{align*}
for every $\varphi\in C_c^\infty(\Omega;\mathbb R)$.
We now show that this forces $F-c=0$ pointwise. Suppose there exists $x_0\in\Omega$ with $F(x_0)>c$. By continuity, there exist $r>0$ and $\alpha>0$ such that $B(x_0,r)\subset\Omega$ and
\begin{align*}
F(x)-c\ge \alpha
\end{align*}
for every $x\in B(x_0,r)$. Choose a nonnegative function $\varphi\in C_c^\infty(B(x_0,r);\mathbb R)$ with
\begin{align*}
\int_{B(x_0,r)}\varphi(x)\,d\mathcal L^n(x)>0.
\end{align*}
Then
\begin{align*}
0
=
\int_\Omega (F(x)-c)\varphi(x)\,d\mathcal L^n(x)
=
\int_{B(x_0,r)} (F(x)-c)\varphi(x)\,d\mathcal L^n(x)
\ge
\alpha\int_{B(x_0,r)}\varphi(x)\,d\mathcal L^n(x)>0,
\end{align*}
a contradiction. The case $F(x_0)<c$ is identical after replacing $\varphi$ by a nonnegative test function supported where $c-F$ is positive. Therefore $F(x)=c$ for every $x\in\Omega$.
[/proof]
Applying the claim to the continuous function
\begin{align*}
F(x)=U'(\rho_\infty(x))+V(x)+\Phi_\infty(x)
\end{align*}
gives a constant $c\in\mathbb R$ such that
\begin{align*}
U'(\rho_\infty(x))+V(x)+\Phi_\infty(x)=c
\end{align*}
for every $x\in\Omega$. This is the desired Euler-Lagrange condition.
[/step]
