Counteracting dephasing in Molecular Nanomagnets by optimized qudit encodings

Molecular Nanomagnets may enable the implementation of qudit-based quantum error-correction codes which exploit the many spin levels naturally embedded in a single molecule, a promising step towards scalable quantum processors. To fully realize the potential of this approach, a microscopic understanding of the errors corrupting the quantum information encoded in a molecular qudit is essential, together with the development of tailor-made quantum error correction strategies. We address these central points by first studying dephasing effects on the molecular spin qudit produced by the interaction with surrounding nuclear spins, which are the dominant source of errors at low temperatures. Numerical quantum error correction codes are then constructed, by means of a systematic optimization procedure based on simulations of the coupled system-bath dynamics, that provide a striking enhancement of the coherence time of the molecular computational unit. The sequence of pulses needed for the experimental implementation of the codes is finally proposed.


INTRODUCTION
Reliable quantum computation demands the adoption of strategies to protect the information being processed from external noise, i.e., of quantum error correction (QEC) schemes 1 . At the same time, while the ultimate quantum computer is expected to host QEC routines based on abstract, systemindependent error models, the modern pioneering era of noisy intermediate-scale quantum devices 2,3 calls for strategies tailored for the specific physical hardware utilized.
In their essence, QEC algorithms achieve information protection by suitably encoding an elementary two-state computational unit, a logical qubit, into a larger Hilbert space. This permits one to distinguish, and thus detect, different error occurrences without corrupting the information so that it is then possible to retrieve it 4 .
While traditional QEC approaches realize the extra space by exploiting registers of many physical two-level systems, alternative routes to QEC have emerged wherein a logical qubit is hosted by a single many-level system (multi-level or qudit encodings) [5][6][7][8][9][10][11][12][13][14][15][16] . The first advantage of the latter strategy is to prevent an overhead of physical units necessary to implement the code. Also, the manipulation of single or of pairs of logical qubits can be much easier, since they do not require controlling multiple physical objects at once 17 . Moreover, the same multi-level object can provide protection against different types of errors 12,15 .
The most important error in molecular spin systems at low temperature is given by pure dephasing, that is, the suppression of coherence between different spin states. Such a decoherence mechanism originates principally from the hyperfine interaction of the central (electronic or nuclear) molecular spin with the bath of surrounding nuclear spins [45][46][47] . Except from specific situations 48 , decoherence of a central spin induced by a nuclear spin bath is known to produce non-exponential decay behaviour 19,45,47,[49][50][51][52][53] . This is due to many factors, such as non-negligible entanglement building up between the spin and an evolving bath, the limited number of nuclear spins (~10 2 ) usually surrounding the molecular spin S, the slow relaxation timescales of the bath relative to the motion of the central spin. Although mandatory to design targeted codes, a QEC framework that takes into account both the multi-level nature of a spin S larger than 1/2 and the explicit structure and dynamics of the nuclear spin bath is still missing.
In this work, we develop a class of numerical spin qudit codes which are designed based on a detailed microscopic structure of the environment responsible for errors and which provide a strongly enhanced correction efficiency. Moreover, the sequence of control pulses and measurements needed for an experimental implementation of such codes is discussed. While the advantage of these codes as compared to a simple spin 1/2 is evident already using small S, the performance strikingly improves as qudits with larger spin are used, thus positively exploiting more and more available levels in the molecular spectrum.
The codes are derived by first analyzing the decoherence effects experienced by a qudit spin S > 1/2 embedded into a realistic nuclear spin bath, by means of numerical simulations of the coupled qudit-bath dynamics through a cluster-correlation expansion (CCE) 54 . A systematic procedure is then put forward to capture the spin-dephasing process by means of error operators acting on the system, which is then used to derive optimized codewords for QEC. Thanks to the flexibility of the procedure, the numerical codes can be optimized taking into account the specific timescale of free evolution admitted between two subsequent QEC cycles, thus allowing one to largely reduce the number of correction steps sufficient to ensure a desired fidelity. As such, they are an optimal candidate for realizations in near-term devices, in which the implementation of the QEC can be noisy and can take relatively long times.

