Spin-Vibronic Dynamics in Open-Shell Systems beyond the Spin Hamiltonian Formalism

Vibronic coupling has a dramatic influence over a large number of molecular processes, ranging from photochemistry to spin relaxation and electronic transport. The simulation of vibronic coupling with multireference wave function methods has been largely applied to organic compounds, and only early efforts are available for open-shell systems such as transition metal and lanthanide complexes. In this work, we derive a numerical strategy to differentiate the molecular electronic Hamiltonian in the context of multireference ab initio methods and inclusive of spin–orbit coupling effects. We then provide a formulation of open quantum system dynamics able to predict the time evolution of the electron density matrix under the influence of a Markovian phonon bath up to fourth-order perturbation theory. We apply our method to Co(II) and Dy(III) molecular complexes exhibiting long spin relaxation times and successfully validate our strategy against the use of an effective spin Hamiltonian. Our study sheds light on the nature of vibronic coupling, the importance of electronic excited states in spin relaxation, and the need for high-level computational chemistry to quantify it.


INTRODUCTION
The understanding of the interaction among the nuclear degrees of freedom and electronic states, namely the vibronic coupling, is key to understanding the behavior of open-shell coordination compounds and their application in several fields such as photochemistry, [1] photocatalysis [2], ultrafast spectroscopy, [3] and molecular magnetism.[4] For instance, vibronic coupling, also called spin-phonon coupling when restricted to a single spin multiple of states, plays a significant role in the thermally induced relaxation of the magnetic moments of single-molecule magnets (SMMs).The latter exhibit long-lived electronic spin states but relaxation processes due to vibronic coupling hinder their potential applications as quantum bits, memory units, and spintronics elements.[5][6][7][8] To date, a successful description of the thermalization process in SMMs has been achieved using an ab initio theory of open quantum systems.[4,9] This theory allows for the treatment of system-environment interactions on a first principles ground, thus accounting for the coupling between the magnetic ion's spin and molecular vibrations and lattice phonons.Ab initio theory of open quantum systems has already provided important insights into the relaxation dynamics and the resulting loss of spin polarization in SMMs.[9][10][11][12][13][14][15][16][17][18][19] However, the application of this method has been limited to either systems where a mapping of the lowest electronic states into an effective spin Hamiltonian is possible [4,9] , or to systems where the electronic structure is well described by single-reference methods.[20][21][22] Here we aim to remove such restrictions and enable a description of relaxation processes where molecular electronic states are treated explicitly and fully accounting for their multi-reference nature.* lunghia@tcd.ieWhen the spin-Hamiltonian formalism is replaced by a more general approach that uses the full electronic Hamiltonian, the problem can be addressed by adopting the same strategies and concepts employed in the description of excited-state dynamics.[3,[23][24][25] In this context, the investigation of ultrafast relaxation mechanisms revolves around a theoretical exploration of the dynamics of the excited states in the limit where the Born-Oppenheimer (BO) approximation breaks down.[26] A general approach to treat vibronically interacting electronic states would require the use of a diabatic basis where electronic states are no longer parameterized by the nuclear coordinates.Unfortunately, the construction of a diabatic basis for molecular systems is in general not possible [27].Nevertheless, when the vibronic coupling matrix elements change slowly with the nuclear coordinates, an effective method to avoid explicitly constructing a diabatic basis is to expand the vibronic coupling matrix around a chosen reference geometry.The linear vibronic coupling (LVC) model, which involves truncation at the first order, has found widespread application in the investigation of medium to large systems.[28] In this work, we apply ab initio open quantum system theory to the study of angular momentum dynamics of two prototypical SMM complexes, i.e. [Co(C 3 S 5 ) 2 ](Ph 4 P) 2 [29] and [Dy(bbpen)Cl] [30], where H 2 bbpen= N,N'-bis(2-hydroxybenzyl)-N,N'-bis(2methylpyridyl)ethylenediamine).Vibronic coupling matrix elements within the LVC approximation are extracted from the full Hamiltonian of the system and used to evaluate the total relaxation time τ considering both Orbach and Raman relaxation mechanisms.The numerical evaluation of these contributions is discussed and the results are compared with the spin-Hamiltonian approach.We demonstrate that accounting for vibronic coupling among a sufficiently large number of ab initio states can result in an alteration of the computed relaxation times compared to the simple treatment of the lowest angular momentum multiplet.Furthermore, we examine the impact of spin-orbit coupling on the simulated dynamics and analyze its effects, thus providing an unprecedently detailed description of spin-phonon relaxation at the quantum mechanical level.

