Introduction

Electronic spins in solids are natural candidates to encode qubits, the basic units of future quantum computers1. Quantized spin projections give rise to a discrete set of states that can be coherently manipulated by external magnetic field pulses. In addition, they are insensitive to electric noise and can therefore show longer coherence times than qubits based on electric (e.g. superconducting) circuits. Outstanding examples are magnetic defects in semiconductor hosts, such as NV centres in diamond2 or P+ donors in silicon3. These systems are however not easy to tune, and wiring them up into a scalable architecture remains very challenging. A promising alternative to address these questions relies on the use of electron spins hosted by artificial molecules, i.e., organic radicals or transition-metal complexes4,5,6. This chemistry-based approach affords the synthesis of macroscopic numbers of identical qubits. The design of the metal coordination provides some control over the spin levels, which determine the qubit frequencies7,8. Furthermore, the spin coherence can be enhanced by either reducing the number of neighbouring nuclear spins or by engineering qubit states that are insensitive to magnetic field fluctuations. As a result, molecular spin qubit coherence times TM have shown very significant improvements over the past years8,9,10,11,12,13.

But arguably the most appealing trait of this approach is that it offers the possibility of scaling up computational resources within each molecular unit. The most straightforward strategy is to incorporate several (two or more) magnetic centres, each of them realizing a qubit. Several spin carriers can be assembled in a rational way, such that their nature, disposition and interaction result in a set of addressable transitions that allow coherently manipulating all qubit states6. At the synthetic level, this implies a necessary dissymmetry between the different magnetic centres. Early successful examples in either homometallic14,15 or heterometallic16,17 lanthanide dimers have shown that conditions to realize both CNOT and √SWAP two-qubit gates can be fulfilled18, and allowed measuring the quantum coherence of a CNOT gate with TM of about 400 ns19,20,21. Also, modular supramolecular approaches able to incorporate multiple weakly coupled Cr7Ni qubits have been developed22 and used to implement C-PHASE two-qubit gates23. A second strategy is to profit from the internal spin levels of each magnetic centre to create d-dimensional (d > 2) quantum systems, that is, qudits. Examples include the coherent manipulation of the multiple nuclear spin states of 159TbPc2 (I = 3/2, d = 4)24,25, Et4N[163DyPc2] (I = 5/2, d = 6)26 and 173Yb(trensal) (I = 5/2, d = 6)27 mononuclear complexes. The former of these enabled the first realization of a three-level Grover quantum search algorithm in a single molecule, using tools of molecular electronics to read-out the qudit state28. Electronic spin qudits can also be realized in some cases, for instance by using the S = 7/2 spin manifold of a Gd(III) ion provided that it has a sufficiently low magnetic anisotropy to make all levels experimentally accessible29.

The integration of multiple levels in well-defined complexes increases the density of quantum information handled by these systems and could allow embedding specific functionalities, say quantum error correction or some simple quantum algorithms and simulations27,30,31,32, at the molecular scale. These applications demand increasing the dimension of the computational space beyond three qubits while retaining the ability to perform universal quantum operations30,33. Meeting both conditions in a molecule still represents a daunting challenge.

Here we show that it is possible to design and synthesize such building blocks via the combination of the two strategies outlined above. The rational design of a coordination complex with multiple and magnetically inequivalent Gd(III) ions would allow to scale up to six or more addressable qubits within one single molecule, with the potential of realizing different quantum gate operations. We focus here on the realization of this scaling-up in a dissymmetric [Gd2] molecule in which the two Gd(III) ions have a slightly different coordination14,15. The isolated properties of each Gd(III) d = 8 qudit are first determined through the isolation and study of the [LaGd] and [GdLu] analogue molecules, in which the Gd(III) ions occupy either of the two coordination sites present in [Gd2]. The results show that the [Gd2] molecule holds the promise to act as a six-qubit quantum processor or as a d = 64 electronic spin qudit.

Results

Synthesis and structures

We previously showed that the ligand H3L forms a series of dissymmetric homometallic dilanthanide complexes of formula (Hpy)[Ln2(HL)3(NO3)(py)(H2O)], where H3L = 6-(3-oxo-3-(2-hydroxyphenyl)propionyl)-pyridine-2-carboxylic acid and py = pyridine, among which the [Gd2] compound studied here14,15. The two different Ln environments resulting from the different coordination pockets created by the three HL2− ligands (HL2− = double deprotonated H3L) and their disposition give rise to two metal sites with markedly different Ln–O and Ln–N bond lengths, one position being systematically larger than the other. This allows forming in a controlled manner heterometallic dilanthanide molecules, provided their ionic radii are sufficiently different16,17,34. This unique synthetic strategy was first used to isolate the [CeEr], [LaEr] and [CeY] analogues19. Here, in order to study each Gd(III) ion present in [Gd2], albeit isolated from the other, the [LaGd] and [GdLu] analogues are prepared in a similar manner. Note that the differences in ionic radii Δri, respectively of 9.4 and 7.7 pm34, are sufficient to expect the formation of homogeneous heterometallic single phases in both cases16,17. Thus, equimolar amounts of Gd(NO3)3 and either La(NO3)3 or Lu(NO3)3 were made to react with H3L in pyridine. The layering of the resulting clear yellow solutions with Et2O resulted in homogeneous phases of yellow crystals that were found by single-crystal X-ray diffraction to consist of (Hpy)[LaGd(HL)3(NO3)(py)(H2O)]·5py ([LaGd]) and (Hpy)[GdLu(HL)3(NO3)(py)(H2O)]·5py ([GdLu]).

