Introduction

The paradigmatic bases of information in Quantum Information Processing (QIP) are qubits: two-level individually addressable quantum systems. However, several QIP platforms have recently been proposed that instead make use of d-level systems, referred to as qudits1,2,3,4,5,6. In its infancy, classical computing did experiment with ternary, quaternary, or higher-dimensional, bases of information, before eventually settling on the simplest (bits), when near-zero error rates and easy scalability were attained7. Analogously, it could be argued that quantum computing is likely to follow a similar trend in the long-term; as fault-tolerant platforms emerge and technologies mature, the industry could indeed fully settle on multiqubit systems. However, QIP research is currently not in the noise-free regime but near it, and in order to reach significant quantum supremacy8, increasing the total Hilbert space dimension of the physical platform is a primordial requirement. As such, there is a current race to increase the number n of coupled qubits (d = 2) with superconducting platforms leading the way with n = 519 or n = 43310. While in general the Hilbert space dimension increases exponentially in the number of sites, the relatively slow 2n scaling of qubits, compared to dn for qudits, is proving challenging, necessitating ever-more robust systems and complex control mechanisms.

Given these current technical challenges, most qudit-based platforms that have been physically implemented argue for near-term advantages over equivalent multiqubit implementations. Thus, the principal motivations of qudit platforms over qubits include: (i) the underlying physical systems having lower decoherence rates11, (ii) using the redundancy in additional levels for quantum error correction12,13, (iii) the higher density of information per physical system (site)14, (iv) the reduced number of nonlocal, hence more decoherence-sensitive, operations15 or (v) more robust flying quantum memories16,17. Furthermore, qudits present fundamental theoretical advantages, enabling QIP capabilities offered by ⨂SU(d) vs. ⨂SU(2) of qubits18 such as simplifying some quantum algorithms19, and therefore a fault-tolerant qudit quantum computer indeed remains conceivable. Hence, qudits provide an alternative scaling solution by linearly increasing d, instead of scaling up n, the number of sites, as well as increasing efficiency through single qudit gates operating on larger computational subspaces4,9,20,21. However, one of the disadvantages raised for qudits is the larger number of error channels compared to multiple qubits22. In this context, a study of the near-term viability of qudits is needed to investigate the interplay between computational efficiency and noise error rates in higher dimensions.

In this work, we consider one single qudit versus multiqubit systems, in the context of near noise-free implementations. We undertake an inquiry to determine under what conditions on the applied gates a single qudit system does not lose more computational information than an equivalent multiqubit system, even when the qudit system initially presents more potential error channels. For this purpose, a standard measure to quantify the loss of computational information, that we study, is the Average Gate Infidelities (AGI), as defined by Nielsen23 where the average is over the Haar measure. The choice of the AGI ensures the calculated fidelity is not dependent on the input state and therefore remains relevant even if the gate is applied in later stages of a quantum algorithm. We conduct an in-depth analysis in which we compare the computational fidelity of a single qudit and n-qubit systems, both with identical dimensions of Hilbert space, undergoing arbitrary unitary transformations, and evolving under the influence of comparable noisy conditions. Our benchmark for successful analysis is defined by a lower first-order response of the AGI to the environmental noise, providing a measure of computational fidelity independent of initial states. We investigate, for increasing values of d, the respective growth rates of the AGIs with respect to error rates γ and dimensionless gate time γt. The latter quantifies the gate efficiency by indicating how time-efficient operations on these systems are relative to decoherence timescales, therefore this paper presents a study of the first-order connection between the AGI and this time-efficiency γt.

In other words, this study aims to investigate how the AGI scales proportionally with both the error rate and the speed at which gate operations are performed, as well as the dimension d of the qudit. Additionally, this study aims to provide a benchmarking tool to decide if a qudit platform, for a given (γ,t,d) specification, can compensate for its greater number of error channels by leveraging advantageous decoherence times and gate speeds. Both of those quantities depend intrinsically on the physical platform implementing the single qudit or multiqubit system, in particular their coupling to the environment, the mathematical form of the control pulses’ Hamiltonian, and the addressing speed. In particular, given a single qudit platform and a multiqubit platform with equivalent Hilbert space dimensionality, and specifying a fixed pair of parameters (γ,t), one could conduct a comparative analysis to determine if the qudit platform exhibits sufficiently low decoherence and sufficiently rapid gate time to achieve computational fidelities that are competitive with the multiqubit platform. Or, similarly, since increasing d on a single site in a given qudit platform is a prevailing goal for some platforms14, assuming γt remains of the same order of magnitude, this study would also allow setting theoretical upper limits on the value of d in order to remain advantageous.

In the first part, a gate-independent formula is presented for the first-order response in γt of the AGI to Markovian noise in the Lindblad formalism. The first-order formalism corresponds to the quasi-errorless regime of near-term QIP systems. Expressions for the linear dependency of the AGI on γt for a single qudit, multiqubits and also multiqudits are derived for an arbitrary collapse operator. A comparison is then made between the rate of increase of the AGI of a single qudit vs. equivalent multiple qubits.

This is then followed by numerical simulations, performed with the Python package QuTiP24, that complement and illustrate the analytical results. Discussions of the applicability and limits of the linear response formalism for AGI are given and the following aspects are studied: (i) the applicable range of γt and its dependency on the dimension of the qudit; (ii) the extent of the gate-independence of the result; (iii) the applicability to noise models other than pure dephasing; and finally (iv) the conditions on gate times for which either qudits or multiple qubits are advantageous. This latter aspect is then examined in more detail with respect to existing platforms by taking into account their respective decoherence rate and gate operation time.

Results and Discussion

Fluctuation-dissipation relation for a perturbed pure state

Consider a qudit, a d-level quantum system whose dynamics are governed by the Lindblad master equation25:

$$\frac{\,{{\mbox{d}}}\,{{{\rm{\rho }}}}}{\,{{\mbox{d}}}\,t}=-{{{\rm{i}}}}\left[H,\rho \right]+\mathop{\sum }\limits_{k=1}^{K}{\gamma }_{k}\left({L}_{k}\rho {L}_{k}^{{\dagger} }-\frac{1}{2}\left\{{L}_{k}^{{\dagger} }{L}_{k},\rho \right\}\right),$$
(1)

