[proofplan]
We use the standard formal Otto-calculus description of the Levi-Civita connection on the manifold of smooth positive probability densities. A smooth curve has velocity potential $\phi_t$ precisely when its density solves the continuity equation with velocity field $\nabla\phi_t$. The formal geodesic condition is vanishing covariant acceleration, which in gradient coordinates says that the gradient part of $\partial_t\nabla\phi_t+\nabla_{\nabla\phi_t}\nabla\phi_t$ vanishes. The Riemannian identity $\nabla_{\nabla\phi_t}\nabla\phi_t=\nabla(\frac12|\nabla\phi_t|_g^2)$ reduces this to a Hamilton-Jacobi equation up to an additive componentwise time-dependent gauge, and the converse follows by reversing these implications.
[/proofplan]
[step:Recall the formal Otto velocity and acceleration conventions]
Define the smooth positive probability-density manifold
\begin{align*}
\mathcal P^\infty(M)=\left\{\sigma\in C^\infty(M;(0,\infty)):\int_M \sigma\,d\operatorname{vol}_g=1\right\}.
\end{align*}
For $\sigma\in\mathcal P^\infty(M)$, the formal tangent space is
\begin{align*}
T_\sigma\mathcal P^\infty(M)=\left\{h\in C^\infty(M;\mathbb R):\int_M h\,d\operatorname{vol}_g=0\right\}.
\end{align*}
In Otto's formal Riemannian structure, a tangent vector $h\in T_\sigma\mathcal P^\infty(M)$ is represented by a potential $\psi\in C^\infty(M;\mathbb R)$, unique up to an additive function constant on each connected component of $M$, through
\begin{align*}
h=-\operatorname{div}_g(\sigma\nabla\psi).
\end{align*}
Thus the curve $t\mapsto \rho_t\operatorname{vol}_g$ has velocity field $v_t=\nabla\phi_t$ exactly when
\begin{align*}
\partial_t\rho_t+\operatorname{div}_g(\rho_t\nabla\phi_t)=0.
\end{align*}
By the formal Otto Levi-Civita convention specified in the theorem statement, the acceleration of the velocity field $\nabla\phi_t$ is represented by the gradient projection of
\begin{align*}
\partial_t\nabla\phi_t+\nabla^M_{\nabla\phi_t}\nabla\phi_t,
\end{align*}
where $\nabla^M$ denotes the Levi-Civita connection of $(M,g)$. Therefore the curve is a formal $W_2$-geodesic with velocity potential $\phi_t$ precisely when the continuity equation holds and
\begin{align*}
\partial_t\nabla\phi_t+\nabla^M_{\nabla\phi_t}\nabla\phi_t=0
\end{align*}
as a gradient vector field, equivalently up to the kernel of the gradient map on each connected component.
[/step]
[step:Convert vanishing formal acceleration into the Hamilton-Jacobi equation]
Assume first that $t\mapsto \rho_t\operatorname{vol}_g$ is a formal $W_2$-geodesic with velocity field $\nabla\phi_t$. By the preceding formal Otto convention, the continuity equation holds. It remains to identify the scalar equation for $\phi_t$.
For each $t\in(0,1)$, define
\begin{align*}
F_t:M&\to\mathbb R
\end{align*}
\begin{align*}
x&\mapsto \partial_t\phi_t(x)+\frac{1}{2}|\nabla\phi_t(x)|_g^2.
\end{align*}
Since $\phi$ is smooth and $g$ is smooth, $F_t\in C^\infty(M;\mathbb R)$. The Riemannian identity for gradients gives
\begin{align*}
\nabla^M_{\nabla\phi_t}\nabla\phi_t=\nabla\left(\frac{1}{2}|\nabla\phi_t|_g^2\right).
\end{align*}
Also, because $\phi$ is smooth on $[0,1]\times M$, time differentiation commutes with the spatial gradient:
\begin{align*}
\partial_t\nabla\phi_t=\nabla(\partial_t\phi_t).
\end{align*}
Hence vanishing formal acceleration is equivalent to
\begin{align*}
\nabla F_t=0.
\end{align*}
Therefore, on each connected component $M_\alpha$ of $M$, there exists a real number $c_\alpha(t)$ such that
\begin{align*}
F_t(x)=c_\alpha(t)
\end{align*}
for every $x\in M_\alpha$. To record the dependence on $t$, choose a point $x_\alpha\in M_\alpha$ and define $c_\alpha:(0,1)\to\mathbb R$ by $c_\alpha(t)=F_t(x_\alpha)$. Since $\phi$ and $g$ are smooth, the function $(t,x)\mapsto F_t(x)$ is smooth on $(0,1)\times M$, and hence each $c_\alpha$ is smooth on $(0,1)$.
[guided]
We now translate the vector equation for geodesics into a scalar equation for the potential. Fix $t\in(0,1)$ and define the smooth function
\begin{align*}
F_t:M&\to\mathbb R
\end{align*}
\begin{align*}
x&\mapsto \partial_t\phi_t(x)+\frac{1}{2}|\nabla\phi_t(x)|_g^2.
\end{align*}
This function is the natural candidate because the nonlinear acceleration term of a gradient vector field is itself a gradient. Indeed, for every smooth function $u:M\to\mathbb R$, the Levi-Civita connection gives the identity
\begin{align*}
\nabla^M_{\nabla u}\nabla u=\nabla\left(\frac{1}{2}|\nabla u|_g^2\right).
\end{align*}
Applying this with $u=\phi_t$ gives
\begin{align*}
\nabla^M_{\nabla\phi_t}\nabla\phi_t=\nabla\left(\frac{1}{2}|\nabla\phi_t|_g^2\right).
\end{align*}
Since $\phi\in C^\infty([0,1]\times M;\mathbb R)$, differentiating in $t$ commutes with taking the spatial gradient, so
\begin{align*}
\partial_t\nabla\phi_t=\nabla(\partial_t\phi_t).
\end{align*}
Adding the two identities, the formal acceleration vector field becomes
\begin{align*}
\partial_t\nabla\phi_t+\nabla^M_{\nabla\phi_t}\nabla\phi_t
=
\nabla\left(\partial_t\phi_t+\frac{1}{2}|\nabla\phi_t|_g^2\right)
=
\nabla F_t.
\end{align*}
The geodesic condition says that this acceleration vanishes in Otto's gradient-coordinate description. Thus $\nabla F_t=0$. A smooth function with zero gradient is constant on each connected component of $M$, so for every connected component $M_\alpha$ there is a real number $c_\alpha(t)$ satisfying
\begin{align*}
F_t(x)=c_\alpha(t)
\end{align*}
for all $x\in M_\alpha$. To see that this constant depends smoothly on time, choose a point $x_\alpha\in M_\alpha$ and set $c_\alpha(t)=F_t(x_\alpha)$. The map $(t,x)\mapsto F_t(x)$ is smooth because $\phi$ and $g$ are smooth, so $c_\alpha:(0,1)\to\mathbb R$ is smooth. This is exactly the Hamilton-Jacobi equation, except for the additive function of time that potentials cannot detect through their gradients.
[/guided]
[/step]
[step:Remove the componentwise additive gauge]
Let $\{M_\alpha\}_{\alpha\in A}$ be the finite family of connected components of the compact manifold $M$. Let $\mathcal L^1$ denote one-dimensional Lebesgue measure on $\mathbb R$. For each $\alpha\in A$, define
\begin{align*}
a_\alpha:(0,1)&\to\mathbb R
\end{align*}
\begin{align*}
t&\mapsto \int_{1/2}^{t} c_\alpha(s)\,d\mathcal L^1(s),
\end{align*}
and extend smoothly to $[0,1]$ after choosing any smooth endpoint extension. Define the gauge function
\begin{align*}
a:[0,1]\times M&\to\mathbb R
\end{align*}
\begin{align*}
(t,x)&\mapsto a_\alpha(t)\quad\text{if }x\in M_\alpha.
\end{align*}
Then $a(t,\cdot)$ is constant on each connected component, so
\begin{align*}
\nabla a(t,\cdot)=0.
\end{align*}
Set
\begin{align*}
\widetilde\phi:[0,1]\times M&\to\mathbb R
\end{align*}
\begin{align*}
(t,x)&\mapsto \phi(t,x)-a(t,x).
\end{align*}
Then $\nabla\widetilde\phi_t=\nabla\phi_t$, so the velocity field and the continuity equation are unchanged. Moreover, for $x\in M_\alpha$ and $t\in(0,1)$,
\begin{align*}
\partial_t\widetilde\phi_t(x)+\frac{1}{2}|\nabla\widetilde\phi_t(x)|_g^2
=
\partial_t\phi_t(x)-c_\alpha(t)+\frac{1}{2}|\nabla\phi_t(x)|_g^2
=
0.
\end{align*}
Thus, after the allowed additive gauge choice, the Hamilton-Jacobi equation holds on $M$.
[/step]
[step:Reverse the argument to obtain vanishing formal acceleration]
Conversely, suppose that after an allowed additive gauge choice the pair $(\rho_t,\phi_t)$ satisfies
\begin{align*}
\partial_t\rho_t+\operatorname{div}_g(\rho_t\nabla\phi_t)=0
\end{align*}
and
\begin{align*}
\partial_t\phi_t+\frac{1}{2}|\nabla\phi_t|_g^2=0
\end{align*}
for every $t\in(0,1)$ as smooth functions on $M$. The first equation says, by the definition of Otto velocity potentials, that $\nabla\phi_t$ is the velocity field of the density curve.
Taking the spatial gradient of the second equation gives
\begin{align*}
\nabla(\partial_t\phi_t)+\nabla\left(\frac{1}{2}|\nabla\phi_t|_g^2\right)=0.
\end{align*}
Using again the commutation $\nabla(\partial_t\phi_t)=\partial_t\nabla\phi_t$ and the Riemannian identity
\begin{align*}
\nabla\left(\frac{1}{2}|\nabla\phi_t|_g^2\right)=\nabla^M_{\nabla\phi_t}\nabla\phi_t,
\end{align*}
we obtain
\begin{align*}
\partial_t\nabla\phi_t+\nabla^M_{\nabla\phi_t}\nabla\phi_t=0.
\end{align*}
Hence the formal Otto covariant acceleration vanishes. Therefore $t\mapsto\rho_t\operatorname{vol}_g$ is a formal $W_2$-geodesic with velocity field $v_t=\nabla\phi_t$.
[/step]
