[proofplan] The proof has two parts. First, for a smooth minimizer of the constrained action, the smooth Lagrange multiplier principle for the no-flux continuity equation gives a multiplier $\phi$; velocity variations force the minimising velocity to be $\nabla\phi$, while density variations force the Hamilton-Jacobi equation up to a time-dependent constant. This constant is removed by the freedom to add a function of time to $\phi$. Conversely, the Hamilton-Jacobi inequality gives the quadratic-cost dual lower bound, while pairing the continuity equation with $\phi$ gives equality in that bound. The constant kinetic energy hypothesis then converts optimality on the full interval into the constant-speed geodesic identity on every subinterval. [/proofplan]

[step:Use the stated smooth constrained first-order condition for the no-flux action] By the regularity hypothesis in the theorem statement, there exists a multiplier $\psi\in C^1([0,1]\times\overline U;\mathbb R)$ such that $(\rho,v,\psi)$ is stationary for \begin{align*} \mathscr L[\rho,v,\psi] := \int_0^1\int_U \rho_t |v_t|^2\,d\mathcal L^n\,d\mathcal L^1(t) + 2\int_0^1\int_U \psi_t\left(\partial_t\rho_t+\nabla\cdot(\rho_t v_t)\right)\,d\mathcal L^n\,d\mathcal L^1(t), \end{align*} under all $C^1$ density variations with zero endpoint values and zero spatial mean at each time, and all $C^1$ velocity variations tangent to $\partial U$. Here $\psi_t:U\to\mathbb R$ denotes $x\mapsto\psi(t,x)$. The factor $2$ is only a normalization convention and is chosen so that the velocity variation gives $v_t=\nabla\psi_t$ without an additional factor. [/step]

[step:Vary the velocity to obtain a potential representation]Let \begin{align*} w:[0,1]\times\overline U&\to\mathbb R^n \end{align*} be a $C^1$ vector field such that $w_t\cdot\nu=0$ on $\partial U$ for every $t\in[0,1]$. For $\varepsilon\in\mathbb R$ sufficiently small, define the velocity variation \begin{align*} v^\varepsilon_t:U&\to\mathbb R^n, \end{align*} \begin{align*} x&\mapsto v_t(x)+\varepsilon w_t(x). \end{align*} Stationarity of $\mathscr L$ in the velocity variable gives \begin{align*} 0 = \frac{d}{d\varepsilon}\bigg|_{\varepsilon=0}\mathscr L[\rho,v^\varepsilon,\psi]. \end{align*} Differentiating under the integral sign gives \begin{align*} 0 = 2\int_0^1\int_U \rho_t v_t\cdot w_t\,d\mathcal L^n\,d\mathcal L^1(t) + 2\int_0^1\int_U \psi_t\nabla\cdot(\rho_t w_t)\,d\mathcal L^n\,d\mathcal L^1(t). \end{align*} For each fixed $t$, the divergence theorem on the $C^1$ domain $U$ gives, where $\mathcal H^{n-1}$ denotes the $(n-1)$-dimensional Hausdorff measure on $\partial U$, \begin{align*} \int_U \psi_t\nabla\cdot(\rho_t w_t)\,d\mathcal L^n = -\int_U \rho_t\nabla\psi_t\cdot w_t\,d\mathcal L^n + \int_{\partial U}\psi_t\rho_t w_t\cdot\nu\,d\mathcal H^{n-1}. \end{align*} The boundary integral vanishes because $w_t\cdot\nu=0$ on $\partial U$. Hence \begin{align*} \int_0^1\int_U \rho_t (v_t-\nabla\psi_t)\cdot w_t\,d\mathcal L^n\,d\mathcal L^1(t)=0 \end{align*} for every such $w$. Since $\rho_t>0$ on $U$ and $w$ is arbitrary in the interior, the fundamental lemma of the calculus of variations gives \begin{align*} v_t=\nabla\psi_t \end{align*} on $U$ for every $t\in[0,1]$.[/step]

