Unconventional two-step spin relaxation dynamics of [Re(CO)3(im)(phen)]+ in aqueous solution

Full-dimensional excited-state dynamics simulations including explicit solvation show an unprecedented two-step intersystem crossing mechanism with electronic- and nuclear-driven components in [Re(CO)3(imidazole)(phenanthroline)]+.


S1.1 Initial condition generation
The procedure of sampling the ground state coordinates and velocities is reported elsewhere, 1 and only briefly presented here. Supplementary Fig. S1 gives a graphical overview over the steps involved.
First, classical molecular dynamics simulations were carried out. The force field employed for the complex was taken from the generalized AMBER force field, 6,7 where the missing bonded parameters involving Re were parametrized based on a few-ps QM/MM trajectory. The partial charges were obtained from a RESP fit 8 to electron density obtained at the B3LYP/6-31G* (LANL2DZ for Re 9 ) level of theory, using the Gaussian09 software. 10 The Van-der-Waals radius of Re was set to 1.47 Å. 11 The complex was solvated in a 12 Å truncated octahedron box of 1054 TIP3P 12 H 2 O molecules and a chloride ion was added to neutralize the charge of the complex. After minimization, the system was thermalized to 300 K (NVT ensemble) for 20 ps, followed by equilibration to 1 bar for 200 ps (NPT ensemble). A 10 ns production run (NPT ensemble) was used to obtain 500 equally spaced geometry/velocity snapshots. These steps are represented by the horizontal arrows in Supplementary Fig. S1. All MD simulations employed a time step of 0.5 fs, which was necessary in order to avoid constraining the C-H and N-H bond vibrations of the complex. Note that the O-H bond lengths of the water molecules were also not constrained.
As discussed in Refs. 1,13-15 classical MD at 300 K will bestow too little internal energy to a molecule like [Re(CO) 3 (Im)(Phen)] + , compared to the zero-point energy. In order to provide a more appropriate amount of internal energy, each of the 500 snapshots was heated to 600 K for 20 ps with all solvent molecules frozen (blue arrows in Supplementary Fig. S1). Subsequently, all atoms were imaged into the force field frozen solvent force field QM/MM Excited Figure S1: Preparation of initial conditions. The main steps were classical, force-field-based MD simulations to sample 500 molecular snapshots, local reheating of the solute, QM/MM relaxation in the ground state, and excited-state dynamics simulations. In the excited state, less trajectories were run, because not every initial condition was selected for excitation in the initial state selection scheme.
primary cell, periodic boundary conditions, barostat and thermostat turned off, and the system propagated for 100 fs in the NVE ensemble (black arrows in Supplementary Fig. S1). The resulting endpoints of the 500 trajectories were collected and converted into SHARC initial condition format. During this conversion, we accounted for the fact that AMBER uses a leapfrog-like algorithm (which stores coordinates R(t) and velocities v(t − ∆t 2 )) whereas SHARC uses the velocity Verlet algorithm (coordinates R(τ ) and velocities v(τ )). To compensate this, we computed where ∆t was the time step employed in the AMBER simulations (0.5 fs).
Subsequently, we ran short ground state QM/MM trajectories (level of theory in Supplementary Section S1.2) in order to reduce any effects from imperfect force field parameters. According to Ref. 1, this short relaxation at QM/MM level shifts all bond length and angle distributions much closer to the ab initio optimized values, and shifts the absorption spectrum (see below) by 0.2 eV. In order to avoid that coherent motion originating from the switch from force field to ab initio is carried to the excited state, the short QM/MM trajectories were run for randomized simulation times (between 50 fs and 100 fs).
The end points of these trajectories were then collected and constitute the actual starting points of the excited-state simulations.
Besides initial coordinates and velocities, the initial active state and electronic wave function coefficients are also required to launch the excited-state trajectories. These were obtained for each geometry from a vertical excitation calculation including 30 singlet states at the level of theory given in Supplementary Section S1.2. The absorption spectrum, obtained by summing over all 500 snapshots and convoluting with Gaussians (full width at half maximum (FWHM) of 0.1 eV), is shown in Supplementary   Fig. S2. The most relevant part of the spectrum is located around the experimental excitation energy of [Re(CO) 3 (Im)(Phen)] + , which is 400 nm. 16 Vertical excitation at this wave length can populate the adiabatic states S 1 to S 5 . Whether an initial condition is accepted and which state is actually the initial state for each of the 500 geometries was decided by a stochastic algorithm. 17 Here, in order to receive a sufficient number of initial states from the 500 geometries, we chose an excitation window of 2.8 eV to 3.2 eV, which contains the experimental excitation energy and only the states S 1 to S 5 . Inside this window, the stochastic algorithm accepted 151 initial conditions, with 43 starting in S 1 (28%), 53 in S 2 (35%), 25 in S 3 (17%), 25 in S 4 (17%), and 5 in S 5 (3%). Out of these 151 initial conditions, 100 were actually simulated. From these, 6 were unusable due to network errors that lead to corrupted restart files. Hence, the results presented here and in the main text are based on 94 trajectories, out of which 29 started in S 1 (31%), 32 in S 2 (34%), 17 in S 3 (18%), 14 in S 4 (15%), and 2 in S 5 (2%). As can be seen no initial state was underrepresented in the final set of trajectories. The initial wave function coefficients were simply set to c spin-free i = δ iα , where α is the initial active spin-free state.  Figure S2: Simulated absorption spectrum of [Re(CO) 3 (Im)(Phen)] + in aqueous solution. The total spectrum is given by the thick black line, the contributions of the individual spin-free adiabatic states are shown in color. The excitation energy window from where initial states were drawn is shown as a light grey box. molecules and chloride ion the MM region. Excited states were computed with the Tamm-Dancoff approximation (TDA), and SOCs were computed through a perturbational ZORA formalism. 23 Mixedquality integration grids were used. 1

