We study the scaling behavior of two-dimensional (2D) crystalline membranes in the flat phase by a renormalization group (RG) method and an ε-expansion. Generalization of the problem to non-integer dimensions, necessary to control the ε-expansion, is achieved by dimensional continuation of a well-known effective theory describing out-of-plane fluctuations coupled to phonon-mediated interactions via a scalar composite field, equivalent for small deformations to the local Gaussian curvature. The effective theory, which will be referred to as Gaussian curvature interaction (GCI) model, is equivalent to theories of elastic D-dimensional manifolds fluctuating in a (D + d c )-dimensional embedding space in the physical case D = 2 for arbitrary d c . For D = 2, instead, the GCI model is not equivalent to a direct dimensional continuation of elastic membrane theory and it defines an alternative generalization to generic internal dimension D. After decoupling interactions through a Hubbard-Stratonovich transformation, we study the GCI model by perturbative field-theoretic RG within the framework of an expansion in ε = (4 − D). We calculate explicitly RG functions at two-loop order and determine the exponent η characterizing the long-wavelength scaling of correlation functions to order ε 2 . The value of η at this order is shown to be insensitive to Feynman diagrams involving vertex corrections. As a consequence, the self-consistent screening approximation for the GCI model is shown to be exact to O(ε 2 ). In the physical case of a single out-of-plane displacement field, d c = 1, the O(ε 2 ) correction is suppressed by a small numerical prefactor. As a result, despite the large value of ε = 2, extrapolation of the first and second order results to D = 2 leads to very close numbers, η = 0.8 and η 0.795. The calculated exponent values are close to earlier reference results obtained by nonperturbative renormalization group, the self-consistent screening approximation and numerical simulations. These indications suggest that a perturbative analysis of the GCI model could provide an useful framework for accurate quantitative predictions of the scaling exponent even at D = 2.