[guided]The purpose of the velocity variation is to identify the part of $v_t$ that actually contributes to the minimum. We introduce an arbitrary smooth vector field \begin{align*} w:[0,1]\times\overline U&\to\mathbb R^n \end{align*} with $w_t\cdot\nu=0$ on $\partial U$. This tangency condition is exactly what is needed so that the varied velocity \begin{align*} v^\varepsilon_t:U&\to\mathbb R^n, \end{align*} \begin{align*} x&\mapsto v_t(x)+\varepsilon w_t(x) \end{align*} still satisfies the boundary condition $v^\varepsilon_t\cdot\nu=0$ on $\partial U$. Stationarity in the velocity variable says \begin{align*} 0 = \frac{d}{d\varepsilon}\bigg|_{\varepsilon=0}\mathscr L[\rho,v^\varepsilon,\psi]. \end{align*} The derivative of the kinetic term is \begin{align*} \frac{d}{d\varepsilon}\bigg|_{\varepsilon=0} \int_0^1\int_U \rho_t|v_t+\varepsilon w_t|^2\,d\mathcal L^n\,d\mathcal L^1(t) = 2\int_0^1\int_U \rho_t v_t\cdot w_t\,d\mathcal L^n\,d\mathcal L^1(t). \end{align*} The derivative of the constraint term is \begin{align*} 2\int_0^1\int_U \psi_t\nabla\cdot(\rho_t w_t)\,d\mathcal L^n\,d\mathcal L^1(t). \end{align*} Thus \begin{align*} 0 = 2\int_0^1\int_U \rho_t v_t\cdot w_t\,d\mathcal L^n\,d\mathcal L^1(t) + 2\int_0^1\int_U \psi_t\nabla\cdot(\rho_t w_t)\,d\mathcal L^n\,d\mathcal L^1(t). \end{align*} We now integrate by parts in the spatial variable. For each fixed $t$, the map $x\mapsto\rho_t(x)w_t(x)$ is $C^1$ on $\overline U$, and $U$ has $C^1$ boundary, so the divergence theorem applies; here $\mathcal H^{n-1}$ denotes the $(n-1)$-dimensional Hausdorff measure on $\partial U$: \begin{align*} \int_U \psi_t\nabla\cdot(\rho_t w_t)\,d\mathcal L^n = -\int_U \rho_t\nabla\psi_t\cdot w_t\,d\mathcal L^n + \int_{\partial U}\psi_t\rho_t w_t\cdot\nu\,d\mathcal H^{n-1}. \end{align*} The boundary term is zero because $w_t\cdot\nu=0$ on $\partial U$. Substituting this into the stationarity identity gives \begin{align*} \int_0^1\int_U \rho_t (v_t-\nabla\psi_t)\cdot w_t\,d\mathcal L^n\,d\mathcal L^1(t)=0. \end{align*} Because the variation $w$ is arbitrary in the interior of $U$ and because $\rho_t(x)>0$ for every $x\in U$, the only possible conclusion is \begin{align*} v_t=\nabla\psi_t \end{align*} on $U$ for every $t\in[0,1]$. Positivity of $\rho_t$ is essential here: it lets us divide out the weight in the variational identity and obtain a pointwise statement about the velocity field itself.[/guided]

