Asymmetric conformation of the high-spin state of iron(II)-tris(2,2-bipyridine): Time-resolved x-ray absorption and ultraviolet circular dichroism

Spin crossover (SCO) complexes are molecules containing transition metals, most frequently Fe(II) ions, which undergo a ΔS = 2 change of spin from the singlet low-spin (LS, S = 0) ¹A₁ ground state to the lowest high-spin ⁵T₂ (HS, S = 2) quintet excited state under the effect of external stimuli such as temperature, pressure, light and/or engineering of the ligands.^1,2 This makes them attractive for several applications, in particular, magnetic data storage. The discovery of light-induced spin state trapping (LIESST) at cryogenic temperatures in 1984 represented a significant development in this respect,³ with some complexes reaching HS state lifetimes of several days depending on the organic ligand and the temperature. In the near-octahedral ligand field of Fe(II) SCO complexes, the d-orbitals split into lower t_2g bonding and upper antibonding orbitals. In the ground state, all six electrons are in the former while in the HS state, two of these end up in the upper e.g., orbitals leading to an elongation of the Fe-N bond length by about 0.2 Å, which is nearly identical throughout the family of Fe(II)-polypyridine complexes.^4,5

In recent years, Fe(II)-tris-bipyridine ([Fe(bpy)₃]²⁺) has emerged as the prototypical system for the description of light-induced SCO dynamics. Indeed, the relatively high energy difference between the LS and HS states leaves photoexcitation as the only stimulus for SCO. In the visible and ultraviolet, photoexcitation leads to rapid decay with unity quantum yield of the electronic states to the lowest electronic state, which is the metastable high-spin state (HS, S = 2).^1,6–18 Visible excitation of the singlet metal-to-ligand-charge-transfer (¹MLCT) state is more commonly used. Thereafter, decay occurs via the intermediate ³MLCT⁷ and the metal-centered (MC) states¹³ into the HS ⁵T₂ state. Various theoretical interpretations have been proposed for the mechanism and the pathway(s) leading from the LS to the HS state.^19–25 At room temperature and in an aqueous environment, the back-SCO proceeds via non-radiative HS state decay into the LS ground state in ∼0.6 ns.^5,26 This is a relatively fast decay considering the ΔS = 2 character of the transition and it is indeed the shortest for the entire family of Fe(II) polypyridine complexes.

The HS-LS back-SCO in Fe(II)-polypyridine complexes has been rationalized in the framework of the energy gap law^2,27 assuming a single reaction coordinate (the Fe-N bond) characterized by a symmetric displacement of all Fe-N bonds, and a second order spin–orbit coupling matrix element between the HS and LS states. In this picture, the low-temperature relaxation is attributed to tunneling from the HS to the LS states due to their largely different Fe-N equilibrium distances. However, already 40 years ago, theoretical works suggested that in addition to the Fe-N radial expansion, twisting or torsional distortions may also play a role in the back-SCO mechanism of unconstrained tris-chelate complexes, such as [Fe(bpy)₃]²⁺, in the sense that the entropy-driven racemization of enantiopure or enantioenriched solutions at room temperature (RT) is accompanied by a change of spin-state.^28–35 In particular, low-frequency torsional twisting modes associated with racemization result in a surface crossing of the ¹A₁ and ⁵T₂ states, while stabilizing the ³T₁ state and thus increasing the associated spin–orbit coupling term. In recent years, theoretical studies of the SCO have been carried out by density functional theory (DFT) (e.g., see Ref. 36 and references therein) often comparing the iron atom with chemically similar metals (for instance Ru)³⁷ or by time-dependent DFT (TD-DFT).^19–21,^38–40 However, several of these approaches relied on the assumption of a symmetric radial displacement of all N atoms along the Fe-N bonds, which is considered to be accompanied by substantial torsional motion, converting the molecule between an octahedral structure (D₃) and a trigonal-prismatic one (D_3h or C_2v) from the LS to the HS states.³⁴ 

Regarding the geometric structure of the LS state of [Fe(bpy)₃]²⁺, x-ray diffraction⁴¹ and x-ray absorption spectroscopy (XANES and EXAFS) data are available,^5,42 complemented by infrared³⁷ and Raman^43,44 spectra. For the HS state, XANES and EXAFS data are available,^5,42 as well as low-frequency vibrational wave packet dynamics.^{8,12,15,17,18} Very recently, ultrafast transient IR spectra of the HS have also been reported.⁴⁵ However, it should be mentioned that for several Fe(II)-polypyridine SCO systems with somewhat different ligands, there is a rich literature of steady-state and time-resolved IR and Raman spectroscopy whose information content can, to a certain extent, be applied to [Fe(bpy)₃]²⁺.^26,43–48

Ultrafast optical spectroscopy studies showed that the HS state is reached in ≪100 fs, upon excitation of the ¹MLCT state, often yielding a rich pattern of wave packet oscillations in the HS state.^8,9,13,15 This trend (ultrafast SCO and generation of wave packets in the HS state surface) is not specific to [Fe(bpy)₃]²⁺ but was observed in other Fe(II) SCO complexes,^18,49 including non-polypyridine ones.^50–54 Most of these studies reported observation of low-frequency modes, which are very important since they concern the Fe–ligand bonds. However, it is not to be excluded that higher frequency modes could coherently be excited had a higher temporal resolution been used in the above-cited experiments.

