Here's an elementary type argument I learned, based on unique splitting $T=S+A$ of linear operators $T: \mathbb{R}^n \to \mathbb{R}^n$. 

(i) every linear operator $T$ can be uniquely decomposed $T=S+A$ where $S$ is symmetric ${}^tS=S$ and $A$ is antisymmetric ${}^tA=-A$. 

The meaning of symmetric and antisymmetric operators has a physical interpretation, where every force $T$ is factored as a radial force along the orthogonal eigenspaces of $S$, and the curl-type rotational force $A$. 

(ii) Now with respect euclidean quadratic cost $c(x,y)=-x\cdot y$, the $c$-optimal transportation plans are those plans which have zero rotational component, i.e. $T=S=\nabla \phi$ for a function $\phi$. It is elementary that $T$ is positive (semi) definite iff $\phi$ is (semi) convex. That quadratic-optimal $\pi$ are supported on the graph of gradients of convex potentials is proved by Kantorovich duality, and the form of the dual max program, which I express as: 

$$\min_{\pi} \int_{X\times Y} c(x,y) d\pi(x,y) \leq -\int_X \phi(x) d\sigma(x) + \int_Y \psi(y) d\tau(y), $$ where $\pi$ ranges over all couplings from a source measure $\sigma$ to target measure $\tau$, which we typically assume are absolutely continuous with respect to some background geometric measure. The potentials $\phi$, $\psi$ are convex and concave, respectively, and satisfy the pointwise inequality $-\phi(x)+\psi(y)\leq c(x,y)=-x.y$. 



In *physical* terms, I understand Brenier-McCann's polar factorization theorem as saying that quadratically optimal transport plans move particles along rotation-free trajectories. Every source particle $dx$ has a target $dy$ and is direction by the quadratic-optimal transport plan $\pi$ to "Keep straight, and don't twist or rotate until you arrive at your destination." 

Now you wanted a proof:

The cost $c(x,y)=-x.y$ satisfies the best hypotheses from optimal transport, especially the injectivity of the map $Y\to T_x X$ defined by $y\mapsto\nabla_x c(x,y)=-y$, for every $x\in X$. This is the (Twist) condition following standard terminology of McCann, Villani, etc.. From the Legendre Fenchel inequality $-\phi(x)+\psi(y) \leq c(x,y)$, we find equality holds iff $x$ belongs to the $c$-subdifferential of $\psi=\phi^c$ at $y$, where $\phi$ is a $c$-convex function. And differentiating the equality we find $-\nabla_x \phi(x)=\nabla_x c(x,y)$, which yields $y=\nabla_x \phi(x)$. Therefore the $c$-optimal transport map has the form $y=T(x)=\nabla_x \phi(x)$, where $\phi$ is the $c$-convex potential maximizing Kantorovich's dual program. Differentiating the equality case of $-\phi(x)+\psi(y) \leq c(x,y)$ yields $$D^2_{xx}\phi(x)=D(\nabla_x \phi(x))\geq 0.$$ And thus we find Brenier-McCann's theorem that quadratic-optimal transport maps have the form $y=T(x)=\nabla_x \phi(x)$ with respect to some convex potential $\phi$.