Compounds [LaGd] and [GdLu] crystallize in the monoclinic space group P21/n, with the asymmetric unit coinciding with their formula and Z = 4 (see Supplementary Table 1 and the CIF files available as Supplementary Data 1, for [LaGd], and Supplementary Data 2, for [GdLu]). The pyridinium cation forms a hydrogen bond with one of the carboxylic oxygen atoms from one of the HL2− ligands. Two of the five lattice pyridine molecules are hydrogen bonded to the coordinated water molecule (Supplementary Fig. 1, Supplementary Table 2). The metal complexes (Fig. 1) show the structure observed consistently within this series of compounds, with two lanthanide ions bridged and chelated by three HL2− ligands in two opposite orientations with respect to the La···Gd or Gd···Lu vector. Thus two different coordination environments are present, with two/one (O,N,O) pockets and one/two (O,O) chelates respectively for site 1/2. The coordination sphere is completed by a nitrate ion on site 1 (La in [LaGd] and Gd in [GdLu]), and by one pyridine molecule and a water molecule on site 2 (Gd in [LaGd] and Lu in [GdLu]). In [LaGd], the coordinated nitrate is bidentate, resulting in a coordination number CN of 10 for the largest La(III) ion, while the rest have CN of 9. The intermetallic separations are respectively 3.8941(3) and 3.761(1) Å for La···Gd and Gd···Lu (Supplementary Table 3), to be compared with the intermediate 3.804(1) Å Gd···Gd separation in [Gd2]. The average Ln–O bonds to HL2− donors are longer at site 1, 2.544 Å (La in [LaGd]) and 2.432 Å (Gd in [GdLu]), than at site 2, 2.405 Å (Gd in [LaGd]) and 2.342 Å (Lu in [GdLu]). The differences in average bond distances ΔO between both sites are thus respectively 0.139 and 0.090 Å, larger than those observed in the homometallic analogues [La2] (ΔO = 0.045 Å) and [Gd2] (ΔO = 0.026 Å)14, an indication that the lanthanide ions have taken selectively the optimal relative positions following their respective ionic radii16,17.

Fig. 1: Molecular structures of [LaGd] and [GdLu] spin qudits.
figure 1

For clarity, hydrogen atoms not involved in H-bonding and lattice pyridine molecules are omitted. Atomic sites involved in the coordination are labelled. Colour code: O, red; N, light blue; C, grey.

The coordination sphere at both Gd sites also differs in shape, as indicated by continuous shape measures (Supplementary Fig. 2)35,36. Thus, the Gd2 site in [LaGd] is closest to an ideal capped antiprism (C4v symmetry) and to an ideal tricapped trigonal prism (D3h symmetry). The Gd1 site in [GdLu] is more irregular as indicated by distances to any ideal polyhedron 3–4 times larger. These differences are identical to those found between the two Gd sites in [Gd2]. The metal site composition is confirmed for both compounds by the unreasonable relative displacement parameters and worse final refinement obtained for any other distribution of metals within the molecular structure.

The molecular metal composition observed in the solid-state is confirmed by Electrospray Ionization Mass Spectrometry (ESI-MS), which also allows ensuring the integrity of the molecular structure in solution. Thus, the spectra obtained from solutions of [LaGd], [Gd2] and [GdLu] each exhibit a prominent peak with respectively m/z of 1146.98, 1166.01 and 1183.02 and isotopic distribution that agree with the corresponding [LnLn’(H2L)(HL)2]+ fragment. These fragments correspond to the molecular moiety in the absence of the terminal ligands, which are the only species detected by the technique of the dissociation equilibrium. In the case of [GdLu], a small impurity of [Gd2] is observed indicating a process of partial scrambling. Indeed, the dissociation of the terminal ligands may conduce to a relaxation of the structure and therefore to a decrease of the selectivity as it was shown by previous ESI-MS and density functional theory studies with compounds of the same family16,17. Nevertheless, the selectivity remains important and follows the relative Δri, which for [LaGd] and [GdLu], is relatively high. While the ESI-MS technique does not allow a quantitative analysis, since the [Gd2] impurity arises from the dissociated fraction (expected to be marginal) and the selectivity of the full molecule is that of the solid state, we are confident the amount of impurity in solution remains negligible. Therefore, the heterometallic molecular composition seen in the solid-state structure is maintained in solution, which is relevant for the studies of coherent spin dynamics done on diluted solutions. These studies confirm that any fraction of [Gd2] in the solution of [GdLu] is not detectable since this would have an effect on the quantum coherence in the opposite direction as that observed (see below).

Isolated qudits in [LaGd] and [GdLu]: magnetic dissymmetry

A crucial requirement to properly define a qudit, or N qubits, using a system with multiple energy levels d = 2 N is that the energy spectrum has some nonlinearity, i.e., that the levels are not simply equidistant as those of a harmonic oscillator30. In the case of Gd(III), with its ground state S = 7/2, this condition relies on the existence of a finite magnetic anisotropy7,29. An advantage of this ion is that, because of its L = 0 configuration, the intrinsic anisotropy of the free ion is negligible. Any zero-field splitting of the d = 8 spin levels necessarily arises from small distortions that the coordination environment induces on the close to spherical 4f electronic shell. This property makes Gd(III) a kind of model crystal field probe and allows modifying the magnetic anisotropy via changes in the local coordination. Often, this anisotropy is also quite small, much smaller than those typically found for other lanthanides with L ≠ 04,5,37,38,39, thereby making these levels accessible via conventional magnetic spectroscopy techniques.

The continuous-wave electron paramagnetic resonance (cw-EPR) spectra of [LaGd] and [GdLu] are shown in Fig. 2. The experiments have been performed on powdered samples using both X-band (frequency ω/2π =  9.886 GHz) and Q-band (ω/2π = 33.33 GHz) spectrometers. Spectra measured on each sample at different frequencies do not simply scale with H/ω, where H is the external magnetic field. This shows already that Gd(III) acquires a net magnetic anisotropy in the two possible coordination sites 1 and 2 (Fig. 1). Besides, the spectra of the two samples are also different, thus showing that the magnetic anisotropy depends on the coordination environment (Fig. 2e).

Fig. 2: Magnetic spectroscopy of [LaGd] and [GdLu] spin qudits.
figure 2

X-band and Q-band cw-EPR spectra measured on powdered samples of [LaGd] (a, b) and of [GdLu] (c, d) at T = 6 K. The red solid lines are experimental data, whereas the blue solid lines are simulations obtained with the spin Hamiltonian (1) and the parameters given in the text. The latter have been shifted down for clarity. Panel e compares the splitting of the Gd(III) spin levels at zero magnetic fields in both molecules, which correspond to the coordination sites 1 and 2 in the structures of Fig. 1. The wave functions of all level doublets are also indicated in this panel.

