Vibrational wavepacket dynamics in Fe carbene photosensitizer determined with femtosecond X-ray emission and scattering

The non-equilibrium dynamics of electrons and nuclei govern the function of photoactive materials. Disentangling these dynamics remains a critical goal for understanding photoactive materials. Here we investigate the photoinduced dynamics of the [Fe(bmip)2]2+ photosensitizer, where bmip = 2,6-bis(3-methyl-imidazole-1-ylidine)-pyridine, with simultaneous femtosecond-resolution Fe Kα and Kβ X-ray emission spectroscopy (XES) and X-ray solution scattering (XSS). This measurement shows temporal oscillations in the XES and XSS difference signals with the same 278 fs period oscillation. These oscillations originate from an Fe-ligand stretching vibrational wavepacket on a triplet metal-centered (3MC) excited state surface. This 3MC state is populated with a 110 fs time constant by 40% of the excited molecules while the rest relax to a 3MLCT excited state. The sensitivity of the Kα XES to molecular structure results from a 0.7% average Fe-ligand bond length shift between the 1 s and 2p core-ionized states surfaces.


Supplementary Note 2: Correction of Kβ XES spectra measured at the LCLS
Misalignment of one of the four analyzer crystals of the multicrystal spectrometer during the LCLS experiment resulted in distorted Kβ XES spectra. In order to remove this distortion, we express the measured (uncorrected) spectrum S uncorr as a sum of two correct spectra S corr that are shifted with respect to each other: ( ) = (1 − ) ( ) + ( − ) Here, r is the relative intensity of the misaligned (shifted) spectrum and s is the shift. Following from the properties of the δ-function, we can express S uncorr as a convolution ( ) = ( ) * ( ) where S corr is convolved with a function D describing the misalignment: ( ) = (1 − ) ( ) + ( − ) Misalignment correction is therefore equivalent to a deconvolution problem and can be readily solved by applying the convolution theorem. Correct spectrum can be thus calculated from a measured (uncorrected) spectrum with a following formula Direct application of the convolution theorem to carry out the deconvolution is possible here because D is a sum of δ-functions and is noise free. Therefore, singularities do not appear in the division with { }. In general, deconvolution of noisy data requires regularization (e.g. by imposing some smoothness criterium to the result), this is fortunately not necessary here (due to the specific form of D). The correction procedure includes finding appropriate relative intensity r and shift s of the two delta peaks in D. This is achieved by comparison to the known [Fe(bmip) 2 ] 2+ ground state spectrum measured at the SSRL with an identical spectrometer. In Supplementary Figure 3 we demonstrate the application of this procedure to the ground state Kβ XES spectrum. We found that the misaligned portion of the spectrum is shifted by 20 pixels and it corresponds to 28% of the total intensity. Exactly the same D as shown in Supplementary Figure 3 was then applied to correct all the time-resolved XES spectra at different time delays. This is valid because the correction procedure defined above is a linear operation and does not depend on the shape of S uncorr . Therefore, it can be applied to spectra with any shape (i.e. not only the ground state spectrum, but also to all the transient spectra). Also, it is evident from the power spectrum of D, { } * { } in Supplementary  Figure 3 that this procedure does not involve any filtering (e.g. low pass) or smoothing of the data. We demonstrate latter by comparing the residuals of S corr taken with respect to the reference spectrum and the residuals of S uncorr taken with respect to a "distorted" reference spectrum (Supplementary Figure 4).

Supplementary Note 3: Kβ XES singular value decomposition
Singular value decomposition (SVD) of the [Fe(bmip) 2 ] 2+ time-resolved Kβ XES data is displayed in Supplementary Figure 5. It is evident that the first SVD component includes all the time-dependent signal, whereas higher SVD components only include the noise. Therefore, there is no changes in the shape of the time-resolved Kβ XES differences signal. Comparison of the data and the reconstruction based on the first SVD component is shown Supplementary Figure 6.

