  The identification of relevant collective coordinates is crucial for the
interpretation of coherent nonlinear spectroscopies of complex molecules and
liquids. Using an $\hbar$ expansion of Liouville space generating functions, we
show how to factorize multitime nonlinear response functions into products of
lower-order correlation functions of collective coordinates, and derive closed
expressions for linear, second and third order response functions. In addition
to providing systematic quantum corrections, $\hbar$ offers a convenient
bookkeeping device even for the purely classical response, since including
quantum fluctuations allows to circumvent the expensive computation of
stability matrices which is a major bottleneck in Molecular Dynamics (MD)
simulations. The existing classical simulation strategies, including
Mode-Coupling in $\mathbf{k}$ space and in real-space, Langevin equations, and
Instantaneous Normal Modes are compared from a unified viewpoint.