THEORETICAL METHODS
Vibronic Coupling.We start considering the Molecular Coulomb Hamiltonian (MCH) within the BO approximation which contains the kinetic energy of the electrons (T e ) and the Coulombic interactions electronelectron (V e−e ), electrons-nuclei (V e−n ), and nuclei-nuclei The spin-free MCH wave functions and we assume here that the MCH eigenfunctions are non-degenerate.The MCH Hamiltonian of equation 1 alone is not suitable for describing the magnetic properties of SMM systems, for which it is essential to consider the spin-orbit (SO) interaction that allows the mixing of orbital and spin degrees of freedom.Within the oneelectron effective approximation [31,32], the SO Hamiltonian, ĤSO , is added to the spin-free MCH one as We will refer to H MCH SO as the total Hamiltonian Ĥ matrix elements in the MCH basis of Eq. 2. In general, ĤSO couples the spin-free states |ψ MCH i ⟩ and H MCH SO is nondiagonal in this representation.
By diagonalizing H MCH SO , we obtain the SO-corrected eigenvalues, E i , and the corresponding eigensstates, |ψ i ⟩ The matrix U, which contains as columns the coefficients of the SO eigenfunctions with respect to the MCH basis, is such that where H is the diagonal representation of Ĥ and Within BO approximation, the nuclear degrees of freedom are kept frozen at their equilibrium geometry R 0 and the interaction between electrons and nuclei is purely electrostatic in nature.This assumption prevents a correct description of all those phenomena in which the nuclear motion happens on a short time-scale leading to significant mixing of nuclear and electronic degrees of freedom.Such mixing is responsible for population transfer between different adiabatic states |ψ i ⟩ and it plays a crucial role in describing relaxation processes in molecular systems.[33] We can partly lift the BO approximation considering the LVC expansion in the nuclear displacement of the total Hamiltonian Ĥ where the parametric dependence of Ĥ on the atomic coordinates of the equilibrium geometry R 0 was made explicit.The dependency of Ĥ on the vibrational degrees of freedom is at the origin of the vibronic coupling between the states |ψ i ⟩ and the nuclei.The choice of employing Cartesian coordinates to evaluate the vibronic coupling matrix elements is primarily driven by practical considerations.However, it is always possible to transform these elements into derivatives with respect to normal modes by using where m i is the i-th atomic mass and L αi the Hessian matrix eigenvectors.
Electronic dynamics.The time evolution of the system is determined by the Liouville equation which can be solved within the Born-Markov approximation to obtain the transition rate W ij between two states |ψ i ⟩ and |ψ j ⟩.Here, we describe the nuclei as a phonon bath characterized by normal modes Q α and associated energies ℏω α .Contributions from both one-and two-phonon processes are involved in the relaxation of molecular Kramers systems exhibiting significant magnetic anisotropy.Considering one-phonon processes, the transition rate Ŵ where ℏω ji = E j − E i and the term ∂ Ĥ/∂Q α provides the intensity of the vibronic coupling between the electrons and the α-phonon Q α .The function where nα = (e ℏωα/kBT − 1) −1 is the Bose-Einstein distribution accounting for the phonons' thermal population, k B is the Boltzmann constant, and the Dirac delta functions enforce energy conservation during the absorption and emission of phonons by the spin system, respectively.The Orbach relaxation mechanism, described by Eq. 9, considers the transfer of population through the absorption and emission of a single phonon.An alternative pathway for relaxation towards equilibrium is due to twophonon processes.i.e. the Raman relaxation mechanism.
We model two-phonon spin-phonon transitions, W 2−ph ji , as where the terms ) involve the contribution of all the spin states |ψ k ⟩ at the same time, often referred to as a virtual state.The function G 2−ph fulfills a similar role as G 1−ph for onephonon processes and includes contributions from the Bose-Einstein distribution and imposes energy conservation.G 2−ph accounts for all two-phonon processes, i.e. absorption of two phonons, emission of two phonons or absorption of one phonon and emission of a second one.The latter process is the one that determines the Raman relaxation rate, and in this case G 2−ph reads Once all the matrix elements W n−ph ji have been computed, the relaxation time τ −1 can be predicted by simply diagonalizing W n−ph ji and taking the smallest non-zero eigenvalue.
The study of W 1−ph provides the Orbach contribution to the relaxation rate, τ −1 Orbach , while W 2−ph provides the Raman contribution, τ −1 Raman .The total relaxation time is thus computed as Raman .The theoretical framework just introduced has been successfully applied in the study of SMMs from first principles.[9,19] In these works, lattice harmonic frequencies ω α /2π and normal modes Q α are computed by finite differentiation at density functional theory (DFT) level of theory, while the Hamiltonian Ĥ employed to extract the magnetic properties is the spin Hamiltonian (vide infra).Although this choice is often well justified when studying SMMs, a general approach to the problem must take into account the full ab initio Hamiltonian Ĥ.This choice has the main advantage of preventing any possible loss of information during the construction of the spin Hamiltonian.Moreover, this approach can be applied in situations where it is challenging to determine a priori the specific electronic states that contribute to the relaxation process at a given temperature.
In this paper, our aim is to explore different levels of theory in the evaluation of ∂ Ĥ/∂Q α terms and their effect on the calculated relaxation time.The starting point of our analysis is the evaluation of the Hamiltonian derivative ∇ Ĥ in the ab initio wave function basis, i.e.
∇ Ĥ contains the partial derivative of the Hamiltonian with respect to the 3N Cartesian degrees of freedom, where N is the number of atoms in the system.In the next section the details of the evaluation of ∇ Ĥ are discussed.
Hamiltonian derivative.Vibronic coupling effects are introduced at the first perturbative order in nuclear displacement by considering the gradient of Ĥ, whose matrix elements can be expressed as follows using the Hellmann-Feynman theorem Several theoretical approaches rely on the evaluation of Eq. 15 to unravel the electronic excited-state dynamics beyond the BO approximation.[3] Thus, different electronic structure codes offer the possibility to evaluate forces on atoms and non-adiabatic couplings (NACs) employing both numerical and analytical methods.Nevertheless, the availability of such features is often restricted to specific ab initio methods and rarely with the inclusion of SO effects.Hence, our aim is to build a general framework that allows us to evaluate Eq. 15 numerically, computing the derivative coupling elements ⟨ψ j (R 0 )|∇ψ i (R 0 )⟩ using the wave function overlap between states at different geometries and with the inclusion of SO effects.The need to use the wavefunction overlap to evaluate derivative coupling becomes clear by writing the derivative of ψ i (R 0 ) along R a as central finite differentiation (16) The evaluation of ⟨ψ j (R 0 )|∇ψ i (R 0 )⟩ reduces to the computation of the overlaps ⟨ψ j (R 0 )|ψ i ( R)⟩, where R = R 0 ± ∆R a and a denotes one of the 3N displacement direction used in the numerical differentiation.For this last purpose, it is more practical to use Eq. 6 and work with the MCH wave functions since their composition in terms of Slater determinant coefficients, molecular orbitals (MOs), and atomic orbitals (AOs) can be obtained directly from ab initio quantum chemistry calculations.The derivative coupling between the SO states i and j can be rewritten as where K SO ij (R 0 ) and K U ij (R 0 ), obtained using the product rule of differentiation, represent the derivative couplings between MCH states in the SO basis and the variation of the rotational matrix U respectively, i.e. and It is important to notice that in order to evaluate Eq. 15 special attention must be paid to the calculation of the terms K SO ij (R 0 ) and K U ij (R 0 ) to avoid phase inconsistency between wave functions calculated at different geometries.

