<p>This paper is concerned with the inverse elastic scattering problem to determine the shape and location of an elastic obstacle with the Neumann boundary condition in two or three dimensions. By establishing a one-to-one correspondence between the Herglotz wave function and its kernel, we introduce the far-field operator which is crucial in the factorization method. We present a theoretical factorization of the far-field operator and rigorously prove the properties of its associated operators involved in the factorization. Unlike the Dirichlet problem where the boundary integral operator of the single-layer potential involved in the factorization of the far-field operator is weakly singular, the boundary integral operator of the conormal derivative of the double-layer potential involved in the factorization of the far-field operator with Neumann boundary conditions is hypersingular, which forces us to prove that this operator is isomorphic using Fredholm’s theorem. Meanwhile, we present theoretical analyses of the factorization method for various illumination and measurement cases, including compression-wave illumination and compression-wave measurement, shear-wave illumination and shear-wave measurement, and full-wave illumination and full-wave measurement. In addition, we also consider the limited aperture problem and provide a rigorous theoretical analysis of the factorization method in this case. Numerous numerical experiments are carried out to demonstrate the effectiveness of the proposed method, and to analyze the influence of various factors, such as polarization direction, frequency, and multi-scale scatterers on the reconstructed results.</p>