Although well-documented theoretically,^{22,28,30–32,35} experimental evidence for the above-mentioned torsional distortions is scarce.³⁵ Ideally, this would be provided by a structural analysis of the HS state. X ray diffraction has largely been used in the case of pressure-, temperature- and light-induced SCO⁴ complexes that can be crystallized and in which the HS state is long-lived,^50,51 but the torsional distortion was not reported, probably due to limited resolution. The non-symmetric HS geometry was attributed by Pàpai et al. to Jahn–Teller (JT) effects arising from the degenerate electronic configuration of the quintet state multiplicity due to the occupation of two e_g^* type antibonding orbitals. Symmetry breaking is not only due to JT distortions;²²  ab initio molecular dynamics studies of an aqueous solution of [Fe(bpy)₃]²⁺ with counterions show that solute–solvent interactions—with water molecules intercalated between the bpy ligands^24,55—might play a key role in modulating the conformation of such an organometallic compound at room temperature.

In a recent article, Oppermann et al. combined ultrafast transient absorption, transient anisotropy and transient circular dichroism (CD) spectroscopy to the study of photoexcited enantiomers of a di-methylated form of the [Fe(bpy)₃]²⁺complex (denoted [Fe(dm-bpy)₃]²⁺).⁵⁶ To our knowledge, this work established, for the first time, an experimental verification of the above-mentioned decade long predictions of a departure of the HS state from the D₃ symmetry. Indeed, the transient CD results show that the HS state adopts the so-called Ray–Dutt twisting coordinate,³⁴ leading to an enhanced spin–orbit coupling of the HS state with the LS ground state, thus facilitating its non-radiative decay and here, we model the structural information on the HS state based on the transient CD data.

For a short-lived HS SCO complexes such as [Fe(bpy)₃]²⁺, time-resolved x-ray⁵⁷ and electron⁵⁸ diffraction have been implemented, but aside from being in crystals, the structural analysis of the HS state assumed a symmetric displacement of the ligands. For solutions, time-resolved x-ray absorption spectroscopy (XAS)^59,60 is the structural tool of choice, as it delivers detailed short-range structural information via the analysis of the x-ray near-edge absorption structure (XANES)⁶¹ and the extended x-ray absorption fine structure (EXAFS).^42,62 In time-resolved XAS, the measured signal is the difference between the absorption of the sample after and before photoexcitation.⁵⁹ 

Here we computationally investigate the molecular structure of the [Fe(bpy)₃]²⁺ complex in the LS and HS states by identifying the stationary points of the geometry landscape through the study of ab initio molecular vibrations by normal mode analysis within the QM-DFT approach in vacuo and in the solvent, using the popular Conductor-like Polarizable Continuum Model (C-PCM) method.⁶³ These methods confirm the already observed trends that while the LS state is indeed a symmetric D₃ structure, the HS state exhibits non-symmetric Fe–N distortions leading to a C₁ symmetry. These conclusions are supported by different experimental observables, namely,: the steady-state and time-resolved UV-visible absorption spectra of [Fe(bpy)₃]²⁺ in water,^8,15 its steady-state and time-resolved XANES spectra,^5,42,64 and the steady-state and time-resolved CD spectra.⁵⁶ Further to this, we calculated the vibrational modes and spectra of the system, which we benchmark in the low-frequency vibrational modes against the vibrational wave packets observed in the HS state,^{8,15,17,18,45} while the high frequency modes are calculated and benchmarked against transient IR data in Ref. 45.

All quantum mechanical (QM) computations were performed in vacuo and with the Conductor-like Polarizable Continuum Model (C-PCM)⁴⁷ method using the Gaussian package.⁶⁵ The work by Daku et al.³⁶ was a seed to investigate the performance of density functional theory (DFT) methods for computing reference molecular geometries for the fitting of the XANES spectra. To this end, considering the electronic states investigated there and starting from x-ray crystallography⁴¹ and the XANES fitted data⁵ the molecular geometry of [Fe(bpy)₃)]²⁺ was extracted and its symmetric D₃ structure preliminarily optimized at the BHLYP/6–311+G^** level for the ¹A₁ singlet spin (Low-Spin – LS) ground state. A comprehensive refinement of the ground state molecular structure computed using the hybrid BHLYP^66,67 was then performed by adopting the same 6–311+G^** basis set using the GGA (generalized gradient approximation) PBE⁶⁸ functional with the inclusion (10%) of the exact (Hartree−Fock) exchange (hereafter termed as PBE^*). This choice was actually driven by results presented by Daku et al.;³⁶ in particular the authors concluded that hybrid functionals tend to give longer metal–ligand bond lengths than the GGAs, and that the more important the exact-exchange contribution the longer the metal–ligand bond. The admixture of 25% exact-exchange to the PBE functional—as in the PBE1PBE formulation—was reported to be in connection with an increase in the Fe–N distance by 0.022 (from 1.971 to 1.993 Å with 6–311+G(d,p) atomic basis set). For this reason, always remaining in the field of the hybrid and GGA functionals for numerical reasons, we have thought to introduce in our work and for the first time, a relatively small fraction (10%) of the exact-exchange in the PBE functional. Such a sort of “contamination” while resembling that already applied in the case of the B3LYP functional (termed as B3LYP^* in the literature)³⁶ allow us to better reproduce the UV/Vis absorption signals in the region between 350 and 450 nm in gas-phase. The same procedure was adopted for the HS ⁵T₂ open-shell excited state to investigate its lowest energy conformation treated in an unrestricted Kohn–Sham (KS) approach, after releasing the constraint of the D₃ molecular symmetry. The calculations were performed using an ultrafine grid of Gaussian, that is 99 590 points, for the integral calculation and to minimize the numerical noise, the “VeryTight” criterion was adopted for all geometry optimization calculations (RMS force threshold set to 10⁻⁶ Hartrees/in bohr). All stationary points characterized by calculations of the first and second derivatives were found to be a local minimum at the best level of method adopted, namely, PBE^*/6–311+G^** while a further energy refinement was carried out at PBE^*/6–311++G(3df,3pd)//PBE^*/6–311+G^** level of calculation. On the optimized geometries, restrained electrostatic potential (RESP)⁶⁹ charges were calculated using the AMBER package for all atoms with optimized ab initio distances as equilibrium values.⁶⁶ 

