Analytic Linear Vibronic Coupling Method for First-Principles Spin-Dynamics Calculations in Single-Molecule Magnets

Abstract

Accurate modeling of vibronically driven magnetic relaxation from ab initio calculations is of paramount importance to the design of next-generation single-molecule magnets (SMMs). Previous theoretical studies have been relying on numerical differentiation to obtain spin-phonon couplings in the form of derivatives of spin Hamiltonian parameters. In this work, we introduce a novel approach to obtain these derivatives fully analytically by combining the linear vibronic coupling (LVC) approach with analytic complete active space self-consistent field derivatives and nonadiabatic couplings computed at the equilibrium geometry with a single electronic structure calculation. We apply our analytic approach to the computation of Orbach and Raman relaxation rates for a bis-cyclobutadienyl Dy(III) sandwich complex in the frozen-solution phase, where the solution environment is modeled by electrostatic multipole expansions, and benchmark our findings against results obtained using conventional numerical derivatives and a fully electronic description of the whole system. We demonstrate that our LVC approach exhibits high accuracy over a wide range of coupling strengths and enables significant computational savings due to its analytic, “single-shot” nature. Evidently, this offers great potential for advancing the simulation of a wide range of vibronic coupling phenomena in magnetism and spectroscopy, ultimately aiding the design of high-performance SMMs. Considering different environmental representations, we find that a point charge model shows the best agreement with the reference calculation, including the full electronic environment, but note that the main source of discrepancies observed in the magnetic relaxation rates originates from the approximate equilibrium electronic structure computed using the electrostatic environment models rather than from the couplings.

Impact of Finite Phonon Lifetimes on the Spin Dynamics of Single-Molecule Magnets

“…Hamiltonian acting in the (2J+1)-dimensional |mJ⟩ basis. 50 The spin-phonon coupling Hamiltonian for each mode is evaluated using a linear vibronic coupling model, where we first determine the first derivatives of the CF parameters for each atomic Cartesian degree of freedom, then convert these into the normal mode basis using the appropriate linearcombination of atomic displacements from the normal mode vector. 50 This is done for all phonon modes, at all q-points in the sampling mesh.…”

“…Magnetic relaxation rates are determined using our code Tau, considering the Orbach and Raman-I rate expressions derived from perturbation theory, 14,22 which are given by Equations 1-6 in the supporting information of Ref. 50 . Calculation of the Raman rates using these perturbative expressions requires an energy gap between the lowest Kramers doublet, which occurs in experiments due to the presence of a dipolar magnetic field and/or the driving AC magnetic field in the experiment, and thus we apply a magnetic field of 2 Oe along the main magnetic axis of the molecule, splitting the ground doublet by ca.…”

“…The Raman mechanism involves two phonons, with no restrictions on the energies and wavevectors, and the rate is obtained as the double-integral over pairs of modes. 22,50 We restrict the domain of the phonon energies to 0 ≤ 𝜔 ≲ 𝜔 GHI , where 𝜔 GHI = 267 cm -1 is chosen as the minimum in the DoS above the low-energy pseudo-acoustic peak. The cut-off is applied to avoid divergences in the in the Raman-I rates (where the denominator of Equations 3-6 in the supporting information of Ref.…”

“…The cut-off is applied to avoid divergences in the in the Raman-I rates (where the denominator of Equations 3-6 in the supporting information of Ref. 50 would go zero when ℏ𝜔 is resonant with a CF excitation), and is sufficiently smaller than the first crystal field excitation of 421 cm -1 . The double-integral is transformed into a one-dimensional integral due to the conservation of energy via the Dirac delta function, and is performed over anti-Lorentzian phonon lineshapes (Equation 7 in the supporting information of Ref.…”

“…The double-integral is transformed into a one-dimensional integral due to the conservation of energy via the Dirac delta function, and is performed over anti-Lorentzian phonon lineshapes (Equation 7 in the supporting information of Ref. 50 ) to an equivalent range of 𝜇 ± 2𝜎 (95%) using the trapezoidal method with 40 equidistant steps; anti-Lorentzian lineshapes are used to ensure that the DoS goes exactly to zero at zero energy.…”