Open Access Article

This Open Access Article is licensed under a

Creative Commons Attribution 3.0 Unported Licence

Luis
Escalera-Moreno†
^{a},
José J.
Baldoví†
^{b},
Alejandro
Gaita-Ariño
*^{a} and
Eugenio
Coronado
*^{a}
^{a}Instituto de Ciencia Molecular (ICMol), Univ. de Valencia, C/Catedrático Beltrán 2, E-46980 Paterna, Spain. E-mail: alejandro.gaita@uv.es; eugenio.coronado@uv.es
^{b}Max Planck Institute for the Structure and Dynamics of Matter, Luruper Chaussee 149, D-22761 Hamburg, Germany

Received
26th December 2017
, Accepted 7th March 2018

First published on 7th March 2018

Very recently the closely related fields of molecular spin qubits, single ion magnets and single atom magnets have been shaken by unexpected results. We have witnessed a jump in the phase memory times of spin qubits from a few microseconds to almost a millisecond in a vanadium complex, magnetic hysteresis up to 60 K in a dysprosium-based magnetic molecule and magnetic memory up to 30 K in a holmium atom deposited on a surface. With single-molecule magnets being more than two decades old, this rapid improvement in the physical properties is surprising and its explanation deserves urgent attention. The general assumption of focusing uniquely on the energy barrier is clearly insufficient to model magnetic relaxation. Other factors, such as vibrations that couple to spin states, need to be taken into account. In fact, this coupling is currently recognised to be the key factor that accounts for the slow relaxation of magnetisation at higher temperatures. Herein we will present a critical perspective of the recent advances in molecular nanomagnetism towards the goal of integrating spin–phonon interactions into the current computational methodologies of spin relaxation. This presentation will be placed in the context of the well-known models developed in solid state physics, which, as we will explain, are severely limited for molecular systems.

These experimental records have been supported by advances in theoretical modelling, but a precise description of spin dynamics at the nanoscale is still extremely challenging. For many years, modelling of slow magnetic dynamics in nano-objects, such as single-molecule magnets and single-ion magnets, relied mostly on the Orbach mechanism. The effective barrier for the reversal of magnetization is now routinely estimated from first principles,^{20,21} which allows a rational design of these nanomagnets.^{22–24} In contrast, Raman processes are very often taken into account only parametrically. This evidences that control over spin dynamics in molecular nanomagnets is still an open problem, which requires the modelling of vibrations and of their coupling to the spin energy levels from first principles.

Herein, we discuss the current difficulties in the search for a relationship between the molecular structure and spin dynamics. In order to get an appropriate perspective, we will start from the achievements and drawbacks of a static picture that aims to correlate chemical structures with the spectroscopic and magnetic properties of molecular nanomagnets; then, we will pass through the problematic focus on what we argue is the first stage of this problem – the energy barrier – and, finally, we will review what is being nowadays recognised as the current stage of this problem: the role of vibrations.

In contrast, the magnetic properties of SIMs and mononuclear spin qubits are largely determined by the magnetic anisotropy of a single ion, which results from the combination of spin–orbit coupling and the crystal field.^{28} The relative strength of such electronic interactions relies on the electronic configuration of the magnetic centre, with remarkable differences between d-block (ligand field > spin–orbit coupling) and f-block element ions (ligand field < spin–orbit coupling).^{29} In the latter we can also distinguish between lanthanides and actinides.

Crystal field theory is key for the description of the energy level scheme. This frequently requires the determination of a large number – up to 27 – of crystal field parameters (CFPs). The estimation of CFPs can be done based on a few alternatives. The first one is the rationalisation of the experimental features of complexes already synthesised and characterised empirically. This has traditionally been the default option of spectroscopists and consists in the direct fit of spectroscopic data while varying a non-negligible set of CFPs. The non-vanishing CFPs depend on the point group of symmetry of the molecule.^{30} An accurate description of magnetic properties following a phenomenological approach has also included thermodynamic techniques, such as powder and single-crystal magnetic susceptibility and torque magnetometry.^{31} The second option is using a computational approach to obtain CFPs and then modelling magnetic properties, or even to predict them using the real chemical structure of the coordination complex as an input. In this direction, there are mainly two alternatives that have proven to be useful in molecular magnetism, namely the electrostatic crystal-field approach, which considers a point-charge distribution around the central ion,^{32–35} and the more expensive ab initio calculation, based on the complete active space self-consistent field (CASSCF).^{36–40} Comparisons between these two approaches have been made elsewhere.^{41,42}

Partly guided by theoretical efforts, a rich community of experimental chemists and physicists has worked for decades to increase the energy barrier with the objective of achieving high blocking temperatures, ideally up to room temperature. Popular chemical families include beta-diketonates, aromatic rings, polyoxometalates and phthalocyaninato anions as ligands, and the consensus seems to favour Kramers ions with an oblate f-shell distribution (and notably the Dy^{3+} ion), with an axially elongated coordination environment, rigid polyhapto ligands and diamagnetically diluted samples. There are now record barriers (assuming an Orbach mechanism) that are at least an order of magnitude higher than those reported in cluster-type SMMs. At the same time, SIMs working at room temperature are still a distant dream. One of the main reasons is that the employed energy-barrier framework is an oversimplification. Thus, while recognizing the important victory of being able to systematically design and prepare systems with higher effective thermal barriers, we need to put this into perspective.

Three figures of merit have been frequently employed in the analysis of the dynamical magnetic properties of SMMs:

