[proofplan]
We first fix a compactly supported zero-mass perturbation $\sigma$ and choose a small interval of perturbation parameters on which $\rho+\varepsilon\sigma$ stays positive. The internal energy variation is localised to $\operatorname{supp}\sigma$, where smoothness of $U$ and compactness justify differentiating under the integral sign by dominated convergence. The potential energy is affine in the perturbation. For the interaction energy, we use the quadratic growth of $W$ and the finite second moment of $\rho$ to justify all double integrals, expand the expression algebraically, and use the evenness of $W$ to identify the two cross terms.
[/proofplan]
[step:Choose perturbation parameters preserving positivity and mass]
Let $\sigma\in C_c^\infty(\mathbb R^n;\mathbb R)$ satisfy
\begin{align*}
\int_{\mathbb R^n}\sigma(x)\,d\mathcal L^n(x)=0.
\end{align*}
Let $K=\operatorname{supp}\sigma$, which is a compact subset of $\mathbb R^n$. Since $\rho$ is continuous and strictly positive, the number
\begin{align*}
m=\min_{x\in K}\rho(x)
\end{align*}
is positive if $K\neq\varnothing$. If $K=\varnothing$, then $\sigma=0$ and all three formulas hold, so assume $K\neq\varnothing$. Let
\begin{align*}
M_\sigma=\sup_{x\in K}|\sigma(x)|.
\end{align*}
If $M_\sigma=0$, again $\sigma=0$. Otherwise set
\begin{align*}
\varepsilon_0=\frac{m}{2M_\sigma}.
\end{align*}
For every $\varepsilon\in(-\varepsilon_0,\varepsilon_0)$ and every $x\in K$,
\begin{align*}
\rho(x)+\varepsilon\sigma(x)\ge m-|\varepsilon|M_\sigma> \frac{m}{2}>0.
\end{align*}
For $x\notin K$, $\sigma(x)=0$, so $\rho(x)+\varepsilon\sigma(x)=\rho(x)>0$. Thus $\rho+\varepsilon\sigma$ is positive. Moreover,
\begin{align*}
\int_{\mathbb R^n}(\rho(x)+\varepsilon\sigma(x))\,d\mathcal L^n(x)=1+\varepsilon\int_{\mathbb R^n}\sigma(x)\,d\mathcal L^n(x)=1.
\end{align*}
Hence these perturbations are mass-preserving probability densities.
[/step]
[step:Differentiate the internal energy on the compact support of the perturbation]
Define the internal-energy perturbation map
\begin{align*}
I:(-\varepsilon_0,\varepsilon_0)\to\mathbb R,\qquad \varepsilon\mapsto \mathcal U[\rho+\varepsilon\sigma].
\end{align*}
Since $\sigma=0$ on $\mathbb R^n\setminus K$, we have
\begin{align*}
I(\varepsilon)-I(0)=\int_K\left(U(\rho(x)+\varepsilon\sigma(x))-U(\rho(x))\right)\,d\mathcal L^n(x).
\end{align*}
Let
\begin{align*}
R=\left[\frac{m}{2},\sup_{x\in K}\rho(x)+\varepsilon_0M_\sigma\right]\subset(0,\infty).
\end{align*}
The derivative $U'$ is continuous on the compact interval $R$, so
\begin{align*}
A=\sup_{r\in R}|U'(r)|<\infty.
\end{align*}
For $0<|\varepsilon|<\varepsilon_0$, the mean value theorem applied to the function $r\mapsto U(r)$ gives
\begin{align*}
\left|\frac{U(\rho(x)+\varepsilon\sigma(x))-U(\rho(x))}{\varepsilon}\right|\le A|\sigma(x)|
\end{align*}
for every $x\in K$. The function $x\mapsto A|\sigma(x)|$ belongs to $L^1(K,\mathcal L^n)$ because $\sigma$ is smooth with compact support. For each $x\in K$, ordinary one-variable differentiation gives
\begin{align*}
\lim_{\varepsilon\to0}\frac{U(\rho(x)+\varepsilon\sigma(x))-U(\rho(x))}{\varepsilon}=U'(\rho(x))\sigma(x).
\end{align*}
By dominated convergence,
\begin{align*}
\delta\mathcal U[\rho;\sigma]=I'(0)=\int_K U'(\rho(x))\sigma(x)\,d\mathcal L^n(x).
\end{align*}
Since $\sigma=0$ on $\mathbb R^n\setminus K$, this is
\begin{align*}
\delta\mathcal U[\rho;\sigma]=\int_{\mathbb R^n}U'(\rho(x))\sigma(x)\,d\mathcal L^n(x).
\end{align*}
[guided]
The only possible difficulty in the internal energy is that the integral is over all of $\mathbb R^n$, while the assumptions on $U(\rho)$ do not give a global domination for the differentiated integrand. The compact support of $\sigma$ removes that problem. Define
\begin{align*}
I:(-\varepsilon_0,\varepsilon_0)\to\mathbb R,\qquad \varepsilon\mapsto \mathcal U[\rho+\varepsilon\sigma].
\end{align*}
Because $\sigma(x)=0$ for $x\notin K$, the difference $I(\varepsilon)-I(0)$ is supported only on $K$:
\begin{align*}
I(\varepsilon)-I(0)=\int_K\left(U(\rho(x)+\varepsilon\sigma(x))-U(\rho(x))\right)\,d\mathcal L^n(x).
\end{align*}
Thus we only need domination on the compact set $K$.
The perturbation parameter $\varepsilon_0$ was chosen so that $\rho(x)+\varepsilon\sigma(x)\ge m/2$ on $K$. Also, for $x\in K$ and $|\varepsilon|<\varepsilon_0$,
\begin{align*}
\rho(x)+\varepsilon\sigma(x)\le \sup_{z\in K}\rho(z)+\varepsilon_0M_\sigma.
\end{align*}
Hence all values of $\rho(x)+\varepsilon\sigma(x)$ lie in the compact interval
\begin{align*}
R=\left[\frac{m}{2},\sup_{z\in K}\rho(z)+\varepsilon_0M_\sigma\right]\subset(0,\infty).
\end{align*}
Since $U\in C^\infty((0,\infty);\mathbb R)$, the derivative $U'$ is continuous, and therefore bounded on $R$. Define
\begin{align*}
A=\sup_{r\in R}|U'(r)|.
\end{align*}
For $0<|\varepsilon|<\varepsilon_0$, the one-variable mean value theorem applied between $\rho(x)$ and $\rho(x)+\varepsilon\sigma(x)$ gives
\begin{align*}
\left|\frac{U(\rho(x)+\varepsilon\sigma(x))-U(\rho(x))}{\varepsilon}\right|\le A|\sigma(x)|.
\end{align*}
The right-hand side is integrable over $K$ because $\sigma$ is smooth and compactly supported. Pointwise in $x$, the same one-variable differentiability gives
\begin{align*}
\lim_{\varepsilon\to0}\frac{U(\rho(x)+\varepsilon\sigma(x))-U(\rho(x))}{\varepsilon}=U'(\rho(x))\sigma(x).
\end{align*}
The dominated convergence theorem now applies to the difference quotients on the finite-measure compact set $K$, yielding
\begin{align*}
I'(0)=\int_K U'(\rho(x))\sigma(x)\,d\mathcal L^n(x).
\end{align*}
Finally, since $\sigma$ vanishes outside $K$, extending the domain of integration from $K$ to $\mathbb R^n$ does not change the integral:
\begin{align*}
\delta\mathcal U[\rho;\sigma]=\int_{\mathbb R^n}U'(\rho(x))\sigma(x)\,d\mathcal L^n(x).
\end{align*}
[/guided]
[/step]
[step:Use linearity to compute the potential energy variation]
Define the potential-energy perturbation map
\begin{align*}
J:(-\varepsilon_0,\varepsilon_0)\to\mathbb R,\qquad \varepsilon\mapsto \mathcal V[\rho+\varepsilon\sigma].
\end{align*}
Since $V$ has at most quadratic growth, there is a constant $C_V>0$ such that
\begin{align*}
|V(x)|\le C_V(1+|x|^2)
\end{align*}
for all $x\in\mathbb R^n$. The finite second moment of $\rho$ gives $V\rho\in L^1(\mathbb R^n,\mathcal L^n)$, and compact support of $\sigma$ gives $V\sigma\in L^1(\mathbb R^n,\mathcal L^n)$. Therefore
\begin{align*}
J(\varepsilon)=\int_{\mathbb R^n}V(x)\rho(x)\,d\mathcal L^n(x)+\varepsilon\int_{\mathbb R^n}V(x)\sigma(x)\,d\mathcal L^n(x).
\end{align*}
Differentiating this affine function of $\varepsilon$ at $0$ gives
\begin{align*}
\delta\mathcal V[\rho;\sigma]=\int_{\mathbb R^n}V(x)\sigma(x)\,d\mathcal L^n(x).
\end{align*}
[/step]
[step:Verify the interaction integrals are finite]
Since $W$ has at most quadratic growth, there is a constant $C_W>0$ such that
\begin{align*}
|W(z)|\le C_W(1+|z|^2)
\end{align*}
for every $z\in\mathbb R^n$. For $x,y\in\mathbb R^n$,
\begin{align*}
|x-y|^2\le 2|x|^2+2|y|^2.
\end{align*}
Hence
\begin{align*}
|W(x-y)|\le C_W(1+2|x|^2+2|y|^2).
\end{align*}
Using $\int_{\mathbb R^n}\rho\,d\mathcal L^n=1$ and the finite second moment of $\rho$, we obtain
\begin{align*}
\int_{\mathbb R^n}\int_{\mathbb R^n}|W(x-y)|\rho(x)\rho(y)\,d\mathcal L^n(y)\,d\mathcal L^n(x)<\infty.
\end{align*}
Thus $\mathcal W[\rho]$ is finite.
For each fixed $x\in\mathbb R^n$, the same estimate gives
\begin{align*}
\int_{\mathbb R^n}|W(x-y)|\rho(y)\,d\mathcal L^n(y)<\infty,
\end{align*}
so the convolution
\begin{align*}
W*\rho:\mathbb R^n\to\mathbb R,\qquad x\mapsto \int_{\mathbb R^n}W(x-y)\rho(y)\,d\mathcal L^n(y)
\end{align*}
is well-defined. Since $\sigma$ is bounded and supported in the compact set $K$, the estimates above also imply
\begin{align*}
\int_{\mathbb R^n}\int_{\mathbb R^n}|W(x-y)|\,|\sigma(x)|\rho(y)\,d\mathcal L^n(y)\,d\mathcal L^n(x)<\infty.
\end{align*}
The analogous integrals with $\rho(x)|\sigma(y)|$ and with $|\sigma(x)||\sigma(y)|$ are finite as well. Therefore all terms in the interaction-energy expansion below are absolutely integrable, and Fubini's theorem applies to them.
[/step]
[step:Expand the interaction energy and identify the two cross terms]
Define the interaction-energy perturbation map
\begin{align*}
K_W:(-\varepsilon_0,\varepsilon_0)\to\mathbb R,\qquad \varepsilon\mapsto \mathcal W[\rho+\varepsilon\sigma].
\end{align*}
By bilinearity of multiplication and the absolute integrability verified above,
\begin{align*}
K_W(\varepsilon)=\mathcal W[\rho]+\frac{\varepsilon}{2}A+\frac{\varepsilon}{2}B+\frac{\varepsilon^2}{2}C,
\end{align*}
where
\begin{align*}
A=\int_{\mathbb R^n}\int_{\mathbb R^n}W(x-y)\sigma(x)\rho(y)\,d\mathcal L^n(y)\,d\mathcal L^n(x),
\end{align*}
\begin{align*}
B=\int_{\mathbb R^n}\int_{\mathbb R^n}W(x-y)\rho(x)\sigma(y)\,d\mathcal L^n(y)\,d\mathcal L^n(x),
\end{align*}
and
\begin{align*}
C=\int_{\mathbb R^n}\int_{\mathbb R^n}W(x-y)\sigma(x)\sigma(y)\,d\mathcal L^n(y)\,d\mathcal L^n(x).
\end{align*}
Using Fubini's theorem to interchange the variables $x$ and $y$ in $B$, and using the evenness of $W$, we get
\begin{align*}
B=\int_{\mathbb R^n}\int_{\mathbb R^n}W(y-x)\sigma(x)\rho(y)\,d\mathcal L^n(y)\,d\mathcal L^n(x)=A.
\end{align*}
Therefore
\begin{align*}
K_W(\varepsilon)=\mathcal W[\rho]+\varepsilon A+\frac{\varepsilon^2}{2}C.
\end{align*}
Differentiating this quadratic polynomial in $\varepsilon$ at $0$ gives
\begin{align*}
\delta\mathcal W[\rho;\sigma]=A.
\end{align*}
By the definition of $W*\rho$ and Fubini's theorem,
\begin{align*}
A=\int_{\mathbb R^n}(W*\rho)(x)\sigma(x)\,d\mathcal L^n(x).
\end{align*}
Thus
\begin{align*}
\delta\mathcal W[\rho;\sigma]=\int_{\mathbb R^n}(W*\rho)(x)\sigma(x)\,d\mathcal L^n(x).
\end{align*}
[guided]
The interaction energy is quadratic in the density, so its first variation comes from the cross terms. The first point is to expand the perturbed product:
\begin{align*}
(\rho(x)+\varepsilon\sigma(x))(\rho(y)+\varepsilon\sigma(y))=\rho(x)\rho(y)+\varepsilon\sigma(x)\rho(y)+\varepsilon\rho(x)\sigma(y)+\varepsilon^2\sigma(x)\sigma(y).
\end{align*}
The preceding step verified absolute integrability for each term after multiplication by $W(x-y)$. Hence we may integrate this algebraic identity term by term. Define
\begin{align*}
A=\int_{\mathbb R^n}\int_{\mathbb R^n}W(x-y)\sigma(x)\rho(y)\,d\mathcal L^n(y)\,d\mathcal L^n(x),
\end{align*}
\begin{align*}
B=\int_{\mathbb R^n}\int_{\mathbb R^n}W(x-y)\rho(x)\sigma(y)\,d\mathcal L^n(y)\,d\mathcal L^n(x),
\end{align*}
and
\begin{align*}
C=\int_{\mathbb R^n}\int_{\mathbb R^n}W(x-y)\sigma(x)\sigma(y)\,d\mathcal L^n(y)\,d\mathcal L^n(x).
\end{align*}
Then
\begin{align*}
\mathcal W[\rho+\varepsilon\sigma]=\mathcal W[\rho]+\frac{\varepsilon}{2}A+\frac{\varepsilon}{2}B+\frac{\varepsilon^2}{2}C.
\end{align*}
Why do the two linear terms match? The term $B$ has $\sigma$ in the $y$ variable, whereas $A$ has $\sigma$ in the $x$ variable. Since $B$ is absolutely integrable, Fubini's theorem permits us to interchange the two variables:
\begin{align*}
B=\int_{\mathbb R^n}\int_{\mathbb R^n}W(y-x)\sigma(x)\rho(y)\,d\mathcal L^n(y)\,d\mathcal L^n(x).
\end{align*}
The function $W$ is even, so $W(y-x)=W(-(x-y))=W(x-y)$. Therefore
\begin{align*}
B=A.
\end{align*}
This is exactly where the factor $1/2$ in the definition of $\mathcal W$ is used: the two identical linear terms combine to one copy of $A$. Thus
\begin{align*}
\mathcal W[\rho+\varepsilon\sigma]=\mathcal W[\rho]+\varepsilon A+\frac{\varepsilon^2}{2}C.
\end{align*}
Differentiating this polynomial at $\varepsilon=0$ gives
\begin{align*}
\delta\mathcal W[\rho;\sigma]=A.
\end{align*}
Finally, the convolution was defined by
\begin{align*}
(W*\rho)(x)=\int_{\mathbb R^n}W(x-y)\rho(y)\,d\mathcal L^n(y).
\end{align*}
Substituting this definition into $A$ gives
\begin{align*}
A=\int_{\mathbb R^n}(W*\rho)(x)\sigma(x)\,d\mathcal L^n(x).
\end{align*}
Hence
\begin{align*}
\delta\mathcal W[\rho;\sigma]=\int_{\mathbb R^n}(W*\rho)(x)\sigma(x)\,d\mathcal L^n(x).
\end{align*}
[/guided]
[/step]
[step:Pass from directional formulas to representatives modulo constants]
The three computations show that, for every $\sigma\in C_c^\infty(\mathbb R^n;\mathbb R)$ with zero integral,
\begin{align*}
\delta\mathcal U[\rho;\sigma]=\int_{\mathbb R^n}U'(\rho(x))\sigma(x)\,d\mathcal L^n(x),
\end{align*}
\begin{align*}
\delta\mathcal V[\rho;\sigma]=\int_{\mathbb R^n}V(x)\sigma(x)\,d\mathcal L^n(x),
\end{align*}
and
\begin{align*}
\delta\mathcal W[\rho;\sigma]=\int_{\mathbb R^n}(W*\rho)(x)\sigma(x)\,d\mathcal L^n(x).
\end{align*}
If $q:\mathbb R^n\to\mathbb R$ is any one of the functions $U'(\rho)$, $V$, or $W*\rho$, and if $c\in\mathbb R$, then
\begin{align*}
\int_{\mathbb R^n}(q(x)+c)\sigma(x)\,d\mathcal L^n(x)=\int_{\mathbb R^n}q(x)\sigma(x)\,d\mathcal L^n(x)+c\int_{\mathbb R^n}\sigma(x)\,d\mathcal L^n(x)=\int_{\mathbb R^n}q(x)\sigma(x)\,d\mathcal L^n(x).
\end{align*}
Thus mass-preserving perturbations determine first-variation representatives only up to additive constants. Consequently,
\begin{align*}
\frac{\delta\mathcal U}{\delta\rho}=U'(\rho),\qquad \frac{\delta\mathcal V}{\delta\rho}=V,\qquad \frac{\delta\mathcal W}{\delta\rho}=W*\rho
\end{align*}
are valid representatives modulo additive constants. This completes the proof.
[/step]