In order to render these arguments quantitative, we have performed fits of the experimental spectra using the EPR simulation package EasySpin40. Taking into account the low symmetry of the Gd coordination sites in these molecules, we have considered the simplest lowest-order spin Hamiltonian

$${\cal{H}} = - DS_z^2 + E\left( {S_x^2 - S_y^2} \right) - g\mu _B{\mathbf{HS}},$$
(1)

where D and E are zero-field splitting parameters and g is the electron spin g-factor. The results of these fits are compared to the experimental results in Fig. 2. A reasonably good agreement was obtained by setting D = 2 GHz, ǀEǀ = 0.67 GHz and g = 1.99 for [LaGd] and D = 2.4 GHz, ǀEǀ = 0.8 GHz and g = 1.99 for [GdLu]. The fact that the orthorombicity ǀEǀ/D is close to the maximum value (1/3) probably arises from the low symmetry of both coordination sites. The line broadenings suggest that there are sizeable distributions in D and E, of about 60% for both compounds.

Another experimental technique that provides information on the structure of electronic energy levels is heat capacity. Results obtained for the two molecular “monomers” are shown in Fig. 3. Above 10 K, the specific heat cP/R is dominated by excitations of vibrational modes. This lattice contribution had been measured directly on the diamagnetic derivative [YLa]19, and agrees very well with the high-temperature behaviour. The additional anomaly observed, at H = 0, for both [LaGd] and [GdLu] below 2 K must then be associated with the finite splitting of the Gd(III) electronic spin levels. A nice feature of these results is that the position of the specific heat maximum provides a direct measure of the overall zero-field splitting41.

Fig. 3: Field-dependent specific heat of [GdLu] and [LaGd] spin qudits.
figure 3

Curves were measured at zero and different applied magnetic fields on powdered samples of [GdLu] (a) and [LaGd] (b). Data measured for the diamagnetic complex [YLa]19, which arise from vibrational excitations only, are also shown in both panels. The solid lines are numerical calculations of the magnetic heat capacity, derived from Eq. (1) with the same parameters, and their distributions, that account for the cw-EPR data of Fig. 2 to which the nonmagnetic contribution has been added.

The comparison between results measured for [LaGd] and [GdLu] confirms that the magnetic anisotropy is slightly stronger for the latter. The Schottky-like broad maximum is indeed well accounted for by numerical calculations based on Eq. (1) using the same D and E parameters derived from EPR experiments (and with the same D- and E- strains). A good agreement is also found for data measured under non-zero magnetic fields. Results for the two compounds become progressively closer to each other as H increases, due to the relatively weak magnetic anisotropy of Gd(III) and to the predominance of the Zeeman term in Eq. (1) for sufficiently strong H.

In conclusion, the results shown in this section demonstrate that Gd(III) ions coordinated in the molecular structures of [LaGd] and [GdLu] have low-lying electronic energy level schemes that provide a basis for two different spin qudits.

Exchange coupling in [Gd2]

In this section, we turn our attention to the molecular dimer [Gd2]. This complex hosts the two magnetic ions, in different coordination sites, whose properties in isolation have been discussed in the previous section. This discussion suggests that this molecule, with (2S + 1) × (2S+1) = 64 = 26 unequally spaced levels, provides a proper implementation of a d = 64 qudit or of six qubits. However, an additional necessary ingredient, related to the condition of universality that we discuss below, is the existence of a net coupling between the two spins.

In order to get information of the spin-spin interaction within the molecule, we have compared the magnetic, spectroscopic and thermal properties of the [Gd2] complex with those measured on [LaGd] and [GdLu]. Figure 4 shows the magnetic response of the three samples. While the χT plots of [LaGd] and [GdLu] agree with the predictions for isolated Gd(III) ions (an almost temperature-independent value, in agreement with Curie’s law), the data for [Gd2] show a decrease below ~4 K. This behaviour is compatible with the existence of a weak isotropic coupling described by the following spin Hamiltonian

$${\cal{H}} = {\cal{H}}_1 + {\cal{H}}_2 - J{\mathbf{S}}_1{\mathbf{S}}_2$$
(2)

where H1 and H2 are the spin Hamiltonians of each isolated Gd(III) ion, given by Eq. (1) with the appropriate parameters, and J is the spin-spin coupling constant. As a further simplification, the anisotropy axes at sites 1 and 2 have been taken as collinear. A reasonably good agreement with the experimental results is found for an antiferromagnetic interaction with J =  −0.42 GHz (see Fig. 4). This interaction is likely dominated by intramolecular dipole-dipole interactions, which would give rise to −0.92 GHz ≤ J ≤ +2.73 GHz, depending on the orientation of the main anisotropy axis z. Therefore, Eq. (2) must be regarded as a simplified description, with an effective isotropic J, of the energy level scheme and the overall magnetic moment of [Gd2].

Fig. 4: Magnetic coupling in the [Gd2] dimer.
figure 4

Product of the molar susceptibility χ times temperature measured on powder samples of [LaGd], [GdLu] and [Gd2]. The lines are numerical calculations. In the case of the molecules hosting only one magnetic ion, they are derived from the spin Hamiltonian (1) with the parameters determined independently from EPR experiments (Fig. 2). In the case of [Gd2], the calculations include also a spin-spin coupling, as in Eq. (2), with two values for the interaction constant J.

The results of EPR measurements performed on [Gd2] are also compatible with the existence of magnetic interaction between the two Gd ions. Figure 5 shows the results measured at X-band and Q-band. These spectra are not simple superpositions of those measured on the [LaGd] and [GdLu] monomers (see Figure S3). The strength of J is, however, quite small as compared with that of the single-ion anisotropies. Whereas the latter give rise to zero-field splittings of order 20–30 GHz (Fig. 2) the energy scale of the JS1S2 term is about 4 GHz. For this reason, it is not possible to accurately determine J solely from EPR data. Yet, as Fig. 5 shows, the results are compatible with calculations performed with the same J inferred from magnetic measurements, albeit with an additional broadening that might point to the existence of some ‘J-strain’ or, as it might be expected from the discussion above, some anisotropy in the coupling between the two Gd spins.

