Intuition/elegant reason for why Langevin diffusion converges to $\exp(-U)$? Given a potential function $U: \mathbb{R}^n \to \mathbb{R}$, Langevin diffusion is gradient descent plus a Brownian motion term: $X' = -\nabla U(X) + \sqrt{2} \text{ }dW$.
It happens that the stationary distribution of Langevin diffusion is very nice: it is proportional to $\exp(-U)$. I can verify this by plugging $\exp(-U)$ into a PDE and checking that some second-order terms miraculously cancel. But seeing as $\exp(-U)$ is such a simple solution, I expect an elegant reason for it... or at the very least, some intuition for why the stationary distribution at $x$ should depend only on $U(x)$. Is there any intuition to be found here?
 A: The reason a Langevin diffusion leaves $\nu(x)=e^{-U(x)}$ invariant is because it is symmetric or reversible with respect to $\nu$.  In comparison to general diffusion processes, the ergodic properties of symmetric diffusions are not only easier to analyze, but it is often possible to obtain an explicit formula for the stationary density of a symmetric diffusion.  This symmetry property is a bit more transparent when the Langevin equation is written as $$
dX_t = \nabla \log \nu(X_t) dt + \sqrt{2} dW_t \;. \tag{$\star$}
$$ The symmetry property can then be easily verified by checking that the infinitesimal generator of the Langevin diffusion is symmetric in an $L^2$-inner product weighted by $\nu$.  This property of the Langevin diffusion is analogous to the detailed balance condition of a Markov chain.  It is also the basis of standard MCMC methods based on Langevin diffusions. In fact, many random-walk based MCMC methods (like random walk Metropolis) with target probability density proportional to $\nu(x)$, can weakly approximate the corresponding Langevin diffusion ($\star$).  This is remarkable because, e.g., the random walk Metropolis algorithm does not involve the gradient of $U$.
To read more about this connection, see the introduction and Theorem 5.2 of this paper.
A: Here is one way to think about this. Let $\rho(X,t)$ denote the transient joint density. Consider the free energy functional
$$F(\rho) := \mathbb{E}\left[U + \log \rho\right] = D_{\rm{KL}}\left(\rho||\exp(-U)\right) \geq 0,$$
where $D_{\rm{KL}}$ denotes the Kullback-Leibler divergence. This shows that $\exp(-U)$ is the minimizer of $F$.
Let $\zeta := 1+\log\rho + U$. The evolution of $\rho(X,t)$ is governed by the Fokker-Planck or Kolmogorov's forward PDE
$$\frac{\partial\rho}{\partial t} = \nabla\cdot\left(\rho\nabla\zeta\right).$$
Notice that
$$\frac{{\rm{d}}F}{{\rm{d}}t} = -\mathbb{E}\left[\|\nabla\zeta\|^{2}\right],$$
where the expectation operation $\mathbb{E}$ is with respect to the joint density $\rho(X,t)$. Therefore, $F$ is a Lyapunov functional for the above PDE, i.e., $F$ decreases along the solution trajectory $\rho(X,t)$, and is equal to zero at the stationary density $\propto \exp(-U)$. Another way to say the same is that the above PDE is gradient flow of $F$ with respect to the 2-Wasserstein metric; see e.g., Ch. 8.3 in Villani's book "Topics in Optimal Transportation".