(1) The effective barrier U_{eff}, which can be quantified by fitting the variation of the ac susceptibility signal with frequency and temperature to the Arrhenius equation.

(2) The first excited magnetic energy level E_{1}, which is ideally determined by spectroscopy, but frequently just estimated by theoretical calculations.

(3) Hysteresis loops, which can be observed below a certain blocking temperature. This is actually the critical parameter that bars the gate for devices and applications, although in the field it is common to see “blocking temperature” used in relation to the ac magnetometry signal, since not all complexes display hysteresis at 2 K.

It has often been assumed that U_{eff} and E_{1} are identical and directly control the blocking temperature. However, individual studies have repeatedly shown that this is not the case,^{43} with U_{eff} being significantly lower than E_{1} (which has recently been attributed to the presence of off-resonance phonons due to the finite phonon lifetimes which offer a wider energy window),^{44} or with both U_{eff} and E_{1} being two orders of magnitude higher than the blocking temperature.^{45} What is happening? It is likely that there is no single answer, but it seems clear that all relevant physical processes – including Orbach (Or), Raman (Ra) and direct (Di) mechanisms – should be taken into account in each case. This is conceptually not so different from a simple electric problem, which we shall use for illustration purposes. Let us picture two electric circuits (see Fig. 2), one in series and the other in parallel, for which we want to minimise the current flow, just like we want to minimise spin relaxation in our molecular magnets. What is the simplest systematic strategy to increase the overall resistance in a simple circuit? In the series circuit, we can just pick any resistor, say R_{Or}, and raise its resistance, and R_{series} will escalate with no limit. In the parallel circuit the situation is different: when R_{Or} rises over a certain threshold the current flows exclusively through R_{Ra} and R_{Di}—the paths of least resistance—making R_{Or} an irrelevant part of the circuit.

Fig. 2 Series circuit (left) versus parallel circuit (right). Modified with permission from Mets501 (CC by-sa 3.0) series circuit and parallel circuit. |

Back to molecular magnetism, we qualitatively have a similar situation: given several relaxation pathways, the spin will most commonly relax via the fastest one. It is therefore easy to understand that, after the thermal barrier has risen over a certain threshold, the spin will just use a different relaxation mechanism. So, further raising the barrier and thus blocking the path of most resistance will be irrelevant for all practical purposes. In the case of SIMs, the community has already done a good job in raising the barrier and is now starting to admit that molecular vibrations are the next pathway that needs to be blocked. Actually, Liddle and van Slageren, in a tutorial review published in 2015, already pointed out that the magnitude of the crystal field splitting is not the only factor governing the slow relaxation of molecular nanomagnets. They highlighted the importance of Raman processes and explicitly indicated the necessity, for energy dissipation, of transferring the energy, via phonons, from the spin system to the thermal bath.^{29} For the study of these phenomena, we need to take advantage of the tools, methods and concepts that were developed by physicists working on the thermal dependence of the crystal field Hamiltonian. Thus, we will now overview such reports, summarizing the most crucial equations that relate vibrations and spin energy levels, as well as the open problems that still need to be addressed.

(i) The calculation of the static contribution B_{stat}

(ii) The use of the Debye model to obtain the contribution B_{ac} of the acoustic branches of the phonon spectrum

(iii) The approximation of optical branches using a single-phonon model to obtain their collective contribution B_{op}.

Following this scheme, the thermal dependence B(T) is decomposed into three terms:

B(T) = B_{stat}(T) + B_{ac}(T) + B_{op}(T) | (1) |

Regarding step (i), the calculation of the static contribution was generally achieved either by means of thermal expansion coefficients—either by a complete diagonalisation or by using perturbation formulae—or by using point-charge models.

The following step (ii) is central to our interests and relies on the Debye model (explained in Fig. 3), which is known for being the first approach to correctly reproduce the behaviour of the specific heat in simple solids. Within this model, one starts with the generalised coordinates that describe each atomic motion in a solid according to a given phonon k:^{50}

(2) |

The main difficulty in eqn (2) usually comes from the evaluation of the exponential factor. An approximated solution is possible by combining the assumption of cubic symmetry with the so-called long wavelength approximation (i.e. |·| ≪ 1). Thus, the phase factor e^{i}^{·}^{} can be approximated as ·. Here one finds the first remarkable limitation for our purposes, as complex molecular crystals generally lack cubic symmetry. Moreover, systems where the LWA fails can be encountered in the literature, especially when the working temperature is of the order of the Debye temperature or higher.^{55} Indeed, as the temperature is raised, phonons of short wavelength are also excited and thus the integral in eqn (4) cannot properly describe a correct temperature dependence. In the context of molecular crystals, which present Debye temperatures of the order of tens of kelvin,^{56} this failure is expected to appear much below room temperature.

An expansion of B_{ac} and B_{op} in terms of these coordinates gives rise to expressions that depend on expectation values and so on. Generally, under an anharmonic phonon model for atomic displacements, may be non-zero. In contrast, if the model is harmonic, this expectation value is identically zero, while would be the first term in the expansion different from zero. The expressions derived by Shrivastava are truncated at second order and consider harmonic phonons. Thus, only the effect from quadratic atomic displacements is incorporated,^{50,57} and the phonon-induced modulation of B_{ac} is proportional to . Whereas this approach has also been successfully recovered by the models recently proposed for magnetic molecular crystals,^{48} we need to point out that there is a second important limitation. There are relevant anharmonic effects, such as lattice spacing or phonon–phonon interactions, especially at high temperatures, which cannot always be safely ignored.^{60} Recently, the relevance of anharmonicity in phonons for spin dynamics has already been the subject of study in the context of magnetic molecules.^{44}

