# Nonlinearity of the Fidelity in Open Qudit Systems: Gate and Noise Dependence in High-dimensional Quantum Computing

Jean-Gabriel Hartmann, Denis Janković, Rémi Pasquier, Mario Ruben, Paul-Antoine Hervieux

► To cite this version:

Jean-Gabriel Hartmann, Denis Janković, Rémi Pasquier, Mario Ruben, Paul-Antoine Hervieux. Nonlinearity of the Fidelity in Open Qudit Systems: Gate and Noise Dependence in High-dimensional Quantum Computing. Quantum, 2025, 9, pp.1690. 10.22331/q-2025-04-07-1690. hal-05090542

HAL Id: hal-05090542

https://hal.science/hal-05090542v1

Submitted on 30 May 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Nonlinearity of the Fidelity in Open Qudit Systems: Gate and Noise Dependence in High-dimensional Quantum Computing

Jean-Gabriel Hartmann<fnref n="1" />, Denis Janković<fnref n="1" />,<fnref n="2" />, Rémi Pasquier<fnref n="1" />, Mario Ruben<fnref n="2" />,<fnref n="3" />,<fnref n="4" />, and Paul-Antoine Hervieux<fnref n="1" />

High-dimensional quantum computing has generated significant interest due to its potential to address scalability and error correction challenges faced by traditional qubit-based systems. This paper investigates the Average Gate Fidelity (AGF) of single qudit systems under Markovian noise in the Lindblad formalism, extending previous work by developing a comprehensive theoretical framework for the calculation of higher-order correction terms. We derive general expressions for the perturbative expansion of the Average Gate Infidelity (AGI) in terms of the environmental coupling coefficient and validate these with extensive numerical simulations, emphasizing the transition from linear to nonlinear behaviour in the strong coupling regime. Our findings highlight the dependence of AGI on qudit dimensionality, quantum gate choice, and noise strength, providing critical insights for optimising quantum gate design and error correction protocols. Additionally, we utilise our framework to identify universal bounds for the AGI in the strong coupling regime and explore the practical implications for enhancing the performance of near-term qudit architectures. This study offers a robust foundation for future research and development in high-dimensional quantum computing, contributing to the advancement of robust, high-fidelity quantum operations.

## 1 Introduction

High-dimensional quantum computing (QC) has generated remarkable scientific interest of late, introducing a shift from traditional computing paradigms. While qubit-based quantum information processing

Jean-Gabriel Hartmann: jeangabriel.hartmann@ipcms.unistra.fr

Denis Janković: denis.jankovic@ipcms.unistra.fr

Rémi Pasquier: rpradeep@outlook.fr

Mario Ruben: mario.ruben@kit.edu

Paul-Antoine Hervieux: hervieux@unistra.fr

(QIP) platforms, particularly superconducting qubits, have the highest technological maturity, they are faced with near-term technical challenges of scalability and error correction [1, 2]. On the other hand, recent advancements have highlighted the potential of qudits — quantum systems with *d* levels — as powerful alternatives for novel QC architectures [3, 4, 5, 6].

Indeed, qudits offer several advantages, including (i) lower decoherence rates in certain physical systems [7], (ii) enhanced quantum error correction through additional levels [8, 9, 10] and stabiliser codes [11], as well as (iii) higher information density for reducing circuit complexity and enabling novel algorithm design [12]. They also promise more robust flying quantum memories [13, 14].

While classical computing ultimately settled on bits once sufficient scalability and fault-tolerance were achieved, early platforms did experiment with multi-level systems [15]. Analogously, it can be argued that quantum computing's infancy stage could benefit from exploring higher-dimensional bases to address these current issues. In the push towards universal QC, increasing the total Hilbert space dimension of quantum systems is critical [16]. This dimension is determined by *d*ⁿ, where *d* is the dimensionality and *n* the number of qudits. Despite impressive advancements in superconducting platforms [17], scaling the number of qubits continues to pose significant challenges. Thus, qudit-based approaches may provide an avenue for increasing the Hilbert space with fewer physical units [17, 6, 18, 19].

On the other hand, with the increased number of excited states utilised to implement higher-dimensional states in physical systems, qudits may introduce a greater number of error channels compared to qubits. This could lead to increased sensitivity to environmental noise, affecting coherence times and complicating error correction processes [20]. In a previous paper [21], we studied this scenario in detail by comparing the effects of Markovian noise on multi-qubit and single-qudit systems of equal Hilbert space dimension. We developed a theoretical model

<footnote n="1">¹Université de Strasbourg, CNRS, Institut de Physique et Chimie des Matériaux de Strasbourg, UMR7504, F-67000 Strasbourg, France</footnote>

<footnote n="2">²Institute of Nanotechnology (INT), Karlsruhe Institute of Technology, P.O. Box 3640, 76021 Karlsruhe, Germany</footnote>

<footnote n="3">³Institute for Quantum Materials and Technologies (IQMT), Karlsruhe Institute of Technology, P.O. Box 3640, 76021 Karlsruhe, Germany</footnote>

<footnote n="4">⁴Université de Strasbourg, CNRS, Centre Européen de Science Quantique - Institut de Science et d'Ingénierie Supramoléculaires, UMR 7006, F-67000 Strasbourg, France</footnote>
of the first-order effects of the noise in the Lindblad formalism, demonstrating how the noise impacts the performance of quantum gates in these different systems, and supported these results with numerical simulations.

The fundamental quantities investigated were the average gate fidelity (AGF, $\bar{F}$) - equivalently infidelity (AGI, $\bar{I}$) - and the figure-of-merit ($\tau := t_{\text{gate}}/T_2$). As the name suggests, the AGF is useful in that it integrates out any specific features due to a particular choice of initial state. Therefore, unlike the state fidelity or process fidelity, which compute the target-to-output overlaps of quantum state transfers and unitary gates, respectively, the AGF is fundamental as a platform-agnostic measure of the quality of a system's interaction with its environment [22, 23]. Similarly, the figure-of-merit is important for quantifying realistically the circuit depth that can be achieved in a system, based on the gate time, $t_{\text{gate}}$, and decoherence time, $T_2$.

Those results showed that in the quasi-error-free weak coupling regime, the first-order response of the AGI (a linear function of the dimensionless coupling strength $\gamma t_{\text{gate}}$) was sufficient to characterise the behaviour of the qudit or multi-qubit system. Specifically, a scaling law was derived, showing how the AGI is affected by qudit dimension and the conditions on the figure-of-merit under which qudits may be competitive with multi-qubit systems. These conditions lead to stricter requirements on the figure-of-merit as the qudit dimension increases. However, first-order approximations are insufficient for fully understanding the behavior of qudits under realistic, near-term, conditions where higher-order noise effects are expected to become significant, either through stronger coupling or longer gate times (or circuit depths) approaching the decoherence limit.

Thus, this paper aims to extend the understanding of the AGF of single qudit systems under Markovian noise conditions through a general perturbative expansion. Specifically, we focus on higher-order correction terms and their implications on quantum gate performance in the strong coupling regime. Building upon our previous work, we develop a comprehensive theoretical framework that includes these nonlinear effects, with a constructive method generalising the AGF to arbitrary order in $\gamma t$. This is supported by detailed numerical simulations, with an emphasis on the differences between the first and second order terms in the case of pure dephasing. By doing so, we aim to capture the nuanced impact of noise on the fidelity of quantum gates more precisely.

The key research question addressed in our study is: How do higher-order correction terms and noise coupling strength influence the AGI of single qudit systems, and what implications do these have for the selection and design of quantum gates? Our findings provide important insights on: (i) setting bench-

marks, as well as limits, for the performance of noisy quantum systems of arbitrary dimension by identifying gate-dependent effects, (ii) methodologies for improving performance through optimising basis gates that have favourable decoherence characteristics, and (iii) advancing error correction protocols by enabling the cancellation of higher-order noise effects.

The paper is structured as follows: We begin in Sec. 2.1 with a review of the relevant theoretical background concerning the Lindblad formalism for qudit open quantum systems and the superoperator representation of noisy quantum channels, as well as generalised quantum gates on qudits. In Sec. 2.2, we present our general result for the perturbative expansion of the AGI in noise strength, deriving the general expressions to arbitrary order in $\gamma t$ and discussing the implications of its structure. Section 2.3 examines the nonlinearity of the AGI in the strong coupling regime of pure dephasing through numerical simulations that elucidate key characteristics of this behaviour and motivate the proceeding investigations that make the link with the theoretical results. Subsequently, in Sec. 2.4 we explore the gate- and noise-dependence of the AGI, with particular emphasis on how the choice of quantum gate can influence the performance of the quantum channel. Finally, Sec. 2.5 presents a detailed analytical and numerical study of the first- and second-order correction terms to the AGI, extending the results of the previous paper and placing them on a more rigorous theoretical foundation that serves future investigations.

# 2 Results and Discussion

## 2.1 Mathematical Foundations of the Average Gate Fidelity

We begin with the study of a single qudit, of arbitrary dimension $d$, whose state is represented by a $d \times d$ density matrix $\rho(t)$, and evolves according to the Gorini-Kossakowski-Sudarshan-Lindblad (GKSL) master equation [24, 25], which may be expressed in superoperator form as [26]

$$ \partial_t \rho = S[\rho] + D[\rho]. \quad (1) $$

Here, setting $\hbar = 1$, the Liouvillian superoperator $S$ represents the unitary evolution operator of the von Neumann equation, and $D$ the dissipation superoperator coupling the qudit to the environment,

$$ S[\rho] := -i[H, \rho], \quad (2) $$

$$ D[\rho] := \sum_{k=1}^{K} \gamma_k \mathcal{L}_k[\rho] \quad (3) $$

$$ = \sum_{k=1}^{K} \gamma_k \left( L_k \rho L_k^{\dagger} - \frac{1}{2} \{ L_k^{\dagger} L_k, \rho \} \right), \quad (4) $$
where $H$ is the time-independent interaction Hamiltonian, and $L_k$ one of $K$ possible collapse operators characterising the Markovian noise with coupling coefficient $\gamma_k$ and $\mathcal{L}_k$ the matrix superoperator form [27].

The solution to the master equation for $\rho(0) = \rho_0$ is therefore given by the completely-positive and trace-preserving (CPTP) quantum channel $\mathcal{E}$ as [28]

$$ \rho(t) = \mathcal{E}[\rho_0] \qquad (5) $$

$$ = e^{(S+D)t} \rho_0. \qquad (6) $$

Considering the application of a quantum gate $\mathcal{U}[\rho] = U\rho U^\dagger$ to the system then corresponds to a composition operation, that in the matrix superoperator representation reduces to the matrix product, from which we may express the AGF, $\bar{\mathcal{F}}(\mathcal{E}, \mathcal{U})$, of implementing a unitary operator over the quantum channel in the presence of a noisy environment [23]:

$$ \bar{\mathcal{F}}(\mathcal{E}, \mathcal{U}) := 1 - \bar{\mathcal{I}}(\mathcal{E}, \mathcal{U}) \qquad (7) $$

$$ = \int_{\mathcal{H}} \langle U^\dagger \mathcal{E} [\rho_0] U \rangle_0 d\rho_0 \qquad (8) $$

$$ = \int_{\mathcal{H}} \langle (U^\dagger \circ \mathcal{E}) [\rho_0] \rangle_0 d\rho_0, \qquad (9) $$

where the matrix product has been rewritten in the superoperator representation. The integral is taken uniformly over the Haar measure $\mathcal{H}$ of the state space, with the expectation value $\langle A \rangle_0 = \mathrm{Tr}\{A\rho_0\}$ representing the average over all initial states. For convenience, we shall use $\bar{\mathcal{F}}$ and $\mathcal{F}$ interchangeably for the AGF, and likewise for the AGI.

In this study we shall consider pure initial states, $\mathrm{Tr}\{\rho_0^2\} = 1$. Then the state space is reduced to the complex-projective space $\mathbb{CP}^{d-1}$ such that the Haar measure (over all $\rho_0$) then induces the Fubini-Study measure, the integral of which is normalised, $\int_{\mathcal{H}} d\rho_0 = 1$, [29] and may be calculated analytically (see Appendix A.3). Furthermore, unless specified, we shall also assume that: (i) the dissipator term is dominated by a single noise channel, such that $\mathcal{D} = \gamma\mathcal{L}$, (ii) the associated superoperators $\mathcal{S}$ and $\mathcal{L}$ are time-independent and (iii) we are working in the interaction picture, such that the free evolution Hamiltonian $H_0 = \mathbb{0}_d$, and therefore $\mathcal{S}_0 = \mathbb{0}_{d \times d}$ such that the unitary evolution simplifies to $\mathcal{U} = e^{(S_0+S_c)t} = e^{St}$ with $\mathcal{S} = \mathcal{S}_c$ corresponding to only the control (driving) Hamiltonian $H = H_c$ implementing the quantum gate. This allows us to rewrite Eqs. (5) and (6) as

$$ \mathcal{E}[\rho_0] = e^{(S+\gamma\mathcal{L})t} \rho_0. \qquad (10) $$

We note that the $(d \times d) \times (d \times d)$ superoperators acting on the $\mathcal{H} \otimes \mathcal{H}$ Hilbert space of the quantum dynamical semigroup (QDS) may be written explicitly in Liouville matrix form in terms of the standard

operator form as [30]:

$$ \mathcal{U}[\cdot] := (U^* \otimes U) \mathrm{vec}(\cdot), \qquad (11) $$

$$ \mathcal{S}[\cdot] := -i (\mathbb{1}_d \otimes H - H^T \otimes \mathbb{1}_d) \mathrm{vec}(\cdot), \qquad (12) $$

$$ \mathcal{L}[\cdot] := \left( L^* \otimes L - \frac{1}{2} \left( L^\dagger L \otimes \mathbb{1}_d + (\mathbb{1}_d \otimes L^\dagger L)^T \right) \right) \mathrm{vec}(\cdot), \qquad (13) $$

where we identify the adjoint ($\dagger$), complex-conjugate ($\ast$) and transpose ($T$) operations, $\mathbb{1}_d$ the $d \times d$ identity matrix, $\otimes$ the tensor, or equivalently here for linear maps, the Kronecker product, and $\mathrm{vec}(\cdot)$ the vectorisation operation that, when applied to a $d \times d$ density matrix $\rho$, produces a $1 \times d^2$ column vector by stacking vertically each column of the matrix from left to right. This corresponds to the explicit mapping

$$ \mathrm{vec}(\rho) := \sum_{i,j} \rho_{ij} |i\rangle\langle j| \rightarrow \sum_{i,j} \rho_{ij} |j\rangle \otimes |i\rangle, \qquad (14) $$

of the quantum state from the Hilbert to the Fock-Liouville space [27]. Furthermore, the quantum gate $U_g$ and associated control Hamiltonian $H_c$ (and their respective superoperator representations $\mathcal{U}$ and $\mathcal{S}$) are related by the matrix exponent (and, inversely) logarithm operations. Therefore, given some quantum gate, it may be implemented by the control term,

$$ H_c t_g = i \log U_g, \qquad (15) $$

where the now time-independent $H_c$ represents a single, ideal control pulse that implements the gate in time $t_g$, which can be modulated by the absolute amplitude of the pulse $\|H_c\|$. Existence and uniqueness of the matrix logarithm map are well established under the conditions that $U_g$ is unitary and $-\pi < H_c t_g \le \pi$ [31]. This is an idealisation that necessitates simultaneous control over all pulse amplitudes, phases and detunings on each possible transition between all $d$ states. Platform-dependent physical constraints on the control pulses may preclude experimental realisation of such a control term. Nevertheless it allows for a time-optimal multi-chromatic pulse for studying the robustness of quantum gates that is independent of constraints imposed by the physical platform or pulse control technique [32].

We now precisely define the gate and collapse operators that will be of interest in this study. We begin with the following set of single-quitit quantum gates $U_g \in \{\mathbb{1}, X, Z, F\}$ that generalise the single-qubit identity, Pauli-x, Pauli-z and Hadamard gates to $d$ dimensions, and generate the $d$-dimensional gen-
eralised Clifford algebra (GCA) [33]

$$1_d := \sum_{j=0}^{d-1} |j\rangle\langle j|, \qquad (16)$$

$$X_d := \sum_{j=0}^{d-1} |(j+1) \mod d\rangle\langle j|, \qquad (17)$$

$$Z_d := \sum_{j=0}^{d-1} \omega^j |j\rangle\langle j|, \qquad (18)$$

$$F_d := \frac{1}{\sqrt{d}} \sum_{j=0}^{d-1} \sum_{k=0}^{d-1} \omega^{jk} |j\rangle\langle k|, \qquad (19)$$

where $\omega := e^{\frac{2\pi i}{d}}$ is the d-th root of unity. The subscript d is used to specify the dimension of the operator, and may be omitted in cases where it is clear from context or when referring to the general operator. The X, Z and F gates retain unitarity (and tracelessness for X and Z) as well as the relations $X^d = Z^d = 1_d$ and $X = FZF^\dagger$ of involution and change-of-basis, respectively. However, it is worth noting that they are no longer Hermitian for $d > 2$. Furthermore, in the generalised forms, it can be observed that (i) X is extended from a bit-flip (NOT) gate to a cyclic permutation (INC/ SHIFT) gate that increments the qudit state, (ii) Z is extended from a phase-flip to a phase-shift (CLOCK) gate of each state over the d roots of unity, and (iii) the Walsh-Hadamard gate is extended to the matrix form of the Quantum Fourier Transform (QFT) (equivalently referred to as the Discrete Fourier Transform, Sylvester or Chrestenson matrix).

In addition to this discrete set of Clifford gates, we can also consider non-Clifford gates that may be useful for quantum error correction (QEC) protocols for universal QC. In general, since the unitary group $U$ is not finite, any unitary matrix sampled randomly from, for example, the Haar measure over $U$ will almost always be non-Clifford. Therefore, by including in our investigations the Haar-random gates enables comparison of the differences in behaviour of the Clifford and non-Clifford gates. More specifically, one noteworthy non-Clifford gate is the (generalised) $\pi/8$ gate [34],

$$T_d := \sum_{j=0}^{d-1} e^{i \frac{j \pi}{8}} |j\rangle\langle j|, \qquad (20)$$

which when considered together with the $Z_d$ gate, we may extend to a generalised-PHASE gate

$$\text{PHASE}_d(\phi) := \sum_{j=0}^{d-1} e^{ij\phi} |j\rangle\langle j|. \qquad (21)$$

Setting $\phi$ to $\frac{\pi}{8}$, $\frac{2\pi}{d}$ or $\frac{\pi}{4}$ recovers $T_d$, $Z_d$ or the $S_d = \sqrt{Z_d}$ gates, respectively, where we may note that

special cases including $Z_d$ and $S_d$ are Clifford, while in general for arbitrary choices of phase, including $\phi = \frac{\pi}{8}$, will likely be non-Clifford (see Appendix B.1 for further details on the gate definitions).

Regarding the choice of Markovian collapse operators $L \in \{J_z, J_x, J_\perp\}$ representing the coupling of the qudit to the environment in the master equation, we shall primarily consider the effect of pure dephasing by the operator $J_z$, and of secondary consideration are the bit-flip error and spin relaxation effected by the operators $J_x$ and $J_\perp$, respectively. These higher-order spin operators from spin-$\frac{1}{2}$ to spin-$\frac{d-1}{2}$ can be obtained through the generalised Gell-Mann matrices, and written as