Phase correction.
The electronic MCH wavefunctions are obtained by solving the eigenvalue problem 1 parametrized with the nuclear coordinates of the system.Consequently, if |ψ MCH i (R)⟩ is a valid solution to the problem for a given geometry R so is e iϕi(R) |ψ MCH i (R)⟩.The real number ϕ i (R) ultimately depends on the specific implementation of the eigensolver.
This arbitrariness of the phase factor value must be taken into account when calculating the overlap terms ⟨ψ MCH α (R 0 )|ψ MCH β ( R)⟩ that enter Eq. 18. [34,35] In order to fix the phase, the wave functions at displaced geometries are transformed as The phase-correction factor f i is defined as and it is chosen such that the diagonal elements of the overlap matrix are real and positive.The vector ⃗ f then contains the values of the phase-correction factors f i for each of the MCH wavefunction |ψ MCH i ( R)⟩ as elements, and it is used to transform all the phase-dependent matrices expressed in the MCH basis, such as H MCH SO and the angular momentum operators J MCH i , L MCH i , and S MCH i (i = x, y, z).A further phase-correction step is required when H MCH SO is diagonalized to obtain the SO eigenfunctions.More precisely, we are interested in obtaining phase consistency between the rotational matrices U at displaced coordinates as these enter the expression of ∇ Ĥ through K U .When dealing with Kramers systems, the degeneracy of the SO states has also to be considered in addition to the arbitrary phase introduced by solving the eigenvalue problem.The phase tracking algorithm proposed by S. Mai et al. [36] is used to correct the phase of the rotational matrix U( R) with respect to U(R 0 ).Firstly, the overlap matrix U(R 0 ) † U( R) is computed and put into block-diagonal form by setting to zero all the matrix elements belonging to non-degenerate eigenstates.We refer to such a matrix as M and, since the eigenvalues of H MCH SO correspond to the energies of the Kramers doublets, each diagonal block is two dimensional.Subsequently, M is Löwdin orthonormalized to extract the matrix Φ such that the transformation U( R) → U( R)Φ † gives rotational matrices with the correct phase factors.In our implementation, Φ is obtained by computing the single value decomposition (SVD) of M, i.e. and Spin Hamiltonian.The most common way of approaching the problem of calculating relaxation times in SMMs is to go through the spin-Hamiltonian formalism that, by incorporating the relevant degrees of freedom and interactions, provides a simplified description that can effectively reproduce experimental observations.[37][38][39][40][41] Within this theoretical framework, it is possible to evaluate Eqs. ( 9)-( 13) in a reduced Hilbert space spanned by the ground J-multiplet of the bare magnetic ion of the system, where J denotes the total angular momentum quantum number.[9,19] The mapping of the total Hamiltonian Ĥ into the model space gives the generalised spin Hamiltonian Ĥs which can be expressed as where the operators Ôl m are a tesseral function of the total angular momentum operator Ĵ of rank l and order m.To obtain the spin Hamiltonian coefficients from first principles, the lowest 2J + 1 ab initio wavefunctions are put in a one-to-one correspondence with the magnetic ground state of the ion.To do so, the operator Ĵz in the basis of the lowest 2J + 1 ab initio wavefunctions is diagonalized to obtain the basis states of the spin Hamiltonian | J, mj ⟩.Note that this basis corresponds to Zeeman states |J, m j ⟩ (−J < m j < J) only if the 2J + 1 lowest ab initio wavefunctions are completely decoupled from the other excited states, i.e.only if the operator Ĵz in the ab initio basis is block diagonal.Thus, in practical computations, it is crucial to verify whether the eigenvalues of Ĵz in the initial ab initio basis closely resemble the expected values that would be derived for purely Zeeman states.If this correspondence holds true, the ab initio Hamiltonian is expressed in this new basis and rewritten using Steven's operators.[42] Finally the set of coefficients B l m are adjusted to reproduce the Hamiltonian matrix elements of the lowest 2J + 1 ab initio states In a similar way, we can define the operator which describes the coupling of spin and atomic displacements.We will refer to the elements ∂B l m /∂R a as spinphonon coupling coefficients and to ∇ Ĥ as vibronic coupling operator.The spin-phonon coupling coefficients can computed in two different ways, either by numerically differentiating the parameters B l m , or by projecting the vibronic coupling operator onto the lowest 2J + 1 ab initio states similarly to the static case, [43]