where ρ(t) is the density matrix of the system at time t, H the Hamiltonian of the system, Lk the so-called collapse operators characterizing the Markovian noise, and γk the decay parameters for each of the K noise processes. H = H0 + Hc(t) where H0 models the free evolution of the physical system and encompassing its internal interactions, and Hc(t) is a time-dependent pulse Hamiltonian allowing the controlled evolution. Moreover, the interactions of H with the collapse operators determine relevant timescales such as the gate-time t and the decoherence time T2 that are thus inherent to the physical realization under consideration.

The aim is to study the effect of a single collapse operator \(\sqrt{{\gamma }_{1}}{L}_{1}=\sqrt{\gamma }L\) on short timescales and under small-amplitude noise, i.e., γt ≪ 1. Under these assumptions, one can consider an ansatz of the form:

$$\rho (t)={\rho }^{{{{\rm{* }}}}}-\gamma tM+{{{\mathcal{O}}}}({(\gamma t)}^{2}),$$
(2)

with ρ* the noiseless target state, which is the solution of \(\dot{\rho }=-{{{\rm{i}}}}\left[H,\rho \right]\) after time t, and M the perturbation matrix resulting from the presence of a small-amplitude noise. Terms in \({{{\mathcal{O}}}}({(\gamma t)}^{2})\) include terms whose prefactor is of the form \({\left({\gamma }^{l}{t}^{k}\right)}_{l+k\ge 3}\), a more in-depth discussion is available in Supplementary Notes 2b.

One can easily see that the use of (2) in (1) leads to

$$M=\frac{1}{2}\left\{{L}^{{\dagger} }L,{\rho }^{{{{\rm{* }}}}}\right\}-L{\rho }^{{{{\rm{* }}}}}{L}^{{\dagger} }.$$
(3)

Consider now a quantum operation bringing the initial state ρ0 to a final state ρ(t) at time t. One define the fidelity \({{{\mathcal{F}}}}\) of this final state relative to some target state ρ*26 as

$${{{\mathcal{F}}}}(\rho (t),{\rho }^{{{{\rm{* }}}}})\equiv {\left[{{{\rm{Tr}}}}\left(\sqrt{\sqrt{\rho (t)}{\rho }^{{{{\rm{* }}}}}\sqrt{\rho (t)}}\right)\right]}^{2}.$$
(4)

Subsequently, the infidelity, is then defined as

$${{{\mathscr{E}}}}\equiv 1-{{{\mathcal{F}}}}.$$
(5)

Since ρ* is a pure state (\({\rho }^{{{{\rm{* }}}}}=\left\vert {\varphi }^{{{{\rm{* }}}}}\right\rangle \left\langle {\varphi }^{{{{\rm{* }}}}}\right\vert\)) Eq. (4) simplifies to26

$${{{\mathcal{F}}}}(\rho (t),{\rho }^{{{{\rm{* }}}}})={{{\rm{Tr}}}}\left(\rho (t){\rho }^{{{{\rm{* }}}}}\right).$$
(6)

Finally, substituting Eq. (2) into Eq. (6) leads to (see Supplementary Notes 1a)

$${{{\mathscr{E}}}}({\rho }^{{{{\rm{* }}}}})=\gamma t{\Delta }_{{{{\rm{* }}}}}L+{{{\mathcal{O}}}}({(\gamma t)}^{2}).$$
(7)

where \({\Delta }_{{{{\rm{* }}}}}L={\langle {L}^{{\dagger} }L\rangle }_{{{{\rm{* }}}}}-{\langle {L}^{{\dagger} }\rangle }_{{{{\rm{* }}}}}{\langle L\rangle }_{{{{\rm{* }}}}}\) with \({\langle {L}^{{\dagger} }L\rangle }_{{{{\rm{* }}}}}\equiv {{{\rm{Tr}}}}\left({\rho }^{{{{\rm{* }}}}}{L}^{{\dagger} }L\right)\).

Average Gate Fidelity of a single qudit

Only \({{{\mathscr{E}}}}({\rho }^{{{{\rm{* }}}}})\) for a specific ρ* was obtained in the previous subsection. However, is there a state-independent approach to obtaining the infidelity of a quantum gate under small-amplitude noise? One defines the quantum gate U applied during a time duration t, whose resulting operation brings all initial states ρ0 to all corresponding ρ* = Uρ0U†. There is then a definition of the average gate fidelity of a quantum channel \({{{\mathcal{E}}}}\), attempting to carry the unitary operation U despite a noisy environment, which reads as follows23

$$\begin{array}{l}\bar{{{{\mathcal{F}}}}}({{{\mathcal{E}}}},U)\;=\;{\displaystyle\int}\,d{\rho }_{0}\,{{{\mathcal{F}}}}(\rho (t),{\rho }^{{{{\rm{* }}}}})\\ \qquad\qquad=\,{\displaystyle\int}\,d{\rho }_{0}\,{\left\langle {U}^{{\dagger} }{{{\mathcal{E}}}}[{\rho }_{0}]U\right\rangle }_{0}={\displaystyle\int}\,d{\rho }_{0}\,{\left\langle \left({{{{\mathcal{U}}}}}^{{\dagger} }\,\circ\, \,{{{\mathcal{E}}}}\right)[{\rho }_{0}]\right\rangle }_{0},\end{array}$$
(8)

where the normalized integral is over the Fubini-Study measure on pure states (sometimes called the Haar measure)27, \({{{{\mathcal{U}}}}}^{{\dagger} }[\rho ]\equiv {U}^{{\dagger} }\rho U\) and \({{{\mathcal{E}}}}[{\rho }_{0}]=\rho (t)\).

Introducting \({\tilde{E}}_{k}={E}_{k}U\) the Kraus operators such that

$$\rho (t)=\left({{{{\mathcal{U}}}}}^{{\dagger} }\,\circ\, {{{\mathcal{E}}}}\right)[{\rho }_{0}]=\mathop{\sum}\limits_{k}{\tilde{E}}_{k}{\rho }_{0}{\tilde{E}}_{k}^{{\dagger} }=\mathop{\sum}\limits_{k}{E}_{k}{\rho }^{{{{\rm{* }}}}}{E}_{k}^{{\dagger} },$$
(9)

the Average Gate Fidelity \(\bar{{{{\mathcal{F}}}}}\) given in (8) can be rewritten as28

