The parabolic approximation, applied to the propagation of time harmonic scalar waves, replaces the governing Helmholtz equation with a SchrOdinger equation. The primary computational advantage of this approximation is that it is first order in the range coordinate. The validity of the approximation limits variations in an inhomogeneous wave number field to be both small, over the total range of the experiment, and slow, as measured on a length scale determined by an averaged wavelength. Thus, the validity of the approximation can be severely strained in a number of applications; e.g., it is difficult to justify its application in the presence of regions of rapid, even discontinuous, changes near a surface that runs in the range direction. In this paper we consider the derivation of extended parabolic wave theories which retain the feature that they are first order in the range coordinate. Distinguishing the extended theories from the ordinary one, i.e. the SchrOdinger equation, is the manner in which the cross-range coordinates enter. In the ordinary theory, they appear via the two-dimensional Laplacian, a differential or local operator; in the extended theories they appear via non-local operators. Obtaining explicit forms for specific nonlocal cross-range operators requires the solution of factor-ordering problems which have their counterpart in the correspondence between quantized theories governing the evolution of dynamical systems and their related classical formulations. The operator construction can be related to the construction of both coordinate and phase space path integrals for the appropriate wave propagator and further provides for the formulation and exact solution of a multidimensional nonlinear inverse problem appropriate for ocean acoustic and seismic modelling.