S1.3 SHARC excited-state dynamics
We considered 7 singlet states (ground state plus 6 excited singlets) and 8 triplet states, i.e., a total of 31 spin-mixed states (7+3×8). In order to construct the gradient of the current active spin-mixed state, we mixed 5 the gradients of all singlet and triplet states that are closer than 0.3 eV to the active state; in this way, we reduce the amount of gradient computations while keeping accurate spin-mixed gradients.
The simulations were run until 250 fs, with a nuclear time step of 0.5 fs. The electronic wave function was propagated with a 0.02 fs time step with the local diabatization algorithm. 24 The necessary wave function overlaps were computed with the WFoverlap code 25 based on auxiliary wave functions generated from the TDA transition density matrices. 26 In order to speed up these overlap calculations, the auxiliary wave function vectors were truncated 25 to 99.97% of the norm, keeping only the largest vector elements.
During surface hops, the velocity vectors of all atoms of [Re(CO) 3 (Im)(Phen)] + were rescaled, but not the velocities of MM atoms. Similarly, in the energy-based decoherence correction 27 instead of the full kinetic energy of the system we only considered the kinetic energy of the [Re(CO) 3 (Im)(Phen)] + atoms.

S1.4 SHARC frozen-nuclei dynamics
The frozen-nuclei dynamics was carried out with the same settings mentioned above, except that the nuclei were replaced by a single dummy atom and the nuclear velocity was set to zero (which also renders the decoherence correction inactive). The quantum-chemical data (Hamiltonian matrix, transition dipole moment matrices) from the first time step of each trajectory was fed to SHARC at all time steps. Gradients were set to zero and the overlap matrix to a unit matrix. In this way, 94 trajectories were propagated until 50 fs.

S1.5 Frozen-nuclei dynamics with explicit laser fields
The frozen-nuclei dynamics including explicit laser fields was done for a single frozen-nuclei trajectory, using three different laser fields. In all cases, the same randomly chosen polarization vector, a phase of 0, and a central energy of 3.3397 eV (corresponding to the energy of a bright spin-mixed state for that geometry) were chosen. The envelope was of sin 2 type: the laser was off for the first 1 fs, then the envelope increased between 1 fs and the pulse center as sin 2 , and then decreased over the same time span, giving a symmetric envelope. The simulations were continued for 5 fs after the pulse finished. The maximum field strengths were chosen to transfer about 75% of the ground state population to the excited states. Note that the simulations only considered the transition dipole moments between S 0 and the excited singlet states; all other static and transition dipole moments were omitted from the dipole moment matrices.

S2.1 Electronic populations and representations
In SHARC, the actual nuclear dynamics is propagated in the so-called "diagonal" representation of the electronic states. 5 The diagonal states are computed as the eigenstates of the full molecular Hamiltonian including spin-orbit couplings (also called "spin-mixed states" in the main text). Because directly computing these eigenstates is very demanding, we apply a two-step procedure. 23,28 First, a number of relevant eigenstates of the (spin-free) molecular Coulomb Hamiltonian (MCH) are obtained as usual with electronic structure code (in this case, ADF), and subsequently the full Hamiltonian matrix in the basis of these MCH eigenstates is constructed (including the MCH eigenenergies and spin-orbit couplings) and diagonalized to produce the diagonal states. A validation of this two-step spin-orbit treatment is presented in Supplementary Section S3.
Due to the two-step procedure, in SHARC the MCH-diagonal transformation matrix is known and it is possible to either analyze the state populations in the diagonal or the MCH representation. The electronic populations in the MCH (i.e., spin-free) representation can be computed in two different ways.
On one hand, one can compute "classical" populations, which requires that the active states in the diagonal representation of all trajectories and all time steps is mapped to MCH states, which can then be counted.
This provides statements like "At time step t, 80% of trajectories are moving in active states that are predominantly S 2 ." However, due to strong mixing in the present case, this mapping is non-trivial. On the other hand, one can compute "quantum" populations by incoherent averaging over the quantum amplitudes. The latter can easily be transformed between diagonal and MCH representation: where i is an MCH state, α is a diagonal state, and U is the MCH-diagonal transformation matrix.
Note that here we employ a decoherence correction that makes the classical and quantum populations approximately equal. We also note that we obtain identical results when computing the electronic populations with the quasi-Wigner protocol of Subotnik and coworkers, 29 showing that the "quantum" populations are appropriate.
We also note that the obtained MCH populations give no indication about the electronic character (i.e., they are not diabatic). The electronic character can be analyzed through the charge transfer analysis (see Supplementary Section S2.4).