As phonon energies are close enough to describe a continuum, the series is usually converted into an integral, which describes the overall contribution of the acoustic phonon spectrum to B for harmonic quadratic atomic displacements. This integral introduces a third limitation: it impedes determining which phonon modes in contribute the most to modulate B_{ac}. Since one of the main relaxation pathways in molecular spin qubits and SMMs can be via spin–vibration coupling, it would be desirable to be able to check each individual mode contribution in order to rationally design these molecular systems and slow down this relaxation. Fortunately, this is trivial to do, simply by keeping the series expression instead of switching to the integral.

Before is converted into an integral, the sum over the square atomic locations involved in each collective motion is substituted by a mean value R^{2}, with R being the lattice nearest-neighbour distance in the considered crystal of cubic symmetry. The conversion of the series into an integral then gives:

(3) |

(4) |

The Debye model presents practical and fundamental limitations. Experimentally determined values of Θ_{D} sometimes differ by tens or even hundreds of Kelvin depending on the technique.^{62–64} In other cases, the thermal dependence of Θ_{D} in eqn (4) is employed as a last resource to fit room-temperature data.^{64–66} It may work and provide practical applications,^{60} but makes Θ_{D} unphysical since eqn (4) should only be used as long as the LWA is fulfilled, i.e., at not too high temperatures.

Note that the Debye model is useful for simple solids but has a limited applicability in molecular solids. First, this model assumes a specific phonon spectrum, which could fail in complex crystals of a rather general symmetry where the paramagnetic entities cannot be considered zero-dimensional anymore. Second, the dispersion relation ω = v·k has been used to express the integral in eqn (4). Although this relation is frequently employed, it might not work for some systems and should be consequently replaced by another one depending on the specific structure and properties of the crystal. Third, there is a well-known danger of using the Debye model at high temperatures,^{46} which is now the most interesting regime for the communities of molecular magnetism and spin qubits.^{15,67,68} Already in 1969 it was claimed that theoretically probing the high temperature regime would only be possible whenever non-Debye calculations were available, which should be point-to-point calculations.^{46} These calculations should consider explicitly the exact phonon spectrum of each particular crystal, and possible angular^{54} and thermal dependencies of sound velocities in a given crystal. Others have elaborated on this point, discussing about the replacement of the Debye phonon density by the real one.^{52,64,66} Over the last few years, it has been pointed out that fitting temperature dependencies of spin–lattice relaxation times sometimes requires using the real phonon density.^{69,70}

Finally, in step (iii), the contribution of the optical branches is accounted for by all the optical modes:

(5) |

The phonon optical branch is described using a single-mode harmonic model, with effective frequency ω_{eff} and distortion coordinate Q. This constitutes the fourth main limitation in this procedure, which is usually justified by stating that optical phonons have been usually found in narrow frequency ranges. In the case of molecular solids, this is no longer the case, since the frequencies of molecular vibrations span over two orders of magnitude.

The calculation of 〈Q^{2}〉 results in:

(6) |

(7) |

Even with these four main limitations, this procedure has produced useful insights into the theoretical rationalization of spin dynamics in magnetic molecules. For illustration, we can comment on an electron spin relaxation study by Eaton and co-workers that considered a series of Cu(II) complexes in a wide temperature window.^{72} The semi-empirical model they used invoked contributions from several relaxation processes achieving an excellent reproduction of the thermal dependence of the spin–lattice relaxation time, T_{1}. An almost temperature-independent direct process was found to be significant below 20 K, Raman processes dominated between 20 K and 60 K, and local modes of energies around 300 K (200 cm^{−1}) were found to be very significant already at temperatures of 60 K and above. Since no low-lying electronic states are expected for Cu(II) complexes, these authors did not even consider Orbach processes, in contrast to the previously mentioned excessive focus on the barrier that has been so pervasive in the SMM community.

The same authors demonstrated how detailed experimental information can be useful to find the correct relaxation mechanism. In the case of bis(diethyldithiocarbamato) copper(II), Cu(dtc)_{2} (chemically diluted into a diamagnetic analogue), T_{1} was found to be frequency independent. This ruled out a mechanism involving a thermally activated process and instead indicated that relaxation proceeds via a local mode. In general, distinguishing between a local mode and a thermally activated process requires experimental data at temperatures up to or beyond the temperature corresponding to this characteristic energy or relaxation measurements with at least two different microwave frequencies.

Besides all the above mentioned limitations, this whole methodology is semi-empirical.^{57,63,73–77} The relevant parameters are extracted by means of fittings to experimental data, from independent experiments or tabulated values, or simply estimated. Without an independent predictive capability, this means that these models cannot facilitate a rational molecular design.