The UV-Visible (200–700 nm) absorption and circular dichroism (CD) spectra were computed using the time-dependent DFT⁷⁰ and within the continuum polarizable PCM methods by re-optimizing the gas phase molecular geometry in water, dichloromethane (DCM) and trichloromethane (TCM), using the same functional and basis set indicated so far, with a number of transition states from 100 for the LS state to 150 for the HS state to compute their absorbances. In order to compare computed spectra with the experimental ones, calculated oscillator strengths were convoluted with Gaussian functions and computed with the aid of GaussSum⁷¹ code with a full-width half-maximum (FWHM) of 0.40 eV.

The MXAN software procedure employs the full multiple scattering (MS) approach within the muffin-tin approximation for the shape of the potential to calculate the XANES part of the x-ray absorption cross section. The total charge density needed to calculate the whole potential is obtained as the superposition of the charge density of each single atomic species obtained through a self-consistent Dirac–Hartree–Fock procedure. In this way, we can use neutral or non-neutral atoms, once we have obtained information on the net charge in each individual atomic species from other computational programs, as done in this work.^61,72 The MXAN procedure has also the possibility to fit transient spectra, i.e., the difference between the XANES spectra of a laser excited sample minus the unexcited one (LS state). This method reduces the influence of possible systematic errors in the experiment and calculations and at the same time strongly increases the sensitivity of the data to small changes.^61,72

The latest version of MXAN also has the ability to optimize the two nonstructural parameters of the theory, which are the muffin-tin radii and the interstitial potential (see Ref. 61 for details). This is particularly important when comparing fits where the geometric structure has been derived from quantum mechanical calculations, as in this case, or from structures derived from molecular dynamics simulations. In this case, there is no further geometric adjustment of the derived structure as the fit is done only on the two nonstructural parameters. Therefore, this type of fit has always the same number of parameters.

In the following, we will present the results of the calculations and benchmark them against experimental observables for both the ground LS and the HS states.

The PBE^*/6–311+G^** optimized geometry of LS [Fe(bpy)₃)]²⁺ in water dilute solution is shown in Fig. 1, along with the relevant atoms to define its molecular structure parameters. Details about distances and angles are given in Table I where they are compared with the values obtained by the C-PCM molecular geometry optimization and those from x-ray crystallography.⁴¹ The minimum energy optimization of the singlet spin state generates a D₃ molecular geometry at all levels of computations, which was preserved throughout the optimization steps without imposing any geometrical constraints. As a consequence, all Fe–N distances, in plane (N5–Fe–N6, Fe–N5–C2, Fe–N5–C6) and out-of-plane (N3–Fe–N5–C2) angles are equal for all bipyridine ligands. These results agree with the QM/MM simulations of the last ten years,^39,73,74 showing that a DFT approach in the free-molecule approximation, provides a fairly good agreement with the crystallographic structure (Table I), mostly evident on the Fe–N distance but overall, on the value of the Root-Mean-Square Deviation (RMSD) of 0.203 Å over all bond distances of the molecules (gas phase vs x-ray). Optimized D₃ geometries of the isolated complex in the LS state (¹A₁) reported Fe-N distances of 1.997 Å using the hybrid BHLYP functional;⁷⁴ with the same hybrid functional, the same authors reported an average value for the Fe-N distance—computed considering a nanoseconds timescale Car-Parrinello Molecular Dynamics (CPMD) in water dilute solution in a cubic box—of 1.982 Å with a standard deviation (σ) of 0.056 Å. Our simulations by using an all-electron approach report a Fe–N value of 1.966 Å in aqueous environment and of 1.971 Å in gas phase (see Table I). Moreover, the C2-C2a (see Fig. 1) distances of almost 1.466 Å also result in agreement with the statistical distribution extracted from CPMD-based sampling (average value = 1.466 Å, σ = 0.030 Å). In addition, the distribution of the N_(1,3,5)-Fe-N_(2,4,6) angle being centered at 81.5° with a σ = 2.1° well agrees with the detected value of 81.7° at PBE^*/6–311+G^** level of theory for the [Fe(bpy)₃)]²⁺complex.