S2.2 Densities of states
The densities of states were obtained by collecting all energy differences between the S 0 and the computed singlet and triplet states and convoluting this data with Gaussians (FWHM of 0.05 eV). The active state density was likewise obtained from the energy differences between S 0 and the current active spin-mixed state. For Figure 1f-g, these convoluted densities of state were averaged over the time intervals 8-12 fs ( Figure 1f) and 200-250 fs (Figure 1g), respectively.

S2.3 Time-resolved emission spectrum
The raw (not temporally broadened) emission spectrum was obtained by: where E n and (f osc ) n are the energy difference and oscillator strength for the active→ground state transition and W = (0.25 eV) 2 4 ln 2 defines the width of the Gaussians used to broaden the spectrum on the energy axis (0.25 eV FWHM). To obtain the temporally broadened emission spectrum, the raw spectrum was additionally convoluted in the time domain with a Gaussian: with β = (100 fs) 2 4 ln 2 defining the temporal width (100 fs FWHM). In order to obtain an estimate of the luminescence decay time, the broadened emission spectrum was integrated over all energies and fitted with a bi-exponential decay function f (t) convoluted with G(t).
The bi-exponential fitting function is: Here, A is an overall intensity factor, τ 1 is the fast time constant, τ 2 is the slow time constant, and r is the ratio of the contributions of fast and slow channel.

