Exploiting clock transitions for the chemical design of resilient molecular spin qubits

Abstract

Molecular spin qubits are chemical nanoobjects with promising applications that are so far hampered by the rapid loss of quantum information, a process known as decoherence. A strategy to improve this situation involves employing so-called Clock Transitions (CTs), which arise at anticrossings between spin energy levels. At CTs, the spin states are protected from magnetic noise and present an enhanced quantum coherence. Unfortunately, these optimal points are intrinsically hard to control since their transition energy cannot be tuned by an external magnetic field; moreover, their resilience towards geometric distortions has not yet been analyzed. Here we employ a python-based computational tool for the systematic theoretical analysis and chemical optimization of CTs. We compare three relevant case studies with increasingly complex ground states. First, we start with vanadium(IV)-based spin qubits, where the avoided crossings are controlled by hyperfine interaction and find that these S = 1/2 systems are very promising, in particular in the case of vanadyl complexes in an L-band pulsed EPR setup. Second, we proceed with a study of the effect of symmetry distortions in a holmium polyoxotungstate of formula [Ho(W₅O₁₈)₂]⁹⁻ where CTs had already been experimentally demonstrated. Here we determine the relative importance of the different structural distortions that causes the anticrossings. Third, we study the most complicated case, a polyoxopalladate cube [HoPd₁₂(AsPh)₈O₃₂]⁵⁻ which presents an unusually rich ground spin multiplet. This system allows us to find uniquely favorable CTs that could nevertheless be accessible with standard pulsed EPR equipment (X-band or Q-band) after a suitable chemical distortion to break the perfect cubic symmetry. Since anticrossings and CTs constitute a rich source of physical phenomena in very different kinds of quantum systems, the generalization of this study is expected to have impact not only in molecular spin science but also in other related fields such as molecular photophysics and photochemistry.

---

Introduction

Quantum two-level systems based on spin states, or simply “spin qubits”, are promising candidates as building blocks for quantum technologies.^1–3 The potential of a magnetic system for its application as a qubit relies on the capacity to generate, control and read out different quantum superpositions of two spin states (|0> and |1>). These requirements can be achieved within essentially any two-level quantum system, such as the simplest case of a spin S = 1/2. However, quantum information is fragile and is quickly lost for most quantum systems. The physical processes that lead to this information loss are collectively known as decoherence and can be grouped into two major sources: magnetic noise and vibrations. Different degrees of resilience against the two decoherence sources can be found in the many possible different spins systems.

Molecular magnetism offers an attractive approach to design spin qubits.^4,5 In molecular spin qubits, the spin typically resides on a magnetic metal ion, and chemical design of the coordination complex offers possibilities for a rational optimization of the physical properties, making them a promising pathway to achieve coherent qubits. The local coordination of the ion strongly influences the qubit properties by determining the wave functions of its |0> and |1> states. In comparison with solid-state approaches, coordination chemistry offers a tremendous range of possible environments and, therefore, of designing qubits with suitable characteristics, conferring a competitive advantage to magnetic molecules over other promising candidates such as NV centers in diamond^6–8 or phosphorus impurities in silicon.^2,9 At the same time, molecular spin qubits allow comparable survival time of quantum spin coherence T₂, the characteristic time of the exponential decay of spin echoes, typically associated to spin–spin relaxation processes.³ The optimal conditions to obtain the maximum values of T₂ for all solid state approaches include extreme dilution of the spin in a diamagnetic matrix. The reason for this requirement is the linear dependence of the transition energy between two states of a spin qubit as a function of the local field, due to the Zeeman effect, which means that qubits are strongly affected by the magnetic noise generated by neighboring spins.

Chemistry provides a promising method, complementary to dilution, to effectively “engineer” the molecular spin states to reduce this noise and thereby obtain longer quantum coherence times T₂. The origin of this quantum stability is the special wavefunction mixing happening when the spin states |↑> and |↓> experience an avoided level crossing, giving rise to a tunneling splitting Δ_CT.¹⁰ At these avoided crossings the two spin states present zero Zeeman slope, making the transition frequency Δ between these two states insensitive to small changes in the magnetic field. This results in optimal operating points, known as Clock Transitions (CTs) and are characterized by the magnitude of their tunneling splitting Δ_CT (see Fig. 1 up). Since this protection is only absolute at a particular field B_CT, one can compare the protection offered by CTs in different systems by quantifying the curvature at the anticrossing, or, equivalently, by the sensitivity to the magnetic field ∂Δ/∂B in the vicinity of the anticrossing (see Fig. 1, down).

Fig. 1 (Up) Anticrossings between pairs of states with equal Zeeman slopes dE/dB but either large (blue) or a small (black) tunneling splitting Δ_CT. (Down) Sensitivity to the magnetic field |∂Δ/∂B|, near the anticrossings: in the black curve the protection is lost quickly whereas a larger Δ_CT means a relatively wide window of protection for the blue curve.

However, working with CTs possesses particular challenges. In most EPR experiments, the resonance condition between the magnetic energy level splitting and the microwave frequency can always be met just by sweeping the magnetic field, achieving a continuous increase of the splitting due to the Zeeman effect. In contrast, each CT happens at a defined magnetic field B_CT and has a unique value of energy Δ_CT given by the spin Hamiltonian and which cannot be controlled by the Zeeman term. Since the microwave frequencies are typically difficult to control, this possesses an intrinsic limitation in experimentally characterizing and exploiting CTs. As a consequence, in order to exploit CTs one requires a certain control over the energy level structure; this can be achieved by designing either the crystal-field around the magnetic ion and/or the electro-nuclear hyperfine interactions.^9,11 The former merely requires to choose the appropriate combination of metal ion and coordination environment, but the hyperfine interaction may also be tuned in molecules, for example, through s–d orbital mixing.¹² Rational chemical design in principle offers a path forward, and indeed some of us demonstrated this in the [Ho(W₅O₁₈)₂]⁹⁻ polyoxometalate (in short, HoW₁₀ POM) that presents CTs because of a combination of the local symmetry and the nature of the ground doublet.¹³ These POMs can then be packed in the solid state at unusually high concentrations while retaining desirably long coherence times. Since this initial contribution, a few further works have employed the same kind of strategy to obtain coherence in molecular spin qubits at unusually high concentrations.^12,14,15 However, this has so far mostly been a matter of trial-and-error.

