Generalized master equations valid for the third order response of an optically driven multi-level electronic system are derived within Zwanzig projection formalism. Each of three time intervals of the response function is found to require specific master equation and projection operator.
Exact cumulant response functions for the harmonic profiles of the potential energy surfaces leading to GAUSSIAN spectral diffusion are reproduced. The proposed method accounts for the nonequilibrium state of bath at the border of intervals, and can be used to improve calculations of ultrafast non-linear spectra of energy transferring systems.