Fig. 5: Magnetic spectroscopy of the [Gd2] dimer.
figure 5

X-band (a) and Q-band (b) cw-EPR spectra measured on a powder sample of [Gd2] at T = 6 K. The curves at the top are experimental results, whereas those at the bottom are results of numerical calculations derived from Eq. (2) with J = −0.42 GHz.

Similar considerations apply to the results of heat capacity experiments, which are shown in Fig. 6. The Schottky anomaly associated with the magnetic anisotropy of both ions dominate data measured above 0.35 K. Still, these data are compatible, within the uncertainties of the experiment and the underlying model, with the predictions derived from Eq. (2) for J =  −0.42 GHz. For sufficiently strong magnetic fields (see Fig. 6b), the differences between the specific heat data of isolated and coupled spins (and even between those of [LaGd] and [GdLu]) tend to vanish, as expected.

Fig. 6: Field-dependent specific heat of the [Gd2] dimer.
figure 6

Specific heat of [Gd2], compared to those of monomers [LaGd] and [GdLu] and of the diamagnetic [YLa]19, measured at zero magnetic fields (a) and under a strong 3 T magnetic field (b). The dashed line represents the case in which the two spins in [Gd2] would not interact. It was obtained by adding the results measured on [LaGd] and [GdLu]. The solid line shows the results of a numerical calculation based on Eq. (2) with J = −0.42 GHz.

Magnetic energy level scheme: six-qubit encoding and universal operation

The results described in previous sections show that Eq. (2) provides a reasonably good account of the low-lying magnetic energy levels of [Gd2] and of all experimental quantities that derive from it. We next discuss, on the basis of this description, the potential of this molecular dimer to implement multiple qubits. The energy scheme of this molecule in a magnetic field, shown in Fig. 7a, consists of a set of 64 unequally spaced levels, with level separations between adjacent levels of a few GHz. This scheme obviously admits a labelling of the levels in terms of the states of a qudit (from |n = 1〉 for the ground level to |n = 64〉 for the highest excited one) or in terms of the states of six qubits (say, from |000000〉 to |111111〉). However, this condition, i.e., that the space dimension is large enough, is not sufficient to ensure that such a small “processor” could perform universal quantum operations. The condition of universality implies that any quantum superposition of the basis states (the 64 mentioned above) can be generated starting from any initial state, e.g., starting from the system initialized in its ground state. It implies that there exist non-forbidden transitions between different levels, which can be univocally addressed by setting the frequency of a microwave pulse or the external magnetic field, and which form a complete set, in the sense, defined above, of allowing “visiting” all possible spin states.

Fig. 7: Six-qubit encoding and universality in the [Gd2] dimer.
figure 7

a Energy level scheme of [Gd2] calculated with Eq. (2) and the parameters given in the text. The magnetic field is oriented along the diagonal between the magnetic axes x, y and z of the two Gd(III) ions, which are taken as collinear for simplicity. Possible labelling of the levels in terms of the basis states of six qubits is shown, in which adjacent levels differ only in the state of one qubit. b Colour map of the Rabi frequencies, calculated at 0.5 T, linking different levels, numbered from 1, for the ground level, up to 64, for the highest energy one. c Map of transitions accessible by concatenating different allowed resonant transitions at the same field, showing that the [Gd2] system allows universal operations.

Figure 7b shows the Rabi frequencies for photon-induced resonant transitions between any two levels. This plot shows a dense map of allowed transitions, with rates exceeding 10–15 MHz mT−1, which link the ground state with any other basis state. In order to perform a more rigorous demonstration, we have adapted the mathematical proof for universality42, based on the formalism of Lie algebras, to this particular situation. Figure 7c shows that any two states of the basis can be connected by concatenation of resonant transitions. This question has received relatively little attention in connection with molecular qubits, but it is not as simple to achieve as it might appear at first sight. Its relationship with the formal structure of the underlying spin Hamiltonian can be better understood by comparing the situation of the two constituent Gd ions, realized independently in [LaGd] and [GdLu] complexes, with that of [Gd2]. Applying the same mathematical method, we find that the Hamiltonian given by Eq. (1) affords a complete set of operations for each qudit (see Supplementary Fig. 4). However, Eq. (2), which describes the two qudits in [Gd2], is universal only when J ≠ 0 (compare Fig. 7 with Supplementary Fig. 5, which illustrates the situation for J = 0). The reason for this condition, which we anticipated above, is that conditional operations between states of the two Gd(III) ions can be implemented with simple resonant pulses if and only if each of them has some anisotropy and the two are magnetically coupled.

Quantum spin coherence and relaxation

The previous section shows that [Gd2] provides a suitable platform for a six-qubit quantum processor, in the sense of having the right Hilbert space and a sufficient number of allowed transitions. However, it also sets a final stringent condition, namely that transitions between these levels can be implemented coherently. The latter depends on the spin-coherence and spin-relaxation times, which we discuss in this section. The spin dynamics has been experimentally studied, at 6 K, on diluted solutions of [LaGd], [GdLu] and [Gd2] in MeOH-d4:EtOH-d6, with concentrations in the range 0.38–0.61 mmol L−1. As is well known9, the use of deuterated solvents reduces the hyperfine couplings with the solvent nuclear spins and enhances the electronic spin coherence. The comparison with results measured for a sample of [GdLu] diluted in a non-deuterated solvent mixture, shown in Supplementary Fig. 6a, confirms this enhancement. For all three molecular systems, Electron Spin-Echo (ESE) signals are observed over the entire 10–810 mT range of magnetic fields used, either using 2- or 3-pulse sequences (Supplementary Figs. 69).

Therefore, both the isolated qudits in [LaGd] and [GdLu] and the exchange-coupled pair in [Gd2] present measurable quantum coherence over the full magnetic field range, thus supporting the picture discussed in the previous section. The echo intensity varies with H, which allows obtaining echo-induced EPR spectra. These are in good agreement with the cw spectra, although the presence of a strong modulation of the echo decay (see below) results in variations with the time interval between pulses (Supplementary Figs. 7 and 8).