Generally speaking, CTs appear for every avoided crossing, which is to say in any and every magnetic molecule presenting extradiagonal elements in their spin Hamiltonians, independently of whether these arise from the crystal field or from the hyperfine coupling. From the theoretical point of view, if the form of the spin Hamiltonian is known, it is in principle possible to obtain either exact or approximate analytical solutions that describe at which magnetic field and excitation energy the CTs will appear, as well as their curvature. However, there is no systematic procedure for the exploration of chemical structures that allows the optimization of CTs. Furthermore, there are also perturbations in the molecular geometry or in the parameters of the spin Hamiltonian due for example to crystalline defects or to thermal molecular vibrations. Depending on the nature of the spin Hamiltonian, CTs can be more or less sensitive to this source of noise. Still, there is currently no developed methodology to deal with this problem.

In this work we characterize the behavior of anticrossings in molecular spin qubits. The robustness against magnetic noise is quantified via the curvature of the levels participating in the CT, while the robustness versus molecular distortions is quantified as the relative change in the transition frequency at reasonable degrees of distortions. We perform this kind of analysis for different representative kinds of spin Hamiltonians: from a simple S = 1/2 transition metal, well isolated from its excited states, to a lanthanide in cubic symmetry presenting a highly degenerate ground multiplet. Our final purpose here is to find some insights into the chemical design of “shallow” CTs where a spin qubit can be most resilient. For these goals, in this work we develop and employ a computational tool that assists in the automated exploration the influence of different molecular parameters.

Methods: algorithmic detection of avoided crossings

Herein we present the numerical method and the software implementation which we developed in this work and which is capable of getting an autonomous numerical characterization of the coordinates and curvature of any anticrossing in discrete {x, y₁, y₂, …, y_n}-type datasets. In the case of magnetic molecules this means that, starting from files containing Zeeman diagrams for a set of energy levels {B, E₁, E₂, …, E_n}, the program automatically obtains, for every possible CT, its magnetic field intensity (B_CT), its tunneling splitting (Δ_CT) and more crucially also the curvature at the CT, which allows to quantify the robustness of each CT vs. magnetic field noise.

The program starts by translating the set of discrete energy points into continuous functions, a task which involves automatically distinguishing between crossings and anticrossings (see Fig. 2 for illustration and ESI Section S1 for details†). It processes energy differences only up to a specified frequency threshold, which is defined by the user, e.g. 15 GHz (0.5 cm⁻¹). Within the chosen frequency window, the program then obtains for each anticrossing the second derivative ∂²Δ/∂B²-which is a very good approximation for the curvature (k) – (see ESI Section S1 for details†). A systematic application of the code to a data set obtained by diagonalization of a spin Hamiltonian with a controlled parameter variation allows a fast and systematic study of the influence of each parameter on the frequency and curvature of the anticrossings.

Fig. 2 Ordering of the energy levels between pre-processing (numerical data, symbols) and post-processing (analytical fit, lines) is the same in anticrossings but altered in crossings.

We have applied this methodology to study representative cases in order to deepen our understanding of the nature of the parameters governing the behavior. In some cases, we rely on the computational package SIMPRE to obtain the spin Hamiltonian parameters from a controlled variation in the coordination environment.^11,16 SIMPRE is based on the following Hamiltonian:^17,18

where k is the order and q the operator range which varies between −k and +k of the Stevens operator equivalents Ô_k^q as defined by Ryabov.¹⁹ The crystal field parameters (CF parameters) can be expressed as:

B_k^q = a_k(1 − σ_k)A_k^q〈r^k〉, (2)

where a_k = α, β, γ are the multiplicative Stevens coefficients for k = 2, 4, 6, σ_k are the Sternheimer shielding parameters of the f shell and 〈r^k〉 are the expectation values of the radii. The CF parameters expressed as A_k^q can then be estimated by an effective electrostatic model of N point charges around a rare earth ion, which parameterizes the electric field effect produced by the surrounding ligands using the following relations:where Z_i is the effective charge in units of e, the elementary charge, R_i, θ_i, ϕ_i are the effective spherical coordinates associated with the i-th donor atom, Z_kq are the tesseral harmonics and p_kq are the prefactors of the spherical harmonics.

Application to three case studies

We compared the robustness of CT, in terms of curvature and stability of the transition energy, for three metal-based complexes with three different representative kinds of spin Hamiltonians. We start with a vanadyl complex (a well isolated S = 1/2 ground doublet) where we study the influence of parallel and perpendicular hyperfine couplings (A_∥, A_⊥); a well-known example of this kind is provided by the [VO(C₃S₅)₂] complex (Fig. 3a and b).²⁰ We subsequently applied our methodology to study the HoW₁₀ POM where one can find a large tunneling splitting in the electronic ground state, Δ_CT = 9 GHz (0.3 cm⁻¹), due to the effect of the antiprismatic crystal field around the Ho³⁺, with the rest of the electronic spin states being higher in energy (Fig. 3c and d). For this second system we determine the variation of the CTs with regard to different molecular distortions. As a final case study, we focus on the polyoxopalladate [HoPd₁₂(AsPh)₈O₃₂]⁵⁻ (in short HoPd₁₂). This POM has an exceptionally complicated ground state as a result of the cubic coordination symmetry around the Ho³⁺, which leads to a spin doublet and a spin triplet in near degeneracy (Fig. 3e and f). The study of the variation of this complicated energy level scheme with respect to molecular distortions of the cube allows to gain some novel perspective on the relative robustness towards both magnetic noise and thermal noise in this system.