$$J_z := \sum_{j=1}^{d} \frac{1}{2} (d + 1 - 2j) |j\rangle\langle j|, \qquad (22)$$

$$J_x := \sum_{j=1}^{d-1} \frac{1}{2} \sqrt{j(d-j)} (|j+1\rangle\langle j| + |j\rangle\langle j+1|), \qquad (23)$$

$$J_\perp := \sum_{j=1}^{d-1} \sqrt{j(d-j)} (|j+1\rangle\langle j|), \qquad (24)$$

where $J_\perp = J_x - iJ_y$. It is also possible to consider configurations of linear combinations of these operators, which can be used to model more complex and realistic environmental interactions. However, for the purposes of this study, we are interested in understanding their individual interaction with the control Hamiltonian and consequent effect on the AGI in a platform-agnostic way and thus shall only consider each of the operators separately.

## 2.2 Perturbative Expansion of the AGI

Given the simplified expression for the time-evolution of the quantum channel $\mathcal{E}$ in (10), we want to study the effects on the AGF, not only from the choice of $S$ or $\mathcal{L}$, but also from the coupling strength $\gamma$. Therefore, performing a Taylor series expansion of the AGF in powers of $\gamma$,

$$\bar{\mathcal{F}}(\mathcal{E}, \mathcal{U}) = 1 - \sum_{m=1}^{\infty} \gamma^m F^{(m)}(t). \qquad (25)$$

To calculate explicitly the correction terms $F^{(m)}$, we begin by expanding the solution to the master equation,

$$\rho(t) = \rho^{(0)} + \sum_{m=1}^{\infty} \gamma^m \rho^{(m)}, \qquad (26)$$

where the m-th order correction to the non-perturbed solution of the master equation can be written as [35, 36]
$$ \rho^{(m)} = e^{\mathcal{S}t} \left( \int_0^t \int_0^{t_1} \cdots \int_0^{t_{m-1}} \left( \prod_{i=1}^m e^{-\mathcal{S}t_i} \mathcal{L} e^{\mathcal{S}t_i} \right) dt_m \cdots dt_2 dt_1 \right) \rho_0, \quad (27) $$

and we separate the zero order term $\rho^{(0)}$ from the summation, since it corresponds to only unitary evolution through the von Neumann equation in the absence of environmental coupling.

We may then define a new quantity $\tilde{M}^{(m)}(t)$ containing the integral terms, such that $\rho^{(m)} := e^{\mathcal{S}t} (\tilde{M}^{(m)}(t)) \rho_0$, and therefore the m-th order of the solution to the master equation in terms of the quantum channel reduces to $\mathcal{E}^{(m)} = \mathcal{U} \circ M^{(m)}(t)$. Now, in order to render this term $\tilde{M}^{(m)}(t)$ into a more computationally tractable form, we make use of the eponymous lemma of Campbell resulting from the Baker-Campbell-Hausdorff (BCH) Formula [37, 38],

$$ e^{X}Y e^{-X} = \sum_{n=0}^{\infty} \frac{[(X)^{n}, Y]}{n!}, \quad (28) $$

utilising the iterated commutator defined by the recursion relation $[(X)^n, Y] = [X, [(X)^{n-1}, Y]]$ with the halting condition $[(X)^0, Y] = Y$. With $\mathcal{S}$ and $\mathcal{L}$ time-independent, $\tilde{M}^{(m)}(t)$ simplifies to

$$ \tilde{M}^{(m)}(t) = t^m \sum_{n_1=0}^{\infty} \sum_{n_2=0}^{\infty} \cdots \sum_{n_m=0}^{\infty} \left( \prod_{i=1}^{m} \frac{(-t)^{n_i} [(\mathcal{S})^{n_i}, \mathcal{L}]}{n_i! \sum_{j=i}^{m} (n_j+1)} \right) \quad (29) $$

$$ = t^m \frac{\mathcal{L}^m}{m!} + t^m \sum_{n_1=1}^{\infty} \sum_{n_2=1}^{\infty} \cdots \sum_{n_m=1}^{\infty} \left( \prod_{i=1}^{m} \frac{(-t)^{n_i} [(\mathcal{S})^{n_i}, \mathcal{L}]}{n_i! \sum_{j=i}^{m} (n_j+1)} \right). \quad (30) $$

For the sake of future convenience, we define the new quantity $\tilde{M}^{(m)}(t) := t^m M^{(m)}(t)$, from which we obtain the following general form for the correction terms of the AGF in Eq. (25) in terms of Eq. (10):

$$ F^{(m)}(t) = - \int_{\mathcal{H}} \langle \mathcal{U}^{\dagger} \circ \mathcal{U} \circ (t^m M^{(m)}(t)) [\rho_0] \rangle_0 d\rho_0 \quad (31) $$

$$ = -t^m \int_{\mathcal{H}} \mathrm{Tr} \{ M^{(m)} [\rho_0] \rho_0 \} d\rho_0. \quad (32) $$

This integral over the Fubini-Study measure of pure initial states may now be evaluated using results from Weingarten calculus [39] (see Appendix A.3). From here we obtain our main result for the general perturbative expansion of the AGI, $\mathcal{I} = 1 - \mathcal{F} = \sum_{m=1}^{\infty} \mathcal{I}^{(m)}$:

$$ \mathcal{I} = \frac{-1}{d(d+1)} \sum_{m=1}^{\infty} (\gamma t)^m \mathrm{Tr}\{M^{(m)}(t)\}, \quad (33) $$

where $d$ is the dimensionality of the qudit system, and the powers of $(\gamma t)$ have been grouped to form a single, dimensionless quantity.

As an illustrative example of the utility of this expression, in Fig. 1 we present the cumulative correction terms up to fourth order in $\gamma t$ for the QFT gate $F_4$ acting on a qudit of $d=4$. With $t=1$ normalised by the amplitude of the control Hamiltonian $H_c$, the pure dephasing operator $L = J_z$ was applied over the

couplings $\gamma t \in [0, 0.5]$. While the first-order approximation shows good agreement with the simulated AGI for small $\gamma t \ll 1$, the deviation becomes significant as the coupling grows large, necessitating the addition of the higher-order corrections.

Looking closer at this general result for $\mathcal{I}$, let us consider the case of no coherent driving. This is implemented through a zero control Hamiltonian, $H = 0_d$, giving $\mathcal{S} = \mathbb{0}_{d \times d}$ and corresponding to evolution by the identity gate $U_g = 1_d$. Then, the system evolves only through decay due to the collapse operator $\mathcal{L}$, since only the iterated commutator terms $[(\mathcal{S})^{n_i}, \mathcal{L}]$ for $n_i = 0$ in $M^{(m)}(t)$ survive, resulting in

$$ \mathcal{I}_0^{(m)} = \frac{-(\gamma t)^m \mathrm{Tr}\{\mathcal{L}^m\}}{d(d+1) m!}, \quad (34) $$

which is precisely the first term in Eq. (30), with the subscript 0 used to refer to the gate-independent correction term. Therefore, generalising to nonzero $\mathcal{S}$, it is clear that the m-th order correction can always be separated into the sum of two parts as in Eq. (30): These gate-independent terms $\mathcal{O}((\gamma t)^m)$ (when all $n_i = 0$) and the product of the gate-dependent iterated commutators $\mathcal{O}((\gamma t)^m \prod_{i=1}^m t^{n_i})$.

Combining the summation over $m$ in Eq. (33) and the m-th order expression in Eq. (34) allows us to rewrite the gate-independent part as a Taylor expansion
<figure><img src="image_2.png" /><figcaption>Figure 1: Comparison of the first four correction terms of the AGI to numerical simulations. AGIs for the QFT gate applied to a $d = 4$ qudit undergoing pure dephasing ($L = J_z$) at couplings $\gamma t \in [0, 0.5]$. The discrete points (red) represent the simulated AGI values, while the successive dashed-dotted lines represent the correction terms from first to fourth order.</figcaption></figure>

sion,

$$ \mathcal{I}_0 = \frac{\mathrm{Tr}\{\mathbb{1}_{d \times d} - e^{\gamma t \mathcal{L}}\}}{d(d+1)}. \qquad (35) $$

Furthermore, as long as the collapse operator is real and symmetric, obeying $L = L^* = L^T = L^\dagger$, (as for $J_z$, for example), it is possible to use a multinomial expansion on the operator form of $\mathcal{L}$ in Eq. (13) to obtain explicit expressions to arbitrary order $(\gamma t)^m$ for the gate-independent correction in Eq. (34):

$$ \mathrm{Tr}\{\mathcal{L}^m\} = \mathrm{Tr}\left\{\left(L^* \otimes L - \frac{1}{2}L^\dagger L \otimes 1 - \frac{1}{2}1 \otimes L^\dagger L\right)^m\right\} \quad (36) $$

$$ = m! \sum_{k_1+k_2+k_3=m} \frac{\mathrm{Tr}\{L^{k_1+2k_2}\} \mathrm{Tr}\{L^{k_1+2k_3}\}}{k_1!k_2!k_3!(-2)^{k_2+k_3}}. \quad (37) $$

In the case of a pure (spin) dephasing channel $L = J_z$, considering that $\mathrm{Tr}\{J_z^m\} \sim d^{m+1}$, it can be easily seen that this expression for the $m$-th-order correction scales as

$$ \mathcal{O}(\mathrm{Tr}\{\mathcal{L}^m\}) \sim d^{2(m+1)} \quad (38) $$

$$ \Rightarrow \mathcal{O}(\mathcal{I}_0^{(m)}) \sim (\gamma t)^m d^{2m}. \quad (39) $$

Now, given this scaling of the gate-independent AGI, let us refer to the work of Chiesa et al. [40], where the authors proposed the use of a physical qudit to encode a logical qubit. Using their binomial encoding and error-correction protocol would result in a logical qubit robust to pure dephasing up to order $(\gamma t)^m$ where $2m+1 = d$. Effectively, the AGI of such a logical qubit would be of the same order as $\mathcal{I}_0^{(m)}$. In turn, to maintain $\mathcal{I}_0 \ll 1$, this sets a limit on the maximum

dimension (equivalently, order $m$ of noise correction), such that $d < \frac{1}{\sqrt{\gamma t}}$. Such investigations for protecting against decoherence have already been studied with alternative logical qubit encodings and a discrete error evolution protocol in the work of Miyahara et al. [10].

In what follows in the remainder of this study, we shall focus on the first- and second-order correction terms to the AGF. We note that the preceding analytical approach provides a framework for computing any correction term to arbitrary order. However, explicit calculations up to second-order serve as an extension of our prior work [21], in which the first-order correction to the AGF was investigated and the consequent conditions on gate efficiency derived in terms of the figure-of-merit as a function of qudit dimension. Thus, we study in detail the effects due to the second-order correction terms as they relate to our particular system of interest, particularly with regard to the choice of quantum gate and collapse operator.

## 2.3 Nonlinearity of the AGI in the Strong Coupling Regime

Let us now motivate our investigation into the AGI in the regime of strong coupling between the qudit and the environment. As identified in Fig. 1, the linear first-order approximation may be sufficient in the weak coupling regime, but nonlinear effects become significant above $\gamma t \sim 10^{-1}$. It is therefore necessary to identify a clear transition or separation between these two regimes, and how the AGI behaves in the limit of strong coupling.

To investigate, we perform calculations of the AGI over a large range of scales for $\gamma t$. Therefore, in Fig. 2, we show numerical simulations of the AGIs for qudits under pure dephasing ($L = J_z$) where the coupling strengths varied logarithmically in $\gamma t \in [10^{-5}, 10^3]$. Evolutions were repeated on qudits with dimensions $d \in \{2, 4, 8, 16\}$ for 100 Haar-random quantum gates (see Appendix C.1 for a discussion on methods for generating random quantum gates), with the inset showing the linear regime for small $\gamma t \in [10^{-5}, 10^{-2}]$ and the main figure showing the transition region and strong-coupling regime. The dashed-dotted straight lines give the first-order correction terms at each qudit dimension, defined by the linear function [21]

$$ \mathcal{I}^{(1)} = \frac{d(d-1)}{12}\gamma t, \quad (40) $$

showing the $d$-dependence in the weak coupling regime.

Further effects of the system dimension on the AGI are notable for larger couplings: The start of deviation from linearity, the plateau values $\mathcal{I}^*$ that the AGI curves converge to in the large $\gamma t$ limit, and the associated saturation points $(\gamma t)^*$ (both denoted by the superscript *) at which the curves reach their respective plateau values. Each of these effects confirms
<figure><img src="image_3.png" /></figure>

**Figure 2:** Deviation from linearity of the AGIs for qudits under strong dephasing. Simulations of the AGI against the strength of pure dephasing ($L = J_z$) are plotted on log-log axes over $\gamma t \in [10^{-2}, 10^3]$. The simulations were performed on a set of qudit dimensions, $d \in \{2, 4, 8, 16\}$ each for 100 Haar-random quantum gates. The dashed lines represent the linear regime given by the first-order correction to the AGI at each dimension. The inset shows the linear regime for smaller values of $\gamma t \in [10^{-5}, 10^{-2}]$. For stronger noise the AGIs exhibit a transition from the linear to nonlinear regime that then saturates at a stable plateau value. These plateau values vary for qudit dimension and also for different quantum gates, highlighting the gate-dependence. The saturation points $(\gamma t)^*$ of the stable regime are also dependent on both the qudit dimension and gate type, with higher-dimensional qudits deviating from linearity and saturating earlier.

an overall decrease in robustness as the system dimension increases: Compared to the behaviour for qubits, the AGIs (for arbitrary gates) of higher-dimensional qudits begin to deviate from the linear regime, and reach the higher plateau limits, at lower noise thresholds. These quantities shall be studied in more detail in Sec. 2.4.

Finally, it is interesting to note the gate dependence of the $I^*$ and $(\gamma t)^*$, which is more pronounced for smaller $d$. Given the sample size of 100 gates, it is apparent that, at each dimension, there exists an upper and lower bound to the AGI plateau values, with all the sampled gates converging within this range. It is also noticeable (particularly so for $d=2$) that the $(\gamma t)^*$ and $I^*$ are positively correlated, as curves with lower plateau infidelity appear to reach their plateau values at lower $(\gamma t)^*$ than those with higher plateau infidelities.

Studying these effects of the large-$\gamma t$ behaviour in more detail, in Fig. 3 we present simulations of a sample of 1 million Haar-random qubit ($d=2$) gates over $\gamma t \in [10^{-2}, 10^4]$. This larger sample size provides stronger evidence that the $I^*$ are indeed bounded above and below, as shown by the red and black horizontal lines, respectively. It was observed that the randomly generated quantum gates fall into one of two groups, indicated by the red and blue shaded curves.

<figure><img src="image_4.png" /></figure>

**Figure 3:** Large-$\gamma t$ behaviour of the AGIs for Haar-random qubit gates under pure dephasing. Simulations of the AGI were performed for a sample of one million (50 displayed) Haar-random qubit ($d=2$) gates under pure dephasing ($L=J_z$) for $\gamma t \in [10^{-2}, 10^4]$. For each gate, the AGI curves fall into one of two groups, exhibiting different behaviour near their respective saturation points ($(\gamma t)^*$). The curves in red ($I_{\text{mono}}$) approach their plateau values monotonically from below. The blue curves ($I_{\text{over}}$) rapidly approach their plateau values, before overshooting once with a single turning point, and then converging monotonically from above. The degree of blue shading for each curve indicates the degree of overshoot above the plateau value. All sampled gates converged to their $I^*$ plateau values within the range [$I_{\text{min}}^*, I_{\text{max}}^*$] indicated by the black and red horizontal lines. In the figure and insets, the lighter and darker shaded regions indicate the vertical extent of the sampled gates and their standard deviations, with the solid and dashed lines indicating the mean ($\tilde{I}^*$) and median ($\tilde{I}^*$) AGI values, respectively. The inset figures show each group of gates in isolation. The red group span the full bounding range, but are weighted upwards by outliers near the upper limit. The blue group do not span the full range, but are weighted towards the lower limit. Only the average AGI taken over both groups converges towards the expected mean value of 0.5.

The red shaded curves ($I_{\text{mono}}$) approach their plateau values monotonically from below. The blue shaded curves ($I_{\text{over}}$) all experience an overshoot above their plateau value, with a single turning point, before settling to $I^*$ monotonically from above. From the full sample, approximately 73% of all gates fell in the red, monotonic group, and the remaining 27% in the blue, overshoot group. The degree of blue shading represents the amount of overshoot, defined by the difference $I_{\text{max}} - I^*$ between the height of the peak and the plateau value. Furthermore, there appears to be some correlation between this quantity and the saturation point; the larger the overshoot the longer the curve takes to settle. Gates that experience this overshoot may therefore perform worse at intermediate noise strengths over the strong coupling limit.

The insets show each category plotted separately, from which it may be seen that, unlike the blue curves,
the red curves span the range between the upper and lower bounds, indicated by the lightly-shaded regions. The darker shaded regions, show the standard deviation about the mean, $\bar{\mathcal{I}}^*$, with the dashed line being the median value, $\bar{\mathcal{I}}^*$. We can see from these quantities that, despite spanning the range, the red curves are not distributed uniformly, but instead weighted upwards towards the upper bound, while the blue curves are weighted towards the lower bound. Both groups have significant numbers of outliers at either extreme, shifting the medians away from the means. In fact, it is only when averaging over the combined samples that the distribution becomes centered about the mean of $\bar{\mathcal{I}}^* = 0.5$ (see Appendix C.2 for further investigations on the distributions of these groups of curves within the bounding region).

## 2.4 Gate- and Noise-Dependence of the Fidelity

Having identified particular aspects of interest in the AGI behaviour in the strong coupling regime, let us investigate these characteristics in further detail.

Repeating the simulations presented in Fig. 3 for larger dimensions, we observed that the lower and upper bounds, $\mathcal{I}_{\min}^*$ and $\mathcal{I}_{\max}^*$ respectively, and the mean value, $\bar{\mathcal{I}}^*$, depend only on the dimension $d$ of the system, and are given precisely by the following expressions, for any $d$:

$$ \mathcal{I}_{\max}^* = 1 - \frac{1}{d+1}, \quad (41) $$

$$ \bar{\mathcal{I}}^* = 1 - \frac{1}{d}, \quad (42) $$

$$ \mathcal{I}_{\min}^* = 1 - \frac{2}{d+1}. \quad (43) $$

In Appendix C.2, we present these further examples of the distributions of $\mathcal{I}^*$ for ensembles of Haar random gates with dimensions $d = 2, 4, 8, 16$, which follow these expressions. It has previously been shown [41] that, for the case of a qudit in the singlet state evolving in a noisy quantum channel, the infidelity is bounded by $1 - \frac{2}{d+1} \le \mathcal{I}^* \le 1 - \frac{1}{d}$. This result matches with our case of a single qudit in the large-$d$ limit due to convergence of the mean and upper bound. However, for our single qudit system we shall use our existing framework to prove these limits for arbitrary $d$.

Writing the quantum channel as $\mathcal{E} = \mathcal{U} \circ \mathcal{M}$, and substituting into the expression for the AGF in Eq. (9),

$$ \mathcal{F} = \int_{\mathcal{H}} \mathrm{Tr}\{M[\rho_0]\rho_0\}d\rho_0 \quad (44) $$

$$ = \frac{\mathrm{Tr}\{M\} + \mathrm{Tr}\{M[1]\}}{d(d+1)}. \quad (45) $$