Phase memory times TM have been obtained by measuring the spin-echo decay following a Hahn sequence of π/2 and π pulses (typically 16 ns and 32 ns long, respectively) separated by a varying interval τ. For all three compounds, the ESE intensity decreases in a similar manner with increasing τ, the decay being slightly more rapid in the case of [Gd2]. The ESE decay exhibits a strong modulation with a frequency that increases with H, independently of the compound. To quantify these effects, the ESE decays were modelled (Supplementary Fig. 10) through the equation

$$y\left( \tau \right) = y_0 + A_{2p}e^{ - 2\tau /T_{\mathrm{M}}}\left\{ {1 + ke^{ - \lambda \tau }\cos \left( {2\pi \nu \tau + \phi } \right)} \right\},$$
(3)

in which y0 is a constant background, accounting for the floating level of the electronic output signal, A2p is the initial amplitude, k the relative amplitude of the modulated signal, ν the average frequency of the oscillating component and λ its decay rate, and ϕ the non-zero phase due to the detector dead-time. The magnetic field dependence of A2p exhibits a maximum in the range 300–450 mT, slightly broader for [GdLu] than [GdLa] (Supplementary Fig. 12). The [Gd2] data are further broadened and cannot be obtained as a combination of those of [LaGd] and [GdLu], in agreement with the cw-EPR results (see Fig. 5 and Supplementary Figs. 3 and 9). The phase coherence times of the isolated qubits in [LaGd] and [GdLu] are similar, although slightly larger for the later (Fig. 8). In both cases, TM decreases with H, respectively from 1.05/1.31 μs at 300 mT to 0.73/0.85 μs at 750 mT. On the contrary, TM of [Gd2] remains virtually constant at ca. 0.73 μs for all applied fields. Remarkably, the phase memory times derived for the isolated qubits in [LaGd] and [GdLu] are about twice longer than those found in a 1% magnetically diluted crystal of the polyoxometallate [GdW30]29. The significantly shorter coherent times found for [Gd2] are in agreement with it being more sensitive to decoherence. This effect can be tentatively associated with the coupling between the two spins, which results in the presence of more possible excitations contributing to decoherence. Within this picture, states of the two coupled Gd(III) ions feel decoherence sources affecting either site 1 or site 2. As it happens with other properties (see e.g. Fig. 6), the differences between TM of the monomers and of [Gd2] are also progressively reduced as the magnetic field increases, likely because the Zeeman interaction begins then to dominate over the spin–spin coupling.

Fig. 8: Spin coherence of the molecular qudits.
figure 8

Comparison of the ESE-detected EPR spectra (a) and the field dependence of the phase memory times TM (b), measured on diluted MeOH-d4:EtOH-d6 solutions of [LaGd], [GdLu] and [Gd2] at 6 K.

Spin-lattice relaxation times T1 have been determined through spin-recovery measurements after a π – T – π/2 – τ – π pulse sequence, with a varying length of the first interval T and a fixed τ. In all three compounds, the intensity increases with T. These curves were modelled (Supplementary Fig. 11) with a single exponential function to determine T1. The field-dependences of the inversion recovery amplitude AIR (Supplementary Fig. 13) resemble those of A2p, again with broader maxima for [Gd2], that cannot be obtained as a combination those of [LaGd] and [GdLu]. Values of T1 for [LaGd] and [GdLu] lie in the range 2.5–3.5 μs (Supplementary Fig. 13) and depend only weakly on H. As for TM, values of T1 for the exchange-coupled [Gd2] are significantly lower, ca. 1.5–2 μs.

These experiments give also a hint on possible sources of decoherence. The frequency of the ESE modulation in the 2-pulse experiments is very close to the Larmor frequency of the 2H nucleus over the whole range of magnetic fields and for the three compounds (Supplementary Fig. 12). The modulation can therefore be ascribed to the interaction with the “remote” deuterons of the solvents. The coupling to each solvent nucleus is very weak, but the sum of many contributions gives rise to the large modulation observed. The origin of the modulation was confirmed by measuring a solution of [GdLu] in non-deuterated solvents (Supplementary Fig. 6a) for which the frequency modulation corresponds to the larger Larmor frequency of 1H. This suggests that in these experiments hyperfine interactions with the solvent are still very relevant and, therefore, that the maximally attainable TM for isolated molecules can be even higher than the values shown in Fig. 8.

Coherent control: Rabi oscillations

Spin nutation experiments provide information on the Rabi frequencies, to be compared with TM, and serve to illustrate the difficulties that experiments on frozen solutions face when trying to implement quantum gates in spin qudits. Experiments were performed at 6 K on the same diluted MeOH-d4:EtOH-d6 solutions of [LaGd], [GdLu] and [Gd2] described above. They involve the measurement of the ESE generated by a variable duration pulse (0.1 ≤ tp ≤ 1.6 μs) refocused by a π pulse, with τ intervals (100 ≤ τ ≤ 200 ns) between the two pulses and between the refocusing pulse and ESE detection. Representative results are shown in Fig. 9, for [Gd2], and in the supplementary material (Supplementary Figs. 1416), for [LaGd] and [GdLu]. In all cases, a strongly damped coherent oscillation is observed. The decay is much faster than what one would expect for a single coherent transition between states with decoherence times TM ≈ 0.8–1.2 μs, even if one considers possible inhomogeneities in the amplitude B1 of the micro-wave field pulse43,44. The fast decay arises instead from the destructive interference of Rabi oscillations between different states, with also different Rabi frequencies, that the initial microwave pulse is able to excite in a sample of randomly oriented molecules. A Fourier transformation of the time-dependent signals shows indeed sizeable distributions of Rabi frequencies in the three samples. At 410 mT and an attenuation of 10 dB, which corresponds to B1 ≤ 0.275 mT29, we find average ΩR ≈ 12 MHz for [LaGd], ΩR ≈ 17 MHz for [GdLu] and ΩR ≈ 20 MHz for [Gd2], which lie within the range of values expected for these systems (see Fig. 7 and Supplementary Fig. 4).

Fig. 9: Spin nutation for the [Gd2] dimer.
figure 9