Fig. 3 Molecular structure of the three case studies, and simplified scheme of their electronic structures, emphasizing the multiplicity of their low-energy spin states and the gap to the rest of the spectrum: (a and b) [VO(C₃S₅)₂] (c and d) [Ho(W₅O₁₈)₂]⁹⁻, in short HoW₁₀ (e and f) [HoPd₁₂(AsPh)₈O₃₂]⁵⁻, in short HoPd₁₂.

Vanadyl dithiolate complexes

Behavior of the CTs

Vanadium(IV) complexes have been proven to be exceptionally promising as spin qubits, given their records in phase memory times, T₂ (ref. 21 and 22) and the existing proposals for quantum gates and spin filters based on them.^23–25 Beyond that, vanadium(IV) complexes can offer the advantage of studying a well isolated S = 1/2 system (see for example Fig. 3, left). This is ideal from the point of view of the theory as it serves as a perfect test bed for our computational tool since an analytical solution is available. We employed the following Hamiltonian:

Ĥ = μ_BB × g × Ŝ + Ŝ × A × Î (6)

where μ_B is the Bohr magneton, B is the magnetic field, g is the electronic Landé g-tensor, Ŝ is the electronic spin, Î is the nuclear spin and A is the hyperfine tensor. Owing to the characteristics of this Hamiltonian, the CTs will be fully determined by the hyperfine term. For all practical purposes, the nuclear Zeeman term is negligible compared to the electronic Zeeman term since μ_N ≪ μ_B. For vanadium(IV), we study the whole electro-nuclear spin manifold i.e. the states defined by: M_S = ±1/2, M_I = {±7/2, ±5/2, ±3/2, ±1/2}.

Note that vanadyl complexes and their vanadium analogues, as for example [VO(C_xS_y)₂]²⁻vs. [V(C_xS_x)₃]²⁻, tend to present very different A_∥/A_⊥ ratios:^3,20A_∥ is much higher in vanadyl complexes, while A_⊥ presents higher values in vanadium complexes. Moreover, the dispersion of the hyperfine parameters between different complexes tends to be lower for the vanadyl complexes, because of the major role played by the closely bonded oxygen in the spatial distribution of the vanadium d-electron density. Therefore, as case study we will focus vanadyl complexes.

For our study, we employed the averaged values of the four vanadyl dithiolates complexes referenced in ref. 20 and 21: ; . We varied them by ±15%, a reasonable range given their typically narrow parametric dispersion.²⁶ We limit the study to the variation of each one of the two parameters, keeping the other in its typical value. In simple anticrossings such as this, the curvature can be analytically estimated as k = γ_z²/2Δ, with γ_z = g × ΔM_S × μ_B, where g = 2 for vanadyl, and ΔM_S = 1. The numerical processing of the 658 calculated curves, of which we represent a selection in Fig. 4, coincides exactly with the analytical solutions (see ESI Section S2a†). In particular, one can see that the fields at which CTs appear are exclusively a linear function of A_∥ (Fig. 4a), while Δ and k are exclusively a linear function of A_⊥ (Fig. 4b). A deeper wavefunction analysis demonstrates that these CTs are only allowed in EPR parallel operating mode (see details in ESI Section S6†).

Fig. 4 (a) Zeeman diagrams (B_z) showing the effects of a ±15% variation in A_∥; the full graphs (see Fig. S7†) presents vertical and horizontal symmetry planes, but for clarity here the effects of the A_∥ are only shown for the bottom right part of the graph. Analytically-derived magnetic field positions (B_CT) of CTs 1–4 as a function of A_∥ are included as annotations. Curves in the same color correspond to the same spin levels along the variation. (b) Zeeman diagrams (B_z) showing the effects of −15% (red) and +15% (blue) variation in A_⊥, together with the analytically-derived transition energies Δ_CT of the different CTs. The analytically-derived curvatures (k) of the different CTs as a function of Δ_CT.

In our case example, CTs present the frequencies Δ₁ = 0.53 GHz (0.01768 cm⁻¹), Δ₂ = 0.51 GHz (0.01701 cm⁻¹), Δ₃ = 0.45 GHz (0.01501 cm⁻¹), Δ₄ = 0.35 GHz (0.01167 cm⁻¹), with curvatures in the range 25–35 cm⁻¹ T⁻². Table S1† summarizes the characteristic relevant parameters of this compound; more details can be found in the ESI Section S2.†

Although the tunneling splittings, Δ_i, in this kind of complexes are relatively small, their curvatures k_i, are remarkably low due to the low value of S. A further advantage of CTs derived from hyperfine coupling is the fact that for a given vanadyl complex there are CTs at four different transition frequencies, with the widest tunneling splitting being 50% larger than the narrowest. Since the values of the hyperfine coupling can vary from complex to complex, more opportunities for obtaining a resonance at the appropriate microwave source (L-band, 0.8–1.4 GHz, ∼0.03–0.045 cm⁻¹) are possible,^22,27,28 as well as novel possibilities to operate a multi-qubit system by applying local magnetic fields.²⁹

Influence of hyperfine coupling on the CTs

Although estimating the variation of A with distortions is computationally expensive for heavy atoms,³⁰ the effect of a parameter strain on these CTs can be estimated from Fig. 4. A relatively large deviation of 15% in the value of either A_⊥ or A_∥ has an effect on the transition frequency that is equivalent to a magnetic noise of about 10 mT, meaning the resilience towards magnetic noise is expected to be in tandem with a resilience towards vibrational noise. Therefore, high resilience is expected to be obtained using this hyperfine-driven approach. The price to pay for this is the equipment needed to access the very low energies that characterize these CTs, which will make them more challenging to detect and use. The generalization of this case study would be relevant to other hyperfine-based spin qubits such as the flip-flop qubit scheme based on ³¹P impurities in a ²⁸Si Matrix.³¹

Antiprismatic holmium POM complex

Behavior of the CTs

