Theoretical Design of Optimal Molecular Qudits for Quantum Error Correction

We pinpoint the key ingredients ruling decoherence in multispin clusters, and we engineer the system Hamiltonian to design optimal molecules embedding quantum error correction. These are antiferromagnetically coupled systems with competing exchange interactions, characterized by many low-energy states in which decoherence is dramatically suppressed and does not increase with the system size. This feature allows us to derive optimized code words, enhancing the power of the quantum error correction code by orders of magnitude. We demonstrate this by a complete simulation of the system dynamics, including the effect of decoherence driven by a nuclear spin bath and the full sequence of pulses to implement error correction and logical gates between protected states.

S = 1/2 (red) exchange multiplets as a function of J/J ′ . For J/J ′ ≳ 1.2 we find eight doublets below the first S = 3/2. The effect of such a DMI is also a re-normalization of the energy gaps between different multiplets, evident also in zero field. We stress that a different form of the DMI (including also x or y components) would modify the energies and the inter-multiplet matrix elements without altering our conclusions.
We finally note that addressing the energy gaps between the lowest multiplets in Fig. S1-S-2 (b) by EPR requires frequencies in the W-band range. This choice was motivated by the need to resolve all the different energy gaps. However, it is important to stress that proper pulse-shaping techniques could be employed to reduce leakage even in presence of closer energy levels, obtained for instance by an overall reduction of all the couplings in the spin Hamiltonian.

Derivation of the master equation
In this section, we derive the master equation used to describe decoherence given by the S-3 where N i = 7 is the number of ions and N s the number of nuclear spins in the bath. The operators s α j and I α n (with α = x, y, z) are spin operators for the ions and the nuclear spins, respectively. The dipolar coefficients d αβ jn are given explicitly in the main text. In the basis of the eigenstates |µ⟩ of the system Hamiltonian, defined by H S = µ E µ |µ⟩⟨µ|, the interaction In this expression, we have denoted with I β n (t) = e iH B t I β n e −iH B t the interaction picture coupling operators for the nuclear spins, we have used the notation L µµ ′ = |µ⟩⟨µ ′ |, and the coupling coefficients readd In the following, we consider only diagonal operators on the system (α = z in Eq. S1-S3), because transitions between different eigenstates of the system induced by the bath are forbidden by the large mismatch between system (E µ − E µ ′ with µ ̸ = µ ′ ) and bath energy gaps, compared to the couplings d αβ jn . Hence, terms ∝ L µµ ′ in Eq. (S2) can be neglected in the secular approximation yielding H SB (t) = Ns n=1 µ,µ ′ β=x,y,zd β n,µµ L µµ ⊗ I β n (t) .
Note thus that the structure of the system eigenstates will be contained ind β n,µµ ′ ≡d zβ n,µµ ′ through the matrix elements of the spin operators s z j . In interaction picture and under Born and Markov approximations, the reduced dynamics of the system can be described by the Markovian quantum master equation S1 S-4 Here, ρ(t) is the reduced density matrix of the system, obtained from the full system-bath density matrix ρ SB (t) by tracing over the bath degrees of freedom, ρ(t) = tr B [ρ SB (t)]. Inserting Eq. (S4) into (S5) and using the shortcut notation E µµ ′ = E µ − E µ ′ , one arrives at the expression where we have also introduced the bath noise spectrã as the (one-sided) Fourier transform of the bath correlation functions The correlation functions carry the information about the structure of the bath. In general, they are not accessible for interacting spin baths and can at best be approximated with numerical methods. S2 Introducing the rates equation (S6) can then be brought into the form S-5 where D(L)ρ is a Lindblad dissipator, S1 D(L)ρ = 2LρL † − L † Lρ − ρLL † . The term in the second line of Eq. (S12) is apparently not in Lindblad form, but it can be brought in such a form by diagonalizing the matrix Γ µν . S1 The Hamiltonian H LS = µ S µ |µ⟩⟨µ| is diagonal in the energy basis and only describes weak energy shifts of the system eigenstates. It will thus be neglected in the numerical simulations. The master equation (S12) will be used (with additional Hamiltonian terms describing the resonant pulses) in the numerical simulations.
While diagonal elements of ρ S (t) are preserved, off-diagonal elements in the energy basis as given in the main text. In the last line, Eq. (S14), we have introduced . Solving (S13) thus gives the exponential decay ρ µν (t) = e −γµν t ρ µν (0) for the off-diagonal elements, with decay rate γ µν = −2Γ µν + Γ µµ + Γ νν . In the following, since our focus in on how the structure of system eigenstates, rather than the structure of the bath, impacts decoherence, we will assume that the bath noise spectrum is constant at zero energy, ζ ββ ′ nn ′ (0) = ζ 0 . We have checked, by considering a random distribution of ζ ββ ′ nn ′ (0), that this choice does alter our conclusions.