Physical system and decoherence mechanisms
The system analyzed in the following is a molecular electronic spin S (sketched in Fig. 1), described by the Hamiltonian Here, fŜ x ;Ŝ y ;Ŝ z g are spin S > 1 operators, with eigenstates ofŜ z being defined byŜ z m j i ¼ m m j i. The parameter D indicates the zero-field splitting (assumed to be axial, for simplicity) and Ω = g z μ B B z , with g z the longitudinal g-factor and μ B the Bohr magneton, characterizes the Zeeman interaction with a static magnetic field along the z-direction. The analysis developed here can apply both to a qudit given by a single spin-S ion and to a giant spin originating from the strong exchange interactions between different ions 55 . Also, while we focus here on the case of an electronic qudit, the same treatment can also be applied to a nuclear qudit with small adaptations commented in "Methods".
For a molecular electronic spin S at low temperature, the hyperfine coupling with the surrounding nuclear spin bath is the dominant source of decoherence. Indeed, as typically done in quantum computing platforms, we assume to work at temperatures much smaller than the relevant energy scales of the qudit (Ω; Dk À1 B $ K), such that the thermal population of the excited states is negligible. In these conditions, phonon absorption is very unlikely. At the same time, the qudit energy gaps are much smaller than the Debye energy (typically in the ≳50 K range), thus making also resonant phonon emission (whose probability scales as the third power of the gap) negligible 56 . In general, the effect of spin-phonon coupling on the system dynamics can be calculated from diagonalization of the rate matrix 56 , yielding a decay of both diagonal and off-diagonal elements of the system density matrix (associated to relaxation and dephasing, respectively) on related time-scales. In particular, phonon-induced dephasing is limited by the relaxation time T 1 . A body of experiments on molecular spin qubits and qudits 29,40,57 demonstrate that this is not the case at low temperature, where T 1 increases exponentially and becomes several orders of magnitude longer than the dephasing time. Hence, at low temperature phonons are not responsible for pure dephasing on the spin system and other mechanisms come into play.
Spin dipole−dipole interactions between electronic spins can have an important effect in concentrated samples, but this is strongly reduced if the magnetic centres are diluted in a diamagnetic matrix 21 and are not relevant here because we consider a perspective architecture consisting of a single (or a few) molecular unit 46,58 .
We, therefore, focus on the bath B of nuclear spins surrounding S. The number of nuclear spins in the bath may range from a few tens to a few hundreds in realistic molecular complexes, thus being rather far from the infinite-bath limit underpinning typical Markovianity approximations. By working in the so-called 'coherence window' 59 , in which the system energy gaps are much larger than the gaps of the nuclear spin bath, off-diagonal operators inducing population transfer on the system can be neglected. The system-bath dynamics can be described in this regime by effective spin Hamiltonians featuring only a diagonal coupling between S and the bath, which are derived via perturbation theory. This type of Hamiltonians has been studied in the context of a (pseudo)spin S = 1/2 interacting with a nuclear spin bath 6,45,60 . In "Methods", we deduce an effective Hamiltonian for the dynamics of a generic spin S > 1/2. In interaction picture with respect toĤ S and to first order in Ω −1 , this Hamiltonian is of the form Both the intrinsic and the qudit-conditioned Hamiltonians H ðkÞ B of the bath can be written in the general form where fÎ In a free-decay experiment, the main decoherence process is given by inhomogeneous broadening 49 . The system-bath diagonal coupling a ð1Þ nŜzÎ z k has the effect of a classical random magnetic field-the Overhauser field-on S. Uncertainty in the actual bath state then produces, for the density matrix ρ S ðtÞ of the qudit, a Gaussian decay for the single transition coherence, with (see "Methods"), over timescales much faster than those set by the nuclear spin−spin interaction. This is shown in Fig. 2 The dramatic effect of inhomogeneous broadening on the spin coherence is routinely compensated for in experiments by means of spin-echo/refocusing schemes, whose basic example is the Hahn echo. For different qudit spins, the echo dynamics is shown in Fig. 2b. The realization of the echo transformations is further described in 'Practical Implementation'. The coherence decays Fig. 1 Model system. A spin S larger than 1/2, whose many-level structure is exploited for performing multi-level encodings, interacts with the bath B of surrounding nuclear spins. Entanglement between S and B induces spin dephasing, which is counteracted through quantum error correction. Nuclear spins are plotted from a sample configuration of 100 nuclear spins used in the simulations, generated randomly within a sphere of radius 15 Å and with a minimal distance of 3 Å between spins (see "Methods"). over timescales of tens-to-hundreds of microseconds, signalling that the effect of inhomogeneous broadening is removed to a large extent. The decay is now due to the quantum dynamics of the bath, and is mainly determined by the contribution given by intra-bath interactions in H Optimized qudit encoding While increasingly sophisticated echo pulse sequences can recover the effect of inhomogeneous broadening to a better and better degree, the spin coherence remains irremediably affected by the interacting quantum dynamics of the bath. We derive in the following qudit QEC codes as a means to protect the system from these effects. In particular, we develop a framework for designing optimal numerical codes, which are based on the detailed description of the system-bath dynamics adopted in this work.
A QEC code can be defined by following two fundamental steps. The first step is to identify error operatorsÊ k which describe the effect of the noise source on the system, i.e., such that the state of S at time t can be related to the initial state througĥ The second step is the derivation of computational states 0 L j i and 1 L j i that satisfy Knill−Laflamme conditions for QEC 62 , namely, for all k and j, These conditions demand that, when errorsÊ k affect an initial state ψ L j i ¼ α 0 L j i þ β 1 L j i, the codewords are modified but the corresponding coefficients α and β are not, and the information they carry is thus preserved. Moreover, errors do not create ambiguity between 0 L j i and 1 L j i: error wordsÊ k w L j i span subspaces that are orthogonal to each other for different w = 0, 1.
If these conditions are fulfilled, then different errorsÊ k can be distinguished, detected, and corrected. In practice, a measurement is devised whose outcome discriminates the error occurrence and a recovery operation restores the initial state. Importantly, the 2S + 1 levels of a spin S offer enough state space to detect and correct a number N corr = ⌊S⌋ of error operatorsÊ k 13 , where ⌊S⌋ indicates the largest integer smaller or equal than S. Given that a decomposition of the form (4) involves, in general, a larger number ofÊ k , it is essential to identify the N corr error operators which have a stronger effect, such that the code can be tailored for them ensuring optimal correction.
For the spin-dephasing scenario considered here, a decomposition of the form (4) with exact error operators is not known, thus preventing a derivation of adequate codewords for this type of noise. In order to overcome this limitation, we introduce an iterative numerical optimization procedure which, givenρ S ð0Þ and ρ S ðtÞ computed through CCE, aims at determining a number N corr of operatorsÊ k by decreasing contribution toρ S ðtÞ. Starting from ρ At the n-th step, the distance kρ S ðtÞ À ρ ðnÀ1Þ S k (here, ∥ ⋅ ∥ is the Hilbert−Schmidt norm) is numerically minimized in order to find optimal parameters for a parametrized form orÊ n (specified in the following), and the outcome is used for the subsequent step of the iteration. If a hierarchy ofÊ k exists, a successful optimization will return error operators which give less and less contribution, in the norm, to the density matrix. In this sense, the numerical procedure then leads to an optimal decomposition ofρ S ðtÞ in the form (4).
From the structure of the system-bath HamiltonianĤ of Eq. (1), it follows that ρ S ðtÞ can be generically written in the form (4) with error operators that are diagonal in the basis of states m j i. Given that the Hilbert space of S is finite with dimension 2S + 1, the error operatorsÊ k can be expanded onto a basis fD m g of 2S + 1 diagonal operators. Indeed, a diagonal matrix in this space can have at most 2S + 1 linearly independent entries. This justifies an expansion for the error operators of the form The coefficients E k,m are the free parameters for the numerical optimization. The basis fD m g is chosen in the following to be given by the projectorsD m ¼ m j i m h j over the m j i states. Once the error operators are found, the codewords enabling their QEC are determined by imposing Knill−Laflamme's conditions (5) numerically, as detailed in "Methods". The codewords obtained from this procedure are depicted in Fig. 3a for values of spin from S = 3/2 to S = 9/2. By construction, 0 L j i and 1 L j i have support on different subsets of m j i states in an alternate fashion. This automatically guarantees the fulfilment of Knill−Laflamme's condition (6), reducing the number of free parameters for the numerical search required to impose the condition (5).
The performance of the optimized qudit codes is remarkable, as shown in Fig. 3, where the fidelity after the QEC is reported, for a set of codewords corresponding to different qudit spin S. In particular, for each time t, Fig. 3b represents the squared fidelity of the recovered state with respect to the encoded state, achieved by performing an instantaneous QEC at time t. In panel 3c, we report instead the infidelity 1 À F 2 S ðtÞ in log−log scale for an inset of panel (b). One can observe that, while a squared fidelity above 0.9 is maintained for a spin 1/2 only up to~30 μs, the recovered fidelity is above that value for as long as~40, 65, 240, and 300 μs for qudit spin S = 3/2, 5/2, 7/2 and 9/2, respectively. Similarly, the same spin values guarantee a recovered fidelity above 0.99 for up to~15, 20, 29, and 37 μs, well longer than the spin 1/2,~10 μs. The possibility to recover high fidelity even after rather long evolution times is a , where 0 L j i and 1 L j i are spin-binomial codewords corresponding to spin S 13 (see also "Methods"); the results have been averaged over 2 6 initial spatial configurations of the nuclear spins, as explained in "Methods". b Decay of the echo squared fidelity for the same initialization of (a). For each time point t, a generalized echo transformation of the form e ÀiπŜx at t/2 is understood.
crucial resource for near-term implementations. Indeed, if the rate at which subsequent QEC cycles need to be done is too large, the advantage of the correction may get lost in a realistic implementation because of the non-negligible time necessary to implement all the measurements and recovery operations of the QEC step.
The substantial advantage in increasing the spin of the qudit is further emphasized in Fig. 4, where the gain with respect to the spin 1/2, is reported as a function of time, for different values of the spin S. A remarkable maximal gain, larger than 10, is attained, e.g., for a spin 7/2 at around 60 μs, and a maximal gain around 15 is attained for S = 9/2. Depending on the time at which the optimization is performed, rather different numerical codewords can be obtained, reflecting the interplay of different interaction scales in the system-bath Hamiltonian. However, broad temporal windows can be recognized, in which the codewords maintain essentially the same structure, while being quite different in two different regimes. Therefore, a given set of codewords maintains a stable performance if the QEC is implemented in a rather broad time interval around the optimization time. These features can be observed in Fig. 5a, where the gain as a function of time is shown for three different codeword pairs obtained by optimizing at times t opt = 10, 50, 100 μs for S = 3/2. As expected, while a unique pair giving the largest gain at all times cannot be found, codewords optimized at a given time provide very good performance in a rather large region around that time. We have finally checked that the performance of the numerically optimized codewords does not critically depend on the initial state chosen for our procedure, as demonstrated in Fig. 5b for the exemplary case of S = 3/2. The squared fidelity (colour scale) as a function of time (radial scale) is reported for different values of the angle θ (angular scale) characterizing an initial state of the form cosðθÞ 0 L j i þ i sinðθÞ 1 L j i. Large fidelities are attained for all initial states, with fidelity increasing as one departs from the state with equal weights [θ = π/4], that is, the most decoherence-prone state, used in the other reported simulations.