The complex HoW₁₀ behaves as a molecular spin qubit where, as mentioned above, certain quantum resonances were demonstrated to be exceptionally robust against a highly concentrated electronic spin bath (coherence times T₂ = 8 μs at 1% concentration in an isostructural YW₁₀ matrix).¹³ This behavior is due to the large tunneling splitting in the ground multiplet anticrossings (Fig. 5a), which ultimately can be related with the molecular structure (Fig. 3c) and its distortions (Fig. 5b). This tunneling splitting has been experimentally determined to be Δ_CT = 9 GHz (0.3 cm⁻¹).³² A recent effort has been made to rationalize the limit of the T₂ divergence at the CTs in terms of the nuclear spin bath,³³ but a study in terms of molecular structure is still lacking.

Fig. 5 (a) Zeeman diagrams (B_z) for the ground multiplet of HoW₁₀ with the CF parameters reported in the ESI Section S3† (black lines) and with B₄⁴ parameters that are 2σ larger (blue lines) or 2σ smaller (red lines), see text for details. Energy origin has been chosen at the center of the anticrossing for convenience. (b) Coordination sphere of HoW₁₀, emphasizing the deviation from the ideal D_4d symmetry: a skew angle ϕ ≠ 45° and the off-center position of Ho h ≠ h′. (c) Evolution of the calculated tunneling splitting Δ_CT for rising values of angular distortion α. The average crystallographic value is marked with a vertical line. (d) Evolution of Δ_CT for rising values of the vertical displacement d. The averaged crystallographic value is marked with a vertical line and the splitting determined by EPR is marked with a horizontal line.

To describe this system, one needs to extend the Hamiltonian in eqn (6) by including Crystal Field parameters:

where B_k^q, Ô_k^q are, respectively, the crystal field (CF) parameters and CF operators in the convention of the extended Stevens operators.^17,18 The double summation parameterizes the CF interaction, which for an ideal D_4d symmetry would contain only the axial terms B₂⁰Ô₂⁰, B₄⁰Ô₄⁰ and B₆⁰Ô₆⁰. Additionally, the effects of deviations from exact symmetry were recovered in previous works by introducing an extradiagonal B₄⁴Ô₄⁴ term, compatible with a C₄ symmetry. We started by employing our computational tool to automatically find the CTs that had been characterized experimentally, using pulsed EPR X-band employing the reported parameters (see Fig. 5a).³⁴ We studied the influence of the CF parameters on the anticrossings, both in terms of energy and curvature. Finally, we studied the variation of the CTs in terms of the main structural parameters: the skew angle ϕ and the displacement of the Ho³⁺ atom with respect to the center of the antiprism, (h − h′)/2 (see Fig. 5b).

In terms of the analysis of the CF parameters, we focused on the B₄⁴ parameter, since the M_J = ±4 ground doublet together with the C₄ coordination symmetry in HoW₁₀ mean that Ô₄⁴ (or, equivalently, Ô₆⁴) is the extradiagonal operator that governs the anticrossings in this system (see ESI Section S3†). We used the reported average value of B₄⁴ = 3.14 × 10⁻³ cm⁻¹; as a variation range, ΔB₄⁴ = B₄⁴ ± 2σ, we chose twice the empirical width for the Gaussian distribution responsible for the observed inhomogeneous broadening.¹³

We solved the system for increasing values of B₄⁴ in 43 steps and processed each case automatically with our script to extract magnetic fields, B_CT, transition energies, Δ_i, and curvatures, k_i, of each of the 43 × 8 = 344 anticrossings in the desired energy window (in this case, E < 15 GHz ≈ 0.5 cm⁻¹). Our numerical determination of Δ_i and k_i at the anticrossings here serves as a confirmation of an approximate but excellent analytical solution of a problem that was already discussed in ref. 13 and that was also employed in the V(IV) above: the curvature is analytically estimated as k = γ_z²/2Δ, with γ_z = g × ΔM_J × μ_B, where g is 1.25 for Ho³⁺ and ΔM_J is equal to 8. These parameters are summarized in Table 1.

Table 1 Relevant parameters associated with the CTs of the three case studies. Note that *¹ is referred to energy levels appearing in Fig. 6c colored in orange, and *² for the levels colored in blue

System	Δ CT (tunneling splitting) GHz	k (curvature) cm^−1 T^−2	g	ΔMJ	Ref.
VO(C3S5)2	0.35–0.53	∼24–37	2	1	20
HoW10	9	∼36	1.25	8	13 and 32
HoPd12	0–70	(*^1) ∼ 0, (*^2) ∼ 70	1.25	N.A.	39 and 40

---

Control of CTs by crystal field

Both the transition energies and the curvatures of the anticrossings vary with the extradiagonal parameter B₄⁴. From the point of view of the chemical structure, the presence of any extradiagonal parameters in this structure is related with the distortion of the coordination sphere with respect of a perfect D_4d. This has commonly been attributed, in absence of calculations, to the C₄-preserving skewing distortion characterized by an angle ϕ (see Fig. 5b). This skew angle has been structurally characterized in at least four POM families with LnO₈ similar coordination environments: LnW₁₀, [Ln(β₂-SiW₁₁O₃₉)₂]¹³⁻ (ref. 35) [Ln(β-Mo₈O₂₆)₂]⁵⁻ and [Ln{Mo₅O₁₃(OCH₃)₄NNC₆H₄-p-NO₂}₂]³⁻.³⁶ Depending on the Ln ion, skew angles in the range 39° < ϕ < 44.5° have been found but no rationalization has been done in terms of the molecular structure.

In general, the dependence of the extradiagonal parameters with the distortions are neither linear nor simple. To test the most commonly invoked hypothesis, we studied an idealized version of HoW₁₀ where the coordination sphere is a perfect square antiprism except for the skew angle ϕ ≠ 45°. The Δ vs. ϕ dependence shows that the experimental average skew angle ϕ = 44.2° in HoW₁₀ only accounts for a small fraction of the experimental tunneling splitting (Fig. 5c). Undoubtedly, other molecular distortions play a significant role in the appearance of the large Δ in this case. Among these are surely the off-center position of Ho³⁺, which is 2d = h − h′ = 0.05 Å closer to one of the [W₅O₁₈]⁶⁻ moieties, and the non-coplanarity of the planes containing the two O₄ squares, which in this case form an angle β = 1.11°.