Supplementary Note 4: Kβ XES oscillatory dynamics
We fitted independently the h delay traces in order to quantitatively evaluate the Kβ XES oscillatory signal. In this analysis the time constants of the Kβ delay trace were not constrained by the Kα data and the MLCT* and the 3 MLCT intensities were allowed to be different. This allows most optimal description of the non-oscillatory Kβ signal and extraction of the oscillatory component. Results of the fits are shown in Supplementary Figure 7. Fit of the Kβ oscillatory signal yields a period of 280±12 fs and a damping constant of 266±110 fs. We find that within the experimental uncertainties, the relative amplitude (~1% signal modulation of the total XES signal) and the period of the oscillations in Kα and Kβ XES are the same. However, there is a discrepancy between the damping constants. We believe that this could be explained by the changing noise level in the Kβ delay trace that seems to be higher after 0.5 ps. Note that overall, the Kβ XES noise level is σ = 0.4%, which is nearly an order-of-magnitude higher than the σ = 0.05% in Kα XES experiment. Although the higher Kβ XES noise level does not allow drawing definite conclusions, we propose that the origin of the Kβ oscillatory signal is the same as for the Kα oscillatory signal. That is, linear shifting of the Kβ XES emission energies due to a small relative displacement of the 1s and 3p core-ionized states PESs.  Figure 8. The data was measured with the same experimental set-up as in Ref. [2]. We find no significant fluence dependent changes in the observed dynamics. Differences in the shape of TA spectra at 0-0.1 ps is likely due to a presence of the cross-phase modulation (CPM) signal. Difference intensities at 0.5-2 ps show linear dependence with the laser fluence (Supplementary Figure  9). A thorough analysis of the TA data is presented in the next section.   Figure  10-S12). This data has been measured with three different laser fluences: 1.5 mJ/cm 2 , 10 mJ/cm 2 and 25 mJ/cm 2 (0.5 kHz, 200 µm pump diameter; same experimental set-up as in Ref. [2]). Close inspection of the 2D TA data reveals that the dynamics is not a single exponential. Fitting of delay traces at different wavelength regions with a biexponential parallel kinetic model revels time constant from 1 to 20 ps. In order to consistently describe these dynamics over the whole detected wavelength range we carried out a global analysis of the TA data. The excited state lifetimes were fixed to the values retrieved from the XES. The 3 MC and 3 MLCT lifetimes are thus 1.5 ps and 9 ps, respectively. Additionally, we included a component corresponding to a hot ground state. Results of this global analysis for the three different laser fluences are shown in Supplementary Figure 10-S12. We find that this model successfully describes the data at all three laser fluences. In Supplementary Figure 13 we display a comparison of a hot ground state difference spectrum extracted from the fit and a difference spectrum of a simulated hot ground state spectrum. Good agreement between these confirms the assignment of the ~20 ps signal to the hot ground state species. We therefore conclude that XES and XSS data measured with 45 mJ/cm 2 laser fluence at the LCLS is consistent with the lower fluence TA experiments. In addition to the excited state dynamics, we find that TA is also sensitive to the temperature of the recovered ground state [Fe(bmip) 2 ] 2+ molecules.

Supplementary Note 7: UV-visible transient absorption oscillatory dynamics
The oscillatory dynamics in [Fe(bmip) 2 ] 2+ was first reported by Liu et al. [2] in UV-visible transient absorption with 485 nm pump (Supplementary Information of Ref. [2]). In Supplementary Figure 14 is a comparison of the UV-visible transient absorption (TA) and the Kα XES kinetics. TA experiment was carried out at 400 nm pump and 520 nm probe wavelengths (same experimental set-up as in Ref. [2]). Both the TA and the Kα XES show oscillatory dynamics. Comparison and fitting of these dynamics is shown in Supplementary  Figure 15. Fit of the TA oscillatory signal yields a period of 280±3 fs and a damping constant of 1480±580 fs. The period is same in both TA and X-ray experiment and can be therefore assigned to the same underlying dynamics. Note that the TA fit shows deviations from the data in delay region from -0.3 ps to 0.4 ps (resulting in overestimation of the fitted damping constant). This is most likely due to a cross phase modulation that is often observed in TA experiments.