Experiments were performed at T = 6 K and for a 410 mT magnetic field, that is, close to the maximum observed in the echo-induced spectrum, and for three different microwave pulse powers (attenuations of 10, 14 and 20 dB corresponding to B1 ≈ 0.275 mT, 0.16 mT and 0.092 mT, respectively), on a diluted MeOH-d4:EtOH-d6 solution of [Gd2]. For each set of data, the Fourier transform is shown as inset, evidencing oscillations with characteristic frequencies that correspond to the Larmor frequencies of 2H and 1H, as indicated, in addition to the main Rabi frequency.

In addition, these Fourier plots reveal also two additional oscillations. Their frequencies are independent of B1, unlike ΩR, but proportional to the external magnetic field, and agree well with the Larmor frequencies of the 2H and 1H nuclei (Fig. 9 and Supplementary Figs. 1416). They likely result from a cross-polarization of these nuclei by the Gd spin, known as Hartmann-Hahn effect45, and confirm that the Gd spin is coupled to the surrounding deuterium and proton nuclear spins, contributing to the quantum spin decoherence.

Discussion

In this work, we have designed a molecular structure that can host two Gd(III) ions in distinct coordination environments. Combined with the intrinsic properties of this ion, in particular, its S = 7/2 and zero orbital moment, this design enables realizing two different spin qudits, each with d  = 8 or, equivalently, three qubits. These qudits can be studied in isolation, in either [LaGd] and [GdLu] molecules in which they are accompanied by a diamagnetic counterpart or can be both integrated in the same molecular unit. In the latter case, we have shown that the molecular dimer [Gd2] meets all conditions needed to realize six qubits (or a d = 64 qudit).

The integration of a high number of addressable quantum states in well-defined molecular units represents one of the distinctive traits of the chemical approach to quantum technologies and can provide a number of competitive advantages in this field. Such molecules can locally act as tiny quantum processors themselves, the equivalent of NISQS (Noisy Intermediate Size Quantum Systems) in other schemes46, thus implementing simple algorithms. Of special relevance would be the correction of errors, e.g. phase errors that are likely to be most disturbing in the case of spin qubits, at the molecular scale27,30. For this, it is not even necessary that the dimension of the Hilbert space defined by the available spin states be a power of two. These molecular NISQS can also serve as a suitable basis where to map quantum simulations of simple systems, such as small molecules33,47. An advantage of using a single physical unit to perform these simple computations is that it reduces the number of non-local operations, which require switching on and off interactions between different parts of a quantum circuit, and that are often more sensitive to decoherence and, therefore, error31,32.

The fulfilment of this potential defines also a series of important challenges. A prominent one is how to individually address each of the relevant transitions. This requires working on well-oriented molecules, to avoid the difficulties discussed above in connection with nutation experiments and decreasing significantly the linewidth of each resonance. Even though coherence times are reasonably long in these systems, the resonance lines show a relatively large inhomogeneous broadening, which mainly arises from the random orientation of molecules in frozen solutions and from the existence of some distributions in the parameters that define the spin Hamiltonian (exchange and anisotropy constants). These difficulties can possibly be overcome by synthesizing well-organized molecular frameworks, in which all molecules need to be oriented in the same manner and, at the same time, magnetically diluted or sufficiently far away from each other48,49,50. In the end, the ideal situation would be to explore the response of individual molecules, either via the application of single-molecule electronics24,25,28 or by enhancing the current sensitivity of magnetic spectroscopic techniques51,52,53,54.

A second challenge is how to move forward beyond the level of a few qubits. As it becomes clear by inspecting the energy level scheme of a d = 64 qudit, the resonant transitions that provide the basic set of operations suffer already from a “frequency crowding”, which would seriously hinder addressing individual resonances for qudits of increasing dimension. Therefore, scaling beyond this point must proceed by linking different units via coherent mediators. Implementing such switchable couplings within molecules (e.g. via chemical linkers sensitive to some external stimuli, such as light) remains a challenging goal22,55,56. The alternative is to combine a chemical, bottom up, approach with the integration of functional molecules into solid-state circuits. A recent proposal for a scalable architecture considers the possibility of applying superconducting circuits, namely on-chip resonators, to control, read-out and communicate magnetic molecules57,58. A further advantage of these devices, and of dealing with electronic spins having sizeable energy splittings, is that they are compatible with working at very low temperatures, which are necessary to initialize the qudit to its ground state (its population becomes larger than 99.99% for T < 0.1 K and a magnetic field of 0.5 T). Within this scheme, medium size molecular qudits, able to embed some critical functionalities such as phase error correction, would provide very attractive building blocks to reach computational performances difficult to match by other platforms.

Methods

Synthesis

6-(3-oxo-3-(2-hydroxyphenyl)propionyl)pyridine-2-carboxylic acid (H3L). The ligand H3L was synthesized as described previously14.

(Hpy)[Gd2(HL)3(NO3)(py)(H2O)]·5py ([Gd2]). Compound [Gd2] was synthesized as described previously15. Purity was checked by elemental analysis, mass spectrometry and multiple cell determinations on single crystals, including full structure determination.

(Hpy)[LaGd(HL)3(NO3)(py)(H2O)]·5py ([LaGd]). A 10 mL pyridine solution of ligand H3L (30.0 mg, 0.105 mmol) is added slowly to a 10 mL pyridine solution of La(NO3)3·6H2O (15.2 mg, 0.035 mmol) and Gd(NO3)3·6H2O (15.8 mg, 0.035 mmol). The mixture is stirred at room temperature for one hour, yielding a clear shiny yellow solution that is layered with diethylether. Crystals of [LaGd] form slowly and are recovered after one weak in a 63 % yield. ESI-MS: [LaGd(HL)2(H2L)]+ m/z = 1146.98. Elemental analysis found (calc.) for [LaGd(H2L)(HL)2(NO3)(py)(H2O)]·3H2O: C 44.11 (44.16); H 2.77 (3.04); N 4.87 (5.15). IR (KBr, ʋ/cm−1): 3410 mb, 1618s, 1583 s, 1560 m, 1528 s, 1463 m, 1400 s, 1384 s, 1298 m, 1239w, 1204w, 1148w, 1120w, 1059 m, 949 m, 891w, 760 m, 706 m, 663 m, 635w, 568w.