Given that the quantum channel is CPTP, and that $\mathcal{U}$ is a unitary operator, we know that the superoperator $M$ must also be trace-preserving and thus

<figure><img src="image_5.png" /><figcaption>Figure 4: Gate-dependence of the AGIs and qudit state evolution as functions of noise strength. (Data points) AGI simulations of a $d=4$ qudit under pure dephasing were time-evolved for each of the listed quantum gates: the (blue circles) QFT gate, and (other points) interpolated $X^n$ gates for $\eta \in [0, 1]$ as shown. Additionally shown for each gate: the (solid horizontal lines) plateau value $\mathcal{I}^*$ of the AGI, (dashed curves) state purity ($\mathrm{Tr}(\rho^2)$), and (dashed-dotted curves) mean coherences $\mathbb{E}[\rho_{ij}] = \frac{2}{d(d-1)} \sum_{i>j} |\rho_{i,j}|$ for $i > j$. Note that: (i) the AGI data points of the QFT and $X^n$ gate for $\eta = \frac{3d-2}{4d}$ overlap and converge to the same plateau value (for any $d$), (ii) the $\eta$ values were chosen as they produce equally spaced AGI plateau values (for any $d$), (iii) the final state purity and coherences for each gate were averaged over a set of 2400 random pure Hermitian initial states, and (iv) their curves converge to $\frac{1}{d-1}$ and 0 respectively, following the behaviour of the AGIs as they transition from the linear to nonlinear to plateau-like regimes.</figcaption></figure>

$\mathrm{Tr}\{M[1]\} = d$. Furthermore, in the strong coupling limit where the collapse operator $L = J_z$ describes pure dephasing, it is clear that the action of $M$ on a state $\rho$ is to transform it to a diagonal density matrix where the coherences (off-diagonal elements) have all decayed to zero. Expressing $M = \sum_{k=1}^K E_k^* \otimes E_k$ in the Kraus representation, the trace can be written as

$$ \mathrm{Tr}\{M\} = \sum_{k=1}^{K} \sum_{i=0}^{d-1} |E_{k, ii}|^2. \quad (46) $$

Furthermore, it can be shown that (see Appendix A.4 for a detailed derivation) each of the $d$ squared-elements are bounded by $0 \le |E_{k, ii}|^2 \le 1$ and hence in the large noise regime the AGI is bounded by

$$ 1 - \frac{2}{d+1} \le \mathcal{I}^* \le 1 - \frac{1}{d+1}. \quad (47) $$

Given the significance of this result, it is now of interest to consider how the behaviour is affected by the choice of quantum gate. In other words, we want to identify any particular unitary operators $U_g$, or properties thereof, whose AGIs $\mathcal{I}(\mathcal{E}, U_g)$ when $(\gamma t) \gg 1$ are precisely the $\mathcal{I}_{\min}^*$, $\mathcal{I}_{\max}^*$ or $\bar{\mathcal{I}}^*$.

In Fig. 4, we present simulations of the AGI, as a function of $\gamma t$, on a $d=4$ qudit for a set of gates that
do indeed reflect this behaviour when undergoing pure dephasing. The AGI values (discrete data points) are accompanied by dashed and dotted curves representing the state purity and mean coherences, while the horizontal lines show the $\mathcal{I}^*$ values of the associated gates. The purities and coherences were calculated over a set of $n = 2400$ randomly generated pure Hermitian density matrices. The purity $\text{Tr}\{\rho^2\}$ decays from 1 to a constant $\frac{1}{d-1}$, while the averages of the $\frac{d(d-1)}{2}$ coherences decay to zero, confirming the dephasing action of the $J_z$ collapse operator.

Regarding the set of gates, in addition to the QFT gate, we introduce the interpolated-$X$ gate (see Appendix B), $X^\eta$, where the matrix is raised to the $\eta$-th power for $\eta \in [0, 1]$. Clearly, $\eta = 1$ or 0 produces the $X$ or $\mathbb{1}_d$ gate, respectively, while for $0 < \eta < 1$ the action interpolates between these two extremes.

We find that

$$ \mathcal{I}^*(\mathcal{E}, \mathbb{1}_d) = \mathcal{I}_{\min}^*, \qquad (48) $$

$$ \mathcal{I}^*(\mathcal{E}, F_d) = \bar{\mathcal{I}}^*, \qquad (49) $$

$$ \mathcal{I}^*(\mathcal{E}, X_d) = \mathcal{I}_{\max}^*. \qquad (50) $$

Specifically, the $\mathcal{I}^*$ are extremised above and below by the $X$ and identity gates, respectively, while the mean value $\bar{\mathcal{I}}^*$ indeed corresponds to the action of the QFT gate that creates a superposition state. Furthermore, by careful choice of the values of $\eta$, we can obtain intermediate values of $\mathcal{I}^*$. For example, the $\eta$ values chosen and shown in Fig. 4 lead to curves with plateau values that are equally spaced within the bounding region. We observed that for the value $\eta = \frac{3d-2}{4d}$, the $\mathcal{I}^*$ matches that of the QFT gate, also converging to the mean.

The effect of these gates on the AGI can be understood by following the arguments in Appendix A.4 regarding the operator $M$ and its form based on the choice of quantum gate. It is clear that for the identity gate, or more generally any diagonal unitary matrix, $\text{Tr}\{M\} = d$ which minimises $\mathcal{I}^* = 1 - \frac{2}{d+1}$. Indeed, this result applies directly to the set of generalised-PHASE gates, including Clifford CLOCK gate $Z_d$ and the non-Clifford $\pi/8$ gate $T_d$, whose AGI curves $\mathcal{I}(\gamma t)$ follow that of the identity gate and plateau at $\mathcal{I}_{\min}^*$.

Similarly, the QFT gate produces a superposition state with each of the $d$ basis elements of $\rho$ being $\frac{1}{d}$, thus $\text{Tr}\{M\} = 1$ and $\mathcal{I}^* = 1 - \frac{1}{d}$. Finally the $X$ gate increments all states circularly, and, having no diagonal elements, results in traceless $\text{Tr}\{M\} = 0$ and thus $\mathcal{I}^* = 1 - \frac{1}{d+1}$.

Let us now consider how the choice of quantum gate affects the saturation point $(\gamma t)^*$ at which the AGI attains the plateau value $\mathcal{I}^*$. Figure 5 presents simulations of the plateau saturation points as a function of the qudit dimension for set of $X^\eta$ gates over a range of $\eta$ values from 0 to 1. In addition, the calculations were repeated over a set of $n = 4800$ Haar random gates at each dimension from $d = 2$ to 12.

<figure><img src="image_6.png" /><figcaption>Figure 5: Saturation points $((\gamma t)^*)$ of the AGIs for interpolated $X^\eta$ gates as functions of qudit dimension. AGIs were calculated for a set of 12 interpolated $X^\eta$ gates for $0 \le \eta \le 1$, at dimensions $d \in [2, 12]$. For each dimension and gate, the AGI curve over $\gamma t$ was interpolated by a cubic spline. The saturation points $(\gamma t)^*$ were identified by root-finding algorithm of the points at which they converged (within $\varepsilon < 10^{-8}$) to their AGI plateau values. The plotted data points, on log-linear axes, indicate these root values as functions of $d$ for each gate, and the dashed lines represent their associated power-law fits. This analysis was repeated for a set of $n = 4800$ Haar-random gates at each dimension. The light shaded areas are bounded by the maximum and minimum saturation points found at each dimension. The darker shaded areas represent a $1\sigma$ deviation about the mean (solid gray line).</figcaption></figure>

The grey curve represents the mean value of the saturation point at each dimension, while the dark grey and light gray shaded regions indicate the standard deviations and max-min values, respectively.

For each of the simulated $X^\eta$ gates, the dashed-dotted lines represent power-law fits to the discrete data points at each dimension. The power-law model is given by

$$ (\gamma t)^* = \alpha(d - \beta)^\delta, \qquad (51) $$

and, fitted by nonlinear least-squares regression, resulted in $R^2 > 0.999$ for each fit. It is interesting to note that for the case of the identity gate, the saturation point remains constant independent of dimension.

Continuing the investigation, we performed regression analysis of the fitted model parameters $\alpha$, $\beta$ and $\delta$ across the set of $\eta$ values. The fits are presented in Appendix C.3. Parameter $\alpha$ was best described by a sigmoid model while $\beta$ and $\delta$ by exponential models,

$$ \alpha(\eta) = \frac{\alpha_0}{(1 + e^{\alpha_1(\eta - \alpha_2)})} + \alpha_3, \qquad (52) $$

$$ \beta(\eta) = \beta_0 e^{\beta_1(\eta - \beta_2)} + \beta_3, \qquad (53) $$

$$ \delta(\eta) = \delta_0 e^{\delta_1(\eta - \delta_2)} + \delta_3. \qquad (54) $$

Thus far in our investigation, we have restricted analysis to the effect of pure dephasing by the operator $J_z$ on a qudit. Let us conclude this section by
<figure><img src="image_7.png" /><figcaption>Figure 6: Effect of bit-flip errors and relaxation processes on the AGI. X^η and QFT gate simulations were repeated on qudits for d = 4 undergoing either (circles) bit-flip (L = J_x) errors or (squares) spin relaxation (L = J_-), rather than pure dephasing noise. The following gates were simulated: 1 = X^0 (black, η = 0), X^η (green, η = 3d/2), X (red, η = 1) and QFT (blue). The solid horizontal lines show the upper (1 - 1/d, red) and lower (1 - 2/d, black) bounds of the AGIs, as well as the mean (1 - 1/d, blue). For L = J_-, the AGIs for all matrix powers of the X^η gate behaved almost exactly as the QFT gate, while for L = J_x the plateau AGI values for the non-zero matrix powers were observed to deviate upwards from the identity gate, but not to the same extent as for the case of pure dephasing. On the other hand, the QFT gate response was the same for all types of noise.</figcaption></figure>

studying the effects of the $L = J_x$ and $L = J_-$ operators, corresponding to bit-flip errors and relaxation processes. In Fig. 6, we present simulations of 4 gates: the QFT gate F and 3 X^η gates for $η = 0, 3d/2, 1$, for $d = 4$ qudits interacting with the environment via each of the collapse operators separately. We may make the following remarks based on these results: (1) the plateau values $\mathcal{I}^*$ for the gates remain bounded between $\mathcal{I}_{\min}^*$ and $\mathcal{I}_{\max}^*$, (2) the QFT gate plateaus at the same mean value $\mathcal{I}^*$ for all $L \in \{J_z, J_x, J_\_\}$, (3) for bit-flip errors the behaviour of the identity gate remains the same, while the X^η gates plateau below the QFT gate, and (4) for relaxation processes all gates have nearly the same behaviour as the QFT gate, and all plateau at $\mathcal{I}^*$.

## 2.5 First- and Second-Order Correction Terms to the AGI

Section 2.2 presented our general result for the perturbative expansion of the AGI to arbitrary order m in Eq. (33). Now, let us investigate in further detail the first and second-order correction terms.

Beginning with the first-order correction, at $m = 1$, it can be seen that (see Appendix A.5) the trace over

$M^{(1)}$ simplifies to

$$ \mathcal{I}^{(1)} = - \frac{(\gamma t)}{d(d+1)} \mathrm{Tr}\{\mathcal{L}\}, \qquad (55) $$

which confirms the gate-independent nature of the first-order term as the AGI depends only on the noise superoperator $\mathcal{L}$. Additionally, this trace can be expressed in terms of the regular operator L as

$$ \mathrm{Tr}\{\mathcal{L}\} = |\mathrm{Tr}\{L\}|^2 - d \mathrm{Tr}\{L^\dagger L\}, \qquad (56) $$

reproducing the results of [21]. For the case of pure dephasing, we obtain the familiar result in Eq. (40)

$$ \mathcal{I}^{(1)} = \frac{\gamma t}{12} d(d-1), \qquad (57) $$

using the element-wise definition of the dephasing operator in Eq. (22)

$$ (J_z)_{ii} = \frac{d - (2i - 1)}{2}. \qquad (58) $$

Now, if we compare this analytical result with numerical simulations, we observe in Fig. 7 the relative error of the first-order correction to the AGI for dimensions $d = 2, 4$ over $\gamma t \in [0, 0.1]$. This relative error is expressed by

$$ \varepsilon_{\mathcal{I}}^{(1)} = \frac{\mathcal{I} - \mathcal{I}^{(1)}}{\mathcal{I}}. \qquad (59) $$

The simulations were performed for the QFT, identity, X and X^η gates for $η = 3d/2$. The dashed lines show the expected analytic results for each dimension. The dark and light shaded regions represent the standard deviations and max-min values taken over a sample of Haar random gates

Clearly, the gate independence of the first-order correction term matches precisely the identity gate. The remaining gates show the expected deviation from linearity, with the analytic value underestimating the AGI values. We observe that the relative error for $d = 2$ is approximately one order of magnitude smaller than for $d = 4$.

Moving on to the second-order correction term, we have from Eq. (33)

$$ \mathcal{I}^{(2)} = \frac{-(\gamma t)^2}{d(d+1)} \times \sum_{n_1, n_2=0}^{\infty} \frac{(-t)^{n_1+n_2} \mathrm{Tr}\{[(S)^{n_1}, \mathcal{L}][(S)^{n_2}, \mathcal{L}]\}}{n_1!(n_2+1)!(n_1+n_2+2)}, \qquad (60) $$

which we reduce (see Appendix A.6) to a single iterated commutator using the relation $n_1 + n_2 + 2 = s$,

$$ \mathcal{I}^{(2)} = - \frac{(\gamma t)^2}{d(d+1)} \sum_{s=0}^{\infty} \mathrm{Tr}\{\mathcal{L}[(S)^s, \mathcal{L}]\} \frac{(-t)^s}{(s+2)!}. \qquad (61) $$
<figure><img src="image_8.png" /><figcaption>Figure 7: Relative error of the first-order correction to the AGI, and gate-dependent deviation from linearity. AGIs were simulated for 4 gates acting on qubits (circles) and $d = 4$ qubits (squares) undergoing $L = J_z$ pure dephasing. The relative error $\epsilon_I^{(1)} = \frac{I-I^{(1)}}{I}$ was calculated for the first-order analytical correction. The dashed-dotted and dashed lines represent the first-order term of the perturbative expansion of the AGI, for $d = 2$ and $d = 4$, respectively. These gate-independent terms correspond to the behaviour of the identity gate, $|1\rangle = X^0$ (black, $\eta = 0$). The $X^\eta$ (green, $\eta = \frac{3d-2}{4d}$), $X^1$ (red, $\eta = 1$) and $F$ (blue, QFT) gates show significant deviation from the gate-independent term, and for $d = 4$ are roughly 10 times larger than for $d = 2$. The solid and dashed grey lines represent the mean and median relative errors over a sample of $n = 4800$ Haar-random gates, with the light and dark shaded areas indicating the min-max bounds and standard deviations.</figcaption></figure>

For $s = 0$, the trace term reduces further to $\text{Tr}\{\mathcal{L}^2\}$ by definition of the iterated commutator. Additionally, we show (see Appendix A.7) that, using the property [42],

$$[(A)^s, B] = \sum_{k=0}^{s} (-1)^k \binom{s}{k} A^{s-k} BA^k, \quad (62)$$

the trace over the iterated commutator is zero for all odd values of $s$. Therefore,

$$\mathcal{I}^{(2)} = -\frac{(\gamma t)^2}{d(d+1)} \times \left( \frac{\text{Tr}\{\mathcal{L}^2\}}{2} + \sum_{s=2}^{\infty} \frac{\text{Tr}\{\mathcal{L}[(S)^s, \mathcal{L}](-t)^s\}}{(s+2)!} \right). \quad (63)$$

This shows us that the AGI remains gate-independent at order $\gamma^2 t^2$. Indeed the first appearance of gate-dependence through the $S$ operator occurs at order $\gamma^2 t^4$, when $s = 2$ and the iterated commutator term contains $S$ and is non-zero.

Since the summation of $s$ is over all even integers, we can introduce a maximum cutoff value to render this calculation computationally tractable. Indeed, for the correction term to not be unbound, we must have that the sum is convergent, allowing us to iden-

<figure><img src="image_9.png" /><figcaption>Figure 8: Iteration-dependent convergence of the second-order AGI correction against qudit dimension. The critical order $s_\epsilon$ was calculated as the number of iterations over the summation term in the second-order AGI correction such that the absolute difference $|f_2(s_\epsilon) - f_2(s_\epsilon - 1)| < \epsilon$ is less than the error threshold $\epsilon = 1 \times 10^{-8}$. Calculations were performed for a set of interpolated $X^\eta$ gates and the QFT gate, as well as $n = 2500$ Haar-random gates, over dimensions $d \in [2, 64]$. The identity gate ($\eta = 0$) only required the minimal $s = 2$ iterations, independent of dimension, since the identity commutes with the collapse operator. Of the interpolated gates, the number of iterations required scales with the interpolation parameter $\eta$, while the QFT gate shows marked variation over $d$, without exceeding the critical order for $\eta = 1$. The mean of the Haar-random gates (dark grey) also lies below the $\eta = 1$ curve, while the min-max interval, (light grey), shows gates at lower dimension requiring a larger number of iterations to converge.</figcaption></figure>

tify such a value. We can understand this computationally by identifying the $s$-th order of the summation,

$$f(s) = \text{Tr}\{\mathcal{L}[(S)^s, \mathcal{L}]\} \frac{(-t)^s}{(s+2)!}. \quad (64)$$

Then, for a given error threshold $\epsilon$, the cutoff order $s_\epsilon$ is simply the $s$ at which the absolute difference from the preceding order $s_\epsilon - 2$ is below this threshold,

$$|f(s_\epsilon) - f(s_\epsilon - 2)| < \epsilon. \quad (65)$$

In Fig. 8 we present this iteration-dependent convergence of the second-order AGI correction term. Our set of interpolated $X^\eta$ gates and the QFT gate were simulated over dimensions $d \in [2, 64]$ for a uniform gate time of $t = 1$. We calculated for each of them the cutoff order $s_\epsilon$ within an error threshold $\epsilon = 1 \times 10^{-8}$. This analysis was repeated for a set of $n = 2500$ Haar random gates, and shown by the grey curve (mean value) and shaded region (bounds). It is useful to note that the number of iterations required grows approximately logarithmically with $d$. Furthermore we have a fairly clear range of validity on $s_\epsilon$ for the gates,
<figure><img src="image_10.png" /><figcaption>Figure 9: Time derivative of the convergence of the second-order AGI correction over qudit dimension. Linear gradients of the critical order $s_ε$ as a function of $t$ were fitted by linear regression for each of a set of gates at each dimension $d \in [2, 48]$. The gate set included the QFT, identity ($X^0$), $X$, and uniformly interpolated $X^\eta$ gates and $n = 200$ Haar-random gates. For the identity gate this is constant at zero-gradient, while the gradients for the interpolated gates show rapid initial growth that levels off as $d$ increases. The gradients of the QFT gates for $d > 10$ appear to populate two distinct regimes around 25 and 33. The sampled Haar-random gates exhibit behaviour similar to that of the $X$ gate within one standard deviation of the mean, and the large variation at lower dimensions stabilises as $d$ increases.</figcaption></figure>

although the $X$ gate does not give a consistent upper bound here for smaller dimensions. On the other hand, the identity gate $1 = X^0$ always gives the lower bound of $s_ε = 2$, since then $[1, \mathcal{L}] = 0$.