$$\bar{{{{\mathcal{F}}}}}({{{\mathcal{E}}}},U)=\frac{d+{\sum }_{k}{\left\vert {{{\rm{Tr}}}}\left({\tilde{E}}_{k}{U}^{{\dagger} }\right)\right\vert }^{2}}{d(d+1)}=\frac{d+{\sum }_{k}{\left\vert {{{\rm{Tr}}}}\left({E}_{k}\right)\right\vert }^{2}}{d(d+1)}\,.$$
(10)

Using Eq. (3), one seeks the sets of Kraus operators \(\left\{{\tilde{E}}_{k}\right\}\) or \(\left\{{E}_{k}\right\}\) such that, to \({{{\mathcal{O}}}}({(\gamma t)}^{2})\), ∀ ρ0, ρ*, 

$$\begin{array}{l}\mathop{\sum}\limits_{k}{\tilde{E}}_{k}{\rho }_{0}{\tilde{E}}_{k}^{{\dagger} }\,=\,U{\rho }_{0}{U}^{{\dagger} }-\gamma t\frac{1}{2}\left\{{L}^{{\dagger} }L,U{\rho }_{0}{U}^{{\dagger} }\right\}+\gamma tLU{\rho }_{0}{U}^{{\dagger} }{L}^{{\dagger} }\\ \mathop{\sum}\limits_{k}{E}_{k}{\rho }^{{{{\rm{* }}}}}{E}_{k}^{{\dagger} }\,=\,{\rho }^{{{{\rm{* }}}}}-\gamma t\frac{1}{2}\left\{{L}^{{\dagger} }L,{\rho }^{{{{\rm{* }}}}}\right\}+\gamma tL{\rho }^{{{{\rm{* }}}}}{L}^{{\dagger} }.\end{array}$$
(11)

One can see that the following two sets would work up to the first order in γt

$$\begin{array}{ll}{\tilde{E}}_{0}=\left({{\mathbb{1}}}_{d}-\frac{\gamma t}{2}{L}^{{\dagger} }L\right)U,&\tilde{{E}_{1}}=\sqrt{\gamma t}LU,\\ {E}_{0}={{\mathbb{1}}}_{d}-\frac{\gamma t}{2}{L}^{{\dagger} }L,&{E}_{1}=\sqrt{\gamma t}L.\end{array}$$
(12)

In order to use Eq. (12) in Eq. (8), it is necessary to calculate the traces of the operators. Let us consider a pure dephasing channel of a qudit coupled to a thermal environment through the operator Jz (In general the coupling of a qudit, or qubit, to a thermal environment can be represented by a linear combination, or mixture, of collapse operators, though a pure dephasing channel will typically be present and can be represented by the operator Jz. As a toy model, let us consider a coupling term dominated by a pure dephasing channel.) i.e. \({{{{\mathcal{E}}}}}_{z}\) with L = Jz29. One obtains a gate- (and Hamiltonian-) independent result for the Average Gate Fidelity which reads as (see Supplementary Notes 1b)

$$\overline{{{{\mathcal{F}}}}}({{{{\mathcal{E}}}}}_{z})=1-\frac{\gamma t}{12}d(d-1)+{{{\mathcal{O}}}}({(\gamma t)}^{2}).$$
(13)

In other words the Average Gate Infidelity (AGI) is given by

$$\overline{{{{\mathscr{E}}}}}({{{{\mathcal{E}}}}}_{z})=\frac{\gamma t}{12}d(d-1)+{{{\mathcal{O}}}}({(\gamma t)}^{2}),$$
(14)

or more generally for an arbitrary quantum channel \({{{\mathcal{X}}}}\) with collapse operator L

$$\overline{{{{\mathscr{E}}}}}({{{\mathcal{X}}}})=\frac{\gamma t}{d+1}\left({{{\rm{Tr}}}}({L}^{{\dagger} }L)-\frac{1}{d}{\left\vert {{{\rm{Tr}}}}(L)\right\vert }^{2}\right)+{{{\mathcal{O}}}}({(\gamma t)}^{2}).$$
(15)

Note that it is always possible to find a traceless collapse operator L emulating \({{{\mathcal{X}}}}\)25, so the previous expression can be, in this case, simplified as follows

$$\overline{{{{\mathscr{E}}}}}({{{\mathcal{X}}}})=\frac{\gamma t}{d+1}{{{\rm{Tr}}}}({L}^{{\dagger} }L)+{{{\mathcal{O}}}}({(\gamma t)}^{2}).$$
(16)

It follows from (16) that, if L were independent of d, increasing the dimension d of the Hilbert space would also increase the robustness of qudit gates to a dimension-independent quantum channel.

AGI of qudits vs. qubits

Now let us apply the same technique as described above to another system: an ensemble of n identical dephasing qubits (Hilbert space of dimension d = 2n). In order to compare it with the qudit analysis in the previous subsection, each individual qubit decoheres with the same rate (has the same type, and strength, of environmental coupling) through its spin operator Sz in the same way as the individual qudit (with d = 2). Considering any additional coupling mechanism to the environment arising from inter-qubit interactions would only further disadvantage the multi-qubit implementation. Our considerations then provide a best-case scenario for comparable qubits. This yields the master equation

$$\frac{\,{{\mbox{d}}}\rho }{{{\mbox{d}}}\,t}=-{{{\rm{i}}}}\left[H,\rho \right]+\mathop{\sum }\limits_{k=1}^{n}{L}_{k}\rho {L}_{k}^{{\dagger} }-\frac{1}{2}\mathop{\sum }\limits_{k=1}^{n}\left\{{L}_{k}^{{\dagger} }{L}_{k},\rho \right\}\,,$$
(17)

with

$$L_k = \underbrace{{\mathbb{1}}_2^{(1)} \otimes \cdots \otimes {\mathbb{1}}_2^{(k-1)}}_{k-1}\otimes S_z^{(k)} \otimes \underbrace{{\mathbb{1}}^{(k+1)}_2 \otimes \cdots \otimes {\mathbb{1}}^{(n)}_2}_{n-k}$$
(18)

for k ∈ ⟦1, n⟧.

Using the same reasoning as for dephasing qudits one obtains n + 1 Kraus operators to first order in γt

$${E}_{0}={{\mathbb{1}}}_{{2}^{n}}-\mathop{\sum }\limits_{k=1}^{n}\frac{\gamma t}{2}{L}_{k}^{{\dagger} }{L}_{k},\quad {E}_{k}=\sqrt{\gamma t}{L}_{k}.$$
(19)