COMPUTATIONAL METHODS
Two complexes were chosen to conduct this study: geometry optimization as well as simulations of Γ-point phonons for the two compounds under investigation.[19] They utilized the Perdew-Burke-Ernzerhof (PBE) functional along with DFT-D3 dispersion corrections.[44,45] In the present work, we have reused the same dataset.
ORCA [46] had been used to compute molecular electronic properties using the state-averaged (SA) complete active space self-consistent field (CASSCF) method.For 1, the active space was built with seven electrons in the five 3d orbitals and all states with multiplicity six were considered.For 2, an active space of nine electrons in the seven 4f orbitals was considered and all states with multiplicity six were considered.The RIJCOSX approximation for the Coulomb and Exchange integral was used for both systems.The basis sets DKH-def2-QZVPP for Co atoms, DKH-def2-SVP for H and SARC2-DKH-QZVP for Dy atoms were used.DKH-def2-TZVPP basis set has been used for the rest of the atoms present in the systems.
An in-house Python code has been developed to evaluate ∇ Ĥ.This program reads AOs overlap, MOs, and the CI composition of the MCH wavefunctions from ORCA's output files and generates a valid input file for the program WFOVERLAP.[47,48] The latter is subsequently used to compute the overlap matrix between MCH states at different geometries needed to compute the non-adiabatic couplings ⟨ψ MCH i |∇ψ MCH j ⟩.Our program allows to evaluate ∇ Ĥ in the SO basis while correcting the phase of the electronic wavefunctions and of the rotational matrices.The numerical differentiation is performed using central differentiation around the equilibrium geometry with a step of 0.01 Å.
Second-and fourth-order time-dependent perturbation theories have been used to simulate one-and twophonon processes, respectively.The software MolForge is used for these simulations and it is freely available at github.com/LunghiGroup/MolForge[9].As discussed elsewhere, the simulation of Kramers systems in zero external field requires the use of the non-diagonal secular approximation, where population and coherence terms of the density matrix are not independent of one another.This is achieved by simulating the dynamics of the entire density matrix for one-phonon processes [9,49].An equation that accounts for the dynamics of the entire density matrix under the effect of two-phonon processes resulting from fourth-order time-dependent perturbation theory is not yet available.However, it is possible to remove the coupling between population and coherence terms by applying a small magnetic field along the magnetic easy-axis to break Kramers degeneracy [9].Here we employ the latter strategy to simulate Raman relaxation.

