[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)$.[/step]

[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: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]