A comprehensive theory for interpreting two-dimensional Fourier transform (2D-FT) electron spin resonance (ESR) experiments that is based on the stochastic Liouville equation is presented. It encompasses the full range of motional rates from fast through very slow motions, and it also provides for microscopic as well as macroscopic molecular ordering. In these respects it is as sophisticated in its treatment of molecular dynamics as the theory currently employed for analyzing cw ESR spectra. The general properties of the pulse propagator superoperator, which describes the microwave pulses in Liouville space, are analyzed in terms of the coherence transfer pathways appropriate for COSY (correlation spectroscopy), SECSY (spin-echo correlation spectroscopy), and 2D-ELDOR (electron-electron double resonance) sequences wherein either the free-induction decay (FID) or echo decay is sampled. Important distinctions are made among the sources of inhomogeneous broadening, which include (a) incomplete spectral averaging in the slow-motional regime, (b) unresolved superhyperfine structure and related sources, and (c) microscopic molecular ordering but macroscopic disorder (MOMD). The differing effects these sources of inhomogeneous broadening have on the two mirror image coherence pathways observed in the dual quadrature 2D experiments, as well as on the auto vs crosspeaks of 2D-ELDOR, is described. The theory is applied to simulate experiments of nitroxide spin labels in complex fluids such as membrane vesicles, where the MOMD model applies and these distinctions are particularly relevant, in order to extract dynamic and ordering parameters. The recovery of homogeneous linewidths from FID-based COSY experiments on complex fluids with significant inhomogeneous broadening is also described. The theory is applied to the ultraslow motional regime, and a simple method is developed to determine rotational rates from the broadening of the autopeaks of the 2D-ELDOR spectra as a function of the mixing time, which is due to the development of "motional crosspeaks." The application of this method to recent experiments with nitroxide probes illustrates that rotational correlation times as slow as milliseconds may be measured. It is shown how 2D-ELDOR can be useful to distinguish between the cases of very slow motional (SM) rates with little or no ordering and of very high ordering (HO) but substantial motional rates even though the cw ESR spectra are virtually the same. The effects of motion and of microscopic ordering on the nuclear modulation patterns in 2D-FT-ESR are compared, and it is suggested that these effects could be utilized to further distinguish between SM and HO cases. Key aspects of the challenging computational problems are discussed, and algorithms are described which lead to significant reductions in computation time as needed to permit nonlinear least-squares fitting of the theory to experiments.