We explored the influence of the off-center vertical displacement in the position of Ho³⁺, expressed as d (see Fig. 5b). We found a much stronger effect of this experimental distortion in the crystal structure compared with the torsion angle. The experimental splitting, is completely recovered by the combination of the experimental (averaged) skew angle ϕ = 44.2 and d = 0.0225 Å, remarkably close to the experimental (averaged) value d = 0.025 Å (see Fig. 5d), with the latter having the responsibility for most of the tunneling splitting and the former playing a comparatively minor role. This is expected to change from case to case along the different series of near-D_4d POMs, given the different distortions that they present (accounted by ϕ and d values).

Of course, to understand the behavior beyond this static picture that employ averaged crystallographic positions at high temperature, one would need to consider the thermal displacements due to vibrations at any finite temperatures. The degree of structural variation that one might expect from vibrations in the solid state, and in particular the magnitude of the real-time distortions of Ln-ligand distances in comparatively soft lanthanide complexes were estimated recently by some of us from the study of sub-picosecond dynamics of spin energy levels in magnetic biomolecules.³⁷ This allows an upper bound to the instantaneous distortions expected at low temperature to be extracted: instantaneous distortions of up to 0.01–0.02 Å in the coordination bond length are calculated in the femtosecond time scale, with averaged distortions decreasing with the averaging time following a square root law (see ESI Section S5 for details†).

While a complete analysis of the static and dynamic distortions in the Ho³⁺ coordination sphere HoO₈ in the different series of POMs is outside the scope of the present work, the theoretical framework presented herein, combined with a recently developed theoretical methodology,³⁸ would allow the systematic analysis and rationalization of the resilience of CTs of the different POMs against thermal noise.

Let us rationalize the comparison between anticrossings in HoW₁₀ and in vanadyl complexes. The differences in Δ_CT, g and in M_J result in curvatures for HoW₁₀ of approximately 35 cm⁻¹ T⁻², meaning the protection against magnetic noise is comparable or worse than in the vanadyl case. In terms of resilience towards vibrational noise, one can consider that HoW₁₀ molecules which in the timescale of the pulsed EPR experiment present an extradiagonal parameter that deviates from the average value of B₄⁴ present a larger (or smaller) transition frequency; this effect is comparable to a local magnetic field of over 25 mT, again similar or worse than in the vanadyl case. In terms of molecular approaches to CTs, HoW₁₀ serves here as an example of CT based on the crystal field in lanthanide ions. In comparison with the strategy of employing S = 1/2 systems with CTs based on hyperfine coupling such as vanadyl systems, employing lanthanide ions offer the possibility of the chemical design of molecules with a certain control over crystal field, as a way of tuning Δ_CT. In the case of HoW₁₀, Δ_CT = 9 GHz, allowing to access the CTs with commercial X-band equipment.

Cubic holmium POM complex

Behavior of the CTs

As a last case study case we chose a Ho³⁺ ion in a cubic or near-cubic environment, where the high symmetry can produce a highly degenerate low-energy spectrum with 40 electronuclear spin levels within a ∼2 cm⁻¹ window. This is the case of a cubic, in short, HoPd₁₂ POM (Fig. 3e and f).^39,40 On a preliminary study on this system, we employed a minimalistic theoretical treatment that neglected the nuclear spin to suggest that the avoided crossings could result in exceptionally high coherence times.⁴¹ Herein, we extend that simplified Hamiltonian to include the hyperfine interaction (eqn (7)).

What in the original study appeared as a ground spin doublet and a closely lying triplet (Fig. 3f) is after consideration of hyperfine coupling, which we assumed equal to the value for A = 0.83 GHz (0.0277 cm⁻¹), transformed into a much more complicated structure. We have now five manifolds displaying an intricated set of anticrossings (see ESI for a detailed analysis†). There are of course an enormous number of possible CTs in this system, of which we have selected just two for illustration. In absence of distortion, there is a promising region between the third and the fourth manifolds, where transitions can be found near the X-band energy window and with rather small curvature, which result from the competition between the effective energy level repulsion exerted by the upper and lower manifolds.

Control of CTs by crystal field

Starting from this energy level scheme, the application of a moderate structural axial compression was studied in order to estimate the variation on the relative energies upon distortion of the exact cubic symmetry (Fig. 6c). Note that any possible distortion will create a comparable effect, and in this respect note the recent experimental results on quantum coherent spin-electric control in molecular nanomagnets, where a quantum coherent manipulation of spin states was achieved relying on minute distortions created by an external electric field.⁵⁹ Compression rate is defined as c = ((h₀ − h_c)/h₀) × 100. Calculations up to c = 1% show that the main effect of increasing the axial compression is to rise the energy of the “light blue” manifold (see Fig. 6d–f for three snapshots of the energy level scheme at different compressions).

Fig. 6 (a) Coordination sphere of HoPd₁₂, emphasizing the deviation we create from the ideal O_h cubic symmetry via an axial compression c. The distance h_c between the Ho³⁺ ion and the upper (or lower) O₄ planes is progressively decreased from its ideal value h_c = h₀ (c = 0%) down to h_c = 0.99 × h₀ (c = 1%). (b) Zeeman diagrams (B_z) for the low energy region with the average crystal field parameters. Each color corresponds to an octuplet with consistent Zeeman behavior at high fields. For convenience, the energy origin has been chosen at the center of the fundamental multiplets in the non-distorted cube. (c) Zoom of the upper and lower region and the avoided crossings with Δ in the range: (up) 9–30 GHz (0.3–1 cm⁻¹) and (down) 0.8–15 GHz (0.03–0.5 cm⁻¹). Minimum and maximum value of the curvature of the highlighted CTs, in cm⁻¹ T⁻²: (up) 0.0–70 and (down) 2.5–26. (d) 0.06% of compression (e) 0.65% of compression (f) 0.95% of compression. More details in ESI Section S4b.†

While this cubic system is in many aspects analogous to HoW₁₀, the novel ingredient in this case is the high degeneracy in the spin states, which gives rise to a competition of anticrossings. Our calculations points towards a very high probability of finding suitable CTs in this compound, both in the ideal cubic geometry and also if the experimental reality of the crystallization and/or cooling processes create minor distortions.

