Atomistic characterization of the active-site solvation dynamics of a model photocatalyst

The interactions between the reactive excited state of molecular photocatalysts and surrounding solvent dictate reaction mechanisms and pathways, but are not readily accessible to conventional optical spectroscopic techniques. Here we report an investigation of the structural and solvation dynamics following excitation of a model photocatalytic molecular system [Ir2(dimen)4]2+, where dimen is para-diisocyanomenthane. The time-dependent structural changes in this model photocatalyst, as well as the changes in the solvation shell structure, have been measured with ultrafast diffuse X-ray scattering and simulated with Born-Oppenheimer Molecular Dynamics. Both methods provide direct access to the solute–solvent pair distribution function, enabling the solvation dynamics around the catalytically active iridium sites to be robustly characterized. Our results provide evidence for the coordination of the iridium atoms by the acetonitrile solvent and demonstrate the viability of using diffuse X-ray scattering at free-electron laser sources for studying the dynamics of photocatalysis.

a, shows the cage term extracted from the difference scattering data at 500 fs and 3 ps time delay. b, shows the solvent cage term simulated from the BOMD at the same time delays, scaled by the 13% excitation fraction. c and d, simulated signal for the two primary contributions to the solvent cage dynamics. c, shows the simulated difference scattering signal when increasing the distance between the Ir2(dimen)4 2+ and an acetontrile molecule placed at 3.6 Å distance along the Ir-Ir axis, and moving the acetonitrile to 4.2 in steps of 0.1 Å. d, shows the simulated difference scattering signal when rotating an acetonitrile molecule placed along the Ir-Ir axis at 3.6 Å distance from 90 o to 30 o with respect to the Ir-Ir axis. Inserts sketch the movements giving rise to the difference scattering signals (gray is carbon, blue is nitrogen, white is hydrogen, pink is Ir and green is the rest of the dimen ligand. The two out-of-plane dimen ligands have been omitted for clarity).

Supplementary Figure 16 | Temporal evolution of the RDFs between Ir and each end of ACN. a, shows
the Ir-Me RDF where r for the peak increases with time, in conjunction with a decreasing g(r) -value. b, shows the Ir-N RDF where the grow-in of a peak at short distances occur after roughly 2 ps. a b

Data Extraction and Reduction
The XDS data was measured in scans consisting of many individual sets of pump/probe events measured at different time delays. Typically a scan would consist of 30.000 individual images, with 600 images at each nominal time delay being recorded over the course of ~5 seconds. To facilitate the substantial data reduction process described below, every 4th pump/probe event was recorded without laser excitation ("laser drop-shots" or "unpumped" images) and every 5th event was recorded without x-rays ("x-ray drop-shots") resulting in every 20 th image with neither the laser excitation pulse nor the x-ray probe pulse arriving at the sample position ("dark images"). Difference scattering images are then constructed by subtracting unpumped images from pumped images. Based on the newly implemented timing tool 2 the data could subsequently be sorted and binned according to actual time delay between x-ray and laser, yielding a data set with a ~130 fs time resolution governed by the combined IRF given by the length of the pump-and probe pulses and the temporal widening due to velocity mismatch in the sample rather than the ~500 fs of the xray/laser jitter observed for this set of experiments The dataflow for the data treatment and reduction is illustrated in Supplementary Figure 3 and is described in detail in the following sections.
Measure XDS: The x-ray diffuse scattering is individually measured for a large number of single X-ray pulses produced by the XFEL at 120Hz. The scattering arising from the interaction between each pulse and the sample is measured on the CSPAD area detector and individually stored. All the diagnostic information about the laser and x-ray pulse along with the corrected pump/probe delay determined by the timing tool is also stored for each pulse to allow sorting and filtering based on the individual shot characteristics and the overall machine/instrument performance on a shot-to-shot basis.
Dark correction: A separate "dark" measurement is conducted where the detector is read out at 120 Hz without incident x-rays for two minutes and averaged to form the "dark image". This allows the electronic background noise of the detector to be determined, and as Supplementary Figure 3 shows, this "dark" contribution to the signal varies significantly across the detector. This dark contribution is observed to be constant over the duration of the measurements, and is subtracted from each image prior to all other corrections. The magnitude of the dark measurement is ~1300 ADU/pixel.

