3D-SDFs do not indicate strong hydrogen bonds or other structure in the solvation shell around the bipyridyl carbon atoms. Within 100 fs after excitation, the solvation shell starts to deplete around all nitrogen atoms, particularly in the vicinity of the bipyridyl nitrogen atoms. At the equatorial cyanide ligands, this decrease is only observable on the outward-facing side. At the axial cyanides, the oxygen atom occurrence is visibly shifted towards the C-C bridge of the bipyridyl ligand (Fig. ). This trend continues beyond t = 200 fs. Around 500 fs, the oxygen feature at the C-C bridge splits and moves to accumulate above the two bipyridyl para carbon atoms, forming notable hydrogen bonds to the bipyridyl aromatic system that is negatively charged in the MLCT state. For later times until t=5000 fs, the developed features-rings of diminished solvent density around the cyanides and the newly formed hydrogen bonds-intensify at a somewhat slower pace. The solvation shell dynamics as represented by the TD-3D-SDFs provides a much more holistic depiction of the dynamics, compared to the radial distribution functions (RDFs) that are collected in Section S2.3 and Figs. and. Nonetheless, the RDFs agree on the finding of a fast inertial solvent response leading to a reduction in hydrogen bonding to the cyanide groups. A monoexponential fit of the first temporal component of the singular value decomposition of the time-dependent RDF provides time constants of 50-75 fs for the initial solvent response. This fully agrees with a 75 fs time constant obtained from the singular value decomposition of the 3D-SDF, and is also consistent with other computed and experimental time constants for such responses. In the time-dependent RDFs (Figs. and), one can additionally observe several oscillatory features, e.g., in the distribution of the C bpy -H sol , C CN -H sol , or N-H sol distances. Such oscillations are of significant interest, because they can be observed in X-ray solvent scattering experiments and provide insight into relevant coherent vibrational modes. Here, based on the 3D-SDFs, we find that for [Fe(CN) 4 (bipy)] 2-these oscillations arise almost exclusively from intramolecular solute vibrations that modulate the solutesolvent distances, rather than from coherent fluctuations within the solvent. Making this distinction is only possible with the 3D-SDFs, because it is based on the actual three-dimensional coordinates of all solute and solvent atoms, rather than on inter-atomic distances that can only describe relative motion. We note that Section S2.3 also provides a comparison of the RDFs obtained in the present work with the ones obtained in on-the-fly TDDFT/MM simulations. Here, we obtain excellent agreement within the limits of the TDDFT/MM simulations (lower statistics, shorter simulation times), where all extrema and shoulders of the RDFs and their temporal behavior are reproduced in the VC/MM simulations. We conclude from this comparison to the TDDFT/MM data and from the good agreement with experimental solvent response time constants that the overall solvent dynamics of [Fe(CN) 4 (bipy)] 2-is accurately captured in our present simulations.