(i) One starts by relaxing the relevant geometry and calculates its vibrational spectrum. Depending on the case, this geometry may only involve atoms of the magnetic complex or, additionally, atoms of the rest of the unit cell. This is commonly undertaken via DFT (Density Functional Theory), although recent approaches based on general force fields are also available.^{79} However, DFT does not necessarily guarantee a decrease in the overall system energy as the atomic orbital basis is enlarged. Thus, depending on previous experience, finding a systematic method to relax the geometry may become a hard task. In this context, there are three methodological issues that need to be discussed. First, one should not guide the relaxation process aiming for a perfect match between the X-ray structure—usually extracted at T ≥ 100 K—and the relaxed geometry, which is at the absolute energy minimum of the chosen theoretical method.^{29,48,67} This can be solved either by low-temperature crystallography or by correcting high temperature effects such as libration in the ≥100 K experimental geometry. Second, if the steric pressure by the environment is crucial for the molecular structure, it needs to be taken into account by including, for example, a set of frozen nearest counter-ions during the geometrical optimisation.^{48} Only occasionally can this be dropped and perform a relaxation in a vacuum. Third, usually the calculation of phonons is computationally demanding; thus, only one or a few directions in the reciprocal space are chosen, which limits the physical value of the results. For instance, if only the unit cell gamma-point is taken into account, only vibrational modes restricted to a single unit cell can be taken into account, meaning that the vast majority of intermolecular modes are neglected.^{44,78} Note that the more modes are included in the model, the more likely it will be to find all the potential spin relaxation channels. Whenever the number of modes becomes too large, one will have to sense and select only those ones that could contribute the most to relaxation.

(ii) In a second step, the relaxed geometry is distorted following each vibrational mode, generating a finite set of distorted geometries. For each mode, the selection of the lower and upper bounds of the distortion coordinate is not unique, but one criterion may be distorting the geometry until reaching the energy of the first excited vibrational state. Thus, one can safely use the energies of a harmonic oscillator. Likewise, the criterion used to choose those discrete values that each distortion coordinate must take is not unique yet either.^{33,48}

(iii) Once the relevant magnetic anisotropy parameters affecting the spin relaxation are identified, they can be extracted from each distorted geometry by means of either ab initio or DFT point calculations. For instance, among these parameters one can find the g factor in spin–1/2 molecular spin qubits,^{48} or crystal field parameters such as the ZFS.^{15,44,78} A critical approximation at this step arises when periodic boundary conditions are not incorporated into these point calculations, so they are performed on a single isolated molecule. This may impose a severe limitation on the calculation quality, as long-range effects derived from the presence of charged species in the crystal are being completely removed. One can partially overcome this problem by placing near point charges simulating the outer electrostatic shells of the molecule and testing whether these effects are important or not.^{44}

(iv) Finally, the spin–vibration coupling is introduced as the modulation that each vibrational coordinate exerts on the relevant parameters. This coupling is characterised via derivatives of these parameters with respect to the vibrational coordinates, which are calculated either analytically or numerically and employing the results of the above point calculations. At this point, it is still important to develop and agree on a robust procedure to calculate these derivatives. This can be an important source of numerical error depending on their quality.^{33,48,67} The remarkable achievement of this methodology is that assumptions derived from the Debye model are no longer required as spin–phonon coupling coefficients are individually and explicitly evaluated.

The last stage that completes this process and connects with measurable magnitudes like magnetic relaxation times or spin decoherence times is the inclusion of the calculated spin–phonon coefficients in an appropriate master equation.

Initially, we need to focus on the record value for the spin–spin relaxation time T_{2} = 675 μs (at T = 10 K) which was set in 2015 by the vanadium complex [V(C_{8}S_{8})_{3}]^{2−} (Fig. 5(d)) with perdeuterated tetraphenylphosphonium counter-ions and in a diluted frozen CS_{2} solution, avoiding nuclear spins.^{9} At low temperature, T_{2} is governed by temperature-independent interactions with the spin bath, which in this case is unusually low. However, as the number of active phonons increases with rising temperature, T_{1} decreases and eventually limits T_{2}. In this case, T_{2} decreases by about an order of magnitude for every rise in temperature of 30–40 K, signalling a strong coupling between spin states and vibrations.

In the opposite extreme, one finds vanadyl phthalocyanine VOPc,^{68} with an almost constant T_{2} = 1 μs between 5 K and 300 K (as usual, at high dilutions). For VOPc, a preliminary analysis attributed the high values of T_{1} and T_{2} at high temperatures to the rigidity of the vanadyl moiety, and such rigidity was also shown to be important in another related study.^{67} If the V–O vibration is the only one that couples with the spin state, its frequency would govern the temperature at which T_{1} starts to be short. It would be tempting to speculate on whether the marked difference between [V(C_{8}S_{8})_{3}]^{2−} and VOPc in the thermal dependence of T_{1} is related to the vanadyl moiety or to their very different environments—a crystal compared with a frozen solution—but since this is a multifactorial problem, calculations are required before jumping to conclusions.

Another well-studied case in this series is [Cu(mnt)_{2}]^{2−},^{13,48} which displays T_{2} = 68 μs at low temperature and T_{2} = 600 ns at 300 K. We identified the two modes with the strongest spin–vibration coupling in the range 100–300 K; these involve distortions outside the molecular plane.^{59} In fact, this molecule is planar in the crystal, but theoretically relaxes in a vacuum to a skewed structure. One can extrapolate that increasing the chemical pressure driving the planarity of the complex will raise the vibrational frequency of these modes, thus decreasing their availability as spin relaxation paths. This could result in the survival up to higher temperatures of long spin–lattice relaxation times T_{1}, and, indirectly, in long spin–spin decoherence times T_{2} up to higher temperatures. In contrast, frozen solution samples permit T_{1} measurements in the absence of crystalline pressure and, indeed, these measurements revealed shorter T_{1} times for this complex.^{80}