[step:Vary the density to obtain the Hamilton-Jacobi equation up to a time function] Let \begin{align*} \sigma:[0,1]\times\overline U&\to\mathbb R \end{align*} be a $C^1$ function satisfying $\sigma_0=0$, $\sigma_1=0$, and \begin{align*} \int_U\sigma_t\,d\mathcal L^n=0 \end{align*} for every $t\in[0,1]$. For sufficiently small $\varepsilon\in\mathbb R$, define \begin{align*} \rho^\varepsilon_t:U&\to(0,\infty), \end{align*} \begin{align*} x&\mapsto\rho_t(x)+\varepsilon\sigma_t(x). \end{align*} Stationarity in the density variable gives \begin{align*} 0 = \frac{d}{d\varepsilon}\bigg|_{\varepsilon=0}\mathscr L[\rho^\varepsilon,v,\psi]. \end{align*} Differentiating gives \begin{align*} 0 = \int_0^1\int_U \sigma_t|v_t|^2\,d\mathcal L^n\,d\mathcal L^1(t) + 2\int_0^1\int_U \psi_t\partial_t\sigma_t\,d\mathcal L^n\,d\mathcal L^1(t) + 2\int_0^1\int_U \psi_t\nabla\cdot(\sigma_t v_t)\,d\mathcal L^n\,d\mathcal L^1(t). \end{align*} Integration by parts in time has no endpoint contribution because $\sigma_0=\sigma_1=0$, and integration by parts in space has no boundary contribution because $v_t\cdot\nu=0$ on $\partial U$. Hence \begin{align*} 0 = \int_0^1\int_U \sigma_t\left( |v_t|^2-2\partial_t\psi_t-2v_t\cdot\nabla\psi_t \right)\,d\mathcal L^n\,d\mathcal L^1(t). \end{align*} Using $v_t=\nabla\psi_t$, this becomes \begin{align*} 0 = -\int_0^1\int_U \sigma_t\left( 2\partial_t\psi_t+|\nabla\psi_t|^2 \right)\,d\mathcal L^n\,d\mathcal L^1(t). \end{align*} Define \begin{align*} g:[0,1]\times U&\to\mathbb R, \end{align*} \begin{align*} (t,x)&\mapsto 2\partial_t\psi_t(x)+|\nabla\psi_t(x)|^2. \end{align*} The preceding identity says \begin{align*} \int_0^1\int_U \sigma_t g_t\,d\mathcal L^n\,d\mathcal L^1(t)=0 \end{align*} for every $C^1$ function $\sigma$ with zero endpoint values and zero spatial mean at each time. Applying the fundamental lemma first with variations supported in an arbitrary time subinterval of $(0,1)$ and then with arbitrary zero-mean spatial variations gives that, for each $t\in(0,1)$, the function $g_t$ is constant on $U$. Since $g$ is continuous on $[0,1]\times U$, the same conclusion extends to $t=0$ and $t=1$ by taking one-sided limits. Thus there exists a continuous function \begin{align*} c:[0,1]&\to\mathbb R \end{align*} such that \begin{align*} \partial_t\psi_t+\frac12|\nabla\psi_t|^2=c(t) \end{align*} on $U$ for every $t\in[0,1]$. [/step]

[step:Normalize the potential by adding a function of time] Define \begin{align*} a:[0,1]&\to\mathbb R, \end{align*} \begin{align*} t&\mapsto -\int_0^t c(r)\,d\mathcal L^1(r), \end{align*} and define \begin{align*} \phi:[0,1]\times\overline U&\to\mathbb R, \end{align*} \begin{align*} (t,x)&\mapsto \psi(t,x)+a(t). \end{align*} Then $\nabla\phi_t=\nabla\psi_t=v_t$ on $U$. Also, \begin{align*} \partial_t\phi_t+\frac12|\nabla\phi_t|^2 = \partial_t\psi_t+a'(t)+\frac12|\nabla\psi_t|^2 = c(t)-c(t) = 0 \end{align*} on $U$ for every $t\in[0,1]$. If $\widetilde\phi\in C^1([0,1]\times\overline U;\mathbb R)$ is another potential satisfying $\nabla\widetilde\phi_t=v_t=\nabla\phi_t$ on $U$, then \begin{align*} \nabla(\widetilde\phi_t-\phi_t)=0 \end{align*} on $U$. Since $U$ is convex, it is connected. Therefore $\widetilde\phi_t-\phi_t$ is constant on $U$ for each fixed $t$. This proves uniqueness up to addition of a function of $t$. [/step]

[step:Use Hamilton-Jacobi admissibility to get the dynamic dual lower bound] We now prove the converse. Let $0\le s<t\le1$. Define \begin{align*} K_{s,t}:\overline U\times\overline U&\to[0,\infty), \end{align*} \begin{align*} (x,y)&\mapsto \frac{|x-y|^2}{2(t-s)}. \end{align*} The global Hamilton-Jacobi admissibility condition says \begin{align*} \phi_t(x)-\phi_s(y)\le K_{s,t}(x,y) \end{align*} for every $x,y\in\overline U$. Here $\Pi(\mu_t,\mu_s)$ denotes the set of Borel probability measures on $\overline U\times\overline U$ whose first marginal is $\mu_t$ and whose second marginal is $\mu_s$. Therefore, for every coupling $\pi\in\Pi(\mu_t,\mu_s)$, \begin{align*} \int_U\phi_t\,d\mu_t-\int_U\phi_s\,d\mu_s = \int_{\overline U\times\overline U}(\phi_t(x)-\phi_s(y))\,d\pi(x,y) \le \int_{\overline U\times\overline U}K_{s,t}(x,y)\,d\pi(x,y). \end{align*} Taking the infimum over all $\pi\in\Pi(\mu_t,\mu_s)$ gives \begin{align*} \int_U\phi_t\rho_t\,d\mathcal L^n-\int_U\phi_s\rho_s\,d\mathcal L^n \le \frac{1}{2(t-s)}W_2^2(\mu_s,\mu_t). \end{align*} Equivalently, \begin{align*} W_2^2(\mu_s,\mu_t) \ge 2(t-s)\left( \int_U\phi_t\rho_t\,d\mathcal L^n-\int_U\phi_s\rho_s\,d\mathcal L^n \right). \end{align*} [/step]