(Hpy)[GdLu(HL)3(NO3)(py)(H2O)]·5py ([GdLu]). A 10 mL pyridine solution of ligand H3L (30.0 mg, 0.105 mmol) is added slowly to a 10 mL pyridine solution of Gd(NO3)3·6H2O (15.8 mg, 0.035 mmol) and Lu(NO3)3·6H2O (16.4 mg, 0.035 mmol). The mixture is stirred at room temperature for one hour, yielding a clear yellow solution to which 200 μL conc. HNO3 is added. After 10 more min stirring, the solution is layered with diethylether. Crystals of [GdLu] are recovered after 2 weeks in a 53% yield. ESI-MS: [GdLu(HL)2(H2L)]+ m/z = 1183.02. Elemental analysis found (calc.) for [GdLu(H2L)(HL)2(NO3)(py)(H2O)]·6.7H2O: C 40.33 (41.06); H 2.56 (3.34); N 5.50 (4.79). IR (KBr, ʋ/cm−1): 3422 mb, 1619s, 1585 s, 1560 m, 1528 s, 1465 m, 1401 s, 1384 s, 1325 m, 1299 m, 1239w, 1208w, 1147w, 1121w, 1058 m, 950 m, 892w, 756w, 708 m, 666w, 637w, 569w.

Magnetic measurements

Magnetic measurements were performed using a Quantum Design SQUID MPMS-XL magnetometer through the Physical Measurements unit of the Servicio de Apoyo a la Investigación-SAI, Universidad de Zaragoza. All data were corrected for the sample holders and grease contributions, determined empirically as well as for the intrinsic diamagnetism of the sample, estimated using Pascal constants.

Heat capacity experiments

Heat capacity data were measured, down to T = 0.35 K, with a commercial physical property measurement system (PPMS, Physical Measurements unit of the Servicio de Apoyo a la Investigación-SAI, Universidad de Zaragoza) that makes use of the relaxation method. The samples, in powder form, were pressed into pellets and placed onto the calorimeter on top of a thin layer of Apiezon N grease that fixes the sample and improves the thermal contact. The raw data were corrected from the known contributions arising from the empty calorimeter and the grease.

Electron paramagnetic resonance experiments

Continuous-wave (cw) EPR measurements were performed with a Bruker Biospin ELEXSYS E-580 spectrometer operating in the X-band and Q-band. Solid-state cw-EPR measurements were performed at RT on polycrystalline samples placed in quartz tubes.

In addition, pulsed time domain (TD) measurements were performed at X-band frequencies. In these experiments, the typical widths of the π/2 and π pulses were 16 and 32 ns, respectively. In order to avoid unwanted echoes, a 2-step (4-step) phase cycle was used in the 2-pulse (inversion recovery and three-pulse) experiments. The high power microwave excitation was obtained by using a TWT amplifier. A dielectric low Q cavity from Bruker was used as resonator. TD-EPR measurements were done on frozen solutions, for which gas-flow Helium cryostats were used. The polycrystalline solids were dissolved in a 1:1 mixture of fully deuterated methanol and ethanol. The use of deuterated solvents is intended to limit decoherence due to protons. The TD-EPR experiments were performed at 6 K on 0.61 ([LaGd]), 0.38 ([GdLu]) and 0.44 ([Gd2]) mmol L−1 solutions. Measurements on 0.13 and 0.08 mmol L−1 solutions of [GdLu] were also performed to discard any variation of T1 and TM in the range of concentrations used. Additional TD-EPR measurements were performed on a 0.38 mmol/L solution of [GdLu] in a 1:1 mixture of non-deuterated methanol and ethanol, in order to show the effect of the solvent nuclear spins on the spin coherence times and on the modulation of the ESE signals measured after 2-pulse and 3-pulse sequences. Finally, cw-EPR measurements were also done on the same frozen solutions, albeit in the range 20–80 K due to signs of saturation at 6 K.

Single-crystal X-ray diffraction

Data for compounds [LaGd] and [GdLu] were collected at 100 K on a Bruker APEX II QUAZAR diffractometer equipped with a microfocus multilayer monochromator with MoKα radiation (λ = 0.71073 Å), respectively on a yellow lath and a yellow needle of dimensions 0.86 × 0.18 × 0.12 and 0.40 × 0.03 × 0.03 mm3. Data reduction and absorption corrections were performed with SAINT and SADABS59, respectively. Both structures were solved by intrinsic phasing with SHELXT60 and refined by full-matrix least-squares on F2 with SHELXL-201461. In the case of [GdLu], remaining voids in the structure with no significant electron density peaks were analysed and taken into account with PLATON SQUEEZE62 that recovered a total of 224 electrons per cell, over four equivalent voids of 193 Å3. These figures being reasonable for one diffuse pyridine solvent molecule per void, i.e., one per formula unit, this has been reflected in the reported formula. For both structures, the metal site composition was confirmed by refining the structure with the two homometallic situations as well as with the Ln sites inverted. These resulted in relatively poorer agreement factors and most importantly in unrealistic combinations of Ueq values at the metal sites.

Theory: universality test for molecular spin qudits

A universal quantum processor must be able to implement any unitary operation within the computational states embedded in a Hilbert space. In the particular case discussed here, i.e., qudits simulating qubits: any transition between two molecular states must be driven by applying an external field. To show if a particular molecular qudit is universal (in this sense) or not, we write the spin Hamiltonian as,

$${\cal{H}} = {\cal{H}}_s - g\mu _B{\mathbf{B}}_1({\boldsymbol{t}}){\mathbf{S}}.$$
(4)

Here \(\vec S\) is the total spin operator and B1(t) is a time-dependent external magnetic field that induces resonant transitions between different eigenstates of the spin Hamiltonian \({\cal{H}}_s\). These eigenstates form the computational basis. Since the dimension d of the Hilbert spaces of all spin qudits considered in this work is d = 2N, with N = 3 for [LaGd] and [GdLu], and N = 6 for [Gd2], these eigenstates may be denoted as either |n〉 (with n  = 1, …, d) or as |00…0〉 to |11…1〉. If all of these states are accessible from any other one, the qudit realizes a N-qubit processor.