In this case (see Supplementary Notes 1c),

$$\overline{{{{\mathscr{E}}}}}({{{{\mathcal{E}}}}}_{z})=\frac{\gamma t}{4}\frac{n{2}^{n}}{{2}^{n}+1}+{{{\mathcal{O}}}}({(\gamma t)}^{2})=\frac{\gamma t}{4}\frac{{\log }_{2}(d)d}{d+1}+{{{\mathcal{O}}}}({(\gamma t)}^{2})\,.$$
(20)

Let us stress that (20) yields the same result as Abad et al.30 in the case of identically dephasing qubits with no energy relaxation.

The analytical expressions (16) and (20) are one of the main results of this work.

Following those last two results, two expressions for the AGI have been found: for a single qudit, one finds an infidelity that scales as d2: \(\overline{{{{{\mathscr{E}}}}}_{d}}({{{{\mathcal{E}}}}}_{z})={c}_{d}\gamma t\) in (14) and for an ensemble of n qubits one finds an infidelity that scales as \({\log }_{2}(d)\): \(\overline{{{{{\mathscr{E}}}}}_{{{{\rm{b}}}},n}}({{{{\mathcal{E}}}}}_{z})={c}_{{{{\rm{b}}}},n}\gamma t\) in (20). Moreover, in the case of pure dephasing, one can define the T2,d dephasing time between two energy-adjacent levels for a qudit. It then shares the same expression (in terms of γ) as the typical T2,b dephasing time of a single qubit, namely \(\frac{1}{{T}_{2}}=\frac{\gamma }{2}\).

The ratio between two average gate infidelities, of duration td and tb,n for qudit and n qubits, respectively, becomes

$$\frac{\overline{{{{{\mathscr{E}}}}}_{d}}({{{{\mathcal{E}}}}}_{z})}{\overline{{{{{\mathscr{E}}}}}_{{{{\rm{b}}}},n}}({{{{\mathcal{E}}}}}_{z})}=\frac{{c}_{d}{t}_{d}/{T}_{2,d}}{{c}_{{{{\rm{b}}}},n}{t}_{{{{\rm{b}}}},n}/{T}_{2,{{{\rm{b}}}}}}.$$
(21)

Therefore, in order for a single qudit (d = 2n) to outperform an ensemble of n qubits in noise-robustness, i.e., to have a smaller AGI, the following inequality must hold true

$$\frac{{t}_{{{{\rm{b}}}},n}/{T}_{2,{{{\rm{b}}}}}}{{t}_{d}/{T}_{2,d}} > \frac{{c}_{d}}{{c}_{{{{\rm{b}}}},n}}=\frac{{d}^{2}-1}{3{\log }_{2}(d)}=\frac{{4}^{n}-1}{3n}.$$
(22)

This expression quantifies the requirements on the figure of merit which is the gate time in units of decoherence time τd = td/T2,d relative to τb,n = tb,n/T2,b in order for the qudit to yield higher-fidelity gates. Moreover, it confirms that the infidelity of an ensemble of n identical qubits and the infidelity of a single qudit will generally not have the same linear behaviour in γt even if they have the same T2, thus simply having τd < τb,n is not sufficient to guarantee a more noise-resilient qudit. Moreover, Eq. (22) provides a more precise condition on the ratio of figure of merits than a simple qualitative result such as \(\frac{{d}^{2}}{{\log }_{2}(d)}\), while maintaining the expected \(O({d}^{2}/{\log }_{2}(d))\) behaviour as d → ∞. In particular, see Supplementary Notes 1b for the full analytical calculations, including the derivation of the non-trivial factor \(\frac{1}{3}\).

On a side note, the previous calculations could also be applied to an ensemble of N qudits under identical pure dephasing, in which case we have

$$\overline{{{{{\mathscr{E}}}}}_{d,N}}({{{{\mathcal{E}}}}}_{z})=\frac{\gamma t}{12}\frac{N{d}^{N}}{{d}^{N}+1}({d}^{2}-1)+{{{\mathcal{O}}}}({(\gamma t)}^{2}),$$
(23)

and for 2n = dN one obtains

$$\frac{{t}_{{{{\rm{b}}}},n}/{T}_{2,{{{\rm{b}}}}}}{{t}_{d,N}/{T}_{2,d}} > \frac{{c}_{d,N}}{{c}_{{{{\rm{b}}}},n}}=\frac{{d}^{2}-1}{3{\log }_{2}(d)}.$$
(24)

Let us also note that for L arbitrary (23) yields

$$\overline{{{{{\mathscr{E}}}}}_{d,N}}({{{\mathcal{X}}}})=\gamma tN\frac{{d}^{N-1}}{{d}^{N}+1}\left({{{\rm{Tr}}}}({L}^{{\dagger} }L)-\frac{1}{d}{\left\vert {{{\rm{Tr}}}}(L)\right\vert }^{2}\right)+{{{\mathcal{O}}}}({(\gamma t)}^{2}).$$
(25)

This equation encompasses both scenarios under investigation in this paper until this point. We recall that we construct two systems of equivalent Hilbert Space dimension: a single qudit of dimension d and a system of N qubits. (25) reflects two different scalings of the AGI in N and d respectively. In N, it is linear, and, in d, it scales as \(\left({{{\rm{Tr}}}}({L}^{{\dagger} }L)-\frac{1}{d}{\left\vert {{{\rm{Tr}}}}(L)\right\vert }^{2}\right)\) since the dimension affects the definition of L. In Supplementary Notes 1b we computed that spin-based L lead to a quadratic scaling in d. Therefore, in the qudit subsection, for a single qudit, we study this quantity for fixed N = 1 and varying d, while in the multiqubits subsection, it is for fixed d = 2, but varying N. The subtlety in the latter case being that, by construction \(N:= {\log }_{2}(d)\), with d being that of the single qudit, hence the different scalings in d in (14) and (20). See Fig. 1 for a visual summary of the results.

Fig. 1: Visual synthesis of the infidelity scaling in multiple qubits and a single qudit.
figure 1

Summary diagram illustrating the selected collapse operators and the associated analytically derived expected infidelity scalings as functions of the Hilbert Space dimension, as derived from (20) and (14). This is depicted for two distinct systems: multiple qubits (left) and a single qudit (right). The term 'infidelity scaling' here refers to the slopes of the first-order-in-γt AGIs, denoted as c in (21).

