In this paper we study Willmore Legendrian surfaces (that is Legendrian surfaces which are critical points of the Willmore functional). We use an equality proved in Luo ( arXiv:1211.4227v6 ) to get a relation between Willmore Legendrian surfaces and contact stationary Legendrian surfaces in $$\mathbb {S}^5$$ , and then we use this relation to prove a classification result for Willmore Legendrian spheres in $$\mathbb {S}^5$$ . We also get an integral inequality for Willmore Legendrian surfaces and in particular we prove that if the square length of the second fundamental form of a Willmore Legendrian surface in $$\mathbb {S}^5$$ belongs to [0, 2], then it must be 0 and L is totally geodesic or 2 and L is a flat minimal Legendrian tori, which generalizes the result of Yamaguchi et al. (Proc Am Math Soc 54:276–280, 1976). We also study variation of the Willmore functional among Legendrian surfaces in 5-dimensional Sasakian manifolds. Let $$\Sigma $$ be a closed surface and $$(M,\alpha ,g_\alpha ,J)$$ a 5-dimensional Sasakian manifold with a contact form $$\alpha $$ , an associated metric $$g_\alpha $$ and an almost complex structure J. Assume that $$f:\Sigma \mapsto M$$ is a Legendrian immersion. Then f is called a contact stationary Legendrian Willmore surface (in short, a csL Willmore surface) if it is a critical point of the Willmore functional under contact deformations. To investigate the existence of csL Willmore surfaces we introduce a higher order flow which preserves the Legendre condition and decreases the Willmore energy. As a first step we prove that this flow is well posed if $$(M,\alpha ,g_\alpha ,J)$$ is a Sasakian Einstein manifold, in particular $$\mathbb {S}^5$$ .