A departure from Gaussian statistics is quickly met in the realistic case of mismatched ground-and excited state frequencies in harmonic vibrational modes (as present in the GBOM Hamiltonian), as well as in the presence of any anharmonicity. It can be shown analytically (see SI Secs. I and II) or numerically (SI Sec. III), that the introduction of higher order cumulants leads to an improved agreement with the exact response function in the short timescale limit, but that higher order terms become increasingly more volatile and divergent in the long timescale limit. In principle, a complete summation of the cumulant expansion to infinite order would lead to a cancellation of divergent terms and well-behaved response functions for all timescales. However, a truncation at finite order can yield unphysical lineshapes, as can be demonstrated even in the case of the harmonic GBOM Hamiltonian (SI Sec. II and Ref. 39), and only the second order cumulant approximation is guaranteed to yield linear response functions that do not diverge as t → ∞ for arbitrary Hamiltonians. This divergent behavior makes the benefit of evaluating higher order cumulants in realistic systems, such as by introducing a third order correction that can be computed in condensed phase systems from MD, questionable. In previous work, we have demonstrated numerically that monitoring the deviations of the energy gap fluctuations from Gaussian statistics, by calculating the skewness µ (3) of the distribution, can serve as an indicator of the importance of higher order cumulant contributions. It was found that for skewness values beyond 0.3, the third order cumulant correction becomes numerically unstable and unreliable due to the neglecting of higher order terms; In these cases the numerically stable, albeit inaccurate, second order cumulant approximation should be used instead.