Lattice vibrations were equally critical in explaining the behavior of dysprosocenium,^{15} a record-setting single-ion magnet based on a single Dy(III) cation sandwiched between two (1,2,3)tri-tert-butyl pentacene anions (Fig. 5(a)). This system presents a crystal field splitting of about 1500 cm^{−1} (Fig. 5(b)), magnetic hysteresis at temperatures of up to 60 K and an effective barrier U_{eff} = 1223 cm^{−1} (1760K), something that could naïvely be identified with an extremely uniaxial coordination environment in an ideal geometry. In fact, this was claimed in a parallel discovery of the same record SIM.^{16} This claim does not correspond to the reality of the molecular structure: the complex has a bent shape and bears no correspondence with any ideal symmetry, pseudoaxial or not. Instead, its unique spin dynamics were related to an equally unique spin–phonon coupling of the constrained metal–ligand vibrational modes, intrinsic to the bis-η^{5}-Cp^{ttt} coordination geometry. In particular, four modes have been identified as detrimental, in the sense of coupling strongly to spin states that participate in spin relaxation. These modes involve motion of the two C–H groups on each aromatic ring. This moved the authors to suggest the substitution of these groups by heavier analogues. In cases like this, deuteration would have a minor practical effect compared with halogenation or substitution by an organic group R, but at the same time it would allow a cleaner theoretical analysis since the static crystal field effect would be essentially intact.

Let us finally focus on single atom magnets, in which a neutral atom sits on the top of an insulator. In the first studied example, the neutral atom is Ho and a thin MgO layer separates it from an Ag substrate (Fig. 5(c)).^{14} When physicists described this minimalistic system they highlighted the role of the MgO layer as “protecting the quantum magnet from scattering with electrons and phonons of the substrate” or, in other words, decoupling it from the phonon and conduction electron baths. In this case, instead of blocking the detrimental modes that couple with the electron spin as we suggest in the chemical approach, this is achieved by choosing a system with a low phonon density of states such as MgO. This is simple and effective, but apparently precludes the possibility of a progressive chemical optimization.

A molecular-based analogue of this construction employed a whole monolayer of terbium bis-phthalocyanine complexes TbPc_{2} on MgO/Ag(100), rather than a single Ho atom.^{81} It was reported that magnetic remanence and hysteresis opening obtained with TbPc_{2} on MgO tunnel barriers outperform the ones of any other surface adsorbed SMM as well as those of bulk TbPc_{2}. However, hysteresis disappears above 8 K in this molecular monolayer, in contrast to a blocking temperature of 30 K for single Ho atoms. Whereas the phonon spectrum of MgO is equally poor in both cases, the difference might lie in the richer phonon spectrum available to a compact monolayer of TbPc_{2} molecules compared to isolated and thus “cleaner” Ho atoms.

Finally, let us discuss a last issue from the point of view of calculations. As previously stated in eqn (1), the dynamical contribution to the temperature dependence of the magnetic anisotropy consists of two parts: the contribution of acoustic phonons and that of optical phonons. Over time, the importance of including both acoustic and optical phonons in spin dynamics has been stressed.^{44,53,74–77,82} Indeed, some studies had to be revisited for not properly including the effects of both of them.^{49,83} Thus, any theoretical model to be developed should first consider the importance that both acoustic and optical phonons could have on spin relaxation before arbitrarily neglecting either of them.

A well-known problem which is nevertheless not routinely considered in molecular spin dynamics is the first static contribution in eqn (1), despite having been repeatedly proven to be crucial in other contexts.^{48,53,74–77} Recently, in the relevant experimental regime, magnetic relaxation times (extensible to spin–lattice decoherence times) are being calculated by assuming that relaxation is mainly phonon-induced.^{15,44} But, this static contribution can also play a key role in the thermal modulation of spin Hamiltonian parameters, so its effects should be considered in calculating these relaxation times. Moreover, such a static contribution may become important at temperatures higher than the nitrogen boiling point and thus would have to be considered in this thermal regime. Although ab initio calculations on spin dynamics are close to recovering the experimental order of magnitude of relaxation rates, discrepancies like temperature independent shifts between experimental and calculated relaxation times still remain.^{15} Already in 1973, it was discussed that the inclusion of the exact density of phonon states instead of the Debye ω^{2} value can automatically incorporate the effect of the lattice thermal distortion.^{48,60} Thus, if this static effect is finally proven to be important in calculating magnetic relaxation rates, its proper inclusion could be undertaken by considering the exact phonon spectrum or, alternatively, its most relevant parts.

Nowadays, the exact phonon spectrum in molecular systems is not considered in the calculations of the spin dynamics. For example, the relaxation dynamics of [Dy(Cp^{ttt})_{2}]^{+} were consistent between the crystalline phase and the amorphous frozen solution. Thus, localized molecular vibrations were assumed to govern the spin dynamics, and, therefore, only the gas-phase vibrational modes were considered.^{15} In the discussion, the authors pointed out that this oversimplification was a possible cause behind the fact that the temperature dependence of their calculated Raman mechanism deviates considerably from the experiment. A second example is provided by the molecule [(tpa^{Ph})Fe]^{−}. In this case, acoustic phonons were included in the modelling,^{44} but only in a minimal expression, namely the unit cell gamma-point normal modes. Again, this simplification might be behind the order-of-magnitude deviation between predicted and experimental relaxation times.

Notice that given the importance of an exact phonon dispersion in determining spin relaxation at elevated temperatures, it will be crucial to obtain insight into the environmental phonons that can contribute to these relaxation processes. Indeed, one needs to recall that phonons are an essential part of the dissipation pathway towards the thermal bath.^{29} In this pathway, local vibrations play the role of a link between spin states and phonons.^{84} Thus, the relationship between molecular vibrational modes and the lattice phonons that contribute to spin relaxation should be investigated, which can only be done if we have access to a good estimate of the phonon dispersion. This would allow us to understand which structural features of the magnetic molecules control their relation to the bulk environment and thereby govern spin relaxation.