Using this geometry, we compute the ground state TD-DFT optical absorption spectrum at PBE^*/6–311+G^** level, which is compared in Fig. 2(a) with the experimental absorption spectrum for [Fe(bpy)₃)]²⁺ system at room temperature⁸ and with the one computed with C-PCM method in trichloromethane (TCM) solvent (see Fig. S1 of the supplementary material for comparison with DCM and TCM solvents) after geometry re-optimization. The molecular orbitals involved in the vertical excitation processes for LS multiplicity state are displayed in Fig. 3. Our computed spectra are in line with the experimental data. In fact, the bands at ∼350 and ∼500 nm in the experimental spectrum are at 364 nm (oscillator strength, f = 0.062) and 505 nm (f = 0.066), respectively, in the computed one. Moreover, as expected, these excitations are accompanied by a significant charge transfer (CT) character. The former mainly involves a linear combination of the following single excitations at electronic level: HOMO-2→LUMO + 4, HOMO→LUMO + 7, HOMO-1→LUMO + 4 and HOMO→LUMO + 8, while the latter (known for its ¹MLCT character) shows a multiple excitation character involving—with different weights—the HOMO-2, HOMO-1, LUMO and LUMO + 1 orbitals. Looking at the molecular orbitals in Fig. 3 (upper panel), it is obvious that the unoccupied molecular orbitals reflect a pure valence π^* character. Last, the strongest absorption band experimentally detected around 300 nm in a TCM solution with a shoulder at shorter wavelengths is theoretically estimated as a single band showing a maximum at ∼291 nm. A deeper analysis actually reveals that, this band is mainly composed of absorption features at 288 nm (f = 0.59) and 294 nm (f = 0.17). The former signal is essentially composed of HOMO-5→LUMO, HOMO-4→LUMO + 2 and HOMO-3→LUMO + 1 excitations, while the latter of HOMO-5→LUMO + 2 and HOMO-5→LUMO + 1 single character excitations. In summary, we get a remarkable agreement between experimental and computed spectrum in the 200–700 nm region, not only in terms of band energies but also of relative intensities.

Having benchmarked our computational approach against the LS molecular structure and its UV-Visible absorption spectrum, we calculate the HS state most stable molecular structure using the same PBE^*/6–311+G^** model and computational procedure. From the initial steps of the optimization, the geometry of this state quickly diverges from the D₃ symmetry and ends up in a non-symmetric molecular structure after releasing the major structural parameters around the central Fe²⁺ ion. The resulting optimized geometry parameters are also given in Table I, and Fig. 4 shows the differences between the HS $[ 5 A ( t 2 g 4 e g 2 − like ) ]$ and LS $[ 1 A 1 ( t 2 g 6 e g 0 ) ]$ structures, referred here (and in the following discussion) to the PBE^*/6–311+G^** level of calculation in the gas phase. In the same figure, we report—for the HS state of the [Fe(bpy)₃)]²⁺complex—the Electronic Spin Density by Difference (ESD-D) surface which essentially allow us to derive that the two unpaired electrons characterizing the quintet multiplicity result localized over the central iron (II) atom. For completeness, as for the LS state, the computational pipeline for the HS state has been carried out also in the solvent using C-PCM and for a quick comparison the most relevant (in water solvent) optimized structural parameters are also reported in Table I. The resulting HS structure shows three sets of different Fe–N bond distances with mean values of 2.182, 2.192, and 2.167 Å and slightly distorted in-plane and out-of-plane angles leading to a C₁ symmetry. In addition, we also considered the HS state in a C₂ symmetric sub-space. Such a simulation allows us to identify an optimized structure for such a species featuring interatomic distances resembling those already appreciated by Sousa and coworkers in Ref. 21 (for instance, the Fe–N₁, Fe–N₃ and Fe–N₅ previously estimated in the gas-phase at 2.190, 2.210, and 2.198 Å at PBE0/def2-TZVP level are found lying at 2.171, 2.189, and 2.183 Å in DCM). Furthermore, the C₂ symmetric structure is found at a higher potential energy by 0.19 kcal/mol when compared with the predicted asymmetric counterpart. It is therefore evident that thermal fluctuations at finite-temperature conditions will lead to a population ratio suggesting the coexistence of both symmetric (C₂) and asymmetric HS conformations. Moreover, such an energy gap is found to be in line with the lowest energy bending normal frequencies of the bpy ligands in the HS state multiplicity.

The difference between LS and HS structures can be appreciated by focusing on the most relevant parameters characterizing the iron-bipyridine interaction in Table I, namely, the Fe–N distance and the in-plane and N–Fe–N–C dihedral angle. The latter expresses the extent to which the iron atom is in the plane defined by the aromatic ligands. Evidently, for the LS ground state only one value is reported common to all ligands, because of the D₃ symmetry. In the HS state, the Fe–N interactions are weakened (their distances increase by about 0.2 Å) as a result of the antibonding character of the e_g [SOMO-1(α), SOMO(α)] orbitals depicted in Fig. 3, lower panel. In Table I, we see that the Fe–N distances in ideal gas-phase split into three sets, namely, 2.181–2.183, 2.192–2.193, and 2.166–2.169 Å (with an average distance of 2.181 Å and a standard deviation of 0.011 Å) and correspondingly the bipyridine in plane N-Fe-N angle of 175.2° in D₃ symmetry of the GS splits into three values: 167.2°, 167.8°, and 171.2°, conferring a low symmetry final structure to the HS state. The complex appears as an axially elongated octahedral complex in which the equatorial bonds are, however, equivalent in pairs, both of them (the shorter ones) being located on the same ligand and the other two (the longer ones) on the remaining two ligands. In the water solution, modeled via C-PCM approach, the average value of the Fe–N interatomic distance is equal to 2.174, with a standard deviation of 0.010 Å (see data collected in Table I), i.e., shorter. The observed distortions agree with data available in the literature using pseudopotentials which imply in water solutions changes of 0.01–0.02 Å in the Fe−N bond lengths and of ca. 4° in the N−Fe−N^* bending angle for complex, upon a LS → HS change.⁷⁴ It is also interesting to observe that a connection exists between shorter (and thus stronger) Fe–N bonds and the Fe–N–C–C dihedral angles close to 180°. The nitrogen lone pairs coordinating the central iron atom lie in the planes of the aromatic rings and the dihedral planarity permits to maximize the superposition with the receiving orbitals of the Fe atom. The same considerations obviously hold for a π-back bonding interaction with the π^* bipyridine orbitals.