Full size image

And if each qudit k has a different set of noise parameters (γk, Lk), an even more general formula arises :

$$\overline{{{{{\mathscr{E}}}}}_{d,N}}({{{\mathcal{X}}}})=\frac{{d}^{N-1}}{{d}^{N}+1}\mathop{\sum }\limits_{k=1}^{N}{\gamma }_{k}t\left({{{\rm{Tr}}}}({L}_{k}^{{\dagger} }{L}_{k})-\frac{1}{d}{\left\vert {{{\rm{Tr}}}}({L}_{k})\right\vert }^{2}\right)+{{{\mathcal{O}}}}({({\gamma }_{k}t)}^{2}).$$
(26)

This formula through its general form, can be applied to any qudits whose physical implementation implies different collapse operators from the ones considered in this paper.

Process fidelity & averaged fluctuation-dissipation relation

One may link the fluctuation-dissipation relation obtained in (7) with the results regarding average gate infidelities from Eq. (15),

$$\overline{{{{\mathscr{E}}}}}({{{\mathcal{X}}}})=\gamma t\int\,d{\rho }^{{{{\rm{* }}}}}{{{\mathscr{E}}}}({\rho }^{{{{\rm{* }}}}})=\gamma t\int\,d{\rho }^{{{{\rm{* }}}}}{\Delta }_{{{{\rm{* }}}}}L+{{{\mathcal{O}}}}({(\gamma t)}^{2}).$$
(27)

This integral over the Fubini-Study measure can formally be computed using Weingarten calculus methods31 (see Supplementary Notes 1d) which can be expressed as

$$\int\,d\rho \Delta L=\frac{1}{d+1}{{{\rm{Tr}}}}({L}^{{\dagger} }L)-\frac{1}{d(d+1)}{\left\vert {{{\rm{Tr}}}}(L)\right\vert }^{2},$$
(28)

leading to (15).

In contrast to this formal approach, a more physically-informed approach to obtain the same result was proposed in the previous subsections.

Furthermore, it is possible to express all the computed average gate infidelities as process/entanglement infidelities \({{{{\mathscr{E}}}}}^{{{{\rm{(p)}}}}}\) making use of the relation \(D{{{{\mathscr{E}}}}}^{{{{\rm{(p)}}}}}=(D+1)\overline{{{{\mathscr{E}}}}}\), with D = d, 2n or dN, the dimension of the Hilbert space32. This yields the expression

$${{{{\mathscr{E}}}}}_{d,n}^{{{{\rm{(p)}}}}}({{{\mathcal{X}}}})=\gamma t\frac{n}{d}\left({{{\rm{Tr}}}}({L}^{{\dagger} }L)-\frac{1}{d}{\left\vert {{{\rm{Tr}}}}(L)\right\vert }^{2}\right)+{{{\mathcal{O}}}}({(\gamma t)}^{2}),$$
(29)

which is linear in the number of subsystems n. Likewise we have

$${{{{\mathscr{E}}}}}_{d}^{{{{\rm{(p)}}}}}({{{{\mathcal{E}}}}}_{z})=\frac{\gamma t}{12}({d}^{2}-1)+{{{\mathcal{O}}}}({(\gamma t)}^{2}),$$
(30)
$${{{{\mathscr{E}}}}}_{{{{\rm{b}}}},n}^{{{{\rm{(p)}}}}}({{{{\mathcal{E}}}}}_{z})=\frac{\gamma t}{4}n+{{{\mathcal{O}}}}({(\gamma t)}^{2}).$$
(31)

Note that (31) has been verified experimentally, for example by Ozaeta and McMahon33.

Fit and deviation from the linear behaviour

Using the procedures described in the Methods section, we simulated single qudits of dimension d under pure dephasing, with \(H={{\mathbb{0}}}_{d}\) and small γt ∈ [0, 10−4] (γ ~ 10−4 in some nuclear spins in molecular magnets such as in the experiments of Godfrin et al.20). For example, the simulations were performed for even dimensions d ∈ ⟦2, 22⟧. Fitting the AGIs \(\overline{{{{{\mathscr{E}}}}}_{d}}({{{{\mathcal{E}}}}}_{z})={c}_{d}\gamma t\) as a function of γt yielded the slopes cd that are shown in Fig. 2, along with their analytical expression as a function of d predicted in (14).

Fig. 2: Rate of increase of \(\overline{{{{{\mathscr{E}}}}}_{d}}({{{{\mathcal{E}}}}}_{z})={c}_{d}({J}_{z})\gamma t\) as a function of qudit dimension.
figure 2

The data shown was generated for \(H={{\mathbb{0}}}_{d}\) and γt ∈ [0, 10−4]. The circled dots show the numerical results. The solid curve presents the expected analytical result given by (14).

Full size image

The same simulations were repeated for larger values of γt ∈ [5 × 10−4, 1 × 10−2] and \(H={{\mathbb{0}}}_{d}\). The AGIs were then computed and are shown in Fig. 3 alongside the linear infidelity predicted in (14).

Fig. 3: Average gate infidelities as a function of γt for \(H={{\mathbb{0}}}_{d}\).
figure 3

The data points show the computed values. The solid lines represent the linear theoretical behaviour from (13). Each colour/marker pair corresponds to a different value of d.

Full size image

For more insight, Fig. 4 shows the relative deviation of the computed infidelities from the expected first-order linear behaviour for a broader range of γt up to 5 × 10−2.

Fig. 4: Relative deviation \(1-\frac{{\overline{{{{{\mathscr{E}}}}}_{d}}}^{{{{\rm{sim}}}}}}{{\overline{{{{{\mathscr{E}}}}}_{d}}}^{{{{\rm{th}}}}}}\) as a function of γt for \(H={{\mathbb{0}}}_{d}\).
figure 4

The quantities \({\overline{{{{{\mathscr{E}}}}}_{d}}}^{{{{\rm{sim}}}}}\) and \({\overline{{{{{\mathscr{E}}}}}_{d}}}^{{{{\rm{th}}}}}\) were obtained from numerical computations and (14) respectively. Each marker corresponds to a different value of d.

Full size image

