Summary: Optimization on Riemannian manifolds widely arises in eigenvalue computation, density functional theory, Bose-Einstein condensates, low rank nearest correlation, image registration, signal processing, and so on. We propose an adaptive quadratically regularized Newton method which approximates the original objective function by the second-order Taylor expansion in Euclidean space but keeps the Riemannian manifold constraints. The regularization term in the objective function of the subproblem enables us to utilize a Cauchy-point-like condition as the standard trust-region method for proving global convergence. The subproblem can be solved inexactly either by first-order methods or by performing corresponding Riemannian Newton-type steps. In the latter case, we can further take advantage of negative curvature directions. Both global convergence and superlinear local convergence are guaranteed under mild conditions. Extensive computational experiments and comparisons with other state-of-the-art methods indicate that the proposed algorithm is very promising.