I am learning the idea of "gradient" of a functional in Otto's calculus. It is defined as follows.

Suppose the space we are thinking about is $(\mathcal{P}_{2,AC}(\mathbb{R}^d),W_2)$, the space of probability measures with finite second moment that is absolutely continuous w.r.t. Lebesgue measure, and equipped with 2-Wasserstein distance. So later on we use the density $\rho$ instead of measure in this space. The "tangent space" at $\rho$ is defined to be $T_{\rho}\mathcal{P}_{2,AC}(\mathbb{R}^d)=\{f:\int fdx=0\}$. A vector $v=\nabla\phi$ is said to be "coupled" with some $f\in T_{\rho}\mathcal{P}_{2,AC}(\mathbb{R}^d)$ if they satisfy the equation $$-\nabla\cdot(\rho\nabla\phi)=f.$$

This definition of "couple" is consistent with the idea of "continuity equation" in optimal transport. The Otto's metric tensor is then defined as $$<f,f'>_{\rho}:=\int \rho\nabla\phi\cdot\nabla\phi'dx,$$

with $\phi,\phi'$ coupled with $f,f'\in T_{\rho}\mathcal{P}_{2,AC}(\mathbb{R}^d)$ resp.

Having these in hand we can definethe gradient of a functional $F(\rho)$: $\mathop{grad} F(\rho)\in T_{\rho}\mathcal{P}_{2,AC}(\mathbb{R}^d)$ is the function such that $$<\mathop{grad}F(\rho),f>_{\rho}=d_{\rho}F(f)=\frac{d}{dt}\mid_{t=0}F(\rho_t),$$

for every $f\in T_{\rho}\mathcal{P}_{2,AC}(\mathbb{R}^d)$, and $\rho_t$ is any curve such that $\rho_0=\rho$ and $\rho'(0)=f$.

By this definition we can easily calculate out that eg. when $F(\rho)=\int\rho\log\rho dx$ being the entropy functional, $\mathop{grad}F(\rho)=-\Delta\rho$.

My question is now what is the gradient of the functional $F(\rho)=W_2^2(\rho dx,\eta dx)$ for any given $\eta dx\in \mathcal{P}_{2,AC}(\mathbb{R}^d)$. My guess is that, if $\phi_*$ is the **unique** Kantorovich potential(the solution to the dual Kantorovich problem between $\rho dx$ and $\eta dx$ with quadratic cost), then we may have $$\mathop{grad}W_2^2(\rho,\eta)=-\nabla\cdot(\rho\nabla\phi_*),$$

since if we represent the Waserstein distance in terms of dual problem, then $$W_2^2(\rho,\eta)=\max_{\phi(x)+\phi^c(y)\leq |x-y|^2}\int\phi\rho dx+\int\phi^c\eta dx.$$

If there is no maximum in the formula(which means $\phi$ is fixed), then the above formula surely holds. Now I am thinking that whether we can use the stability of optimal transport map or something to prove that the guess is correct, but I don't know how to make that work.