Average gate infidelities linear in γt with gradients \(\frac{d(d-1)}{12}\) were expected in the case of a single qudit under pure dephasing, according to (14). Figure 2 supports this for small values of γt: a least-squares fit of the computed gradients yields the expected relationship with 1 − R2 < 10−5. Simulations for larger values of γt (Fig. 3) highlight deviations from this linear behaviour. These originate from \({{{\mathcal{O}}}}({(\gamma t)}^{2})\) terms of the form (γt)k>1 (see Supplementary Notes 2b). Moreover, for fixed values of γt, as d increases, the amplitude of this deviation is observed to increase (Fig. 4). This implies that the range of γt values for which the AGI can be treated linearly diminishes with increasing qudit dimension. Assuming a prefactor of the order d4 for the (γt)2 term in the AGI series expansion (as Supplementary Notes 2b hints), this provides an estimate of the range for which the deviation from linearity is negligible: γt ≪ 1 and \(\frac{1}{{d}^{2}}\).

Gate dependence

While the linearity of the AGI does not scale well with d, Eq. (14) has another important characteristic that deserves to be studied: the gate independence of the AGI. This was investigated over a large number of random gates for a given dimension d. Random unitary quantum gates in U(d) were sampled from the circular unitary ensemble, which represents a uniform distribution over the unitary square matrices of dimension d, also known as the Haar measure on the unitary group U(d), and implemented on a qudit through a Hamiltonian obtained by gradient-ascent methods. We decided to model qudits as ladder systems, with one pulse per transition between adjacent levels as considered for example in the experiments of Godfrin et al.20, for a single-molecule magnet (TbPc2, qudit with d = 4), the d − 1 pulses are then each represented by a control Hamiltonian in the interaction picture. More details are discussed in the Methods section. There are a large number of parameters that can influence the results under consideration, such as the free-evolution Hamiltonian or the matrix form of the control pulses, both of which are inherent to the physical realization. Therefore other physical implementations and reference frames for the pulses can be considered, and the deviation from linearity they cause needs to be studied in more detail. The AGIs were then computed for γt ∈ [10−5, 10−3], which lie in the typical ranges observed in current platforms as seen in Table 2, and their rate of increase as a function of γt was fitted. Figure 5 shows the statistical distributions of the relative deviations from the linear behaviour of the obtained rates.

Fig. 5: Hamiltonian-induced deviation of infidelity gradients from linear behavior.
figure 5

Statistical distributions of the relative deviation from the linear behaviour in (14) of the numerically obtained infidelity gradients cd for Ng = 5000 gates for γt ∈ [10−5, 10−3], as a function of the dimension d ∈ ⟦3, 8⟧. The candlestick bar chart should be interpreted as indicated in the upper right, with σ denoting the standard deviation, and the error bars denoting the extremal values. The lower right inset shows the same results for d ∈ ⟦2, 4⟧.

Full size image

Considering Fig. 5, the relative deviation from the linear behaviour for different random gates seems to exceed 1‰ rarely and was not observed outside the < 1% range. Moreover, the range of deviations decreases as the dimension d increases. The inset highlights a noticeable irregularity for d = 2, a single qubit, where the relative deviation is of the order of 1%. Note that for the \(H={{\mathbb{0}}}_{d}\) case simulated in Fig. 2, this shift remained < 1‱, coinciding with the dashed line in Fig. 5, including the case d = 2. Beginning at d = 2, the gradient distributions appear broad and off-centred from the \(H={{\mathbb{0}}}_{d}\) case. As d increases further, the distributions become progressively concentrated around 0. The gate-dependence also arises from \({{{\mathcal{O}}}}({(\gamma t)}^{2})\) terms, with γt2 being dominant in the γt2∣∣H∣∣ ≳ 1 regime (see Supplementary Notes 2b). Therefore, the AGI can only be considered gate-independent when γt ≪ 1 and \(\gamma t\ll \frac{1}{| | H| | t}\). An informative figure showing the deviation from linearity and gate-dependence at higher values of γt is available in Supplementary Notes 2a.

Other cases than pure dephasing

Figure 6 shows AGI rates of increase for channels different from pure dephasing namely: \(\overline{{{{{\mathscr{E}}}}}_{d}}({{{{\mathcal{E}}}}}_{x})\), \(\overline{{{{{\mathscr{E}}}}}_{d}}({{{{\mathcal{E}}}}}_{{{{\rm{+}}}}})\), and \(\overline{{{{{\mathscr{E}}}}}_{d}}({{{{\mathcal{E}}}}}_{x,y,z})\) corresponding to bit-flip, amplitude damping and depolarizing channels respectively. The simulations were performed again with \(H={{\mathbb{0}}}_{d}\), small γt ∈ [0, 10−4] and for even dimensions d ∈ ⟦2, 22⟧.

Fig. 6: Rate of increase of \(\overline{{{{{\mathscr{E}}}}}_{d}}({{{\mathcal{X}}}})={c}_{d}(J)\gamma t\) as a function of qudit dimension.
figure 6

The data shown was generated for \(H={{\mathbb{0}}}_{d}\) and γt ∈ [0, 10−4]. The markers show the numerical results. The solid curves represent the expected linear responses according to (16). Each marker/colour pair corresponds to a different error channel \({{{\mathcal{X}}}}\), with collapse operators J specified in the legend.

Full size image

Consider \({{{\mathcal{R}}}}\), the unitary transformation representing a change of basis, such as a 3D real-space rotation. The average gate fidelity defined in (8) is invariant under the transformation \(\rho \to {{{{\mathcal{R}}}}}^{{\dagger} }\rho {{{\mathcal{R}}}}\). This is supported by a comparison of the results for L = Jz and L = Jx in Figs. 2 and 6, respectively, since the two gradients appear to share the same dependency in d. Moreover, let {lk} be an ensemble of traceless collapse operators with corresponding error channels {ek}, then define L = ∑klk and associated error channel \({{{\mathcal{E}}}}\). From (10) and (12), as long as \({{{\rm{Tr}}}}({l}_{k}^{{\dagger} }{l}_{j})=0,\,\forall j\,\ne\, k\) then \(\overline{{{{\mathscr{E}}}}}({{{\mathcal{E}}}})={\sum }_{k}\overline{{{{\mathscr{E}}}}}({e}_{k})\). Figure 6 again supports such behaviour since simulations with the collapse operator L = J+ ≡ Jx + iJy yield gradients twice as large as, and L = Jx + Jy + Jz yield gradients three times as large as, the L = Jz case.

