[proofplan]
We differentiate the potential energy twice along the smooth compactly supported geodesic. The continuity equation converts the first derivative into the pairing of $\nabla V$ with the velocity field $\nabla\phi_t$. Differentiating this expression again gives one term from $\partial_t\nabla\phi_t$ and one term from $\partial_t\rho_t$; the Hamilton-Jacobi equation and compact-support integration by parts make the derivative-of-$\phi_t$ terms cancel. The only remaining contribution is the Hessian quadratic form of $V$, and positivity of this form gives the formal convexity conclusion.
[/proofplan]
[step:Differentiate the potential energy and use the continuity equation]
For $t\in I_0$, define the map $E:I_0\to\mathbb R$ by $E(t)=\mathcal V[\rho_t]$.
Since $\operatorname{supp}\rho_t\subset K$ for $t\in I_0$, and $V$ and $\rho$ are smooth, differentiating under the integral sign is justified on compact subintervals of $I_0$. Hence
\begin{align*}
E'(t)=\int_{\mathbb R^n}V(x)\partial_t\rho_t(x)\,d\mathcal L^n(x).
\end{align*}
Using the continuity equation,
\begin{align*}
\partial_t\rho_t=-\nabla\cdot(\rho_t\nabla\phi_t),
\end{align*}
we obtain
\begin{align*}
E'(t)=-\int_{\mathbb R^n}V(x)\nabla\cdot(\rho_t(x)\nabla\phi_t(x))\,d\mathcal L^n(x).
\end{align*}
The vector field $x\mapsto V(x)\rho_t(x)\nabla\phi_t(x)$ has compact support contained in $K$. Therefore compact-support integration by parts gives
\begin{align*}
E'(t)=\int_{\mathbb R^n}\nabla V(x)\cdot\nabla\phi_t(x)\rho_t(x)\,d\mathcal L^n(x).
\end{align*}
[guided]
We first isolate the quantity whose second derivative is the formal Hessian. Define the map $E:I_0\to\mathbb R$ by $E(t)=\mathcal V[\rho_t]$.
Because all densities $\rho_t$ with $t\in I_0$ are supported in the same compact set $K$, every integral over $\mathbb R^n$ is actually an integral over $K$. The functions $V$ and $\rho$ are smooth, so the map $t\mapsto V(x)\rho_t(x)$ is smoothly differentiable and compactly supported in $x$ uniformly for $t$ in compact subintervals of $I_0$. Thus differentiating under the integral sign gives
\begin{align*}
E'(t)=\int_{\mathbb R^n}V(x)\partial_t\rho_t(x)\,d\mathcal L^n(x).
\end{align*}
The continuity equation says that the time variation of the density is produced by transport with velocity $\nabla\phi_t$:
\begin{align*}
\partial_t\rho_t=-\nabla\cdot(\rho_t\nabla\phi_t).
\end{align*}
Substituting this identity into the derivative gives
\begin{align*}
E'(t)=-\int_{\mathbb R^n}V(x)\nabla\cdot(\rho_t(x)\nabla\phi_t(x))\,d\mathcal L^n(x).
\end{align*}
Now we move the divergence from $\rho_t\nabla\phi_t$ onto $V$. This produces no boundary term because $\rho_t$ has compact support contained in $K$, so the vector field $x\mapsto V(x)\rho_t(x)\nabla\phi_t(x)$ is compactly supported. The compact-support integration by parts identity therefore gives
\begin{align*}
-\int_{\mathbb R^n}V(x)\nabla\cdot(\rho_t(x)\nabla\phi_t(x))\,d\mathcal L^n(x)
=\int_{\mathbb R^n}\nabla V(x)\cdot\nabla\phi_t(x)\rho_t(x)\,d\mathcal L^n(x).
\end{align*}
Hence
\begin{align*}
E'(t)=\int_{\mathbb R^n}\nabla V(x)\cdot\nabla\phi_t(x)\rho_t(x)\,d\mathcal L^n(x).
\end{align*}
[/guided]
[/step]
[step:Differentiate the velocity pairing a second time]
Define the map $A:I_0\times\mathbb R^n\to\mathbb R$ by $A(t,x)=\nabla V(x)\cdot\nabla\phi_t(x)$.
From the previous step,
\begin{align*}
E'(t)=\int_{\mathbb R^n}A(t,x)\rho_t(x)\,d\mathcal L^n(x).
\end{align*}
Again using the common compact support of $\rho_t$, differentiation under the integral sign gives
\begin{align*}
E''(t)=\int_{\mathbb R^n}\nabla V(x)\cdot\nabla(\partial_t\phi_t)(x)\rho_t(x)\,d\mathcal L^n(x)+\int_{\mathbb R^n}\nabla V(x)\cdot\nabla\phi_t(x)\partial_t\rho_t(x)\,d\mathcal L^n(x).
\end{align*}
Substituting the continuity equation into the second term yields
\begin{align*}
E''(t)=\int_{\mathbb R^n}\nabla V(x)\cdot\nabla(\partial_t\phi_t)(x)\rho_t(x)\,d\mathcal L^n(x)-\int_{\mathbb R^n}A(t,x)\nabla\cdot(\rho_t(x)\nabla\phi_t(x))\,d\mathcal L^n(x).
\end{align*}
Since $x\mapsto A(t,x)\rho_t(x)\nabla\phi_t(x)$ has compact support, integration by parts gives
\begin{align*}
E''(t)=\int_{\mathbb R^n}\nabla V(x)\cdot\nabla(\partial_t\phi_t)(x)\rho_t(x)\,d\mathcal L^n(x)+\int_{\mathbb R^n}\nabla A(t,x)\cdot\nabla\phi_t(x)\rho_t(x)\,d\mathcal L^n(x).
\end{align*}
[/step]
[step:Use the geodesic equation to cancel the potential-derivative terms]
The Hamilton-Jacobi equation gives
\begin{align*}
\partial_t\phi_t+\frac{1}{2}|\nabla\phi_t|^2=c(t).
\end{align*}
Taking the spatial gradient, and using that $c(t)$ has no $x$-dependence, we obtain
\begin{align*}
\nabla(\partial_t\phi_t)(x)=-\nabla\left(\frac{1}{2}|\nabla\phi_t(x)|^2\right).
\end{align*}
Since $\phi_t$ is smooth,
\begin{align*}
\nabla\left(\frac{1}{2}|\nabla\phi_t(x)|^2\right)=D^2\phi_t(x)\nabla\phi_t(x).
\end{align*}
Thus
\begin{align*}
\nabla(\partial_t\phi_t)(x)=-D^2\phi_t(x)\nabla\phi_t(x).
\end{align*}
Also,
\begin{align*}
\nabla A(t,x)=D^2V(x)\nabla\phi_t(x)+D^2\phi_t(x)\nabla V(x).
\end{align*}
Substituting these two identities into the formula for $E''(t)$ gives
\begin{align*}
E''(t)=\int_{\mathbb R^n}\nabla V(x)\cdot\bigl(-D^2\phi_t(x)\nabla\phi_t(x)\bigr)\rho_t(x)\,d\mathcal L^n(x)+\int_{\mathbb R^n}\bigl(D^2V(x)\nabla\phi_t(x)+D^2\phi_t(x)\nabla V(x)\bigr)\cdot\nabla\phi_t(x)\rho_t(x)\,d\mathcal L^n(x).
\end{align*}
Because $D^2\phi_t(x)$ is symmetric, the two terms involving $D^2\phi_t(x)$ cancel pointwise. Therefore
\begin{align*}
E''(t)=\int_{\mathbb R^n}\nabla\phi_t(x)\cdot\bigl(D^2V(x)\nabla\phi_t(x)\bigr)\rho_t(x)\,d\mathcal L^n(x).
\end{align*}
[/step]
[step:Evaluate at $t=0$ and identify the formal Hessian]
Evaluating the preceding identity at $t=0$ gives
\begin{align*}
\left.\frac{d^2}{dt^2}\right|_{t=0}\mathcal V[\rho_t]
=\int_{\mathbb R^n}\nabla\phi_0(x)\cdot\bigl(D^2V(x)\nabla\phi_0(x)\bigr)\rho_0(x)\,d\mathcal L^n(x).
\end{align*}
By the formal Wasserstein Hessian convention for smooth Wasserstein geodesic variations, the left-hand side is
\begin{align*}
\operatorname{Hess}_{W_2}\mathcal V(\rho_0)[\nabla\phi_0,\nabla\phi_0].
\end{align*}
This proves the Hessian formula.
[/step]
[step:Deduce formal displacement convexity from positive semidefiniteness]
Assume that $D^2V(x)$ is positive semidefinite for every $x\in\mathbb R^n$. Then for every $x\in\mathbb R^n$,
\begin{align*}
\nabla\phi_0(x)\cdot\bigl(D^2V(x)\nabla\phi_0(x)\bigr)\ge 0.
\end{align*}
Since $\rho_0(x)\ge 0$, the Hessian formula gives
\begin{align*}
\operatorname{Hess}_{W_2}\mathcal V(\rho_0)[\nabla\phi_0,\nabla\phi_0]\ge 0.
\end{align*}
The same computation applies at any time along any smooth compactly supported Wasserstein geodesic in the stated class. Hence the second derivative of $t\mapsto\mathcal V[\rho_t]$ is nonnegative along such geodesics, which is precisely formal displacement convexity on this smooth compactly supported class.
[/step]