One way to check universality is as follows. We use the fact that the energy spectrum has some nonlinearity, i.e., that the levels are not simply equidistant as those of a harmonic oscillator. This arises from the combination of single-ion anisotropy, the dissymmetry between the two Gd(III) coordination sites and the mutual interaction between the two Gd(III) spins. This means that we can address a transition between any two states, say |n〉 and |m〉, by making the frequency ω of the driving magnetic field B1(t) resonant with the transition frequency \(\omega _{nm} = (E_n - E_m)/\hbar\), provided that the matrix element \(\left\langle n \right|{\mathbf{B}}_1\left( {\boldsymbol{t}} \right){\mathbf{S}}\left| m \right\rangle \ne 0\). If this happens, it is fair to say that we can implement any unitary operation of the form \(e^{iH_{nm}t}\) with \(H_{nm} = \left| n \right\rangle \left\langle m \right| + h.c.\) The allowed transitions define a set \(L = \{ H_{nm},H_{n\prime m\prime }, \ldots \}\).

It is however expected that L will not cover all the necessary transitions. The natural question is then which extra transitions can be implemented by concatenating the different elements of L. The formal answer42 to this question is that the accessible transitions are those belonging to the Lie algebra \({\cal{L}}\) generated by L. This is a natural consequence of the Lie-Trotter formulas, i.e.

$$e^{i\left[ {H_{nm},H_{n\prime m\prime }} \right]t} = \mathop {{\lim }}\limits_{C \to \infty } \left( {e^{\frac{{H_{nm}t}}{{\sqrt C }}}e^{\frac{{H_{n\prime m\prime }t}}{{\sqrt C }}}e^{ - \frac{{H_{nm}t}}{{\sqrt C }}}e^{ - \frac{{H_{n\prime m\prime }t}}{{\sqrt C }}}} \right)^C.$$
(5)

The different commutators of the elements of L are also elements of the corresponding Lie group, \({\cal{L}}\). In practice, to find the allowed transitions, we may compute the Hnm by expressing BB1(t)S in the basis of eigenstates of \({\cal{H}}_s\) and compute all possible commutators. By doing so, we can check if \({\cal{L}}\) covers the full Hilbert space. In a finite dimensional basis, this is easily computed since the non-zero commutators are of the form:

$$i\left[ {H_{n\,m},H_{n\,m\prime }} \right] = i\left| m \right\rangle \left\langle {m\prime } \right| + h.c.$$
(6)

We finally notice that the proof says nothing about how to perform operations. It only checks if they are possible.

Universality of [GdLu] and [LaGd] three-qubit systems (or qu8its)

These molecules have a spin Hamiltonian \({\cal{H}}_s\) given by Eq. (1). It is rather simple to prove that this Hamiltonian leads to universal operations. It turns out that for intermediate magnetic fields, say \(\mu _0H_z \sim 0.5 - 1\,{\mathrm{T}}\), the level spectrum consists of non-equidistant levels and that \(\left\langle n \right|S_x\left| {n + 1} \right\rangle \ne 0\) for all eigenstates |n〉 (see Supplementary Fig. 4). From the point of view of the universal operation, this is rather favourable since \(i\left[ {H_{n\,n + 1},H_{n + 1\,n + 2}} \right] = i\left| n \right\rangle \left\langle {n + 2} \right| + h.c.\), and so on. Therefore, every transition can be performed by concatenating commutators.

Universality of the [Gd2] six-qubit system (or qu64it)

Putting together several qudits, they must interact as a prerequisite to be universal. The [Gd2] dimer is described by Eq. (2), which includes a weak, antiferromagnetic interaction, thus it fulfils this condition. Besides, from the previous discussion it follows that if some symmetry is shared by \({\cal{H}}_s\) and B1(t)S, the system is presumably not universal. For example, a parity selection rule can impose that only transitions\(\left\langle n \right|{\mathbf{B}}_1\left( {\boldsymbol{t}} \right){\mathbf{S}}\left| {n + 2} \right\rangle \ne 0\). Then, this system becomes not universal: starting from an even labelled state |n〉, the odd states cannot be accessed, and viceversa. This condition might introduce difficulties in situations that are closer to that described by the spin Hamiltonian of [Gd2]. The isotropic Heisenberg model, \({\cal{H}} = - J{\mathbf{S}}_1{\mathbf{S}}_2 - g\mu _{\mathrm{B}}{\mathbf{H}}({\mathbf{S}}_1 + {\mathbf{S}}_2)\) conserves the total spin S. Therefore, a spin dimer described by this model is definitely not universal for an external driving BB1(t)S, with S = S1 + S2. Transitions between different total spin states are forbidden. Fortunately, in [Gd2] the different anisotropy terms acting on each spin break this symmetry, besides providing a non-equidistant level spectrum. However, contrary to what happens in the single qudit case, the dc magnetic field H cannot be much higher that J, D and E. The reason is twofold. If H is strong enough, the two Gd(III) ions are effectively decoupled. At the same time, the anisotropy terms become less important and the system is closer to an isotropic Heisenberg system. Our numerical calculations, whose results are shown in Fig. 7 and in Supplementary Fig. 5, confirm this. Figure 7b and the top panels of Supplementary Fig. 5 show maps of the Rabi frequencies for resonant transitions between any pair of states. Using transitions with Rabi frequencies larger than 0.2 MHz mT−1, we construct contour plots (Fig. 7c and bottom panels of Supplementary Fig. 5) of those transitions that are feasible by concatenation of all possible commutation operators. The system is universal if the matrix built in this way spans all points (except the diagonal). Figure 7 shows the universal character of the dimer of two coupled Gd(III) ions in a not too strong magnetic dc-field. Then, Supplementary Fig. 5 shows how this property is broken down by either setting J = 0 or increasing the dc-external field. The conclusion is that [Gd2] can, under the appropriate conditions, simulate either a six-qubit processor or two decoupled three-qubit systems.