Non-linear corrections:
In the second step of the data reduction applied here, the CSPAD images were corrected for X-ray energy-and intensity non-linearity of the individual pixels as described in detail 3 . These signal variations arising from the non-linearity of the detector responses of the same order of magnitude as the difference signals that are the starting point for the data analysis and are therefore removed in order to facilitate the analysis of the difference scattering signals.
The non-linear corrections were determined for each data acquisition scan based on the unpumped images (see above) in each scan, thus avoiding effects from drifts and changes in pulse characteristics and detector response over the time scale of the entire experiment (5 days).
Calculated corrections: Following the corrections for detector nonlinearity, standard detector corrections were applied to every image. These corrections compensate for variations in solid angle coverage by each pixel, for the x-ray polarization as well as for the x-ray absorption through the slightly tilted liquid sheet 4,5 . As a final step in the correction procedure, the 2D images were corrected by binning the individual pixels as a function of 2θ and adjusting a gain factor for each pixels such that a radially isotropic scattering signal of the corrected unpumped XDS data was enforced.
Quantitative scaling: The individual scattering images were put on an absolute scale of electron units/liquid unit-cell by following the 1D scaling-procedure described in detail in previous work 6 . The measured scattering was scaled to a simulated scattering curve at Q > 3 Å -1 . The simulated scattering curve was calculated from the sum of the coherent and incoherent scattering from an ensemble of molecules representing the stoichiometry of the sample.
Difference scattering : Before calculating the difference signals on which the structural analysis rely, outlier detection and removal was accomplished through a procedure relying on comparing each (azimuthally integrated) image to the median of the rest, a method closely analogous to the one used in our previous work and discussed in detail in the SI of 7 , except with the Chauvenet criterion replaced by a median-comparison . Typically ~10% of the images were rejected in this step.
In order to isolate the signal arising from the structural changes induced in the sample by the pump laser pulse, difference scattering images were constructed by subtracting the average of the nearest two (bracketing) unpumped images from each pumped image. This converts the data from each individual acquisition scan to a (large, ~20000 individual images) set of 2D difference scattering images, each of which is "time stamped" with the actual delay time (~10 fs accuracy) between the laser-and X-ray pulse for that particular difference image.
Temporal binning: To improve the S/N-ratio to a level where quantitative analysis of the difference signals is feasible, the data from 26 acquisition scans was combined and sorted according to time delay. This allows an averaging in "temporal bins", producing averaged difference scattering images in predefined time delay intervals, thus reducing the massive amount of data further significantly increasing the signal to noise ratio of the averaged images.
The data presented here was extracted from a set of 26 consecutive scans comprising ~700.000 individual pump-probe events. To obtain a uniform S/N ratio over the entire time delay range the temporal binning was done for a fixed number of images (200) rather than for a fixed range of time delays. As described above, the data acquisition protocol was designed to ensure almost-uniform temporal coverage in the -2 ps to 4ps region of central interest to these experiments, and the average temporal bin width in this region was approximately 8 femtoseconds 2D data reduction: Due to the short duration of the pulses and short delays studied in these experiments, the difference scattering images contain some anisotropy that will not have had time to rotationally average. Only a subset of solute molecules with a favorable alignment of the transition dipole moment are excited. This preferential excitation results in a population of excited-state molecules where the structural change takes place along the laser polarization . This in turn leads to anisotropic difference scattering signals 8,9 . The anisotropic 2D scattering images can, however, be fully described by two contributions which are a function scattering vector Q only, but which are weighted by the angle with respect to the laser polarization axis. By following the procedure presented in 8 , this allows the robust extraction of the isotropic scattering described by the Debye equation and in turn allows all subsequent analysis to be carried out within the general framework applied in our previous work 6,7,10 .

