  We present a theoretical framework for higher-order correlation functions
involving multiple times and multiple points in a classical, many-body system.
Such higher-order correlation functions have attracted much interest recently
in the context of various forms of multi-dimensional spectroscopy, and have
found an intriguing application as proposed measures of dynamical
heterogeneities in structural glasses. The theoretical formalism is based upon
projection operator techniques which are used to isolate the slow time
evolution of dynamical variables by expanding the slowly-evolving component of
arbitrary variables in an infinite, ``multi-linear'' basis composed of the
products of slow variables of the system. Using the formalism, a formally exact
mode coupling theory is derived for multi-point and multi-time correlation
functions. The resulting expressions for higher-order correlation functions are
made tractable by applying a rigorous perturbation scheme, called the
N-ordering method, which is exact for systems with finite correlation lengths
in the thermodynamic limit. The theory is contrasted with standard mode
coupling theories in which the noise or fluctuating force appearing in the
generalized Langevin equation is assumed to be Gaussian, and it is demonstrated
that the non-Gaussian nature of the fluctuating forces leads to important
contributions to higher-order correlation functions. Finally, the higher-order
correlation functions are evaluated analytically for an ideal gas system for
which it is shown that the mode coupling theory is exact.