In Ref. 56, the absorption spectrum of the HS state was extracted for the 260–360 nm region. The C-PCM (in TCM) UV-Vis absorption spectrum of the HS state simulated via TD-DFT, is compared to the experimental one in Fig. 2(b) (see Fig. S2 for a comparison with the gas phase computation). It is estimated for a mixture of C₁(57.9%) +C₂(42.1%) structures, using a population ratio derived from the energy difference of 0.19 kcal/mol at 298 K; moreover, given the slight difference observed between the two structures, the corresponding vertical UV/Vis absorptions are indeed very similar (differences are confined within a few nanometers). The simulated HS spectrum shows a dominant ⁵π→π^* ligand charge (LC) band in the 300 nm region and confirms the absence of intense absorption bands in the >350 nm region in agreement with the experimental spectrum, but also with that of other Fe(II) SCO complexes, whose HS states can be stabilized at low temperatures.² Such a prominent peak mainly arises from three quintet electronic transitions with oscillator strengths of 0.12, 0.19 e 0.22 au. with the following composition: SOMO-9(β)→LUMO(β), SOMO-8(β)→LUMO(β), SOMO-2(β)→LUMO + 4(β), SOMO-2(α)→LUMO + 2(α), SOMO-1(β)→LUMO + 3(β), see the molecular orbitals depicted in Fig. 3, lower panel. We also note that the theory predicts a slightly blue-shifted LC band compared to the experiment.⁵⁶ The remaining low-intensity bands (with oscillator strengths of ∼0.01 a.u.) at longer wavelengths in the 350–600 nm range fall at 362 and 445 nm; these two signals resulting substantially populated by SOMO(β)→LUMO + 8(β), SOMO(β)→LUMO + 4(β), SOMO-5(α)→LUMO + 2(α), and SOMO(α)→LUMO + 2(α) may be then assigned as ⁵MLCT transitions.^2,75 In addition, we would like to highlight that our TD-DFT computations estimate two low-lying quintet excited state at 5705 and 4787 nm resulting from SOMO(β)→LUMO(β) and SOMO(β)→LUMO + 1(β) permutation, respectively (see Fig. 3, lower panel). Such a character in the excitation patterns qualitative agrees with the analysis recently reported by Lee and co-workers with a symmetric counterpart (i.e., ⁵T₂) of the [Fe(bpy)₃)]²⁺ system at DFT level in acetonitrile solution (see Table S11 in Ref. 76).

To summarize this part, both the LS and HS optical absorption spectra are well reproduced by the present computation. The above discussion deals with stabilized structures. However, strictly speaking ensemble averaging would be needed given the floppy nature of the molecule.^22,74 This is especially important for the torsional angles of the Fe(II) complex in solution, for which the distribution of conformations is expected to be broad, as experimentally confirmed through time-resolved UV anisotropy and CD measurements.⁵⁶ From a computational point of view, structural averaging the geometries of LS/HS conformers is a particularly involved task, which is beyond the scope of the present paper. Rather, the stabilized structures obtained above are considered the most probable ones for the respective states, although a thorough consideration of the distribution of conformations would be needed in the future.

This section is dedicated to the comparison between the LS ground state and HS excited state XANES experimental data with the calculations obtained using the stabilized structures obtained in Sec. III A. This is achieved through the MXAN procedure, where the fits are carried out only on nonstructural parameters, as specified in Sec. II B.

We first benchmark our theoretical results on the LS XANES spectrum. The comparison between the experimental data (red dots) and the PBE^* derived geometric structures at the best-fit conditions (see Sec. II B) is shown in Fig. 5. The fit obtained using the PBE^* geometry yields R_sq = 3.18. The agreement between the experimental data and the theoretical calculation is quite good over the whole energy range, with a small discrepancy due to the presence of an extra peak in the energy region about 10 eV after the edge. This good agreement is consistent with the PBE^* Fe–N value being close to that of the crystallographic data,⁴¹ the fit of the XANES spectrum done by moving the geometric structure,⁵ and the fit of the EXAFS spectrum,⁴² all of which yielded identical values (2 ± 0.02 Å).

We now turn to the case of the HS state. In previous analyses,^5,42 the fit of the differential XANES spectrum (excited minus unexcited XANES spectrum at 80 ps time delay) was carried out assuming a D₃ symmetry also in the HS state. The resulting Fe–N bond increase in the HS state was found to be ∼0.2 Å. However, we now see through the PBE^* calculations that a departure from the LS D₃ symmetry occurs in the HS state, with a non-negligible distortion due to inequivalent Fe–N bond elongations (Table I), although the average Fe–N distance remains ∼2.18 Å, in agreement with Refs. 5 and 42.

