In this article, we propose and analyze an effective Hessian recovery strategy for the Lagrangian finite element method of arbitrary order. We prove that the proposed Hessian recovery method preserves polynomials of degree
on general unstructured meshes and superconverges at a rate of
on mildly structured meshes. In addition, the method is proved to be ultraconvergent (two orders higher) for the translation invariant finite element space of any order. Numerical examples are presented to support our theoretical results.