Supplementary Note 8: Excitation yield determination
In order to determine the photoexcitation yield of [Fe(bmip) 2 ] 2+ solute molecules during the time-resolved experiment at the LCLS we measured also time-resolved Kβ XES of a reference sample. Latter measurement was carried out immediately after the [Fe(bmip) 2 ] 2+ experiment to guarantee consistent laser excitation conditions. Reference sample was [Fe(btbip) 2 ] 2+ [btbip = 2,6-bis(3-tert-butyl-imidazole-1-ylidine)-pyridine]. This sample was recently investigated with XSS [3]. After a MLCT excitation this complex relaxes ultrafast to a long-lived 5 MC state and this allows it to use for a robust excitation yield determination. Global analysis of the [Fe(btbip) 2 ] 2+ difference Kβ XES kinetics is shown in Supplementary  Figure 16. Comparison of the data and the fit are in Supplementary Figure 17. We find 0.2 ps MLCT lifetime, 279±33 ps 5 MC lifetime and a complete GS recovery. The 5 MC lifetime agrees with the 260 ps lifetime reported in Ref. [3]. In Supplementary Figure 18  Secondly, because we cannot assume a linear excitation regime, we take into account the lowest order non-linear absorption effects. These are: 1) ground state depletion, 2) stimulated emission (SE), and 3) excited state absorption (ESA). SE cross section of 1 MLCT state at 400 nm is equal to the ground state absorption at 400 nm. However, SE is significantly suppressed due to ultrafast vibrational relaxation and 1 MLCT-> 3 MLCT intersystem crossing during the duration of the pump laser pulse (45 fs FWHM). Latter has been measured to be faster than 20 fs in [Fe(bpy) 3 ] 2+ [4]. We therefore introduce 10 fs time constant that depopulates 1 MLCT and populates 3 MLCT. We assume that the 3 MLCT cannot be stimulated to the ground state with 400 nm light. ESA cross section we estimate from the short time-scale UV-visible transient absorption spectra in Supplementary Figure 19. By assuming that the ESA cross section is a slowly varying function the wavelength, we estimate that the ESA cross section is very likely <50% of the ground state cross section. In Supplementary Figure 20 we have simulated the excitation yield dependence from the laser fluence for different ESA cross sections (assumed to be same for 1 MLCT and 3 MLCT). We find that ESA cross section influences the total excitation yield only a small amount and the dominant non-linear effect is ground state depletion that results in a saturated absorption profile. Given the number of assumptions and experimental uncertainties in this estimation, we conclude that the excitation yield of [Fe(bmip) 2 ] 2+ is very likely 84±10%. Therefore, the XES population kinetics analysis was carried out with excitation yields in a range from 74% to 94%. As presented in the main manuscript, variation of the excitation yield within this range does not change the conclusion we draw.

Supplementary Note 9: Comparison of different kinetic models
Different kinetic models were tested to fit the Kα and Kβ XES data. Below we present the results from four kinetic models that were test. The Kα and Kβ XES delay traces were fitted simultaneously and in all the kinetic models the excitation yield was fixed to 84%. All the models succeed in capturing time-dependent dynamics, except the model displayed in Supplementary Figure 21. There the ground state recovery is modeled by a single time constant. It is therefore clear that at least two excited states are required to model the observed decay. A three-state sequential model in Supplementary Figure 22 captures the time-dependent intensity changes. However, it requires the long-lived 3 MLCT state to have too small difference amplitudes, also significantly smaller from the short-lived MLCT* state.
Alternative to a sequential model is presented in Supplementary Figure 23. This model is the same as discussed in the main manuscript and it includes branching of the MLCT* into 3 MC and 3 MLCT that subsequently decay parallel to the ground state. Possible variation of this kinetics is shown in Supplementary Figure 24. There we have tested the likely scenario that the 3 MLCT does not decay to the ground state directly, but via the 3 MC. Because the 3 MLCT lifetime is significantly longer than the 3 MC lifetime, then no considerable 3 MC population builds-up at late delays. We find that both models in Supplementary Figure 23