Comparison between F and C cases
In order to compare the effect of ferromagnetic (F) vs competing (C) exchange interactions, we consider the molecular structure reported in Fig. 1-(a) of the main text with either S-6 negative (F) or positive (C) J i,j . For F, one can easily compute the elements of Γ mm ′ within the ground S = 9/2 multiplet: includes the effect of the coupling between system and bath spins. By assuming isotropic g j , one gets d αβ the hydrogen nuclear spin spectroscopic factor, µ N is the nuclear magneton and r α jn is the α component of the distance between the j-th electronic spin and the n-th nuclear spin. ζ 0 is the bath spectral function, assumed to be a constant.
Here we have exploited the Wigner-Eckart theorem within the S multiplet ⟨m|s z j |m⟩ = p j ⟨m|S z |m⟩. Then, neglecting small anisotropy effects, ⟨m|S z |m⟩ = m. We have introduced here a phenomonological coherence time (as usually done in experiments), defined as . This can be substituted for comparison in the analogous expression for the antiferromagnetic case, thus getting rid of the unknown bath spectral function ζ 0 . S-7 Impact of the structure of the eigenstates on decoherence In the following, we argue that the essential dependence on the system eigenstates that rules decoherence at short times, as outlined above, is confirmed also when avoiding to make Markov approximations for the bath. In particular, by restarting from the microscopic system-bath interaction Eq. (S1) and neglecting off-diagonal couplings (which is justified since the gaps in the system are much larger than the coupling to the bath), the interactionpicture system-bath Hamiltonian in the basis of system eigenstates |µ⟩ can be written in the and coefficients h µ n = j d zz jn ⟨µ| s z j |µ⟩ . We will indicate with U ± µ (t) the propagator generated by ±H µ (t). Although U ± µ (t) do have the form of a simple exponential of ±H µ (t) because of the time dependence, for short times (compared to the inverse of the system-bath couplings) it can be estimated with a Magnus expansion, where we have used Eq. (S22) and kept terms up to second-order in system-bath coupling.
The free evolution of the system coherences (in interaction picture for the system) can then be written like S-8 where ρ SB is the system-bath density matrix. Under Born approximation for the initial state, ρ SB (0) = ρ(0) ⊗ ρ B (0), this gives the coherence decay ρ µν (t) = L µν (t)ρ µν (0) with Assuming that the time t is small as compared to the energy scales in H gives, for µ ̸ = ν, Keeping only the leading term up to first order in t gives L µν (t) = n cos[(h ν n − h µ n )t/2] . This is well approximated at short times by the Gaussian decay L µν (t) ≈ e −Γt 2 , with Hence, even without the Markov approximation, the short-time decay rate of the system coherences is related to the spin structure of the systems by a squared sum of weighted differences between expectation values of local spin operators.