A single qudit vs an ensemble of qubits

An ensemble of n qubits were simulated under identical pure dephasing, with \(H={{\mathbb{0}}}_{{2}^{n}}\) and small γt ∈ [0, 10−4]. The simulations were performed for n ∈ ⟦1, 7⟧. Fitting the AGI \(\overline{{{{{\mathscr{E}}}}}_{{{{\rm{b}}}},n}}({{{{\mathcal{E}}}}}_{z})={c}_{{{{\rm{b}}}},n}\gamma t\) as a function of γt yielded the slopes cb,n that are shown in Fig. 7, along with their analytical expression as a function of d given in (20).

Fig. 7: Rate of increase of \(\overline{{{{{\mathscr{E}}}}}_{{{{\rm{b}}}},n}}({{{{\mathcal{E}}}}}_{z})={c}_{{{{\rm{b}}}},n}(\{{L}_{k}\})\gamma t\) as a function of Hilbert Space dimension d = 2n.
figure 7

The data shown was generated for \(H={{\mathbb{0}}}_{{2}^{n}}\) and γt ∈ [0, 10−4]. The {Lk} collapse operators are the ones defined in (18). The circled dots show the numerical results. The solid curve presents the expected theoretical result according to (20). The dashed line shows \({{{{\mathscr{E}}}}}_{{{{\rm{b}}}},n}^{{{{\rm{(p)}}}}}({{{{\mathcal{E}}}}}_{z})\) given in (31) which is linear in \(n={\log }_{2}(d)\).

Full size image

The same simulations were performed on a single qudit with dimension d = 2n and Fig. 8 shows the ratios \(\frac{{c}_{d}}{{c}_{{{{\rm{b}}}},n}}\) for n ∈ ⟦1, 6⟧ as well as the theoretical curve provided by Eq. (22) on which the points should be falling. According to the same Eq. (22), this curve also highlights the critical values of τb/τd, denoting the figure of merit τk = tk/T2,k = γktk/2, with respect to qudit/qubits advantage in terms of the rate of increase of the AGI.

Fig. 8: Potential range for the ratio of figure of merits τb/τd.
figure 8

The rounded circles show the numerical values obtained for cd/cb,n. The solid curve comes from (22) and highlights the theoretical critical values of T2,d/T2,b.

Full size image

The AGI gradients obtained for an ensemble of n qubits under identical pure dephasing were expected to follow a \(\frac{d{\log }_{2}(d)}{4(d+1)}\) relationship as a function of d(20). Figure 7 justifies this for small values of γt with the least-square fit now yielding 1 − R2 < 10−7. Finally, Fig. 8 provides quantitative data for the ratio of the decoherence times of a single qudit vs an ensemble of qubits. Some values of interest are summarised in Table 1. For example, in order for a qu-8-it (qudit with d = 8) to present a computational fidelity advantage over 3 qubits for a fixed gate time, the qudit platform needs a coherence time at least 7 times longer than the multiqubit platform. Note that an intuitive scaling such as \(\frac{{d}^{2}}{{\log }_{2}(d)}\) would indicate a much more demanding constraint of 21.5.

Table 1 Ratios of gate times in units of decoherence times between qubits and qudits for specific values of n and d (τb/τd needs to be larger than the critical values in order for a single qudit to be advantageous vs an equivalent ensemble of n qubits)
Full size table

From Table 2, state-of-the-art single qudit platforms, such as trapped ions4, present coherence times of the order of 100 ms for a single qu-7-it, orders of magnitude longer than superconducting qubits33,34,35,36. Trapped ions present γt ≈ 10−3, while γt ≈ 10−2 for superconducting qubits; this ratio of 10 would allow qudits with d ≲ 10 to still be advantageous i.e., according to (22), the single qu-7-it would still maintain a higher average gate fidelity over one gate acting on the whole Hilbert Space than the multiqubit platform. Another comparison with superconducting qubits could be molecular nuclear spin qudits, where some proposals put γt ≈ 10−4 (see Moreno-Pineda et al.1), and whose coherence times are ~ 6 − 7 times larger than the superconducting qubit case. With figure of merits τ ~ 100 larger than superconduction platforms, single molecular nuclear spin qudits with d ≲ 40 are still advantageous over equivalent superconducting qubits, i.e. n ~ 5. Such high-d qudit platforms can still be conceivable, given that some specific quantum operations on d = 52 have already successfully been implemented on, for example, Rydberg atoms37. However, it remains to be seen if universal quantum gate generation will become easily achievable in practice with such high d.

Table 2 Decoherence times (T2) and gate times (tn) of different qubit/qudit platforms
Full size table

Finally, one can compare (22) and (24) to discuss conditions on N qudits outperforming \(N{\log }_{2}(d)\) qubits. From this, if a single qudit outperforms \({\log }_{2}(d)\) qubits, the advantage remains conserved as long as the multiqudit gate time scales slower from 1 qudit to N qudits than the multiqubit from \({\log }_{2}(d)\) to \(N{\log }_{2}(d)\) qubits.

Perspectives on scalability and pathways beyond the first order

Given the rapid development of quantum computing platforms with very different physical properties, such as decoherence time or Hilbert space dimension (see Table 2), there is a growing need for detailed elaboration of the tradeoffs between their information density and noise error rates. By combining analytical results and numerical simulations, we have performed a comparative study of gate efficiency for systems composed of sets of qubits or qudits. A fluctuation-dissipation-like relation for the gate infidelity of an operation on a pure state was derived. We then put forward a physically-informed method to obtain the first-order effect of Markovian noise on the average gate infidelity (AGI). A connection was made between the latter and the first gate-independent result. The rate of increase of the AGI of a single qudit vs equivalent multiple qubits under pure dephasing was compared. This yielded a critical curve of the ratio of their respective gate times in units of decoherence time, a quantity indicating how time-efficient operations on a particular system are. Values on either side of the curve specify which of the two systems had a higher rate of increase of the AGI. To compete in terms of gate fidelity, as the dimension increases, the efficiency of qudit gates must not simply always be larger than the multiqubit one by a factor \(O({d}^{2}/{\log }_{2}(d))\), but precisely by a factor \(\frac{{d}^{2}-1}{3{\log }_{2}(d)}\), which makes a significant difference for lower values of d for which it provides less demanding constraints. Additionally, analytical expressions of linear response for arbitrary collapse operators and a general multiqudit system were presented (see (16) and (25)). They may be useful to those working in the field of quantum computing e.g., to, as mentioned in the introduction, benchmark qudit platforms either in terms of maximal practical d or in terms of conditions on the figure of merit to compensate for the greater noise scaling, in comparison with current state-of-the-art multiqubit platforms.

