Ultrafast carbon monoxide photolysis and heme spin-crossover in myoglobin via nonadiabatic quantum dynamics

Light absorption of myoglobin triggers diatomic ligand photolysis and a spin crossover transition of iron(II) that initiate protein conformational change. The photolysis and spin crossover reactions happen concurrently on a femtosecond timescale. The microscopic origin of these reactions remains controversial. Here, we apply quantum wavepacket dynamics to elucidate the ultrafast photochemical mechanism for a heme–carbon monoxide (heme–CO) complex. We observe coherent oscillations of the Fe–CO bond distance with a period of 42 fs and an amplitude of ∼1 Å. These nuclear motions induce pronounced geometric reorganization, which makes the CO dissociation irreversible. The reaction is initially dominated by symmetry breaking vibrations inducing an electron transfer from porphyrin to iron. Subsequently, the wavepacket relaxes to the triplet manifold in ∼75 fs and to the quintet manifold in ∼430 fs. Our results highlight the central role of nuclear vibrations at the origin of the ultrafast photodynamics of organometallic complexes.

H emeproteins containing a porphyrin-iron complex (heme) play a major role in the storage and transport of diatomic molecules 1 . The diatomic ligand dissociation occurs concurrently to a spin-crossover (SCO) transition from low-spin (LS) to high-spin (HS) of the Fe II center 2 . The heme complex with carbon monoxide (CO) in myoglobin is one of the most studied hemeproteins of this kind [3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18] . Heme-CO is initially in a singlet LS state (S = 0), which transforms to a quintet HS state (S = 2) upon CO dissociation 10 . The unbound CO in the myoglobin cavity along with the motions of the remaining heme initiate a "protein quake" that opens a channel for the diatomic molecule to escape 3 . Over a century ago, Halden and Lorrain discovered that the CO dissociation can be initiated photochemically 19 , by a mechanism which is still under debate. The photolysis can be initiated upon absorption to the lowest singlet porphyrin band ( 1 Q) or the second porphyrin band ( 1 B) 20 . Here, we focus on the photochemistry of the heme-CO complex upon excitation to the lowest Q-band of porphyrin.
The heme-CO photolysis is an ultrafast process. Recent pumpprobe X-ray experiments of myoglobin, with an initial pump pulse exciting the 1 Q state, seem to agree on a two-step kinetic reaction: (i) a first step taking <50-70 fs, attributed to both CO photolysis and partial SCO, and (ii) a completion of the spin transition to the HS state in~300-400 fs 5,6,20 . Despite the numerous studies on heme-CO photolysis, the kinetics and mechanism of dissociation are still under debate, notably regarding the ultrafast nature of the reaction, and the spin and character of the photolytic state. The rate of photolysis has never been experimentally reported, apart from the upper bound of 50-70 fs. As far as the photolytic state is concerned, the most widely accepted hypothesis is that dissociation occurs from a metal-to-ligand charge-transfer (MLCT) state [15][16][17][18]20 . Still, experiments and theory do not provide a unified picture to date. On the one hand, Franzen, Martin et al., based on time-resolved absorption and Raman experiments, described a "rapid spin state change that must precede photolysis" 20 . They consider that photolysis is occurring from a triplet metal→porphyrin ring transfer ( 3 MLCT dissociation) 20,21 . In the model by Franzen et al., the ultrafast reaction is due to valence tautomerism. This mechanism implies a rapid interconversion of several quasidegenerate electronic states involving d-transitions, some of which are dissociative for the CO bond. On the other hand, Head-Gordon et al., employing time-dependent density-functional theory, proposed a dissociation in a singlet metal → σ Ã CO state ( 1 MLCT dissociation) [15][16][17][18] . In the model of Head-Gordon et al., CO photolysis occurs in a Marcus-like process, in which the initial 1 Q population is transferred to the photolytic state after crossing a barrier of <0.2 eV. However, this hypothesis seems in contradiction with experimental observations of absence of fluorescence emission, an ultrafast reaction and a unity quantum yield for the heme-CO complex 22 . The theoretical method used provided an incomplete description of ultrafast heme-CO photolysis, since the coupled electron-nuclear motion was neglected. Recent theoretical studies of similar iron complexes indicate the fundamental role of nuclear motions in explaining the ultrafast nature of intersystem crossing (ISC) [23][24][25] . Such electronnuclear strong coupling has been experimentally observed in time-resolved X-ray absorption in a Fe(II) tris(2,2′-bipyridine) complex 26,27 and also in the ferricyanide ion 28 .
Here, we resolve the mechanism of ultrafast heme-CO photolysis by means of quantum wavepacket dynamics, which accounts for nuclear and electronic coherent motion 29 . We find that CO photolysis occurs in around 20 fs in the 1 MLCT band, prior to the spin transition. The ultrafast nature of the photolysis is due to strong vibronic couplings and a band of quasicontinuous states. Upon dissociation, the SCO process induces a sequential transition of the remaining heme first into the triplet 3 MLCT and second to the HS quintet manifold 5 MLCT.

