[proofplan]
Fix an arbitrary compactly supported smooth formal Wasserstein geodesic $(\rho_t)_{t\in I}$ with velocity field $\nabla\phi_t$. The formal second-variation identity along a formal $W_2$ geodesic identifies the ordinary second derivative of $t\mapsto\mathcal F[\rho_t]$ with the formal Wasserstein Hessian evaluated on the instantaneous velocity. Since each $\phi_t$ is smooth and compactly supported, the assumed Hessian inequality applies at every interior time. The remaining point is the one-variable calculus fact that a continuous function whose restriction to the interior has nonnegative second derivative is convex on the whole interval.
[/proofplan]
[step:Apply the formal geodesic second-variation identity]
Let
\begin{align*}
f:I\to\mathbb R,\qquad t\mapsto \mathcal F[\rho_t].
\end{align*}
Because $\mathcal F$ is twice differentiable in the formal $W_2$ calculus and $(\rho_t)_{t\in I}$ is a smooth formal $W_2$ geodesic with velocity field $\nabla\phi_t$, the formal geodesic second-variation formula gives, for every $t$ in the interior $I^\circ$ of $I$,
\begin{align*}
f''(t)=\operatorname{Hess}_{W_2}\mathcal F(\rho_t)[\nabla\phi_t,\nabla\phi_t].
\end{align*}
Here the geodesic hypothesis is precisely what removes the formal covariant-acceleration term from the general second-variation formula.
[guided]
Define the scalar function
\begin{align*}
f:I\to\mathbb R,\qquad t\mapsto \mathcal F[\rho_t].
\end{align*}
The goal is to prove convexity of this one-variable function. To connect convexity to the Wasserstein hypothesis, we use the formal second-variation rule in the $W_2$ calculus.
For a general smooth curve in the formal Wasserstein manifold, differentiating $\mathcal F[\rho_t]$ twice has two contributions: the Hessian of $\mathcal F$ applied to the velocity and a term involving the formal covariant acceleration of the curve. In the present theorem, $(\rho_t)_{t\in I}$ is assumed to be a formal $W_2$ geodesic with velocity field $\nabla\phi_t$. By the definition of a formal geodesic in this calculus, the formal covariant acceleration vanishes. Therefore the second derivative reduces to the Hessian term:
\begin{align*}
f''(t)=\operatorname{Hess}_{W_2}\mathcal F(\rho_t)[\nabla\phi_t,\nabla\phi_t]
\end{align*}
for every $t\in I^\circ$.
This is the only place where the geodesic assumption is used. Without it, the sign of the Hessian alone would not control the full second derivative of $\mathcal F[\rho_t]$ along an arbitrary curve.
[/guided]
[/step]
[step:Use compact support of the potentials to get a nonnegative scalar second derivative]
Fix $t\in I^\circ$. By hypothesis, $\rho_t\in\mathcal P_{2}^{\infty}(\mathbb R^n)$. Also, by the compact-support assumption on the geodesic potentials, $\phi_t\in C_c^\infty(\mathbb R^n)$. Applying the assumed Hessian inequality with $\rho=\rho_t$ and $\phi=\phi_t$ gives
\begin{align*}
\operatorname{Hess}_{W_2}\mathcal F(\rho_t)[\nabla\phi_t,\nabla\phi_t]\ge 0.
\end{align*}
Combining this with the second-variation identity yields
\begin{align*}
f''(t)\ge 0
\end{align*}
for every $t\in I^\circ$.
[/step]
[step:Convert nonnegativity of the second derivative into convexity on the interval]
We first prove convexity on $I^\circ$. Let $a,b,c\in I^\circ$ with $a<b<c$. Since $f''\ge0$ on $I^\circ$, the derivative $f'$ is nondecreasing on $I^\circ$. Hence, for every $r\in[a,b]$ and every $s\in[b,c]$,
\begin{align*}
f'(r)\le f'(s).
\end{align*}
Using the one-dimensional fundamental theorem of calculus on the compact intervals $[a,b]$ and $[b,c]$, we obtain
\begin{align*}
\frac{f(b)-f(a)}{b-a}\le \frac{f(c)-f(b)}{c-b}.
\end{align*}
Indeed, the left-hand side is the average of $f'$ over $[a,b]$, and the right-hand side is the average of $f'$ over $[b,c]$; monotonicity of $f'$ gives the displayed inequality. Rearranging this slope inequality gives
\begin{align*}
f(b)\le \frac{c-b}{c-a}f(a)+\frac{b-a}{c-a}f(c).
\end{align*}
Thus $f$ is convex on $I^\circ$.
It remains only to include possible endpoints of $I$. The smoothness of the curve and the formal differentiability of $\mathcal F$ imply that $f$ is continuous on $I$. If $x,y\in I$ and $\lambda\in[0,1]$, choose sequences $(x_k)_{k\in\mathbb N}$ and $(y_k)_{k\in\mathbb N}$ in $I^\circ$ such that $x_k\to x$ and $y_k\to y$, whenever an endpoint approximation is needed. Applying convexity on $I^\circ$ to $x_k$ and $y_k$ and then passing to the limit by continuity gives
\begin{align*}
f((1-\lambda)x+\lambda y)\le (1-\lambda)f(x)+\lambda f(y).
\end{align*}
Therefore $f$ is convex on all of $I$.
[guided]
We now reduce the conclusion to ordinary one-variable calculus. We have already proved that
\begin{align*}
f''(t)\ge 0
\end{align*}
for every $t\in I^\circ$. The first consequence is that $f'$ is nondecreasing on $I^\circ$: if $r<s$ are points of $I^\circ$, then the fundamental theorem of calculus gives
\begin{align*}
f'(s)-f'(r)=\int_r^s f''(u)\,d\mathcal L^1(u)\ge 0.
\end{align*}
Thus $f'(r)\le f'(s)$ whenever $r<s$.
Now take $a,b,c\in I^\circ$ with $a<b<c$. Since $f'$ is nondecreasing, every value of $f'$ on $[a,b]$ is at most every value of $f'$ on $[b,c]$. Averaging this inequality over the two intervals gives
\begin{align*}
\frac{1}{b-a}\int_a^b f'(r)\,d\mathcal L^1(r)\le \frac{1}{c-b}\int_b^c f'(s)\,d\mathcal L^1(s).
\end{align*}
Applying the fundamental theorem of calculus to the two integrals yields
\begin{align*}
\frac{f(b)-f(a)}{b-a}\le \frac{f(c)-f(b)}{c-b}.
\end{align*}
This is exactly the monotonicity of secant slopes. Multiplying by the positive number $(b-a)(c-b)$ and rearranging gives
\begin{align*}
(c-a)f(b)\le (c-b)f(a)+(b-a)f(c).
\end{align*}
Dividing by $c-a>0$ gives
\begin{align*}
f(b)\le \frac{c-b}{c-a}f(a)+\frac{b-a}{c-a}f(c).
\end{align*}
This is the three-point form of convexity on $I^\circ$.
Finally, if $I$ has endpoints, the second derivative statement was only used on the interior $I^\circ$. The function $f$ is continuous on $I$ because $(\rho_t)_{t\in I}$ is smooth and $\mathcal F$ is differentiable along such curves in the formal calculus. Given $x,y\in I$ and $\lambda\in[0,1]$, approximate any endpoint among $x$ and $y$ by interior points. For sequences $(x_k)_{k\in\mathbb N}$ and $(y_k)_{k\in\mathbb N}$ in $I^\circ$ with $x_k\to x$ and $y_k\to y$, convexity on $I^\circ$ gives
\begin{align*}
f((1-\lambda)x_k+\lambda y_k)\le (1-\lambda)f(x_k)+\lambda f(y_k).
\end{align*}
Passing to the limit using continuity of $f$ gives
\begin{align*}
f((1-\lambda)x+\lambda y)\le (1-\lambda)f(x)+\lambda f(y).
\end{align*}
So $f$ is convex on the full interval $I$.
[/guided]
[/step]
[step:Identify the convex function with the functional along the geodesic]
By definition of $f$, the convexity just proved says exactly that
\begin{align*}
\mathcal F[\rho_{(1-\lambda)s+\lambda t}]\le (1-\lambda)\mathcal F[\rho_s]+\lambda\mathcal F[\rho_t]
\end{align*}
for all $s,t\in I$ and all $\lambda\in[0,1]$. Hence $t\mapsto\mathcal F[\rho_t]$ is convex on $I$. Since the compactly supported smooth formal Wasserstein geodesic was arbitrary, $\mathcal F$ is formally displacement convex along every such geodesic.
[/step]