Numerical simulations contributed to the discussion on the validity and limits of the linear response assumption. This further restricted the ranges of possible γt ≪ 1 accounting for qudit dimension, gate, and noise type. For example, the larger the dimension, the lower the relative gate-dependent response. Finally, after simulations supported the analytical critical curve, different current platforms were studied with respect to this condition on gate time efficiency. Given equivalent Hilbert space dimensions, viable qudit platforms (leveraging advantageous decoherence times and gate speeds to compensate for the higher rate of increase in AGI) capable of outperforming equivalent state-of-the-art multiqubit ones in gate fidelity have been found for pure dephasing. Moreover, this performance could be extended to qudits with d as large as ~ 40 in the case of nuclear spins in molecular magnets, for example. Some multiqubit platforms still outperform any existing qudit platform regarding scalability in the number of subsystems. However, it is conceivable that some scalable qudit platforms continue to outperform equivalent multiqubit systems in terms of attainable fidelity. Further study of how multiqudit and multiqubit gate times scale with the number of subsystems is needed. Moreover, this study was limited to first-order noise responses. However, using the notation “qu-j-it" for a qudit of dimension d = j, through carefully chosen quantum error correction schemes, it was recently demonstrated it is possible to entirely remove the first-order response of logical qu-k-its embedded in physical qu-d-its (k < d) through carefully chosen encodings12,13. Of particular interest in the authors’ future work is the study of the Hamiltonian-dependent response of the dimension-dependent AGI. More generally, an additional study of higher-order responses of logical qudits vs. physical qubits would therefore also be required to assess the viability of these logical error-resilient qudits. This could elucidate if the presence of single qudit advantage, as quantified by this paper, is robust to system scaling and if qudits will remain useful beyond the NISQ era.

Methods

Numerical noisy qudit/multiqubit simulation

The data and source code used for generating the results shown in this paper are publicly available on the RADAR4KIT Repository38. All simulations were done using the Python package QuTiP24 version 4.7, SciPy version 1.7.3, and NumPy version 1.21.5. This subsection aims to present the modus operandi for obtaining the “numerical results" referenced in the different figures: the AGIs (\(\overline{{{{\mathscr{E}}}}}\)) for Figs. 3 and 4 and the slope of the AGIs (c) in the other figures.

Standard packages

Essential functions for simulating quantum dynamics and fitting curves to data are provided by the QuTiP library, including functions for propagator calculation in superoperator form (qt.propagator) and gate fidelity evaluation (average_gate_fidelity from qutip.metrics), along with the curve_fit function from scipy.optimize.

Parameters

We define the system’s dimension d, the decay parameter γ, and the collapse operators {Lk} under consideration. The collapse operators are QObj instances characterized by their matrix form in the canonical basis. Additionally, we generate a list of time points for simulating the system’s evolution. Considering the quantity of interest in this study is γt, γ is chosen as fixed, and the range of γt is then given by the range of the time points.

Time evolution

The simulation of the quantum system’s time evolution is facilitated by computing the propagator using the system’s Hamiltonian, the list of time points, and the collapse operators multiplied by \(\sqrt{\gamma }\). This generates a time-dependent propagator in the form of a list of QObj superoperators for different values of γt. The system’s Hamiltonian will be discussed in further detail in the following subsection. However, apart from Fig. 5 studying the gate/Hamiltonian-dependence, the other figures report simulations done with a vanishing Hamiltonian \(H={{\mathbb{0}}}_{d}\) since the quantities under consideration are considered Hamiltonian-independent.

Fidelity calculation

At each γt the average gate fidelity is computed relative to a target gate, in the case of \(H={{\mathbb{0}}}_{d}\): the identity matrix. This is the quantity displayed in Figs. 3 and 4.

Curve fitting

A curve is fitted to the calculated fidelities over the range of γt using the curve_fit function. This process involves fitting the function 1 − cγt for the parameter c. The obtained slopes c({Lk}) are then the ones displayed in the Figs. 2, 5, 6, 7 and 8. Moreover the least-square fit parameter R given by the fitting functions is the one reported in this study.

Random gate and pulse Hamiltonian generation

Gate generation

In the study of the gate-dependent deviation from the analytical results of this manuscript, for each dimension d under consideration, a set of Ng = 5000 gates have been randomly generated with the Bristol39 package in Python. The gates have been drawn from the circular unitary ensemble, and are thus considered to be uniformly distributed over the Haar measure. Subsequently, to generate an associated set of pulses for each gate, we have used the optimize_pulse_unitary function from the pulse optimization module (control.pulse_optim) of QuTiP.

Pulse generation

The pulse generation is done through gradient-ascent methods using the GRAPE algorithm40 and was run in parallel for each gate using a high-performance cluster. The numerical optimizer used by default is the L-BFGS-B method. Assuming the control Hamiltonian, as discussed after (1), takes the form

$$H{{{\rm{c}}}}(t)=\mathop{\sum }\limits_{k=0}^{N}{u}_{k}(t){H}_{k},$$
(32)

with Hk being a basis set of controls and uk(t) representing the time-dependent control amplitudes, the optimization process involves finding the set of uk(t) that best approximates the target gate.

Choice of Hamiltonian

For the simulations reported in this paper, we decided to model qudits as ladder systems, with one pulse per transition between adjacent levels as considered for example in the experiments of Godfrin et al.20, for a single-molecule magnet (TbPc2, qudit with d = 4), the d − 1 pulses are then each represented by two control Hamiltonians in the interaction picture. More explicitly, the basis set of controls is chosen to be the ensemble of pairs \(\left\vert k\right\rangle \left\langle k+1\right\vert +\left\vert k+1\right\rangle \left\langle k\right\vert\) and \({{{\rm{i}}}}(\left\vert k\right\rangle \left\langle k+1\right\vert -\left\vert k+1\right\rangle \left\langle k\right\vert )\), with k running from 1 to d − 1. Moreover, H0, the free-evolution, is chosen to be vanishing since we consider the interaction reference frame.