- R. Sessoli, D. Gatteschi, A. Caneschi and M. A. Novak, Nature, 1993, 365, 141–143 CrossRef CAS.
- R. Sessoli, H. L. Tsai, A. R. Schake, S. Wang, J. B. Vincent, K. Folting, D. Gatteschi, G. Christou and D. N. Hendrickson, J. Am. Chem. Soc., 1993, 115, 1804–1816 CrossRef CAS.
- M. Mannini, F. Pineider, P. Sainctavit, C. Danieli, E. Otero, C. Sciancalepore, A. M. Talarico, M.-A. Arrio, A. Cornia, D. Gatteschi and R. Sessoli, Nat. Mater., 2009, 8, 194–197 CrossRef CAS.
- M. Evangelisti and E. K. Brechin, Dalton Trans., 2010, 39, 4672 RSC.
- S. Sanvito, Chem. Soc. Rev., 2011, 40, 3336–3355 RSC.
- M. F. Crommie, Science, 2005, 309, 1501–1502 CrossRef CAS.
- M. N. Leuenberger and D. Loss, Nature, 2001, 410, 789–793 CrossRef CAS.
- S. Thiele, F. Balestro, R. Ballou, S. Klyatskaya, M. Ruben and W. Wernsdorfer, Science, 2014, 344, 1135–1138 CrossRef CAS.
- J. M. Zadrozny, J. Niklas, O. G. Poluektov and D. E. Freedman, ACS Cent. Sci., 2015, 1, 488–492 CrossRef CAS.
- M. Shiddiq, D. Komijani, Y. Duan, A. Gaita-Ariño, E. Coronado and S. Hill, Nature, 2016, 531, 348–351 CrossRef CAS.
- W. P. Schleich, Nature, 2000, 403, 256–257 CrossRef CAS.
- R. Sessoli, ACS Cent. Sci., 2015, 1, 473–474 CrossRef CAS.
- K. Bader, D. Dengler, S. Lenz, B. Endeward, S.-D. Jiang, P. Neugebauer and J. van Slageren, Nat. Commun., 2014, 5, 5304 CrossRef CAS.
- F. Donati, S. Rusponi, S. Stepanow, C. Wackerlin, A. Singha, L. Persichetti, R. Baltic, K. Diller, F. Patthey, E. Fernandes, J. Dreiser, Ž. Šljivančanin, K. Kummer, C. Nistor, P. Gambardella and H. Brune, Science, 2016, 352, 318–321 CrossRef CAS.
- C. A. P. Goodwin, F. Ortu, D. Reta, N. F. Chilton and D. P. Mills, Nature, 2017, 548, 439–442 CrossRef CAS PubMed.
- F. S. Guo, B. M. Day, Y. C. Chen, M. L. Tong, A. Mansikkamäki and R. A. Layfield, Angew. Chem., Int. Ed., 2017, 56, 11445–11449 CrossRef CAS PubMed.
- S. K. Gupta, T. Rajeshkumar, G. Rajaraman and R. Murugavel, Chem. Sci., 2016, 7, 5181–5191 RSC.
- Y. Chen, F. Ma, X. Chen, B. Dong, K. Wang, S. Jiang, C. Wang, X. Chen, D. Qi, H. Sun, B. Wang, S. Gao and J. Jiang, Inorg. Chem., 2017, 56, 13889–13896 CrossRef CAS PubMed.
- R. Sessoli, Nature, 2017, 548, 400–401 CrossRef CAS PubMed.
- S. Cardona-Serra, L. Escalera-Moreno, J. J. Baldoví, A. Gaita-Ariño, J. M. Clemente-Juan and E. Coronado, J. Comput. Chem., 2016, 37, 1238–1244 CrossRef CAS PubMed.
- R. J. Blagg, L. Ungur, F. Tuna, J. Speak, P. Comar, D. Collison, W. Wernsdorfer, E. J. L. McInnes, L. F. Chibotaru and R. E. P. Winpenny, Nat. Chem., 2013, 5, 673–678 CrossRef CAS.
- J. J. Baldoví, S. Cardona-Serra, J. M. Clemente-Juan, E. Coronado, A. Gaita-Ariño and A. Palii, Inorg. Chem., 2012, 51, 12565–12574 CrossRef.
- N. F. Chilton, Inorg. Chem., 2015, 54, 2097–2099 CrossRef CAS PubMed.
- L. Ungur, J. J. Le Roy, I. Korobkov, M. Murugesu and L. F. Chibotaru, Angew. Chem., Int. Ed., 2014, 53, 4413–4417 CrossRef CAS.
- A. Gaita-Ariño, H. Prima-García, S. Cardona-Serra, L. Escalera-Moreno, L. E. Rosaleny and J. J. Baldoví, Inorg. Chem. Front., 2016, 3, 568–577 RSC.
- L. Thomas, F. Lionti, R. Ballou, D. Gatteschi, R. Sessoli and B. Barbara, Nature, 1996, 383, 145–147 CrossRef CAS.
- D. Gatteschi and R. Sessoli, Angew. Chem., Int. Ed., 2003, 42, 268–297 CrossRef CAS PubMed.
- R. Sessoli and A. K. Powell, Coord. Chem. Rev., 2009, 253, 2328–2341 CrossRef CAS.
- S. T. Liddle and J. van Slageren, Chem. Soc. Rev., 2015, 44, 6655–6669 RSC.
- S. Alvarez, P. Alemany, D. Casanova, J. Cirera, M. Llunell and D. Avnir, Coord. Chem. Rev., 2005, 249, 1693–1708 CrossRef CAS.
- M. Perfetti, E. Lucaccini, L. Sorace, J. P. Costes and R. Sessoli, Inorg. Chem., 2015, 54, 3090–3092 CrossRef CAS PubMed.
- H. Bethe, Ann. Phys., 1929, 395, 133–208 CrossRef.
- J. D. Rinehart and J. R. Long, Chem. Sci., 2011, 2, 2078 RSC.
- J. J. Baldoví, S. Cardona-Serra, J. M. Clemente-Juan, E. Coronado, A. Gaita-Ariño and A. Palii, J. Comput. Chem., 2013, 34, 1961–1967 CrossRef.
- J. J. Baldoví, J. J. Borrás-Almenar, J. M. Clemente-Juan, E. Coronado and A. Gaita-Ariño, Dalton Trans., 2012, 41, 13705 RSC.
- N. F. Chilton, D. Collison, E. J. L. McInnes, R. E. P. Winpenny and A. Soncini, Nat. Commun., 2013, 4, 2551 Search PubMed.
- L. F. Chibotaru and L. Ungur, J. Chem. Phys., 2012, 137, 64112 CrossRef CAS.
- F. Aquilante, J. Autschbach, R. K. Carlson, L. F. Chibotaru, M. G. Delcey, L. De Vico, I. Fdez. Galván, N. Ferré, L. M. Frutos, L. Gagliardi, M. Garavelli, A. Giussani, C. E. Hoyer, G. Li Manni, H. Lischka, D. Ma, P. Å. Malmqvist, T. Müller, A. Nenov, M. Olivucci, T. B. Pedersen, D. Peng, F. Plasser, B. Pritchard, M. Reiher, I. Rivalta, I. Schapiro, J. Segarra-Martí, M. Stenrup, D. G. Truhlar, L. Ungur, A. Valentini, S. Vancoillie, V. Veryazov, V. P. Vysotskiy, O. Weingart, F. Zapata and R. Lindh, J. Comput. Chem., 2016, 37, 506–541 CrossRef CAS PubMed.
- D. Aravena, M. Atanasov and F. Neese, Inorg. Chem., 2016, 55, 4457–4469 CrossRef CAS.
- D. Aravena and E. Ruiz, Inorg. Chem., 2013, 52, 13770–13778 CrossRef CAS PubMed.
- J. J. Baldoví, Y. Duan, R. Morales, A. Gaita-Ariño, E. Ruiz and E. Coronado, Chem.–Eur. J., 2016, 22, 13532–13539 CrossRef PubMed.
- J. van Leusen, M. Speldrich, H. Schilder and P. Kögerler, Coord. Chem. Rev., 2015, 289–290, 137–148 CrossRef CAS.
- J. M. Zadrozny, D. J. Xiao, M. Atanasov, G. J. Long, F. Grandjean, F. Neese and J. R. Long, Nat. Chem., 2013, 5, 577–581 CrossRef CAS PubMed.
- A. Lunghi, F. Totti, R. Sessoli and S. Sanvito, Nat. Commun., 2017, 8, 14620 CrossRef PubMed.
- K. S. Pedersen, J. Dreiser, H. Weihe, R. Sibille, H. V. Johannesen, M. A. Sørensen, B. E. Nielsen, M. Sigrist, H. Mutka, S. Rols, J. Bendix and S. Piligkos, Inorg. Chem., 2015, 54, 7600–7606 CrossRef CAS PubMed.
- K. N. Shrivastava, Phys. Rev., 1969, 187, 446–450 CrossRef CAS.
- K. N. Shrivastava, Phys. Lett. A, 1970, 31, 454–455 CrossRef CAS.
- K. N. Shrivastava, Chem. Phys. Lett., 1973, 22, 622–624 CrossRef CAS.
- A. Manoogian, J. Magn. Reson., 1979, 36, 1–6 CAS.
- K. N. Shrivastava, J. Phys. C: Solid State Phys., 1982, 15, 3869–3876 CrossRef CAS.
- R. L. Hartman, J. S. Bennett and J. G. Castle, Phys. Rev. B: Solid State, 1970, 1, 1946–1953 CrossRef.
- J. H. Van Vleck, Phys. Rev, 1940, 57, 426–447 CrossRef CAS.
- M. Dong-Ping, C. Ju-Rong and M. Ning, Commun. Theor. Phys., 2001, 36, 357–364 CrossRef CAS.
- M. Vanhaelst, P. Matthys and E. Boesman, Phys. Status Solidi A, 1978, 85, 639–644 CrossRef CAS.
- A. Singh and K. N. Shrivastava, Phys. Status Solidi A, 1980, 101, 51–56 CrossRef CAS.
- L. J. Berliner, G. R. Eaton and S. S. Eaton, Distance measurements in biological systems by EPR, Kluwer Academic, 2002 Search PubMed.
- K. N. Shrivastava, R. S. Rubins, S. Hutton and J. E. Drumheller, J. Chem. Phys., 1988, 88, 705–706 CrossRef CAS.
- L. Escalera-Moreno, N. Suaud, A. Gaita-Ariño and E. Coronado, 2015, arXiv:1512.05690v1.
- L. Escalera-Moreno, N. Suaud, A. Gaita-Ariño and E. Coronado, J. Phys. Chem. Lett., 2017, 8, 1695–1700 CrossRef CAS PubMed.
- K. N. Shrivastava, Chem. Phys. Lett., 1973, 20, 106–107 CrossRef.
- K. N. Shrivastava, Phys. Rep., 1975, 20, 137–227 CrossRef.
- J. S. M. Harvey, Chem. Phys. Lett., 1971, 10, 62–64 CrossRef CAS.
- L. E. Misiak and M. Subotowicz, Solid State Commun., 1991, 80, 761–766 CrossRef CAS.
- K. N. Shrivastava, G. V. Rubenacker, S. L. Hutton, J. E. Drumheller and R. S. Rubins, J. Chem. Phys., 1988, 88, 634–636 CrossRef CAS.
- R. S. Rubins and S. K. Jani, J. Chem. Phys., 1977, 66, 3297–3299 CrossRef CAS.
- A. M. Karo and J. R. Hardy, Phys. Rev., 1963, 129, 2024–2036 CrossRef CAS.
- M. Atzori, E. Morra, L. Tesi, A. Albino, M. Chiesa, L. Sorace and R. Sessoli, J. Am. Chem. Soc., 2016, 138, 11234–11244 CrossRef CAS PubMed.
- M. Atzori, L. Tesi, E. Morra, M. Chiesa, L. Sorace and R. Sessoli, J. Am. Chem. Soc., 2016, 138, 2154–2157 CrossRef CAS.
- S. K. Hoffmann and S. Lijewski, J. Magn. Reson., 2013, 227, 51–56 CrossRef CAS PubMed.
- S. K. Hoffmann and S. Lijewski, J. Magn. Reson., 2015, 252, 49–54 CrossRef CAS.
- W. M. Walsh, Phys. Rev., 1959, 114, 1473–1485 CrossRef CAS.
- A. J. Fielding, S. Fox, G. L. Millhauser, M. Chattopadhyay, P. M. H. Kroneck, G. Fritz, G. R. Eaton and S. S. Eaton, J. Magn. Reson., 2006, 179, 92–104 CrossRef CAS.
- K. L. Wan, S. L. Hutton, J. E. Drumheller and R. S. Rubins, J. Chem. Phys., 1987, 86, 3801–3803 CrossRef CAS.
- W.-C. Zheng and S.-Y. Wu, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 226, 409–412 CrossRef CAS.
- W.-C. Zheng and S.-Y. Wu, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 54, 1117–1122 CrossRef CAS.
- R. S. Rubins, J. Phys. Chem. Solids, 2005, 66, 1675–1682 CrossRef CAS.
- W.-C. Zheng, G.-M. Jia, L. He and W.-Q. Yang, Spectrochim. Acta, Part A, 2011, 78, 818–820 CrossRef PubMed.
- A. Lunghi, F. Totti, S. Sanvito and R. Sessoli, Chem. Sci., 2017, 8, 6051–6059 RSC.
- J. K. Bristow, J. M. Skelton, K. L. Svane, A. Walsh and J. D. Gale, Phys. Chem. Chem. Phys., 2016, 18, 29316–29329 RSC.
- S. Lenz, K. Bader, H. Bamberger and J. van Slageren, Chem. Commun., 2017, 53, 4477–4480 RSC.
- C. Wäckerlin, F. Donati, A. Singha, R. Baltic, S. Rusponi, K. Diller, F. Patthey, M. Pivetta, Y. Lan, S. Klyatskaya, M. Ruben, H. Brune and J. Dreiser, Adv. Mater., 2016, 28, 5195–5199 CrossRef PubMed.
- P. G. Klemens, Phys. Rev., 1962, 125, 1795–1798 CrossRef.
- H. Klein, U. Scherz, M. Schulz, H. Setyono and K. Wisznewski, Z. Phys. B: Condens. Matter Quanta, 1977, 28, 149–157 CrossRef CAS.
- A. Palii, S. Ostrovsky, O. Reu, B. Tsukerblat, S. Decurtins, S.-X. Liu and S. Klokishner, J. Chem. Phys., 2015, 143, 84502 CrossRef PubMed.