Supplementary Note 10: Anisotropic XSS signal
Time-resolved XSS difference signal was decomposed into isotropic and anisotropic components following the procedure described in Ref. [5]. The anisotropic XSS signal is shown in Supplementary Figure 25. We found that the anisotropic XSS can be fully described by the optical Kerr-effect (OKE) signal induced in the acetonitrile. The time-dependence of the OKE signal is in an agreement with the dynamics observed in Ref. [6]. The fast initial dynamics <200 fs is slightly non-exponential, however the data can be described with a very good accuracy by two parallel exponential decays with 110 fs and 1.4 ps time constants (in agreement with Ref. [6]). Fit with this kinetics is shown in Fig, S25C. Decay time constants were kept fixed in the fit and we used the OKE signal to check the time-zero (t 0 ) and the width of the Gaussian instrument response function (IRF). The best fit values are t 0 = 20±10 fs and FWHM=110±20 fs.

Supplementary Note 11: Global analysis of XSS
Isotropic XSS data is shown in Supplementary Figure 26. Q-axis was calibrated using a reference acetonitrile heat signal (Supplementary Figure 27A) [7]. Comparison in Supplementary Figure 27A also shows that at time scales >50 ps only solvent heat signal remains. The time-dependence of the solvent heat amplitude was analyzed by fitting the solvent heat signal at each time delay (Supplementary Figure 27B). We found that solvent heating can be well described with two exponential rise times corresponding to 0.35 ps and 13. MLCT*, 3 MLCT and 3 MC populations (from Kα/Kβ XES analysis), oscillatory component (from Kα XES fit), solvent heating and a delay gaussian function that describes the residual ultrafast dynamics that is not captured by the other components. Presence of this residual ultrafast non-exponential component is a result highly non-equilibrium (ballistic) dynamics at early time scales (<200 fs), possibly related to intra-molecular and/or solvation motions within the MLCT manifold. It is also possible that the isotropic XSS signal is contaminated by a small amount with the ultrafast OKE signal due to imperfect anisotropic/isotropic decomposition. However, we emphasize that these uncertainties in the early dynamics do not influence the interpretation of the solute dynamics at later times in any way.
Comparison of the 2D XSS signal reconstructed from the three SVD components and the 2D XSS signal reconstructed from the global fitting shows a very good agreement (Supplementary Figure 30).

Supplementary Note 12: Comparison of different damping mechanisms
Both the solute XSS signal and the Kα XES intensity dependence from the average Fe-ligand bond length R is linear. However, because the wavepacket follows a broad distribution of distances g(R,t), then this can influence the observed signals (e.g. broadening of the Kα XES spectrum in addition to the shift). In order to evaluate these effects we calculated the observed signals for three different g(R,t), but with the same ensemble average <R(t)> (Supplementary Figure 32). We found no effect to the simulated signals depending on whether the damping mechanism is inelastic (cooling, Supplementary Figure 32A-C) or elastic (Supplementary Figure 32D-F). In the third case where g(R,t)=δ(R-<R(t)>), we found only very small differences from the two previous cases at early time delays <0.3 ps (Supplementary Figure  32G-I). Therefore, the observed signals depend only on the <R(t)> and any effects due to the distribution of bond lengths can be omitted. Left column: simulation with inelastic damping mechanism of g(R,t) (vibrational cooling). Middle column: simulation with elastic damping mechanism of g(R,t) (dephasing). Right column: simulation with only the ensemble average Fe-L distances <R(t)>. In all the simulations oscillation period is 278 fs and the exponential damping time constant is 500 fs. Note that g(R,t) also decays with the 1.5 ps lifetime of the 3 MC state.

