µSQUID-EPR mapping of tunnelling gaps
The multilevel character of Et4N[160GdPc2] (or [160GdPc2]−) (S = 7/2, L = 0, I = 0), the sizable anisotropy, and the crystal packing, make this system an excellent test bed for QPI effects. In contrast to 3d-based MMs, where QPI is observable via µSQUID studies, the exceedingly large tunnelling gaps in a 4f-MM require a different approach. We, henceforth, exploit the resonant absorption observed in the field-orientation dependent M(H) curves upon microwave absorption in the frequency range of ν = 0.1–20 GHz employing µSQUID-EPR technique (Fig. 1a). For this purpose, a single micro-crystal of [160GdPc2]− was mounted on the µSQUID-EPR chip and measured at a base temperature of 30 mK. The crystal was aligned so that the easy axis and preferably also the hard axis (see below) of the MM are parallel to the applied fields, that is, µ0H|| along the easy axis and µ0Htr along the hard axis or at least within the hard–medium plane (Fig. 1a). This alignment is often challenging and typically requires multiple placement attempts, as the external vector field is strictly confined to the µSQUID plane (see “Methods”). Figure 1c shows the µSQUID-EPR “frequency map” (see “Methods”) and the corresponding Zeeman diagram for [160GdPc2]− with µ0H|| precisely along the easy axis and µ0Htr = 0. As a pre-requisite for this study, the µSQUID-EPR allows precise determination of spin Hamiltonian parameters40.
Fig. 1: Frequency-dependent µSQUID-EPR Investigation.
a Schematic figure of µSQUID-EPR chip comprising the µSQUID loops and the coplanar waveguide for microwave irradiation. b Crystallographic orientation of the [160GdPc2]− complex with the easy axis (blue arrow), medium (green), and hard (red) axes, and crystal packing of the molecules, showing the orientation of the molecule in the unit cell along the c-crystallographic axis (right panel). Colour code: Gd, purple; N, cyan; C, grey. Hydrogens and the Et4N+ counter cation are omitted for clarity. c Zeeman diagram for the [160GdPc2]− complex as determined via µSQUID-EPR studies and parameters as described in the text. The three marked tunnel splitting can be identified as Δ1: (ms =) −5/2 → + 1/2, Δ2: −5/2 → + 3/2 and Δ3: −7/2 → + 1/2. The numbers from 0 to 7 correspond to the energy levels from the lowest to the highest in energy. d Frequency map (ΔM(H||,ν)) with the easy axis applied along the easy axes of the crystal. The frequencies shown in this panel correspond to the energy positions for the investigation of spin interference effects. e Zoomed region for each explored tunnel gap.
The field-dependent frequency map (ΔM(H||,ν)) exhibits several resonant absorption peaks, all of which can be fitted employing by the following axial parameters: \({B}_{2}^{0}\) = −680.3(5) MHz, \({B}_{4}^{0}\) = −1.45(1) MHz and g = 2.0 (see Figs. 1c, d and S3), for a Hamiltonian of the form (1):
$${{{H}}}_{{{\rm{Gd}}}}=g{\mu }_{B}{\mu }_{0}{{\bf{H}}}{{\boldsymbol{.}}}\hat{{{\bf{S}}}}+{\sum }_{n=1}^{3}{B}_{2n}^{0}{O}_{2n}^{0}+\left({B}_{2}^{2}{O}_{2}^{2}+{B}_{4}^{4}{O}_{4}^{4}\right).$$
(1)
here, the first term represents the electronic Zeeman interaction, while \({O}_{k}^{q}\) and \({B}_{k}^{q}\) are the Extended Stevens operators42,43. The second term in the equation comprises the axial ligand field parameters, while the third term consists of the two transverse ligand field parameters. The term \({B}_{6}^{0}\) was excluded from the fitting as it tends to zero. Notably, transitions originating from excited states (e.g., (2→3) in Fig. 1d) are observed, consistent with partial thermal population of these levels. Although the bath temperature is 30 mK, the relative populations of these states indicate an effective spin temperature exceeding 500 mK under microwave irradiation, with the precise value depending on the applied microwave power, pulse width, and delay (see SI Section 2).
The µSQUID-EPR setup also allows the application of the field along any direction (θ, with respect to the easy axis of [160GdPc2]−) within the µSQUID plane (Fig. 1b). The presence of transverse fields results in the observation of EPR forbidden transitions with (Δms ≠ 1), which carry detailed information regarding the spin Hamiltonian of the system. The fitting of the µSQUID-EPR “angular map” (ΔM(H,θ), see “Methods”), likewise, provides access to the transverse ligand field parameters of [160GdPc2]−, i.e., \({B}_{2}^{2}\) = -273(3) and \({B}_{4}^{4}\) = 3.0(3) MHz (See Fig. S4). These values are consistent with previously determined parameters40. The presence of \({B}_{2}^{2}\) is a consequence of reduced symmetry from the ideal D4d. A similar effect has been also observed in the archetypal [Mn12] complex44. Note that although the fourth‑order parameters appear small compared to the second‑order terms, a meaningful comparison should consider the scaled quantities such as \({B}_{2}^{0}{S}^{2}\) and \({B}_{4}^{0}{S}^{2}\).
The frequency map (ΔM(H||,ν) in Fig. 1a, and insets obtained with higher frequency resolution permit the resolution of the avoided crossings, i.e., hybridisation that quantifies the resonant QTM (or tunnel splitting, Δi) between different |ms〉 states. The three tunnel splittings (six, including both polarities of the longitudinal field) with the highest visibilities are indicated by the circles in Fig. 1c, d. These three gaps correspond to Δ1 = |−5/2〉 →|+1/2〉, Δ2 = |−5/2〉 →|+3/2〉 and Δ3 = |−7/2〉 →|+1/2〉. We denote these as: −ms → + ms–n, hence, for Δ1: −ms = −5/2 → ms–n = ½ (i.e., n = 2 (even parity)), and similarly for Δ2: n = 1 (odd Parity) and Δ3: n = 3 (odd Parity), as also described in earlier works14,39. Here, the ms are the diabatic state labels, the true eigenstates are mixed and do not correspond uniquely to these ms values. Denoting the change of spin at each transition as Δmi, we have Δm1 = 3 at Δ1 and Δm2,3 = 4 at Δ2,3.
A major advance with respect to our previous work40, however, is the direct experimental observation of how the entire spin manifold (Zeeman diagram) of the MM, including Δi, reacts to transverse fields. Frequency maps at different Htr were collected with Htr aligned nearly along the hard axis of the MM. This is presented as an animation (SI. V1) comprising ~ 72 h of continuous data collected with sweep rate for H|| < 20 mT/s (adiabatic sweep), 0.1 GHz steps in frequency, and 4 mT steps in Htr. Some of the frames, captured in Fig. 2 at different constant Htr, indicate the non-trivial oscillations in Δi(Htr). At certain Htr values, the gap Δ1 (even n) closes, while Δ2 (odd n) tends to approach its maxima, evidencing the direct signature of the parity effect (see below). The animation and Fig. 2 also reveal that the entire spin manifold’s reaction to small Htr is mostly contributed by these oscillations of various tunnel gaps.
Fig. 2: Transverse field frequency-map variation.
Frequency map (ΔM(H||,ν)) variation upon transverse field (Htr) application, highlighting the oscillating behaviour of Δ1–3 at Htr a −22 mT, b −10 mT, c +2 mT, d +14 mT, e +26 mT, f +38 mT, g +50 mT, and h +62 mT. The frequency maps were collected with the field (H||) applied along the easy axes of the crystal.
Once the oscillating gaps are detected (confirming the direction of Htr), further investigation can be achieved by fixing the microwave radiation frequency (i.e., ν = 6.05 GHz, 8.90 GHz in Fig. 1d) and varying Htr, i.e., the absorption peak intensities (ΔM) plotted with H|| and Htr. Figure 3 shows the ΔM(H||, Htr) maps with H|| varied (at <20 mT/s) along the easy axis, in the presence of different constant Htr values. The separations between two neighbouring absorption peaks (δH||) associated with a Δi are approximately monotonic functions of the corresponding Δi (see Fig. 1 inset and “Methods”); hence, a study of δH||(Htr) gives access to Δi(Htr). Notably, the δH||(Htr) clearly shows oscillatory features with the indication that the minima in Δ1 coincide with the maxima in Δ2,3 (see the enlarged sections in Fig. 3a, b).
Fig. 3: Transverse field study at fixed frequency.
Transverse field dependent absorption maps (ΔM(H||, Htr)) between Htr = ±120 mT at a fixed frequency of a 6.05 GHz and b 8.90 GHz. Corresponding exact numerical simulations employing Eq. (1) are shown in (c, d), respectively. The thick coloured lines and the zoomed regions highlight the oscillating behaviour of the tunnel splitting Δi upon Htr variation.
In addition to the oscillations of δH||(Htr) corresponding to Δ1–3, near H|| ~ 0, the maps exhibit a pronounced tiling pattern, i.e., periodic features along the Y axis, due to the topological effects or QPI on the whole manifold. The oscillations of zero-longitudinal-field QTM gaps (−ms → + ms) plausibly cause this periodic feature, as these gaps have a significant role in offsetting different states in the Zeeman diagram. These features can be quantitatively addressed by exact numerical diagonalisation of (1) under the influence of transverse fields (Fig. 3c, d). The simulations capture the periodic (tiling) features, while their enlarged sections show the oscillating gaps from the pair transitions43,45 (see Section 6 in SI and animated Zeeman diagrams: SI V2).
A closer inspection of Fig. 3a, b reveals that each oscillating pair of lines is accompanied by a second pair whose separation increases steadily with transverse field, indicating that the prominent features consist of one pair with an oscillating spacing δH||(Htr) and another with a monotonically increasing spacing. The same effect causes each of the resonant absorption lines in Fig. 2 to develop into two lines at higher transverse fields (see Fig. 2g), where one pair exhibits an oscillating tunnel gap and the other a monotonic increase with Htr. This behaviour is due to two molecular orientations of the unit cell, leading to a hard–medium plane alignment, with parallel easy axis arrangement (Fig. S1). At 5% dilution, it is expected that the applied Htr aligns with the hard axis and medium axes of the statistically distributed molecules at both sites within the unit cell (see Fig. S1). However, their independent responses remain distinguishable and are advantageous in practice, as they offer additional constraints to uniquely determine the direction of Htr. For simplicity, however, we simulate and discuss one molecular orientation at a time.
Topological quenching and parity effect
Corroboration of the topological quenching of the tunnelling gaps can be gained by rotating the crystal to align Htr at different angles in the hard-medium plane of the MMs while maintaining the easy axis aligned with H|| (Fig. 4). When Htr is aligned nearly along the hard axis (φ = 100° where φ denotes the angle between the medium-axis and Htr), the most prominent oscillations are observed (solid points in Fig. 4a) in Δ1–3(Htr), as extracted (see “Methods”) from the corresponding δH||(Htr) in Fig. 3a, b. On the other hand, when Htr is not along the hard axis (φ = 10°), the oscillations diminish (solid points in Fig. 4b showing Δ1(Htr) for different crystal orientations and in Fig. S6 for Δ2,3(Htr)).
Fig. 4: Spin parity effects and influence of transverse terms on DPs.
a Measured tunnel splittings Δ1,2,3 (solid points) as a function of transverse fields, compared to the numerical simulations (lines) for φ = 100° (angle between Htr and medium axis) and \({B}_{4}^{4}\) = 3.9 MHz; b Measured tunnel splittings Δ1-vs-transverse fields (solid points) for different crystal orientations, compared to the numerical simulations (lines) for \({B}_{4}^{4}\) = 3.9 MHz and φ = 100°, 55° and 10°, respectively; c Simulated tunnel splittings Δ1,2 as a function of transverse fields for different \({B}_{4}^{4}\) with a fixed φ = 100° and d different φ with a fixed \({B}_{4}^{4}\); e, f Anti-crossing energy levels simulated as a function of Htr in the Hard-medium plane and corresponding differences (Δi) as contour maps. The white, red arrows in (f) represent the position of the DPs with and without \({B}_{4}^{4}\) contribution. Panels a and b also show their associated error bars.
Likewise, it can be noted that the 2nd minima of Δ1(Htr) are nearly aligned with the maxima in Δ2,3(Htr), a sign of topological quenching of tunnel gaps and the parity effect. A period of ~ 70 mT is evident in Δ1(Htr). Only three minima are observed in Δ1, as only ~3 diabolic points are expected in Δ1 (see later). In contrast to Mn12 and Fe8 (integer S), the Δ1(Htr) with even-n (in −ms → + ms–n) shows a minimum at zero transverse field, while Δ2,3(Htr) for odd-n exhibit maxima. This behaviour aligns with Kramer’s degeneracy for a half-integer spin system (S = 7/2), where even-n tunnel gaps must vanish at zero field.
Another difference to [Mn12]39 and [Fe8]14 is that the expected maxima in Δ2,3(Htr) at zero Htr (Fig. 4a) seem to be diminished, with the minima moving towards smaller Htr. This observation is particularly compelling as topologically quenched tunnelling emerges at relatively low transverse fields, making their presence experimentally undeniable and strongly motivating a deeper investigation into the role of the 4th-order transverse ligand field parameter in spin systems.
The minima in the observed oscillations in Fig. 4 are the DPs. Encircling such a point in the transverse-field space of either of the two intersecting levels accumulates a phase of \(\pi\), whereas paths that do not enclose the DP yield an accumulated phase of 0. This quantized 0/\(\pi\) behavior is the Longuet–Higgins46 subcase of the Berry phase and serves as the topological invariant. In spin systems, several DPs (for example, 84 DPs exist for a S = 7/2 system26) for each pair of levels arise due to QPI or, in other words, topological quenching, driven by transverse magnetic fields, and transverse ligand field parameters (mainly E or \({B}_{2}^{2}\)) that create anisotropy in the lateral plane and two dominant tunnelling paths. These paths accumulate different Berry phases, leading to interference patterns modulated by transverse fields. DPs can be modelled through methods such as the Feynman path integral (Instanton)15,47,48,49 and the discrete WKB50,51,52. The former offers some intuition on the system’s topology as it requires a classical analogue energy surface. We, therefore, employ the classical analogous energy diagram for [160GdPc2]−, with different signs and hypothetical magnitudes of 4th order anisotropic parameter (\({B}_{4}^{4}\)), by allowing the continuous orientation of the spin in Eq. (1). Figure S7 indicates that \({B}_{4}^{4}\) leaves its imprints on the instantons (classical least action paths15,47,48,49), especially its large values can lead to four instantons instead of two, as predicted for fourfold symmetry. Classical energy diagrams, however, fail to capture details of the QPI or topological quenching.
Accounting for the quantum (or semi-classical) calculations, the tunnel splittings are very much more sensitive to small \({B}_{4}^{4}\) values. The exact analytical calculations for such spin-Hamiltonians can be convoluted, as exemplified in simpler systems52. The predicted period of oscillations if only one axial and one transverse parameter in the Hamiltonian is considered is:
$$\Delta {H}_{1}=4.8\times 10^{-5}({2k}_{B}/g{{{\upmu }}}_{B})\sqrt{2{B}_{2}^{2}({B}_{2}^{2}+{3B}_{2}^{0})},$$
(2)
where \({k}_{B}\) is the Boltzmann constant, \(g\) ≈ 2.00; as derived using Eq. 10 in ref. 15, or Eq. 2 in ref. 14, replacing anisotropy terms with \({B}_{2}^{0}\), \({B}_{2}^{2}\) parameters. Note that this is an incomplete theoretical description of the experimental system as the observations have indicated significant effect of 4th order parameters. For [160GdPc2]− this expression yields a period ≈ 80 mT for \({B}_{2}^{0}\,=\,-680.3\,{{\rm{MHz}}},\,{B}_{2}^{2}\,=\,-273\,{{\rm{MHz}}}\), i.e., comparable to the experimentally observed periods ~ 70 mT. A \({B}_{4}^{4}\) dependent shift (toward the medium axis) of certain DPs at large transverse fields was predicted for non-zero \({B}_{2}^{0}\), \({B}_{2}^{2}\), \({B}_{4}^{4}\) (see ref. 53). In contrast, our observations suggest an opposite trend, i.e., the shift of certain DPs at small transverse fields. This highlights a sign-dependent role of \({B}_{4}^{4}\) not fully explored (see section 9 in SI for more details). Although several of the studies have focused14,39 on integer‑spin systems, the theoretical framework that describes DPs in anisotropic spin Hamiltonians is formulated for arbitrary spin values and is not intrinsically limited to integer spins. Consequently, the qualitative behaviour of DPs extends seamlessly to half‑integer systems. Except for the distinction associated with the Kramers-degeneracy—the tunnel-gaps that vanish at Htr = 0 have even n for half-integer spin systems and odd n for the integer spin systems—the evolution of DPs follows the same symmetry principles in both classes, depending on the transverse and axial anisotropy terms. Aside from this difference, the evolution and symmetry of DPs in both classes of systems are governed by the same underlying axial and transverse anisotropy terms. Thus, our half‑integer S = 7/2 findings—especially the \({B}_{4}^{4}\)-dependent shifts—may be compared qualitatively with the earlier works14,39, even though exact numerical agreement is not anticipated.
The interpretation of the experimental results of [160GdPc₂]−, hence, must consider a Hamiltonian with non-zero \({B}_{2}^{0}\), \({B}_{2}^{2}\), \({B}_{4}^{0}\), \({B}_{4}^{4}\), H|| taking into account the sign dependence of \({B}_{4}^{4}\) (relative to \({B}_{2}^{2}\)). Here, we examine the DP shifts induced by the fourth‑order transverse anisotropy terms (SI, Section 9), whereas a detailed perturbative framework that includes spin‑parity and path‑dependent corrections will be provided in a subsequent theoretical analysis. In addition, to allow the inclusion of all parameters and modelling of DPs under 2D Htr, here we use the exact numerical diagonalisation. Tunnel gaps for the chosen DPs were extracted from the simulated Zeeman diagrams at different Htr (see Section 6 in SI and SI.V2).
It is worth noting that the DPs discussed here are distinct from clock transitions54,55, even though both involve field-dependent extrema in energy levels and relate to decoherence resilience in different ways. Clock transitions occur when the first-order field derivative of a transition energy vanishes over a broad field-range, thereby reducing sensitivity to magnetic field-noise. By contrast, DPs are true degeneracies created by destructive quantum interference and are central to designing geometric quantum gates in molecular spin systems. Whether the clock-transition gap lies exactly between two DPs (i.e., the constructive interference point), or elsewhere, depends on its transverse-field (and transverse-parameter) dependent landscape, an aspect that merits further investigation.
Numerical simulation of fourth-order transverse parameter effects & Berry phase
For the system studied here, numerical analysis shows that axial (\({B}_{2}^{0}\)) and transverse (\({B}_{2}^{2}\)) ligand field parameters alone cannot fully account for the observed Δ(Htr). To reach the optimum fitting, we simulate Δ1,2 for different hypothetical \({B}_{4}^{4}\) values (including zero) with other parameters fixed from angular and frequency maps (Fig. 4c). The transverse field angle φ = 100°, i.e., close to the hard axis, was chosen to match experimental data and reveals parity-dependent oscillations in Δ1,2. At Htr = 0, the gap Δ1 (Δmᵢ = 4) is quenched while Δ2,3 (Δmᵢ = 3) has maxima modulated by \({B}_{4}^{4}\). Transitions with Δmᵢ = 4 are more sensitive to \({B}_{4}^{4}\) than that with Δmᵢ = 3, since the 4th order terms remain in the matrix (transition) elements involving Δmi = 4. Increasing \({B}_{4}^{4}\) from negative to positive values suppresses the Δ₂ maximum, indicating that a large \({B}_{4}^{4}\) with the opposite sign to \({B}_{2}^{2}\) disrupts constructive interference (at Htr = 0) between the QTM pathways for odd n.
Figure 4d shows simulated Δ1(Htr) and Δ2(Htr), for various angles (φ) of Htr in the hard-medium plane, using a \({B}_{4}^{4}\) value consistent with experimental data. Note that φ is defined relative to the medium axis; hence, φ = 90° shows the QPI (oscillations) most prominently, while φ = 0° shows more classic behaviour of Δ1,2(Htr). Eventually, the unique set of parameters used to fit the experimental findings in Fig. 4a and b are: g = 2.00 and \({B}_{2}^{0}\) = −680.3, \({B}_{4}^{0}\) = −1.45, \({B}_{2}^{2}\) = −273, \({B}_{4}^{4}\) = 3.9 MHz. Some of the pre-obtained values, acquired from angular/frequency map fitting, carried significant uncertainty. In contrast, the fittings of Δ1,2,3(Htr) offer excellent precision in transverse parameters, especially \({B}_{4}^{4}\), surpassing conventional orientation mapping techniques. To investigate the role of \({B}_{4}^{4}\) on topological quenching or DPs, and particularly the reason behind the contrast with a theoretical work53, we simulated the relevant DPs in the 2D space Htr (hard-medium plane) and the extracted gaps between intersecting states, as shown in Fig. 4e, f. As we vary \({B}_{4}^{4}\) up to and beyond the estimated values in [160GdPc2]−, see Fig. S8, we find that the DPs for Δ1 (Δmi = 3) do not move with changing \({B}_{4}^{4}\), while the DPs at Δ2 (Δmi = 4) at small Htr, first move inward towards each other and then away towards the medium axis (Fig. S8m–o), in contrast to the predictions53. We notice that the opposite sign of \({B}_{4}^{4}\) (compared to that estimated in [160GdPc2]−), indeed matches the prediction in ref. 53 (Fig. S8p–r), i.e., the DPs at large Htr move towards the medium axis. To clarify this sign-dependent aspect further, we simulated a tunnel gap at zero H|| for the extreme hypothetical scenarios (\({B}_{4}^{4}\) > 0 and <0 with \({B}_{2}^{2}=0\)) and discussed how, in the presence of non-zero \({B}_{2}^{2}\), the shift of DPs must depend on the sign of \({B}_{4}^{4}\) (see section 9 in SI). Hence, our observations point to a sign-dependent role of \({B}_{4}^{4}\), complementing the predictions in ref. 53, and opening the possibility of shifting DPs toward lower transverse fields Htr, thereby facilitating experimental accessibility.
Finally, to numerically confirm the topological nature of the gap minima, we evaluated the Berry phase acquired by the two energy levels near a degeneracy under a closed circuit in the transverse magnetic field. The accumulated phase was computed by numerically diagonalizing the Hamiltonian and determining the phase associated with the selected eigenstate along a circular trajectory in parameter space (see https://doi.org/10.5281/zenodo.19451754). When the path encircles a gap minimum, the system acquires a Berry phase of π, a result observed for several level pairs. This quantized response confirms that the gap minima correspond to true DPs displaying the Longuet–Higgins form of the Berry phase46.
In conclusion, we have demonstrated a magneto-spectroscopic method for the direct determination of QPI and parity effects in 4f-MMs. The technique, exploiting the resonant absorption of the M(H) loops and the ability to apply a transverse field along any direction in the x-y plane, allows not just the precise determination of the spin Hamiltonian parameters but also the controlled quenching of QTM at precisely determined transverse fields. Furthermore, the non-trivial evolution of the DPs and their dependence upon \({B}_{4}^{4}\) term can, in principle, be exploited for the design of robust MMs towards QTM, making these systems more attractive towards quantum technologies. Although the experimental techniques are insufficient to distinguish between Abelian and non-Abelian holonomies, the presence of genuine DPs suggests that more advanced parameter-space control could, in principle, enable access to non-Abelian holonomies for HQC—an aspirational yet not unreachable prospect. Hence, identifying such degeneracies in MMs is central both to simple geometric quantum gate designs and to the development of more advanced holonomic architectures in the future. Moreover, the observed motion of DPs toward experimentally accessible transverse fields suggests new synthetic strategies for chemists: by tailoring spin centre arrangements and local field environments, it may be possible to harness topologically quenched tunnelling rates. Together, these results not only deepen our understanding of spin dynamics but also chart a path forward for the rational design of next-generation quantum materials.