Practical implementation
Having analyzed the ideal efficiency of the optimized qudit codes, we now turn to a discussion of how to implement in practice all the steps of the QEC procedure, namely encoding, detection, and recovery, whose formal description is given in "Methods". The manipulation of a spin S > 1/2 system requires of course more complex control sequences compared to a spin S = 1/2. Nonetheless, it can be realized in a total time sufficiently short so that it does not significantly impact the efficiency of the ideal QEC, as detailed in the following.
The population transfers required among spin states can be realized using sequences of resonant microwave/radiofrequency pulses. These are described in the Hamiltonian by time-dependent control fields of the form g y μ B B y ðtÞ cosðωtÞŜ y where the envelope B y (t) is typically rectangular or Gaussian. These pulses induce transitions from a spin state m j i to m ± 1 j i when ω is set at the corresponding transition frequency, implementing a two-state unitary rotation Y m ðθÞ ¼ exp½θ=2ð m þ 1 j im h j À m j i m þ 1 h jÞ between states m j i and m þ 1 j i of an arbitrary angle θ. All the and numerical codewords 0 L j i and 1 L j i optimized at t = 5 μs. c Infidelity 1 À F 2 S for an inset of panel (b) representing the region of F 2 S ! 0:9. Fig. 4 Gain, defined in Eq. (9), with respect to non-corrected spin 1/2 for the data in Fig. 3. The curves indicate the average of G S ðtÞ over the different nuclear spin spatial configurations (generated as explained in "Methods"), while the shaded areas mark the region included between G S ðtÞ ± σ S ðtÞ, where σ S ðtÞ ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi hG 2 S ðtÞi À hG S ðtÞi 2 q is the standard deviation and 〈 ⋅ 〉 denotes averaging. The shaded areas show that a large gain is obtained for all spatial configurations at given spin S. The standard deviation σ S (t) increases with S since the denominators in Eq. (9) become smaller and smaller for increasing S. A comparison with spin-binomial codes is further given in "Methods".  steps of the QEC are illustrated in the following, and the explicit realization for a S = 5/2 qudit code is depicted in Fig. 6. The implementation proposed here generalizes the one proposed for spin-binomial codes 13 .
Encoding. We assume that the information, i.e., the coefficients α and β, is initially encoded in a simple state such as α À1=2 j iþ β 1=2 j i. The preparation of the logical state α 0 L j i þ β 1 L j i for arbitrary α and β is then realized by alternating pulses Y m (θ) distributing population among the different m j i states and π-pulses Y m ( ± π) which rearrange the different populations in the correct order. The angles θ of the two-level rotations are fully determined by the components of the code-words on each m j i state. The explicit sequence for S = 5/2 is given in panel 'encoding' of Fig. 6a with the angles given in Table 1. Once the state has been codified, the system is left to decay freely for a time t/2, then a spin-echo pulse sequence is performed, and the QEC finally takes place at time t starting with the error detection. These pulse sequences can also be used to perform single-qubit gates between the code-words, for instance by first re-mapping the code-words to, e.g., ± 1=2 j i states, performing the desired twostate operation, and re-encoding.
Spin echo. The Hahn echo for a spin 1/2 can be understood as a magnetic pulse along x or y which effectively flips the spin. Then, the spin can be viewed as effectively evolving withĤ for a time t/2 and with the same Hamiltonian but withŜ z changed to ÀŜ z for an equal time t/2, where t is the time at which the QEC is performed. Similarly, the echo scheme is extended here to a larger spin S > 1/2 by considering a 'generalized-pulse' transformation of the form U echo ¼ e ÀiπŜx at time t/2 which inverts the spin, sending state m j i to Àm j i. This transformation can be realized with a sequence of S + 1/2 (for half-integer S) resonant π-rotations along x or y, coupling pairs of m j i $ Àm j istates, followed by Y m (±π) pulses to rearrange populations. For instance, in the case of a spin 5/2, U echo is obtained by three independent rotations 5=2 j i $ À5=2 j i, 3=2 j i $ À3=2 j i and 1=2 j i $ À1=2 j i. Due to the lack of a direct matrix element between m j i and Àm j i in the architecture considered here, each of these Δm > 1 transformations needs to be decomposed into Δm = ±1 transitions. This is done, for instance, using the strategy discussed in "Methods-Pulse sequences".
Detection. To realize the error detection, two additional ingredients are introduced in the system considered until now. The first one is a weak coupling of the qudit to a spin s A = 1/2 ancilla. The ancilla is described by adding to the Hamiltonian of Eq. (12) the terms where fŝ x A ;ŝ y A ;ŝ z A g are spin-1/2 operators for the ancilla. The first term in Eq. (10) describes the Zeeman coupling of the ancilla to the static magnetic field whilst the second one describes the ancilla-qubit coupling parametrized by the tensor J. For J x;y ( jΩ À Ω A j, such that states of qudit and ancilla remain essentially factorized, the excitation energies of the ancilla Δ . Blue and orange colours represent amplitudes associated to α and β (i.e., 0 L j i and 1 L j i), respectively. To exemplify the effect of the basis rotation and the outcome of the measurement during the detection, the first-error case (related to the error operatorÊ 1 of Eq. (8)) is shown in (b). Table 1. Angles for the pulse sequence realizing the QEC for S = 5/2 as depicted in Fig. 6.
J zŝ z AŜ z only, i.e., Δ By irradiating the ancilla with a resonant magnetic pulse at angular frequency Δ ðmÞ A it is thus possible to selectively excite the ancilla only if the qudit is in state m j i. A subsequent measurement of the state of the ancilla then reveals whether the qudit state has support on m j i or not. This mechanism will be exploited in the following to detect the different possible errors without corrupting α and β. Apart from this selective excitation immediately followed by a measurement, the ancilla is always in its ground state, and thus it does not affect the previously developed treatment of the qudit incoherent dynamics. For this reason, its coupling to the nuclear spins is also irrelevant for the present discussion.
The second ingredient is a coupling of the magnetic molecule to a microwave resonator. Crucial steps towards achieving the strong coupling between magnetic molecules and a resonator have been experimentally demonstrated recently 63 . This coupling can then be exploited to measure the ancilla, building on techniques well developed in the field of circuit quantum-electrodynamics [64][65][66] and adapted to Molecular Nanomagnets 67,68 . The coupling of the molecule to the resonator induces a shift of the resonance frequency of the resonator which depends on the ancilla-qudit state. As explained below, this can be exploited to measure the state of the ancilla without collapsing the qudit state.
The error detection is described in an abstract setting by a projective measurement with the projector operators given in Eq. (25) of "Methods". Given the difficulty to implement similar operators, that project into complex superpositions of system eigenstates, the detection is divided into two steps. In the first step, a sequence of pulses is performed which rotates the full basis of error words into the basis of m j i states in both error spaces corresponding to 0 L j i and 1 L j i (see panel 'basis rotation' of Fig. 6 for S = 5/2) 13 . This operation thus converts the detection of the projectors (25) into an easier-toimplement measurement in the m j i basis, and the unitarity of the transformation ensures that α and β are preserved. At this point, every possible post-error state is of the form α m j i þ β m 0 j i for different pairs of states ð m j i; m 0 j iÞ. We now aim to induce an excitation of the ancilla only for a superposition state of the qudit with components on ð m j i; m 0 j iÞ, without collapsing the superposition. We achieve this by applying a two-tone two-photon drive at frequencies Δ  69 . Then, in the dispersive limit, the coupling G between resonator and ancilla induces a shift of the cavity angular frequency ω c of ±G 2 /δ m , with δ m ¼ Δ ðmÞ A À ω c and the sign of the shift depending on the state of the ancilla 64 . Since here we need to measure the ancilla irrespective of the state of the qudit in the subspace ð m j i; m 0 j iÞ, a frequency-independent measurement of the state of the ancilla must be performed. Two different approaches to solve this same issue, by detecting the amplitude (but not the frequency) of the output field, are described in 69 . The ancilla is then measured by exploiting its coupling to the resonator and the qudit wavefunction is projected onto such states. The sequence of measurements is then repeated probing each (mutually exclusive) pair of ð m j i; m 0 j iÞ states sequentially, returning a yes/no answer at each step if the system is found in the corresponding error state, and stopping if a positive outcome is obtained. Hence, there will be at most ⌊S⌋ measurements given that the number of possible errors for a qudit of spin S is ⌈S⌉, where ⌈S⌉ (⌊S⌋) indicates the smallest integer larger (largest integer smaller) or equal than S.
Recovery. After detection, the system has been projected into a superposition state of the form α m j i þ β m 0 j i with known m and m 0 . The simplest way to restore the encoded state α 0 L j i þ β 1 L j i is then to first use a few π-pulses to send m j i ! 1=2 j i and m 0 j i ! À1=2 j i, and then to repeat the pulse sequence which implements the encoding (panel 'recovery' in Fig. 6). Alternatively, one can save a few pulses by redesigning the encoding sequence starting from each possible pair m j i, m 0 j i resulting from detection.
Impact on performance. The non-instantaneous duration of the QEC procedure in a realistic implementation (during which information is not protected), together with related potential imperfections, may cause a loss of efficiency in the correction. We thus here discuss to what extent such effects can reduce the expected performance.
The operations described to implement the QEC involve sequences of resonant pulses (and ancilla measurements) only. In electronic spin systems, a single π-pulse requires less than 10 ns for achieving a state transfer with high fidelity, and this time could be further reduced by pulse-shaping techniques [70][71][72] . The measurement time for the ancilla readout through a microwave resonator can be roughly estimated from the field of circuit QED to be of 40−100 ns with fidelity above 0.98 65,73 . Then, for the spins S ≤ 9/2 considered here the QEC procedure requires a total time ranging from a few hundreds of nanoseconds to few microseconds at most, and is hence much shorter than the decay time that can be allowed by the optimized code-words while ensuring a recovery fidelity above 0.99. Indeed, the latter can be of tens of microseconds, as visible from Figs. 3 and 5.
The practical implementation of the QEC is thus expected not to significantly reduce the correction performance for the qudits studied here. However, one could also predict that the growth of the complexity of the implementation for very large spins will eventually set a tradeoff between gain and duration of the QEC favouring the use of moderately large spins, similarly to what was observed for spin-binomial codes 13 . Importantly, this limitation can be mitigated in the present scheme by optimizing the codewords at larger times. Moreover, it should be noted that the bottleneck in the specific experimental implementation proposed here for large spins is related to the rapid scaling of the number of pulses required with S. This, in turn, is due to the low connectivity of the 2S + 1 spin levels that permits resonant state transfers only between states with Δm = ±1. A possible way around this problem is to consider magnetic molecules with competing interactions 46,[74][75][76][77] , for which the multi-level structure used for the qudit encoding is given by low-energy states belonging to different multiplets that can provide larger state connectivity.

DISCUSSION
We have investigated decoherence effects produced by a realistic nuclear spin bath on a spin qudit S > 1/2 in Molecular Nanomagnets, by simulating the coupled system-bath dynamics via a CCE. Building on this analysis, we have developed a versatile numerical strategy to construct optimal QEC codes tailored for the specific spin-dephasing errors induced by the bath, thus bridging the gap between traditional general-purpose correction algorithms and the necessity of hardware-specific strategies meeting current experimental capabilities. The resulting qudit codes achieve a remarkable performance, and can be optimized by taking into account constraints on the time interval between subsequent QECs. Moreover, the increase in performance with the increase of the qudit spin is striking, signalling that the codes exploit the available levels of the molecular system as a resource very efficiently. The proposed codes can be implemented experimentally using standard sequences of resonant control pulses. Such sequences are explicitly designed and discussed, and their practical realization is predicted not to set important limits on the efficiency of the codes. Given these results, the proposed codes are a promising candidate for realizing error-protected quantum computational units embedded at the single-molecule level, a central building block for implementing reliable quantum information processing on short-term molecular devices.
Recent works 47 point out that the CCE method used here correctly reproduces the phenomenology of coherence enhancement due to the existence of a nuclear diffusion barrier 53 . An interesting perspective is the integration of the framework developed in this work with chemical engineering techniques for achieving an even longer lifetime of the error-corrected logical qubit through this mechanism. The synergy of tailored QEC codes as investigated here with engineered nuclear spin distributions may pave the way towards a class of highly coherent molecular qubits.
The framework developed in this work is general, and can be applied to a wide landscape of molecular systems and also to other individual spin systems such as impurities in semiconductors, in order to design a proof-of-principle experiment to demonstrate the effectiveness of the QEC code. In addition, it can be extended, in the future, to investigate decoherence effects affecting superpositions of spin states belonging to different spin multiplets 74,75 . This would be interesting since, on the one hand, the use of many low-m spin states belonging to different multiplets may allow one to increase the number of levels available for error correction without exasperating dephasing effects given by large m À m 0 transitions. On the other hand, it would enable a thorough study of standard block encodings embedded in a single molecule, wherein a register of qubits is achieved through many effective spin-1/2 systems selected from different spin multiplets.

Derivation of the effective Hamiltonian
The spin Hamiltonian describing the interacting evolution of the molecular spin S and the bath B of N nuclear spins iŝ The tensors D n contain dipole−dipole couplings between S and B, while the tensors E n;m contain nuclear−nuclear dipolar couplings. The elements of D n satisfy and the same holds for the corresponding elements of E n;m . A canonical perturbative (Schrieffer-Wolff) transformation generated by with h.c. indicating the hermitian conjugate, is constructed such that H G ¼ e GĤ e ÀG , within first order in Ω −1 , does not contain off-diagonal couplings between S and B with respect to the states m j i 45,48,60,78 . As detailed in the following, G is proportional to Ω −1 and leading orders inĤ G can thus be computed from the Baker−Campbell−Hausdorf expansion The coefficients G αβ k are determined explicitly by imposing the cancellation of the off-diagonal couplings between S and B to first order. This results in the relation The coefficients G αβ k then read These expressions are indeed of order Ω −1 and the transformation generated by G is thus perturbative, such that its effect on the initial factorized qudit-nuclei state is neglected. By keeping terms in Eq. (15) to first order in Ω −1 only and neglecting energy-non-conserving terms, the effective Hamiltonian of Eqs. (1) and (2) is obtained with coefficients n;m ¼ E þÀ n;m ; d ð0Þ n;m ¼ E zz n;m =2; a ð1Þ n ¼ D zz n ; b ð1Þ n ¼ 2 Ω jD þz n j 2 À jD þþ n j 2 À D þÀ and b ð2Þ n ¼ c ð2Þ n ¼ d ð2Þ n ¼ 0. The energy of S is also renormalized according tõ but this is absorbed into the interaction picture in Eq. (1).

Simulations and dephasing timescales
The configuration of nuclear spins in space is generated randomly within a sphere of radius 15 Å from the spin S, as sketched in Fig. 1. Further, a minimal distance of 3 Å is considered between nuclear spins and between each nuclear spin and S. In all the simulations presented in this work, configurations of N = 100 nuclear spins are considered, whose initial state is taken to be thermal at temperature T = 2 K. Moreover, a static magnetic field B z = 1 T along z is assumed, for achieving the regime of large Zeeman energy of S. The decoherence function, is computed by means of a CCE 54,61 . This expansion decomposes L nm (t) as a product of contributions originating from clusters involving an increasing number of nuclear spins, and is formally equivalent to a perturbative expansion in small intra-bath effective couplings. Clusters involving more and more spins contribute smaller and smaller corrections to L nm (t), justifying a truncation of the expansion to few-spin clusters for practical applications. For inclusion up to n-size clusters, we call this truncation CCEn, which yields the truncated function L ðnÞ nm ðtÞ. The effect of inhomogeneous broadening is well captured by CCE-1 61,79 . Given that nuclear gaps are of the order of millikelvin in magnitude, an initial thermal state of the nuclear spin-bath at temperatures T of a few kelvins is to a good approximation a fully unpolarized state Under this approximation, the CCE-1 can be solved analytically also for S > 1/2 giving Here, D zz k is the hyperfine coupling /Ŝ zÎ z k between S and the k-th nuclear spin. For small D zz k t, this is well approximated by e ÀðnÀmÞ 2 Γ 2 t 2 with Γ 2 ¼ P k ðD zz k Þ 2 =4 as given in the main text. For the echo dynamics, we find that CCE-2 gives essentially converged results, as tested by including also larger clusters. For this reason, the numerical results reported here are obtained using CCE-2. This convergence confirms that intrinsic nuclear flip-flop processes give the dominant contribution to spin dephasing after echo. Indeed, an inspection of the magnitude of the coefficients a ðkÞ n;m , b ðkÞ n;m , c ðkÞ n;m , d ðkÞ n;m reveals three fundamental interaction scales at play, which ordered by decreasing strength, are associated to (i) the diagonal coupling between S and nuclear spins (terms /Ŝ zÎ z k which are compensated for by the echo), (ii) the intrinsic interacting evolution (terms /Î α nÎ β m ), (iii) the S-conditioned interacting evolution (terms /Ŝ zÎ α mÎ β n ). These different energy scales are responsible for contributions to decoherence over different timescales, with (ii) being dominant in the echo decay.

Nuclear spin qudit
The quantitative analyses presented in this work are focused on the case of an electronic spin qudit. Nevertheless, the theoretical framework applies also to the case of a nuclear spin qudit. The crucial approximation underlying the physical model studied here is that the qudit energy gaps are much larger than the gaps of surrounding nuclear spins of the bath. For an electronic qudit, these energy differences can be made intrinsically large by using a sufficiently large static magnetic field. For a nuclear spin qudit of a magnetic ion (coupled to an electronic spin), whose interaction with surrounding nuclear spins is mediated by virtual excitations of the electronic spin, the necessary energy difference mainly originates from contact hyperfine interaction between the nuclear and electronic spin. The construction of the effective Hamiltonian, and the hierarchy of the interaction scales at play, then follows as discussed above.

Derivation of numerical qudit codes
The iteration defined by Eq. (7) is first used to determine the error operators given in Eq. (8). The numerical codewords 0 L j i and 1 L j i are found by starting from the following ansatz, inspired by spin-binomial codes 13 , This ansatz permits one to impose Eq. (6) by construction and hence to reduce the number of free parameters, thanks to 0 L j i and 1 L j i having nonzero components on different sets of m j i states. Knill−Laflamme conditions (5) are finally enforced on the coefficients γ ð0Þ=ð1Þ ℓ by numerically minimizing the function

Detection and recovery
The abstract detection and recovery operations follow the general QEC theory of ref. 62 . Once a set of operators fÊ k g k¼0; ;bSc to be corrected with a spin S is identified, the two error subspaces corresponding to 0 L j i and 1 L j i are first determined. These two subspaces are defined as the linear span of the statesÊ k 0 L j i andÊ k 1 L j i for all k, respectively. For each of these subspaces, a basis f e that outcome to the computational states 0 L j i and 1 L j i, respectively. Since the coefficients α and β of the encoded state have been preserved, this operation then fully restores the logical state ψ L j i. The recovery is formalized by a set of transformations fR k g such that for c = 0, 1. Here, each transformationR k is constructed as a rotation in the two-dimensional space spanned by e ðcÞ k E and c L j i.

Pulse sequences
To systematically convert a generic transformation U acting on the state space of a spin S > 1/2 into a sequence of resonant Y m (θ) pulses, one can exploit known decomposition strategies from quantum control theory 13,80 .
In a first step, the unitary U can be decomposed into a sequence of twostate planar rotations using the algorithm given in ref. 80 . This gives a sequence of two-state transfers which involve states m j i and m 0 j i also with jm À m 0 j>1. To further convert such two-state rotations into rotations Y m (θ), i.e., with jm À m 0 j ¼ 1, one finally exploits π-pulses to bring the population of m 0 close to m and back. For instance, defining Y m;m 0 ðθÞ ¼ exp θ=2 À m j i m 0 h j ð þm 0 j i m h j ½ (27) for m 0 > m, one can iteratively exploit the properties Y m;mþ2 ðθÞ ¼ Y m;mþ1 ðπÞY mþ1;mþ2 ðÀθÞY m;mþ1 ðÀπÞ: ¼ Y m;mþ1 ðÀπÞY mþ1;mþ2 ðθÞY m;mþ1 ðπÞ: As an example, the pulse sequence depicted in panel 'basis rotation' of Fig.   6a, which realizes a transformation mapping the basis of error words e ðcÞ k E to the basis of m j i states, is obtained with the procedure sketched here. The resulting angles are given in Table 1.

Spin-binomial codes
In order to compare the numerical codes with other qudit approaches to spin dephasing, we test recently-proposed spin-binomial codes 13 . These codes are based on a description of spin dephasing as produced by a Markovian bath which couples to the system via operator ffiffi ffi γ pŜ z . The latter model can only describe an exponential decay of coherence with rate γ and contributions of order (γt) n to the density matrix can be computed exactly and are determined by powers up toŜ n z . We find that, while spinbinomial codes still give an interesting performance, they are largely overwhelmed by the numerical codes. This can be appreciated by comparing Figs. 3 and 4 with Fig. 7, both in terms of fidelity and gain. The fact that spin-binomial codes still give an increasing gain for increasing qudit spin despite being designed for a simpler dephasing mechanism, suggests that small powers of the coupling operatorsŜ z play a fundamental role in the decoherence process also in the present scenario, over short timescales with respect to the intra-bath interaction strength.

DATA AVAILABILITY
The datasets generated during and/or analyzed during the current study are available from the corresponding author on reasonable request. Fig. 7 Performance of spin-binomial codes. a Average squared fidelity as a function of time for different values of the qudit spin S, for initial state ψ L j i ¼ ð 0 L j i þ i 1 L j iÞ= ffiffi ffi 2 p and spin-binomial codewords 0 L j i and 1 L j i. The results are averaged over 2 6 nuclear spin configurations. b Gain, defined in Eq. (9), with respect to a spin 1/2, for the same data of panel (a). The shaded area are constructed as in Fig. 4.