## Footnotes |

† These authors contributed equally to this work. |

‡ Phonons are lattice vibrations that can be imagined as particles carrying a quantum of vibrational energy. |

§ Using the definition of the Debye temperature (see Fig. 3), k_{B}Θ_{D} is equal to the phonon energy of maximum frequency. Thus, the Debye temperature can be interpreted as the temperature at which the highest-frequency vibration (hence, every one of them) is excited. Macroscopically, the Debye temperature can be regarded as a measure of the hardness of the crystal. Typical Debye temperatures range from 38 K for cesium to 2230 K for carbon. |

¶ Let us picture a vibration that propagates in the direction that is perpendicular to a given crystallographic plane. This will be a longitudinal phonon if the stretching and compression happen between successive planes, so that the geometric distortion is parallel to the direction of propagation. Conversely, it will be a transverse phonon if there is a lateral displacement between successive planes, so that the geometric distortion is perpendicular to the direction of propagation. In a one-dimensional solid, atoms are restricted to move along a given straight line, so phonons corresponded to longitudinal waves. In three-dimensional solids, atoms are not restricted anymore to the direction of propagation, and can also vibrate up and down, producing transverse waves. An effective sound velocity is commonly used to describe the speed propagation of a phonon, where one can distinguish a longitudinal and a transverse velocity, respectively. This effective sound velocity is also related to the hardness of the crystal. |

This journal is © The Royal Society of Chemistry 2018 |