Additionally, we observe that neither the statistical region sampled by the GBOMs, nor the spline fit to Φ(µ (3) , µ (4) ) are symmetric around the origin of µ (3) = 0, µ (4) = 0. Thus, the volatility of the third order cumulant approximation does not depend solely on the magnitude of non-Gaussian behavior (|µ (3) |, |µ (4) |), and specifically the sign of skewness of the energy gap fluctuations imparts a substantial effect on the value of the dampening factor Φ. This finding can be rationalized by examining the relative performance of the undampened (Φ=0) third order cumulant approximation against the second order cumulant approximation. We define the relative metric κ rel = (κ 3 -κ 2 )/κ 2 , where κ 2 and κ 3 correspond to the metric of Eqn. 23 evaluated for the second order and undampened third order cumulant approximation respectively. A plot of κ rel across the range of sampled GBOMs can be found in Fig. . As can be seen, systems with energy gap fluctuation statistics presenting with a negative skew always yield a third order cumulant spectrum that outperforms the second order cumulant approximation. It may be demonstrated analytically (see SI Sec. II. A) for a system with uncoupled GBOM modes (corresponding to the Duschinsky rotation matrix being equal to the identity matrix, [J] = [I]), that a negative skew is indicative of a relaxation in vibrational frequency upon excitation (ω e < ω g ). This causes a divergence proportional to t 2 in the real part of the third order cumulant lineshape function, leading to a dampening of the overall response function and a physical linear absorption spectrum. Conversely, as ω e > ω g a positive skew occurs in the simple model system. As ω e /ω g becomes increasingly large, the lineshape function diverges proportionally to -t 2 , causing an unphysical divergent response function. Additionally, oscillatory divergences with linearly growing amplitudes appear more dominantly in the short timescale relevant to linear optical response and these terms can cause unphysical features in the resulting spectral lineshape. Interestingly, our calculations indicate that an analogous statement holds for GBOMs with couplings between individual modes as described by the Duschinsky rotation: A negative skew is indicative of a stable third order cumulant approximation that improves over the second order cumulant lineshape. For positive skewness values, truncating a cumulant expansion at third order almost universally leads to a deterioration of the spectrum in comparison to the second order approximation.