Data Analysis
The difference scattering signal -S(Q,t) -contains contributions from all the structural changes that occur in the probed sample volume following laser excitation. It is usually expressed as the sum of 3 terms: i) The solute termSsolute(Q,t), arising from structural changes in the solute ii) The solvent term Ssolvent(Q,t), arising from structural changes of the bulk solvent iii) The solute-solvent cross-term (cage term) Scage(Q,t), arising from structural changes in the environment local to the solute, such as in the solvent-shell surrounding the solute.
The solute term -Ssolute(Q,t) The solute term arise from changes in the scattering signal caused by changes in the solute structure. Since two ground state structures have been identified for Ir2(dimen)4 2+ , the solute term needs to account for the depletion of both the depopulated ground states, the structural dynamics of the excited states excited being populated from each of the ground state structures, as well as the fact that the two ground state conformers of Ir2(dimen)4 2+ interconvert on the relevant time scale of the experiment. This can be written as: Where α is the concentration of excited states and β(t) and (1-β(t)) are the relative fractions of depleted ground state structures. SGS1(Q) and SGS2(Q) are the simulated scattering signals from the two ground state structures, SES1(Q,t), SES2(Q,t) are the simulated scattering signals (see below) signal from the excited state population excited from each of the two ground state configurations. In this representation, the time evolution of the difference signal arises from both the excited state structure as well as from the ground state depletion distribution. Thus ) The scattering signal from the ground and excited state structures (SGS1(Q), SGS2(Q), SES1(Q,t) and SES2(Q,t)) was calculated from DFT-derived molecular geometries of ground and excited state Ir2(dimen)4 2+ , by employing the Debye equation: Where k, l and m runs over all atoms in the molecule, fk,l(Q) are the atomic form factors, and rkl is the interatomic distances between every atom pair.
The fitting of measured difference scattering signals utilizing sets of such simulated signals from a range of putative molecular geometries is described in detail elsewhere 6 , and the computational details of the DFT calculations are described in Section 'Supplementary Methods: Computational Details -DFT'. The Ir-Ir distances of the ground state DFT calculations, were constrained to the Ir-Ir distances identified in 10 .
The excited state structures were parameterized in terms of their most significant structural dynamics according to 11 ; namely the Ir-Ir distance (dIr-Ir) and the ligand dihedral twist (DN-Ir-Ir-N), i.e.
For the range of excited state structures used in the analysis, the Ir-Ir distance and ligand dihedral twist were simultaneously constrained in a series of excited-state DFT calculations, such that a matrix of excited state molecular geometries was generated. The Ir-Ir distance was varied between 2.5 Å and 5 Å in steps of 0.5 Å and the ligand dihedral twist was varied between 0 and 50 degrees in steps of 5 degrees. Within the boundaries of this DFT-optimized 'excited state structure matrix' the excited state structure of any combination of Ir-Ir distance and dihedral twist spanned by the matrix could be constructed as a linear interpolation between the nearest structures and from each of these structures the scattering S(Q) could be calculated through the Debye equation introduced above.

The solvent term -Ssolvent(Q,t)
The solvent term arises from structural changes of the bulk solvent structure associated with changes in the hydrodynamic parameters (temperature, pressure, density) of the solvent. Assuming a classical continuum description, the equilibrated state of the solvent can be expressed as a function of two independent hydrodynamical variables, here chosen as the temperature (T) and the density (). The Ssolvent that originates from the bulk-solvent response can be described as a function of their variations ΔT(t) and Δρ(t). Investigations at free-electron laser and synchrotron sources 12,13,14 have demonstrated that for the laser fluence used in these experiments, a first order treatment is adequate to model the response of the solvent from the few picosecond timescale to the hundreds of milliseconds time scale. Within this framework the solvent term is quantified through the following linear combination: are the difference scattering signals arising from a change in temperature, T, at constant density ,, and from a change in density, , at constant temperature, T, respectively. The reference difference scattering signals, were acquired independently during a dedicated study at ID09b, ESRF 12 . It has been shown that for times t such that t < d/vs, where d is the FWHM of the laser spot and vs is the speed of sound in the liquid, no thermal expansion of the solvent has yet taken place 14 . With the present experimental conditions (d = 255 µm FWHM and vs = 1280 m/s for MeCN), thermal expansion is expected to happen on the 200 ns time scale, which is well beyond the temporal window probed in this experiment. Thus, Ssolvent reduces to the contribution from impulsive solvent heating: The validity of ignoring the density response was further verified by trying to include the density term in the fitting procedure. This returned a best-fit value for density change of -0.02 kg/m 3 with an uncertainty of ±0.07 kg/m 3 and we thus conclude any contribution from density changes to be below this 0.001% detection threshold for this parameter.