NUMERICAL RESULTS
Electronic structure and spin-phonon coupling.The computed energy spectra of the two compounds under investigation are reported in Figure 2. The energies of the lowest 2J + 1 states used to build the spin Hamiltonian span 280.2 cm −1 for 1 and 821 cm −1 for 2. It is noteworthy that the energy gap between the ground state J multiplet and the next first excited state is significantly smaller in the case of complex 1 compared to complex 2. In the former case, such energy separation amounts to 376.4 cm −1 and is comparable to the energy spacing between the first two Kramer doublets.On the other hand, for complex 2, the lowest 2J + 1 states are well separated from the first next excited state by an energy of 2424.5 cm −1 .
As a first step in our analysis, we benchmark the quality of the computed ∇ Ĥ matrix elements by comparing the derivative of the SH coefficients ∂B l m /∂R a obtained i) by differentiating B l m from equation 24 and ii) by using directly equation 27.For the latter, the ab initio eigenfunctions are written in terms of (pseudo) angular momentum states as and used to write the total Hamiltonian in the basis of the lowest 2J + 1 eigenfunctions as Finally, vibronic coupling matrix elements are computed using equations 15, 18, and 20.Because of the definition 28, the terms K SO ij (R 0 ) and K U ij (R 0 ) are now expressed with respect to the J-multiplet | J, mα ⟩ and the rotation matrix W(R 0 ).Since the Zeeman states | J, mα ⟩ contain only angular momentum degrees of freedom, the terms K SO ij (R 0 ) vanish and the vibronic coupling only depends on the energies E n of the ab initio eigenfunctions and on the rotation matrix W(R 0 ).The results are shown in Figure 3 where the coefficients ∂B l m /∂R a computed with the two approaches are compared.As a measure of the linear relationship between the two sets of values, we calculate the Pearson's Product Moment Correlation Coefficient (PPMCC) and the Root Mean Square Error (RMSE).For complex 1 the PPMCC and the RMSE are 0.9999 and 0.008 cm −1 Å−1 respectively, while for 2 the computed PPMCC is 0.9958 and the RMSE is 0.036 cm −1 Å−1 .
Relaxation time with the full Hilbert space.Analysis of the electron dynamics within ab initio open quantum system theory was conducted to compare the computed relaxation times τ obtained through both the proposed full Hamiltonian method and within the spin-Hamiltonian theory.Under the full Hamiltonian framework, the computation of matrix elements Ŵ 1−ph and Ŵ 2−ph is not confined to the subset of the Hilbert space defined by the ground state J-multiplet and Figure 4 shows Raman relaxation obtained using various Hilbert space sizes.We do not display relaxation times associated with the Orbach mechanism because they remain unchanged when expanding the dimension of the Hilbert space due to the lack of resonant phonons with the highenergy excited states introduced.In the case of complex 2, convergence is achieved immediately, and the Raman relaxation times remain consistent with those obtained using a minimal 2J + 1-dimensional Hilbert space.However, a different behavior is observed for complex 1, where convergence is only achieved when utilizing Hilbert spaces with more than twenty states.This difference is consistent with the spectra reported in Figure 2, where for complex 1, the energy of the first excited state outside the ground-state J-multiplet and the density of states is significantly lower compared to compound 2.
Comparison with the spin Hamiltonian.Next, we compare the converged relaxation time values obtained within the full Hamiltonian framework to those obtained using the spin-Hamiltonian approximation.Results are reported in Figure 5.The decomposition of the computed τ in terms of Orbach and Raman contributions and experimental extracted values of τ are also reported.At elevated temperatures, the dominance of Orbach relaxation is observed, whereas the influence of the Raman mechanism emerges only at lower temperatures.This behavior is mainly due to the temperature dependency of the phonon population.Orbach processes involve the absorption of a single phonon in resonance with the excited electronic states, if resonant phonons are thermally activated then the Orbach mechanism dominates the relaxation process.
As temperature decreases, resonant phonons are no longer populated and relaxation happens through absorption/emission of a pair of out-of-resonance low-energy phonons, i.e. the Raman mechanism.
The Orbach relaxation times, determined using both methodologies, generally display similar orders of magnitude, with the most significant deviation being a factor of two observed for compound 1.In contrast, substantial differences emerge when comparing the computed Raman relaxation times for both complexes.In the case of 1 the relaxation time obtained using the converged Hilbert space is three times shorter than the one predicted with the use of the spin Hamiltonian over the entire temperature range.This phenomenon is expected due to the growing number of excited states considered in the spin dynamics, which can facilitate coupling between the ground-state Kramers doublet m j = ±J, even at relatively low temperatures.Overall,   better agreement with experimental results is achieved when utilizing the full Hamiltonian framework.Turning our attention to complex 2, significant differences between the two methods emerge in the Raman relaxation mechanism at low temperatures.Within the full Hamiltonian approach, the relaxation times are FIG.6. Orbach (green, gray) and Raman (blue, black) relaxation times τ as a function of temperature computed using the full Hamiltonian approach with and without the inclusion of SO derivatives.In orange experimentally extracted date are reported.Left panel shows results for complex 1 and right panel for complex 2.
more than an order of magnitude smaller compared to both experimental results and the outcomes obtained using the spin-Hamiltonian framework.This substantial deviation persists even when considering only the lowest 2J + 1 ab initio states.One could be tempted to blame this deviation on some technical aspects of the calculations, primarily the displacement step utilized in the evaluation of the non-adiabatic coupling matrix elements.However, a strong correlation exists when comparing the coefficients ∂B l m /∂R a obtained with the two methodologies (see.Figure 3).This leaves us with two options: i) the process of constructing a spin Hamiltonian is less susceptible to numerical noise or ii) the spin Hamiltonian is not in this case able to fully capture the physics behind spin-vibronic coupling.
Spin-orbit coupling effect.So far, SO coupling and its dependency on the atomic coordinates have been fully accounted for through the evaluation of the terms K SO ij (R 0 ) and K U ij (R 0 ).Here, our intent is to show how considering the variation of the SO operator with respect to nuclear coordinates in the evaluation of vibronic coupling matrix elements can affect the computed relaxation times τ .To achieve this, we assume the term K U ij (R 0 ) to be negligible in Eq. 17 and proceed to reevaluate Ŵ 1−ph and Ŵ 2−ph .This approximation is equivalent to assuming that the SO coupling matrix elements remain constant at their equilibrium values when a small displacement in atomic coordinates is applied.In Figure 6, we show Orbach and Raman relaxation times calculated within the full Hamiltonian framework, both with and without the inclusion of SO derivatives, in comparison to experimental results.For both complexes 1 and 2, we observe that Orbach relaxation times are slower when the matrix elements K SO ij (R 0 ) are excluded from the calculation.In the case of 1, omitting the SO derivative leads to Orbach relaxation times that are over two orders of magnitude larger.However, for compound 2, the effect of the SO derivative only marginally impacts the Orbach relaxation times in the high-temperature regime.Regarding the Raman relaxation mechanism, we note dis-tinct behavior for the two systems under investigation.In 1, removing the K SO ij (R 0 ) matrix elements results in relaxation times that are up to three times larger compared to the full Hamiltonian picture.Conversely, in complex 2, the Raman relaxation times are consistently smaller when the SO derivative is excluded, with a difference of one order of magnitude in the low-temperature regime.Furthermore, in all cases, significant changes between Raman and Orbach relaxation times are observed in the high-temperature limit.

