[proofplan] We prove the formula directly, without appealing to a general dynamical-plan theorem. First we check that the pushed-forward measures $\mu_t$ lie in $\mathcal{P}_p(\mathbb{R}^n)$ and have the prescribed endpoints. Then, for $0 \le s \le t \le 1$, we push the same optimal plan $\pi$ forward by $(T_s,T_t)$ to obtain a coupling of $\mu_s$ and $\mu_t$, whose transport cost is exactly $|t-s|^p$ times the original optimal cost. This gives the upper bound on $W_p(\mu_s,\mu_t)$, and the triangle inequality along the chain $\mu_0,\mu_s,\mu_t,\mu_1$ forces equality. [/proofplan]

[step:Verify that the interpolated measures lie in $\mathcal{P}_p(\mathbb{R}^n)$ and have the correct endpoints] For each $t \in [0,1]$, the map $T_t: \mathbb{R}^n \times \mathbb{R}^n \to \mathbb{R}^n$ is continuous, hence Borel measurable. Therefore the pushforward measure $\mu_t := (T_t)_\#\pi$ is a Borel probability measure on $\mathbb{R}^n$. We now verify the finite $p$-moment condition. Since $t \in [0,1]$, the convexity inequality $|a+b|^p \le 2^{p-1}(|a|^p+|b|^p)$ gives, for every $(x,y) \in \mathbb{R}^n \times \mathbb{R}^n$, \begin{align*} |T_t(x,y)|^p \le 2^{p-1}\bigl((1-t)^p |x|^p + t^p |y|^p\bigr). \end{align*} Since $(\operatorname{pr}_1)_\#\pi=\mu_0$ and $(\operatorname{pr}_2)_\#\pi=\mu_1$, we obtain \begin{align*} \int_{\mathbb{R}^n} |z|^p \, d\mu_t(z) = \int_{\mathbb{R}^n \times \mathbb{R}^n} |T_t(x,y)|^p \, d\pi(x,y). \end{align*} Thus \begin{align*} \int_{\mathbb{R}^n} |z|^p \, d\mu_t(z) \le 2^{p-1}\int_{\mathbb{R}^n \times \mathbb{R}^n} \bigl(|x|^p+|y|^p\bigr) \, d\pi(x,y). \end{align*} Using the marginal identities again, \begin{align*} \int_{\mathbb{R}^n} |z|^p \, d\mu_t(z) \le 2^{p-1}\left(\int_{\mathbb{R}^n} |x|^p \, d\mu_0(x) + \int_{\mathbb{R}^n} |y|^p \, d\mu_1(y)\right) < \infty. \end{align*} Hence $\mu_t \in \mathcal{P}_p(\mathbb{R}^n)$. At the endpoints, $T_0=\operatorname{pr}_1$ and $T_1=\operatorname{pr}_2$, so \begin{align*} \mu_0 = (T_0)_\#\pi \end{align*} and \begin{align*} \mu_1 = (T_1)_\#\pi. \end{align*} Thus the curve has the required endpoints. [/step]

[step:Push the optimal plan forward to couple two intermediate measures]Fix $s,t \in [0,1]$ with $s \le t$. Define the two-time interpolation map $S_{s,t}: \mathbb{R}^n \times \mathbb{R}^n \to \mathbb{R}^n \times \mathbb{R}^n$ by \begin{align*} S_{s,t}(x,y) = (T_s(x,y),T_t(x,y)). \end{align*} Let \begin{align*} \pi_{s,t} := (S_{s,t})_\#\pi. \end{align*} Then $\pi_{s,t}$ is a Borel probability measure on $\mathbb{R}^n \times \mathbb{R}^n$. Let $\operatorname{pr}_1,\operatorname{pr}_2: \mathbb{R}^n \times \mathbb{R}^n \to \mathbb{R}^n$ denote the coordinate projections. Since $\operatorname{pr}_1 \circ S_{s,t}=T_s$ and $\operatorname{pr}_2 \circ S_{s,t}=T_t$, the functoriality of pushforward measures gives \begin{align*} (\operatorname{pr}_1)_\#\pi_{s,t} = (T_s)_\#\pi = \mu_s. \end{align*} Similarly, \begin{align*} (\operatorname{pr}_2)_\#\pi_{s,t} = (T_t)_\#\pi = \mu_t. \end{align*} Therefore $\pi_{s,t} \in \Pi(\mu_s,\mu_t)$.[/step]

