Summary: A new higher order finite element method for elliptic partial differential equations on a stationary smooth surface \(\Gamma\) is introduced and analyzed. We assume that \(\Gamma\) is characterized as the zero level of a level set function \(\phi\) and only a finite element approximation \(\phi_h\) (of degree \(k\geq 1\)) of \(\phi\) is known. For the discretization of the partial differential equation, finite elements (of degree \(m\geq 1\)) on a piecewise linear approximation of \(\Gamma\) are used. The discretization is lifted to \(\Gamma_h\), which denotes the zero level of \(\phi_h\), using a quasi-orthogonal coordinate system that is constructed by applying a gradient recovery technique to \(\phi_h\). A complete discretization error analysis is presented in which the error is split into a geometric error, a quadrature error, and a finite element approximation error. The main result is a \(H^1(\Gamma)\)-error bound of the form \(c(h^m+h^{k+1})\). Results of numerical experiments illustrate the higher order convergence of this method.