DISCUSSION
The presented numerical method allows the prediction of Orbach and Raman relaxation times in SMMs using ab initio wavefunctions from electronic structure methods avoiding the construction of a spin Hamiltonian, thus generalizing previously proposed strategies [4].
Overall, the results are in good agreement with experiments, but interestingly the inclusion of the full electronic Hilbert space improves over spin Hamiltonian results for [Co(C 3 S 5 ) 2 ] − but worsens them for [Dy(bbpen)Cl].We believe that multiple effects are at play.On the one hand, the results for [Co(C 3 S 5 ) 2 ] − , clearly show that multiple electronic excited states can contribute to the Raman relaxation, going beyond the ground-state spin multiplet.This is only possible with the proposed approach and it is an effect that must be fully accounted for going forward.On the other hand, the use of a full Hilbert space leads to the largest deviation with respect to experiments in [Dy(bbpen)Cl] despite being the most accurate method in principles.This points to a potential presence of error cancellation effects in the use of the spin Hamiltonian in [Dy(bbpen)Cl] and highlights the need for extra care in comparing simulations and experiments.In general, discrepancies between theory and experiments up to one order of magnitude are not unprecedented for this kind of simulation.These differences might be in part ascribed to the presence of relaxation mechanisms that extend beyond the scope of this study, particularly the absence of spin-spin dipolar cross-relaxation.[50,51] Fluctuations in the computed relaxation times can also be attributed to the specific limitations of the ab initio methods.In this regard, three main points deserve further analysis: i) phonon calculations are limited to the Γ-point, neglecting any effects arising from acoustic phonons and the dispersion of optical modes.As shown previously, this effect might lead to substantial deviations at low temperature; ii) whilst the CASSCF electronic structure method is able to capture the multiconfigurational nature of SMMs, it is not efficient in accounting for dynamical correlation effects, which could be significant in this context; iii) the evaluation of wavefunction overlap through finite differences is a general and powerful approach, but a careful convergence of NAC vectors with respect to the differentiation step size is required.
Going beyond methodological considerations, the results presented here shed light on the origin of spin-phonon coupling at the quantum mechanical level.For instance, the proposed method is able to isolate the key role of SO coupling derivatives in the study of spin-phonon dynamics in SMMs.To the best of our knowledge, the importance of this term has not been considered in previous studies, [43].Moreover, the study of the Dy compound shows evidence that the spin Hamiltonian approximation might not always be fully justified and that the full wavefunction contains additional information.We envision that the extension of our method to other molecular complexes will further help unravel the contributions to spin relaxation of pure electronic origin.The present study also paves the way for a systematic exploration of the role of electronic excited-state dynamics in open-shell transition metal and lanthanide-based coordination compounds and its interplay with spin dynamics.[52][53][54] Indeed, effects such as inter-system crossing and internal convention also find their origin in vibronic coupling and this work strongly supports the possibility of extending these simulations to account for those processes.
In conclusion, we have presented a novel computational method able to accurately describe the dynamics of the ground state magnetization of open-shell systems under the effect of vibronic coupling up to two-phonon relaxation processes.Our method generalizes previous approaches based on the effective spin Hamiltonian and shows the importance of including the effect of electronic excited states.Moving forward, this work represents a pivotal point toward delivering a complete and accurate picture of spin dynamics from first principles.

FIG. 2 .
FIG. 2. Energy eigenvalues in cm −1 computed for 1 (left) and 2 (right) up to 8000 cm −1 .The eigenvalues corresponding to the lowest 2J +1 states are highlighted in red.The scale's zero point is aligned with the energy of the first Kramers doublet.

FIG. 3 .
FIG. 3. Comparison between the computed matrix elements ∂B l m /∂Ra obtained by numerical differentiation of B l m (Method 1) and by projecting NACs between ab initio wavefunctions into the spin Hamiltonian subspace (Method 2).Left panel shows results for complex 1 and right panel for complex 2.

FIG. 4 .
FIG. 4. Raman Spin-phonon relaxation times τ as a function of temperature for different sizes of the Hilbert space used to compute the two-phonon transition rate Ŵ 2−ph .In orange experimentally extracted data are reported.The left panel shows results for complex 1 and the right panel for complex 2.

FIG. 5 .
FIG.5.Orbach (green) and Raman (blue) relaxation times τ as a function of temperature computed within the spin-Hamiltonian framework (Method 1, empty dots) and by employing the full Hamiltonian approach (Method 2, filled dots).In orange experimentally extracted date are reported.Left panel shows results for complex 1 and right panel for complex 2.