[guided]We need an admissible transport plan between $\mu_s$ and $\mu_t$. The construction should remember that both $\mu_s$ and $\mu_t$ were produced from the same pair $(x,y)$ sampled according to $\pi$. This is why we define a map that records both intermediate positions at once. For fixed $s,t \in [0,1]$ with $s \le t$, define the map $S_{s,t}: \mathbb{R}^n \times \mathbb{R}^n \to \mathbb{R}^n \times \mathbb{R}^n$ by \begin{align*} S_{s,t}(x,y) = (T_s(x,y),T_t(x,y)). \end{align*} This map is Borel measurable because $T_s$ and $T_t$ are continuous maps from $\mathbb{R}^n \times \mathbb{R}^n$ to $\mathbb{R}^n$. We define \begin{align*} \pi_{s,t} := (S_{s,t})_\#\pi. \end{align*} Thus $\pi_{s,t}$ is the law of the pair of random points obtained by starting at $x$, ending at $y$, and observing the straight-line interpolation at times $s$ and $t$. To prove that $\pi_{s,t}$ is a coupling of $\mu_s$ and $\mu_t$, we check its two marginals. Let $\operatorname{pr}_1,\operatorname{pr}_2: \mathbb{R}^n \times \mathbb{R}^n \to \mathbb{R}^n$ denote the coordinate projections. The first coordinate of $S_{s,t}(x,y)$ is $T_s(x,y)$, so \begin{align*} \operatorname{pr}_1 \circ S_{s,t} = T_s. \end{align*} Using the identity $(f \circ g)_\#\nu=f_\#(g_\#\nu)$ for pushforward measures, we get \begin{align*} (\operatorname{pr}_1)_\#\pi_{s,t} = (\operatorname{pr}_1)_\#(S_{s,t})_\#\pi = (\operatorname{pr}_1 \circ S_{s,t})_\#\pi = (T_s)_\#\pi = \mu_s. \end{align*} The same computation for the second coordinate gives \begin{align*} (\operatorname{pr}_2)_\#\pi_{s,t} = (\operatorname{pr}_2)_\#(S_{s,t})_\#\pi = (\operatorname{pr}_2 \circ S_{s,t})_\#\pi = (T_t)_\#\pi = \mu_t. \end{align*} Hence $\pi_{s,t}$ has first marginal $\mu_s$ and second marginal $\mu_t$, which is precisely the statement that $\pi_{s,t} \in \Pi(\mu_s,\mu_t)$.[/guided]

[step:Use the pushed-forward coupling to obtain the constant-speed upper bound] Because $\pi_{s,t}$ is a coupling of $\mu_s$ and $\mu_t$, the definition of $W_p$ gives \begin{align*} W_p(\mu_s,\mu_t)^p \le \int_{\mathbb{R}^n \times \mathbb{R}^n} |u-v|^p \, d\pi_{s,t}(u,v). \end{align*} By the definition of $\pi_{s,t}$ as the pushforward of $\pi$ under $S_{s,t}$, \begin{align*} \int_{\mathbb{R}^n \times \mathbb{R}^n} |u-v|^p \, d\pi_{s,t}(u,v) = \int_{\mathbb{R}^n \times \mathbb{R}^n} |T_s(x,y)-T_t(x,y)|^p \, d\pi(x,y). \end{align*} For every $(x,y) \in \mathbb{R}^n \times \mathbb{R}^n$, \begin{align*} T_t(x,y)-T_s(x,y) = (t-s)(y-x). \end{align*} Therefore \begin{align*} |T_t(x,y)-T_s(x,y)|^p = |t-s|^p |x-y|^p. \end{align*} Combining this with the optimality of $\pi$, \begin{align*} W_p(\mu_s,\mu_t)^p \le |t-s|^p \int_{\mathbb{R}^n \times \mathbb{R}^n} |x-y|^p \, d\pi(x,y). \end{align*} Since $\pi$ is optimal, \begin{align*} W_p(\mu_s,\mu_t)^p \le |t-s|^p W_p(\mu_0,\mu_1)^p. \end{align*} Taking $p$-th roots gives \begin{align*} W_p(\mu_s,\mu_t) \le |t-s| W_p(\mu_0,\mu_1). \end{align*} [/step]

[step:Apply the triangle inequality to force equality] For $0 \le s \le t \le 1$, apply the preceding upper bound to the pairs $(0,s)$, $(s,t)$, and $(t,1)$. This gives \begin{align*} W_p(\mu_0,\mu_s) \le s W_p(\mu_0,\mu_1). \end{align*} Also, \begin{align*} W_p(\mu_s,\mu_t) \le (t-s)W_p(\mu_0,\mu_1). \end{align*} And \begin{align*} W_p(\mu_t,\mu_1) \le (1-t)W_p(\mu_0,\mu_1). \end{align*} Since $W_p$ is the Wasserstein metric on $\mathcal{P}_p(\mathbb{R}^n)$, its triangle inequality gives \begin{align*} W_p(\mu_0,\mu_1) \le W_p(\mu_0,\mu_s)+W_p(\mu_s,\mu_t)+W_p(\mu_t,\mu_1). \end{align*} Substituting the three upper bounds yields \begin{align*} W_p(\mu_0,\mu_1) \le \bigl(s+(t-s)+(1-t)\bigr)W_p(\mu_0,\mu_1). \end{align*} Since $s+(t-s)+(1-t)=1$, every inequality in this chain must be an equality. In particular, \begin{align*} W_p(\mu_s,\mu_t) = (t-s)W_p(\mu_0,\mu_1). \end{align*} For arbitrary $s,t \in [0,1]$, exchanging $s$ and $t$ if necessary gives \begin{align*} W_p(\mu_s,\mu_t) = |t-s|W_p(\mu_0,\mu_1). \end{align*} Thus $(\mu_t)_{t \in [0,1]}$ is a constant-speed geodesic from $\mu_0$ to $\mu_1$, with interpolation formula $\mu_t=(T_t)_\#\pi=((1-t)\operatorname{pr}_1+t\operatorname{pr}_2)_\#\pi$. [/step]