In order to test the results of the PBE^* structural determination for the HS state, we simulate the differential transient data using the optimized HS geometry following the same procedure used for the LS state that is, at fixed molecular geometry without any further molecular structure refinement to improve the overall fit. Figure 6(a) compares the transient experimental data⁶⁴ (black diamonds) and two calculations made with the PBE^* geometry and with (green dashed line) and without (black solid line) RESP charges. As for the ground state, the use of RESP charges slightly improves the quality of the fit (the error function R_sq decreases from 1.38 to 1.26). In both calculations, the value of R_sq is much lower than 2.97, which is the one found with the global fit (structural in a D₃ symmetry and nonstructural parameters) procedure used in Ref. 5. The value of 2.97 is obtained using the more recent experimental transient data with a higher signal to noise ratio.⁶⁴ In other words, the PBE^* structure reproduces the transient experimental data better than the ones in Ref. 5. To further highlight the improved agreement (with respect to the fit done imposing the D₃ symmetry⁵) between the experimental data and the calculations, using PBE^* nonsymmetrical geometry with and without charges, Fig. 6(b) shows the comparison between the integral of the f_e(E) function for the three fits. The calculations releasing the D₃ symmetry constraint show a significant improvement and among these the one using the PBE^* geometry with RESP charges produces a lower integral in a wide energy range, confirming—as already reported in the literature using pseudopotentials—the preference for the distorted geometry and an asymmetric charge distribution.

In order to complement the above analysis of the XANES, we also carried out simulations of the recently reported steady-state and time-resolved deep-UV CD spectra.⁵⁶ In this work, the molecules were the very similar iron(II)tris(4,4′-dimethyl-2,2′-bipyridine), associated for stereocontrol with the enantiopure Δ or Λ-enantiomer of tris(3,4,5,6-tetrachlorobenzene-1,2-diolato-κ2O1,O2)phosphorus (V) (P(O₂C₆Cl₄)^3– or TRISPHAT) anions.⁷⁸ As seen in Fig. 2, the LS absorption is dominated by the intense LC band around 300 nm. This band is sensitive to the chiral conformation of the molecule, and the CD signal stems from the coupling of the three dipoles that are oriented along the long axis of the bpy ligands.^56,79 Figure 7(a) shows the experimental ground state CD spectrum (in TCM solvent) in the region of the LC band, and Fig. 7(b) shows the simulated one, using the same structural parameters and methods as those used for the UV-Visible absorption and the XANES spectra, optimized for the TCM solvent. Figure S3 shows the simulated LS CD spectra for both water, DCM and TCM solvents. The water spectrum is slightly blue shifted with respect to the other two solvents, reflecting the same trends as in the linear absorption spectrum (Fig. S2). The agreement in Fig. 7 is remarkable, bearing in mind that the deviations that occur have the same origin as those in the UV-Visible absorption spectrum (Fig. 2). In particular, the dissymmetry between the red positive and blue negative lobes of the signal is well reproduced.

We now turn to the case of the HS state. Figure 8(a) shows the experimental CD spectrum of the HS state retrieved from the pump-probe spectrum at 4 ps⁵⁶ in DCM and Fig. 8(b) shows the calculated one in the same solvent, using the same electronic and structural parameters that were used to simulate the UV-Visible absorption (C₁+C₂) and the XANES spectra. Figure S3 compares the simulated HS CD spectra in various solvents. As for the LS state, we also observed a slight blue shift of the water spectrum with respect to the other two solvents. Figure S4 presents the same comparison as in Fig. 8 but using TCM C-PCM modeled calculations. Concerning Figs. 8 and S4, one notes that peak positions of the HS state are well reproduced by the simulation, as are the relative intensities of the two lobes. There are also some deviations, in particular in the experimental spectrum a shoulder appears on the blue positive lobe just above 290 nm, which is not reproduced in the simulations, as was the case with the UV-Visible absorption spectrum (Fig. 2). It may be due to an additional state, not accounted for in our calculations, or to a vibronic band of the LC transition. In summary, our simulated CD spectra agree very well with the experimental ones both in the case of the LS state and the HS state. Most importantly, this agreement is achieved using the same electronic and structural parameters as used for the UV-Visible and the XANES spectra. Therefore, the electronic and structural characteristics of the system in the HS state are supported by two very different types of observables.

One of the goals in the study of the vibrational spectrum of the two electronic states of the [Fe(bpy)₃]²⁺ complex is to characterize the stationary point(s) optimized at the PBE^*/6–311+G^** level as local minimum(a) of the potential energy surfaces (PES). Ultrafast transient absorption spectroscopy have revealed a number of (low-frequency) vibrational wave packets generated in the forward SCO process,^8,15 while transient IR spectroscopy has identified a number of high frequency modes in the 1200–1600 cm⁻¹ region. In the latter, the experimentally derived modes were compared to DFT calculated ones for the ground state, the intermediate ³T state and the ⁵T state.