S-9
Process tomography and derivation of the code words We consider a subset of the system Hilbert space of dimension d. Starting from the solution of the Lindblad equation for the pure dephasing mechanism described in this work, we look for a set of operators {E k } (such that k E k E † k = I) providing a Kraus decomposition for the system density matrix at a given time t at which the error correction will be performed, i.e.
The operators E k represent the different quantum errors occurring on the system. From Eq.(S32) we note that the E k operators must be diagonal. Hence, only d linearly independent E k exist. These can in general be determined through quantum process tomography, as described in Ref. S3 (Sec. 8.4.2), which is simplified here given the diagonal structure of the E k . We first decompose the E k on a basis of diagonal operatorsẼ i , and we re-write Eq. (S33) as with χ ij = k ϵ ik ϵ * jk . By choosing the operators {Ẽ i } as projectors onto the eigenstates E µ = |µ⟩⟨µ| and comparing with Eq. (S32) one easily finds χ ij = e −γ ij t .
Bi diagonalizing χ, one then finds the Kraus operators. Indeed, being V † the matrix that diagonalizes χ and ξ the vector of the eigenvalues, we can rewrite S-10 from which one obtains Detection of a number n of different errors is possible provided the Hilbert space has at least dimension 2n: indeed, one needs at least 2n orthogonal states for mapping the effect on each error on the two logical states, E k |0 L ⟩ and E k |1 L ⟩, to distinguishable states. In the systems considered here, one cannot in general correct all d operators E k , but only a subset of size d/2. It is thus important to identify the error operators that give largest contribution to reconstruct ρ(t). To this end, we observe that the squared norm tr(E † k E k ) of E k is proportional to the eigenvalues ξ k of χ. Hence, by sorting the eigenvalues of the χ matrix in decreasing order, the required hierarchy for the Kraus operators E k is obtained.
Once the dominant error operators E k are determined, a quantum error correction scheme is formulated by finding logical states |0⟩ L and |1⟩ L fulfilling Knill-Laflamme conditions (KLc), Eq. 4 of the main text.
To this end, we consider code words defined on orthogonal superpositions of the the system eigenstates: where A and B are disjoint subsets of d/2 levels. With this choice, |0 L ⟩ and |1 L ⟩ are not mixed by diagonal E k operators by construction, thus automatically fulfilling the second Exploiting the decomposition of E k into projectors on the eigenstates (Eq. S37), the first KLc can be re-written as S-11 This set of conditions for j, k = 0, ..., d/2 − 1, together with the normalisation conditions µ |a µ | 2 = µ |b µ | 2 = 1, form a linear system of equations in the variables |a µ | 2 , |b µ | 2 . This system is in general overdetermined, but numerically underdetermined in practical cases.
Hence, at least one solution can be found, either by matrix pseudo-inversion or by numerical optimisation.

Actual implementation of Quantum Error Correction
To identify the different errors E k we generalise the procedure usually exploited in qubitbased quantum error correction codes to measure a set of stabilizers. S5

Error detection
The scheme we follow is summarised in Fig. 3(a) of the main text. Here we detail each step and follow the corresponding evolution of system state. To do this, we first introduce an orthonormal basis set of so-called error words {|e 0 k ⟩ , |e 1 k ⟩}, k = 0, ...d/2 − 1 spanning the error spaces defined as span{E k |0 L ⟩} k=1,...,d/2 and span{E k |1 L ⟩} k=1,...,d/2 , respectively. For performing the error detection, we need to realize a measurement distinguishing different error words, described by projectors P k = |e 0 k ⟩⟨e 0 k | + |e 1 k ⟩⟨e 1 k |. In this way, the system is projected post-measurement into a known error word in both error spaces, without collapsing the stored information. This then allows one to implement a corresponding k-dependent recovery operation. Implementing these projectors onto superposition states of the system eigenstates directly can be hard experimentally. Therefore, in our scheme we first perform a transformation that entangles each of the d/2 error words α |e 0 k ⟩ + β |e 1 k ⟩ with a different eigenstate of an ancillary d/2-level qudit. The desired projectors on the system are then obtained by measuring the ancilla in its eigenbasis.
The sequence of operations (see Fig. 3 of the main text) is the following: 0. Start with the qudit in a generic logical state |ψ⟩ (on which pure dephasing has possibly S-12 occurred) and a d/2−levels ancillary qudit in its ground state |0⟩, i.e. the combined ancilla-qudit state is the tensor product |Φ⟩ = |0⟩ ⊗ |ψ⟩.