The cage term -Scage(Q,t)
In the present case of Ir2(dimen)4 2+ the cage term arises predominantly from changes in the structure of the innermost solvent shells of the solute. The Born-Oppenheimer Molecular Dynamics (BOMD) simulations described in the 'Supplementary Methods: Computational Details -BOMD' section shows that two separate processes contribute to the solvation dynamics of [Ir2(dimen)4] 2+ (1) an initial desolvation of acetonitrile methyl groups solvating the axial sites along the Ir-Ir axis of the Ir2(dimen)4 2+ ground state, followed by (2) a coordination of these axial sites of the excited state Ir2(dimen)4 2+ by the nitrogens of the acetonitrile solvent molecules. In the following, the difference scattering signal associated with these two processes is referred to as ) (Q S n desolvatio  and ) (Q S on coordinati  respectively. Thus, the cage term can be expressed as: Where a(t) and b(t) are scaling factors.
As described in 15 , the cage term can be simulated directly from the BOMD trajectories by simulating the scattering of the pairwise radial distribution functions between solute and solvent atoms (ignoring solvent-solvent and solute-solute atom pairs).
This cage scattering was simulated for both the 20 ps ground state trajectory and for the 40 excited state trajectories (running for 3.5 ps each), and by subtracting simulated ground state cage scattering from that of the excited state, the cage term (difference scattering signal) was simulated. As the desolvation and coordination processes are consecutive, the signal from each process was estimated from the BOMD simulations as the difference scattering signal simulated on early (300 fs) and late (3 ps) time delays respectively (see 'Supplementary Methods: Data Analysis -Data Fitting' for more information).

SXDS(Q,t) from the sample
Summarizing the previous sections, S can be expressed as:  Figure 4b) indicate that six components (colored stars) are necessary to describe the data above the noise-level (black stars), which is described by the magnitude of the singular values decreasing in a linear fashion on a semilog plot. The corresponding first six singular vectors are presented in Supplementary Figure 4c and their time evolution is presented in Supplementary Figure 4d. The time-evolution of each of the six singular vectors has significant dynamic around time-zero of the experiment. While there is not necessarily a correspondence between singular vectors and physically meaningful components, the singular value decomposition of the dataset does show that at least six distinct time-resolved processes can be distinguished in the recorded difference scattering data.

Data Fitting
The above equation (Supplementary Equation 1) describes the parameterization of the model used to fit the recorded XDS difference signals. The difference signal for each time delay in the acquired data set is fitted independently by optimizing the following seven fit parameters, which are expressed in terms of six time-dependent parameters: a) Long ground state conformer depletion - These parameters are assumed to describe the full evolution of the data set. However, even though the acquired set of difference signals hold enough information content for a seven-parameter fit to be feasible 6 , it is well known that several of these parameters (In particular the excitation fraction and the Ir-Ir distance in the excited state) can be strongly correlated and that this affects the robustness and reliability of the fit to a significant degree 6 . We therefore implemented the following five-step fit procedure: 1) Identifying the presence of a cage term in the acquired data. 2) Determine the cage term from the data. 3) Fit the data at all time delays with unconstrained parameters. 4) Use the results of step (3) to quantify the excitation fraction and ground state interconversion on the t>2 ps time scale. 5) Fit the dynamics of the remaining parameters on the t<2 ps time scale using the ground state depletion dynamics determined in step (4).
For all of these steps, the fitting was done using the constrained-minimization fmincon routine implemented in Matlab® with a uniform noise distribution (Q) estimated from the standard deviation of the 'laser-off' shots.

