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) schemes1. At the same time, while the ultimate quantum computer is expected to host QEC routines based on abstract, system-independent error models, the modern pioneering era of noisy intermediate-scale quantum devices2,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 it4.

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 once17. Moreover, the same multi-level object can provide protection against different types of errors12,15.

A very promising architecture for the implementation of multi-level encodings is given by molecular nanomagnets13,16,18. Indeed, these highly coherent systems19,20,21,22,23,24,25,26,27 offer many accessible spin levels28,29,30,31, which can be manipulated with high accuracy through electromagnetic pulses, and they can be chemically engineered to meet desired requirements32,33,34,35,36,37,38,39,40,41,42,43,44.

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 spins45,46,47. Except from specific situations48, decoherence of a central spin induced by a nuclear spin bath is known to produce non-exponential decay behaviour19,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 (~102) usually surrounding the molecular spin \({{{\mathcal{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.

Results

Physical system and decoherence mechanisms

The system analyzed in the following is a molecular electronic spin \({{{\mathcal{S}}}}\) (sketched in Fig. 1), described by the Hamiltonian \({\hat{H}}_{{{{\mathcal{S}}}}}=D{\hat{S}}_{z}^{2}+{{\Omega }}{\hat{S}}_{z}\). Here, \(\{{\hat{S}}_{x},{\hat{S}}_{y},{\hat{S}}_{z}\}\) are spin S > 1 operators, with eigenstates of \({\hat{S}}_{z}\) being defined by \({\hat{S}}_{z}\left|m\right\rangle =m\left|m\right\rangle\). The parameter D indicates the zero-field splitting (assumed to be axial, for simplicity) and Ω = gzμBBz, with gz 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 ions55. 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”.

Fig. 1: Model system.
figure 1

A spin \({{{\mathcal{S}}}}\) larger than 1/2, whose many-level structure is exploited for performing multi-level encodings, interacts with the bath \({{{\mathcal{B}}}}\) of surrounding nuclear spins. Entanglement between \({{{\mathcal{S}}}}\) and \({{{\mathcal{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”).

Full size image

For a molecular electronic spin \({{{\mathcal{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 (\({{\Omega }},D{k}_{{{{\rm{B}}}}}^{-1} \sim\) 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) negligible56. In general, the effect of spin-phonon coupling on the system dynamics can be calculated from diagonalization of the rate matrix56, 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 T1. A body of experiments on molecular spin qubits and qudits29,40,57 demonstrate that this is not the case at low temperature, where T1 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 matrix21 and are not relevant here because we consider a perspective architecture consisting of a single (or a few) molecular unit46,58.

We, therefore, focus on the bath \({{{\mathcal{B}}}}\) of nuclear spins surrounding \({{{\mathcal{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 \({{{\mathcal{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 bath6,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 \({\hat{H}}_{{{{\mathcal{S}}}}}\) and to first order in Ω−1, this Hamiltonian is of the form

$$\hat{H}={\hat{H}}_{{{{\mathcal{B}}}}}^{(0)}+{\hat{S}}_{z}{\hat{H}}_{{{{\mathcal{B}}}}}^{(1)}+({\hat{S}}^{2}-{\hat{S}}_{z}^{2}){\hat{H}}_{{{{\mathcal{B}}}}}^{(2)}.$$
(1)

Both the intrinsic and the qudit-conditioned Hamiltonians \({H}_{{{{\mathcal{B}}}}}^{(k)}\) of the bath can be written in the general form

$${H}_{{{{\mathcal{B}}}}}^{(k)}=\mathop{\sum }\limits_{n=1}^{N}\left({a}_{n}^{(k)}{\hat{I}}_{n}^{z}+{b}_{n}^{(k)}{({\hat{I}}_{n}^{z})}^{2}\right)+\mathop{\sum }\limits_{\begin{array}{c}{n,\!m\!=\!1}\\ m\ne n\end{array}}^{N}\left({c}_{n,m}^{(k)}{\hat{I}}_{n}^{+}{\hat{I}}_{m}^{-}+{d}_{n,m}^{(k)}{\hat{I}}_{n}^{z}{\hat{I}}_{m}^{z}\right),$$
(2)

where \(\{{\hat{I}}_{k}^{x},{\hat{I}}_{k}^{y},{\hat{I}}_{k}^{z}\}\) are spin operators for the k-th nuclear spin of the bath and \({\hat{I}}_{k}^{\pm }={\hat{I}}_{k}^{x}\pm {{{\rm{i}}}}{\hat{I}}_{k}^{y}\). The coefficients \({a}_{n,m}^{(k)},{b}_{n,m}^{(k)},{c}_{n,m}^{(k)},{d}_{n,m}^{(k)}\) are a function of the hyperfine couplings between \({{{\mathcal{S}}}}\) and \({{{\mathcal{B}}}}\), the nuclear spin−spin dipolar couplings, and Ω. In the following, nuclear spins are assumed to be protons (I = 1/2), since hydrogen nuclei typically represent the major decoherence source. Other relevant parameters and details on the simulations are given in “Methods”.

System-bath entanglement, generated by the Hamiltonian of Eq. (1), can be interpreted in terms of ‘which-way information’ accumulated in the state of the nuclear spins: depending on the state \(\left|m\right\rangle\) of \({{{\mathcal{S}}}}\), the bath undergoes different interacting evolutions described by Hamiltonians \({\hat{H}}_{{{{\mathcal{B}}}},m}=\left\langle m\right|\hat{H}\left|m\right\rangle\). These conditioned bath evolutions result in a decay of coherences in the system density matrix \({\rho }_{{{{\mathcal{S}}}}}(t)\), according to \(\left\langle n\right|{\rho }_{{{{\mathcal{S}}}}}(t)\left|m\right\rangle ={L}_{nm}(t)\left\langle n\right|{\rho }_{{{{\mathcal{S}}}}}(0)\left|m\right\rangle\). The function \({L}_{nm}(t)={{{\mbox{tr}}}}_{{{{\mathcal{B}}}}}\left[{{{{\rm{e}}}}}^{-{{{\rm{i}}}}{\hat{H}}_{{{{\mathcal{B}}}},n}t}{\rho }_{{{{\mathcal{B}}}}}(0){{{{\rm{e}}}}}^{{{{\rm{i}}}}{\hat{H}}_{{{{\mathcal{B}}}},m}t}\right]\), with \({\rho }_{{{{\mathcal{B}}}}}(0)\) the initial bath state, is computed numerically in the following through a CCE54,61.

In a free-decay experiment, the main decoherence process is given by inhomogeneous broadening49. The system-bath diagonal coupling \({a}_{n}^{(1)}{\hat{S}}_{z}{\hat{I}}_{k}^{z}\) has the effect of a classical random magnetic field—the Overhauser field—on \({{{\mathcal{S}}}}\). Uncertainty in the actual bath state then produces, for the density matrix \({\rho }_{{{{\mathcal{S}}}}}(t)\) of the qudit, a Gaussian decay for the single transition coherence,

$$\left\langle m\right|{\rho }_{{{{\mathcal{S}}}}}(t)\left|n\right\rangle \approx {{{{\rm{e}}}}}^{-{\left(n-m\right)}^{2}{{{\Gamma }}}^{2}{t}^{2}}\left\langle m\right|{\rho }_{{{{\mathcal{S}}}}}(0)\left|n\right\rangle ,$$
(3)

with \({{{\Gamma }}}^{2}={\sum }_{k}{({a}_{n}^{(1)})}^{2}/4\) (see “Methods”), over timescales much faster than those set by the nuclear spin−spin interaction. This is shown in Fig. 2, where the squared fidelity \({{{{\mathcal{F}}}}}_{S}^{2}(t)\) with respect to the initial state, with \({{{{\mathcal{F}}}}}_{S}(t)={{{\mbox{tr}}}}_{{{{\mathcal{S}}}}}\sqrt{\sqrt{{\rho }_{{{{\mathcal{S}}}}}(t)}{\rho }_{{{{\mathcal{S}}}}}(0)\sqrt{{\rho }_{{{{\mathcal{S}}}}}(t)}}\)4, is depicted for different spin S. The fidelity decays over timescales of hundreds of nanoseconds.

Fig. 2: Inhomogeneous broadening and spin echo.
figure 2

a Decay of the squared fidelity with respect to the initial state under inhomogeneous broadening for different spins S initialized in a state \(\left|{\psi }_{L}\right\rangle =(\left|{0}_{L}\right\rangle +{{{\rm{i}}}}\left|{1}_{L}\right\rangle )/\sqrt{(2)}\), where \(\left|{0}_{L}\right\rangle\) and \(\left|{1}_{L}\right\rangle\) are spin-binomial codewords corresponding to spin S13 (see also “Methods”); the results have been averaged over 26 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 \({{{{\rm{e}}}}}^{-{{{\rm{i}}}}\pi {\hat{S}}_{x}}\) at t/2 is understood.

Full size image

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 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}_{{{{\mathcal{B}}}}}^{(0)}\), of the form \({c}_{n,m}^{(0)}{\hat{I}}_{n}^{+}{\hat{I}}_{m}^{-}\). If the latter were absent, the echo would recover unit fidelity independently from S over timescales of hundreds of microseconds, until virtual flip-flops described by the terms of type \({c}_{n,m}^{(1)}{\hat{S}}_{z}{\hat{I}}_{n}^{+}{\hat{I}}_{m}^{-}\) set in. This effect is still partially visible in Fig. 2b in the fact that for short timescales with respect to interactions in \({H}_{{{{\mathcal{B}}}}}^{(0)}\), the fidelity exhibits almost overlapping decay for different S.

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 \({\hat{E}}_{k}\) which describe the effect of the noise source on the system, i.e., such that the state of \({{{\mathcal{S}}}}\) at time t can be related to the initial state through

$${\hat{\rho }}_{{{{\mathcal{S}}}}}(t)=\mathop{\sum}\limits_{k}{\hat{E}}_{k}{\hat{\rho }}_{{{{\mathcal{S}}}}}(0){\hat{E}}_{k}^{{\dagger} }.$$
(4)

The second step is the derivation of computational states \(\left|{0}_{L}\right\rangle\) and \(\left|{1}_{L}\right\rangle\) that satisfy Knill−Laflamme conditions for QEC62, namely, for all k and j,

$$\left\langle {0}_{L}\right|{\hat{E}}_{k}^{{\dagger} }{\hat{E}}_{j}\left|{0}_{L}\right\rangle =\left\langle {1}_{L}\right|{\hat{E}}_{k}^{{\dagger} }{\hat{E}}_{j}\left|{1}_{L}\right\rangle ,$$
(5)
$$\left\langle {0}_{L}\right|{\hat{E}}_{k}^{{\dagger} }{\hat{E}}_{j}\left|{1}_{L}\right\rangle =0.$$
(6)

These conditions demand that, when errors \({\hat{E}}_{k}\) affect an initial state \(\left|{\psi }_{L}\right\rangle =\alpha \left|{0}_{L}\right\rangle +\beta \left|{1}_{L}\right\rangle\), 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 \(\left|{0}_{L}\right\rangle\) and \(\left|{1}_{L}\right\rangle\): error words \({\hat{E}}_{k}\left|{w}_{L}\right\rangle\) span subspaces that are orthogonal to each other for different w = 0, 1.

If these conditions are fulfilled, then different errors \({\hat{E}}_{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 Ncorr = ⌊S⌋ of error operators \({\hat{E}}_{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 \({\hat{E}}_{k}\), it is essential to identify the Ncorr 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 \({\hat{\rho }}_{{{{\mathcal{S}}}}}(0)\) and \({\hat{\rho }}_{{{{\mathcal{S}}}}}(t)\) computed through CCE, aims at determining a number Ncorr of operators \({\hat{E}}_{k}\) by decreasing contribution to \({\hat{\rho }}_{{{{\mathcal{S}}}}}(t)\). Starting from \({\hat{\rho }}_{{{{\mathcal{S}}}}}^{(0)}\equiv {\hat{E}}_{0}{\hat{\rho }}_{{{{\mathcal{S}}}}}(0){\hat{E}}_{0}^{{\dagger} }\), the n-step estimate \({\hat{\rho }}_{{{{\mathcal{S}}}}}^{(n)}\) to \({\hat{\rho }}_{{{{\mathcal{S}}}}}(t)\) of the iteration is defined according to

$${\hat{\rho }}_{{{{\mathcal{S}}}}}^{(n)}={\hat{\rho }}_{{{{\mathcal{S}}}}}^{(n-1)}+{\hat{E}}_{n}{\hat{\rho }}_{{{{\mathcal{S}}}}}(0){\hat{E}}_{n}^{{\dagger} }.$$
(7)

At the n-th step, the distance \(\parallel {\rho }_{{{{\mathcal{S}}}}}(t)-{\rho }_{{{{\mathcal{S}}}}}^{(n-1)}\parallel\) (here, ∥ ⋅ ∥ is the Hilbert−Schmidt norm) is numerically minimized in order to find optimal parameters for a parametrized form or \({\hat{E}}_{n}\) (specified in the following), and the outcome is used for the subsequent step of the iteration. If a hierarchy of \({\hat{E}}_{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 \({\hat{\rho }}_{{{{\mathcal{S}}}}}(t)\) in the form (4).

From the structure of the system-bath Hamiltonian \(\hat{H}\) of Eq. (1), it follows that \({\rho }_{{{{\mathcal{S}}}}}(t)\) can be generically written in the form (4) with error operators that are diagonal in the basis of states \(\left|m\right\rangle\). Given that the Hilbert space of \({{{\mathcal{S}}}}\) is finite with dimension 2S + 1, the error operators \({\hat{E}}_{k}\) can be expanded onto a basis \(\{{\hat{D}}_{m}\}\) 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

$${\hat{E}}_{k}=\mathop{\sum }\limits_{m = 0}^{2S}{E}_{k,m}{\hat{D}}_{m}.$$
(8)

The coefficients Ek,m are the free parameters for the numerical optimization. The basis \(\{{\hat{D}}_{m}\}\) is chosen in the following to be given by the projectors \({\hat{D}}_{m}=\left|m\right\rangle \,\left\langle m\right|\) over the \(\left|m\right\rangle\) 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, \(\left|{0}_{L}\right\rangle\) and \(\left|{1}_{L}\right\rangle\) have support on different subsets of \(\left|m\right\rangle\) 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).

Fig. 3: Numerically optimized qudit codes.
figure 3

a Absolute value of the overlap of the optimized code-words \(\left|{0}_{L}\right\rangle\) (blue) and \(\left|{1}_{L}\right\rangle\) (orange) with each state \(\left|m\right\rangle\) for different spins S. b Squared average fidelity (over nuclear configurations) for different qudit spins S, for initial state \(\left|{\psi }_{L}\right\rangle =[\left|{0}_{L}\right\rangle +{{{\rm{i}}}}\left|{1}_{L}\right\rangle ]\sqrt{2}\) and numerical codewords \(\left|{0}_{L}\right\rangle\) and \(\left|{1}_{L}\right\rangle\) optimized at t = 5 μs. c Infidelity \(1-{{{{\mathcal{F}}}}}_{S}^{2}\) for an inset of panel (b) representing the region of \({{{{\mathcal{F}}}}}_{S}^{2}\ge 0.9\).

Full size image

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-{{{{\mathcal{F}}}}}_{S}^{2}(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 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,

$${{{{\mathcal{G}}}}}_{S}(t)=\frac{1-{{{{\mathcal{F}}}}}_{1/2}^{2}(t)}{1-{{{{\mathcal{F}}}}}_{S}^{2}(t)},$$
(9)

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.

Fig. 4: Gain, defined in Eq. (9), with respect to non-corrected spin 1/2 for the data in Fig. 3.
figure 4

The curves indicate the average of \({{{{\mathcal{G}}}}}_{S}(t)\) over the different nuclear spin spatial configurations (generated as explained in “Methods”), while the shaded areas mark the region included between \({{{{\mathcal{G}}}}}_{S}(t)\pm {\sigma }_{S}(t)\), where \({\sigma }_{S}(t)=\sqrt{\langle {{{{\mathcal{G}}}}}_{S}^{2}(t)\rangle -{\langle {{{{\mathcal{G}}}}}_{S}(t)\rangle }^{2}}\) 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”.

Full size image

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 topt = 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 (\theta )\left|{0}_{L}\right\rangle +{{{\rm{i}}}}\sin (\theta )\left|{1}_{L}\right\rangle\). 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.

Fig. 5: Optimization time and initial states.
figure 5

a Gain as a function of time for a logical state \(\left|{\psi }_{L}\right\rangle =(\left|{0}_{L}\right\rangle +{{{\rm{i}}}}\left|{1}_{L}\right\rangle )/\sqrt{2}\) for numerical codewords optimized at three different times topt = 10, 50, 100 μs and S = 3/2. b Fidelity (colour scale) as a function of time (radial scale) for different initial superposition states \(\cos \theta \left|{0}_{L}\right\rangle +{{{\rm{i}}}}\sin \theta \left|{1}_{L}\right\rangle\) (θ depicted in angular scale) of a set of optimized numerical codewords for S = 3/2.

Full size image

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}{\mu }_{{{{\rm{B}}}}}{B}_{y}(t)\cos (\omega t){\hat{S}}_{y}\) where the envelope By(t) is typically rectangular or Gaussian. These pulses induce transitions from a spin state \(\left|m\right\rangle\) to \(\left|m\pm 1\right\rangle\) when ω is set at the corresponding transition frequency, implementing a two-state unitary rotation \({Y}_{m}(\theta )=\exp [\theta /2(\left|m+1\right\rangle \left\langle m\right|-\left|m\right\rangle \left\langle m+1\right|)]\) between states \(\left|m\right\rangle\) and \(\left|m+1\right\rangle\) of an arbitrary angle θ. All the 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 codes13.

Fig. 6: Pulse sequence for S = 5/2.
figure 6

a Explicit sequence of pulses to implement the QEC for a spin S = 5/2. Each horizontal line represents a spin state \(\left|m\right\rangle\) with time flowing from left to right. Each grey box indicates a pulse implementing a rotation Ym(θ) as described in the text, with the values of the angles θk and ϕk given in Table 1. b Time evolution of the absolute value of the amplitude \(| \left\langle m| \psi (t)| \,\right\rangle |\) of the qudit state \(\left|\psi (t)\right\rangle\) over each \(\left|m\right\rangle\) state at different stages of the control sequence, corresponding to the control pulses depicted in (a). The sequence holds for any choice of α and β, though panel (b) shows an example with \(\alpha =1/\sqrt{2}\), \(\beta ={{{\rm{i}}}}/\sqrt{2}\). Blue and orange colours represent amplitudes associated to α and β (i.e., \(\left|{0}_{L}\right\rangle\) and \(\left|{1}_{L}\right\rangle\)), 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 \({\hat{E}}_{1}\) of Eq. (8)) is shown in (b).

Full size image

Encoding

We assume that the information, i.e., the coefficients α and β, is initially encoded in a simple state such as \(\alpha \left|-1/2\right\rangle +\beta \left|1/2\right\rangle\). The preparation of the logical state \(\alpha \left|{0}_{L}\right\rangle +\beta \left|{1}_{L}\right\rangle\) for arbitrary α and β is then realized by alternating pulses Ym(θ) distributing population among the different \(\left|m\right\rangle\) states and π-pulses Ym( ± π) 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 \(\left|m\right\rangle\) 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., \(\left|\pm 1/2\right\rangle\) states, performing the desired two-state operation, and re-encoding.

Table 1 Angles for the pulse sequence realizing the QEC for S = 5/2 as depicted in Fig. 6.
Full size table

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 \(\hat{H}\) for a time t/2 and with the same Hamiltonian but with \({\hat{S}}_{z}\) changed to \(-\hat{{S}_{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}_{{{{\rm{echo}}}}}={{{{\rm{e}}}}}^{-{{{\rm{i}}}}\pi {\hat{S}}_{x}}\) at time t/2 which inverts the spin, sending state \(\left|m\right\rangle\) to \(\left|-m\right\rangle\). 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 \(\left|m\right\rangle \leftrightarrow \left|-m\right\rangle\) states, followed by Ym(±π) pulses to rearrange populations. For instance, in the case of a spin 5/2, Uecho is obtained by three independent rotations \(\left|5/2\right\rangle \leftrightarrow \left|-5/2\right\rangle\), \(\left|3/2\right\rangle \leftrightarrow \left|-3/2\right\rangle\) and \(\left|1/2\right\rangle \leftrightarrow \left|-1/2\right\rangle\). Due to the lack of a direct matrix element between \(\left|m\right\rangle\) and \(\left|-m\right\rangle\) 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 sA = 1/2 ancilla. The ancilla is described by adding to the Hamiltonian of Eq. (12) the terms

$${{{\Omega }}}_{A}{\hat{s}}_{A}^{z}+\mathop{\sum}\limits_{k=x,y,z}{{\mathbb{J}}}_{k}{\hat{s}}_{A}^{k}{\hat{S}}^{k},$$
(10)

where \(\{{\hat{s}}_{A}^{x},{\hat{s}}_{A}^{y},{\hat{s}}_{A}^{z}\}\) 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 \({\mathbb{J}}\). For \({{\mathbb{J}}}_{x,y}\ll | {{\Omega }}-{{{\Omega }}}_{A}|\), such that states of qudit and ancilla remain essentially factorized, the excitation energies of the ancilla \({{{\Delta }}}_{A}^{(m)}\) depend on the state \(\left|m\right\rangle\) of the qudit via the diagonal coupling \({{\mathbb{J}}}_{z}{\hat{s}}_{A}^{z}{\hat{S}}_{z}\) only, i.e.,

$${{{\Delta }}}_{A}^{(m)}={g}_{A}{\mu }_{{{{\rm{B}}}}}{B}_{z}+{{\mathbb{J}}}_{z}m.$$
(11)

By irradiating the ancilla with a resonant magnetic pulse at angular frequency \({{{\Delta }}}_{A}^{(m)}\) it is thus possible to selectively excite the ancilla only if the qudit is in state \(\left|m\right\rangle\). A subsequent measurement of the state of the ancilla then reveals whether the qudit state has support on \(\left|m\right\rangle\) 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 recently63. This coupling can then be exploited to measure the ancilla, building on techniques well developed in the field of circuit quantum-electrodynamics64,65,66 and adapted to Molecular Nanomagnets67,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 \(\left|m\right\rangle\) states in both error spaces corresponding to \(\left|{0}_{L}\right\rangle\) and \(\left|{1}_{L}\right\rangle\) (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-to-implement measurement in the \(\left|m\right\rangle\) basis, and the unitarity of the transformation ensures that α and β are preserved. At this point, every possible post-error state is of the form \(\alpha \left|m\right\rangle +\beta \left|m^{\prime} \right\rangle\) for different pairs of states \((\left|m\right\rangle ,\left|m^{\prime} \right\rangle )\).

We now aim to induce an excitation of the ancilla only for a superposition state of the qudit with components on \((\left|m\right\rangle ,\left|m^{\prime} \right\rangle )\), without collapsing the superposition. We achieve this by applying a two-tone two-photon drive at frequencies \({{{\Delta }}}_{A}^{(m)}\), \({{{\Delta }}}_{A}^{(m^{\prime} )}\) (panel ‘ancilla measurement’ of Fig. 6)69. Then, in the dispersive limit, the coupling G between resonator and ancilla induces a shift of the cavity angular frequency ωc of ±G2/δm, with \({\delta }^{m}={{{\Delta }}}_{A}^{(m)}-{\omega }_{c}\) and the sign of the shift depending on the state of the ancilla64. Since here we need to measure the ancilla irrespective of the state of the qudit in the subspace \((\left|m\right\rangle ,\left|m^{\prime} \right\rangle )\), 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 in69. 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 \((\left|m\right\rangle ,\left|m^{\prime} \right\rangle )\) 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 \(\alpha \left|m\right\rangle +\beta \left|m^{\prime} \right\rangle\) with known m and \(m^{\prime}\). The simplest way to restore the encoded state \(\alpha \left|{0}_{L}\right\rangle +\beta \left|{1}_{L}\right\rangle\) is then to first use a few π-pulses to send \(\left|m\right\rangle \to \left|1/2\right\rangle\) and \(\left|m^{\prime} \right\rangle \to \left|-1/2\right\rangle\), 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 \(\left|m\right\rangle\), \(\left|m^{\prime} \right\rangle\) 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 techniques70,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.9865,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 codes13. 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 interactions46,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 works47 point out that the CCE method used here correctly reproduces the phenomenology of coherence enhancement due to the existence of a nuclear diffusion barrier53. 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 multiplets74,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^{\prime}\) 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.

Methods

Derivation of the effective Hamiltonian

The spin Hamiltonian describing the interacting evolution of the molecular spin \({{{\mathcal{S}}}}\) and the bath \({{{\mathcal{B}}}}\) of N nuclear spins is

$${\hat{H}}_{{{{\mathcal{S}}}}{{{\mathcal{B}}}}}={\hat{H}}_{{{{\mathcal{S}}}}}+\mathop{\sum }\limits_{n=1}^{N}{\omega }_{n}{\hat{I}}_{n}^{z}+\mathop{\sum }\limits_{n=1}^{N}\hat{{{{\boldsymbol{S}}}}}\cdot {{\mathbb{D}}}_{n}\cdot {\hat{{{{\boldsymbol{I}}}}}}_{n}+\mathop{\sum}\limits_{n\ne m}{\hat{{{{\boldsymbol{I}}}}}}_{n}\cdot {{\mathbb{E}}}_{n,m}\cdot {\hat{{{{\boldsymbol{I}}}}}}_{m},$$
(12)

where \({\hat{H}}_{{{{\mathcal{S}}}}}=D{\hat{S}}_{z}^{2}+{{\Omega }}{\hat{S}}_{z}\), \(\hat{{{{\boldsymbol{S}}}}}=\{{\hat{S}}^{z},{\hat{S}}^{+},{\hat{S}}^{-}\}\), \({\hat{{{{\boldsymbol{I}}}}}}_{n}=\{{\hat{I}}_{n}^{z},{\hat{I}}_{n}^{+},{\hat{I}}_{n}^{-}\}\). The tensors \({{\mathbb{D}}}_{n}\) contain dipole−dipole couplings between \({{{\mathcal{S}}}}\) and \({{{\mathcal{B}}}}\), while the tensors \({{\mathbb{E}}}_{n,m}\) contain nuclear−nuclear dipolar couplings. The elements of \({{\mathbb{D}}}_{n}\) satisfy

$${{\mathbb{D}}}_{n}^{++}={\left({{\mathbb{D}}}_{n}^{--}\right)}^{* },\quad {{\mathbb{D}}}_{n}^{+-}={\left({{\mathbb{D}}}_{n}^{+-}\right)}^{* }={{\mathbb{D}}}_{n}^{-+},{{\mathbb{D}}}_{n}^{+z}={\left({{\mathbb{D}}}_{n}^{-z}\right)}^{* },$$
(13)

and the same holds for the corresponding elements of \({{\mathbb{E}}}_{n,m}\). A canonical perturbative (Schrieffer-Wolff) transformation generated by

$$G=\mathop{\sum}\limits_{\beta \!=\!\pm\! ,z}[{G}_{k}^{+\beta }{\hat{S}}^{+}{\hat{I}}_{n}^{\beta }-\,{{\mbox{h.c.}}}\,],$$
(14)

with h.c. indicating the hermitian conjugate, is constructed such that \({\hat{H}}_{G}={{{{\rm{e}}}}}^{G}\hat{H}{{{{\rm{e}}}}}^{-G}\), within first order in Ω−1, does not contain off-diagonal couplings between \({{{\mathcal{S}}}}\) and \({{{\mathcal{B}}}}\) with respect to the states \(\left|m\right\rangle\)45,48,60,78. As detailed in the following, G is proportional to Ω−1 and leading orders in \({\hat{H}}_{G}\) can thus be computed from the Baker−Campbell−Hausdorf expansion

$${\hat{H}}_{G}=\hat{H}+[G,\hat{H}]+\frac{1}{2!}[G,[G,\hat{H}]]+\ldots \,.$$
(15)

The coefficients \({G}_{k}^{\alpha \beta }\) are determined explicitly by imposing the cancellation of the off-diagonal couplings between \({{{\mathcal{S}}}}\) and \({{{\mathcal{B}}}}\) to first order. This results in the relation

$$\left[G,{\hat{H}}_{{{{\mathcal{S}}}}}+\mathop{\sum }\limits_{n=1}^{N}{\omega }_{n}{\hat{I}}_{n}^{z}\right]=-\mathop{\sum}\limits_{\begin{array}{c}{\alpha\! =\!\pm\! ,}\\ {\beta\! =\!\pm\! ,z}\end{array}}\mathop{\sum }\limits_{n=1}^{N}{{\mathbb{D}}}_{n}^{\alpha \beta }{\hat{S}}^{\alpha }{\hat{I}}_{n}^{\beta }.$$
(16)

The coefficients \({G}_{k}^{\alpha \beta }\) then read

$$\begin{array}{l}{G}_{k}^{++}=\frac{{{\mathbb{D}}}_{k}^{++}}{{{\Omega }}\,+\,D\,+\,{\omega }_{k}},\quad {G}_{k}^{+-}=\frac{{{\mathbb{D}}}_{k}^{+-}}{{{\Omega }}\,+\,D\,-\,{\omega }_{k}},\\ {G}_{k}^{+z}=\frac{{{\mathbb{D}}}_{k}^{+z}}{{{\Omega }}\,+\,D}.\end{array}$$
(17)

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

$$\begin{array}{lll}{a}_{n}^{(0)}&=&{\omega }_{n},\qquad {b}_{n}^{(0)}=0,\qquad {c}_{n,m}^{(0)}={{\mathbb{E}}}_{n,m}^{+-},\\ {d}_{n,m}^{(0)}&=&{{\mathbb{E}}}_{n,m}^{zz}/2,\qquad {a}_{n}^{(1)}={{\mathbb{D}}}_{n}^{zz},\\ {b}_{n}^{(1)}&=&\frac{2}{{{\Omega }}}\left[| {{\mathbb{D}}}_{n}^{+z}{| }^{2}-| {{\mathbb{D}}}_{n}^{++}{| }^{2}-{\left({{\mathbb{D}}}_{n}^{+-}\right)}^{2}\right],\\ {c}_{n,m}^{(1)}&=&\frac{2}{{{\Omega }}}({{\mathbb{D}}}_{n}^{++}{{\mathbb{D}}}_{m}^{--}+{{\mathbb{D}}}_{n}^{-+}{{\mathbb{D}}}_{m}^{+-}),\\ {d}_{n,m}^{(1)}&=&\frac{2}{{{\Omega }}}{{\mathbb{D}}}_{n}^{+z}{{\mathbb{D}}}_{m}^{-z},\\ {a}_{n}^{(2)}&=&\frac{2}{{{\Omega }}}\left[| {{\mathbb{D}}}_{n}^{++}{| }^{2}-{({{\mathbb{D}}}_{n}^{+-})}^{2}\right],\end{array}$$
(18)

and \({b}_{n}^{(2)}={c}_{n}^{(2)}={d}_{n}^{(2)}=0\). The energy of \({{{\mathcal{S}}}}\) is also renormalized according to

$$\tilde{{{\Omega }}}={{\Omega }}+\frac{2I(I+1)}{{{\Omega }}}\mathop{\sum }\limits_{n=1}^{N}\left[| {{\mathbb{D}}}_{n}^{++}{| }^{2}+{({{\mathbb{D}}}_{n}^{+-})}^{2}\right],$$
(19)

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 \({{{\mathcal{S}}}}\), as sketched in Fig. 1. Further, a minimal distance of 3 Å is considered between nuclear spins and between each nuclear spin and \({{{\mathcal{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 Bz = 1 T along z is assumed, for achieving the regime of large Zeeman energy of \({{{\mathcal{S}}}}\). The decoherence function,

$${L}_{nm}(t)={{{\mbox{tr}}}}_{{{{\mathcal{B}}}}}\left[{{{{\rm{e}}}}}^{-{{{\rm{i}}}}{\hat{H}}_{{{{\mathcal{B}}}},n}t}{\rho }_{{{{\mathcal{B}}}}}(0){{{{\rm{e}}}}}^{{{{\rm{i}}}}{\hat{H}}_{{{{\mathcal{B}}}},m}t}\right]\,,$$
(20)

is computed by means of a CCE54,61. This expansion decomposes Lnm(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 Lnm(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 CCE-n, which yields the truncated function \({L}_{nm}^{(n)}(t)\).

The effect of inhomogeneous broadening is well captured by CCE-161,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

$${\rho }_{{{{\mathcal{B}}}}}(0)=\frac{{{{{\rm{e}}}}}^{-{H}_{{{{\mathcal{B}}}}}^{(0)}/{k}_{B}T}}{{{{\mbox{tr}}}}_{{{{\mathcal{B}}}}}[{{{{\rm{e}}}}}^{-{H}_{{{{\mathcal{B}}}}}^{(0)}/{k}_{B}T}]}\simeq {{\mathbb{1}}}_{{{{\mathcal{B}}}}}/{2}^{N}.$$
(21)

Under this approximation, the CCE-1 can be solved analytically also for S > 1/2 giving

$${L}_{nm}^{(1)}(t)=\mathop{\prod }\limits_{k=1}^{N}\cos \left[\frac{(n-m){{\mathbb{D}}}_{k}^{zz}}{2}t\right].$$
(22)

Here, \({{\mathbb{D}}}_{k}^{zz}\) is the hyperfine coupling \(\propto {\hat{S}}_{z}{\hat{I}}_{k}^{z}\) between \({{{\mathcal{S}}}}\) and the k-th nuclear spin. For small \({{\mathbb{D}}}_{k}^{zz}t\), this is well approximated by \({{{{\rm{e}}}}}^{-{(n-m)}^{2}{{{\Gamma }}}^{2}{t}^{2}}\) with \({{{\Gamma }}}^{2}={\sum }_{k}{({D}_{k}^{zz})}^{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}_{n,m}^{(k)}\), \({b}_{n,m}^{(k)}\), \({c}_{n,m}^{(k)}\), \({d}_{n,m}^{(k)}\) reveals three fundamental interaction scales at play, which ordered by decreasing strength, are associated to (i) the diagonal coupling between \({{{\mathcal{S}}}}\) and nuclear spins (terms \(\propto {\hat{S}}_{z}{\hat{I}}_{k}^{z}\) which are compensated for by the echo), (ii) the intrinsic interacting evolution (terms \(\propto {\hat{I}}_{n}^{\alpha }{\hat{I}}_{m}^{\beta }\)), (iii) the \({{{\mathcal{S}}}}\)-conditioned interacting evolution (terms \(\propto {\hat{S}}_{z}{\hat{I}}_{m}^{\alpha }{\hat{I}}_{n}^{\beta }\)). 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 \(\left|{0}_{L}\right\rangle\) and \(\left|{1}_{L}\right\rangle\) are found by starting from the following ansatz, inspired by spin-binomial codes13,

$$\left|{0}_{L}/{1}_{L}\right\rangle =\mathop{\sum }\limits_{\begin{array}{c}{\ell \!=\!0}\\ \ell \ {{{\rm{odd}}}}/{{{\rm{even}}}}\end{array}}^{S+1/2}{\gamma }_{\ell }^{(0)/(1)}\left|-S+\ell \right\rangle .$$
(23)

This ansatz permits one to impose Eq. (6) by construction and hence to reduce the number of free parameters, thanks to \(\left|{0}_{L}\right\rangle\) and \(\left|{1}_{L}\right\rangle\) having non-zero components on different sets of \(\left|m\right\rangle\) states.

Knill−Laflamme conditions (5) are finally enforced on the coefficients \({\gamma }_{\ell }^{(0)/(1)}\) by numerically minimizing the function

$$\mathop{\sum }\limits_{k,j = 0}^{2S}\left|\left\langle {0}_{L}\right|{\hat{E}}_{k}^{{\dagger} }{\hat{E}}_{j}\left|{0}_{L}\right\rangle -\left\langle {1}_{L}\right|{\hat{E}}_{k}^{{\dagger} }{\hat{E}}_{j}\left|{1}_{L}\right\rangle \right|.$$
(24)

Detection and recovery

The abstract detection and recovery operations follow the general QEC theory of ref. 62. Once a set of operators \({\{{\hat{E}}_{k}\}}_{k = 0,\ldots ,\lfloor S\rfloor }\) to be corrected with a spin S is identified, the two error subspaces corresponding to \(\left|{0}_{L}\right\rangle\) and \(\left|{1}_{L}\right\rangle\) are first determined. These two subspaces are defined as the linear span of the states \({\hat{E}}_{k}\left|{0}_{L}\right\rangle\) and \({\hat{E}}_{k}\left|{1}_{L}\right\rangle\) for all k, respectively. For each of these subspaces, a basis \(\{\left|{e}_{k}^{(0)}\right\rangle \}\) (\(\{\left|{e}_{k}^{(1)}\right\rangle \}\)) is selected. This is chosen here to be given by a Gram−Schmidt orthonormalization of states \({\hat{E}}_{k}\left|{0}_{L}\right\rangle\) (\({\hat{E}}_{k}\left|{1}_{L}\right\rangle\)). The detection measurement is then described by the projectors

$${\hat{P}}_{k}=\left|{e}_{k}^{(0)}\right\rangle \left\langle {e}_{k}^{(0)}\right|+\left|{e}_{k}^{(1)}\right\rangle \left\langle {e}_{k}^{(1)}\right|,$$
(25)

with k = 0, …, ⌊S⌋. Finally, the recovery operation, given an outcome j of the measurement, maps back the states \(\left|{e}_{j}^{(0)}\right\rangle\) and \(\left|{e}_{j}^{(1)}\right\rangle\) corresponding to that outcome to the computational states \(\left|{0}_{L}\right\rangle\) and \(\left|{1}_{L}\right\rangle\), respectively. Since the coefficients α and β of the encoded state have been preserved, this operation then fully restores the logical state \(\left|{\psi }_{L}\right\rangle\). The recovery is formalized by a set of transformations \(\{{\hat{R}}_{k}\}\) such that

$${\hat{R}}_{k}\left|{e}_{k}^{(c)}\right\rangle =\left|{c}_{L}\right\rangle ,$$
(26)

for c = 0, 1. Here, each transformation \({\hat{R}}_{k}\) is constructed as a rotation in the two-dimensional space spanned by \(\left|{e}_{k}^{(c)}\right\rangle\) and \(\left|{c}_{L}\right\rangle\).

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 Ym(θ) pulses, one can exploit known decomposition strategies from quantum control theory13,80. In a first step, the unitary U can be decomposed into a sequence of two-state planar rotations using the algorithm given in ref. 80. This gives a sequence of two-state transfers which involve states \(\left|m\right\rangle\) and \(\left|m^{\prime} \right\rangle\) also with \(| m-m^{\prime} | > 1\). To further convert such two-state rotations into rotations Ym(θ), i.e., with \(| m-m^{\prime} | =1\), one finally exploits π-pulses to bring the population of \(m^{\prime}\) close to m and back. For instance, defining

$${Y}_{m,m^{\prime} }(\theta )=\exp \left[\theta/2 \left(-\left|m\right\rangle \left\langle m^{\prime} \right|\right.+\left|m^{\prime} \right\rangle \left\langle m\right|\right]$$
(27)

for \(m^{\prime}\, > \,m\), one can iteratively exploit the properties

$${Y}_{m,m+2}(\theta )={Y}_{m,m+1}(\pi ){Y}_{m+1,m+2}(-\theta ){Y}_{m,m+1}(-\pi ).$$
(28)
$$={Y}_{m,m+1}(-\pi ){Y}_{m+1,m+2}(\theta ){Y}_{m,m+1}(\pi ).$$
(29)

As an example, the pulse sequence depicted in panel ‘basis rotation’ of Fig. 6a, which realizes a transformation mapping the basis of error words \(\left|{e}_{k}^{(c)}\right\rangle\) to the basis of \(\left|m\right\rangle\) 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 codes13. These codes are based on a description of spin dephasing as produced by a Markovian bath which couples to the system via operator \(\sqrt{\gamma }{\hat{S}}_{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 \({\hat{S}}_{z}^{n}\). We find that, while spin-binomial 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 \({\hat{S}}_{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.

Fig. 7: Performance of spin-binomial codes.
figure 7

a Average squared fidelity as a function of time for different values of the qudit spin S, for initial state \(\left|{\psi }_{L}\right\rangle =(\left|{0}_{L}\right\rangle +{{{\rm{i}}}}\left|{1}_{L}\right\rangle )/\sqrt{2}\) and spin-binomial codewords \(\left|{0}_{L}\right\rangle\) and \(\left|{1}_{L}\right\rangle\). The results are averaged over 26 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.

Full size image