This paper establishes the theoretical foundation for statistical applications of an intriguing new type of spatial point processes called critical point processes. These point processes, residing in Euclidean space, consist of the critical points of latent smooth Gaussian random fields or of subsets of critical points like minima, saddle points etc. Despite of the simplicity of their definition, the mathematical analysis of critical point processes is non-trivial involving for example deep results on the geometry of random fields, Sobolev space theory, chaos expansions, and multiple Wiener-It{ô} integrals. We provide explicit expressions for fundamental moment characteristics used in spatial point process statistics like the intensity parameter, the pair correlation function, and higher order intensity functions. The crucial dependence structure (attraction or repulsiveness) of a critical point process is discussed in depth and is in particular related to the dimension of the points and the type of critical points (extrema, saddle points, or all of the critical points). We propose simulation strategies based on spectral methods or smoothing of grid-based simulations and show that resulting approximate critical point process simulations asymptotically converge to the exact critical point process distribution. Finally, under the increasing domain framework, we obtain asymptotic results for linear and bilinear statistics of a critical point process. In particular, we obtain a multivariate central limit theorem for the intensity parameter estimate and a modified version of Ripley's K-function.