Supplementary Note 13: Simulations of the XSS signal components
The solute and solvent cage XSS signals of the 3 MC state were simulated based on the calculated DFT structures from Ref. [8] (included in Supplementary Note 19). Supplementary Figure 33A shows the dependence of the solute XSS signals from the average Fe-ligand distance if solvent is not included. Supplementary Figure 33A shows the solvent cage signal retrieved from molecular dynamics (MD) simulations for the GS and the 3 MC state in their optimal geometries (R = 1.943 Å and R = 2.066 Å, respectively). In the MD simulations, the GS and the 3 MC DFT structures of [Fe(bmip) 2 ] 2+ were solvated in a cubic box (50 Å size) of acetonitrile molecules using the three-site interaction potential derived by Guardia et al. [9]. The complex is frozen at the optimized geometry by applying harmonic positional restraints with force constant of 1000 kcal/mol and the parameters for intermolecular interactions are taken from the OPLS 2005 force field [10]. MD trajectories were calculated with a Nose-Hoover thermostat at 300 K for total time of 5 ns, with a time step of 2 fs. Frames were saved every ps and subsequently used to calculate the RDFs, which were sampled with a radial bin of 0.01 Å. For these simulations, we use the Desmond software package developed at D. E. Shaw Research [11].

Supplementary Note 14: RASSCF orbital covalencies
Metal-ligand covalencies of the nominally Fe 3d eg orbital and the nominally ligand occupied σ orbital retrieved from the RASSCF calculations are shown in Supplementary Table 1.
Covalencies are calculated for the optimal GS and 3 MC geometries and for valence and 1s and 2p core-ionized states. The σ-donation is decreased in the 3 MC state with respect to the GS. Also, the covalency is higher in core hole states than in the valence states.

Supplementary Note 15: RASSCF core-ionized PES
RASSCF energies of the core-ionized states (nominally with respect to the optimal 3 MC state energy), together with the core-ionized PESs retrieved from parabolic fits are shown in Supplementary Figure 34. Force constants of all the PESs is very similar: k 1s = 38.0 eV/Å 2 and k 2p = 37.9 eV/Å 2 (latter corresponds to the average 2p PES in Supplementary Figure 34, solid gray line). However, the 1s and 2p PESs are displaces with respect to each other. Minimum energy position of the PESs is respectively R 1s =2.032 Å and R 2p =2.046 Å (latter corresponds again to the average 2p PES). The displacement between the 1s and 2p coreionized PESs is therefore ∆R core =0.014 Å. Because of this displacement, the energy difference between the states changes with a slope κ core =k core ∆R core , where κ core is first-order vibronic coupling and k core ≡k 1s =k 2p .
Supplementary Figure 34. Calculated core-ionized state PESs. Blue crosses are the calculated RASSCF energies of the 1s core-ionized states. Solid black line is a parabolic fit to these energies. Red crosses are the calculated RASSCF energies of a three 2p core-ionized states that contribute most to the intensity of the Kα 1 XES spectrum. Dashed gray lines are parabolic fits to these energies. Solid gray line is the average of the three 2p core-ionized PES.