Next, we need to investigate the effect of the gate time parameter $t$ on the second-order correction term, since in the summation it is raised to the power of $s$. We choose our set of quantum gates, $U_g \in \{X^\eta, F\}$ with $\eta \in \{0, \frac{3d-2}{8d}, \frac{9d-6}{16d}, \frac{3d-2}{4d}, 1\}$, and simulate them on qubits of dimension $d \in [2, 48]$. Instead of varying the noise parameter $\gamma$ with constant gate time $t = 1$, we fix $\gamma = 1$ and then vary the gate time $t \in [0, 5]$. Note that, in order to implement the full quantum gate during this modified gate time, the amplitude of the control-pulse Hamiltonian matrix must be modulated by the inverse of the gate time, $H_c \to \frac{1}{t}H_c$. We observed (see Appendix C.4 for an example of $s_ε$ against $t$ for the $X$ gate) that for each gate at each dimension, the number of iterations required to reach convergence was roughly linear in time: $s_ε \propto t$. Hence to each such set of data we fitted using linear regression the model

$$s_{\varepsilon} = mt + c, \qquad (66)$$

where $m = \frac{ds_\varepsilon}{dt}$. Figure 9 shows this curve of fitted gradients for each gate as a function of the qudit dimension. This analysis is also extended to a

<figure><img src="image_11.png" /><figcaption>Figure 10: Relative error of the AGI up to second-order. AGIs were simulated for 4 fixed gates and $n = 1000$ Haar-random gates acting on qubits (circles) and $d = 4$ qudits (squares) undergoing $L = J_z$ pure dephasing. The relative error $\varepsilon_I^{(2)} = \frac{\mathcal{I} - \mathcal{I}^{(1)} - \mathcal{I}^{(2)}}{\mathcal{I}}$ was calculated for the second-order perturbative correction. The error for $d = 4$ is roughly 1% while for qubits this is reduced to 0.01%, two orders of magnitude lower than the relative error at first-order.</figcaption></figure>

set of $n = 200$ Haar random gates at each dimension, indicated by the grey curve and shaded regions. The interpolated gates all appear to evolve smoothly to a nearly constant value at large $d$. On the other hand, the QFT gate for $d > 6$ appears to oscillate randomly between two constant regimes near $\frac{ds_\varepsilon}{dt} \approx 25$ and $\frac{ds_\varepsilon}{dt} \approx 33$.

Based on these two analyses, we can identify a safe bound on $s_ε$ for precise calculation of the second-order correction. Recall our standard approach of normalising the Hamiltonian such that the gate time $t = 1$, varying $\gamma$ only, and generally considering $d < 64$. Then, we observe that approximately 50 iterations are sufficient to calculate the second-order correction to a precision of $\varepsilon < 1 \times 10^{-8}$.

We now present in Fig. 10 our results for the relative error of the AGI up to second-order. The results presented follow the same approach as for those up to first-order shown in Fig. 7. Choosing the same set of gates, $\{X^0, X^\eta, X^1, F\}$ with $\eta = \frac{3d-2}{4d}$ and evolving over $\gamma t \in [0, 0.1]$, we show the scaling of the relative error for $d = 2, 4$,

$$\varepsilon_I^{(2)} = \frac{\mathcal{I} - \mathcal{I}^{(1)} - \mathcal{I}^{(2)}}{\mathcal{I}}. \qquad (67)$$

Comparing these results with those of just the first-order in Fig. 7, we see that the error for $d = 4$ is reduced by one order of magnitude to $\sim 1\%$, while for $d = 2$ the error of $\sim 0.01\%$ is lower by two orders of magnitude. Additionally, the sign of the relative error has changed from negative to positive, indicating that despite the increased precision, the correction up to second order will tend to understate the actual AGI
value, with this error growing quadratically in $\gamma t$.

Finally, as for the first-order case, it is interesting to express the trace of the repeated commutator of superoperators $S$ and $\mathcal{L}$ in terms of regular operators $H$ and $L$. This will give a clear representation of the effect of the control Hamiltonian on the AGI, as well as its interaction with the collapse operator. Let us consider separately the two trace terms of $\mathcal{I}^{(2)}$ in Eq. (63): $\mathrm{Tr}\{\mathcal{L}^2\}$ and $\mathrm{Tr}\{\mathcal{L}[(S)^s, \mathcal{L}]\}$. We find (see Appendix A.5) that the gate-independent expression can be written as

$$ \mathrm{Tr}\{\mathcal{L}^2\} = |\mathrm{Tr}\{L^2\}|^2 + \frac{1}{2} \mathrm{Tr}\{L^\dagger L\}^2 + \frac{d}{2} \mathrm{Tr}\{(L^\dagger L)^2\}, \quad (68) $$

which, in the case of pure dephasing, $L = J_z$, reduces to

$$ \mathrm{Tr}\{J_z^2\} = \frac{d^2(3 - 5d^2 + 2d^4)}{120}, \quad (69) $$

using the element-wise definition in Eq. (58). This allows us express the $\gamma^2 t^2$ correction in terms of the first-order correction term from Eq. (57),

$$ -\frac{(\gamma t)^2}{d(d+1)} \frac{\mathrm{Tr}\{J_z^2\}}{2} = -(\gamma t) \frac{2d^2 - 3}{20} \mathcal{I}^{(1)}. \quad (70) $$

However, it is not necessarily the case that any $m$-th order expression can be written recursively in terms of lower orders. Equation (70) for the second-order appears to be a special case.

For the iterated commutator, we find that (see Appendix A.6)

$$ \begin{aligned} \mathrm{Tr}\{\mathcal{L}[(S)^s, \mathcal{L}]\} & = \sum_{k=0}^{s} \sum_{j=0}^{k} \sum_{l=0}^{s-k} \binom{s}{k} \binom{k}{j} \binom{s-k}{l} (-1)^{s+k-j-l} \\ & \times \mathrm{Re} \left[ \mathrm{Tr}\{LH^{k-j}LH^{s-k-l}\}^* \mathrm{Tr}\{LH^jLH^l\} \right. \\ & \qquad \left. + \frac{1}{2} \mathrm{Tr}\{L^T L^* H^{s-j-l}\}^* \mathrm{Tr}\{L^T L^* H^{j+l}\} \right. \\ & \qquad \left. + \frac{1}{2} \mathrm{Tr}\{L^T L^* H^{k-j} L^T L^* H^{s-k-l}\}^* \mathrm{Tr}\{H^{j+l}\} \right. \\ & \qquad \left. - \mathrm{Tr}\{LH^{k-j}L^T L^* H^{s-k-l}\}^* \mathrm{Tr}\{LH^{j+l}\} \right. \\ & \qquad \left. - \mathrm{Tr}\{L^T L^* H^{k-j}LH^{s-k-l}\}^* \mathrm{Tr}\{LH^{j+l}\} \right]. \end{aligned} \quad (71) $$

This gives a complete description of the second-order gate-dependence in terms of the explicit Hamiltonian terms, without resorting to the superoperator representation. The superoperator and operator expressions have complexity $O sd^6)$ and $O((sd)^3)$, respectively. Since we observed in Fig. 8 that the maximal $s_\epsilon$ needed to reach convergence is roughly logarithmic in $d$, the complexities can then be reduced to $O(d^6 \log d)$ and $O((d \log d)^3)$. Therefore, for $d \gg 2$,

it is computationally preferred to use the operator-based expression, while for smaller qudits with $d \sim 2$, the superoperator representation is generally preferable for numerical simulations.

### 3 Conclusions and Outlook

In this work, we have performed a comprehensive analysis of the AGF for single qudit open quantum systems coupled to Markovian noise environments in the Lindblad superoperator formalism. Our primary contributions and findings are as follows:

We have expanded the AGF perturbatively in powers of the dimensionless coupling constant $\gamma t$. From this we derived our main result in Eq. (33) for the general expressions of the correction terms to arbitrary order. The result was expressed in terms of the iterated commutators of the gate and noise superoperators, $S$ and $\mathcal{L}$, respectively, based on the expansion obtained in Eq. (30). It is particularly significant that these corrections were expressed in a form that may be directly implemented in numerical calculations. We found that, at each order above the first, the correction term could be separated into a gate-independent term $O((\gamma t)^m)$ depending only on $\mathcal{L}$, and gate-dependent terms containing $S$ and higher powers of $t$. The first-order correction term, $\mathcal{I}^{(1)}$, is gate-independent and depends solely on the noise superoperator $\mathcal{L}$. The gate-dependence appears for the first time at $\gamma^2 t^4$, in the second-order correction term $\mathcal{I}^{(2)}$, highlighting the significant role of the interaction between the control Hamiltonian $S$ and the noise operator $\mathcal{L}$. Explicit expressions for the first- and second-order correction terms in the operator representation were also derived.

Our numerical simulations under pure dephasing noise revealed a clear transition from linear to nonlinear behaviour in the AGI as the noise coupling strength $\gamma t$ increases. For $\gamma t \sim 1$, the AGI deviates significantly from the linear approximation and eventually reaches a plateau at a stable value. The plateau values, $\mathcal{I}^*$, as well as their corresponding saturation points, $(\gamma t)^*$, were found to be dependent not only on the qudit dimension $d$ but also on the specific quantum gate implemented.

Indeed, utilising our theoretical framework we were able to derive in Eq. (47) a significant result for universal upper and lower bounds of the AGI plateau values in the strong coupling regime, depending only on the dimensionality of the system. Furthermore we were able to identify specific gates that saturate these bounds: The identity gate ($1_d$) and the generalised-NOT ($X$) achieve the lower and upper limits, respectively, while the QFT gate ($F_d$) was found to plateau at the mean AGI value of the Haar measure. Additionally, we introduced the interpolated $X^\eta$ gates to demonstrate the range of AGI values between these bounds. The saturation points of these
interpolated gates were studied numerically, and were
well-described by a power-law model of $(\gamma t)^*$ vs $d$.

We analyzed the convergence behavior of the
second-order correction term, identifying that approx-
imately 50 iterations are sufficient for precise calcu-
lation up to a qudit dimension of $d < 64$ within an
absolute tolerance of $\varepsilon < 1 \times 10^{-8}$. The number of it-
erations to convergence was found to scale logarithmi-
cally with the system dimension $d$, and linearly with
the gate time $t$.

Numerical simulations confirmed the analytical pre-
dictions, demonstrating that the relative error of
the AGI is significantly reduced when including the
second-order correction term. For $d = 4$, the error is
reduced to the order of 1%, and for $d = 2$, it is reduced
to 0.01%, indicating the necessity of higher-order cor-
rections for accurate fidelity modeling, particularly as
the dimension of the system increases.

The detailed insights gained from this study have
several important applications and implications for
the field of quantum computing, particularly given
the growing interest in qudits. Understanding the be-
havior of the AGI and its dependence on the type
of quantum gate and environmental noise facilitates
the optimization of quantum gate design. Our ap-
proach provides a methodology for identifying quan-
tum gates that have favourable fidelity characteristics.
Thus, through appropriate choice or design of basis
gates it may be possible to enhance the performance
of quantum circuits. This could lead to significant im-
provements for optimising the robustness of quantum
algorithms on near-term noisy platforms.

The results underscore the importance of incorpo-
rating higher-order corrections into error correction
protocols. This is particularly relevant for qudits of
higher dimensions, where noise effects are more pro-
nounced. We envisage that the detailed perturbative
expansions provided in this work could be leveraged
to develop advanced error correction techniques, for
example by incorporating the higher-order correction
terms into the cost function of optimal control meth-
ods. In particular, our correction terms can be used
to quantify and mitigate the error in logical qubit em-
bedding protocols up to arbitrary order. Future re-
search is needed to show how this could result in pro-
tection against such errors above first-order, and how
this might incorporate the gate-dependencies identi-
fied.

The derived bounds and behaviours of AGI for
different gates and dimensions also serve as bench-
marks for assessing the performance of quantum sys-
tems. These benchmarks can guide experimental-
ists in evaluating and improving their quantum hard-
ware and gate implementations. Indeed further study
is needed to extend these results from the idealised
single-pulse (time-independent) Hamiltonian to real-
istic pulse-based optimal control techniques for gate
generation.

While this study focused primarily on pure dephas-
ing and bit-flip errors, future research should explore
further the impact of other noise models, such as am-
plitude damping and depolarizing noise. Understand-
ing how different types of noise affect the AGI will
provide a more comprehensive framework for design-
ing noise-resilient quantum systems.

Extending the framework to multi-qudit systems is
a necessity. These introduce additional complexity
due to inter-qudit interactions and correlated noise
effects. Investigating these factors will be crucial for
the scalability of qudit-based architectures, but could
prove fruitful in reducing the circuit complexity of
long range entangling operations.

Finally, it would be beneficial to have experimen-
tal verification of these predictions. Testing the de-
rived correction terms and transition points in real
quantum systems will help to validate and refine the
models presented.

In conclusion, this study provides a detailed the-
oretical foundation for understanding the fidelity of
quantum gates in noisy environments. The findings
contribute to the ongoing efforts to develop robust,
high-fidelity quantum operations, paving the way for
practical and scalable quantum computing.

## Acknowledgments

This work was funded by the French National
Research Agency (ANR) through the Programme
d'Investissement d'Avenir under contract ANR-11-
LABX-0058_NIE and ANR-17-EURE-0024 within
the Investissement d'Avenir program ANR-10-IDEX-
0002-02. D.J. and M.R. gratefully acknowledge fi-
nancial support from the Deutsche Forschungsgemein-
schaft (DFG, German Research Foundation) through
the Collaborative Research Centre "4f for Future"
(CRC 1573, project number 471424360) project B3.
J-G.H. also acknowledges QUSTEC funding from the
European Union's Horizon 2020 research and inno-
vation program under the Marie Skłodowska-Curie
Grant Agreement No. 847471. The authors would
like to acknowledge the High Performance Comput-
ing Center of the University of Strasbourg for sup-
porting this work by providing scientific support and
access to computing resources. Part of the computing
resources were funded by the Equipex Equip@Meso
project (Programme Investissements d'Avenir) and
the CPER Alsac calcul/Big Data.

## Code and Data Availability

All code and data utilised in the simulations, analyses,
and generation of the results presented in this study
are publicly available via Zenodo [43] under a Creative
Commons Attribution (CC BY) license.
# A Complementary Derivations

Here, we present the detailed derivations of the analytical results in the main text.

## A.1 Integration of the m-th Order Solution to the Master Equation

**Theorem A.1.** Consider an open quantum system of a single qudit of dimension $d$, initiated in a pure state $\text{Tr}\{\rho_0^2\} = 1$, under the influence of a time-independent noise superoperator $\mathcal{L}$ with coupling constant $\gamma$, and evolving via a time-independent unitary superoperator $\mathcal{S}$ over time $t$. The solution to the master equation in Eq. (1) is given by the quantum channel in Eq. (10), and can be expressed in a perturbative expansion in $\gamma$ by Eq. (26), where the m-th order correction term is given by Eq. (27). Then, the nested integral can be evaluated, by induction, and results in Eq. (29):

$$ \tilde{M}^{(m)}(t) = \int_{0}^{t} \cdots \int_{0}^{t_{m-1}} \left( \prod_{i=1}^{m} e^{-S t_i} \mathcal{L} e^{S t_i} \right) dt_m \cdots dt_1 \quad (72) $$

$$ = t^m \sum_{n_1, \dots, n_m = 0}^{\infty} \left( \prod_{i=1}^{m} \frac{(-t)^{n_i} [(\mathcal{S})^{n_i}, \mathcal{L}]}{n_i! \sum_{j=i}^{m} (n_j + 1)} \right). \quad (73) $$

*Proof.* In the calculations that follow, we make use of Campbell's lemma, based on the Baker-Campbell-Hausdorff Formula [37, 38], as stated in Eq. (28),

$$ e^X Y e^{-X} = \sum_{n=0}^{\infty} \frac{[(X)^n, Y]}{n!}, \quad (74) $$

with $[(X)^n, Y] = [X, [(X)^{n-1}, Y]]$ and $[(X)^0, Y] = Y$, and proceed with proof by induction.

*Base Case:* To first order, $m = 1$,

$$ \int_{0}^{t} e^{-S t_1} \mathcal{L} e^{S t_1} dt_1 = \int_{0}^{t} \sum_{n_1=0}^{\infty} \frac{[(-S t_1)^{n_1}, \mathcal{L}]}{n_1!} dt_1 \quad (75) $$

$$ = \sum_{n_1=0}^{\infty} \frac{(-1)^{n_1}}{n_1!} [\mathcal{S})^{n_1}, \mathcal{L}] \int_{0}^{t} t_1^{n_1} dt_1 \quad (76) $$

$$ = t \sum_{n_1=0}^{\infty} \frac{(-t)^{n_1} [(\mathcal{S})^{n_1}, \mathcal{L}]}{n_1!(n_1+1)} \quad (77) $$

$$ = \tilde{M}^{(1)}(t). \quad (78) $$

*Induction Step:* Assume that the result in Eq. (73) holds at m-th order, and relabel indices $t, t_1, \dots, t_{m-1} \to t_1, t_2, \dots, t_m$ and $t \to t_1$:

$$ \tilde{M}^{(m)}(t_1) = \int_{0}^{t_1} \cdots \int_{0}^{t_m} \left( \prod_{2=1}^{m+1} e^{-S t_i} \mathcal{L} e^{S t_i} \right) dt_{m+1} \cdots dt_2 \quad (79) $$

$$ = t_1^m \sum_{n_2=0}^{\infty} \cdots \sum_{n_{m+1}=0}^{\infty} \left( \prod_{i=2}^{m+1} \frac{(-t_1)^{n_i} [(\mathcal{S})^{n_i}, \mathcal{L}]}{n_i! \sum_{j=i}^{m+1} (n_j+1)} \right). \quad (80) $$

Then, to order $m+1$:

$$ \int_{0}^{t} \tilde{M}^{(m)}(t_1) (e^{-S t_1} \mathcal{L} e^{S t_1}) dt_1 = \int_{0}^{t} \left( t_1^m \sum_{n_2=0}^{\infty} \cdots \sum_{n_{m+1}=0}^{\infty} \prod_{i=2}^{m+1} \frac{(-t_1)^{n_i} [(\mathcal{S})^{n_i}, \mathcal{L}]}{n_i! \sum_{j=i}^{m+1} (n_j+1)} \right) \left( \sum_{n_1=0}^{\infty} \frac{(-t_1)^{n_1} [(\mathcal{S})^{n_1}, \mathcal{L}]}{n_1!} \right) dt_1 \quad (81) $$

$$ = \sum_{n_2=0}^{\infty} \cdots \sum_{n_{m+1}=0}^{\infty} \frac{(-1)^{n_1} [(\mathcal{S})^{n_1}, \mathcal{L}]}{n_1!} \prod_{i=2}^{m+1} \left( \frac{(-1)^{n_i} [(\mathcal{S})^{n_i}, \mathcal{L}]}{n_i! \sum_{j=i}^{m+1} (n_j+1)} \right) \int_{0}^{t} t_1^{m} t_1^{n_1} \prod_{i=2}^{m+1} t_1^{n_i} dt_1 \quad (82) $$