We present a complete analysis (at the highest level made possible by today's computing resources—DFT with a very extensive basis set) of the vibrational spectrum starting from the calculation of the normal modes of vibration with the study of the Potential Energy Distribution (PED) of both the LS and HS state. In fact, it is thanks to the theoretical study of the vibrational normal modes and the differences between the LS and HS states that we can hypothesize which variations in the internal coordinates between the two states are the most significant in breaking the D₃ symmetry in the HS states. It should also be noted that the use of vibrational normal modes has been a key factor in the interpretation of the time-dependent XANES spectra of Fe(II)-SCO complexes obtained at x-ray Free Electron Lasers.¹⁷ Last but not least, in our opinion this analysis could hopefully pave the way to further theoretical and experimental studies aimed at assessing the key structural features of the strongly coupled vibronic states emerging in the transition from the LS to the HS state in [Fe(bpy)₃]²⁺.

Of the 177 fundamental vibrational modes of the LS D₃ form, 30 are A₁ Raman-active, 29A₂ IR-active and 59 E IR- and Raman-active ones.³⁷ The experimental FTIR and Raman mode frequencies of LS [Fe(bpy)₃]²⁺ were compiled in Ref. 37, based on references therein. The low- to medium-frequency modes here calculated at PBE^*/6–311+G^** level in the gas phase are presented in Table II and compared with the experimental mode frequencies of the LS state. A more detailed PED analysis of the theoretical vibrational spectra of the ligand and of the two electronic states of the complex is provided in Tables S3(a), S4(a), and S5. A further basis for the vibrational assignment of the infrared spectra was provided from the comparison of the theoretical spectra of the complexes with that of the bpy ligand (Table S2). The high frequency (especially above 500 cm⁻¹) vibrational modes are due to the ligand.⁴⁶ The discussion hereafter is restricted to the main differences between the LS and the HS states detected in the calculations, after benchmarking the latter against the low frequency modes observed in the HS state by ultrafast spectroscopy.^8,15 Figure S5 shows the computed IR spectra of LS and HS states. In the following, we discuss separately the low frequency and high frequency regions of the vibrational spectrum.

The low-frequency region: Table S3(a) shows the PED analysis for the LS state vibrations, while Table II compares the calculated vibrational modes for the LS state with the available experimental data from Ref. 37. The agreement is very satisfactory both in the low- and high-frequency regions. The Fe-N stretch mode is around 380 cm⁻¹ for most Fe(II)-polypyridine complexes.⁸⁰ Examining Table S3(a) shows that indeed, a mode at 368 cm⁻¹ dominates this vibration, although other contributions are also present.

For the HS state, as already mentioned, only experimental data in the form of wave packet dynamics in ultrafast pump-probe experiments are available.^8,12,15 The highest temporal resolution was achieved in Ref. 15, which allowed the observation of three modes at 127, 157, and 225 cm⁻¹, with the first being the most prominent and the other two gradually weaker. Table S4(a) shows the PED analysis for the calculated modes of the HS state and Table S4(b) shows the mode frequencies with their assignments.⁸¹ The identification of the modes appearing as wave packets is still a matter of discussion.^15,17 However, the 225 cm⁻¹ mode can safely be attributed to the Fe–N radial stretch vibration as discussed in Ref. 80. In addition, it is nearly identical for all Fe(II)-polypyridine complexes. Table S4(a) (also reported in Table II) indeed shows that the Fe–N stretch is dominated by modes at 257, 253, and 209 cm⁻¹ but also at 121, 114, and 110 cm⁻¹. In fact, the computed harmonic frequencies and normal mode analysis of the group of frequencies at 209, 253, and 257 cm⁻¹ strongly reveal that the asymmetric character of vibrations in this region are all uniquely assignable to Fe–N asymmetric stretching with a minimal contribution of C–N–Fe OOP bending involving the pyridine groups.

Concerning the mode at 157 cm⁻¹, Table S4(b) shows that modes calculated between 190 and 145 cm⁻¹ could all contribute. While a definitive assignment cannot be given, it is important to note that none of the calculated modes contains a significant contribution of Fe–N stretch vibrations but are rather dominated by bending modes and improper angle torsions. In detail, doublet bands at 145/148 cm⁻¹ refer to the bending modes of OOP nitrogen displacements while the doublet at 188/190 cm⁻¹ could safely be assigned to nitrogen in-plane bending modes with a significant displacement of Fe atom. Finally, for the mode at 127 cm⁻¹, the 130 to 110 cm⁻¹ region shows that the Fe–N stretch mode is dominant, although bending and torsional contributions are also present. In fact, the computed harmonic frequencies reveal that normal modes at 110, 114, and 121 cm⁻¹ do actually refer to Fe–N symmetric stretch modes with a partial contribution of out-of-plane (OOP) bending involving C–N–Fe and C–C–C groups. In summary, we can safely attribute the 127 and 225 cm⁻¹ vibrations to predominantly Fe–N stretch modes, while the 157 cm⁻¹ vibration is dominated by bending and torsional modes. These assignments are reported in Table II.

The high frequency region: The modes spanning from 400 to 3000 cm⁻¹ are predominantly due to the organic ligand (Table S2(b)) and the major changes due to the electronic state may be revealed in the CH stretching region and around 1600, 1400, 1300, 700, and below 600 cm⁻¹. This is highlighted in Tables S5 and S6 for both states where the calculated vibrational frequencies (in cm⁻¹), the corresponding IR band intensity (in km mol⁻¹) values and symmetry are provided, and the spectra are shown in Fig. S5. In Table S5, the A₁ symmetry vibrations of the LS state are also included. Although these modes have no IR activity in the D₃ symmetry group, the symmetry breakdown occurring in the HS state is expected to make them IR active and to lift the degeneracy of doubly degenerate modes. Thus, in principle, the vibrational spectrum of the HS state is expected to be more complex than that of the LS state.

The CH stretch region is expected to reflect significant changes of the vibrational structure between the two electronic states. The calculated LS and HS CH stretch frequencies and occasionally, the IR band intensity for the most intense bands is shown in Tables S4 and S5 for the LS and HS states.

From these results, one may easily foresee that several bands of each electronic state are actually expected to overlap. This conclusion is schematically represented in Table S5 reporting the frequencies of the peaks which overlap. Relying on these computations, the most intense CH stretch bands would be produced by the LS state around 3150 cm⁻¹ while the lowest frequency cluster of bands would be evidence of the HS state around 3120 cm⁻¹.

Employing the same approach, let us consider the other bands of the spectra of the two electronic states. As one may note in Fig. S5, the modes in the intermediate region of the HS state for which the three bands centered at 1590 cm⁻¹ and those at 1580 cm⁻¹ produce intense components while the lower intensity band at 1565 cm⁻¹ is the result of the overlap of the three bands marked as (⁠ $* *$⁠) in Table S6. The same reasoning can be employed for the other typical bands, namely, the multiplet at 1400 cm⁻¹ (CC stretching, CH and ring bend modes), the band at 1300 cm⁻¹ (CN and CC stretch modes) and the bands in the 600–800 cm⁻¹ range and those below 600 cm⁻¹, mainly ring modes and CH rocking modes, described in detail in Table S6.

In this way, for the LS ground state it is evident that the 1600 cm⁻¹ closely lying modes would produce a single band whose intensity is largely due to the E symmetry vibrations and the lower intensity band around 1560 cm⁻¹ is due to the overlap between the E and A₂ modes in this range. This region is very different in the HS state where the LS vibrational modes are now replaced by a clear triplet of mode frequencies around 1590, 1580, and 1550 cm⁻¹. The same applies to the bands between 1470 and 1432 cm⁻¹, and the next bands are the three medium-weak intensity bands in the 1300–1200 cm⁻¹ range and in the region of 700 cm⁻¹. The above discussion, along with Fig. S5, provides a recipe for the identification of the distortions in the HS state by means of time-resolved IR or Raman spectroscopy, as recently reported in Ref. 45.

We have investigated the electronic and molecular structure of the low-spin ground state and the high-spin excited state of the prototypical [Fe(bpy)₃]²⁺ molecule, using the PBE^* functional, modeling the solvent interaction via C-PCM method. The calculated molecular structure of the HS state points to a nonequivalent elongation of the six Fe–N bonds in line with previous theoretical predictions,^28,31 and more recent computational reports,^22,34 as well as to a distortion of the bpy ligands, leading to a symmetry lowering of the system from D₃ in the LS state to C₁. The results of the calculations are benchmarked against the available experimental observables, namely, the ultraviolet-visible absorption spectra, the steady-state and transient Fe K-edge XANES spectra and ultraviolet circular dichroism spectra for both the LS and HS states. Bearing in mind that we are dealing with the most probable configuration of a rather floppy molecule, the simulations do provide good agreement with the experiment, correlating these very diverse observables to one set of structural parameters. The observed symmetry reduction along the LS→HS transition is expected to cause the mixing of electronic states and therefore the relatively short lifetime of the HS spin state.

Finally, the stationary points of the potential energy surface were characterized as local minima by calculating harmonic vibration frequencies. The resulting calculated low-frequency vibrational modes are in good agreement with the data from pump-probe experiments and allow us to identify the modes that appear as wave packets in the HS state. However, these are not the modes that cause the asymmetric deformation in the HS. We computed the frequencies and IR absorption spectra from the QM harmonic force field for the LS and HS states, which show substantial differences between them, especially in the mid- to high-frequency region, thus providing a tool for further investigations of the distortions in the HS state using transient IR spectroscopy or time-resolved Raman spectroscopy.

We thank the Department of Chemistry of the University of Rome “Sapienza” and CINECA supercomputing centre for computing resources. This research was supported by the “Departments of Excellence-2018” Program (Dipartimenti di Eccellenza) of the Italian Ministry of Education, University and Research, DIBAF-Department of University of Tuscia, Project “Landscape 4.0 – food, wellbeing and environment”. CZ acknowledges support by the Rome Technopole foundation within the PNRR action in the field of the NextGenerationEU - Sec. 4. “Digital, Industry and Aerospace”. MC acknowledges support of the ERC via the CHIRAX grant (No. 101095012).

The authors have no conflicts to disclose.

Nico Sanna: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal); Supervision (equal); Validation (equal); Writing – original draft (lead); Writing – review & editing (equal). Costantino Zazza: Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal); Validation (equal); Writing – review & editing (equal). Giovanni Chillemi: Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal); Validation (equal). Elisabetta Pace: Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal); Validation (equal). Francesco Cappelluti: Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal). Luigi Bencivenni: Formal analysis (equal); Investigation (equal); Methodology (equal); Resources (equal); Software (equal). Malte Oppermann: Data curation (equal); Validation (equal). Maurizio Benfatto: Conceptualization (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal); Writing – original draft (equal); Writing – review & editing (equal). Majed Chergui: Conceptualization (equal); Formal analysis (equal); Validation (equal); Writing – original draft (equal); Writing – review & editing (equal).

The data that support the findings of this study are available from the corresponding authors upon reasonable request.