Supplementary Note 16: Equations of time-dependent signals
Functions used to describe the time-dependence of Kα XES, Kβ XES, XSS and UV-vis TA data are presented below. System of kinetic rate equations that define the excited state populations is: For the clarity of the presentation we have defined here 1 = MLCT*, 2 = 3 MC, 3 = 3 MLCT and 4 = recovered hot GS, 5 = recovered cold GS. Respective rate constants are k 1,2,3,4,5 . r is the branching ratio between 3 MC and 3 MLCT. The boundary conditions are: where N exc corresponds to the initial photoexcited population. Additionally, 1,2,3,4,5 ( < 0) = 0 . Solutions of the above system of equations are convoluted with a Gaussian instrument response function Where t 0 is time-zero and = /�2√2 2� is standard deviation. Thus, excited state populations are ( ) = ( * )( ) (4) , i=1,2,3,4,5. Photoionized population present during the Kα/Kβ XES and XSS experiment is described by � + 1� (5) Where r ion is the ratio of photoexcited population that is ionized. In this case, other populations are rescaled: , i=1,2,3,4,5. The GS population that is not photoexcited is Note that P are defined as relative populations and therefore ∑ ( ) = 1 at any t. Oscillatory signal is described by a damped harmonic oscillation: Here = 2 / , T is the period and τ damp is the damping time constant of the oscillations. Additionally, ( < 0) = 0 and convolution with the MLCT* population is carried out ( ) = ( 1 * )( ) (9) Here 1 is the normalized MLCT* population 1 ( ) = 1 ( )/ ∫ 1 ( ) . Solvent heat signal present in XSS is described by two exponential rise times: where τ i is the exponential rise time and ( < 0) = 0. The residual ultrafast signal present in the isotropic XSS signal at short time scales is described by a Gaussian function: Row vectors V retrieved from the singular-value decomposition of the difference XSS data are model as is the amplitude associated with the time-dependent component i (different for each row vector V j , j=1,2,3). Note that 6 ( ) is excluded because it very small contribution is effectively included in the solvent heat term ℎ ℎ ( ). Row vectors V retrieved from the singular-value decomposition of the UV-vis TA difference data in Section 6 are model as

Supplementary Note 17: Simulation of the 3 MC wavepacket dynamics
Nuclear wavepacket dynamics on the 3 MC surface is simulated with a single-mode quantum harmonic oscillator. Dynamics is calculated using the Lindblad master equation = − ℏ [ , ] + ∑ �2 + − + − + � (16) Here ρ is the density matrix of the system, H is the Hamiltonian of the system and C j is the collapse operator describing interaction of the system with the environment (j=1,2). Hamiltonian is = ℏ ( + 0.5) (17) Where = + is number operator and a is annihilation operator. ω is angular frequency of the oscillator. Oscillator period is = 2 / = 278 fs (ℏ ≈ 14.9 meV). Elastic dephasing is described with a following collapse operator (18) and vibrational cooling is described with (19) where γ 1,2 are dephasing and cooling rates, respectively. Observed exponential time constants are related to the rates by = 2ℏ/ . Initial state is taken to be a coherent Gaussian where R and P are the position and the momentum of the wavepacket and m is the effective mass of the oscillator. In the simulations presented in this manuscript, initial coordinates are R init = 1.943 Å -2.066 Å = -0.123 Å and P init =0. Value of the effective mass was selected to give the initial energy of the harmonic oscillator = 0.5 2 2 = 0.5 . Latter corresponds to a reasonable reorganization energy of the 3 MC state with respect to the 3 MLCT state [12]. For our purpose of calculating the bond length distribution, m is also relevant because the FWHM of coherent Gaussian wavepacket depends on it: Γ = 2.355�ℏ (2 ) ⁄ . In the simulation Γ = 0.025 Å, which is a reasonable spread of the initial wavepacket (note however that Γ has no effect to <R(t)>). Lindblad master equation was numerically solved on a harmonic oscillator basis using the QuTiP package (http://qutip.org) [13,14]. Number of vibrational basis functions was 50. This value was selected to be higher than | | 2 + 3| |, where | | 2 ≈ 33.7 and | | ≈ 5.8 are the mean and the standard deviation of the number of vibrational quanta in the initial state, respectively. Evolution of bond length distribution was subsequently calculated with 0 ( , ) = ∑ 50 , ( ) ( ) * ( ) (20) where ( ) are vibrational basis functions (harmonic oscillator eigenstates). In order to take into account population and depopulation of the 3 MC state, the 0 ( , ) is appropriately convoluted and weighted with the 3 MC decay (defined in the previous Supplementary Note 16): ( , ) is used to describe the experimentally observed oscillatory wavepacket dynamics (shown in Fig. 5B in the main manuscript and Supplementary Figure  32).