Among the large number and variety of CTs that can be found on this system, we can highlight the transition between levels 5 and 14 at a field 0.0247 T, with a transition energy Δ = 15 GHz (0.50 cm⁻¹). See Fig. 6c up. Not only are the two levels are very flat, with curvatures k₅ = 0.82386 and k₁₄ = 0.39665: the curvature of both are in the same direction, meaning the second derivative of the transition energy is the difference – instead of the addition – of their curvatures.

Moreover, comparing HoPd₁₂ with previous systems, we conclude that HoPd₁₂ offers a unique protection from magnetic noise but an extreme sensitivity to molecular distortions. The reason for this is that the magnitude of Δ_CT in some of the CTs is not solely determined by a tunneling splitting, but rather by the energy difference between two different manifolds. In other words, the change in the diagonal CF terms directly affect the transition energy. Even in absence of dynamical calculations, this is likely to mean a high thermal effect, since any distortion will dramatically affect the transition energies, precisely at the CTs. Comparison with previous studies on the time-averaged effects of vibrations on spin energy levels⁴⁸ allows us to expect important spectral diffusion effects, which should be addressed experimentally by hole-burning. A way of looking at the difference between HoPd₁₂ and HoW₁₀ in the control requirements for the electrical manipulation of CTs is by comparing the effect of a vertical compression, c, for HoPd₁₂ to the effect of a d in HoW₁₀. Linear interpolating from Fig. 5d, shows that a vertical displacement of the Ho atom of the order of 10⁻⁶ Å results in a variation of Δ of about 1 MHz, so in the range of what is achieved in the experiments on HoW₁₀via an electrical field.⁵⁹ A comparable compression onto HoPd₁₂ results, for transition 21–30 (Δ = 13.5 GHz = 0.45 cm⁻¹, k = 5.2 cm⁻¹ T⁻²) (see ESI Section S4a†), in a change of about 200 MHz (see ESI Section S4a†): an effect that is stronger by two orders of magnitude. This means that the extraordinary sensitivity of HoPd₁₂ against molecular distortions creates a unique situation also in terms of electrical control of spin states at the CT. Note that to validate these predictions we are still missing an experimental EPR spectrum of HoPd₁₂, which is predicted to present an extremely complex structure. Most importantly, all relevant transitions are occurring from excited levels, which will be depopulated at low temperature. As an alternative, if the temperature is high enough, T₂ will be limited in practice by a relatively fast T₁, thus limiting the actual applicability of this system.

Conclusions

Anticrossings between quantum states are the source of different kinds of attractive physical phenomena, including clock transitions. At the same time, these optimal working points cannot be externally tuned by an external magnetic field and need to be tuned by chemical means. In this work we have addressed this problem theoretically by developing a python-based computational tool for the systematic analysis and chemical optimization of CTs. We employed the code to study three representative systems from lower to higher complexity on its electronic structure and Hamiltonian.

To analyze the influence of transverse hyperfine coupling in simple S = 1/2 systems, we started by studying CTs in V(IV) complexes. More precisely, in vanadyl complexes the hyperfine coupling is relatively robust against distortions from the coordination environment, a good protection against magnetic noise can be obtained in these systems. Nevertheless, these CTs present very low energies that make them difficult to access experimentally.

Next, by studying HoW₁₀, a paradigmatic system where CTs arise from the ligand field, we numerically rationalized the variation of both the transition frequency and the curvature of CTs with different kinds of molecular distortions from the ideal D_4d symmetry and give insights about the influence of the distortions that originates the anticrossings. Compared with the vanadyl case, we found a large sensitivity towards molecular distortions, since the energy levels are controlled by the crystal field, and a similar protection against magnetic noise. In addition to CTs with quantum tunneling splitting proportional to a most common frequency EPR band, such as X-band. Our final case study was a near-cubic HoPd₁₂ polyoxopalladate, which served to illustrate the cubic high-symmetry limit in which the ground state determined by the CF contains a large admixture of energy levels. In this regime we found new kinds of CTs resulting from different combinations of anticrossing levels, that result in a unique sensitivity to molecular distortions, offering a large range of tunneling splitting frequencies; and an equally unique protection from magnetic noise.

Note that CTs do not only appear in the context of molecular spin qubits. Rather, they are a general kind of quantum transitions that are uniquely robust against environmental noise. CTs are known since their use in the first cesium atomic standards of frequency and time interval,⁴² which gave rise to similar effects in atomic traps or atomic lattices,^43–46 and to the so-called ZEFOZs (for ‘zero first order Zeeman shift’).^47,48 Unique physical effects have been reported in very different systems and fields and employing a variety of terms, including “avoided level crossings”,^17,49 “avoided resonance crossings”⁵⁰ and “Landau-Zener crossings”.^51–55 The method developed and employed in this work is general and therefore can be used to study the evolution of the energy levels and crossings of a quantum system when an external physical stimulus is applied (magnetic field, electric field, pressure, etc.). It can equally be used to automatically explore physical or chemical modifications of any kind of systems presenting anticrossings, independently of the nature of the noise. The model employed herein can be applied to facilitate the detection of anticrossings in a number of fields, from the theoretical molecular design of strong coupling between plasmons and single-molecule excitons^56,57 to the prediction of the electronic and optical properties of twisted bilayer graphene.⁵⁸

As a final remark, the present work serves as a first numerical estimate of the sensitivity to structural or vibrational distortions of different kinds of CTs. Note it is possible to adjust molecular geometries by either physical pressure⁵⁹ or chemically.^60,61 Furthermore, recent experimental progress in the coherent control of spin states employing electrical fields,⁶² where the control is mediated via minimal molecular distortions, are another possible avenue for the applicability of the present methodology. On the other hand, if coupled with recently developed theoretical methodology,³⁸ this could serve to rationalize the influence of molecular vibrations in the effectiveness of CTs, and, eventually, to the chemical engineering of molecular spin qubits where the electron spin is not just shielded from magnetic noise but at the same time also protected from molecular vibrations.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