Step (1) -Showing the need of a solvent cage term in describing the data
To investigate the potential presence of a solute-solvent cage term in the acquired XDS data, the data recorded for time delays around 5 ps (4.7 ps to 5.3 ps) was fitted using parameters (a) through (e) with (c) and (d) (the excited-state structural parameters) locked to the ~100 ps results previously obtained at a synchrotron source 10 .
Data, fit and residuals for the 5ps delay are presented in Supplementary Figure 5a. The residual matches closely what we reported for the investigation of excited state Ir2(dimen)4 2+ on longer time scales 10 .
Supplementary Figure 5 a, shows the residual of the constrained fit directly compared with the solute-solvent cage term calculated from the BOMD simulations described in the 'Supplementary methods: Computational Details -BOMD' section below. Even though the relative amplitudes of the oscillatory features in the signal vary, good qualitative agreement between the experimental data and the BOMD result is observed. This indicates that a cage response is present in the data and should be accounted for in the analysis.
To ensure that potential bias from possible erroneous structure determination in the earlier synchrotron result is not an issue, another fit was conducted without constraining the excited state structure. Supplementary Figure 6 shows a comparison between the residuals with and without constrained structural parameters and again with a direct comparison with the BOMD cage term. Very little difference between the residual obtained via these two procedures is observed.
Supplementary Figure 7 compares the amplitude of the solute and cage components from the data and the BOMD simulations to ensure that the relative amplitude of the cage component is not mistakenly enhanced in the fit. Supplementary Figure 7a shows the averaged data recorded between 1 and 3 ps (black circles) fitted only with the solute structural components (red curve) and the resulting residual assigned to the cage dynamics (blue curve). Supplementary Figure 7b shows the averaged solute (red) and cage (blue) contribution in the same time frame calculated from the BOMD simulations.
Step (2) -Extraction of the cage signal from the difference scattering data.
As described below in 'Supplementary Methods: Computational Details -BOMD', the BOMD simulations show that the solvation dynamics proceed in two steps, (1) an initial desolvation of acetonitrile methyl groups close to the Ir sites of the Ir2(dimen)4 2+ upon the photoinduced Ir-Ir contraction, followed by (2) an excited state coordination to the Ir sites of the Ir2(dimen)4 2+ by the nitrogen groups of the acetonitrile solvent.
To investigate whether the acquired difference signals contain contributions from such solvation dynamics, the data is fitted with parameters (a) through (e) plus the average of components (f) and (g) calculated from the BOMD simulations. The residual between the data and the fit components In summary the solvation dynamics can be extracted from the difference scattering data, and two components are required to describe them. The experimental cage term agrees well with the simulated cage terms from the BOMD simulations. Comparing the experimental data with a simulation of the difference scattering signal arising from translation and rotation of an acetonitrile molecule aligned along the Ir-Ir axis suggests a larger degree of rotational alignment of the coordinating molecules along the Ir-Ir axis than suggested by the BOMD simulations.
Step (3) -Fitting the difference scattering data set with free parameters.
The data is then fitted utilizing the full 7-parameter model to fit the data. The evolution of the 5 scalar components (components a, b, e, f, g) is shown in Supplementary Figure 10. We here note that both the experimentally determined residual cage terms and the simulated cage terms from the BOMD calculations can be used for the analysis. The structural dynamics derived from the analysis is independent on the choice between residual and BOMD cage terms.
The excited state structural parameters (Ir-Ir distance and ligand dihedral twist), determined for time delays exceeding t=3 ps values are almost identical to the 2.9 Å and 15 degrees of the 100 ps synchrotron measurements 10 . However, the ground state depletion kinetics (red and yellow traces) shows non-physical evolution during the first picoseconds (e.g. the long ground state depletion (red trace) grows during the first few hundred femtoseconds, then decay partially towards the 1 ps time delay, and grows in again slowly over the next 2 ps). We ascribe this non-physical behavior to be caused by a strong correlation between the ground state depletion parameters and the excited state structural parameters. As described in e.g. 10 , such strong correlations between parameters can be removed by determining one of them independently and remove it as a free parameter in the remainder of the analysis. In the present case the ground-state depletion dynamics were estimated as described below.

