Figure shows the populations of the ground state S 0 (in red) and of the excited state S 1 (in blue) as functions of time for the simulation starting after the nπ * excitation. The four panels refer to the four flavors of CTTSH implementations presented above. Figure reports the analogous quantities calculated using TSH and TSH-ODC. In Figs. and, the thick lines represent the populations calculated as the fraction of trajectories in a given electronic state m, namely as the ratio F m (t) = N m (t)/N tr , where N m (t) is the number of trajectories for which the active state at time t is m, and N tr is the total number of trajectories; the thin lines represent the electronic population estimated as the average P m (t) over the trajectories of |C α m (t)| 2 , where C α m (t) are the electronic coefficients along the trajectory α. In general F m (t) = P m (t), and the difference is particularly severe if decoherence effects are important and not accounted for, like observed for TSH in Fig. . The decoherence corrections used in TSH-ODC are designed to restore the internal consistency of the TSH procedure in an ad hoc manner, resulting in F m (t) P m (t), as seen in Fig. . Instead, we note a remarkably good agreement across implementations of CTTSH with respect to the electronic populations evaluated as F m (t) or as P m (t). Internal consistency, thus F m (t) P m (t), is correctly achieved with CTTSH when utilizing frozen Gaussians (fixed widths) all along the simulated dynamics, instead, the agreement between F m (t) and P m (t) decreased over long times when thawed Gaussians (variable widths) are employed to calculate quantum momentum. This happens because, as the trajectories become delocalized in space, the widths of the individual Gaussians used for reconstructing the nuclear density increase. Consequently, the quantum momentum, which is evaluated as the spatial derivative of the density, decreases, leading to a reduced effect of decoherence. Nevertheless, we believe that by significantly increasing the number of coupled trajectories, we can attain an accurate representation of the nuclear distribution in configuration space without requiring a substantial increase in Gaussian width. In an ideal scenario, each trajectory would be strongly coupled to several others, making the nuclear density almost insensitive to variations in Gaussian widths within a reasonable range. However, achieving this would necessitate a much larger number of trajectories than the 100-150 (total number of trajectories used in this study) that can currently be accommodated on a single computing node, with the required number increasing with the number of coordinates.