A dissymmetric [Gd2] coordination molecular dimer hosting six addressable spin qubits

Artificial magnetic molecules can host several spin qubits, which could then implement small-scale algorithms. In order to become of practical use, such molecular spin processors need to increase the available computational space and warrant universal operations. Here, we design, synthesize and fully characterize dissymetric molecular dimers hosting either one or two Gadolinium(III) ions. The strong sensitivity of Gadolinium magnetic anisotropy to its local coordination gives rise to different zero-field splittings at each metal site. As a result, the [LaGd] and [GdLu] complexes provide realizations of distinct spin qudits with eight unequally spaced levels. In the [Gd2] dimer, these properties are combined with a Gd-Gd magnetic interaction, sufficiently strong to lift all level degeneracies, yet sufficiently weak to keep all levels within an experimentally accessible energy window. The spin Hamiltonian of this dimer allows a complete set of operations to act as a 64-dimensional all-electron spin qudit, or, equivalently, as six addressable qubits. Electron paramagnetic resonance experiments show that resonant transitions between different spin states can be coherently controlled, with coherence times TM of the order of 1 µs limited by hyperfine interactions. Coordination complexes with embedded quantum functionalities are promising building blocks for quantum computation and simulation hybrid platforms.

E lectronic spins in solids are natural candidates to encode qubits, the basic units of future quantum computers 1 . 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 diamond 2 or P + donors in silicon 3 . 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 complexes [4][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 frequencies 7,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 T M have shown very significant improvements over the past years [8][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 states 6 . At the synthetic level, this implies a necessary dissymmetry between the different magnetic centres.
Early successful examples in either homometallic 14,15 or heterometallic 16,17 lanthanide dimers have shown that conditions to realize both CNOT and √SWAP two-qubit gates can be fulfilled 18 , and allowed measuring the quantum coherence of a CNOT gate with T M of about 400 ns [19][20][21] . Also, modular supramolecular approaches able to incorporate multiple weakly coupled Cr 7 Ni qubits have been developed 22 and used to implement C-PHASE two-qubit gates 23 . A second strategy is to profit from the internal spin levels of each magnetic centre to create ddimensional (d > 2) quantum systems, that is, qudits. Examples include the coherent manipulation of the multiple nuclear spin states of 159 TbPc 2 (I = 3/2, d = 4) 24,25 , Et 4 N[ 163 DyPc 2 ] (I = 5/2, d = 6) 26 and 173 Yb(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 state 28 . 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 accessible 29 .
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 simulations 27,[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 operations 30,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 [Gd 2 ] molecule in which the two Gd(III) ions have a slightly different coordination 14,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 [Gd 2 ]. The results show that the [Gd 2 ] 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 H 3 L forms a series of dissymmetric homometallic dilanthanide complexes of formula (Hpy)[Ln 2 (HL) 3 (NO 3 )(py)(H 2 O)], where H 3 L = 6-(3-oxo-3-(2-hydroxyphenyl)propionyl)-pyridine-2-carboxylic acid and py = pyridine, among which the [Gd 2 ] compound studied here 14,15 . The two different Ln environments resulting from the different coordination pockets created by the three HL 2− ligands (HL 2− = double deprotonated H 3 L) 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 different 16,17,34 14 , an indication that the lanthanide ions have taken selectively the optimal relative positions following their respective ionic radii 16,17 .
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 (C 4v symmetry) and to an ideal tricapped trigonal prism (D 3h 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 [Gd 2 ]. 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], [Gd 2 ] 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'(H 2 L) (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 [Gd 2 ] 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 family 16,17 . Nevertheless, the selectivity remains important and follows the relative Δr i , which for [LaGd] and [GdLu], is relatively high. While the ESI-MS technique does not allow a quantitative analysis, since the [Gd 2 ] 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 [Gd 2 ] 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 oscillator 30 . In the case of Gd(III), with its ground state S = 7/2, this condition relies on the existence of a finite magnetic anisotropy 7,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 ≠ 0 4,5,37-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).
In order to render these arguments quantitative, we have performed fits of the experimental spectra using the EPR simulation package EasySpin 40 . Taking into account the low symmetry of the Gd coordination sites in these molecules, we have considered the simplest lowest-order spin Hamiltonian 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 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 c P /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 splitting 41 .
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 [Gd 2 ]. In this section, we turn our attention to the molecular dimer [Gd 2 ]. 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 = 2 6 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 [Gd 2 ] 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 [Gd 2 ] show a decrease below~4 K. This behaviour is compatible with the existence of a weak isotropic coupling described by the following spin Hamiltonian where H 1 and H 2 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  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.
[GdLu] e) Energy from ground level (GHz) 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 [Gd 2 ].
The results of EPR measurements performed on [Gd 2 ] 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 zerofield splittings of order 20-30 GHz (Fig. 2) the energy scale of the JS 1 S 2 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.
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.
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 [Gd 2 ] 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. 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 universality 42 , 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 [Gd 2 ]. 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 [Gd 2 ], 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 [Gd 2 ] 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 [Gd 2 ] in MeOH-d 4 :EtOH-d 6 , with concentrations in the range 0.38-0.61 mmol L −1 . As is well known 9 , 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 3pulse sequences (Supplementary Figs. 6-9).
Therefore, both the isolated qudits in [LaGd] and [GdLu] and the exchange-coupled pair in [Gd 2 ] 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) Phase memory times T M 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 [Gd 2 ]. 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 in which y 0 is a constant background, accounting for the floating level of the electronic output signal, A 2p 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. As it happens with other properties (see e.g. Fig. 6), the differences between T M of the monomers and of [Gd 2 ] are also progressively reduced as the magnetic field increases, likely because the Zeeman interaction begins then to dominate over the spin-spin coupling. Spin-lattice relaxation times T 1 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 T 1 . The field-dependences of the inversion recovery amplitude A IR (Supplementary Fig. 13 Fig. 13) and depend only weakly on H. 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 2 H 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 1 H. This suggests that in these experiments hyperfine interactions with the solvent are still very relevant and, therefore, that the maximally attainable T M for isolated molecules can be even higher than the values shown in Fig. 8. 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 T M ≈ 0.8-1.2 μs, even if one considers possible inhomogeneities in the amplitude B 1 of the micro-wave field pulse 43,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 B 1 ≤ 0.275 mT 29 , we find average Ω R ≈ 12 MHz for [LaGd], Ω R ≈ 17 MHz for [GdLu] and Ω R ≈ 20 MHz for [Gd 2 ], which lie within the range of values expected for these systems (see Fig. 7 and Supplementary  Fig. 4).
In addition, these Fourier plots reveal also two additional oscillations. Their frequencies are independent of B 1 , unlike Ω R , but proportional to the external magnetic field, and agree well with the Larmor frequencies of the 2 H and 1 H nuclei ( Fig. 9 and Supplementary Figs. 14-16). They likely result from a crosspolarization of these nuclei by the Gd spin, known as Hartmann-Hahn effect 45 , 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 [Gd 2 ] 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 schemes 46 , 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 scale 27,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 molecules 33,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, error 31,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 welloriented 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 other [48][49][50] . In the end, the ideal situation would be to explore the response of individual molecules, either via the application of single-molecule electronics 24,25,28 or by enhancing the current sensitivity of magnetic spectroscopic techniques [51][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 goal 22,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 molecules 57,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.
(Hpy)[Gd 2 (HL) 3 (NO 3 )(py)(H 2 O)]·5py ([Gd 2 ]). Compound [Gd 2 ] was synthesized as described previously 15 . Purity was checked by elemental analysis, mass spectrometry and multiple cell determinations on single crystals, including full structure determination.   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 Xband 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 (4step) 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 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 mm 3 . Data reduction and absorption corrections were performed with SAINT and SADABS 59 , respectively. Both structures were solved by intrinsic phasing with SHELXT 60 and refined by full-matrix least-squares on F 2 with SHELXL-2014 61 . In the case of [GdLu], remaining voids in the structure with no significant electron density peaks were analysed and taken into account with PLATON SQUEEZE 62 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 U eq 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, HereS is the total spin operator and B 1 (t) is a time-dependent external magnetic field that induces resonant transitions between different eigenstates of the spin Hamiltonian 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 = 2 N , with N = 3 for [LaGd] and [GdLu], and N = 6 for [Gd 2 ], 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 B 1 (t) resonant with the transition frequency ω nm ¼ ðE n À E m Þ= h, provided that the matrix element n h jB 1 t ð ÞS m j i≠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 ¼ n j i m h j þ h:c: The allowed transitions define a set L ¼ fH nm ; H n0m0 ; g. 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 answer 42 to this question is that the accessible transitions are those belonging to the Lie algebra L generated by L. This is a natural consequence of the Lie-Trotter formulas, i.e.
The different commutators of the elements of L are also elements of the corresponding Lie group, L. In practice, to find the allowed transitions, we may compute the H nm by expressing gμ B B 1 (t)S in the basis of eigenstates of H s and compute all possible commutators. By doing so, we can check if L covers the full Hilbert space. In a finite dimensional basis, this is easily computed since the nonzero commutators are of the form: 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 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 μ 0 H z $ 0:5 À 1 T, the level spectrum consists of nonequidistant levels and that n h jS x n þ 1 j i≠0 for all eigenstates |n〉 (see Supplementary  Fig. 4). From the point of view of the universal operation, this is rather favourable since i H n nþ1 ; H nþ1 nþ2 Â Ã ¼ i n j i n þ 2 h jþ h:c:, and so on. Therefore, every transition can be performed by concatenating commutators.
Universality of the [Gd 2 ] six-qubit system (or qu64it). Putting together several qudits, they must interact as a prerequisite to be universal. The [Gd 2 ] 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 H s and B 1 (t)S, the system is presumably not universal. For example, a parity selection rule can impose that only transitions n h jB 1 t ð ÞS n þ 2 j i≠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 [Gd 2 ]. The isotropic Heisenberg model, H ¼ ÀJS 1 S 2 À gμ B HðS 1 þ S 2 Þ conserves the total spin S. Therefore, a spin dimer described by this model is definitely not universal for an external driving gμ B B 1 (t)S, with S = S 1 + S 2 . Transitions between different total spin states are forbidden. Fortunately, in [Gd 2 ] 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 [Gd 2 ] can, under the appropriate conditions, simulate either a six-qubit processor or two decoupled three-qubit systems.

Data availability
Raw data sets relevant to this publication are freely available via the FATMOLS community at the ZENODO repository, using the link https://zenodo.org/communities/ fatmols-fet-open-862893/about/. These raw data correspond to Figs. 2  Received: 15 June 2020; Accepted: 19 October 2020; 62. Spek, A. PLATON SQUEEZE: a tool for the calculation of the disordered solvent contribution to the calculated structure factors. Acta Cryst. C. 71, 9-18 (2015).