[step:Pair the continuity equation with the potential]Define \begin{align*} F:[0,1]&\to\mathbb R, \end{align*} \begin{align*} r&\mapsto \int_U\phi_r\rho_r\,d\mathcal L^n. \end{align*} Since $\rho$ and $\phi$ have the stated regularity, differentiating under the integral sign gives \begin{align*} F'(r) = \int_U \partial_t\phi_r\,\rho_r\,d\mathcal L^n + \int_U \phi_r\,\partial_t\rho_r\,d\mathcal L^n. \end{align*} Using the continuity equation $\partial_t\rho_r=-\nabla\cdot(\rho_r\nabla\phi_r)$, we get \begin{align*} F'(r) = \int_U \partial_t\phi_r\,\rho_r\,d\mathcal L^n - \int_U \phi_r\nabla\cdot(\rho_r\nabla\phi_r)\,d\mathcal L^n. \end{align*} The divergence theorem gives, where $\mathcal H^{n-1}$ denotes the $(n-1)$-dimensional Hausdorff measure on $\partial U$, \begin{align*} -\int_U \phi_r\nabla\cdot(\rho_r\nabla\phi_r)\,d\mathcal L^n = \int_U \rho_r|\nabla\phi_r|^2\,d\mathcal L^n - \int_{\partial U}\phi_r\rho_r\nabla\phi_r\cdot\nu\,d\mathcal H^{n-1}. \end{align*} The boundary term is zero by the no-flux condition. Therefore \begin{align*} F'(r) = \int_U\left(\partial_t\phi_r+|\nabla\phi_r|^2\right)\rho_r\,d\mathcal L^n. \end{align*} Using the Hamilton-Jacobi equation, \begin{align*} \partial_t\phi_r+\frac12|\nabla\phi_r|^2=0, \end{align*} we obtain \begin{align*} F'(r) = \frac12\int_U|\nabla\phi_r|^2\rho_r\,d\mathcal L^n = \frac12 E(r). \end{align*}[/step]

[guided]This step explains why the expression involving the endpoint potentials is the correct dual quantity. Define \begin{align*} F:[0,1]&\to\mathbb R, \end{align*} \begin{align*} r&\mapsto \int_U\phi_r\rho_r\,d\mathcal L^n. \end{align*} We differentiate $F$ using the product rule under the integral sign. This is justified by the assumed $C^1$ regularity in time of $\rho$ and $\phi$, together with compactness of $\overline U$. Thus \begin{align*} F'(r) = \int_U \partial_t\phi_r\,\rho_r\,d\mathcal L^n + \int_U \phi_r\,\partial_t\rho_r\,d\mathcal L^n. \end{align*} The continuity equation says that mass is transported by the velocity $\nabla\phi_r$: \begin{align*} \partial_t\rho_r=-\nabla\cdot(\rho_r\nabla\phi_r). \end{align*} Substituting this identity into the derivative of $F$ gives \begin{align*} F'(r) = \int_U \partial_t\phi_r\,\rho_r\,d\mathcal L^n - \int_U \phi_r\nabla\cdot(\rho_r\nabla\phi_r)\,d\mathcal L^n. \end{align*} Now we integrate by parts in space. Since $U$ has $C^1$ boundary and the functions are smooth enough up to the boundary, the divergence theorem gives, where $\mathcal H^{n-1}$ denotes the $(n-1)$-dimensional Hausdorff measure on $\partial U$, \begin{align*} -\int_U \phi_r\nabla\cdot(\rho_r\nabla\phi_r)\,d\mathcal L^n = \int_U \rho_r|\nabla\phi_r|^2\,d\mathcal L^n - \int_{\partial U}\phi_r\rho_r\nabla\phi_r\cdot\nu\,d\mathcal H^{n-1}. \end{align*} The boundary term vanishes exactly because of the no-flux condition \begin{align*} \rho_r\nabla\phi_r\cdot\nu=0 \end{align*} on $\partial U$. Therefore \begin{align*} F'(r) = \int_U\left(\partial_t\phi_r+|\nabla\phi_r|^2\right)\rho_r\,d\mathcal L^n. \end{align*} Finally we use the Hamilton-Jacobi equation \begin{align*} \partial_t\phi_r+\frac12|\nabla\phi_r|^2=0. \end{align*} It converts the integrand into exactly one half of the kinetic energy density: \begin{align*} \partial_t\phi_r+|\nabla\phi_r|^2 = \frac12|\nabla\phi_r|^2. \end{align*} Hence \begin{align*} F'(r) = \frac12\int_U|\nabla\phi_r|^2\rho_r\,d\mathcal L^n = \frac12E(r). \end{align*} This identity is the equality case in the dynamic duality argument: the potential $\phi$ does not merely give a lower bound; along the proposed curve it extracts exactly the kinetic action.[/guided]