$$ = \sum_{n_2=0}^{\infty} \cdots \sum_{n_{m+1}=0}^{\infty} \frac{(-1)^{n_1} [(\mathcal{S})^{n_1}, \mathcal{L}]}{n_1!} \prod_{i=2}^{m+1} \left( \frac{(-1)^{n_i} [(\mathcal{S})^{n_i}, \mathcal{L}]}{n_i! \sum_{j=i}^{m+1} (n_j+1)} \right) \frac{t^{m+1+\sum_{i=1}^{m+1} n_i}}{m+1+\sum_{j=1}^{m+1} n_j} \quad (83) $$

$$ = t^{m+1} \sum_{n_2=0}^{\infty} \cdots \sum_{n_{m+1}=0}^{\infty} \frac{(-t)^{n_1} [(\mathcal{S})^{n_1}, \mathcal{L}]}{n_1! \sum_{j=1}^{m+1} (n_j+1)} \prod_{i=2}^{m+1} \frac{(-t)^{n_i} [(\mathcal{S})^{n_i}, \mathcal{L}]}{n_i! \sum_{j=i}^{m+1} (n_j+1)} \quad (84) $$

$$ = t^{m+1} \sum_{n_2=0}^{\infty} \cdots \sum_{n_{m+1}=0}^{\infty} \prod_{i=1}^{m+1} \frac{(-t)^{n_i} [(\mathcal{S})^{n_i}, \mathcal{L}]}{n_i! \sum_{j=i}^{m+1} (n_j+1)} \quad (85) $$

$$ = \tilde{M}^{(m+1)}(t). \quad (86) $$
□ Weingarten calculus. Specifically, we make use of the following result,

## A.2 Perturbative Expansion of the AGI

We define the Average Gate Infidelity as

$$ \mathcal{I}(\mathcal{E}, \mathcal{U}) = 1 - \mathcal{F}(\mathcal{E}, \mathcal{U}) \quad (87) $$

$$ = 1 - \left( 1 - \sum_{m=1}^{\infty} \gamma^m F^{(m)} \right), \quad (88) $$

where the quantum channel is given by

$$ \mathcal{E} = \mathcal{U} + \sum_{m=1}^{\infty} \gamma^m \mathcal{E}^{(m)} \quad (89) $$

$$ = \mathcal{U} + \sum_{m=1}^{\infty} \gamma^m \mathcal{U} \tilde{M}^{(m)}. \quad (90) $$

Then, from Eq. (9),

$$ \mathcal{F}(\mathcal{E}, \mathcal{U}) = \int_{\mathcal{H}} \langle (\mathcal{U}^{\dagger} \circ \mathcal{E}) [\rho_0] \rangle_0 d\rho_0 \quad (91) $$

$$ = \int_{\mathcal{H}} \left\langle \left( \mathcal{U}^{\dagger} \circ \left( \mathcal{U} + \sum_{m=1}^{\infty} \gamma^m \mathcal{U} \tilde{M}^{(m)} \right) \right) [\rho_0] \right\rangle_0 d\rho_0 \quad (92) $$

$$ = \int_{\mathcal{H}} \left\langle \rho_0 + \sum_{m=1}^{\infty} \gamma^m \tilde{M}^{(m)} [\rho_0] \right\rangle_0 d\rho_0, \quad (93) $$

and evaluating the expectation value over $\rho_0$ of initial states with $\langle A \rangle_0 = \text{Tr}\{A\rho_0\}$,

$$ \begin{aligned} & 1 - \sum_{m=1}^{\infty} \gamma^m F^{(m)} \\ & = \int_{\mathcal{H}} \text{Tr} \left\{ \rho_0^2 + \sum_{m=1}^{\infty} \gamma^m \tilde{M}^{(m)} [\rho_0] \rho_0 \right\} d\rho_0 \end{aligned} \quad (94) $$

$$ = 1 + \sum_{m=1}^{\infty} \gamma^m \int_{\mathcal{H}} \text{Tr} \{ \tilde{M}^{(m)} [\rho_0] \rho_0 \} d\rho_0, \quad (95) $$

since $\text{Tr}\{\rho_0^2\} = 1$ and $\int_{\mathcal{H}} d\rho_0 = 1$, and using Eq. (25). The m-th order term of the AGF is therefore

$$ F^{(m)} = - \int_{\mathcal{H}} \text{Tr} \{ \tilde{M}^{(m)} [\rho_0] \rho_0 \} d\rho_0, \quad (96) $$

and we may write the AGI in terms of $t^m M^{(m)} = \tilde{M}^{(m)}$,

$$ \mathcal{I}(\mathcal{E}, \mathcal{U}) = - \sum_{m=1}^{\infty} (\gamma t)^m \int_{\mathcal{H}} \text{Tr} \{ M^{(m)} [\rho_0] \rho_0 \} d\rho_0. \quad (97) $$

## A.3 Integral over the Fubini-Study Measure

We can evaluate Eq. (97) to arrive at the result in Eq. (33). The integral over the Fubini-Study measure is calculated using methods from [39] relating to

$$ \begin{aligned} & \int_{\mathcal{H}} U_{ka} U_{ic} \bar{U}_{jb} \bar{U}_{rd} dU \\ & = \frac{1}{(d^2 - 1)} \left[ (\delta_{kj} \delta_{ab} \delta_{ir} \delta_{cd} + \delta_{kr} \delta_{ad} \delta_{ij} \delta_{cb}) \right] - \\ & \quad - \frac{1}{d(d^2 - 1)} \left[ (\delta_{kj} \delta_{ad} \delta_{ir} \delta_{cb} + \delta_{kr} \delta_{ab} \delta_{ij} \delta_{cd}) \right], \end{aligned} \quad (98) $$

where, and for the remainder of this section, the indices represent the matrix elements using Einstein sum notation.

Now, beginning with the element-wise definition of the superoperator term,

$$ (M[\rho])_{pi} = M_{kj,pi} \rho_{kj}, \quad (99) $$

where we identify $\rho$ and $\rho_0$ for convenient usage of the element-wise sum notation, this extends to

$$ (M[\rho])_{pq} = (M[\rho])_{pi} \rho_{iq} \quad (100) $$

$$ = M_{kj,pi} \rho_{kj} \rho_{iq}. \quad (101) $$

Taking the trace of this term corresponds to

$$ \text{Tr}\{M[\rho]\rho\} = \delta_{pq} (M[\rho]\rho)_{pq} \quad (102) $$

$$ = M_{kj,ri} \rho_{kj} \rho_{ir}, \quad (103) $$

where we define $r := p = q$. Rewriting the integral of this expression in terms of the associated unitary operator $U$ leads to

$$ \begin{aligned} & \int_{\mathcal{H}} \text{Tr} \{ M [U \rho U^{\dagger}] U \rho U^{\dagger} \} dU \\ & = \int_{\mathcal{H}} M_{kj,ri} (U_{ka} \rho_{ab} \bar{U}_{jb}) (U_{ic} \rho_{cd} \bar{U}_{rd}) dU \end{aligned} \quad (104) $$

$$ = \rho_{ab} M_{kj,ri} \rho_{cd} \int_{\mathcal{H}} (U_{ka} \bar{U}_{jb} U_{ic} \bar{U}_{rd}) dU. \quad (105) $$

If we now substitute the result of this integral from Eq. (98), and evaluate all of the delta function terms, we arrive at

$$ \begin{aligned} & \int_{\mathcal{H}} \text{Tr}\{M[\rho]\rho\}d\rho \\ & = \frac{1}{(d^2 - 1)} [\rho_{aa}\rho_{cc}M_{jj,ii} + \rho_{ac}\rho_{ac}M_{ki,ki}] - \\ & \quad - \frac{1}{d(d^2 - 1)} [\rho_{ac}\rho_{ac}M_{jj,ii} + \rho_{aa}\rho_{cc}M_{ki,ki}] \end{aligned} \quad (106) $$ $$ = \frac{1}{(d^2 - 1)} \left[ (\text{Tr}\{\rho\})^2 \text{Tr}\{M[1]\} + \text{Tr}\{\rho^2\} \text{Tr}\{M\} \right] - \\ & \quad - \frac{1}{d(d^2 - 1)} \left[ \text{Tr}\{\rho^2\} \text{Tr}\{M[1]\} + (\text{Tr}\{\rho\})^2 \text{Tr}\{M\} \right]. \quad (107) $$

Returning to the matrix forms from the Einstein sum notation, it can be seen that:

$$ \rho_{aa} \rho_{cc} = (\text{Tr}\{\rho\})^2, \quad (108) $$

$$ \rho_{ac} \rho_{ac} = \text{Tr}\{\rho^2\}, \quad (109) $$

$$ M_{jj,ii} = \text{Tr}\{M[1]\}, \quad (110) $$

$$ M_{ki,ki} = \text{Tr}\{M\}, \quad (111) $$
which results in

$$
\begin{align}
\int_{\mathcal{H}} \mathrm{Tr}\{M[\rho]\rho\}d\rho & \\
&= \frac{1}{(d^2-1)} \left[ (\mathrm{Tr}\{\rho\})^2 \mathrm{Tr}\{M[1]\} + \mathrm{Tr}\{\rho^2\} \mathrm{Tr}\{M\} \right] - \nonumber \\
&\quad - \frac{1}{d(d^2-1)} \left[ \mathrm{Tr}\{\rho^2\} \mathrm{Tr}\{M[1]\} + (\mathrm{Tr}\{\rho\})^2 \mathrm{Tr}\{M\} \right]. \tag{112}
\end{align}
$$

Finally, noting that Tr{ρ} = 1, Tr{ρ²} = 1 and ρ² = ρ, we arrive at

$$
\begin{align}
\int_{\mathcal{H}} \mathrm{Tr}\{M[\rho]\rho\}d\rho &= \frac{\mathrm{Tr}\{M\} + \mathrm{Tr}\{M[1]\}}{d(d+1)} \tag{113} \\
&= \frac{\mathrm{Tr}\{M\}}{d(d+1)}, \tag{114}
\end{align}
$$

since for the cases where M is composed of a traceless
lindbladian operator L, Tr{M[1]} = 0, and we arrive
at Eq. (33).

A.4 Boundedness of the AGI

We make use of the relation in Eq. (113) to prove that,
in the limit of strong dephasing, the AGI of a single
qudit is bounded above and below, with the bounds
depending only on the dimension of the system.

**Theorem A.2.** Consider an open quantum system of a single qudit of dimension *d*, initiated in a pure state Tr{$\rho_0^2$} = 1, under the influence of a time-independent noise superoperator *L* with coupling constant $\gamma$ and pure dephasing collapse operator *L* = *J*<sub>*z*</sub>, and evolving via an arbitrary time-independent unitary superoperator *S* over time *t*. Then, in the limit of strong coupling where ($\gamma$t) $\gg$ 1, the AGI *I* is bounded in the region

$$
1 - \frac{2}{d+1} \le \mathcal{I} \le 1 - \frac{1}{d+1}. \quad (115)
$$

Proof. In superoperator form, the quantum channel
can be written as $\mathcal{E} = \mathcal{U} \circ M$. Now, since the channel
is completely positive and trace-preserving (CPTP)
and $\mathcal{U}$ is a unitary operator, $M$ must also be CPTP.
Substituting this into the integral equation for the
AGF in Eq. (9) and using Eq. (113), we have

$$
\mathcal{F} = \frac{\mathrm{Tr}\{M\} + \mathrm{Tr}\{M[1]\}}{d(d+1)}. \quad (116)
$$

Therefore, proving the bound on $\mathcal{I}$ is equivalent to
proving that

$$
1 \le (d+1)\mathcal{F} \le 2, \qquad (117)
$$

$$
d \le \mathrm{Tr}\{M\} + \mathrm{Tr}\{M[1]\} \le 2d, \qquad (118)
$$

$$
0 \le \mathrm{Tr}\{M\} \le d, \qquad (119)
$$

since Tr{M[1]} = Tr{1} = d due to the trace-
preservation of M.

Now, to show this, consider the action of *M* on
the state *ρ*. The collapse operator *L* = *J*<sub>*z*</sub> of pure
dephasing is a diagonal matrix which preserves the
populations of *ρ* while causing the coherences to decay
to zero. Thus,

$$
M : \rho \rightarrow \sum_{i=0}^{d-1} \rho'_{ii} |i\rangle\langle i| \quad (120)
$$

$$
\forall \rho, M_{ij}^{kl} \rho_{kl} \propto \delta_{ij} \qquad (121)
$$

$$
M_{ij}^{kl} = \lambda_{i}^{kl} \delta_{ij}. \qquad (122)
$$

In the Kraus representation, using $M = \sum_{k=1}^{K} E_k^* \otimes E_k$, the matrix elements are given, trivially, by

$$
M_{ij}^{nm} = \sum_{k=1}^{K} E_{k,jm}^{*} E_{k,in}, \quad (123)
$$

and, since *M* is trace-preserving,

$$
\forall n, m, \sum_{i=0}^{d-1} M_{ii}^{nm} = \sum_{k=1}^{K} \sum_{i=0}^{d-1} E_{k,im}^{*} E_{k,in} \quad (124)
$$

$$
= \delta_{nm}. \tag{125}
$$

Furthermore,

$$
\mathrm{Tr}\{M\rho\} = \mathrm{Tr}\{\rho\} \qquad (126)
$$

$$
\mathrm{vec}(\mathbf{1})^{\mathrm{T}} \cdot \mathbf{M} \cdot \mathrm{vec}(\rho) = \mathrm{vec}(\mathbf{1})^{\mathrm{T}} \cdot \mathrm{vec}(\rho) \quad (127)
$$

$$
\therefore \sum_{n=0}^{d-1} M_{ij}^{nn} = \delta_{ij}. \qquad (128)
$$

Thus, in particular, if $n = m$:

$$
\sum_{k=1}^{K} \sum_{i=0}^{d-1} |E_{k,in}|^2 = 1 \quad (129)
$$

$$
\Leftrightarrow \sum_{k, i=n} |E_{k,in}|^2 + \sum_{k, i \neq n} |E_{k,in}|^2 = 1 \quad (130)
$$

$$
\therefore \forall n, 0 \le \sum_{k=1}^{K} |E_{k,nn}|^2 \le 1 \quad (131)
$$

Moreover, since

$$
\mathrm{Tr}\{M\} = \sum_{k,ij} E_{k,jj}^* E_{k,ii} \quad (132)
$$

$$
= \sum_k |\mathrm{Tr}\{E_k\}|^2, \qquad (133)
$$

the trace of *M* can be equivalently expressed by

$$
\mathrm{Tr}\{M\} = \sum_{i=0}^{d-1} \sum_{j=0}^{d-1} M_{ij}^{ij} \quad (134)
$$

$$
= \sum_{i=0}^{d-1} M_{ii}^{ii} \qquad (135)
$$

$$
= \sum_{k=1}^{K} \sum_{i=0}^{d-1} |E_{k,ii}|^2 . \qquad (136)
$$
Therefore, combining this last expression with Eq. (131) where each of the d elements $|E_{k,ii}|^2$ are bounded above by 1, it is clear that the sum of all d of them, and thus Tr{$M$}, must be bounded by d,

$$0 \le \sum_{k}^{K} \sum_{i=0}^{d-1} |E_{k,ii}|^2 \le d, \quad (137)$$

proving the boundedness of $\mathcal{F}$ and $\mathcal{I}$. $\square$

## A.5 Gate-Independent Correction Terms

In Sec. 2.2, we found in Eq. (34) that the m-th order AGI correction always contains a gate-independent term proportional to the trace of the m-th power of the collapse operator $\mathcal{L}$ in superoperator form. Specifically, in Sec. 2.5 Eqs. (55) and (63), we presented the first- and second-order gate-independent AGI correction terms.

Since we have in Eq. (13) an expression for $\mathcal{L}$ in terms of the (regular) collapse operators $L$, we can use this to express the correction terms in terms of $L$ instead. To first-order,

$$\mathrm{Tr}\{\mathcal{L}\} = \mathrm{Tr}\left\{L^* \otimes L - \frac{1}{2} \left(L^\dagger L \otimes 1_d + (1_d \otimes L^\dagger L)^T\right)\right\} \quad (138)$$

$$= \mathrm{Tr}\{L^* \otimes L\} - \frac{1}{2} \mathrm{Tr}\{L^\dagger L \otimes 1_d + 1_d \otimes (L^\dagger L)^T\} \quad (139)$$

$$= \mathrm{Tr}\{L^*\} \mathrm{Tr}\{L\} - \frac{d}{2} \left( \mathrm{Tr}\{L^\dagger L\} + \mathrm{Tr}\{(L^\dagger L)^T\} \right) \quad (140)$$

$$= |\mathrm{Tr}\{L\}|^2 - d \mathrm{Tr}\{L^\dagger L\}, \quad (141)$$

giving the result of Eq. (56).

By the same procedure, it is trivial to see that the second-order term in Eq. (68) is given by,

$$\mathrm{Tr}\{\mathcal{L}^2\} = |\mathrm{Tr}\{L^2\}|^2 + \frac{1}{2}|\mathrm{Tr}\{L^\dagger L\}|^2 + \frac{d}{2}\mathrm{Tr}\{(L^\dagger L)^2\}. \quad (142)$$

Furthermore, it is possible to make use of the multinomial expansion to find the m-th power trace to arbitrary order. In general, this produces a non-trivial sum over all L-words of length m. In the specific case of real and symmetric collapse operators $L = L^* = L^T = L^\dagger$, these words can be factorised,

resulting in the triangular sum

$$\mathrm{Tr}\{\mathcal{L}^m\} = \mathrm{Tr}\left\{\left(L^* \otimes L - \frac{1}{2} \left(L^\dagger L \otimes 1_d + (1_d \otimes L^\dagger L)^T\right)\right)^m\right\} \quad (143)$$

$$= m! \sum_{k_1+k_2+k_3=m} \frac{\mathrm{Tr}\{(L^*)^{k_1} (L^\dagger L)^{k_2} \otimes L^{k_1} ((L^\dagger L)^T)^{k_3}\}}{(-2)^{k_2+k_3} k_1! k_2! k_3!} \quad (144)$$

$$= m! \sum_{k_1+k_2+k_3=m} \frac{\mathrm{Tr}\{L^{k_1+2k_2}\} \mathrm{Tr}\{L^{k_1+2k_3}\}}{(-2)^{k_2+k_3} k_1! k_2! k_3!} \quad (145)$$

$$= m! \sum_{k_1=0}^{m} \sum_{k_2=0}^{m-k_1} \frac{\mathrm{Tr}\{L^{k_1+2k_2}\} \mathrm{Tr}\{L^{2m-k_1-2k_2}\}}{(-2)^{m-k_1-1} k_1! k_2! (m-k_1-k_2)!} \quad (146)$$

## A.6 Second-order Gate-Dependent Correction Term

Beginning with the general expression for the second-order term in the AGI expansion in Eq. (33), and combining with the $m=2$ term $M^{(2)}(t)$ in Eq. (30), we obtain Eq. (60) since

$$\mathrm{Tr}\{M^{(2)}(t)\} = \mathrm{Tr}\left\{\frac{\mathcal{L}^2}{2!} + \sum_{n_1=1}^{\infty} \sum_{n_2=1}^{\infty} \left(\prod_{i=1}^{2} \frac{(-t)^{n_i} [(S)^{n_i}, \mathcal{L}]}{n_i! \sum_{j=i}^{2} (n_j + 1)}\right)\right\} \quad (147)$$

