[proofplan]
We compute the first variation of each term in $\mathcal E$ at a smooth positive density $\rho$ along a smooth mass-preserving perturbation. The internal and potential terms give $U'(\rho)$ and $V$, while the evenness of $W$ allows the interaction variation to symmetrise to $W*\rho$. Thus the formal variational derivative is $U'(\rho)+V+W*\rho$. The formal Wasserstein gradient-flow rule then sets the velocity equal to the negative spatial gradient of this variational derivative, and substituting this velocity into the continuity equation gives the displayed PDE.
[/proofplan]
[step:Compute the first variation along a mass-preserving perturbation]
Fix a time $t\in(0,T)$ and abbreviate $\rho=\rho_t$. Let
\begin{align*}
\sigma:\mathbb R^n\to\mathbb R
\end{align*}
be a smooth admissible perturbation with sufficient decay and zero total mass:
\begin{align*}
\int_{\mathbb R^n}\sigma(x)\,d\mathcal L^n(x)=0.
\end{align*}
For $\varepsilon$ in an interval around $0$, define
\begin{align*}
\rho_\varepsilon:\mathbb R^n&\to(0,\infty)
\end{align*}
by $\rho_\varepsilon(x)=\rho(x)+\varepsilon\sigma(x)$, assuming $\rho_\varepsilon$ remains positive and admissible for the formal computation. The zero-mass condition gives
\begin{align*}
\int_{\mathbb R^n}\rho_\varepsilon(x)\,d\mathcal L^n(x)=1,
\end{align*}
so $\rho_\varepsilon$ is a mass-preserving variation.
By the assumed justification for differentiating under the integral sign and by the chain rule applied to $U\in C^2((0,\infty);\mathbb R)$,
\begin{align*}
\left.\frac{d}{d\varepsilon}\right|_{\varepsilon=0}\int_{\mathbb R^n}U(\rho_\varepsilon(x))\,d\mathcal L^n(x)=\int_{\mathbb R^n}U'(\rho(x))\sigma(x)\,d\mathcal L^n(x).
\end{align*}
Similarly, differentiating the potential term gives
\begin{align*}
\left.\frac{d}{d\varepsilon}\right|_{\varepsilon=0}\int_{\mathbb R^n}V(x)\rho_\varepsilon(x)\,d\mathcal L^n(x)=\int_{\mathbb R^n}V(x)\sigma(x)\,d\mathcal L^n(x).
\end{align*}
[guided]
Fix $t\in(0,T)$ and write $\rho=\rho_t$ to focus on the variational calculation at one density. A Wasserstein tangent vector preserves total mass, so we test the first variation against a smooth admissible perturbation
\begin{align*}
\sigma:\mathbb R^n\to\mathbb R
\end{align*}
satisfying
\begin{align*}
\int_{\mathbb R^n}\sigma(x)\,d\mathcal L^n(x)=0.
\end{align*}
For $\varepsilon$ sufficiently close to $0$, define
\begin{align*}
\rho_\varepsilon:\mathbb R^n&\to(0,\infty)
\end{align*}
by $\rho_\varepsilon(x)=\rho(x)+\varepsilon\sigma(x)$, with admissibility and positivity included in the formal hypotheses. Then
\begin{align*}
\int_{\mathbb R^n}\rho_\varepsilon(x)\,d\mathcal L^n(x)=\int_{\mathbb R^n}\rho(x)\,d\mathcal L^n(x)+\varepsilon\int_{\mathbb R^n}\sigma(x)\,d\mathcal L^n(x)=1.
\end{align*}
Thus the variation stays in the affine hyperplane of probability densities.
For the internal energy, the map $\varepsilon\mapsto U(\rho(x)+\varepsilon\sigma(x))$ is differentiable for each $x$ because $U\in C^2((0,\infty);\mathbb R)$ and $\rho_\varepsilon(x)>0$. The theorem assumes enough decay and domination to move the $\varepsilon$-derivative through the integral. Therefore the chain rule gives
\begin{align*}
\left.\frac{d}{d\varepsilon}\right|_{\varepsilon=0}\int_{\mathbb R^n}U(\rho_\varepsilon(x))\,d\mathcal L^n(x)=\int_{\mathbb R^n}U'(\rho(x))\sigma(x)\,d\mathcal L^n(x).
\end{align*}
For the potential term, the integrand depends linearly on $\rho_\varepsilon$. Differentiation under the integral sign gives
\begin{align*}
\left.\frac{d}{d\varepsilon}\right|_{\varepsilon=0}\int_{\mathbb R^n}V(x)\rho_\varepsilon(x)\,d\mathcal L^n(x)=\int_{\mathbb R^n}V(x)\sigma(x)\,d\mathcal L^n(x).
\end{align*}
These two computations identify the internal and confinement contributions to the variational derivative as $U'(\rho)$ and $V$.
[/guided]
[/step]
[step:Use evenness of $W$ to identify the interaction variation]
Define the interaction functional $\mathcal W$ on admissible densities by
\begin{align*}
\mathcal W[\eta]=\frac12\int_{\mathbb R^n}\int_{\mathbb R^n}W(x-y)\eta(x)\eta(y)\,d\mathcal L^n(y)\,d\mathcal L^n(x).
\end{align*}
Using the assumed differentiability under the double integral,
\begin{align*}
\left.\frac{d}{d\varepsilon}\right|_{\varepsilon=0}\mathcal W[\rho_\varepsilon]=\frac12\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)+\frac12\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*}
In the second double integral, interchange the names of the variables $x$ and $y$. Since $\mathcal L^n\otimes\mathcal L^n$ is invariant under this coordinate swap and $W$ is even,
\begin{align*}
W(y-x)=W(-(x-y))=W(x-y).
\end{align*}
Hence the second double integral equals the first. Therefore
\begin{align*}
\left.\frac{d}{d\varepsilon}\right|_{\varepsilon=0}\mathcal W[\rho_\varepsilon]=\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*}
By the definition of convolution,
\begin{align*}
\left.\frac{d}{d\varepsilon}\right|_{\varepsilon=0}\mathcal W[\rho_\varepsilon]=\int_{\mathbb R^n}(W*\rho)(x)\sigma(x)\,d\mathcal L^n(x).
\end{align*}
[/step]
[step:Identify the formal variational derivative of $\mathcal E$]
Combining the three first-variation formulas, we obtain
\begin{align*}
\left.\frac{d}{d\varepsilon}\right|_{\varepsilon=0}\mathcal E[\rho_\varepsilon]=\int_{\mathbb R^n}\left(U'(\rho(x))+V(x)+(W*\rho)(x)\right)\sigma(x)\,d\mathcal L^n(x).
\end{align*}
Thus the formal variational derivative of $\mathcal E$ at $\rho$ is
\begin{align*}
\frac{\delta\mathcal E}{\delta\rho}=U'(\rho)+V+W*\rho.
\end{align*}
Since $t\in(0,T)$ was arbitrary, this gives for the curve $(\rho_t)_{t\in(0,T)}$ the time-dependent potential
\begin{align*}
\Phi_t:\mathbb R^n&\to\mathbb R
\end{align*}
defined by
\begin{align*}
\Phi_t(x)=U'(\rho_t(x))+V(x)+(W*\rho_t)(x).
\end{align*}
By hypothesis, $\Phi_t$ is $C^1$ in the spatial variable.
[/step]
[step:Substitute the Wasserstein gradient velocity into the continuity equation]
In the formal $2$-Wasserstein Riemannian calculus, the gradient-flow velocity associated to the first variation $\Phi_t=\delta\mathcal E/\delta\rho_t$ is
\begin{align*}
v_t:\mathbb R^n&\to\mathbb R^n
\end{align*}
defined by
\begin{align*}
v_t(x)=-\nabla\Phi_t(x)=-\nabla\left(U'(\rho_t(x))+V(x)+(W*\rho_t)(x)\right).
\end{align*}
The formal gradient-flow curve is therefore required to satisfy the continuity equation
\begin{align*}
\partial_t\rho_t+\nabla\cdot(\rho_t v_t)=0.
\end{align*}
Substituting the displayed formula for $v_t$ gives
\begin{align*}
\partial_t\rho_t+\nabla\cdot\left(-\rho_t\nabla\left(U'(\rho_t)+V+W*\rho_t\right)\right)=0.
\end{align*}
Equivalently,
\begin{align*}
\partial_t\rho_t=\nabla\cdot\left(\rho_t\nabla\left(U'(\rho_t)+V+W*\rho_t\right)\right).
\end{align*}
This is exactly the asserted formal aggregation-diffusion Wasserstein gradient-flow equation.
[/step]