Step (4) -Determining the ground state interconversion on the <2 ps time scale
The ground state kinetics were estimated and constrained by fitting the evolution of the two ground state depletion amplitudes in the time-range exceeding 2 ps where the difference signal exhibits only slowly-varying dynamics, and then back-extrapolating the evolution the into the 0-2 ps region of primary interest to the present study. In this procedure, linear extrapolation was applied. The ground state depletion kinetics and resulting fit is shown in Supplementary Figure 11 where the ground state depletion kinetics (red and blue trace) are fitted in the t>2 ps region and backextrapolated to time-zero. The black lines in Supplementary Figure 11 show a linear backextrapolation where the ground state depletion dynamics were constrained such that there is a 1:1 ratio between loss of the (preferentially excited) long ground state conformer and the grow-in of the short ground state conformer. The ground state depletion dynamics described by the combined fit was used in the following analysis.
The back-extrapolation of the ground-state depletion dynamics makes it possible to quantify the absolute concentration of states excited from each of the two ground state structures   , and we find that immediately after photo-excitation the excited-state fraction was 13%, 10% excited from the long GS conformer and 3% excited from the short GS conformer. This matches the expectations from optical spectroscopy that we should primarily excite the long ground state conformer at 480 nm 1 . For the structural analysis the fit of the time evolution and the backextrapolation allows us to fix a and b (and thereby ) (t  ) (t  ) in the analysis to the relationship obtained in the analysis presented above.