$$= \mathrm{Tr}\left\{\frac{\mathcal{L}^2}{2!}\right\} + \sum_{n_1=1}^{\infty} \sum_{n_2=1}^{\infty} \frac{(-t)^{n_1+n_2} \mathrm{Tr}\{[(S)^{n_1}, \mathcal{L}][(S)^{n_2}, \mathcal{L}]\}}{n_1! n_2! (n_1+n_2+2)(n_2+1)} \quad (148)$$

Focusing now on the gate-dependent terms only, we can express the double summation of the trace term using the relation $n_1 + n_2 + 2 = s$

$$\sum_{n_1=1}^{\infty} \sum_{n_2=1}^{\infty} \frac{(-t)^{n_1+n_2} \mathrm{Tr}\{[(S)^{n_1}, \mathcal{L}][(S)^{n_2}, \mathcal{L}]\}}{n_1! n_2! (n_1+n_2+2)(n_2+1)} \quad (149)$$

$$= \sum_{s=0}^{\infty} \sum_{n=0}^{s-2} \frac{(-t)^{s-2} \mathrm{Tr}\{[(S)^n, \mathcal{L}][(S)^{s-2-n}, \mathcal{L}]\}}{n! s(s-1)!} \quad (150)$$

$$= \sum_{s=0}^{\infty} \frac{(-t)^{s-2}}{s(s-1)!} \times \sum_{n=0}^{s-2} \binom{s-1}{n} \mathrm{Tr}\{[(S)^n, \mathcal{L}][(S)^{s-2-n}, \mathcal{L}]\}. \quad (151)$$

**Lemma A.1.** Given operators A and B the following binomial sum of the product of iterated commutators
at order $n$ and $j-n$ is traceless,

$$ \sum_{n=1}^{j} \binom{j+1}{n} \mathrm{Tr}\{[(A)^n, B][(A)^{j-n}, B]\} = 0. \quad (152) $$

*Proof.* Using $(\binom{j+1}{n} = \binom{j+1}{n-1} + \binom{j+1}{n-1})$, we have

$$ \begin{aligned} & \sum_{n=1}^{j} \binom{j+1}{n} \mathrm{Tr}\{[(A)^n, B][(A)^{j-n}, B]\} \\ & = \sum_{n=1}^{j} \binom{j}{n-1} \mathrm{Tr}\{[(A)^n, B][(A)^{j-n}, B]\} + \\ & \qquad + \sum_{n=1}^{j} \binom{j}{n} \mathrm{Tr}\{[(A)^n, B][(A)^{j-n}, B]\}. \end{aligned} \quad (153) $$

Now, for the first term with binomial $(\binom{j+1}{n-1})$ and using the change of index $n \to n+1$,

