Assume $\rho$ is probability density defined on $\mathbb{R}^d\times\mathbb{R}^d$. I am interested in the Wasserstein gradient flow of a functional: \begin{equation*} \mathcal{E}(\rho)=\iint_{\mathbb{R}^d\times\mathbb{R}^d}c(x,y)\rho(x,y)~dxdy+\rm{KL}(\rho_a|\varrho_a)+\rm{KL}(\rho_b|\varrho_b) \end{equation*} Here $\rho_a(x)=\int\rho(x,y)~dy$, $\rho_b(y)=\int\rho(x,y)~dx$ are marginals of $\rho$ and $\varrho_a$,$\varrho_b$ are given probability densities on $\mathbb{R}^d$. I tried to compute the Wasserstein gradient flow of this functional $\mathcal{E}$: $\frac{\partial \rho}{\partial t}=-\rm{grad}_\rho \mathcal{E}(\rho)$,this gives the following Fokker Planck Equation: \begin{equation*} \frac{\partial\rho}{\partial t}=\nabla\cdot\left(\rho\nabla\left(c(x,y)+\log\left(\frac{\rho_a(x)}{\varrho_a(x)}\right)+\log\left(\frac{\rho_b(y)}{\varrho_b(y)}\right)\right)\right) \quad\quad (*) \end{equation*} Here is my question: Is it possible to give a Stochastic Differential Equation whose probability density is exactly the equation (*)? If we directly formulate the SDE like: \begin{equation} d\left[\begin{array}{c}\mathbf{X}_t\\\mathbf{Y}_t\end{array}\right]=\left[\begin{array}{c}\nabla_x(c(\mathbf{X}_t,\mathbf{Y}_t)-\log\varrho_a(\mathbf{X}_t))\\\nabla_y(c(\mathbf{X}_t,\mathbf{Y}_t)-\log\varrho_b(\mathbf{Y}_t))\end{array}\right]dt-\left[\begin{array}{c}\nabla_x\log\rho_a(\mathbf{X}_t)\\\nabla_y\log\rho_b(\mathbf{Y}_t)\end{array}\right]dt \quad\quad (\#) \end{equation} It is clear that the first term \begin{equation*} \left[\begin{array}{c}\nabla_x(c(\mathbf{X}_t,\mathbf{Y}_t)-\log\varrho_a(\mathbf{X}_t))\\\nabla_y(c(\mathbf{X}_t,\mathbf{Y}_t)-\log\varrho_b(\mathbf{Y}_t))\end{array}\right]dt \end{equation*} in (#) can be treated as the drift term, I am now quite confused about the second term in (#). What I want to do is to introduce Brownian Motion to replace the second term. If we consider the case \begin{equation*} d\mathbf{X}_t=-\nabla_x\log\rho(\mathbf{X}_t)dt \end{equation*} where $\rho$ is the law of $\mathbf{X}_t$. Then the corresponding Fokker Planck is $\frac{\partial \rho}{\partial t}=\Delta\rho$. So now if we let $\tilde{\mathbf{X}}_t$ satisfies \begin{equation*} d\tilde{\mathbf{X}}_t=\sqrt{2}~d\mathbf{B}_t \quad \tilde{\mathbf{X}}_0\stackrel{d}{=}\mathbf{X}_0 \end{equation*} Here $\mathbf{B}_t$ is the standard Brownian Motion. We will know $\mathbf{X}_t$ and $\tilde{\mathbf{X}}_t$ have the same law. (i.e.$\mathbf{X}_t\stackrel{d}{=}\tilde{\mathbf{X}}_t$ for any $t\geq 0$) But we can't simply replace the second term in (#) by $\left[\begin{array}{c} \sqrt{2}d\mathbf{B}_t^{(1)} \\ \sqrt{2}d\mathbf{B}_t^{(2)} \end{array}\right]$. Otherwise, the corresponding Fokker Planck becomes \begin{equation*} \frac{\partial\rho}{\partial t}=\nabla\cdot\left(\rho\nabla\left(c(x,y)-\log\varrho_a(x)-\log\varrho_b(y)\right)\right)+\Delta \rho \end{equation*} which is definitely different from (*). Are there any experts who can help me on this?