S2.4 Charge-transfer character analysis
In order to follow the evolution of the electronic wave function, at each time step we computed charge transfer descriptors for all excited states with the TheoDORE program. 30,31 For the analysis, we divided the complex into 9 fragments: (1) Re, (2)-(4) the three CO ligands, (5) imidazole N atoms, (6) imidazole C atoms, (7) phenanthroline N atoms and C atoms in para-position to the N atoms, (8) all other phenanthroline C atoms except for the C=C bridge, and (9) the C=C bridge. H atoms were always included in the fragment of their parent atom. See Supplementary Fig. S3 for a graphical depiction. This fragmentation scheme was chosen based on a correlation analysis as described elsewhere. 31 From the charge transfer analysis, for each of the 94 trajectories for each of the 500 time steps and for each of the 31 states, one obtains a 9 × 9 matrix with elements Ω AB that describes how much contribution of the wave function is given by charge transfer between fragments A and B. From this data, we first compute the charge transfer matrix for the total electronic wave function for that trajectory and step: where Ω i AB is an element of the charge transfer matrix for MCH state i, and the "traj" superscript of U iα , c diag α , and Ω i AB (t) was omitted for clarity. After finding the charge transfer matrices for the total electronic wave function of each trajectory Ω traj AB (t), we average them across all trajectories, leaving one 9×9 matrix for each time step: For a qualitative overview over the electronic character, we sum up multiple matrix elements: Ω phen→phen (t) = A∈phen B∈phen Ω total AB (t), Ω im→phen (t) = A∈im B∈phen Ω total AB (t).
Note that the sum of all elements of Ω total (t) is equal to one. where the full spin-orbit Hamiltonian is diagonalized in the full space spanned by all possible TDA excitations. However, in the SHARC dynamics simulations, we employ a perturbational approach, 23,28 where in a first step the spin-free Hamiltonian is diagonalized in the full excitation space, yielding a number of singlet and triplet states. In the second step, the spin-orbit Hamiltonian is constructed in the subspace spanned by these singlet and triplet states, and diagonalized to yield the spin-orbit energies.
This approach neglects the spin-orbit couplings between the subspace and all higher-lying states, hence the spin-orbit analogue of "dynamical correlation energy" is neglected in the perturbational treatment.
In Supplementary Fig. S4 we show whether the perturbational spin-orbit treatment is sufficient and hence whether the SHARC approach is safely applicable to [Re(CO) 3 Figure S4: Comparison of energies of spin-orbit coupled states using different relativistic approximations. These approximations are: a no SOCs, b perturbational SOC for 6 singlets and 6 triplets, c for 10 singlets and 10 triplets, d for 150 singlets and 150 triplets, and e for full variational 2-component SOC (equivalent to 43,788 singlets and 43,788 triplets).

S4 Supplementary results: Evolution of electronic populations
Supplementary Figure S5 shows the same population data as Figure 1 in the main manuscript, but as a regular line plot. Population This section provides a formal definition of the "spin-orbit wave packet" that is introduced in the main manuscript.
We begin with the definition of the full electronic Hamiltonian: that was already mentioned in passing in Section S3. Importantly, one can either obtain the set of eigenfunctions of the spin-free Hamiltonian: H spin-free |ψ spin-free α = E spin-free α |ψ spin-free α (12) or the set of eigenfunctions of the full Hamiltonian: Note that the eigenfunctions of the latter Hamiltonian do in general not have a well-defined spin multiplicity.
Therefore, it is possible to represent the eigenfunctions of the full Hamiltonian as linear combinations of the spin-free eigenfunctions: Equally, the spin-free eigenfunctions can be represented as linear combination of the full eigenfunctions: Hence, the eigenfunctions of the full Hamiltonian can conveniently be called "spin-mixed". As an example, the states X 1 to X 24 in Figure S4 Unfortunately, the matrix elements ψ spin-free β |Ĥ full |ψ spin-free α do not vanish for α = β, so that the differential equation is not trivially solvable.
Alternatively, one can write the time-dependent electronic wave function in the basis of the eigenfunctions of the full Hamiltonian: In the same way as above, one obtains the equation of motion: but here the Hamiltonian matrix is diagonal, so we can simplify to: with E full j being the eigenenergies of the full Hamiltonian.
The last differential equation can easily be solved to yield the solution of the time-dependent electronic wave function: According to this derivation, the time evolution of the electronic wave function after excitation can best be described as a linear superposition of spin-mixed eigenstates of the full Hamiltonian including SOC. In this basis, only the complex phase of the expansion coefficients depend on time.
The expansion coefficients c full i (0) are determined by projecting the initial wave function |Ψ(0) on the spin-mixed eigenstates: This can be done by inserting the resolution of identity for the spin-free states: Here we assume that the initial electronic wave function is purely singlet-assuming instantaneous, vertical excitation from the singlet ground state. In the simplest such case, c spin-free α = δ ασ , where σ is any single singlet state. Even in this case, multiple c full i (0) will be non-zero. Hence, vertical excitation indeed prepares a coherent linear superposition of the eigenstates of the full Hamiltonian, making the term "spin-orbit wave packet" appropriate.

S6 Supplementary results: Non-Kasha behavior in population flux
The biexponential ISC process is a clear example of non-Kasha behavior. Population transfer from singlet to triplet occurs faster than internal conversion. To illustrate this point, Supplementary Fig. S6 presents the number of net hops between the different electronic states. As can be seen, the number of hops from the S 1 to the triplet states (18 in total) is much smaller than the number of hops from higher singlets to triplets (55 in total). Although the ratio of the net hops shown in Supplementary Fig. S6 are probably not statistically converged (due to including only 94 trajectories), the figure clearly shows the strong non-Kasha behavior of the system. The luminescence decay that is observed experimentally and in the simulated emission spectrum is due to a complex interplay of different effects: (i) ISC from predominantly singlet states to predominantly triplet states, (ii) IC from bright singlet states to less bright singlet states, (iii) modification of brightness of the singlet states by nuclear motion. Additionally, the strong spin-orbit couplings induce strong mixing of the states. Supplementary Fig. S7 summarizes these different effects into one graph. It shows the temporal evolution of the average oscillator strength of each spin-mixed (diagonal) state on the vertical axis. From the vertical axis, one can read how the oscillator strength of the states changes due to the nuclear evolution. Additionally, the figure shows the spin expectation value as color code, and it indicates the most populated states by dot width.
According to the figure, at t = 0 all populated states (large dots) have a relatively large oscillator strength around 0.01. According to the color code, these states are already spin-mixed (spin expectation value around 1), but note that the total electronic wave function at t = 0 is a pure singlet, because the total wave function is a linear combination of the spin-mixed states. Over the course of the simulation, the oscillator strength of some of the strongly populated states decreases (e.g., between 20 and 30 fs from 0.004 to 0.002. More importantly, ISC drives population from the predominantly-singlet states that show high oscillator strengths to the predominantly-triplet states that have lower oscillator strength. Most relevant in this respect is the depopulation of the brightest states (oscillator strength above 0.01), which are largely dominated by S 3 and higher singlet states.