$$ \begin{aligned} & \sum_{n=1}^{j} \binom{j}{n-1} \mathrm{Tr}\{[(A)^n, B][(A)^{j-n}, B]\} \\ & = \sum_{n=0}^{j-1} \binom{j}{n} \mathrm{Tr}\{[(A)^{n+1}, B][(A)^{j-n-1}, B]\} \quad (154) \\ & = \sum_{n=0}^{j-1} \binom{j}{n} \mathrm{Tr}\{[A, [(A)^n, B][(A)^{j-n-1}, B]\} \quad (155) \\ & = - \sum_{n=0}^{j-1} \binom{j}{n} \mathrm{Tr}\{[(A)^n, B][(A)^{j-n}, B]\}, \end{aligned} $$

by cyclic permutation of the trace. This summation is the negation of the second binomial term $(\binom{j+1}{n-1})$ for all terms $1 \le n \le j-1$. Therefore only the $n=0$ and $n=j$ terms remain. Moreover since $(\binom{j}{0}) = \binom{j}{0} = 1$, we have that

$$ \begin{aligned} & \sum_{n=1}^{j} \binom{j+1}{n} \mathrm{Tr}\{[(A)^n, B][(A)^{j-n}, B]\} \\ & = \mathrm{Tr}\{[(A)^j, B]B\} - \mathrm{Tr}\{B[(A)^j, B]\} \quad (157) \\ & = 0. \quad (158) \\ & \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \Box \end{aligned} $$

Now, applying this result to the summation in Eq. (151), where we identify $j=s-2$, we have that all terms for $n \ge 1$ are zero, and therefore only the $n=0$ term survives giving the following relation

$$ \begin{aligned} & \sum_{n=0}^{s-2} \binom{s-1}{n} \mathrm{Tr}\{[(S)^n, \mathcal{L}][(S)^{s-2-n}, \mathcal{L}]\} \\ & = \mathrm{Tr}\{\mathcal{L}[(S)^{s-2}, \mathcal{L}]\}. \end{aligned} \quad (159) $$

From this, with a change of index $s-2=s$, we can rewrite the full second-order correction as

$$ \mathcal{I}^{(2)} = \frac{-(\gamma t)^2}{d(d+1)} \times \left( \frac{\mathrm{Tr}\{\mathcal{L}^2\}}{2} + \sum_{s=2}^{\infty} \frac{(-t)^s \mathrm{Tr}\{\mathcal{L}[(S)^s, \mathcal{L}]\}}{(s+2)!} \right), \quad (160) $$

giving the result in Eq. (63).

Now, as for the gate-independent terms, we wish to express this gate-dependent trace of superoperators explicitly in terms of the operators $H$ and $L$.

Beginning with the binomial expansion of the trace of the iterated commutator [42],

$$ \begin{aligned} & \mathrm{Tr}\{\mathcal{L}[(S)^s, \mathcal{L}]\} \\ & = \sum_{k=0}^{s} (-1)^k \binom{s}{k} \mathrm{Tr}\{\mathcal{L}\mathcal{S}^k \mathcal{L}\mathcal{S}^{s-k}\} \quad (161) \\ & = i^s \sum_{k=0}^{s} \sum_{j=0}^{k} \sum_{l=0}^{s-k} \binom{s}{k} \binom{j}{l} (-1)^{s+k-j-l} \times \\ & \qquad \mathrm{Tr}\{\mathcal{L}((H^*)^{k-j} \otimes H^j) \mathcal{L}((H^*)^{s-k-l} \otimes H^l)\}, \end{aligned} \quad (162) $$

where we used the binomial expansion of the expression for $S$ in terms of $H$ in Eq. (12),

$$ \mathcal{S}^k = i^k \sum_{j=0}^{k} (-1)^{k-j} \binom{k}{j} (H^*)^{k-j} \otimes H^j. \quad (163) $$

Substituting now the expression for $\mathcal{L}$ in terms of $L$ in Eq. (13), and expanding algebraically the products of the operator terms in the trace, we obtain the following nine expressions, the sum of which is equivalent to the superoperator trace (omitting the prefactors and summations for convenience):

$$ \begin{aligned} & \mathrm{Tr}\{LH^{k-j}LH^{s-k-l}\}^* \mathrm{Tr}\{LH^jLH^l\}, \quad (164) \\ & -\frac{1}{2}\mathrm{Tr}\{LH^{k-j}L^T L^*H^{s-k-l}\}^* \mathrm{Tr}\{LH^{j+l}\}, \quad (165) \\ & -\frac{1}{2}\mathrm{Tr}\{LH^{s-j-l}\}^* \mathrm{Tr}\{LH^jL^T L^*H^l\}, \quad (166) \\ & -\frac{1}{2}\mathrm{Tr}\{L^T L^*H^{k-j}LH^{s-k-l}\}^* \mathrm{Tr}\{LH^{j+l}\}, \quad (167) \\ & \frac{1}{4}\mathrm{Tr}\{L^T L^*H^{k-j}L^T L^*H^{s-k-l}\}^* \mathrm{Tr}\{H^{j+l}\}, \quad (168) \\ & \frac{1}{4}\mathrm{Tr}\{L^T L^*H^{s-j-l}\}^* \mathrm{Tr}\{L^T L^*H^{j+l}\}, \quad (169) \\ & -\frac{1}{2}\mathrm{Tr}\{LH^{s-j-l}\}^* \mathrm{Tr}\{L^\dagger LH^j LH^l\}, \quad (170) \\ & \frac{1}{4}\mathrm{Tr}\{L^T L^*H^{s-j-l}\}^* \mathrm{Tr}\{L^\dagger LH^{j+l}\}, \quad (171) \\ & \frac{1}{4}\mathrm{Tr}\{H^{s-j-l}\}^* \mathrm{Tr}\{L^\dagger LH^j L^T L^*H^l\}. \quad (172) \end{aligned} $$

Now, we can simplify these expressions by making the following observations: Firstly, $s$ must be even since (i) heuristically, the trace must be real since we are calculating the fidelity (a real quantity) and therefore the prefactor $i^s$ in Eq. (161) cannot be imaginary, and (ii) formally, we prove this in Appendix A.7. Secondly, by symmetry of the binomial expansions, we are free to make the following change of variables:

$$ \begin{aligned} & j \to k-j, \\ & l \to s-k-l. \end{aligned} $$
For clarity, note that this change of variables does not affect the prefactor

$$(-1)^{s+k-j-l} \rightarrow (-1)^{j+k+l},$$

since, with $s$ even, the ratio of these two terms is

always a power of 2 and therefore equal to 1.

As an example to illustrate the effect of this change of variables, for the first term Eq. (164), the forward and backwards sums are complex conjugates and therefore only the real part can remain (omitting prefactors of the full expression):

$$\sum_{j=0}^{k} \sum_{l=0}^{s-k} \mathrm{Tr}\{LH^{k-j}LH^{s-k-l}\}^{*} \mathrm{Tr}\{LH^{j}LH^{l}\} \rightarrow \sum_{j=0}^{k} \sum_{l=0}^{s-k} \mathrm{Tr}\{LH^{j}LH^{l}\}^{*} \mathrm{Tr}\{LH^{k-j}LH^{s-k-l}\} \quad (173) \\ \Rightarrow \sum_{j=0}^{k} \sum_{l=0}^{s-k} \mathrm{Tr}\{LH^{k-j}LH^{s-k-l}\}^{*} \mathrm{Tr}\{LH^{j}LH^{l}\} = \sum_{j=0}^{k} \sum_{l=0}^{s-k} \mathrm{Re}\{\mathrm{Tr}\{LH^{k-j}LH^{s-k-l}\}^{*} \mathrm{Tr}\{LH^{j}LH^{l}\}\}. \quad (174)$$

Similarly, by comparing the remaining terms in Eqs. (165) to (172) with their counterparts following the change of variables, it is possible to identify the following sums as complex conjugates of one another:

$$\sum_{j=0}^{k} \sum_{l=0}^{s-k} (165 + 166) = \sum_{j=0}^{k} \sum_{l=0}^{s-k} 2 \mathrm{Re}\{165\}, \\ \sum_{j=0}^{k} \sum_{l=0}^{s-k} (167 + 170) = \sum_{j=0}^{k} \sum_{l=0}^{s-k} 2 \mathrm{Re}\{167\}, \\ \sum_{j=0}^{k} \sum_{l=0}^{s-k} (169 + 171) = \sum_{j=0}^{k} \sum_{l=0}^{s-k} 2 \mathrm{Re}\{169\}, \\ \sum_{j=0}^{k} \sum_{l=0}^{s-k} (168 + 172) = \sum_{j=0}^{k} \sum_{l=0}^{s-k} 2 \mathrm{Re}\{168\}.$$

Hence, by combining all of these terms into the full expression for $\mathrm{Tr}\{\mathcal{L}[(S)^s, \mathcal{L}]\}$ in Eq. (161), we obtain the following result of Eq. (71):

$$\mathrm{Tr}\{\mathcal{L}[(S)^s, \mathcal{L}]\} \\ = \sum_{k=0}^{s} \sum_{j=0}^{k} \sum_{l=0}^{s-k} \binom{s}{k} \binom{k}{j} \binom{s-k}{l} (-1)^{s+k-j-l} \times \\ \times \mathrm{Re} \left[ \begin{aligned} & \mathrm{Tr}\{LH^{k-j}LH^{s-k-l}\}^* \mathrm{Tr}\{LH^jLH^l\} \\ & + \frac{1}{2} \mathrm{Tr}\{L^T L^* H^{s-j-l}\}^* \mathrm{Tr}\{L^T L^* H^{j+l}\} \\ & + \frac{1}{2} \mathrm{Tr}\{L^T L^* H^{k-j} L^T L^* H^{s-k-l}\}^* \mathrm{Tr}\{H^{j+l}\} \\ & - \mathrm{Tr}\{LH^{k-j}L^T L^* H^{s-k-l}\}^* \mathrm{Tr}\{LH^{j+l}\} \\ & - \mathrm{Tr}\{L^T L^* H^{k-j} LH^{s-k-l}\}^* \mathrm{Tr}\{LH^{j+l}\} \end{aligned} \right]. \quad (175)$$

## A.7 Tracelessness of the Odd-Order Iterated Commutator

**Theorem A.3.** Given two operators A and B, the product of B with the iterated commutator of order s, $B[(A)^s, B]$, is traceless for all odd s,

$$\forall s = 2k + 1, k \in \mathbb{Z}, \mathrm{Tr}\{B[(A)^s, B]\} = 0. \quad (176)$$

*Proof.* We begin with the trace of the binomial expansion of the iterated commutator [42],

$$B[(A)^s, B] = \sum_{k=0}^{s} (-1)^k \binom{s}{k} BA^{s-k}BA^k \quad (177)$$

$$\therefore \mathrm{Tr}\{B[(A)^s, B]\} = \sum_{k=0}^{s} (-1)^k \binom{s}{k} \mathrm{Tr}\{BA^{s-k}BA^k\}. \quad (178)$$

By symmetry of the binomial expansion $\binom{s}{k} = \binom{s}{s-k}$, and considering $s$ odd $(-1)^{s-k} = -(-1)^k$, we can write the equivalent expansion

$$B[(A)^s, B] = - \sum_{k=0}^{s} (-1)^k \binom{s}{k} BA^k BA^{s-k} \quad (179)$$

$$\therefore \mathrm{Tr}\{B[(A)^s, B]\} = - \sum_{k=0}^{s} (-1)^k \binom{s}{k} \mathrm{Tr}\{BA^k BA^{s-k}\} \quad (180) \\ = - \sum_{k=0}^{s} (-1)^k \binom{s}{k} \mathrm{Tr}\{BA^{s-k}BA^k\}, \quad (181)$$

by cyclic permutation of the trace, $\mathrm{Tr}\{CD\} = \mathrm{Tr}\{DC\} \Rightarrow \mathrm{Tr}\{[C, D]\} = 0$. Hence, we have the result for odd s

$$\mathrm{Tr}\{B[(A)^s, B]\} = - \mathrm{Tr}\{B[(A)^s, B]\} \quad (182)$$

$$\Rightarrow \mathrm{Tr}\{B[(A)^s, B]\} = 0. \quad (183)$$
<figure><img src="image_12.png" /></figure>

# B Gates and Operators

Here we compile additional information and examples regarding the quantum gates and collapse operators used in the main text.

## B.1 Quantum Gates

Given the set of single-qudit gates $U_g \in \{1_d, X, Z, F\}$ defined in Eqs. (16) to (19), we present them in matrix form to illustrate their structure. The identity matrix in $d$ dimensions is trivial, and generated by a null control Hamiltonian $H_c = \mathbb{0}_d$, corresponding to no action on the qudit state. The generalized Pauli-X or SHIFT gate $X = \Sigma_x$ is defined as

$$X = \begin{bmatrix} 0 & 1 & 0 & 0 & \cdots & 0 \\ 0 & 0 & 1 & 0 & \cdots & 0 \\ 0 & 0 & 0 & 1 & \cdots & 0 \\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & 0 & \cdots & 1 \\ 1 & 0 & 0 & 0 & \cdots & 0 \end{bmatrix}, \quad (184)$$

which cyclically permutes the basis states, incrementing them all by 1 level. An equivalent formulation is given by the transpose $X^T$, which cyclically permutes the states downwards,

$$X^T = \begin{bmatrix} 0 & 0 & 0 & \cdots & 0 & 1 \\ 1 & 0 & 0 & \cdots & 0 & 0 \\ 0 & 1 & 0 & \cdots & 0 & 0 \\ 0 & 0 & 1 & \cdots & 0 & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & 0 & \cdots & 1 & 0 \end{bmatrix}. \quad (185)$$

The generalized Pauli-Z or CLOCK gate $Z = \Sigma_z$ is defined as

$$Z = \begin{bmatrix} 1 & 0 & 0 & \cdots & 0 \\ 0 & \omega & 0 & \cdots & 0 \\ 0 & 0 & \omega^2 & \cdots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & 0 & \omega^{d-1} \end{bmatrix}, \quad (186)$$

where the $d$-th roots of unity, $\omega$, are also used to define the matrix form of the Quantum Fourier Transform in $d$ dimensions, being the generalisation of the Walsh-Hadamard (superposition) gate,

$$F = \frac{1}{\sqrt{d}} \times \begin{bmatrix} 1 & 1 & 1 & 1 & \cdots & 1 \\ 1 & \omega & \omega^2 & \omega^3 & \cdots & \omega^{d-1} \\ 1 & \omega^2 & \omega^4 & \omega^6 & \cdots & \omega^{2(d-1)} \\ 1 & \omega^3 & \omega^6 & \omega^9 & \cdots & \omega^{3(d-1)} \\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & \omega^{d-1} & \omega^{2(d-1)} & \omega^{3(d-1)} & \cdots & \omega^{(d-1)(d-1)} \end{bmatrix}. \quad (187)$$

Note that we may consider a further generalisation of the CLOCK gate from a phase shift of the roots of unity to an arbitrary phase $\phi$,

$$\text{PHASE}(\phi) = \begin{bmatrix} 1 & 0 & 0 & \cdots & 0 \\ 0 & e^{i\phi} & 0 & \cdots & 0 \\ 0 & 0 & e^{2i\phi} & \cdots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & 0 & e^{(d-1)i\phi} \end{bmatrix}, \quad (188)$$

where the choice of phase factor allows us to generate not only Clifford gates but also non-Clifford gates, such as for the $T$ gate with $\phi = \pi/8$.

Now, we introduce the interpolated $X$ (SHIFT) and $Z$ (CLOCK) gates, where we raise the $X$ and $Z$ gates to the power $\eta \in [0, 1]$, which allows us to smoothly interpolate the action of the gate(s) between the identity ($\eta = 0$) and the original gate ($\eta = 1$). As an example, consider $\eta = \frac{3d-2}{4d}$ for $d = 2$:

$$X^{0.5} = \frac{1}{\sqrt{2}} \begin{bmatrix} e^{i \frac{\pi}{4}} & e^{-i \frac{\pi}{4}} \\ e^{-i \frac{\pi}{4}} & e^{i \frac{\pi}{4}} \end{bmatrix}. \quad (189)$$

It can be seen that for $0 < \eta < 1$ the $X^\eta$ and $Z^\eta$ gates are non-Clifford.

## B.2 Collapse Operators

We present a brief example for the collapse operators of pure dephasing, bit-flip and relaxation, as defined in Eqs. (22) to (24). In $d = 4$, these are precisely the generalised spin operators $J_z$, $J_x$ and $J_-$ for spin-$\frac{3}{2}$ systems, given by

$$J_z = \frac{1}{2} \begin{bmatrix} 3 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & -1 & 0 \\ 0 & 0 & 0 & -3 \end{bmatrix}, \quad (190)$$

and

$$J_x = \frac{1}{2} \begin{bmatrix} 0 & \sqrt{3} & 0 & 0 \\ \sqrt{3} & 0 & 2 & 0 \\ 0 & 2 & 0 & \sqrt{3} \\ 0 & 0 & \sqrt{3} & 0 \end{bmatrix}, \quad (191)$$

and

$$J_{-} = \begin{bmatrix} 0 & 0 & 0 & 0 \\ \sqrt{3} & 0 & 0 & 0 \\ 0 & 2 & 0 & 0 \\ 0 & 0 & \sqrt{3} & 0 \end{bmatrix}. \quad (192)$$

# C Complementary Results

Here, we present further results and investigations that support the findings of the main text.
<figure><img src="image_13.png" /><figcaption>Figure 11: Method QUTIP，《rand_unitary_haar》(and equivalently CIRQ，《rand_unitary》), based on [44].</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Figure 12: Method QUTIP，《rand_unitary》(and equivalently QUTIP，《rand_herm》).</figcaption></figure>

## C.1 Uniformity of the CUE for Haar-random Unitaries

As noted in Sec. 2.3, the ability to generate ensembles of well-distributed random matrices for arbitrary qudit dimensions is of fundamental importance to the study of the statistical behaviour of the AGI, particularly regarding the aspect of gate-dependence. Given a qudit of dimension $d$, we identify the space of all possible quantum gates with the unitary group $U(d)$ of $d \times d$ unitary matrices, associated to skew-Hamiltonian control terms generated through the matrix exponential map of the Lie algebra $u(d)$. Therefore, the aim is to sample from a uniform distribution over this group.

Since $U$ is a compact, connected Lie group, it has a unique (up to positive scalar multiplication) and operation-invariant measure, referred to as the Haar measure [44]. When normalised to one, it produces a uniform probability measure on $U$, with the associated probability space referred to by the name Circular Unitary Ensemble (CUE), as originally introduced by Dyson [45]. Elements of the CUE are thus typically referred to as Haar-random unitary matrices. The properties of these random matrices have been well-studied, particularly in the standard reference text on random matrix theory by Mehta [46], to which we shall refer. For the specific case of the generation of random unitary matrices, we shall refer to the work of Mezzadri [44].

The spectral statistics form the foundation for the study of random matrices, and in particular we shall concern ourselves with the two primary quantities: (1) the eigenvalue distribution (level density), and (2) the nearest-neighbour eigenvalue ($n$-level, for $n = 1$) spacing distribution.

Since we are studying unitary matrices whose eigen-

values all lie on the unit circle in the complex plane, $\lambda_j = e^{i\theta_j}$, the eigenvalue distribution is thus the distribution of the complex phases $\theta_j$ of the $d$ eigenvalues:

$$ -\pi \le \theta_j < \pi \quad \forall \ 1 \le j \le d. \qquad (193) $$

In the large matrix limit where $d \to \infty$ the eigenvalue phases will tend towards a uniform distribution on the domain $[-\pi, \pi)$ with probability density $\rho(\theta) = \frac{1}{2\pi}$.

Next, the nearest-neighbour eigenvalue spacing distribution $p(s)$ is obtained by: (1) sorting the spectrum of phases $\theta_j$ from smallest to largest, (2) unfolding this spectrum by the factor $\frac{d}{2\pi}$ to normalise the mean spacing to 1, and then (3) taking their successive differences $s_j = \frac{d}{2\pi}(\theta_{j+1} - \theta_j)$, resulting in a list of $d-1$ elements. This distribution for the unitary ensembles is well approximated by the Wigner Surmise $p_w(s)$,

$$ p_w(s) = \frac{32s^2}{\pi^2} e^{-\frac{4s^2}{\pi}}, \qquad (194) $$

where $s$ is taken in the continuous limit of $s_j$. In the original work of Wigner [47] on the Gaussian Unitary Ensemble (GUE) of random Hamiltonian (Hermitian) matrices, this is the exact analytical form for the case of $2 \times 2$ matrices, which applies to very good approximation for large $d$. Additionally since the GUE and CUE are related functionally, these two distributions become equivalent for both ensembles for large $d$.

Given our requirements of sampling random matrices uniformly from the unitary group, we note that there exist numerous (PYTHON-based) packages with methods for generating unitary matrices. Some common, but non-exhaustive, examples include: QUTIP [48, 49] (rand_unitary, rand_unitary_haar, rand_herm), CIRQ [50] (rand_unitary, rand_special_unitary), and BRISTOL [51] (gen_cue).
<figure><img src="image_15.png" /><figcaption>(a) $d = 100$ level density</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>(b) $d = 100$ spacings</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>(c) $d = 2$ level density</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>(d) $d = 2$ spacings</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>(a) $d = 100$ level density</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>(b) $d = 100$ spacings</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>(c) $d = 2$ level density</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>(d) $d = 2$ spacings</figcaption></figure>

<caption>Figure 14: Method CIRQrand special unitary</caption>

Both QU TIP’s **rand unitary haar** and CIRQ’s
**rand unitary** methods implement precisely the al-
gorithm provided by Mezzadri in [44]. The
CIRQ `.rand_special_unitary` method augments the
original `rand_unitary` method by normalising the de-
tterminant of the randomly generated matrix to 1. The
QU TIP `.rand_unitary` method uses the `rand_herm`
method to generate a random Hermitian matrix and
then computes its unitary equivalent by matrix expo-
nentiation.

We tested each of these different methods by sam-
pling one million eigenvalues in each of 2 matrix size
configurations, *d* ∈ {2, 100} for each method. These
two values were chosen to represent the two regimes
of large and small matrix dimensions. The number of
eigenvalues was chosen such that, for the algorithm of
Mezzadri the standard deviation of the mean for the
level density of the eigenvalue phases of the *d* = 100
matrices fell below the threshold of 0.001. The num-
ber of eigenvalues was kept constant at each dimen-
sion by varying the number of quantum gates gen-
erated, specifically 500000 and 10000 for *d* = 2 and
*d* = 4 respectively.

In Figs. 11 to 14, the simulated data of eigen-
value and spacing distributions (blue) are compared
to their analytically expected results (red). We can
see that the method of Mezzadri in Fig. 11 performs
the best for generating uniformly distributed Haar-
random gates. Hence this is our method of choice.
However, even for this method the spacing distribu-
tion does show reduced accuracy for lower dimensions,
despite the level densities remaining uniform. This
could be attributed to the fact that the Wigner Sur-
mise is exact for 2 × 2 Hermitian matrices (but still
very accurate for large *d*), but it is only in the large *d*
limit that the distributions for the CUE are identical
to that of the GUE [46]. Thus some deviation is to
be expected for small unitary matrices.

The QU TIP `.rand_unitary` and QU TIP `.rand herm`
methods in Fig. 12 do not appear suitable for sam-
pling uniformly from the CUE. Indeed, inspection of
the sampled matrices showed that approximately 50%
of the generated unitary matrices were diagonal. On
the other hand, the BRISTOL `.gen cue` in Fig. 13 is
well optimised for high-dimensional systems but is
unsuitable for small dimensions. Finally, it is inter-
esting to note that the special unitary matrices of
CIRQ `.rand_special_unitary` in Fig. 14 are similarly
distributed for large dimension, but show singular be-
haviour for small *d*.

C.2 Distribution of the asymptotic AGIs

We continue our investigations on the different
PYTHON methods for generating Haar-random gates,
namely: BRISTOL.gen cue, QU TIP `.rand_unitary`,
QU TIP `.rand_unitary haar`, CIRQ `.rand_unitary`
and CIRQ `.rand_special_unitary`.

In Sec. 2.3, and Fig. 3 specifically, we discussed the
distribution of the plateau values *I*\* of the AGI in the
large-γ*t* regime. Recall that we refer to these plateau
values *I*\* at the large-γ*t* limit, lim<sub>γ*t*→1</sub> *I* → *I*\*.
We showed that they must be bounded above and below
by Eq. (47) for all gates. We also observed the ap-
pearance of two distinct behaviours (*I*\*<sub>mono</sub> and *I*\*<sub>over</sub>)
in the sampled Haar-random gates: (1) Monotonic
growth of the AGI towards *I*\*, and (2) overshoot with
a single turning point above *I*\*.

We note that, as for all the numerical results in the
main text, the `QU TIP `.rand_unitary_haar` method
was used to generate the data shown in Fig. 3. We
complement this work by repeating these simulations
using the remaining methods (for *d* = 2) and also with
`QU TIP `.rand_unitary_haar` for *d* ∈ {2, 4, 8, 16}. For
each test method and dimension, 100000 gates were
generated, their *I*\* calculated numerically and their
<figure><img src="image_23.png" /><figcaption>Figure 15: Probability distributions of the AGI plateaus for different gate generation methods. 100000 gates were generated using each of the four methods for $d = 2$: BRISTOL.gen_cue in Fig. 15a, QUTIPrand_unitary in Fig. 15b, CIRQ，《rand_unitary》in Fig. 15c and CIRQ，《rand_special》in Fig. 15d. For each gate, the $\mathcal{I}^*$ was calculated and categorised according to whether the curves were monotonic (orange) or overshooting (blue). The distributions of the counts of $\mathcal{I}^*$ in the domain [$\mathcal{I}_{\min}^*, \mathcal{I}_{\max}^*$] (black and red vertical lines, respectively) were binned according to Sturge's formula $n_{bins} = \log_2(n_{gates} + 1)$ and normalised such that the total distribution $p(\mathcal{I}_{total}^*)$ (grey) had unit area. The mean plateau value $\mathbb{E}(\mathcal{I}_{total}^*)$ (cyan, dotted vertical line), was plotted in comparison with the expected mean $\bar{\mathcal{I}}^* = 1 - \frac{1}{d}$ (blue, dashed vertical line). The green line represents the expected uniform distribution $u(\mathcal{I}_{min}^*, \mathcal{I}_{max}^*)$ in the domains.</figcaption></figure>
<figure><img src="image_24.png" /><figcaption>(a) $d=2$</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>(c) $d=8$</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>(b) $d=4$</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>(d) $d=16$</figcaption></figure>

Figure 16: Probability distributions of the AGI plateaus for dimensions $d \in \{2, 4, 8, 16\}$. 100000 gates were generated using the method QUOTIPrand unitary haar of [44] at each dimension: $d=2$ in Fig. 16a, $d=4$ in Fig. 16b, $d=8$ in Fig. 16c, $d=16$ in Fig. 16d. For each gate, the $\mathcal{I}^*$ was calculated and categorised according to whether the curves were monotonic (orange) or overshooting (blue). The distributions of the counts of $\mathcal{I}^*$ in the domain $[\mathcal{I}_{\min}^*, \mathcal{I}_{\max}^*]$ (black and red vertical lines, respectively) were binned according to: Sturge's formula $n_{\text{bins}} = \log_2(n_{\text{gates}} + 1)$ for $d=2$, $n_{\text{bins}} = 100$ for $d=4, 8$ and $n_{\text{bins}} = 200$ for $d=16$, due to the increasing concentration of the counts for the larger dimensions. The counts were normalised such that the total distribution $p(\mathcal{I}_{\text{total}}^*)$ (grey) had unit area. The mean plateau value $\mathbb{E}(\mathcal{I}_{\text{total}}^*)$ (cyan, dotted vertical line), was plotted in comparison with the expected mean $\bar{\mathcal{I}}^* = 1 - \frac{1}{d}$ (blue, dashed vertical line). The green curve represents the fitted Weibull distributions $f(\mathcal{I}^*; \mathcal{I}_0^*, k, \lambda) = \frac{k}{\lambda} \left(\frac{\mathcal{I}^* - \mathcal{I}_0^*}{\lambda}\right)^{k-1} e^{-\left(\frac{\mathcal{I}^* - \mathcal{I}_0^*}{\lambda}\right)^k}$ in the domain.
behaviour categorised as either monotonic or overshooting as shown in Table 1.

Figures 15 and 16 display the data of these experiments by plotting histograms of the probability densities of the $\mathcal{I}^*$ for each category of gate, $\mathcal{I}_{\text{mono}}^*$ in orange and $\mathcal{I}_{\text{over}}^*$ in blue, as well as their combined distribution $\mathcal{I}_{\text{total}}^*$ in grey. The histograms were normalised such that the area under the curve of $p(\mathcal{I}_{\text{total}}^*)$ was equal to 1. The red and black vertical lines show the upper and lower bounds on $\mathcal{I}^*$, respectively. The dashed blue vertical lines show the expected mean value for each dimension.

Looking closer at Fig. 15, the two CIRQ methods in Figs. 15c and 15d appear to have distributions matching most closely with that of QuTiP `.rand.unitary_haar` in Fig. 16a, although the normalisation procedure for the special unitary matrices appears to generate a larger proportion of gates with overshoot. It is interesting to note that the `BRISTOL.gen.cue` method in Fig. 15a generates an increasing number of gates near the expected mean, but only below, with no gates in the range $1 - \frac{1}{d} \le \mathcal{I}^* \le 1 - \frac{1}{d+1}$. As mentioned in Appendix C.1, QuTiP `.rand.unitary` in Fig. 15b generates a large number of diagonal unitary matrices, whose $\mathcal{I}^* = \mathcal{I}_{\text{min}}^* = 1 - \frac{2}{d+1}$ are equivalent to that of the identity matrix $\mathbb{1}_d$.

Figure 16 shows the preferred QuTiP `.rand.unitary_haar` method for different dimensions $d \in \{2, 4, 8, 16\}$. We note that at each dimension the computed mean over all sampled $\mathcal{I}^*$ is equal to (within statistical error) the expected mean shown by the dashed blue vertical line, $\tilde{\mathcal{I}}^* = 1 - \frac{1}{d}$. However, most striking is the fact that for $d > 2$ the total probability density is no longer uniformly distributed between $\mathcal{I}_{\text{min}}^*$ and $\mathcal{I}_{\text{max}}^*$, but becomes increasingly concentrated about the mean.

We remark that this observed behaviour is due to the concentration of measure phenomenon as the dimensionality of the system increases [52, 53]. To understand this, we first demonstrate the following property of the AGI:

**Lemma C.1.** The AGI function $\mathcal{I}: \mathbf{U}(d) \to \mathbb{R}$ is L-Lipschitz continuous with respect to the Frobenius norm $\|A\|_F = \sqrt{\text{Tr}\{A^\dagger A\}}$ over the unitary group $\mathbf{U}(d)$ with Lipschitz constant $L = \frac{2}{d+1}$, such that

$$|\mathcal{I}(U) - \mathcal{I}(U')| \le L \|U - U'\|_F \quad (195)$$

where $U, U' \in \mathbf{U}(d)$.

*Proof.* The AGF $\mathcal{F}$ can be written in terms of the process fidelity $F$ as [41]

$$\mathcal{F} = \frac{dF + 1}{d + 1}, \quad (196)$$

where $F = \frac{1}{d^2} |\text{Tr}\{U^\dagger V\}|^2$ for two unitaries $U, V \in$

$\mathbf{U}(d)$. Then, with $\mathcal{I} = 1 - \mathcal{F}$, we have that

$$\mathcal{I} = \frac{d^2 - |\text{Tr}\{U^\dagger V\}|^2}{d(d+1)}. \quad (197)$$

Considering the difference, we have that

$$\begin{align}
& |\mathcal{I}(U) - \mathcal{I}(U')| \nonumber \\
& = \frac{1}{d(d+1)} \left| \text{Tr}\{U^\dagger V\} \right|^2 - \left| \text{Tr}\{U'^\dagger V\} \right|^2 \tag{198} \\
& = \frac{1}{d(d+1)} \left| \text{Tr}\{U^\dagger V\} \right| - \left| \text{Tr}\{U'^\dagger V\} \right| \times \nonumber \\
& \qquad \times \left| \text{Tr}\{U^\dagger V\} \right| + \left| \text{Tr}\{U'^\dagger V\} \right|, \tag{199}
\end{align}$$

by difference of squares. Now,

$$\begin{align}
& \left| \text{Tr}\{U^\dagger V\} + \text{Tr}\{U'^\dagger V\} \right| \nonumber \\
& \le \left| \text{Tr}\{U^\dagger V\} \right| + \left| \text{Tr}\{U'^\dagger V\} \right| \tag{200} \\
& \le 2d, \tag{201}
\end{align}$$

by the Cauchy-Schwartz inequality and since the trace of a unitary is at most $d$. Similarly,

$$\begin{align}
& \left| \text{Tr}\{U^\dagger V\} - \text{Tr}\{U'^\dagger V\} \right| \nonumber \\
& \le \left| \text{Tr}\{U^\dagger V\} - \text{Tr}\{U'^\dagger V\} \right| \tag{202} \\
& \le \left\| U^\dagger V - U'^\dagger V \right\|_F \tag{203} \\
& \le \left\| (U - U')^\dagger V \right\|_F \tag{204} \\
& \le \|U - U'\|_F, \tag{205}
\end{align}$$

and by unitary invariance of the Frobenius norm. Therefore, substituting these two expressions,

$$|\mathcal{I}(U) - \mathcal{I}(U')| \le \frac{2}{d+1} \|U - U'\|_F, \quad (206)$$

and thus $\mathcal{I}$ is Lipschitz with $L = \frac{2}{d+1}$. $\square$

Now, we can apply this result for $\mathcal{I}$ in conjunction with Levy's Lemma extended to the classical compact group $\mathbf{U}(d)$ [54], which we restate here without proof:

**Theorem C.1.** Let $f: X \to \mathbb{R}$ be L-Lipschitz continuous with Lipschitz constant $L$, and let $X$ be an element of one of the compact classical groups. Then, for each $t > 0$,

$$\mathbb{P} \left[ |f(X) - \mathbb{E}[f(X)]| \ge t \right] \le e^{-\frac{(d-2)t^2}{24L^2}}. \quad (207)$$

Using our result for $\mathcal{I}$, we have that

$$\mathbb{P} \left[ |\mathcal{I}(U) - \mathbb{E}[\mathcal{I}(U)]| \ge t \right] \le e^{-\frac{(d-2)(d+1)^2 t^2}{96}}. \quad (208)$$

This therefore confirms our observation of the concentration of measure phenomenon in Fig. 16; it tells us that the probability of observing a gate whose infidelity is further than $t$ from the expectation value of
<table><caption>**Table 1:** Proportions of Random Gates with and without Turning Points. 100000 gates were generated for each of the method and dimension pairs shown. The $I^*_{mono}$ (%) and $I^*_{over}$ columns show the percentages of gates falling into either of the two categories: (1) monotonic growth, and (2) overshooting behaviour. The column $\epsilon_{rel}$ gives the relative error between the measured mean of the overall data $E(I^*_{total})$ and the expected mean $\bar{I}^* = 1 - \frac{1}{d}$ for the calculated AGI plateau values, given by $\epsilon_{rel} = \left| \frac{E(I^*_{total}) - \bar{I}^*}{\bar{I}^*} \right|$. The quantity $D_{KL}$ is the Kullback-Leibler divergence measuring the relative entropy between the distribution $p(I^*_{total})$ and the reference distribution: $D_{KL}(p(I^*_{total})||u(I^*_{min}, I^*_{max}))$ for $d = 2$ uniform and $D_{KL}(p(I^*_{total})||f(I^*; I^*_0, k, \lambda))$ for $d > 2$ Weibull.</caption><thead><tr><th>d</th><th>Method</th><th>I<sup>*</sup><sub>mono</sub> (%)</th><th>I<sup>*</sup><sub>over</sub> (%)</th><th>ε<sub>rel</sub> (%)</th><th>D<sub>KL</sub></th></tr></thead><tbody><tr><td>2</td><td>BRISTOL.gen_cue</td><td>59.56</td><td>40.44</td><td>0.100</td><td>0.885</td></tr><tr><td>2</td><td>QUdTIPrand_unitary</td><td>100.00</td><td>0.00</td><td>0.243</td><td>0.438</td></tr><tr><td>2</td><td>CIRQrand_unitary</td><td>73.06</td><td>26.94</td><td>1.02 × 10<sup>-3</sup></td><td>5.92 × 10<sup>-5</sup></td></tr><tr><td>2</td><td>CIRQrand_special_unitary</td><td>61.66</td><td>38.34</td><td>6.83 × 10<sup>-4</sup></td><td>1.32 × 10<sup>-4</sup></td></tr><tr><td>2</td><td>QUdTIPrand_unitary_haar</td><td>75.51</td><td>24.49</td><td>2.56 × 10<sup>-4</sup></td><td>7.15 × 10<sup>-5</sup></td></tr><tr><td>4</td><td>QUdTIPrand_unitary_haar</td><td>71.24</td><td>28.76</td><td>4.38 × 10<sup>-5</sup></td><td>2.92 × 10<sup>-3</sup></td></tr><tr><td>8</td><td>QUdTIPrand_unitary_haar</td><td>70.21</td><td>29.79</td><td>1.27 × 10<sup>-5</sup></td><td>1.21 × 10<sup>-3</sup></td></tr><tr><td>16</td><td>QUdTIPrand_unitary_haar</td><td>99.95</td><td>0.05</td><td>7.98 × 10<sup>-4</sup></td><td>1.36 × 10<sup>-3</sup></td></tr></tbody></table>

<figure><img src="image_28.png" /><figcaption>**Figure 17:** Fitted power-law model parameters as functions of $\eta$ for the saturation points $(\gamma t)^*$ of the $X^\eta$ gates. AGIs for $X^\eta$ with $\eta \in (0, 1]$ for $d \in [2, 12]$. Saturation points $(\gamma t)^*$ were calculated by interpolation and root-finding. A power law model was fitted with parameters $\alpha(\eta)$ (blue, circles), $\beta(\eta)$ (red, triangles), and $\delta(\eta)$ (black, squares) were then fitted to sigmoidal and exponential models (dashed lines) as functions of $\eta$.</figcaption></figure>

the AGI decays exponentially in the distance $t^2$ and the dimensionality $(d-2)(d+1)^2$ of the system.

For the case $d=2$, the exponent is zero and thus the probability is constant for all $t$ and so we see in Fig. 16a that the AGI plateau values are uniformly distributed. However, as $d > 2$ for Figs. 16b to 16d, the $I^*$ become exponentially concentrated about the mean $\bar{I}^*$.

This could also account for the fluctuation in the ratios of $I^*_{mono}/I^*_{over}$ at different system sizes. However, further work is needed to elucidate the precise reasons for these observations.

### C.3 Power-law behaviour of the saturation points of the AGI

As shown in Fig. 5 in Sec. 2.4, the saturation points $(\gamma t)^*$ of the AGIs for the interpolated $X^\eta$ gates were observed to follow a power-law model described by Eq. (51). The initial analysis was performed for each gate in a range of values $\eta = (0, 1]$ for $d \in [2, 12]$ by: (1) computing the AGIs numerically for a range of $\gamma t$ values in PYTHON using QUTIP's propagator and process_fidelity methods, (2) interpolating the AGI curves using SciPy.interpolate.CubicSpline, and (3) finding the plateau saturation points $(\gamma t)^*$ by root-finding with SciPy optimize root_scalar.</sub>

Then, this three-parameter power law model in Eq. (51) was fitted on $(\gamma t)^*$ as a function of $d$ for each $\eta$ curve by nonlinear least-squares regression using the SciPy optimize curve_fit method. This produced three sets of data for the fitted parameters $\alpha, \beta, \gamma$ as functions of $\eta$. These three curves were then themselves fitted to their respective sigmoidal and exponential models in Eqs. (52) to (54), using the same SciPy optimize curve_fit method. The models showed good agreement with the data, with $R^2(\alpha) = 0.9992, R^2(\beta) = 0.9992$ and $R^2(\delta) = 0.9989$.

The fitted data points and resulting models are shown in Fig. 17. We note that the outliers at small values of $\eta$ (including $\eta = 0$ for the identity gate with constant $(\gamma t)^*$) are due to the near-constant saturation point values of these gates over the range of $d$, which slowly transition towards more power-law-like behaviour as $\eta$ increases. For example, for $\eta = 0$ the saturation point is constant in $d$, giving a horizontal line. Fitting a power-law curve to this dataset was possible but required unrealistic fit parameters like $\alpha \sim 10^{50}$. Nevertheless, this analysis hints at a potential approach for estimating the saturation points based only on the parameter $\eta$, without requiring the explicit numerical simulation.</sub>
<figure><img src="image_29.png" /><figcaption>Figure 18: Linear regression fits of the critical order $s_ε$ over gate time $t$ for the X gate. The X gate ($η = 1$) was simulated for dimensions $d \in \{2, 4, 8, 16, 32\}$ over gate times $γt \in [0.5, 4.5]$ with $γ = 1$ fixed. The second-order gate-dependent correction term was calculated at each $t$ for increasing order $s$ in the summation until convergence between successive terms was obtained at the critical value $s_ε$. For each $d$, a linear model was fitted to the data as shown by the dotted lines.</figcaption></figure>

### C.4 Time dependence of the second-order correction to the AGI

As alluded to in Fig. 9 in Sec. 2.5, the convergence of the second-order (gate-dependent) AGI correction term as a function of $t$ was observed to be linear in Eq. (66). More specifically, the number of iterations $s_ε$ required to reach convergence of two successive terms in the summation Eq. (63) was calculated according to Eqs. (64) and (65).

To understand the time-dependence of this, for each of the $η$ and $d$ values shown in Fig. 9, the critical $s_ε$ was computed for a range of times $t \in [0, 5]$ with $γ$ fixed at 1. Interestingly, we can see that in Fig. 18 the $s_ε$ does indeed appear linear in $t$ for each dimension shown for the X gate. This behaviour extends to other $η$ values too. Therefore, a linear least-squares regression was performed on these data points using the SciPy `optimize.curve_fit`; the resulting fitted models are shown in Fig. 18 alongside the data points, indicating good agreement. The resulting gradients $m = \frac{ds_ε}{dt}$ for each fitted curve are thus the data shown in Fig. 9.
References

[1] Olivier Ezratty. "Perspective on superconducting qubit quantum computing". The European Physical Journal A 59, 94 (2023).

[2] Nathalie P. De Leon, Kohei M. Itoh, Dohun Kim, Karan K. Mehta, Tracy E. Northup, Hanhee Paik, Benjamin S. Palmer, Nitin Samarth, Sorawis Sangtawesin, and David W. Steuerman. "Materials challenges and opportunities for quantum computing hardware". Science 372, eabb2823 (2021).

[3] Eufemio Moreno-Pineda, Clément Godfrin, Franck Balestro, Wolfgang Wernsdorfer, and Mario Ruben. "Molecular Spin Qudits for Quantum Algorithms". Chemical Society Reviews 47, 501-513 (2018).

[4] Eufemio Moreno-Pineda and Wolfgang Wernsdorfer. "Measuring molecular magnets for quantum technologies". Nature Reviews Physics 3, 645-659 (2021).

[5] Yulin Chi, Jieshan Huang, Zhanchuan Zhang, Jun Mao, Zinan Zhou, Xiaojiong Chen, Chonghao Zhai, Jueming Bao, Tianxiang Dai, Huihong Yuan, Ming Zhang, Daoxin Dai, Bo Tang, Yan Yang, Zhihua Li, Yunhong Ding, Leif K. Oxenløwe, Mark G. Thompson, Jeremy L. O'Brien, Yan Li, Qihuang Gong, and Jianwei Wang. "A programmable qudit-based quantum processor". Nature Communications 13, 1166 (2022).

[6] Martin Ringbauer, Michael Meth, Lukas Postler, Roman Stricker, Rainer Blatt, Philipp Schindler, and Thomas Monz. "A universal qudit quantum processor with trapped ions". Nature Physics 18, 1053-1057 (2022).

[7] Yuchen Wang, Zixuan Hu, Barry C. Sanders, and Sabre Kais. "Qudits and High-Dimensional Quantum Computing". Frontiers in Physics 8, 589504 (2020).

[8] Alessandro Chiesa, Francesco Petiziol, Emilio Macaluso, Sandro Wimberger, Paolo Santini, and Stefano Carretta. "Embedded quantum-error correction and controlled-phase gate for molecular spin qubits". AIP Advances 11, 025134 (2021).

[9] Francesco Petiziol, Alessandro Chiesa, Sandro Wimberger, Paolo Santini, and Stefano Carretta. "Counteracting dephasing in Molecular Nanomagnets by optimized qudit encodings". npj Quantum Information 7, 133 (2021).

[10] Hideyuki Miyahara, Yiyou Chen, Vwani Roychowdhury, and Louis-Serge Bouchard. "Decoherence mitigation by embedding a logical qubit in a qudit". Quantum Information Processing 22, 278 (2023).

[11] Lane G. Gunderman. "Local-dimension-invariant qudit stabilizer codes". Physical Review A 101, 052343 (2020).

[12] Wolfgang Wernsdorfer and Mario Ruben. "Synthetic Hilbert Space Engineering of Molecular Qudits: Isotopologue Chemistry". Advanced Materials 31, 1806687 (2019).

[13] Yunzhe Zheng, Hemant Sharma, and Johannes Borregaard. "Entanglement Distribution with Minimal Memory Requirements Using Time-Bin Photonic Qudits". PRX Quantum 3, 040319 (2022).

[14] Frédéric Bouchard, Robert Fickler, Robert W. Boyd, and Ebrahim Karimi. "High-dimensional quantum cloning and applications to quantum hacking". Science Advances 3, e1601915 (2017).

[15] Nikolay Petrovich Brusentsov and José Ramil Alvarez. "Ternary Computers: The Setun and the Setun 70". In John Impagliazzo and Eduard Proydakov, editors, Perspectives on Soviet and Russian Computing. Volume 357 of IFIP Advances in Information and Communication Technology, pages 74-80. Springer Berlin Heidelberg, Berlin, Heidelberg (2011).

[16] Sergio Boixo, Sergei V. Isakov, Vadim N. Smelyanskiy, Ryan Babbush, Nan Ding, Zhang Jiang, Michael J. Bremner, John M. Martinis, and Hartmut Neven. "Characterizing quantum supremacy in near-term devices". Nature Physics 14, 595-600 (2018).

[17] Frank Arute, Kunal Arya, Ryan Babbush, Dave Bacon, Joseph C. Bardin, Rami Barends, Rupak Biswas, Sergio Boixo, Fernando G. S. L. Brandao, David A. Buell, Brian Burkett, Yu Chen, Zijun Chen, Ben Chiaro, Roberto Collins, William Courtney, Andrew Dunsworth, Edward Farhi, Brooks Foxen, Austin Fowler, Craig Gidney, Marissa Giustina, Rob Graff, Keith Guerin, Steve Habegger, Matthew P. Harrigan, Michael J. Hartmann, Alan Ho, Markus Hoffmann, Trent Huang, Travis S. Humble, Sergei V. Isakov, Evan Jeffrey, Zhang Jiang, Dvir Kafri, Kostyantyn Kechedzhi, Julian Kelly, Paul V. Klimov, Sergey Knysh, Alexander Korotkov, Fedor Kostritsa, David Landhuis, Mike Lindmark, Erik Lucero, Dmitry Lyakh, Salvatore Mandrà, Jarrod R. McClean, Matthew McEwen, Anthony Megrant, Xiao Mi, Kristel Michielsen, Masoud Mohseni, Josh Mutus, Ofer Naaman, Matthew Neeley, Charles Neill, Murphy Yuezhen Niu, Eric Ostby, Andre Petukhov, John C. Platt, Chris Quintana, Eleanor G. Rieffel, Pedram Roushan, Nicholas C. Rubin, Daniel Sank, Kevin J. Satzinger, Vadim Smelyanskiy, Kevin J. Sung, Matthew D. Trevithick, Amit Vainsencher, Benjamin Villalonga, Theodore White, Z. Jamie
Yao, Ping Yeh, Adam Zalcman, Hartmut Neven, and John M. Martinis. "Quantum supremacy using a programmable superconducting processor". Nature 574, 505–510 (2019).

[18] Clément Godfrin, Abdelkarim Ferhat, Rafik Ballou, Svetlana Klyatskaya, Mario Ruben, Wolfgang Wernsdorfer, and Franck Balestro. "Operating Quantum States in Single Magnetic Molecules: Implementation of Grover's Quantum Algorithm". Physical Review Letters 119, 187702 (2017).

[19] Stefan Thiele, Franck Balestro, Rafik Ballou, Svetlana Klyatskaya, Mario Ruben, and Wolfgang Wernsdorfer. "Electrically driven nuclear spin resonance in single-molecule magnets". Science 344, 1135–1138 (2014).

[20] Matthew Otten, Keshav Kapoor, A. Bariş Özgüler, Eric T. Holland, James B. Kowalkowski, Yuri Alexeev, and Adam L. Lyon. "Impacts of noise and structure on quantum information encoded in a quantum memory". Physical Review A 104, 012605 (2021).

[21] Denis Janković, Jean-Gabriel Hartmann, Mario Ruben, and Paul-Antoine Hervieux. "Noisy qudit vs multiple qubits: conditions on gate efficiency for enhancing fidelity". npj Quantum Information 10, 59 (2024).

[22] Michael A. Nielsen and Isaac L. Chuang. "Quantum computation and quantum information". Cambridge University Press. Cambridge, NY (2000).

[23] Michael A. Nielsen. "A simple formula for the average gate fidelity of a quantum dynamical operation". Physics Letters A 303, 249–252 (2002).

[24] Göran Lindblad. "On the generators of quantum dynamical semigroups". Communications in Mathematical Physics 48, 119–130 (1976).

[25] Vittorio Gorini, Andrzej Kossakowski, and E. C. George Sudarshan. "Completely positive dynamical semigroups of $N$-level systems". Journal of Mathematical Physics 17, 821–825 (1976).

[26] Heinz-Peter Breuer and Francesco Petruccione. "The Theory of Open Quantum Systems". Oxford University Press. Oxford (2007). First edition.

[27] Daniel Manzano. "A short introduction to the Lindblad master equation". AIP Advances 10, 025106 (2020).

[28] John Watrous. "The Theory of Quantum Information". Cambridge University Press. Cambridge (2018). First edition.

[29] Ingemar Bengtsson and Karol Zyczkowski. "Geometry of Quantum States: An Introduction to Quantum Entanglement". Cambridge University Press. Cambridge (2006). First edition.

[30] Timothy F. Havel. "Robust procedures for converting among Lindblad, Kraus and matrix representations of quantum dynamical semigroups". Journal of Mathematical Physics 44, 534–557 (2003).

[31] Terry A. Loring. "Computing a logarithm of a unitary matrix with general spectrum". Numerical Linear Algebra with Applications 21, 744–760 (2014).

[32] Maxwell Aifer and Sebastian Deffner. "From quantum speed limits to energy-efficient quantum gates". New Journal of Physics 24, 055002 (2022).

[33] Apostolos Vourdas. "Quantum systems with finite Hilbert space". Reports on Progress in Physics 67, 267–320 (2004).

[34] Mark Howard and Jiri Vala. "Qudit versions of the qubit $\pi/8$ gate". Physical Review A 86, 022316 (2012).

[35] Jose Martínez-Carranza, Francisco Soto-Eguibar, and Héctor M. Moya-Cessa. "Alternative analysis to perturbation theory in quantum mechanics: Dyson series in matrix form". The European Physical Journal D 66, 22 (2012).

[36] Braulio M. Villegas-Martínez, Francisco Soto-Eguibar, and Héctor M. Moya-Cessa. "Application of Perturbation Theory to a Master Equation". Advances in Mathematical Physics 2016, 1–7 (2016).

[37] John E. Campbell. "On a Law of Combination of Operators bearing on the Theory of Continuous Transformation Groups". Proceedings of the London Mathematical Society s1-28, 381–390 (1896).

[38] Brian C. Hall. "Lie groups, Lie algebras, and representations: an elementary introduction". Number 222 in Graduate texts in mathematics. Springer. Cambridge, NY (2015). Second edition.

[39] Benoit Collins, Sho Matsumoto, and Jonathan Novak. "The Weingarten Calculus". Notices of the American Mathematical Society 69, 1 (2022).

[40] Alessandro Chiesa, Emilio Macaluso, Francesco Petiziol, Sandro Wimberger, Paolo Santini, and Stefano Carretta. "Molecular Nanomagnets as Qubits with Embedded Quantum-Error Correction". The Journal of Physical Chemistry Letters 11, 8610–8615 (2020).

[41] Michał Horodecki, Paweł Horodecki, and Ryszard Horodecki. "General teleportation channel, singlet fraction, and quasidistillation". Physical Review A 60, 1888–1898 (1999).
[42] Howard C. Volkin, United States. National Aeronautics, Space Administration, and Lewis Research Center. "Iterated commutators and functions of operators". NASA technical note. National Aeronautics and Space Administration. (1968). url: https://ntrs.nasa.gov/citations/19680027053.

[43] Jean-Gabriel Hartmann, Denis Janković, Rémi Pasquier, Mario Ruben, and Paul-Antoine Hervieux. "Code and Data for "Nonlinearity of the Fidelity in Open Qudit Systems: Gate and Noise Dependence in High-dimensional Quantum Computing"". Zenodo. (2025).

[44] Francesco Mezzadri. "How to generate random matrices from the classical compact groups". Notices of the American Mathematical Society 54, 592-604 (2007).

[45] Freeman J. Dyson. "The Threefold Way. Algebraic Structure of Symmetry Groups and Ensembles in Quantum Mechanics". Journal of Mathematical Physics 3, 1199-1215 (1962).

[46] Madan L. Mehta. "Random matrices". Number 142 in Pure and applied mathematics series. Elsevier. Amsterdam (2004). Third edition.

[47] Eugene P. Wigner. "Characteristic Vectors of Bordered Matrices With Infinite Dimensions". The Annals of Mathematics 62, 548 (1955).

[48] Jesper R. Johansson, Paul D. Nation, and Franco Nori. "QuTiP: An open-source Python framework for the dynamics of open quantum systems". Computer Physics Communications 183, 1760-1772 (2012).

[49] Jesper R. Johansson, Paul D. Nation, and Franco Nori. "QuTiP 2: A Python framework for the dynamics of open quantum systems". Computer Physics Communications 184, 1234-1240 (2013).

[50] Cirq Developers. "Cirq". Google Quantum AI. (2023).

[51] Mehmet Süzen, Cornelius Weber, and Joan J. Cerdà. "Spectral Ergodicity In Deep Learning Architectures Via Surrogate Random Matrices" (2017). arXiv:1704.08303.

[52] Apostolos A. Giannopoulos and Vitali D. Milman. "Concentration Property on Probability Spaces". Advances in Mathematics 156, 77-106 (2000).

[53] Michel Ledoux. "The concentration of measure phenomenon". Number Issue 89 in Mathematical surveys and monographs. American Mathematical Society. Providence, RI (2001).

[54] Elizabeth S. Meckes. "The Random Matrix Theory of the Classical Compact Groups". Cambridge University Press. Cambridge (2019). First edition.