The present work has been funded by the EU (ERC-CoG 647301 DECRESIM, ERC-AdG-788222-Mol2D, COST 15128 Molecular Spintronics Project and the EU-QUANTERA project SUMO) the Spanish MINECO program (grants MAT2017-89993 and CTQ2017-89528 cofinanced by FEDER (EQC2018-004888-P), and Excellence Unit María de Maeztu MDM-2015-0538), and the Generalitat Valenciana (PROMETEU/2017/066, PROMETEU/2019/066) and PO FEDER (IDIFEDER/2018/061). S. C. S. thanks the Spanish MINECO for a Juan de la Cierva Incorporación postdoctoral Fellowship. S. G. S acknowledges the Ministry of Education of Spain (grant PRE2018-083350).

References

1. R. Hanson, V. V. Dobrovitski, A. E. Feiguin, O. Gywat and D. D. Awschalom, Science, 2008, 320, 352.
2. M. Steger, K. Saeedi, M. L. W. Thewalt, J. J. L. Morton, H. Riemann, N. V. Abrosimov, P. Becker and H.-J. Pohl, Science, 2012, 336, 1280.
3. J. M. Zadrozny, J. Niklas, O. G. Poluektov and D. E. Freedman, ACS Cent. Sci., 2015, 1, 488.
4. A. Gaita-Ariño, F. Luis, S. Hill and E. Coronado, Nat. Chem., 2019, 11, 301.
5. E. Coronado, Nat. Rev. Mater., 2020, 5, 87.
6. R. Hanson, V. V. Dobrovitski, A. E. Feiguin, O. Gywat and D. D. Awschalom, Science, 2008, 320, 352.
7. D. Scarabelli, M. Trusheim, O. Gaathon, D. Englund and S. J. Wind, Nano Lett., 2016, 16, 4982.
8. M. W. Doherty, N. B. Manson, P. Delaney, F. Jelezko, J. Wrachtrup and L. C. Hollenberg, Phys. Rep., 2013, 528, 1.
9. J. T. Muhonen, J. P. Dehollain, A. Laucht, F. E. Hudson, R. Kalra, T. Sekiguchi and A. Morello, Nat. Nanotechnol., 2014, 9, 986.
10. C. Zener, Proceedings of the Royal Society of London. Series A, Containing Papers of a Mathematical and Physical Character, 1932, 137, 696.
11. J. J. Baldoví, S. Cardona-Serra, J. M. Clemente-Juan, E. Coronado, A. Gaita-Ariño and A. Palii, Inorg. Chem., 2012, 51, 12565.
12. J. M. Zadrozny, A. T. Gallagher, T. D. Harris and D. E. Freedman, J. Am. Chem. Soc., 2017, 139, 7089.
13. M. Shiddiq, D. Komijani, Y. Duan, A. Gaita-Ariño, E. Coronado and S. Hill, Nature, 2016, 531, 348.
14. C. A. Collett, K. I. Ellers, N. Russo, K. R. Kittilstved, G. A. Timco, R. E. Winpenny and J. R. Friedman, Magnetochemistry, 2019, 5, 4.
15. R. T. Harding, S. Zhou, J. Zhou, T. Lindvall, W. K. Myers, A. Ardavan, G. A. D. Briggs, K. Porfyrakis and E. A. Laird, Phys. Rev. Lett., 2017, 119, 140801.
16. J. J. Baldoví, J. J. Borrás-Almenar, J. M. Clemente-Juan, E. Coronado and A. Gaita-Ariño, Dalton Trans., 2012, 41, 13705.
17. 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.
18. J. J. Baldoví, J. M. Clemente-Juan, E. Coronado, A. Gaita-Ariño and A. Palii, J. Comput. Chem., 2014, 35, 1930.
19. I. D. Ryabov, J. Magn. Reson., 1999, 140, 141.
20. C. J. Yu, M. J. Graham, J. M. Zadrozny, J. Niklas, M. D. Krzyaniak, M. R. Wasielewski, O. G. Poluektov and D. E. Freedman, J. Am. Chem. Soc., 2016, 138, 14678.
21. M. Atzori, L. Tesi, E. Morra, M. Chiesa, L. Sorace and R. Sessoli, J. Am. Chem. Soc., 2016, 138, 2154.
22. J. M. Zadrozny, J. Niklas, O. G. Poluektov and D. E. Freedman, ACS Cent. Sci., 2015, 1, 488.
23. M. Atzori, A. Chiesa, E. Morra, M. Chiesa, L. Sorace, S. Carretta and R. Sessoli, Chem. Sci., 2018, 9, 6183.
24. S. Cardona-Serra, A. Gaita-Ariño, E. Navarro-Moratalla and S. Sanvito, J. Phys. Chem. C, 2018, 122, 6417.
25. S. Cardona-Serra, A. Gaita-Ariño, M. Stamenova and S. Sanvito, J. Phys. Chem. Lett., 2017, 8, 3056.
26. M. Atzori, L. Tesi, E. Morra, M. Chiesa, L. Sorace and R. Sessoli, J. Am. Chem. Soc., 2016, 138, 2154.
27. C. V. Grant, W. Cope, J. A. Ball, G. G. Maresch, B. J. Gaffney, W. Fink and R. D. Britt, J. Phys. Chem. B, 1999, 103, 10627.
28. P. Pietrzyk and Z. Sojka, Appl. Magn. Reson., 2011, 40, 471.
29. M. D. Jenkins, D. Zueco, O. Roubeau, G. Aromí, J. Majer and F. Luis, Dalton Trans., 2016, 45, 16682.
30. F. Neese, Spin-Hamiltonian parameters from first principle calculations: theory and application, in High Resolution EPR, Springer, New York, 2009, pp. 175–229.
31. G. Tosi, F. A. Mohiyaddin, V. Schmitt, S. Tenberg, R. Rahman, G. Klimeck and A. Morello, Nat. Commun., 2017, 8, 1.
32. S. Ghosh, S. Datta, L. Friend, S. Cardona-Serra, A. Gaita-Ariño, E. Coronado and S. Hill, Dalton Trans., 2012, 41, 13697.
33. L. Escalera-Moreno, A. Gaita-Ariño and E. Coronado, Phys. Rev. B, 2019, 100, 064405.
34. S. Cardona-Serra and A. Gaita-Ariño, Dalton Trans., 2018, 47, 5533.
35. J. J. Baldoví, J. M. Clemente-Juan, E. Coronado, Y. Duan, A. Gaita-Ariño and C. Giménez-Saiz, Inorg. Chem., 2014, 53, 9976.
36. J. J. Baldoví, Y. Duan, C. Bustos, S. Cardona-Serra, P. Gouzerh, R. Villanneau, G. Gontard, J. M. Clemente-Juan, A. Gaita-Ariño, C. Giménez-Saiz and A. Proust, Dalton Trans., 2016, 45, 16653.
37. L. E. Rosaleny, K. Zinovjev, I. Tuñón and A. Gaita-Ariño, Phys. Chem. Chem. Phys., 2019, 21, 10908.
38. A. Ullah, J. Cerdá, J. J. Baldoví, S. Varganov, J. Aragó and A. Gaita-Ariño, J. Phys. Chem. Lett., 2019, 10(24), 7678–7683.
39. E. V. Chubarova, M. H. Dickman, B. Keita, L. Nadjo, F. Miserque, M. Mifsud, I. W. C. E. Arends and U. Kortz, Angew. Chem., Int. Ed., 2008, 47, 9542.
40. N. V. Izarova, M. T. Pope and U. Kortz, Angew. Chem., Int. Ed., 2012, 51, 9492.
41. J. J. Baldoví, L. E. Rosaleny, V. Ramachandran, J. Christian, N. S. Dalal, J. M. Clemente-Juan, P. Yang, U. Kortz, A. Gaita-Ariño and E. Coronado, Inorg. Chem. Front., 2015, 2, 893.
42. L. Essen and J. V. Parry, Nature, 1955, 176, 280.
43. M. Takamoto, F. L. Hong, R. Higashi and H. Katori, Nature, 2005, 435, 321.
44. T. Rosenband, P. O. Schmidt, D. B. Hume, W. M. Itano, T. M. Fortier, J. E. Stalnaker, K. Kim, S. A. Diddams, J. C. J. Koelemeij, J. C. Bergquist and D. J. Wineland, Phys. Rev. Lett., 2007, 98, 220801.
45. A. D. Ludlow, M. M. Boyd, T. Zelevinsky, S. M. Foreman, S. Blatt, M. Notcutt, T. Ido and J. Ye, Phys. Rev. Lett., 2006, 96, 033003.
46. A. D. Ludlow, T. Zelevinsky, G. K. Campbell, S. Blatt, M. M. Boyd, M. H. de Miranda, M. J. Martin, J. W. Thomsen, S. M. Foreman, J. Ye, T. M. Fortier, J. E. Stalnaker, S. A. Diddams, Y. Le Coq, Z. W. Barber, N. Poli, N. D. Lemke, K. M. Beck and C. W. Oates, Science, 2008, 319, 1805.
47. G. Heinze, C. Hubrich and T. Halfmann, Phys. Rev. Lett., 2013, 111, 033601.
48. M. Zhong, M. P. Hedges, R. L. Ahlefeldt, J. G. Bartholomew, S. E. Beavan, S. M. Wittig, J. J. Longdell and M. J. Sellars, Nature, 2015, 517, 177.
49. J. P. Pekola, D. S. Golubev and D. V. Averin, Phys. Rev. B, 2016, 93, 024501.
50. J. Wiersig, Phys. Rev. Lett., 2006, 97, 253901.
51. S. N. Shevchenko, S. Ashhab and F. Nori, Phys. Rep., 2010, 492, 1.
52. M. Urdampilleta, S. Klyatskaya, M. Ruben and W. Wernsdorfer, Phys. Rev. B: Condens. Matter Mater. Phys., 2013, 87, 195412.
53. G. Sun, X. Wen, M. Gong, D. W. Zhang, Y. Yu, S. L. Zhu, J. Chen, P. Wu and S. Han, Sci. Rep., 2015, 5, 8463.
54. A. Russomanno, A. Silva and G. E. Santoro, Phys. Rev. Lett., 2012, 109, 257201.
55. A. Palii, B. Tsukerblat, J. M. Clemente-Juan, A. Gaita-Ariño and E. Coronado, Phys. Rev. B: Condens. Matter Mater. Phys., 2011, 84, 184426.
56. R. Chikkaraddy, B. De Nijs, F. Benz, S. J. Barrow, O. A. Scherman, E. Rosta, A. Demetriadou, P. Fox, H. Ortwin and J. J. Baumberg, Nature, 2016, 535, 127.
57. T. Itoh, Y. S. Yamamoto and T. Okamoto, Phys. Rev. B, 2019, 99, 235409.
58. H. Patel, L. Huang, C. J. Kim, J. Park and M. W. Graham, Nature, 2019, 10, 1.
59. J. Hao, F. Hu, J. T. Wang, F. R. Shen, Z. Yu, H. Zhou, H. Wu, Q. Huang, K. Qiao, J. Wang, J. He, L. He, J. R. Sun and B. Shen, Chem. Mater., 2020, 32, 1807.
60. Y. Mizuguchi, E. Paris, T. Sugimoto, A. Iadecola, J. Kajitani, O. Miura, T. Mizokawa and N. L. Saini, Phys. Chem. Chem. Phys., 2015, 17, 22090.
61. H. Kamebuchi, A. Nakamoto, T. Yokoyama and N. Kojima, Bull. Chem. Soc. Jpn., 2015, 88, 419.
62. J. Liu, J. Mrozek, Y. Duan, A. Ullah, J. J. Baldoví, E. Coronado, A. Gaita-Ariño and A. Ardavan, 2020, arXiv:2005.01029 [cond-mat.mes-hall].

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: 10.1039/d0sc01187h

---