Results
Heme-CO model and electronic structure. The model taken as representative of the active center of the myoglobin protein consists of an iron(II) encapsulated in a porphyin ring, and axially ligated to a CO and an imidazole. Hereafter this model is referred to as heme-CO for the complex with carbon monoxide and simply heme for the remaining ligands. The LS state of the heme-CO complex has an electronic configuration (a 2 ) 2 (a 1 ) 2 ðe Ã x ; e Ã y Þ 0 for the porphyrin and (d xy ) 2 (d xz , d yz ) 4 ðd z 2 Þ 0 ðd x 2 Ày 2 Þ 0 for Fe. In its minimum energy structure, the iron and porphyrin are coplanar and the CO is bound to Fe in an upright position perpendicular to the plane. The Fe-C(O) bond is 1.80 Å, around 0.26 Å shorter than the Fe-N bonds with porphyrin or imidazole, which are 2.03 Å and 2.06 Å, respectively. Upon photolysis, the Fe changes into a quintet state, with electronic configuration (d xy ) 1 (d xz , d yz ) 3 d z 2 ð Þ 1 ðd x 2 Ày 2 Þ 1 . In the HS, the Fe-CO bond is unstable. The iron is displaced from the porphyrin plane, distorted to a square-pyramidal structure. The Fe-N bond length is 2.10 Å for the porphyrin nitrogens and 2.17 Å for imidazole.
The full electronic spectrum of heme-CO in the Franck-Condon region is shown in Table 1. This spectrum has been computed in the gas phase, at the CASSCF(10,9)/CASPT2/ANO-RCC-VDZP level of theory for the minimum energy structure of the S 0 ground state (the active space is shown in Supplementary Fig. 1). More details can be found on the effect on the spectrum of the basis set (Supplementary Table 1), the level of theory (Supplementary Table 2) and the model (Supplementary Tables 3-5). At this geometry, the HS state is found 0.7 eV above the LS state. The bright states are the 1 Q states, which correspond to a π → π* transition localized on the porphyrin ring. The 1 Q x and 1 Q y states have an oscillator strength of 1.37⋅10 −2 and 9.13 ⋅ 10 −3 au for Q x and Q y , respectively. In our model, these transitions are found at 2.73-2.79 eV, while experimentally these bands are found at 2.14-2.30 eV 30 . These difference are not unexpected, since on the one hand we do not take into account the electrostatic environment of the protein, and on the other hand the statistical average of several protein conformations. Quantum Mechanics/Molecular Mechanics (QM/MM) calculations decrease the overall spectrum by around 0.84 eV for the lowest part of the spectrum. More importantly, the energetic order and the gap between the states in our model as compared with QM/MM is essentially conserved (Supplementary Table 6). This is essential to guarantee a realistic dynamical treatment. The lowest valence states (up to ca. 2 eV) correspond to dark doubly-degenerate metal-centered (MC) states, involving transitions from the occupied d xy , d xz and d yz orbitals to the unoccupied d z 2 and d x 2 Ày 2 orbitals. The crystal field splitting of the iron orbitals corresponds to a square pyramidal complex. The MLCT band is formed by states of mainly MLCT character, although some MC and LMCT character is also observed (Supplementary Table 7 for a detailed analysis of the state character). These bands of states have strong multiconfigurational character, and were not correctly represented in previous single-reference studies based on density-functional theory 15,18 .
Quantum dynamics. Quantum wavepacket dynamics has been performed using a vibronic model containing the main vibrational coordinates (heme doming, symmetry-breaking, rotational and dissociative modes, see Supplementary Fig. 2 for a plot of the vibrations and Methods section for further details.) The evolution of diabatic state populations during the first 500 fs is shown in Fig. 1a. The results have been obtained by sampling 10 independent quantum dynamics simulations with different initial conditions. The initial conditions are obtained by projecting the heme-CO geometries extracted from molecular dynamics snapshots of a myoglobin protein onto our vibronic model Hamiltonian (Supplementary Table 8.) In the Figure, only the states which are mainly populated are shown ( Supplementary Fig. 3 for the dynamics of all states.) Following the dynamics of each spin manifold, a sequential transfer is found to occur, first from singlet to triplet (S → T) and then from triplet to quintet (T → Q). The time constant for S → T transfer is estimated as 76 ± 15 fs, while the T → Q decay time is 429 ± 70 fs, in perfect agreement with the experimental rates of Cammarata et al. and Franzen et al. 5,6,20 (for the timescales for each initial condition, see Supplementary  Table 9 and Supplementary Note 1.) In the singlet manifold, we observe an abrupt decay of the 1 Q population, which is mainly transferred to the singlet 1 MLCT manifold. The rate for this transfer is estimated to 26 ± 7 fs. The overall reaction is thus sequential 1 Fig. 1b). The initial step is a complete transfer from the 1 Q state to the 1 MLCT state, resulting in a negligible 1 Q population after~100 fs. During the remaining transfer steps the three spin manifolds coexist. This is because the three manifolds are close in energy and strongly mixed through spin-orbit coupling. However, a clearly dominant spin state prevails during different intervals: singlet for t < 75 fs, triplet for 75 < t < 425 fs and quintet at later times.
In Fig. 1c, we show the evolution of the Fe-C(O) and the Fe out-of-plane distances for the first 500 fs (see Supplementary Figs. 4-6 for the corresponding evolution per set of initial conditions). Clearly, we observe damped large amplitude stretching motions for the Fe-C distance, which converges to a quasi-stationary value of 2.2 Å. These motions are coherent, with a period of 42 fs. The Fourier transform of this signal shows a frequency of 800 cm −1 , much faster than the Fe-C(O) stretching frequency of 488 cm −1 . In 21 fs, the Fe-CO distance oscillates between the equilibrium distance (1.7-1.8 Å) and 2.5 Å. These are the typical distances at which the ground-state dissociation occurs 7 . After the first oscillation, the wavepacket is essentially in the MLCT bands, which are repulsive for CO such that the CO does not recombine and undergoes continued oscillations. Relaxation of the structure damps these oscillations until an equilibrium value of 2.2 Å is reached. Both the equilibrium value and the coherent frequency might be disrupted by the presence of the protein, an effect not included in the present model. This would introduce a faster decoherence of the Fe and CO interactions than what we observe in our model. As for the evolution of the Fe out-of-plane motion, we observe oscillations of ±0.2 Å around the initial position, in good agreement with recent time-resolved X-ray crystallography data 4 . As expected, the iron center reacts to the elongation of the Fe-CO distance by following the CO. Very rapidly, it oscillates back to the porphyrin plane, where part of the energy is dissipated to histidine. In the second and subsequent oscillations, the oscillations of Fe and CO are out of phase, indicating that the bond is photolyzed. A large dispersion of the Fe-C distances, which is appearing after 0.2 ps, is due to the delocalization of the wavepacket over the 1,3,5 MLCT bands, indicating an unbound CO molecule. This dispersion is less marked for the Fe out-ofplane distance due to the more localized nature of the Fe atom, bounded strongly to the heme and the proximal histidine.
Experimental evidence shows that the CO is photolyzed within 70 fs 5,6,20 . In order to determine whether dissociation occurs from the singlet or the triplet state, we have considered a reduced model with only the singlet manifold, in which we eliminate the coupling parameters to the triplet and quintet states (Supplementary Figs. 7-9). In this case, an ultrafast transfer from Q → 1 MLCT is still observed. The Fe-C(O) distance again exhibits oscillations with an amplitude of 0.7 Å. The fact that the amplitude is smaller than in the full model can be ascribed to the absence spin-orbit couplings which introduce strong mixing of the 1 MLCT surfaces with 3 MLCT states, that extend the region of oscillations. From these results, we can infer that the transfer to the 1 MLCT manifold is sufficient to dissociate the CO, although the presence of 3 MLCT favors further the photolysis. The dissociation is thus happening between 0.5 and 1.5 periods of the Fe-C(O) oscillation (20-60 fs), when the wavepacket is mainly in the 1 MLCT state.
After photolysis, the SCO mechanism brings the system sequentially to a HS state. After~75 fs, the 3 MLCT band is clearly dominant, although a residual 3 MC participation is also observed ( Supplementary Fig. 1). From the triplet manifold, a population transfer to the HS states is slowly building up, until it becomes dominant at around 400 fs. At this time, the wavepacket is distributed across the 5 MLCT and 5 MC excited states, staying trapped in these states. This is consistent with the slow time constant of~3 ps observed in myoglobin, which has been attributed to a vibrational cooling due to the protein 20,31 . Such excited state trapping has been observed recently in similar iron complexes 27 . This trapping can be understood in terms of the potential energy surfaces depicted in Fig. 2. Upon excitation to the 1 Q band, a sequential ultrafast transfer, first to 1 MLCT and then to 3 MLCT occurs. These transfers are energetically favored and therefore extremely fast. In the 3 MLCT manifold, the wavepacket relaxes in a distribution of 3 MLCT and 3 MC states. The 3 MC state can undergo direct transfer to the 5 MLCT state, which is, however, energetically less favorable (and thus slower) than the singlet → triplet MLCT. The transition from the lowest 5 MLCT to the 5 MC state is energetically unfavoured. On the other hand, the 3 MC states are energetically separated from the rest of the states, such that a transfer to the 5 MC manifold becomes inefficient.
Vibronic effects. The initial structure is close to a C 4v symmetric geometry, in which the CO is in upright position and the porphyrin ring is close to the square planar conformation. The origin of the heme-CO bond is the typical σ-donation π-back-donation mechanism. A partial electron transfer from an occupied σ orbital of CO to an empty d z 2 orbital of Fe occurs simultaneously to a partial back transfer of electron density from the occupied d yz and d xz of Fe to the empty π Ã y;CO and π Ã x;CO orbital, respectively. At this geometry, the π orbitals of porphyrin do not overlap with the dorbitals.
Different types of vibrations can weaken the heme-CO bond, thus liberating the diatomic molecule. In Fig. 3, a schematic representation of the most representative vibrations is shown (for a detailed analysis of all vibrations, see Supplementary Table 10). The most important vibration is the dissociative coordinate (Fig. 3a), which corresponds to a symmetric Fe-CO stretching. This is the main reaction coordinate for CO photolysis. The Fe-CO stretching decreases the overlap between the iron and CO orbitals, thus decreasing both the strength of σ-donation and the π-back-donation. Once the CO is released, the remaining complex is then stabilized in the quintet state. This brings the Fe center out-of-plane, stabilized by π-interactions with the porphyrin antibonding orbitals.
Activation of the dissociative vibration is fundamental for releasing the CO from the heme. This vibration is of a 1 symmetry, as it is a vibration along the principal symmetry axis. However, the initial photon is absorbed exclusively in the porphyrin moiety.  The porphyrin optical Q-state is E-symmetric, thus localized in the perpendicular plane with respect to the principal axis. This implies that the dissociative vibration is not activated in the quasi-C 4v initial geometry. The dissociation of CO thus requires an initial energy transfer from the porphyrin plane to the Fe center, which can only be possible by populating symmetry breaking vibrations. In the excited states, strong symmetry breaking of Jahn-Teller (JT) and pseudo Jahn-Teller (pJT) type occurs 32,33 . The JT effect is the sudden increase of inter-state couplings between an E-degenerate state along symmetry breaking vibrational coordinates. These JT couplings induce a doublewell minimum in the potential energy surfaces at elongated coordinates, while the energy is maximum at the initially symmetric structure. The pJT effect is analogous to the JT effect but applies to non-degenerate states, usually requiring more than a single vibration and more than one electronic state. Both mechanisms are characterized by strong vibronic couplings. In the particular case of the heme-CO complex, the E-symmetry of the Q-band can be broken by b 1 and b 2 vibrations, in the so-called E⊗b JT mechanism 32 . A representative vibration of this type is depicted in Fig. 3d. Such a vibration is not directly coupled to the Fe-CO bond, but is fundamental to induce the interstate coupling between the Q-band and the MLCT band. The relaxation energy induced by such vibrations can be as large as 0.5 eV, which rapidly leads the wavepacket far from the Franck-Condon region. Large amplitude motions have been experimentally observed in similar iron complexes 26 . The rotation of CO (see Fig. 3c) has been discussed at length in the literature, generating a controversy on whether it has an effect on the photochemistry 8,9,16,34 . The rotation of CO slightly weakens the Fe-CO bond due to a less efficient overlap between the densities of each moiety. Furthermore, this rotation induces a Renner-Teller type of coupling (rovibronic coupling), in which e-symmetry bending modes of CO and a 1 -type vibrations couple E and A 1 states, in what is known as the (E⊕A 1 )⊗(e⊕a 1 ) effect. This mechanism can also couple the Q-band (E symmetry) and the MLCT band (A 1 symmetry) to the CO bending (e symmetry) and the dissociative iron-CO stretching (a 1 symmetry). The population evolution of the singlet manifold excluding the JT modes-and thus the rotational modes that are responsible of the Q-MLCT coupling-shows a similar population and Fe-C kinetics. This indicates that any symmetry breaking coordinate, related to protein fluctuations, rotations or symmetry breaking Jahn-Teller effects, will accelerate the coupling to the MLCT band.
The porphyrin doming (a 1 vibration of the porphyrin ring, see modes 8,12 in Supplementary Fig. 3) has been pointed out as an important motion for the photodynamics 5,6,35 . The doming motion is out-of-plane, and has a slow period of 371 fs. It is of the same symmetry as the dissociative mode, to which it is strongly coupled. The potential energy surface cuts along the doming mode indicate a weak reorganization energy and a small vibronic coupling for all potential energy surfaces. Therefore, even though this motion is necessary to stabilize the high-spin state, it plays a passive role and does not significantly affect the short-time dynamics of photolysis (Supplementary Figs. [10][11][12]. Therefore, we conclude that the doming motion is in fact a consequence of the CO dissociation, that is strongly coupled to the Fe-CO dissociation coordinate, and is activated during the second kinetic step of 280 fs 5,6,20 . Spin crossover. The photochemical spin crossover mechanism is an important reaction in many organometallic complexes, in which ligand chromophores allow a usually fast low-to-high spin transition of the complex. Recent theoretical and experimental studies for similar iron complexes have shown that ligands play a fundamental role in modulating the time constants and the electronic pathway during the photoreaction 4-6 . The spin crossover mechanism is mediated by the relativistic spin-orbit coupling (SOC). Although this is a purely electronic transition, the vibrational effects on the SOC seem fundamental to explain the ultrafast nature of spin crossover in metallic complexes. In the particular case of heme-CO, the main SOC strengths are about 200-300 cm −1 , which are of similar size as the JT vibronic couplings. Therefore, internal conversion and intersystem crossing are competing mechanisms in the dynamical relaxation mechanism.
In the heme-CO complex, the initially populated Q-band of porphyrin exhibits a negligible SOC with any other state. In Representation of the nuclear motions (left) and the effect on the Fe-CO orbital interaction (right). The σ CO and d xz /d yz orbitals (in blue) are occupied, while e x/y , π Ã x=y;CO and d z 2 orbitals (in red) are unoccupied. Blue arrows indicate the direction of the partial electron transfer, due to the overlap between an occupied and an unoccupied orbital (solid and dashed arrows indicate strong and weak overlap). Magenta arrows indicate the vibrational motion of the specified atoms. a In the initial structure, the Fe-CO bond is strongest due to partial electron transfer d xz =d yz ! π Ã x =π Ã y (π-back donation) and d z 2 σ CO (σ-donation). b Dissociative vibration, which corresponds to a stretching of the Fe-CO bond. This structure destabilizes the Fe-CO bond and stabilizes by a partial electron transfer d xz =d yz ! e Ã x =e Ã y when the Fe atom is out of plane. c Rotational vibrations, which correspond to the rotation of the CO moiety. This vibration decreases the overlap between Fe and CO orbitals. d Symmetry-breaking vibration, in which the porphyrin is distorted inducing Jahn-Teller and pseudo Jahn-Teller vibronic couplings in the excited state. This vibration has no effect on the Fe-CO bond, but activates the coupling between the porphyrin and the Fe densities NATURE COMMUNICATIONS | DOI: 10.1038/s41467-018-06615-1 ARTICLE NATURE COMMUNICATIONS | (2018) 9:4502 | DOI: 10.1038/s41467-018-06615-1 | www.nature.com/naturecommunications general, SOC is strong when metal density is present in the electronic state. Therefore, the wavepacket is initially dominated by the symmetry breaking vibronic couplings in the singlet manifold, leading to the transition 1 Q → 1 MLCT and 1 Q → 1 MC. Both the 1 MLCT and 1 MC have strong spin-orbit coupling strength with 3 MLCT and 3 MC states, respectively. In general, the SOC is strong when the selection rules ΔS = 1 and ΔL = 1 are obeyed. The latter condition implies that one of the d orbitals changes orientation between the triplet and the singlet state. For example, a strong coupling occurs between the singlet state of character d xy ! d x 2 Ày 2 with the triplet state of character d xz ! d x 2 Ày 2 . We observe that singlet-triplet transitions occur at a faster rate in the MLCT band than within the MC states ( Supplementary Fig. 4). At around 50 fs, the 3 MLCT population becomes dominant, at the expense of the singlet manifold populations. This is because there are a larger number of crossing points between MLCT states of different multiplicity than between MC states (see Fig. 2 and Supplementary Fig. 4). This results in a dominant population of 3 MLCT and 5 MLCT states during the spin-crossover transitions. The 3 Q band triplets, which were pointed out as playing a role in the photophysics of heme-CO, contribute negligibly to the dynamics 20 . The 1 MC state exhibits a residual population after 100 fs, following prevalent transfer to the 3 MC state, with a population which is two times smaller than the 3 MLCT population. Finally, the 3 MLCT → 5 MLCT transfer builds up population in the quintet band. A 3 MC → 5 MC transfer, however, is not observed, due to the fact that the 3 MC and 5 MC states do not cross in an accessible region. Our results show rather that the wavepacket remains trapped in the 5 MLCT band. The trapping occurs in regions of the energy surface where relaxation to lower states (in this case MC states) is slow due to large gaps and small vibronic couplings, whereas a band of electronic states acts as an energy dissipation force, thus localizing the wavepacket in an excited state. Such observation for heme-CO is in line with the valence tautomerism model of Franzen,Martin et al. 20 In summary, the existence of quasi-continuous MLCT bands degenerate to the 1 Q state and strong vibronic couplings are at the origin of the ultrafast photolysis and spin-crossover in the heme-CO complex.

Discussion
The mechanism of heme-CO photolysis and spin crossover has been elucidated by means of high-dimensional quantum dynamical calculations. The results indicate that the vibronic mechanism can be summarized in terms of a sequential 1 Q → 1 MLCT/MC → 3 MLCT → 5 MLCT transfer. During the first step, which is completed within 20-60 fs, the CO is dissociated, prior to the spin-crossover mechanism. The dissociation occurs in the singlet MLCT excited states. For these states, the σ − π bond between Fe and CO bond is weakened. Subsequently, the wavepacket relaxes stepwise to the quintet band under the effect of spin-orbit coupling. We assign the three experimental rates observed for the heme-CO complex in myoglobin, recently reported in Ref. 20 and Refs. 5,6 as follows: (i) the 50-70 fs time constant corresponds to the singlet-triplet spin crossover, in which the population is mainly in the 3 MLCT excited state; (ii) the 300-400 fs time constant corresponds to the triplet-quintet spin crossover, in which the population is mainly in the 5 MLCT excited state; (iii) the 2.4 ps time constant we tentatively assign to the relaxation of the 5 MLCT excited state to the final HS state of the heme.

Methods
Electronic structure. Unconstrained structure optimization and normal-mode calculations of S 0 in gas phase have been performed with density-functional theory (DFT) using the B3LYP exchange-correlation functional and the LAN2DZ basis set [36][37][38][39] . These computations have been performed with Gaussian09 40 . Excited states have been computed with CASSCF(10,9)/CASPT2 41,42 . The active orbitals are the Fe(II) 3d-orbitals and the Gouterman's porphyrin π-orbitals (see Supplementary Fig. 2 and Supplementary Note 2). State-averaged calculations have been performed for 35 singlet states, 50 triplet states, and 30 quintet states. In all calculations, we employed relativistic Douglas-Kroll-Hess (DKH) Hamiltonian with the ANO-RCC-VDZ basis set [43][44][45] . Larger basis sets do not significantly affect the energetics (Supplementary Table 2). Multi-state CASPT2 has been performed with frozen core and 116 frozen virtual orbitals. Similarly to coupled-cluster 46 , deleting virtual orbitals accelerates computations and keeps the essential physics (Supplementary Table 2). In particular, the spin crossover transition is well represented and the DFT minima are also CASPT2 minima. Spin-orbit couplings have been obtained by expanding the DKH Hamiltonian in the basis of 20 singlets, 20 triplets and 20 quintets 45 . All multi-configuration results have been obtained with Molcas 47 .
Vibronic Hamiltonian. The model Hamiltonian for representing the dissociation and the spin crossover is of the form: in which S = 0, 1, 2 is the total spin quantum number, q is the mass-frequencyweighted vibrational coordinate vector,V SO is the spin-orbit coupling. The vibronic Hamiltonian, represented in the basis of ground-state normal modes, is given byĤ containing the harmonic kinetic energy and potential operator. The potential is given byV in which m s is the azimuthal spin and S Ψ m s I is the state wavefunction. The E S I energies are diabatic at the S 0 minimum geometry; at the quasi-C 4v geometry these can be considered quasi-diabatic. The V S IJ ðqÞ potential is derived from onedimensional potential energy surface cuts which were diabatized by fitting to the following functional form, ðq σ Þ are the spin-free adiabatic energies from the model Hamiltonian and the CASPT2 calculation, respectively. The latter are obtained by one-dimensional potential energy surface cuts along each groundstate normal mode direction. The passive normal modes have been kept frozen along the cuts.
The fitting of quasi-continuous bands requires an accurate initial guess. We employ the following protocol: (i) quasi-diabatic surfaces are obtained by maximum overlap of CI coefficients; (ii) pre-fitting is done on the quasi-diabatic surfaces by employing first and second order Taylor expansion or the exponential function, (iii) the adiabatic PT2 surfaces are fitted with the full diagonal Hamiltonian, and (iv) the full Hamiltonian is fitted until the total error is <0.10 eV.
Vibrational modes. Vibrational coordinates have been chosen by performing onedimensional potential energy cuts and selecting the modes with largest reorganization energy. There are a total of 15 explicit vibrational modes in the model: 1 symmetric and 1 asymmetric Fe-CO stretchings, 2 doming modes, 1 CO stretching, 1 in-plane, and 1 out-of-plane porphyrin vibrations, and 4 doubly degenerate modes (Jahn-Teller modes of porphyrin out-of-plane, a porphyrin inplane mode and a pure CO bending.) These modes have been selected in order to reproduce the symmetry breaking that allows the Q → MLCT ultrafast transition, the dissociation and the relaxation of the quintet heme. A schematic plot of these vibrations can be found in Supplementary Fig. 3. Duschinsky rotation effects and anharmonic couplings between vibrational modes have been neglected in our model; in particular, the latter are negligibly small compared to the vibronic couplings (see Supplementary Table 11 and Supplementary Note 3).
Quantum dynamics. The Multi-Layer Multi-Configuration Time-Dependent Hartree (ML-MCTDH) method 48,49 (Heidelberg package, version 8.5.5) was used to perform wavepacket propagation up to 1 picosecond for a vibronic coupling model comprising 179 electronic states of singlet, triplet and quintet multiplicities, and 15 vibrational modes exhibiting strong vibronic coupling. A four-layer representation of the wavefunction was employed, where 8 single-particle functions (spf) were chosen per layer, with a primitive harmonic-oscillator discrete variable representation (DVR) comprising up to 80 DVR points (see Supplementary Note 4 and Supplementary Fig. 13). Initial conditions were constructed in accordance with the relative oscillator strengths of the singlet ligand and metal-centered states at the Franck-Condon geometry.