[step:Identify the full interval action with the Wasserstein distance] Integrating the identity $F'(r)=E(r)/2$ from $0$ to $1$ gives \begin{align*} \int_U\phi_1\rho_1\,d\mathcal L^n-\int_U\phi_0\rho_0\,d\mathcal L^n = \frac12\int_0^1E(r)\,d\mathcal L^1(r). \end{align*} This is also assumed in the statement, and it is consistent with the preceding computation. Applying the dual lower bound from the previous step with $s=0$ and $t=1$ gives \begin{align*} W_2^2(\mu_0,\mu_1) \ge 2\left( \int_U\phi_1\rho_1\,d\mathcal L^n-\int_U\phi_0\rho_0\,d\mathcal L^n \right) = \int_0^1E(r)\,d\mathcal L^1(r). \end{align*} On the other hand, the curve $(\rho_r,\nabla\phi_r)$ is an admissible no-flux continuity-equation pair joining $\mu_0$ to $\mu_1$. Because $U$ is convex, $\overline U$ is geodesically convex for the Euclidean distance, and the no-flux boundary condition makes the continuity equation valid distributionally after viewing the measures as Borel probability measures on $\mathbb R^n$ supported in $\overline U$. The curve is narrowly continuous by the $C^1$ regularity of $\rho$, and its action is finite because $\overline U$ is compact and $\nabla\phi$ is continuous. Applying the Benamou-Brenier dynamic upper bound [citetheorem:9556] to this admissible curve gives \begin{align*} W_2^2(\mu_0,\mu_1) \le \int_0^1E(r)\,d\mathcal L^1(r). \end{align*} Therefore \begin{align*} W_2^2(\mu_0,\mu_1) = \int_0^1E(r)\,d\mathcal L^1(r). \end{align*} [/step]

[step:Use constant kinetic energy to prove constant speed] Since $E$ is constant on $[0,1]$, there exists $e\in[0,\infty)$ such that \begin{align*} E(r)=e \end{align*} for every $r\in[0,1]$. The previous step gives \begin{align*} W_2^2(\mu_0,\mu_1)=e. \end{align*} Now fix $0\le s<t\le1$. Repeating the same dual lower bound and pairing argument on the interval $[s,t]$ gives \begin{align*} W_2^2(\mu_s,\mu_t) \ge 2(t-s)\left(F(t)-F(s)\right). \end{align*} Since $F'(r)=e/2$, we have \begin{align*} F(t)-F(s)=\frac12 e(t-s). \end{align*} Thus \begin{align*} W_2^2(\mu_s,\mu_t)\ge e(t-s)^2. \end{align*} The Benamou-Brenier dynamic upper bound [citetheorem:9556] applied to the restriction of the same admissible curve to $[s,t]$, after the affine time reparametrization of $[s,t]$ onto $[0,1]$, gives \begin{align*} W_2^2(\mu_s,\mu_t) \le (t-s)\int_s^t E(r)\,d\mathcal L^1(r) = e(t-s)^2. \end{align*} Therefore \begin{align*} W_2(\mu_s,\mu_t)=\sqrt e\,|t-s|. \end{align*} Because $\sqrt e=W_2(\mu_0,\mu_1)$, this is \begin{align*} W_2(\mu_s,\mu_t)=|t-s|W_2(\mu_0,\mu_1) \end{align*} for every $0\le s<t\le1$. Hence $(\mu_t)_{t\in[0,1]}$ is a constant-speed geodesic in $(\mathcal P_2(\overline U),W_2)$. [/step]