Step (5) -Fitting the excited state dynamics with the ground state kinetics locked to the results of step (4)
With the analysis steps described above, the full excited state structural and solvation dynamics can now be determined using Supplementary Equation S.eq.1, but restricting the free parameters of the fit to be: , ) (t a , and ) (t b , the dynamics of which are discussed in detail in the main article. The parameters for the species excited from either ground state conformer were constrained such that the relative structural change towards the final value was the same for the two species. This means that at e.g. ~200 fs time delay when the Ir-Ir contraction is 70% of its final value, the Ir-Ir distance of the ensemble excited from the long ground state conformer (4.3 Å) would have contracted to 3.32 Å, while the Ir-Ir distance of the ensemble excited from the short ground state conformer (3.6 Å) would have contracted to 3.11 Å. By independently fitting ) (t a , and ) (t b and constraining the dynamics of the long and short ground states, the degrees of freedom are limited to physically meaningful ranges. This assumes that the absolute number of excited state molecules that are created by the laser pulse and remains constant within the 5 ps. Constraining the relationship between the structural dynamics of the two groundstate structural imposes the assumption that the excited state potential energy surface is harmonic and that the two ground state conformers are excited to the same excited state potential energy surface. This assumption is valid for the excited state dynamics occurring after the initial Ir-Ir contraction, as illustrated by the wavelength independent Ir-Ir stretch vibration 1 . However, the initial ballistic contraction happening on the <200 fs time scale is assumed to occur along an anharmonic potential energy surface. This entails that the analysis results on the 0-200 fs time scale could be smeared out by the contribution from the 20% of the molecules excited from the short ground state conformer, which might be the reason that the initial Ir-Ir stretch oscillation seems weaker in the extracted experimental kinetics than for the simulation. The temperature increase of the bulk solvent is not discussed in the main paper, but is shown in Supplementary Figure 12. The contribution to the signal from the acetonitrile temperature increase is observed to be quite small <0.1 K, and accounts for less than 4% of the total difference signal. This puts this observable at the limit of the resolution of our experiment. This is also seen in the 'free parameter' fit (Supplementary Figure 10) where the solvent heat is responsible for the smallest contribution.
The relatively low heat influx (0.04 K temperature increase at 6 ps time delay) match well with the expected temperature increase from the 13% excitation fraction of a 6 mM Ir2(dimen)4 2+ solution. The upper bound for the excess excitation energy is given by the 0.81 eV energy difference between the 480 nm excitation and 700 nm emission from the S1 state 1 . With a volumetric heat capacity of 63.56 J/(mol *K) for acetonitrile 12 , complete vibrational cooling of all excited S1 states should lead to a temperature increase of 0.050 K The difference signal, fit and fit components at four different time delays are presented in Supplementary Figure 13.

Computational Details -DFT
The scattering from the structures of ground-and excited state Ir2(dimen)4 2+ was simulated via the Debye-expression 15 , utilizing DFT-optimized molecular geometries. These were calculated using the ORCA 2.8 program package 16 using the one parameter hybrid version of the Perdew-Burke-Erzerhoff functional with 25% HF exchange, PBE0 17 , and Ahlrichs type TZVP basis set 18 . The conducting-like screening solvation model (COSMO) 19 was applied with ε = 36.6 D as appropriate dielectric constant for acetonitrile. In addition, effective core potentials (ECP) were used for Ir 20 . Excited state calculations were facilitated by using the lowest lying triplet, since the excited state triplet and first excited singlet potential energy surfaces are known to be very similar 11 . For all calculations the molecular structure was optimized from crystal structures after which Ir-Ir distance and ligand dihedral twist were constrained. For the ground state calculations, Ir-Ir distance and dihedral twist were constrained to the results of our previous study 10 . For the excited state calculations a total of 121 structures were optimized where the Ir-Ir distance and dihedral twist were systematically varied from 2.5 to 5.0 Å in 11 steps and from 0 to 50 degrees in 11 steps respectively.

Computational Details -BOMD
The acetonitrile (ACN) solvent shell response of solvated complex has been investigated using DFT-based Quantum Mechanical / Molecular Mechanics Born-Oppenheimer Molecular Dynamics (BOMD) simulations, as described in 11 , where the results are also compared to previous experimental results, showing good agreement.
Supplementary Figure 14 shows the Ir-solvent pairwise radial distribution functions for 40 ES trajectories binned in 50 fs intervals, and samples the RDF in each time step, resulting in a total of 200 frames used for each ES curve, whereas the GS RDF was sampled over the entire GS trajectory.
In the ground state, the RDF of the ACN methyl groups and the Ir atoms shows a peak around 3.5 Å (Supplenetary Figure 14a). Since this peak has a value of g(r) < 1.0, it is still less likely to find a solvent molecule in this region, compared to the bulk solvent. However, since the probability is still higher in this region than in the adjacent space, this peak represents what can be termed a quasiordering of the methyl ends with respect to the solvent-accessible Ir-regions of Ir2(dimen)4 2+ , parallel to the Ir-Ir axis, There is no analogue to the Ir-Me peak in the Ir-N RDF. For the Cligand RDFs we observe a steeper increase in probability of finding ACN nitrogen atoms, than methyl groups. These observations correspond to the ACN molecules in the Ir2(dimen)4 2+ side-regions having a preferential orientation perpendicular to the Ir-Ir axis, with the N ends oriented towards the ligands, and the Me ends towards the metals. This interpretation is supported by a Bader analysis 21 of a snapshot of the GS trajectory, which assigns roughly half a formal positive charge to each Ir atom, approx. 1 negative formal charge to each of the ligand nitrogen atoms, which are then stabilized by the electro-positive methyl ACN ends.
For the excited state distributions we observe a fast (< 1 ps) decay and ~ 0.7 Å displacement of the gIr-Me(r) peak, and an equally fast grow-in of the peak around r = 9 Å. On longer time scales (> 2 ps), a clear peak in the gIr-N(r) starts to grow in at around 3.5 Å.
The t < 1 ps features are interpreted as arising from the metal atoms contracting, but as seen from the perspective of the solvent, and are thus not related to the solvent cage itself undergoing significant dynamics.
The model of the (changes in) quasi-ordering of the solvent is further underpinned by Supplementary Figure 15. The figure shows the coordination number ratios of Ir-N:Ir-Me (a), and Ligand-N:Ligand-Me (b). A value of 1 is equivalent to a random orientation, since there is an equal amount of nitrogen atoms and methyl groups in ACN. In the GS, the Ir ratio at shorter distances is below 1, meaning that the solvent has a preferred orientation of pointing the nitrogen atoms away from the Ir atoms. This changes significantly in the excited state trajectories, meaning that the preferred orientation inverts after excitation such that the Ir atoms are now preferably solvated by the nitrogen atoms of the acetonitrile.
With time, the ACN-rotation shown in Supplementary Figure 14 settles, resulting in an Ir-N peak growing in around r ~ 3 Å (Supplementary Figure 14). The temporal evolution is visualized in Supplementary Figure 16, which is a 2D plot of cut-outs of the relevant r regions, evolving in time.
In the top plot, the Ir-Ir contraction is again evident. At t > 2 ps, the Ir-N peak starts to grow in upon coordination of the nitrogen groups of the acetonitrile solvent.