§ INTRODUCTION

While quantum processors continue to progress towards practical use, the errors present in current systems are still the most limiting factor. A dominant error in superconductor-based quantum computers is decoherence.
There have been several proposals to mitigate it, either by designing new qubit types[1, 2, 3] or by further optimizing the existing designs, typically focusing on increasing the coherence times of the circuit[4, 5, 6]. However, in the latter case, one often encounters trade-offs between different circuit properties. For example in transmon qubits, which have emerged as the most widely-used qubit type in large-scale experiments[7, 8, 9, 10, 11], the suppression of charge noise comes at the cost of low anharmonicity which sets a lower bound on the duration of single-qubit operations[12].

This illustrates the importance of understanding the different errors affecting the quantum hardware, as well as the fact that an informed design of the circuit parameters requires the acknowledgment of a plethora of possible error sources which are not necessarily limited only to the coherence properties of the circuit, but also include leakage errors in single qubit gates due to low anharmonicity, state preparation errors due to finite-temperature heating effects[13, 14, 15], as well as the parity-switching error presented in this manuscript. To elaborate further, only taking into account the coherence properties of the transmon and the low anharmonicity, one way to achieve better performance is to increase the transmon anharmonicity (in order to suppress potential leakage errors) while keeping the frequency of the transmon fixed. The latter condition, under the assumption of constant quality factors, ensures the coherence properties of the circuit remain unchanged. However this inadvertently leads into the regime where the transmon charge dispersion becomes more significant. A good understanding of the charge-dispersion related errors is therefore needed in order to make informed design choices of the circuit parameters. One error that is exponentially more pronounced in the low $E_J/E_C$ regime is related to the presence of charge-like excitations of the superconducting condensate, referred to as quasiparticles and the charge-parity of the transmon, thus prompting us to analyze these effects further so that a trade-off between the different error sources can be made.

Quasiparticles can be created through several mechanisms, and are known to cause different types of incoherent errors in superconducting qubit realizations [16, 17]. Particularly, quasiparticle tunnelling across the Josephson junction results in energy relaxation and dephasing in superconducting qubits[18, 19, 20, 21, 22, 23].
Such detrimental quasiparticle-induced effects have, in turn, motivated research for finding mitigation strategies such as normal-metal traps [24, 25, 26], band-gap engineering [27, 28, 29] and improved qubit design [30]. There have also been efforts towards designing new types of superconducting qubits that are expected to be intrinsically robust against quasiparticle tunneling[31].

Suppression of charge-noise susceptibility of the transmon is achieved by adding a large shunt capacitor in parallel with a Josephson junction of a Cooper-pair box [12]. However, the energy levels of the transmon are not completely independent of the offset charge on the transmon island, but instead have a weak $2e$-periodic charge dispersion. Since the presence of a quasiparticle shifts the island charge by $e$, the energy spectrum of the transmon can, thus, be divided into two distinct manifolds based on the parity of the number of quasiparticles on the island as shown in Fig. <ref>(a).
The switching of the charge-parity can occur either due to a preexisting quasiparticle tunneling across the Josephson junction on to the transmon island, or due to photon-assisted breaking of a Cooper pair, as pictured in Fig. <ref>(b)[16]. The timescale of these stochastic parity-switches is referred to as the parity-switching time and, as we argue in the Results Section, it is typically much shorter compared to the time needed to obtain meaningful statistics from the quantum computer.

Since the difference between the two parity manifolds is strongly suppressed, quasiparticle effects in transmons have mainly been analyzed in the context of quasiparticle-induced decoherence. However, the parity-dependent energy splitting of the higher-excited states is much larger compared to that of the first-excited and ground states, as shown in Fig. <ref>(a). Therefore, parity-switching in the second excited state can potentially become a notable source of error for example in diabatic CZ (controlled-Z) or more generally CPHASE (controlled-phase) gates realized in the tunable coupler architecture[32, 33, 34, 35, 36, 37, 38, 39, 40, 41].

A schematic representation of the effect of a parity switch on a single transmon and on a two-qubit gate. (a) The energy diagram of a single offset-charge-sensitive transmon with $E_J/E_C = 10$ with two distinct parity manifolds, marked with + and -. While the ground state $|0\rangle$ also comprises of two distinct parity levels, the difference between them is not visible. (b) Illustration of two parity switching mechanisms. The vertical axis represents the energy, with the left and right regions corresponding to the two sides of the Josephson junction with the middle region corresponding to the insulator. The light grey area corresponds to the density of states of a BCS superconductor on both sides of the junction. Orange: A high energy photon breaks a Cooper pair (dark blue) thus generating two quasiparticles (light blue), with one quasiparticle tunneling across the junction. Green: A preexisting quasiparticle tunnels across the junction. (c) The energy level diagram of the states involved in the operation of a diabatic CPHASE gate. The Rabi oscillation between the levels $|101\rangle \leftrightarrow |002\rangle$ is marked with blue arrows, however the larger charge dispersion of the second excited state means that both parity states of the $|002\rangle$ level cannot be considered degenerate anymore. (d) The lumped element model of the tunable-coupler circuit used in the implementation of diabatic CPHASE gates, consisting of two computational transmons (dark and light blue) referred to in the following as Qubits 1 and 2 (with indices $q_1$ and $q_2$), a flux-tunable coupler (green), denoted with the index $c$, and capacitive couplings between the transmons (black). The readout resonators and drive lines for the implementation of single-qubit gates are not included in the schematics.

Besides the ability to perform fast and high-fidelity gate operations, one of the main reasons for the introduction of the tunable coupler is the fact that the static residual-ZZ interaction between the qubits can be completely suppressed by tuning the coupler transmon to a specific frequency[34, 32], thus leading to reduced cross-talk and spectator qubit related decoherence[42]. When determining this frequency, it has been shown that the higher levels of the tunable coupler system are relevant[35, 32], thus implying that quasiparticle dynamics might affect our ability to effectively decouple the qubits.

The tunable-coupler setup is currently one of the most promising platforms for large-scale quantum computing[7, 8, 10, 11] and, thus, it is important to study its susceptibility to parity-switching errors during quantum gates and idling.

In this paper, we first develop an analytical theory of parity switches in a tunable-coupler based architecture, and demonstrate that the effects of a parity switch in a two-qubit gate can be a relevant source of error, even in the transmon regime. Moreover, we show that this previously unidentified error can, in certain regimes, be the dominant quasiparticle-induced error mechanism during a diabatic two-qubit gate, indicated by a comparison to currently achievable parity switching rates observed in superconducting circuits. Furthermore, we demonstrate that the inherent stochastic nature of parity switching events limits the ability to suppress any unwanted longitudinal interactions between the qubits coupled through tunable couplers[43]. We find that the magnitude of the unwanted interactions make this effect relevant as coherence times advance into the millisecond regime[44, 45].

Secondly, we develop an analysis of the optimal circuit parameters which is based on the accurate modeling of these different error sources as well as on a performance metric related to the gate infidelities of a given quantum circuit. We analyze the dependence of the performance metric as a function of transmon qubit design parameters, and we find an optimal range of the parameters for which the performance metric is maximized. This demonstrates how informed design choices can aid the improvement of current quantum computers, by taking into account the different possible error mechanisms. Therefore, we believe our analysis and results can serve as a guideline when designing transmon-based quantum processors.

§ RESULTS

§.§ Modeling the Dynamics

Here, we present how a parity switch affects the tunable coupler architecture and derive analytical results describing the magnitude of the effect. The Hamiltonian of an individual transmon, not taking into account potential higher order contributions to the Josephson energy [46], is given by [12, 47]

\begin{equation}\label{eq:full_transmon_hamiltonian}
\hat{H} = 4E_C\left(\hat{n} - n_g + \frac{P - 1}{4}\right)^2 - E_J \cos{\hat{\phi}},
\end{equation}

where the operator $\hat{n}$ represents the dimensionless charge and $\hat{\phi}$ is the superconducting phase operator across the Josephson junction. The variables $E_C$, $E_J$ and $n_g$ represent the charging energy of an electron (i.e. the energy required to add a single electron of the Cooper-pair to the transmon island), Josephson energy and dimensionless offset charge, respectively. The variables $\hat{n}$ and $\hat{\phi}$ are related via the canonical commutation relation $[\hat{\phi},\hat{n} ] = i$. Additionally we have included a discrete parity variable $P \in \{-1,+1\}$, corresponding to the parity of the number of electrons that have tunneled across the junction. We observe that the parity term has the same effect as a shift of the offset charge by $\Delta n_g = 1/2$.

The eigenenergies of the Hamiltonian in Eq. <ref> have been analytically determined already in Ref. [12] and the low-energy manifold of a single transmon can be approximated in the asymptotic limit of $E_J/E_C \gg 1$ as

\begin{equation}\label{eq:duffing_with_parity}
\hat{H}/\hbar \simeq \left[\omega + \delta\omega(P,n_g) \right] \hat{a}^\dagger \hat{a} + \frac{\alpha + \delta\alpha(P,n_g)}{2}\hat{a}^\dagger\hat{a}^\dagger \hat{a} \hat{a},
\end{equation}

where $\hat{a}$ are bosonic annihilation operators, and $\hbar\omega\simeq \sqrt{8 E_J E_C} - E_C$ and $\hbar\alpha \simeq -E_C$ are the expressions for the transmon (angular) frequency and the anharmonicity in the asymptotic limit, respectively. Here, we have already taken into account the fact that the different parities have almost identical parameters, and the small differences between them are taken into account with the two parameters depending on the parity, $ \delta\omega(P,n_g)$ and $\delta\alpha(P,n_g)$.
The parity $P$ therefore divides the eigenstates of the Hamiltonian in Eq. <ref> into two distinct manifolds, as illustrated in Fig. <ref>(a).

Denoting the eigenenergies of the original transmon Hamiltonian in Eq. <ref> as $E_i$ with $i \in \{0,1,2,3,... \}$, the difference between the energy levels of the different parity states is asymptotically approximated by [12]

\begin{equation}\label{eq:parity_difference}
E_m^+(n_g) -  E_m^-(n_g) \simeq \epsilon_m \cos(2\pi n_g),
\end{equation}

where the superscript refers to the parity and the charge dispersion $\epsilon_m$ is given by

\begin{equation}\label{eq:asymptotic_charge_dipersion}
\epsilon_m \simeq (-1)^m E_C \frac{2^{4m + 5}}{m!}\sqrt{\frac{2}{\pi}}\left( \frac{E_J}{2 E_C}\right)^{\frac{m}{2} + \frac{3}{4}} e^{-\sqrt{8E_J /E_C}}.
\end{equation}

While the exponential suppression of the charge dispersion with the ratio $E_J/E_C$ is well-known and the main reason for the introduction of the transmon, the formula in Eq. <ref> also predicts a significant increase in the charge dispersion of higher excited states. This means that even though the effect of a parity switch may be small in the computational subspace, the effect can be significantly more pronounced if higher-excited states are involved in the operation of two-qubit gates. More specifically $|\epsilon_2/\epsilon_1| \simeq 2^{5/2}\sqrt{E_J/E_C} \sim 40$, for $E_J/E_C \sim 50$. This difference is even more pronounced since certain effects scale with the square of the charge dispersion, as we show in the following. Due to this rapid scaling, we neglect the effect on the first excited (i.e. $\delta\omega = 0$) state and only focus on the second excited state.

The parameters of the Hamiltonian in Eq. <ref> are related to the eigenenergies of the Hamiltonian in Eq. <ref> via the relations

\begin{align}
E_1^P &= \omega + \delta\omega(P,n_g), \\
E_2^P &= 2\omega + \alpha + 2\,\delta\omega(P,n_g) + \delta\alpha(P,n_g),
\end{align}

and by additionally defining the parity-averaged frequency and anharmonicity as $\hbar \omega = [E_1^+(n_g) + E_1^-(n_g)]/2$ and $\hbar (2\omega + \alpha) = [E_2^+(n_g) + E_2^-(n_g)]/2$[22]. Together with Eq. <ref> and by neglecting the first-excited state charge dispersion (i.e. $\epsilon_2 \gg \epsilon_1$) due to the aforementioned reasons, we arrive at $\delta\alpha(P,n_g) = P\epsilon_2\cos(2\pi n_g)/2$.

§.§.§ The Diabatic $CPHASE$ Gate

We consider a non-adiabatic, i.e. diabatic, CPHASE gate based on the two-qubit gate scheme using tunable couplers that was analyzed in Refs. [32, 33, 34, 35, 36] with similar schemes proposed in Refs. [37, 38, 39, 40]. We show the circuit schematics of the tunable-coupler setup in Fig. <ref>(d). Here, the two computational transmons, which we refer to as Qubits 1 and 2 ($q_{1,2}$), are capacitively coupled with each other and to a third, frequency-tunable, transmon which is referred to as the tunable coupler or simply coupler ($c$).  The introduction of the coupler enables on-demand gate operation between the two computational qubits.

Such a system can be modeled with the Hamiltonian [34, 35]

\begin{align}\label{eq:tqg_ham}
\hat{H}/\hbar &= \sum_{i \in \{q_1,c,q_2\}} \omega_i \hat{a}_i^\dagger \hat{a}_i + \frac{\alpha_i}{2} \hat{a}_i^\dagger\hat{a}_i^\dagger \hat{a}_i \hat{a}_i
- \sum_{\substack{i,j \in \{q_1,c,q_2\} \\ i\neq j}} g_{ij} (\hat{a}_i^\dagger - \hat{a}_i) (\hat{a}_j^\dagger - \hat{a}_j).
\end{align}

The dependence of the transmon parameters on the parity and offset charge is not explicitly shown above in order to simplify the notation. However, we note that the whole system has $2^3 = 8$ distinct parity states.

The main operation principle of the diabatic CPHASE gate is shown in Fig. <ref>(c). The CPHASE gate is implemented by tuning the frequency of the coupler using a flux pulse closer to the frequency of the computational transmons. The conditional phase is collected during a Rabi oscillation between the $|101\rangle$ and $|002\rangle$ states of the system, as illustrated in Fig. <ref>(c). Here, the labeling of the states is defined as $|q_1\, c\, q_2\rangle$. Since the second excited state of one of the computational transmons is significantly populated during the gate operation, the non-degeneracy of the two parity levels can have a direct effect, and therefore quasiparticle tunneling and photon-assisted pair breaking in the transmon can become a notable source of error in the gate operation. We emphasise here that the couplings between the transmons $g_{ij}  = \beta_{ij} \sqrt{\omega_i \omega_j}$ also depend on the frequencies, meaning that while $g_{q_1 q_2}$ is constant, $g_{q_{1,2} c}$ is time dependent. The dimensionless prefactors $\beta_{ij}$ depend on the coupling capacitances, as well as self-capacitances of the transmons in the lumped-element circuit model[34].

The computational basis in this scheme is formed by the eigenstates of the Hamiltonian from Eq. <ref> in the idling configuration (defined below), rather than the local (uncoupled) transmon states [35]. Since the couplings act only as a perturbation to the uncoupled states, we identify the full Hamiltonian eigenstates corresponding to the uncoupled states. More specifically, the computational state $|ij\rangle$, $i,j \in \{0,1\}$ is the eigenstate $|\psi\rangle$ of the Hamiltonian in Eq. <ref> with the maximal overlap $|\langle \psi |i_{q_1} 0_c j_{q_2} \rangle|$.
This notation is employed throughout this manuscript and we omit the subscript indices $q_{1,2}$ and $c$ from now on. The kets with three indices (e.g. $|ijk\rangle$) always denote the local (uncoupled) Fock states of the three-transmon system. The kets with only two indices (e.g. $|ij\rangle$) are used to denote the eigenstates of the whole system that are closest to the local (uncoupled) state $|i0j\rangle$. We denote with $\omega_{ij}$ the angular frequency of the computational state $|ij\rangle$. In general, the state $|11\rangle$ accumulates conditional phase with the rate

\begin{equation}\label{eq:zz_interaction}
\zeta_\mathrm{ZZ} = \omega_{11} - \omega_{01} - \omega_{10} + \omega_{00}.
\end{equation}

However, since the coupler frequency is tunable, it is typically possible to find one or two frequencies $\omega_c$ for which $\zeta_{\rm ZZ}=0$[34, 32, 37]. These special frequencies are referred to as the coupler idling frequencies and denoted with $\omega_c^\mathrm{idle}$.

A variety of pulse shapes can be used to implement the gate effectively and without inducing too many unwanted transitions[35, 32, 48]. Even though our analytical analysis makes minimal assumptions about the pulse shape, we need to choose a specific form for the numerical simulations. In our simulations, we use the so-called flattop Gaussian pulse described by the following formula[32, 49]

\begin{equation}\label{eq:flattop_gaussian}
f(t) = \frac{1}{2}\left[\mathrm{erf}\left(\frac{t - \tau_b}{\sqrt{2}\sigma}\right) - \mathrm{erf}\left(\frac{t - \tau_b  - \tau_c}{\sqrt{2}\sigma}\right) \right] - C.
\end{equation}

The flattop Gaussian is obtained by a convolution of a step function with duration $\tau_c$ and a Gaussian with parameter $\sigma$. The reasoning behind this choice is that the convolution of the flattop pulse with a Gaussian strongly suppresses the spectral component of the flattop pulse at higher frequencies, thus reducing the probability for unwanted transitions. An additional rise time of $\tau_b$ is also introduced, which we fix to $\tau_b = 2\sqrt{2}\sigma$ in the remainder of this paper. Since any gate must have a finite duration, we introduce a cut-off at time $T = 2\tau_b + \tau_c$. The constant $C$ is then introduced to correct for the discontinuity at the beginning and the end of a pulse with finite duration.

Consequently, the coupler frequency during the gate operation is given by

\begin{equation}
\omega_c = \omega_c^\mathrm{idle} + A f(t),
\end{equation}

where $A$ is the pulse amplitude. In the simulations, we optimize the pulse amplitude and duration $\tau_c$ such that the process fidelity of the gate is maximized.

§.§.§ Effective Model

The tunable coupler circuit Hamiltonian in Eq. <ref> is difficult to analyze and we typically rely on numerical studies [35]. It is therefore beneficial to introduce an effective Hamiltonian that can approximate the physics of the system.
Similar to Refs. [50, 34], we introduce the Schrieffer-Wolff transformation as a means to decouple the computational states from the coupler, and assume that the decoupled coupler remains in the ground state during the gate operations. Unlike the approach taken in Refs. [50, 34], where the Hilbert space of the local transmons is truncated to the computational subspace, in our case, we must also include the $|02\rangle$ state to account for the Rabi oscillation used to accumulate the conditional phase. More details of the transformation can be found in the Methods section.

Consequently, we consider here the truncated Hilbert space spanned by the states $\{|00\rangle, |01\rangle,|10\rangle,|11\rangle,|02\rangle\}$. In our analytical  considerations, we thus neglect the possibility of leakage, which is a good assumption for high-fidelity gates. We emphasize that the leakage effects are included in our numerical simulations which have been made with the full Hamiltonian defined in Eq. <ref>.

We additionally neglect the counter-rotating terms within the truncated Hilbert space and, thus, obtain the effective Hamiltonian

\begin{equation}\label{eq:effective_hamiltonian}
\hat{H}_\mathrm{eff}/\hbar \hat{=}
\begin{blockarray}{cccccc}
|00\rangle & |01\rangle & |10\rangle & |11\rangle & |02\rangle \\
\begin{block}{(ccc|cc)c}
0 & 0 & 0 & 0 & 0 &\, |00\rangle \\
0 & \tilde{\omega}_{q_2} & \tilde{g}_{01,10} & 0 & 0 &\, |01\rangle \\
0 &  \tilde{g}_{01,10} & \tilde{\omega}_{q_1} & 0 & 0 &\, |10\rangle \\
\cline{1-5}
0 & 0 & 0 & \tilde{\omega}_{q_1} + \tilde{\omega}_{q_2} & \tilde{g}_{11,02} &\, |11\rangle \\
0 & 0 & 0 & \tilde{g}_{11,02} &  2\tilde{\omega}_{q_2} + \tilde{\alpha}_{q_2} &\, |02\rangle \\
\end{block}
\end{blockarray},
\end{equation}

where the parameters with the tilde denote the perturbed parameters of the original full Hamiltonian in Eq. <ref>, which are derived in the Methods section (See Eqs. <ref>-<ref>) up to second order in $g_{q_i c}/(\omega_{q_i} - \omega_{c})$. Without the loss of generality, we have also assumed that the second-excited state of the second qubit ($q_2$) is used to perform the Rabi oscillation.

The pulse shape used to implement a diabatic CPHASE gate can typically be divided into three distinct stages: (i) fast sweep of the coupler frequency to an operation value close to the resonance with the qubits; (ii) long constant-frequency plateau at the operation frequency; (iii) fast sweep back to idling frequency of the coupler.

Since the plateau is typically longer compared to the fast sweeps [32], and the effective coupling strengths $\tilde{g}_{01,10}$ and $\tilde{g}_{11,02}$ are significantly smaller if the coupler is further detuned, we assume that the plateau stage is the only contributor to the dynamics and the effective Hamiltonian from Eq. <ref> is constant in time. We show later that this approximation is valid for realistic gate durations close to 40 ns and longer. If this approximation cannot be made, the form of the effective Hamiltonian remains unchanged. However, the analytical results derived in the remainder of this section become more complex due to the time dependence of the effective Hamiltonian.

The effective Hamiltonian in Eq. <ref> is block-diagonal. Using matrix exponentiation, one can readily obtain an effective unitary time-evolution operator. Assuming that the computational transmons are detuned by approximately one anharmonicity, i.e. $\tilde{\omega}_{q_1} -\tilde{\omega}_{q_2} \approx \tilde{\alpha}_{q_2}$, the Rabi oscillation frequency between the single-excitation states $|01\rangle$ and $|10\rangle$ is lower than that of the transition between $|11\rangle$ and $|02\rangle$ because $|\tilde{\omega}_{q_1} -\tilde{\omega}_{q_2}|\gg \tilde{g}_{01,10}$. Thus, we can neglect the interaction between the single-excitation states.

Since we are interested in the operator acting on the computational subspace of the system, we further truncate the subspace by excluding the non-computational $|02\rangle$ state.

After accounting for the single-qubit phases, which is typically done via virtual $Z$-rotations [51], the effective time-evolution operator $\hat{U}(t)$ is equivalent to

\begin{align}\label{eq:effective_unitary}
\hat{U}(t) \hat{=}
\begin{blockarray}{ccccc}
|00\rangle & |01\rangle & |10\rangle & |11\rangle & \\
\begin{block}{(cccc)c}
1 & 0 & 0 & 0  &\, |00\rangle \\
0 & 1 &  0 & 0 &\, |01\rangle \\
0 & 0 & 1 & 0 &\, |10\rangle \\
0 & 0 & 0  & \sqrt{P_{11}} e^{i \phi(t)} &\, |11\rangle \\
\end{block}
\end{blockarray}.
\end{align}

Here, we have defined the conditional phase $\phi(t)$ and the population of the $|11\rangle$ state $P_{11}$. Note that this operator is not necessarily trace-preserving, as part of the population of the $|11\rangle$ state might remain in the $|02\rangle$ state, due to potential calibration errors.

The equation for $P_{11}$ is derived by recognizing the block-diagonal structure of the effective Hamiltonian in Eq. <ref>. This simplifies the analysis, reducing it to a standard two-level Rabi oscillation, with $P_{11}$ given by[52]

\begin{equation}\label{eq:perturbative_P11}
P_{11}(t) = 1 - \frac{2 \tilde{g}_{11,02}^2}{\Omega^2}\left[ 1 - \cos(\Omega t)\right],
\end{equation}

where we have defined the qubit-qubit detuning and the Rabi frequency of the transition as $|11\rangle \leftrightarrow |02\rangle$ as $\tilde{\Delta} = \tilde{\omega}_{q_1} - \tilde{\omega}_{q_2} $ and $\Omega = \sqrt{(\tilde{\Delta} - \tilde{\alpha}_{q_2})^2 + 4\tilde{g}_{11,02}^2 }$, respectively.

Calculating the conditional phase is less straightforward since the accumulated phase of the $|11\rangle$ state needs to have the single-excitation phases subtracted from it. Nevertheless, we can readily obtain

\begin{align}\label{eq:perturbative_cphase}
\phi(t) &= \frac{1}{2}\left[ (\tilde{\alpha}_{q_2} - \tilde{\Delta})t + \pi \left(1 - \mathrm{sign}\{\cos(\Omega t/2) \} \right) \right] +  \mathrm{arctan}\left( \frac{\tilde{\Delta} - \tilde{\alpha}_{q_2}}{\Omega} \tan(\Omega t/2)\right).
\end{align}

Eqs. <ref> and <ref> can be used to asses the susceptibility of the gate parameters to a small perturbation, such as a parity switch. Due to the larger charger dispersion of the second-excited state, and the fact that the second-excited state of Qubit 2 ($q_2$) is populated during the gate operation, we can assume that the main contribution of the parity switch is the perturbation of the anharmonicity $\alpha_{q_2}$.

By treating the parity-dependent contribution to the anharmonicity $\delta\alpha$ from Eq. <ref> as a small perturbation, we can obtain the parity-dependent expressions for the conditional phase and $|11\rangle$ population

\begin{align}
\phi(t_g,P_{q_2}) &\approx \phi_0 + \frac{\partial \phi}{\partial \alpha_{q_2}}\bigg|_{t = t_g} \delta\alpha_{q_2}(P_{q_2},n_g), \label{eq:cphase_taylor}\\
P_{11}(t_g,P_{q_2}) &\approx 1 + \frac{1}{2}\frac{\partial^2 P_{11} }{\partial \alpha_{q_2}^2}\bigg|_{t = t_g} \left[ \delta\alpha_{q_2}(P_{q_2},n_g) \right]^2,\label{eq:Pee_taylor}
\end{align}

where the Taylor expansion of the optimal gate parameters for a small perturbation of $\alpha_{q_2}$ evaluated at the mean (parity averaged) anharmonicity $\alpha_{q_2}$ from Eq. <ref> was employed. The Taylor expansion performed here assumes that both parity states suffer from an error of the same magnitude (e.g., the error of the conditional phase for both parities has the same absolute value) and we show in the remainder of the text that this assumption results in higher process fidelities of the gate. At this point, we stress again that $\delta\alpha(P,n_g) = P\epsilon_2\cos(2\pi n_g)/2$.

While the above expressions are completely general also in the non-perturbative regime, the relations given in Eqs. <ref> and <ref> can be used to obtain analytical expressions for ${\partial \phi}/{\partial \alpha_{q_2}}$ and ${\partial^2 P_{11} }/{\partial \alpha_{q_2}^2}$.

We realize from Eq. <ref> that the implementation of a high-fidelity gate with an arbitrary conditional phase $\phi_0$ requires $P_{11}(t_g)=1$, otherwise some population remains outside of the computational subspace. From this observation we readily arrive to the condition for the gate time $\Omega t_g = n\cdot 2\pi$, in which $n \in \mathbb{N}$ is an integer number. This condition enables us to further simplify the relations for the susceptibility of the conditional phase to a parity switch, and we obtain up to the leading order,

\begin{equation}\label{eq:simplified_pert_cphase}
\frac{\partial \phi}{\partial \alpha_{q_2}}\bigg|_{t = t_g} \approx  \frac{t_g}{2\hbar},
\end{equation}

\begin{equation}\label{eq:simplified_pert_P11}
\frac{\partial^2 P_{11} }{\partial \alpha_{q_2}^2}\bigg|_{t = t_g} \sim \mathcal{O}(g_{q_{1}c}^2 g_{q_{2}c}^2 ) +\mathcal{O}(g_{q_{1}q_2} g_{q_{1}c} g_{q_{2}c} )  .
\end{equation}

In the derivation of Eq. <ref> we have neglected the terms proportional to $g_{q_1 c}g_{q_2 c}/(\Delta_{q_2 c} + \alpha_{q_2})^3$ and $g^2_{q_2 c}/(\Delta_{q_2 c} + \alpha_{q_2})^3$ and higher orders. Additionally Eq. <ref> only contains the scaling of the result.

To summarize, the derived relations (equations <ref> and <ref>) are therefore valid under the following assumptions:

* The initial assumptions used to derive the effective Hamiltonian in Eq. <ref> are valid, meaning that:

* The second-order perturbation theory used for the Schrieffer-Wolff transformation is valid, i.e. $g_{q_2c}^2/(\omega_{q_2} - \omega_c)^2 \ll 1$.

* The gate has low leakage outside of the considered subspace.

* The rotating wave approximation for the counter-rotating coupling terms is justified.

* The Rabi frequency associated with the $|01\rangle \leftrightarrow |10\rangle$ transition is much longer compared to the Rabi frequency of the $|11\rangle \leftrightarrow |02\rangle$ transition. In the perturbative regime this is fulfilled if $\sqrt{\tilde{\Delta}^2 + 4\tilde{g}_{01,10}^2} \ll \Omega$.

* The coupler frequency in the interaction regime is relatively constant. If this is not the case, the time dependency of the perturbative parameters must be taken into account.

* The gate has low leakage outside of the computational subspace.

* The simplified formula in Eq. <ref> additionally neglects the terms proportional to  $g_{q_1 c}g_{q_2 c}/(\Delta_{q_2 c} + \alpha_{q_2})^3$ and $g^2_{q_2 c}/(\Delta_{q_2 c} + \alpha_{q_2})^3$, and smaller.

* We have assumed that the main contributor to the perturbation is the Qubit 2 ($q_2$) whose second-excited state is populated during the gate operation. However, if the system is designed in such a way that the charge dispersion of any of the other two transmons is significantly larger, their effects might not be negligible anymore.

Note that both assumptions about the leakage are automatically fulfilled if the gate has a high-fidelity.

§.§ Gate Fidelity Limitations

So far we have quantified how a parity switching event can affect the parameters of the gate unitary. In order to describe the gate performance in a quantum circuit, we also need to consider how frequent are parity switching events.

§.§.§ Parity-Switching Lifetimes

Experiments of parity switching lifetimes typically show the parity switching time to lie in the broad range of $T_P \sim 100 \,\mu \mathrm{s}\,-\, 1 \,\mathrm{s}$ [53, 54, 55, 47, 56, 57]. Even though the parity switching lifetime of transmons might increase in the future, e.g., due to better design and improved shielding, it appears that the parity lifetime may be fundamentally upper bounded by the high-energy quasiparticle burst events, which are observed to happen once every $10-50\,$s[58, 59, 60, 61].

§.§.§ Kraus Operator Description

Comparing the realistic range of parity lifetimes of superconducting qubits to the duration of a single two-qubit gate $t_g$, which is typically in the range of tens to hundreds of nanoseconds [37, 40, 32, 33, 48], we observe that $t_g \ll T_P$. However, any meaningful application of a quantum computer will include re-running an algorithm, comprised of a large number of non-parallel gates $N_\mathrm{gates} $, in order to reduce the statistical uncertainty of the observable being evaluated. Denoting the number of runs by  $N_\mathrm{shots}$, the uncertainty of any observable generally scales as $1/\sqrt{N_\mathrm{shots}}$[52] and therefore $N_\mathrm{shots}$ must be large. All together, the time needed to execute a full algorithm $T_\mathrm{alg}$ roughly scales as $T_\mathrm{alg} \propto N_\mathrm{shots} \cdot N_\mathrm{gates} \cdot t_g$, but it may be realistically even longer, due to the time needed to measure and reinitialize the quantum computer, and any possible pulse schedule compilation of the control electronics[62].

As an example, in Refs. [7, 9] the state of the art devices were run for a total of several minutes in order to obtain meaningful results, which implies the following clear separation of timescales

\begin{equation}\label{eq:separation_of_timescales}
t_g \ll T_P \ll T_\mathrm{alg}.
\end{equation}

The left side of Eq. <ref> indicates that the probability for a parity switch occurring during the operation of a single gate is very low, while the right-hand side suggests that a large number of parity switches can occur during an execution of an algorithm. This means that the effect of the charge-parity switch (CPS) on a diabatic CPHASE gate can be described by the following Kraus operators[63] acting on the two-qubit density matrix $\hat{\rho}$

\begin{equation}\label{eq:kraus_operators}
\mathrm{CPHASE}_\mathrm{CPS}[\hat{\rho}] =\hat{U}_{-} \hat{\rho} \hat{U}^\dagger_{-} + \hat{U}_{ +} \hat{\rho} \hat{U}^\dagger_{ +},
\end{equation}

where $\hat{U}_{ \pm}$ are Kraus operators corresponding to the different parity implementations of the two-qubit gate. Eq. <ref> can be interpreted as a stochastic application of two different gate operators; by assuming that the target conditional phase is $\phi_0$, they can be written as

\begin{equation}\label{eq:kraus_op_definition}
\hat{U}_\pm(t) = \frac{1}{\sqrt{2}}
\begin{pmatrix}
1 & 0 & 0 & 0 \\
0 & 1 & 0 & 0 \\
0 & 0 & 1 & 0 \\
0 & 0 & 0 & \sqrt{1 - \frac{\delta P_{11}}{4}} \,e^{ i\phi_0  \pm i \frac{\delta \phi}{2}}\\
\end{pmatrix},
\end{equation}

with $\delta \phi = \partial\phi/\partial \alpha_{q_2}\,\epsilon_2^{q_2}\cos(2\pi n_g^{q_2})$ and $\delta P_{11} = \partial^2 P_{11}/\partial \alpha_{q_2}^2 \,\left[ \epsilon_2^{q_2} \cos(2\pi n_g^{q_2}) \right]^2/2$, which is a result of Eqs. <ref> and <ref>, together with $\delta\alpha(P,n_g) = P\epsilon_2\cos(2\pi n_g)/2$. We have additionally assumed that both parities are equally likely, however this assumption can also be easily relaxed.

We have defined the Kraus operators in Eq. <ref> so that there are small errors associated with each parity state. We show in the following that this corresponds to higher process fidelities, compared to having one parity state with a perfect fidelity and the second parity state with a larger error. This assumption is, therefore, equivalent to optimally calibrating the gate with respect to the parity switching error.

§.§.§ Gate Fidelity

Using the Kraus operator description in Eq. <ref> allows us to make a formal analysis of performance metric of an arbitrary conditional phase gate in the presence of parity switches and resulting phase and leakage errors.

The process fidelity $\mathcal{F}$ can be computed by [64, 65]

\begin{equation}\label{eq:fidelity_definition}
\mathcal{F}[\mathrm{CPHASE}_\mathrm{CPS}] = \frac{ \frac{1}{d}\sum_{i = +,-} \left|\mathrm{tr}\left\{\hat{U}_\mathrm{CPHASE}^\dagger\hat{U}_i\right\}\right|^2 + 1 - L}{d + 1},
\end{equation}

with the leakage parameter $L = 1 -  \mathrm{tr}\left\{\hat{U}_\mathrm{CPHASE}^\dagger \left[ \sum_{i = +,-} \hat{U}_i \hat{U}_i^\dagger\right] \hat{U}_\mathrm{CPHASE}\right\}/d$.
Here $\hat{U}_\mathrm{CPHASE}$ is the unitary operator of an ideal CPHASE gate and $d$ is the dimension of the computational Hilbert space, which in our case is $d = 4$.

The infidelity $1-\mathcal{F}$ of the operation given in Eq. <ref> for small perturbations $\delta \phi$ and $\delta P_{11}$ is given by the following expression

\begin{align}
\mathcal{F} \approx 1 - \frac{3}{80} \left[\frac{\partial \phi}{\partial \alpha_{q_2}}\bigg|_{t = t_g} \epsilon_2^{q_2} \cos(2\pi n_g^{q_2}) \right]^2 + \frac{1}{32} \frac{\partial^2 P_{11} }{\partial \alpha_{q_2}^2}\bigg|_{t = t_g} \left[ \epsilon_2^{q_2} \cos(2\pi n_g^{q_2}) \right]^2. \label{eq:cphase_Pee_infidelity}
\end{align}

Here, we observe that the effect on the fidelity is of second-order in the charge dispersion, due to the coherent nature of a conditional phase error. However, we have shown in Eq. <ref> of the Methods section that the infidelity of a series of $N$ gates is given by $1 -\mathcal{F}_N \approx  N(1 - \mathcal{F})$. This indicates that the error scales linearly with the number of gates, as is typical of incoherent errors, but quadratically in the error parameter $\delta\phi$, as is expected from a coherent error[66, 67].

In deriving Eq. <ref>, one can observe that the numerical prefactors in front of both terms (here $3/80$ and $1/32$) increase, and therefore the process fidelity decreases, if the gate is not calibrated in such a way that the error is equally distributed between both parity states. For this reason, we have chosen the Kraus operators as given in Eq. $\hat{U}_\pm$.

We conclude from Eq. <ref>, together with Eqs. <ref> and <ref>, that the dominant effect of the parity switch event is the shift in the conditional phase, rather than leakage, since $({\partial \phi}/{\partial \alpha_{q_2}})^2 \gg {\partial^2 P_{11} }/{\partial \alpha_{q_2}^2}$. Moreover the magnitude of the shift in conditional phase is given by Eq. <ref> and Eq. <ref>. Further simplifying Eq. <ref> therefore results in

\begin{equation}\label{eq:simplified_infidelity}
1 - \mathcal{F} \approx \frac{3}{320} \left[ \epsilon_2^{q_2}t_g \cos(2\pi n_g^{q_2})/\hbar \right]^2.
\end{equation}

§.§ Numerical simulations

Here, we compare the above results with the numerically exact treatment of full the Hamiltonian in Eq. <ref>. The parity switch effect is taken into account in the numerical experiments by considering the transmon Hamiltonian from Eq. <ref>. In other words, the system is simulated for all $2^3 = 8$ possible parity states and the results are averaged accordingly.

According to Eq. <ref>, the magnitude of the effect also depends on the offset voltage of Qubit 2, which is typically not known. Thus, we assume for simplicity that $n_g^{q_2} = 0$. Alternatively, as long as all the effects remain second order in the charge dispersion, one can also define an average charge dispersion as $\bar{\epsilon}_2^{q_2} = \epsilon_2^{q_2} \sqrt{ \int_0^1 \mathrm{d}n_g \cos^2(2\pi n_g) } = \epsilon_2^{q_2}/\sqrt{2} $.

§.§.§ Parity-Switch Induced Gate Errors

To achieve high-fidelity gate simulations, it is crucial to carefully select the Hamiltonian parameters in Eq. <ref>. We have described how these parameters were chosen so that high-fidelity gates are possible with arbitrary ratios of $E_J/E_C$ of Qubit 2 in Table <ref> of the Methods section.
We also note that in general, we use parameters that closely resemble those in the implementation presented in Ref. [32].

Effects of a parity switch on the fidelity of a CZ gate. (a) Analytical predictions (red circles) of the conditional-phase difference of the CZ gate between the two parities. The theory predictions from Eqs. <ref> and <ref> are compared to the phase extracted from a full numerical simulation of the gate (blue crosses), with the Hamiltonian parameters from Table <ref>. The inset displays the relative absolute error between the analytical predictions and the numerical results. The $x$-axis represents the $E_J/E_C$ ratio (bottom) and second excited state charge dispersion (top) of the qubit with the higher frequency (Qubit 2). The inset displays the relative error of the analytical prediction. (b) Gate infidelity after a parity switch on Qubit 2, where the theory predictions (red circles) only take into account the effect on the conditional phase, according to Eq. <ref>. The green lines are estimates of the upper bound for the infidelity contribution of quasiparticle induced decoherence on the two-qubit gate system, with the parity lifetimes $T_P$ measured by different references: [1] Kurter, C., et al. (2022)[56], Diamond, S., et al. (2022)[53] and [3] Risté, D., et al. (2013)[54]. The grey region is the region in which errors due to unwanted transitions during the gate operation are more prominent and the effect of a parity switch is negligible. It therefore represents the lower bound of the infidelity, which is largely independent of the ratio $E_J/E_C$.

Here, we examine the effect of a parity-switch on the fidelity of a CZ gate and compare the perturbative results to a full simulation of the Hamiltonian in Eq. <ref>. The full numerical treatment includes also the ramping up and down of the pulse, as described in Eq. <ref>. The pulse parameters in Eq.  <ref> are set to $\tau_b = 2\sqrt{2}\sigma $ with $\sigma = 5\,$ns, and the amplitude of the pulse $A$ and $\tau_c$ are optimized such that the gate fidelity is maximal. More details about the numerical simulations are given in Ref. [36]. The parity of the whole system is then switched from one state to the other according to charge dispersion given by Eqs. <ref> and <ref>.

Since in our analytical approach we approximate the realistic pulse shape shown in Eq. <ref> with a square pulse (i.e., we neglect any dynamics during the ramping up and down of the flux pulse), we have to be careful when comparing the results to those obtained with the numerical data. Particularly, the pulse duration $T$ in the numerics is in general different than the gate duration $t_{g}$ we have defined in our analytical derivation, since the full duration $T$ also includes the time neede to ramp the pulse up and down. When we compare our analytical results with numerical simulations, we first extract the number $n$ of Rabi oscillations from the simulated data by monitoring the population of the $|002\rangle$ state during the gate operation. After $n$ is determined, the analytical gate time is adjusted such that $t_g = n\pi/\Omega$. The effective gate time obtained in this manner does not differ significantly from the duration of the flat part of the pulse, typically less than 5 ns.

The conditional phase of the gate in the simulations is obtained from propagating the state $|\psi(t= 0)\rangle = \frac12 (1,1,1,1)^\mathrm{T}$ and extracting the conditional phase of the $|11\rangle$ state. Fig. <ref>(a) compares the analytical results in Eqs. <ref> and <ref> and to the numerically obtained values and confirms that the two approaches agree up to a good accuracy. Since the magnitude of the charge dispersion in the numerical analysis is the same as in the analytical treatment, the error stems completely from the approximations made in evaluating $\partial \phi / \partial \alpha_{q_2}$. Furthermore, Fig. <ref>(a) also clearly shows that the parity-switching-induced shift in the conditional phase scales exponentially with the ratio of $E_J/E_C$. This, in turn, is due to the scaling of the charge dispersion in Eq. <ref>.

In Fig. <ref>(b), we show numerical data of the full gate fidelity after the parity switch and, as a comparison, the corresponding analytical result calculated using Eq. <ref>. The process fidelity in the numerical example is obtained from propagating a number of input states, reconstructing the effective superoperator of the gate from these simulations and subsequently using Eq. <ref> to obtain the process fidelity of the gate. We observe that the infidelity of the numerical simulation flattens for $E_J/E_C\gtrsim 75$, which is due to other errors in the gate implementation, such as leakage transitions during the pulse ramping up and down. In this region (in grey) the effect of a parity switch is not seen since it is simply too small compared to other errors. On the other hand at lower values of $E_J/E_C$, the numerical results overlap with the parity switching error predicted by the shift in the conditional phase. While clearly demonstrating the magnitude of the error, this result also shows that the leakage error contribution in Eq. <ref> is negligible in the perturbative regime. The gate durations in Fig. <ref> are typically 45 ns $ \lesssim t_g  \lesssim $ 60 ns, with $\tau_c \sim t_g$.

Additionally, Fig. <ref>(b) compares the magnitude of the error due to quasiparticle related decoherence to the parity switch induced error described in this work. Since the quasiparticle induced characteristic decay times $T_1$ and $T_\phi$ depend on a large number of parameters and there are two possible quasiparticle generating mechanisms [see Fig. <ref>(b)] [16], we only provide an upper bound based on the parity switching time observed in the references cited in the caption of Fig. <ref>. This upper bound is determined by noting that in the computational subspace[55]

\begin{equation}
\Gamma^{+-}_{00} + \Gamma^{+-}_{11} + \Gamma^{+-}_{01} + \Gamma^{+-}_{10} \approx 2/T_P,
\end{equation}

where the rates $\Gamma_{ij}^{+-}$ represent the transition rates between states $|i^\pm\rangle \rightarrow |j^\mp \rangle$ in different parity manifolds. Since amplitude damping noise has a larger effect on the infidelity compared to pure dephasing (see Table <ref>), we furthermore assume the worst-case scenario in which each quasiparticle parity switching event produces a $T_1$ decay, so that $1/T_1^{qp} = \Gamma^{+-}_{01} + \Gamma^{+-}_{10} \approx 2/T_P$. This expression for $T_1^{qp}$, together with experimentally measured values of $T_P$, provides an approximate upper bound for the magnitude of the effect of the decoherence. Even though each green line in Fig. <ref>(b) uses a constant measured $T_P$, independent of $E_J/E_C$, the plotted infidelity contribution is not constant due to the varying gate duration $T = 2\tau_b + \tau_c$.

Very interestingly we observe from Fig. <ref>(b) that depending on the ratio of $E_J/E_C$ and the parity lifetimes $T_P$, the contribution of parity switching to the infidelity of the two-qubit gate system can dominate the contribution from quasiparticle-induced relaxation. Note that the two-qubit gates are also typically the noisiest building blocks of a quantum algorithm [66, 7, 8, 11]. We have also provided additional numerical results in the Supplementary Information, showing that the effect of a parity switch on the leakage becomes the main contribution to the infidelity at shorter gate times.

We have additionally included a comparison of the infidelity of a two-qubit gate due to $1/f$-type charge noise in the Supplementary Information in order to compare the two error sources whose contributions scale with the charge dispersion of the transmon. However, unlike the charge dispersion of the transmon which is maximized when $\cos(2\pi n_g) = 1$, the low-frequency charge noise decoherence rate is maximal at $\sin(2\pi n_g) = 1$ . This means that these two errors are mutually exclusive, i.e. if we were hypothetically able to tune $n_g$ to the value where the charge dispersion of the transmon is equal to zero, that point corresponds to the maximal decay rate due to low-frequency charge noise[12]. Nonetheless, in the comparison we have assumed in both the charge-parity switching analysis and the low-frequency charge noise analysis that we are at the noise hotspot, i.e. the value of $n_g$ where the effects are maximal, thus slightly overestimating both errors. We observe that at lower $E_J/E_C$ the charge-parity switching error is dominant and vice versa. The crossover between the infidelity due to the charge-noise-induced dephasing in the computational subspace and the parity switching exhibits a crossover at $E_J/E_C \sim 80$, which corresponds to a gate infidelity of $1-\mathcal{F} \sim 10^{-8}$, i.e. the charge-parity switching error is dominant for $E_J/E_C \gtrsim 80$. However, due to the increased charge dispersion in the higher-excited levels of the transmon and the utilization of the second-excited state, we have additionally analyzed the effect of the second-excited state charge-noise-induced dephasing which was found to be much larger compared to the effects in the computational subspace. In this case, the crossover between the charge-parity error, which was again dominant at lower $E_J/E_C$, was found to occur at $E_J/E_C \sim 65$ and infidelities on the order of $1-\mathcal{F}\sim 10^{-6}$. We note here that our analysis overestimates the effect of the charge noise in a realistic scenario by neglecting the time correlations in the noise and should be treated as an upper bound of the magnitude of the effect. These results demonstrate that the parity-switching error is the dominant charge dispersion related error of the two-qubit gate in the regime of $E_J/E_C \lesssim 65$.

§.§.§ Idling Interaction

One of the main benefits of the tunable coupler architecture is the ability to also suppress any residual interactions between the qubits during idling[43, 34, 35], with unwanted ZZ coupling strengths demonstrated to be below 1 kHz [32]. Here, we show that an uncontrolled parity switching event sets a lower bound on the minimum achievable unwanted ZZ interaction in such systems.

While the conditional phase during the gate operation is obtained by transferring the population outside of the computational subspace, the residual interactions during the idling period are a consequence of dispersive shifts of the energy levels of the computational states, as described in Eq. <ref>.
However, in the presence of parity switches in all three transmons of the tunable-coupler setup, the ZZ-coupling rate can be tuned to exactly zero only for one of the eight equally likely parities. By distributing the errors equally across all parities and assuming the charge dispersions are small, the idling coupling strength (depending on the three individual parities) can be written as

\begin{equation}\label{eq:parity_zz_coupling}
\tilde{\zeta}_\mathrm{ZZ}(P_{q_1},P_c,P_{q_2} ) \simeq \sum_{i\in\{q_1 ,c, q_2 \}} \frac{P_i}{2} \left(\frac{\partial \zeta_\mathrm{ZZ}}{\partial \alpha_i} \epsilon_2^i +\frac{\partial \zeta_\mathrm{ZZ}}{\partial \omega_i} \epsilon_1^i \right),
\end{equation}

with $P_i \in \{-1,+1\} \, \forall i$, and the derivatives are evaluated at the idling point. While we have still assumed $\epsilon_2 \gg \epsilon_1$, we have also acknowledged that the first-excited-state charge dispersion and the derivatives ${\partial \zeta_\mathrm{ZZ}}/{\partial \omega_i}$ can in some cases be significantly larger compared to ${\partial \zeta_\mathrm{ZZ}}/{\partial \alpha_i}$. Unlike the gate fidelity results, the above expression is first order in the charge dispersion.

However, in order to evaluate the derivatives in Eq. <ref>, or the general value of the coupling rate $\zeta_\mathrm{ZZ}$ with fixed Hamiltonian parameters (without including parity effects), it is necessary to go to fourth-order perturbation theory in the coupling strengths, as was done in Refs. [32, 35]. While the complete fourth-order expression is impractical for obtaining any analytical insight, by further assuming the hierarchy of the system parameters $\Sigma_{ic}\gg \Delta_{ic}\gg g_{ic} \gg g_{12}$, where $\Sigma_{ij} = \omega_i + \omega_j$ and $\Delta_{ij} = \omega_i - \omega_j$, for $i,j\in\{q_1,c,q_2 \}$, the cumbersome expressions can be significantly simplified into[35]

\begin{align}\label{eq:perturbative_zz_coupling}
\zeta_\mathrm{ZZ} &\approx \frac{2\left[(\alpha_{q_1} + \alpha_{q_2})\tilde{g}_{01,10}^2 - 2\nu\tilde{g}_{01,10} (2\alpha_{q_1}\alpha_{q_2} + (\alpha_{q_1} - \alpha_{q_2})\Delta_{q_1 q_2}) \right]}{(\Delta_{q_1 q_2} + \alpha_{q_1})(\Delta_{q_1 q_2} - \alpha_{q_2})} + 2\nu^2\left[ 4\alpha_c + \frac{(\alpha_{q_1} + \alpha_{q_2})\Delta_{q_1 q_2}^2}{(\Delta_{q_1 q_2} + \alpha_{q_1}) (\Delta_{q_1 q_2} - \alpha_{q_2})}\right],
\end{align}

where $\nu = g_{q_1 c}g_{q_2 c}/(2\Delta_{q_1 c}\Delta_{q_2 c})\sim 10^{-3}$.

Comparison of the analytical formulas with numerical results obtained via exact diagonalization in the idling regime. (a) Largest derivatives in Eq. <ref> obtained numerically (solid lines) and from Eq. <ref> (filled circles) using the parameters from Table <ref> and $\alpha_{q_2} = -270\cdot 2\pi$ GHz. The $x$-axis represents the coupler frequency, with the two idling frequencies denoted with black dashed vertical lines. The filled circles correspond to the derivatives obtained using Eq. <ref>. (b) Comparison of numerical and analytical results for the root-mean-squared coupling strength at the idling point (defined in Eq. <ref>), with the bar graph displaying the numerical values of the coupling for the 8 different parity states for the highlighted data point. The parameters of the simulation are chosen identically as in Fig. <ref> and are listed in the Methods section, more specifically in Table <ref>. The parameters for the bar plot correspond to the parameters in panel (a).

The difference between the derivatives obtained from Eq. <ref> and the numerically exact result is shown in Fig. <ref>(a) from which we observe that a qualitatively good agreement can be found in the vicinity of the two idling frequencies of the system (vertical black dotted lines). We therefore evaluate the error between the analytical formulas and the numerically exact expression in more detail below.

Since the system has eight uncontrolled, rapidly (compared to experimental timescales, as in Eq. <ref>) switching parity states, we define a parity-averaged idling interaction strength by averaging over all parities:

\begin{equation}\label{eq:averaged_parity_zz}
\langle \tilde{\zeta}^2_\mathrm{ZZ} \rangle_\mathrm{CPS} = \frac{1}{2^3}\sum_{P_{q_1},P_c,P_{q_2} \in \{ -1, +1\}} \left[ \tilde{\zeta}_\mathrm{ZZ}(P_{q_1},P_c,P_{q_2}  ) \right]^2.
\end{equation}

This definition is the parity-averaged mean-square of the idling ZZ strength defined in Eq. <ref>. Parity-averaged coupling strengths are shown in Fig. <ref>(b), where we observe a good agreement between the analytical result and the numerical data. The possible values of the quantity $\zeta_\mathrm{ZZ}(P_{q_1},P_c,P_{q_2}) $ for the third data point [also shown in panel (a)] are visualized in the bar graph. The bar graph demonstrates that Qubit 2 ($q_2$) is the dominant contributor to the perturbation, with small corrections due to the parity of Qubit 1 ($q_1$). This can be seen by noting that the residual coupling strength is dependent mostly on the parity of Qubit 2 ($q_2$). The parity of the coupler in this set of parameters is largely irrelevant due to the high ratio of $E_{J_c}/E_{C_c} \approx 250$, but this is not always the case. For example, considering the gate implementation from Ref. [37], the ratio of $E_{J_c}/E_{C_c}$ for the coupler is much smaller and therefore, the coupler parity in this implementation has a much larger effect on the strength of the residual ZZ coupling.

The above results therefore present a fundamental limit on the magnitude of unwanted interactions that can be achieved in the tunable coupler architecture. However, the overall effect on an algorithm is more complex as it depends on the duration of the execution (since the interaction is always "on"), which in turn depends on the coherence times. We estimate that this error becomes relevant if the coherence times are of the order of $1/ \sqrt{\langle \tilde{\zeta}^2_\mathrm{ZZ} \rangle_\mathrm{CPS}}$, since a significant unwanted conditional phase is accumulated if the algorithm duration $\Delta T$ (which is determined by the coherence times) is long enough such that $\Delta T \cdot\sqrt{\langle \zeta^2_\mathrm{ZZ} \rangle_\mathrm{CPS}} \sim 1$. For $E_J/E_C \sim 50$, the coherence time (or algorithm duration) needed to observe the parity switching induced residual coupling strength is on the order of 1 ms. Note that current fabrication processes are indeed already approaching this value[44, 45].

§.§.§ Optimal Qubit Parameters

Having established the magnitude of the parity-induced error on a two-qubit gate, we have shown that this error can be mitigated by increasing the $E_J/E_C$ ratio of the transmon. However, there are other error sources present in such architectures[36], and while increasing the $E_J/E_C$ ratio will suppress the parity-switching errors, it may also increase the contribution of other possible error sources. Therefore, in order to find better parameters for future transmon-based quantum computers, we must evaluate the contributions of all errors affecting the system. In particular, we estimate optimal regions for the qubit parameters $E_J$ and $E_C$, where the errors contributing to the gate and state preparation infidelities are minimized.

We consider a number of different error sources relevant to superconducting qubits:

* $T_1$ decay due to the coupling to a bath of two-level systems [68, 69, 70, 71, 72, 73, 74].

* $T_\phi$ pure dephasing due to the coupling to magnetic flux noise [75, 76, 77, 68].

* Leakage affecting single-qubit gates due to low anharmonicity [78].

* State preparation errors due to finite-temperature heating effects [13, 14, 15], without the presence of active reset.

* Errors in the two-qubit gate operation due to parity switch effects that are analyzed in this manuscript.

In Methods, we show how the above error sources scale with the transmon Hamiltonian parameters $E_J$ and $E_C$.  We have not included any errors related to the control and calibration of the individual gates, as such errors do not explicitly depend on the qubit parameters and their inclusion, therefore, would not significantly alter the presented results. In other words, we assume that perfect control and calibration are possible. Similarly, measurement errors are present, but have no explicit dependence on $E_J$ and $E_C$.

Scaling analysis for the most relevant error sources in an architecture of flux-tunable transmons.

Parameter scaling      Single-qubit gate infidelity     Two-qubit gate infidelity

Charge noise ($T_1$)      $\Gamma_1 = \frac{1}{T_1} \propto E_C^{3/2}E_J^{1/2}$     $ \frac{1}{3}\Gamma_1 t_\mathrm{SQG} $     $ \frac{2}{5}(\Gamma_1^{q_1} + \Gamma_1^{q_2} ) t_\mathrm{TQG} $

Flux noise ($T_\phi$)     $\Gamma_\phi = \frac{1}{T_\phi} \propto E_C E_J$     $ \frac{1}{6}\Gamma_\phi t_\mathrm{SQG} $      $ \frac{1}{5}(\Gamma_\phi^{q_1} + \Gamma_\phi^{q_2} ) t_\mathrm{TQG} $

Leakage     $P_\mathrm{leak.} \propto E_C^{-\gamma}$, $5 \lesssim \gamma \lesssim 6$      $\frac{1}{3}P_\mathrm{leak.} $     /

Parity Switch     $(\delta\phi)^2 \propto e^{-2\sqrt{8 E_J / E_C}} E_C^{1/2} E_J^{7/2} $     /     $\frac{3}{80}\left(\frac{t_\mathrm{TQG}}{2\hbar}\epsilon^{q_2}_2 \right)^2$

Thermal excitation     $P_{|1\rangle}/P_{|0\rangle} = e^{-\beta (\sqrt{8E_J E_C} - E_C)}$     $P_{|1\rangle}$     /

For simplicity, we again consider the same tunable coupler system as in the previous section, but arranged in a square grid, as pictured in Fig. <ref>(a). In order to avoid frequency crowding issues with such a connectivity[79], the qubits in the array are divided into low and high-frequency transmons. In this case, it is therefore sufficient to consider only a single pair of qubits with different parameters. Moreover, we assume that the qubits are always detuned by approximately one anharmonicity and the anharmonicity of the qubits is similar. We have found that these two conditions are sufficient to be able to perform high-fidelity two-qubit gates, as described in Table <ref>. Since the parameters of the transmons are related by this condition, we can parametrize the whole system in terms of $E_J$ and $E_C$ of the higher-frequency transmon. Some examples of these parameter values that were obtained from experimental demonstrations of the diabatic CZ gate are plotted in Fig. <ref>(b).

In order to quantify the performance of an algorithm execution with a specified pair of parameters $E_J$ and $E_C$ in mind, we define a performance metric $\mathcal{P}$ which we will then maximize. We further define  $\mathcal{P}$ as one minus a weighted sum of the infidelity contributions of all the relevant errors listed in Table <ref>. This sum can be written as

\begin{equation}\label{eq:total_infidelity}
1 - \mathcal{P} = \sum_{i =T_1, T_\phi, \mathrm{parity} } w_{\mathrm{TQG},i} (1 - \mathcal{F}_{\mathrm{TQG},i}) + \sum_{i = q_1,q_2} \sum_{j = T_1, T_\phi, \mathrm{leak.}} w_{\mathrm{SQG},i,j} (1 - \mathcal{F}_{\mathrm{SQG},i,j}) + \sum_{i = q_1,q_2}  w_{\mathrm{SP},i} (1 - \mathcal{F}_{\mathrm{SP},i}),
\end{equation}

where the summation runs across all infidelity contributions, or more explicitly for the simple algorithm pictured in Fig. <ref>(c)

\begin{align}\label{eq:total_infidelity_full}
1 - \mathcal{P} &=  \frac{2}{5}(\Gamma_1^{q_1} + \Gamma_1^{q_2} ) t_\mathrm{TQG} + \frac{1}{5}(\Gamma_\phi^{q_1} + \Gamma_\phi^{q_2} ) t_\mathrm{TQG} + \frac{3}{80}\left(\frac{t_\mathrm{TQG}}{2\hbar}\epsilon^{q_2}_2 \right)^2 \nonumber \\
&+ 2\left[\frac{1}{3}\Gamma_1^{q_1} t_\mathrm{SQG} +  \frac{1}{6}\Gamma_\phi^{q_1} t_\mathrm{SQG} \right] +  \frac{1}{3}P_\mathrm{leak.}^{q_1}+ 2\left[\frac{1}{3}\Gamma_1^{q_2} t_\mathrm{SQG} +  \frac{1}{6}\Gamma_\phi^{q_2} t_\mathrm{SQG} \right] +  \frac{1}{3}P_\mathrm{leak.}^{q_2}\nonumber \\
&+  \left[ \frac{10 ( t_\mathrm{TQG} + 2t_\mathrm{SQG})}{ T_{1,0}} \right] \left(  P_{|1\rangle}^{q_1} + P_{|1\rangle}^{q_2} \right).
\end{align}

The first term in Eq. <ref> and first line in Eq. <ref> correspond to the errors of the two-qubit gate, the second term and line to single-qubit gate errors and the last to the state preparation error.

A similar fidelity approximator was defined in Ref. [7]. We have introduced additional weights in the sum, in order to account for the relative number of single and two-qubit gates, and also to correctly take into account the fact that the error in state preparation occurs only once, while the gate error is significantly amplified after a number of applications of the operation. As the determination of the position of maximal $\mathcal{P}$ within the $(E_C,E_J)$ landscape, using which we find the optimal range for parameters, depends exclusively on the relative values of the weights $w_{i,j}$, we proceed by assuming that all error terms associated with the two-qubit gate are assigned weights of $w_\mathrm{TQG} = 1$. As pictured in Fig. <ref>(c), we analyze a circuit where we perform four single-qubit gates per each two-qubit gate, with half those single-qubit gates being $\pi$ rotations which are more susceptible to leakage. Note that this ratio of single to two-qubit gates arises naturally with the introduction of randomized compiling into the algorithm[66]. The weights $w_{\mathrm{SQG},T_{1,\phi}}$ for the decoherence induced infidelity during a single-qubit gate therefore have a value of $2$ (since a $\pi$ and $\pi/2$ rotation are applied) per qubit, while the leakage error $w_{SQG,\mathrm{leak.}}  = 1$ per each qubit. Lastly, the thermal state preparation error is weighted so that $w_\mathrm{SP} = 10 ( t_\mathrm{TQG} + 2t_\mathrm{SQG}) / T_{1,0}$, since this quantity is the inverse of the approximate number of single and two-qubit gates that we can perform within a specified coherence time $T_{1,0}$ at the reference point. Therefore, the performance metric $\mathcal{P}$ is an approximator for the process fidelity of the simple circuit displayed in Fig. <ref>(c).

(a) Schematic representation of the square grid architecture, with high- (dark blue) and low-frequency (light blue) transmons, connected via tunable couplers (green). (b) Values of $E_J$ and $E_C$ of Qubit 2 ($q_2$) from different experimental implementations of the tunable coupler transmon architecture. Both parameters are extracted from the reported $\omega_{q_2}$ and $\alpha_{q_2}$, and therefore the points are only an approximation[46]. The annotations refer to the following references [1] Collodo, M. C., et al. (2020)[49], [2] Xu, Y., et al. (2020)[41] [3] Sung, Y., et al. (2021)[32], [4] Wu, Y., et al. (2021)[11], [5] Xu, H., et al. (2021)[33] and [6] Google Quantum AI (2022)[8]. The black lines represent constant ratios of $E_J/E_C$ and the red lines correspond to contours of constant qubit frequency $\omega_{q_2}$. (c) Schematic representation of the circuit used to infer the weights in Eqs. <ref> and <ref>, with the state preparation pictured on the left (in purple), single-qubit $\pi$ and $\pi/2$-rotations in green and blue respectively, and the two-qubit CZ gate in orange. $N$ is the number of times the pictured circuit (without the state preparation) is repeated before measurement and therefore an integer determined by the reference coherence time. In our case, we consider $N = \left\lfloor \frac{T_{1,0}}{10 ( t_\mathrm{TQG} + 2t_\mathrm{SQG}) } \right\rfloor$.  (d) The function $1 - \mathcal{P}$ defined in Eq. <ref> plotted for different values of the second (higher frequency) transmon $E_J$ and $E_C$. We consider a single-qubit gate implemented with a Gaussian DRAG pulse with a duration of $t_\mathrm{SQG} = 16\,$ns, a two-qubit gate duration of $t_\mathrm{TQG} = 50\,$ns and three different reference coherence times indicated on top of each panel. For all three cases, the reference $T_{\phi,0} = T_{1,0}$ at $E_J = 12\cdot h\,\mathrm{GHz}$ and $E_C = 0.2\cdot h\,\mathrm{GHz}$ for Qubit 2 ($q_2$) with the parameters of Qubit 1 ($q_1$) given in Table <ref>. All the parameters (transmon parameters and decay times) are scaled accordingly to different values of $E_J$ and $E_C$, and for each qubit individually, as described in Table <ref>. The lighter green contour marks the region in the plot with the lowest values of $1 - \mathcal{P}$, defined by the 10th percentile of the plotted values. The darker green contour is obtained with a density matrix simulation of the circuit from panel (c) with the same errors, but instead of evaluating the function $1-\mathcal{P}$, it is obtained by minimizing the infidelity of the state before measurement.

The values of $1 - \mathcal{P}$ from Eq. <ref> for certain parameters are plotted in Fig. <ref>(d). Due to the fact that the maximum of $\mathcal{P}$ is an optimal solution only for the circuit pictured in Fig. <ref>(c), we have plotted an optimal region of parameters defined by the 10th percentile of the points with the lowest value of $1-\mathcal{P}$. This optimal region is also contrasted with the optimal region obtained with a density matrix simulation of the same circuit, and the good agreement with between the curves shows that $\mathcal{P}$ is a valid performance metric.

The color scale corresponding to the $1 - \mathcal{P}$ has four distinct less-favorable regions, corresponding to the error sources in Table <ref>. For small $E_C$ and, consequently, low anharmonicity $\alpha_{q_2}$, the overall fidelities are low due to relatively large errors caused by leakage during single-qubit gates. For small $E_J$ and large $E_C$, $E_J/E_C$ has a low value and therefore the system experiences relatively large errors arising from parity switching events. If both $E_J$ and $E_C$ are large, the coherence times are short and, therefore, the dominating source of infidelity. Moving perpendicularly to the constant frequency contours towards lower values of $E_C$ and $E_J$, results in lowered frequency and, thus increased errors caused by thermal excitations.

On the other hand, when considering near-future transmons with coherence times close to 0.5 ms[44, 45] in Fig. <ref>(d), we observe that the optimal region is shifted towards larger values of $E_C$, and is even limited by the parity-switching effects in the bottom right corner, thus demonstrating the importance of this effect in future QPU design, when two-qubit gate infidelities surpass the limit of $10^{-3}$ infidelity.

By comparing the data in Fig. <ref>(b) and (d), we observe that the implementation by Xu, H., et al. (2021)[33] is close to the parity-switching induced error region. Since the pure dephasing times of the system were not reported we are not able to asses the relative contribution of the parity switching error to this gate implementation. However, using the results shown in Eqs.
<ref> and <ref>, we observe that the parity-switching-error induced infidelity for the parameters is estimated to be $1 - \mathcal{F} \approx 2\cdot10^{-3}$. As a comparison, the reported $\Gamma_1$ decay times (measured in the idling configuration) in the reference were $\Gamma_1^{q_1} = 20.8\,\mu$s and $\Gamma_1^{q_2} = 28.8\,\mu$s[33]. Together with the reported gate time, these decay times correspond to an infidelity contribution of approximately $6\cdot10^{-4}$, as given by the expression in Table <ref>. These results show that for this specific implementation of the diabatic CZ gate, the described effects of the parity switches are comparable, and possibly even greater than the $T_1$ decay induced infidelity.

The shifts of the optimal region with increasing coherence times seen in Fig. <ref>(d) mean that the qubit parameters must be adapted to the currently achievable coherence times. Fig. <ref>(a) shows how the performance metric $\mathcal{P}$ increases with the achievable coherence properties of the system, provided that the parameters $E_J$ and $E_C$ are adapted to the coherence times. On the other hand, Fig. <ref>(b) displays the value of the function $\mathcal{P}$ without changing $E_J$ and $E_C$. Although the performance metric $\mathcal{P}$ on panels (a) and (b) is initially largely coherence limited, the assumption of fixed parameters $E_J$ and $E_C$ in the simulations of Fig. <ref>(b) shows that simply increasing the coherence times does not necessarily give better fidelities if the effects of other sources of error are not taken into account. In this case, the parity switching error becomes the dominant error source, as seen from the upper panel in Fig. <ref>(b).

The relative contributions to the value of $1 - \mathcal{P}$ (in Eq. <ref>) from the error sources listed in Table <ref> (top) and the value of the performance metric $\mathcal{P}$ as a function of coherence times with assuming $T_{\phi,0} = T_{1,0}$ at the reference point (bottom). Here we considered the same reference point as used in Fig. <ref>(d). Panel (a) shows a well-designed system, where the parameters $E_J$ and $E_C$ are optimally adapted so that $\mathcal{P}$ is maximized for each value of the coherence time on the $x$-axis. Panel (b) presents $\mathcal{P}$ at various coherence times while keeping $E_J$ and $E_C$ fixed, corresponding to the values marked as point [5] in Fig. <ref>(b).

It is also important to mention what other aspects of successfully operating a transmon tunable coupler based quantum computer were not included in the presented analysis.

* We are neglecting any possible cross-talk effects between next-nearest neighbors.

* Idling errors were discussed in this work, but not included in this analysis, as their contribution depends significantly on the specific algorithm being implemented. Taking also this effect into account, the area with low $E_J/E_C$ is even less favorable.

* Other sources of decoherence are expected to have smaller contributions and therefore do not significantly affect the findings presented here.

* The TLS environment is random, meaning that the scaling shown in Table <ref> only holds for the average of a large number of qubits. Additionally, the $T_1$ decay rate is also heavily influenced by the design, i.e., the exact geometry of the capacitor pads of the transmon, meaning that the presented results are only valid for a comparison of qubits with similar designs.

* Recently, two-qubit gates with tunable couplers idled below the qubit frequencies have been demonstrated[37, 38]. In such implementations, the presence of thermal coupler excitations limit the gate performance at lower qubit frequencies in a more complex manner. Our analysis in the low frequency regime therefore only holds for the original tunable coupler proposals, where the coupler is idled above the computational transmons.

§ DISCUSSION

We have provided a novel framework for the optimization of circuit parameters that can be used to guide the future design of transmon-based quantum computers. Our findings reveal the presence of a distinct global performance peak within the $E_J$ and $E_C$ parameter space, which has not been identified before. Moreover, our optimization procedure can be straightforwardly extended to more error sources, provided that the scaling of the error as a function of the system parameters is known analytically, or the infidelity contribution can be evaluated numerically. The latter is typically realistic as long as the error is sufficiently local, i.e., it depends only on the parameters of a handful of transmons at most. While we have based our analysis on transmon qubits connected via tunable couplers, the same principles, albeit with different error sources, can be applied to different types of qubits[2] or co-design chips[80]. Additionally, more parameters than just $E_J$ and $E_C$ can be optimized, e.g., also the gate durations can be considered as free parameters since they are realistically easy to adjust in experiments. The limiting factor here are the errors for which the analytical behavior is unknown and numerical interpolation in a large parameter space is too demanding.

Additionally, we have established how parity switching affects the commonly implemented tunable-coupler mediated diabatic CPHASE gate in a transmon based quantum computer, both analytically and numerically. We have shown that the parity switching error can be the main quasiparticle-related error source of the two-qubit gate. Moreover, we have demonstrated that the experimental implementation of the gate presented in Ref.[33] may have a comparable, if not larger, contribution of parity-switching errors compared to all $T_1$ decay mechanisms.
While the tunable-coupler-based diabatic CPHASE gate is more relevant due to its implementation in leading large-scale experiments[7, 8, 10, 11], we believe that the effects described in this manuscript should be considered in any current or future transmon-based quantum gate which utilizes higher-excited states.

One of the primary anticipated advantages of incorporating tunable couplers into the system is the potential for on-demand complete suppression of $ZZ$-type interactions among the qubits. However, our research in this context has revealed that the stochastic nature of parity switches imposes constraints on this proposition, practically establishing a lower bound on the achievable minimum $ZZ$ coupling strength. The magnitude of this "always on" interaction should be an important consideration when running longer algorithms. More specifically, this effect becomes relevant if the algorithm is long enough to accumulate a considerable conditional phase due to the unwanted coupling strengths shown in Fig. <ref>(b). Since current coherence times are approaching the 1 ms limit[44, 45], the residual idling strength can become relevant if the described effects are not taken into account in the design of the transmon parameters.

One way of mitigating the parity-switching effects would be to attempt to tune the offset charge $n_g$ to the point where both parity manifolds are degenerate. However, such a solution is not practical, since the environmental charge noise would result in a drift of the offset charge $n_g$ as was demonstrated in Refs. [81, 82]. We further note that $\cos(2\pi n_g) = 0$ is the low-frequency charge noise hotspot [12] in which the qubit frequency is maximally sensitive to the fluctuations of the offset charge. Therefore, the qubit is expected to have lower coherence times at this particular value for $n_g$.

It is important to acknowledge the potential influence of other noise mechanisms, not explicitly addressed in this study, on the optimal design parameters. Consequently, we have chosen to present a range of optimal parameters, rather than prescribing a single optimal value for $E_J$ and $E_C$. Our findings, as depicted in Figure <ref>(d), reveal that a two-fold variation in the reference coherence time only marginally adjusts the optimal parameter domain. This observation underscores the robustness of our results, suggesting that the presence of additional, potentially sub-leading noise mechanisms omitted from our simulations is unlikely to precipitate a drastic alteration in the presented outcomes.

§ METHODS

§.§ The Schrieffer-Wolff Transformation

The Schrieffer-Wolff (SW) transformation used in this manuscript was first introduced in Ref. [83] and a similar transformation has been applied to the computational subspace of the two-qubit system in Refs.[34, 50]. The aim of the SW transformation in our case is not to diagonalize the system, but rather to decouple the coupler states from the computational transmons, thus enabling us to only study the reduced system, i.e. we want to simplify the full Hamiltonian into a more tractable reduced model, containing only the relevant states (i.e. the states which have a significant population).

The reason for this is the fact that the computational basis of the system is defined by the eigenstates of the system which have the maximum overlap with states of the form $|i_{q_1} 0_c j_{q_2}\rangle$, i.e. with the coupler always in the ground state. Any excitations of the coupler therefore lead to errors, so in order to analyze the ideal gate dynamics, we constrain ourselves only to the energy levels of the computational subspace and the second excited state used in the Rabi oscillation during the CPHASE gate operation.

In general, a SW transformation is obtained by noting that any unitary operator can be written as $\hat{U} = e^{\hat{S}} = \mathds{1} + \hat{S} + \frac{1}{2}\hat{S}^2 + \dots$, where $\hat{S}$ is anti-hermitian, $\hat{S} = -\hat{S}^\dagger$. Consequently a unitary transformation of an arbitrary Hamiltonian $\hat{H}$ can be expanded in terms of $\hat{S}$ as

\begin{equation}\label{eq:basic_SW_transform}
\hat{U}\hat{H}\hat{U}^\dagger = \hat{H} + \left[ \hat{S},\hat{H}\right] + \frac{1}{2}\left[\hat{S}, \left[\hat{S},\hat{H}\right]\right] + \mathcal{O}(\hat{S}^3).
\end{equation}

As is typical in perturbation theory, we introduce the index $\alpha$ for bookkeeping purposes, and split the full Hamiltonian into a diagonal part and two off-diagonal perturbations, so that $\hat{H} = \hat{H}_0 +  \alpha\hat{V}_1+  \alpha^2\hat{V}_2$. Additionally we rewrite the operator $\hat{S}$ as a first order operator $\hat{S} \rightarrow \alpha\hat{S}$, since if the perturbation is small, $\hat{U}$ should be close to identity. In our case, $\hat{H}$ is given in Eq. <ref> and the first-order perturbation $\hat{V}_1 = -\sum_{i = q_1,q_2} \hbar g_{ic}(\hat{a}_i^\dagger - \hat{a}_i)(\hat{a}_c^\dagger - \hat{a}_c)$ corresponds to the capacitive couplings of the two transmons (Qubits 1 and 2, $q_{1,2}$) to the coupler ($c$), the direct coupling between the qubits $\hat{V}_2 = - \hbar g_{q_1 q_2}(\hat{a}_{q_1}^\dagger - \hat{a}_{q_1})(\hat{a}_{q_2}^\dagger - \hat{a}_{q_2})$ is a second-order perturbation, while $\hat{H}_0$ is a sum of the three independent anharmonic oscillator Hamiltonians. This hierarchy is chosen due to the fact that in all practical scenarios $g_{q_1 q_2} \ll g_{q_1 c}, g_{q_2 c}$[34].

By plugging the ansaetze into Eq. <ref>, and grouping the terms with the same order of $\alpha$, we obtain

\begin{equation}
\hat{U}\hat{H}\hat{U}^\dagger = \hat{H}_0 + \alpha \left( \hat{V}_1 + \left[ \hat{S},\hat{H}_0  \right]\right) + \alpha^2 \left(\hat{V}_2 + \left[ \hat{S},\hat{V}_1\right] +\frac{1}{2}\left[ \hat{S} ,\left[\hat{S}, \hat{H}_0 \right]\right] \right) + \mathcal{O}(\alpha^3).
\end{equation}

Looking at the first-order term, it is natural to choose $\hat{S}$ such that $\left[ \hat{S},\hat{H}_0\right] = -\hat{V}_1$, i.e. so that we cancel any couplings to the coupler states up to lowest order. However in order to do so and account for the couplings to the higher state correctly, we generalize the transformation from Refs. [34, 50]

\begin{align}
\hat{S}_i &= \sum_{n_{q_i},n_c \in \{0,1\} } \sqrt{(n_{q_i} + 1)(n_c + 1)} \left[ \frac{g_{q_ ic}}{\Delta_{q_i c} + n_{q_i} \alpha_{q_i} - n_c \alpha_c }\left( \hat{\pi}_{q_i}^{n_{q_i}+1,n_{q_i}} \hat{\pi}_{c}^{n_{c},n_{c}+1} - \hat{\pi}_{q_i}^{n_{q_i},n_{q_i}+1} \hat{\pi}_{c}^{n_{c}+1,n_{c}}\right)  \right.  \nonumber \\
&- \left. \frac{g_{q_i c}}{\Sigma_{q_i c} + n_{q_i} \alpha_{q_i}  + n_c \alpha_c}\left( \hat{\pi}_{q_i}^{n_{q_i}+1,n_{q_i}} \hat{\pi}_{c}^{n_{c}+1,n_{c}} - \hat{\pi}_{q_i}^{n_{q_i},n_{q_i}+1} \hat{\pi}_{c}^{n_{c},n_{c}+1}\right)\right] , \\
\hat{S} &= \hat{S}_1 + \hat{S}_2.
\end{align}

We have additionally defined $\Delta_{q_ic} = \omega_{q_i} - \omega_c $, $\Sigma_{q_i c} = \omega_{q_i} + \omega_c $ and the operators $\hat{\pi}_{k}^{n,m} = |n\rangle \langle m |$, acting in the Hilbert space of $k\in \{ q_1,c,q_2\}$.

Since we have assumed the coupler remains in the ground state at all times, the effective Hamiltonian is defined on the set of states $\{|0_{q_1} 0_c 0_{q_2}\rangle,|0_{q_1} 0_c 1_{q_2}\rangle,|1_{q_1} 0_c 0_{q_2}\rangle,|1_{q_1} 0_c 1_{q_2}\rangle,|0_{q_1} 0_c 2_{q_2}\rangle \}$. Additionally, we neglect any couplings outside of the subspace of interest, however the resulting effective Hamiltonian still contains terms coupling the levels $|000\rangle \leftrightarrow |101\rangle$ and $|000\rangle \leftrightarrow |002\rangle$. These couplings are neglected in the rotating-wave approximation, as these transitions do not conserve the total occupation number.

By additionally setting the energy of the ground state to zero we arrive at the effective subspace Hamiltonian from Eq. <ref>. The perturbative parameter values are given by

\begin{align}
\tilde{\omega}_{q_i} &= \omega_{q_i}  + \frac{g_{q_i c}^2}{\Delta_{q_i c}} + \frac{2g_{q_i c}^2}{\Sigma_{q_i c} + \alpha_{q_i} } + \frac{g_{q_i c}^2}{\Sigma_{q_i c}},\label{eq:perturbative_hamiltonian_parameters_omega}\\
\tilde{\alpha}_{q_i} &= \alpha_{q_i} -\frac{2g_{q_i c}^2}{\Delta_{q_i c} } + \frac{2g_{q_i c}^2}{\Delta_{q_i c} + \alpha_{q_i} } +\frac{4g_{q_i c}^2}{\Sigma_{q_i c} + \alpha_{q_i}} + \frac{g_{q_i c}^2}{\Sigma_{q_i c}} - \frac{3g_{q_i c}^2}{\Sigma_{q_i c} + 2\alpha_{q_i}}, \\
\tilde{g}_{01,10} &= g_{q_1 q_2} + \frac{ g_{q_1 c}  g_{q_2 c} }{2} \left(\frac{1}{\Delta_{q_1 c}} + \frac{1}{\Delta_{q_2 c}}
-  \frac{1}{\Sigma_{q_1 c}} - \frac{1}{\Sigma_{q_2 c}} \right), \\
\tilde{g}_{11,02} &= \sqrt{2}\left[ g_{q_1 q_2} + \frac{ g_{q_1 c}  g_{q_2 c} }{2} \left(\frac{1}{\Delta_{q_1 c}} + \frac{1}{\Delta_{q_2 c} + \alpha_{q_2}}
-  \frac{1}{\Sigma_{q_1 c}} - \frac{1}{\Sigma_{q_2 c} + \alpha_{q_2}} \right) \right]. \label{eq:perturbative_hamiltonian_parameters_g1102}
\end{align}

§.§ Fidelity Scaling

Scaling of the fidelity with the number of gates, i.e. computing the process fidelity of sequential application of $N$ gates, can be found by using Eq. <ref> to first define the map corresponding to a series of gates

\begin{equation}\label{eq:N_gate_kraus_map}
\left(\mathrm{CPHASE}_\mathrm{CPS}\right)^N[\hat{\rho}] = \sum_{k=0}^N {{N}\choose{k}} \hat{U}_-^k \hat{U}_+^{(N-k)}\hat{\rho} \, \hat{U}_-^{(N-k)} \hat{U}_+^k,
\end{equation}

where we have used the fact that the operators in Eq. <ref> are diagonal and therefore commute with each other and also $\hat{U}_+^\dagger = \hat{U}_-$.
Using the effective Kraus operators from Eq. <ref>, and only considering the conditional phase error, we combine this with the fidelity definition from Eq. <ref>, and obtain

\begin{align}\label{eq:infidelity_scaling}
\mathcal{F}[\left(\mathrm{CPHASE}_\mathrm{CPS}\right)^N] &=  \frac{d + \sum_{k = 0}^N {{N}\choose{k}}  \left|\mathrm{tr}\left\{\left(\hat{U}_\mathrm{CPHASE}^\dagger\right)^N\hat{U}_-^k \hat{U}_+^{(N-k)}\right\}\right|^2}{d^2 + d}\\
&= \frac{4 + \frac{1}{2^N}\sum_{k = 0}^N {{N}\choose{k}} \left[10 + 6\cos\left( \frac{\delta\phi}{2}(N - 2k)\right)\right]}{20}\\
&\approx 1 - \frac{3 }{80}(\delta\phi)^2 \frac{1}{2^N}\sum_{k=0}^N {{N}\choose{k}} (N - 2k)^2 \\
&= 1 - \frac{3}{80} N (\delta\phi)^2.
\end{align}

This result indicates that calibrating the gate such that the error is evenly split between the two parities not only increases the single-gate fidelity but also leads to a more generous scaling of the infidelity $\propto N$, compared to the purely coherent error case for which the error scales as $\propto N^2$ [66].

§.§ Gate parameters for high-fidelity simulations

Finding good gate parameters, both for the Hamiltonian as well as for the pulse for high-fidelity simulations, is not a trivial task. Here, we discuss how to find optimal Hamiltonian parameters at different qubit frequencies and anharamonicities. In general, the only prerequisite is that the qubits are detuned by approximately one anharmonicity, which can be seen from Eq. <ref>.

Table of Hamiltonian parameters (Eq. <ref>) with which high-fidelity CZ gates are possible.

$\alpha_{q_1}$     $\alpha_{q_2} + 10\cdot h\,\mathrm{MHz}$

$\alpha_{c}$     $-110\cdot h\,\mathrm{MHz}$

$\omega_{q_2}$     $4.8\cdot h\,\mathrm{GHz}$

$\omega_{q_1}$     $\omega_{q_2} + \alpha_{q_2} + 10\cdot h\,\mathrm{MHz}$

$\beta_{q_1 c}$     0.015

$\beta_{q_2 c}$     0.015

$\beta_{q_1 q_2}$     0.001

As seen from Table <ref>, we keep some of the parameters in the simulation fixed, while others depend on $\alpha_{q_2}$, which is varied. These parameters are based on the experimental values from Ref. [32] and the coupling coefficients from the Hamiltonian in Eq. <ref> are related to $\beta_{ij}$ via the following relation $g_{ij} = \beta_{ij}\sqrt{\omega_i \omega_j}$, as described in the main text. The small perturbation to $\alpha_{q_1}$ is there in order to ensure that the energy levels are significantly non-degenerate for perturbation theory to apply. This is also completely realistic as a typical fabrication procedure results in seemingly random deviations from the designed values.
More specifically, the gates simulated in Fig. <ref>(a,b) are obtained by varying $\alpha_{q_2}$ with values $\alpha_{q_2} \in [195,230,250,270,300]\,h\,$MHz. Additionally, in order to certify that the analyzed effect is not limited to the choice of parameters presented above, the second data point of Fig. <ref>(a,b) at $E_J/E_C \approx 50$, was generated in the same way as in Table <ref>, but with the change to $\omega_{q_2} = 5.1\cdot 2\pi\mathrm{GHz}$.
The coupler frequency in the idling configuration $\omega_c^\mathrm{idle}$, i.e. before and after performing a gate, is determined by diagonalizing the Hamiltonian to fulfill condition Eq. <ref>.
The pulse parameters from Eq. <ref> are obtained by numerically optimizing the fidelity of the gate, with fixed $\sigma = 5\,$ns and $\tau_b = 2\sqrt{2}\sigma$. Typical values of the amplitude are $A\sim 1 - 1.2 \cdot \,h \,\mathrm{GHz}$ and $\tau_c \sim 60\,$ns.

§.§ Scaling of Transmon Error Sources

In this section we derive the scaling of the noise parameters shown in Table <ref> with the transmon energies $E_J$ and $E_C$.

The noise is modeled by appending the appropriate noise channel after the gate unitary, and the calculated infidelity is thus independent of the unitary dynamics, as was shown in Ref. [84].

§.§.§ Charge noise $T_1$

The $T_1$ decay time of transmon devices is believed to be currently limited by the presence of a discrete number of environmental two-level systems (TLSs) which couple to the qubit via their electrical dipole [69, 68, 70, 71, 72, 73, 74].
Therefore assuming that the charge noise in the system is weak enough, the interaction Hamiltonian between a transmon and an environmental TLS can be derived from Eq. <ref>, by replacing $n_g \rightarrow n_g + \delta\hat{n}_g$, as in Refs. [69, 74], so that

\begin{equation}
\hat{H}_\mathrm{q-TLS} = 8 E_C \, \hat{n}\otimes \delta\hat{n}_g = -4\sqrt{2} E_C\left( \frac{E_J}{8E_C}\right)^{\frac{1}{4}} i (\hat{a} - \hat{a}^\dagger)\otimes \delta\hat{n}_g.
\end{equation}

In the above equation, we have used the asymptotic expression for the number operator $\hat{n}$, derived already in Ref. [12]. The operator $\delta\hat{n}_g$ is defined in the Hilbert space of the TLS, and is related to the parameters of the TLS - more specifically, its electrical dipole [72]. More importantly, $\delta\hat{n}_g$ does not explicitly depend on $E_J$ and $E_C$.

Since the majority of the $T_1$ experiments on transmons display exponential decays[85], the dynamics can be captured by the Lindblad equation. As a second-order approximation is assumed in the derivation of the Lindblad equation, the resulting decay rates are found to be proportional to the square of the coupling coefficient of the transmon to the TLS environment [63]. In our case this translates to

\begin{equation}
\Gamma_1 \propto E_C^{3/2} E_J^{1/2}.
\end{equation}

Here we have omitted the noise spectrum of the environment, since current models assume a flat noise spectrum without any dependence on the qubit frequency[72].

§.§.§ Flux noise $T_\phi$

While it is not strictly necessary for the computational transmons to be flux-tunable, it is often desired as flux-tunability additionally enables the implementation of an iSWAP gate with the same architecture [32, 34]. Having slow flux tunability is also desirable in order to avoid resonances with TLSs in the environment which can severely limit the $T_1$ decay time [73, 74], as well as helping with the issue of frequency crowding [79].

However, unlike the TLS environment producing $T_1$ dynamics, the noise spectrum of magnetic-flux noise is typically observed to have a $1/f^\alpha$ frequency dependence [75, 76, 77, 68], with $ \alpha \sim 1$. The large noise-spectrum amplitude at lower frequencies means that the long-time correlation results in non-Markovian dynamics [86]. Assuming that the noise is slow enough, in order for the adiabatic approximation to hold, the interaction Hamiltonian due to a slowly fluctuating magnetic environment $\delta\hat{\Phi}$ is given by [12]

\begin{equation}
\hat{H}_\mathrm{q-flux} = \frac{\partial \hat{H}}{\partial \Phi} \otimes \delta\hat{\Phi} = \hbar \frac{\partial \omega}{\partial \Phi} \hat{a}^\dagger\hat{a} \otimes \delta\hat{\Phi}.
\end{equation}

The form of the flux dispersion $\partial \omega / \partial \Phi$ of a split-junction transmon, with Josephson energies $E_{J_1}$ and $E_{J_2}$, is determined by the relations [12]

\begin{align}
E_J(\Phi) &= E_{J\Sigma}\cos(\pi\Phi)\sqrt{1 + d^2\tan^2(\pi\Phi)},\\
\hbar\omega(\Phi) &= \sqrt{8 E_C E_J(\Phi)} - E_C,
\end{align}

with $E_{J\Sigma} = E_{J_1} + E_{J_2}$ and $d = |E_{J_1} - E_{J_2}|/(E_{J_1} + E_{J_2})$.
Identically as in the $T_1$-decay example, we assume that the environment operator $\delta\hat{\Phi}$ does not depend on the transmon parameters. Realistically, the magnitude of this operator depends on the dot product between the magnetic-dipole operator of the spins and the surface vector of the SQUID loop.

Even though the Markovian approximation does not hold anymore, the decay rate due to this noise still scales with the square of the coupling coefficient [87, 88], which brings us to

\begin{equation}
\Gamma_\phi \propto \left(\frac{\partial \omega}{\partial \Phi} \right)^2 \propto E_C E_{J\Sigma}.
\end{equation}

While we have explicitly acknowledged only the first order flux dispersion $\partial\omega/\partial\Phi$, it is easy to see that the above scaling is valid even if $\partial\omega/\partial\Phi = 0$, and we must take into account the second, or any other order dispersion $\partial^n\omega/\partial\Phi^n$.
Additionally, in the main text we assume that $E_J \approx E_{J\Sigma}$, which is a realistic assumption precisely since large deviations from this condition result in an increased sensitivity to flux noise.

§.§.§ Single-Qubit Gate Leakage

The low anharmonicity of the transmon is a limiting factor in the operation of single-qubit gates [78, 36] as fast operations drive a part of the population from the computational subspace into the second excited state. A straightforward and effective scheme for mitigating this effect while implementing single-qubit gates known as Derivative Removal by Adiabatic Gate (DRAG) was presented in Ref. [78].

Since typical pulse amplitudes used to perform single-qubit operations are typically much lower compared to the qubit frequency, and the drive is resonant with the qubit, we assume that the rotating-wave approximation is accurate. In the frame rotating with the qubit frequency, the effective Hamiltonian depends on the pulse parameters, anharmonicity, and the detuning between the drive and qubit frequencies. This indicates that the amount of leakage does not explicitly depend on the qubit frequency.

Thus, the process fidelity, similar to the one defined in Eq. <ref>, but in this case for a single-qubit rotation around the $x$ or $y$-axes, depends only on the qubit anharmonicity (or charging energy), since $\alpha \simeq -E_C$ in the transmon limit. This means that even though analytical results are not available, the relationship can be determined numerically and interpolated. This relationship, albeit with different parameters, has already been plotted in Ref. [78], and generally follows a dependence of $P_\mathrm{leak.} \propto E_C^{-\gamma}$, $5 \lesssim \gamma \lesssim 6$, with higher exponents observed at lower $E_C$. The independence of the single-qubit gate infidelity of the qubit frequency (within the transmon regime) was also verified numerically.

The single-qubit gate parameters assumed in Fig <ref>(d) and Fig <ref> were a DRAG Gaussian pulse, with a $\sigma = 4\,$ns and a total duration $t_\mathrm{SQG} = 4\sigma$. The amplitudes of both DRAG components are numerically optimized before interpolating the dependence of $P_\mathrm{leak.}$ on $\alpha$, which was used to generate Fig. <ref>(d). The pulse drive frequency is assumed to be resonant with the qubit. More details are available in Ref. [36]. We note here, that most of the SQG infidelity is due to leakage, rather than phase errors.

§.§.§ Thermal Excitation Error

While other error sources affect the gate performance, the thermal-excitation error considered here only affects the state preparation.
The process fidelity, which is defined as the fidelity averaged over Haar random distributed input states, is therefore not applicable, since this error only affects one input state. We therefore replace the process infidelity [64] with a state infidelity. This can also be qualitatively thought of as replacing the Haar random distribution with a delta-like distribution with a peak at the $|0\rangle$ state. By modeling the thermal excitation as a bit-flip channel $\mathcal{E}[\hat{\rho}] = (1 - P_{|1\rangle}) \hat{\rho} + P_{|1\rangle} \hat{\sigma}_x\hat{\rho}\hat{\sigma}_x$ it is straightforward to see that the state infidelity[52]

\begin{equation}
\mathcal{F} = \langle 0 |\mathcal{E}[|0\rangle \langle 0 |]  | 0\rangle = 1 - P_{|1\rangle} =  \frac{1}{1 + e^{-\beta \omega}},
\end{equation}

where we have additionally assumed that the temperature is low enough such that the population of the higher-excited states is negligible.

url<#>1urlprefixURL doiprefixDOI:

[1]
authorGyenis, A. et al.
journaltitleMoving beyond the transmon: Noise-protected superconducting quantum circuits.
PRX Quantum volume2, pages030101, <10.1103/PRXQuantum.2.030101> (year2021).

[2]
authorKjaergaard, M. et al.
journaltitleSuperconducting qubits: Current state of play.
Annual Review of Condensed Matter Physics volume11, pages369–395, <10.1146/annurev-conmatphys-031119-050605> (year2020).

[3]
authorHuang, H.-L., authorWu, D., authorFan, D. & authorZhu, X.
journaltitleSuperconducting quantum computing: a review.
Science China Information Sciences volume63, <10.1007/s11432-020-2881-9> (year2020).

[4]
authorMartinis, J. M.
journaltitleSurface loss calculations and design of a superconducting transmon qubit with tapered wiring.
npj Quantum Information volume8, pages26, <10.1038/s41534-022-00530-6> (year2022).

[5]
authorEun, S., authorPark, S. H., authorSeo, K., authorChoi, K. & authorHahn, S.
titleShape optimization of superconducting transmon qubit for low surface dielectric loss (year2022).

[6]
authorMenke, T. et al.
journaltitleAutomated design of superconducting circuits and its application to 4-local couplers.
npj Quantum Information volume7, pages49, <10.1038/s41534-021-00382-6> (year2021).

[7]
authorArute, F. et al.
journaltitleQuantum supremacy using a programmable superconducting processor.
Nature volume574, pages505–510, <10.1038/s41586-019-1666-5> (year2019).

[8]
authorAcharya, R. et al.
journaltitleSuppressing quantum errors by scaling a surface code logical qubit.
Nature volume614, pages676–681, <10.1038/s41586-022-05434-1> (year2023).

[9]
authorKim, Y. et al.
journaltitleEvidence for the utility of quantum computing before fault tolerance.
Nature volume618, pages500–505, <10.1038/s41586-023-06096-3> (year2023).

[10]
authorCao, S. et al.
journaltitleGeneration of genuine entanglement up to 51 superconducting qubits.
Nature <10.1038/s41586-023-06195-1> (year2023).

[11]
authorWu, Y. et al.
journaltitleStrong quantum computational advantage using a superconducting quantum processor.
Phys. Rev. Lett. volume127, pages180501, <10.1103/PhysRevLett.127.180501> (year2021).

[12]
authorKoch, J. et al.
journaltitleCharge-insensitive qubit design derived from the cooper pair box.
Phys. Rev. A volume76, pages042319, <10.1103/PhysRevA.76.042319> (year2007).

[13]
authorWenner, J. et al.
journaltitleExcitation of superconducting qubits from hot nonequilibrium quasiparticles.
Phys. Rev. Lett. volume110, pages150502, <10.1103/PhysRevLett.110.150502> (year2013).

[14]
authorJin, X. Y. et al.
journaltitleThermal and residual excited-state population in a 3d transmon qubit.
Phys. Rev. Lett. volume114, pages240501, <10.1103/PhysRevLett.114.240501> (year2015).

[15]
authorHeinsoo, J. et al.
journaltitleRapid high-fidelity multiplexed readout of superconducting qubits.
Phys. Rev. Appl. volume10, pages034040, <10.1103/PhysRevApplied.10.034040> (year2018).

[16]
authorGlazman, L. & authorCatelani, G.
journaltitleBogoliubov quasiparticles in superconducting qubits.
SciPost Physics Lecture Notes <10.21468/scipostphyslectnotes.31> (year2021).

[17]
authorCatelani, G. & authorPekola, J. P.
journaltitleUsing materials for quasiparticle engineering.
Materials for Quantum Technology volume2, pages013001, <10.1088/2633-4356/ac4a75> (year2022).

[18]
authorCatelani, G. et al.
journaltitleQuasiparticle relaxation of superconducting qubits in the presence of flux.
Phys. Rev. Lett. volume106, pages077002, <10.1103/PhysRevLett.106.077002> (year2011).

[19]
authorCatelani, G., authorSchoelkopf, R. J., authorDevoret, M. H. & authorGlazman, L. I.
journaltitleRelaxation and frequency shifts induced by quasiparticles in superconducting qubits.
Phys. Rev. B volume84, pages064517, <10.1103/PhysRevB.84.064517> (year2011).

[20]
authorCórcoles, A. D. et al.
journaltitleProtecting superconducting qubits from radiation.
Applied Physics Letters volume99, pages181906, <10.1063/1.3658630> (year2011).

[21]
authorCatelani, G., authorNigg, S. E., authorGirvin, S. M., authorSchoelkopf, R. J. & authorGlazman, L. I.
journaltitleDecoherence of superconducting qubits caused by quasiparticle tunneling.
Phys. Rev. B volume86, pages184514, <10.1103/PhysRevB.86.184514> (year2012).

[22]
authorCatelani, G.
journaltitleParity switching and decoherence by quasiparticles in single-junction transmons.
Phys. Rev. B volume89, pages094522, <10.1103/PhysRevB.89.094522> (year2014).

[23]
authorPop, I. M. et al.
journaltitleCoherent suppression of electromagnetic dissipation due to superconducting quasiparticles.
Nature volume508, pages369–372, <10.1038/nature13017> (year2014).

[24]
authorRiwar, R.-P. et al.
journaltitleNormal-metal quasiparticle traps for superconducting qubits.
Phys. Rev. B volume94, pages104516, <10.1103/PhysRevB.94.104516> (year2016).

[25]
authorHosseinkhani, A., authorRiwar, R.-P., authorSchoelkopf, R. J., authorGlazman, L. I. & authorCatelani, G.
journaltitleOptimal configurations for normal-metal traps in transmon qubits.
Phys. Rev. Appl. volume8, pages064028, <10.1103/PhysRevApplied.8.064028> (year2017).

[26]
authorHosseinkhani, A. & authorCatelani, G.
journaltitleProximity effect in normal-metal quasiparticle traps.
Phys. Rev. B volume97, pages054513, <10.1103/PhysRevB.97.054513> (year2018).

[27]
authorSun, L. et al.
journaltitleMeasurements of quasiparticle tunneling dynamics in a band-gap-engineered transmon qubit.
Phys. Rev. Lett. volume108, pages230509, <10.1103/PhysRevLett.108.230509> (year2012).

[28]
authorRiwar, R.-P. & authorCatelani, G.
journaltitleEfficient quasiparticle traps with low dissipation through gap engineering.
Phys. Rev. B volume100, pages144514, <10.1103/PhysRevB.100.144514> (year2019).

[29]
authorMarchegiani, G., authorAmico, L. & authorCatelani, G.
journaltitleQuasiparticles in superconducting qubits with asymmetric junctions.
PRX Quantum volume3, pages040338, <10.1103/PRXQuantum.3.040338> (year2022).

[30]
authorPan, X. et al.
journaltitleEngineering superconducting qubits to reduce quasiparticles and charge noise.
Nature Communications volume13, pages7196, <10.1038/s41467-022-34727-2> (year2022).

[31]
authorSmith, W. C., authorKou, A., authorXiao, X., authorVool, U. & authorDevoret, M. H.
journaltitleSuperconducting circuit protected by two-cooper-pair tunneling.
npj Quantum Information volume6, pages8, <10.1038/s41534-019-0231-2> (year2020).

[32]
authorSung, Y. et al.
journaltitleRealization of high-fidelity cz and $zz$-free iswap gates with a tunable coupler.
Phys. Rev. X volume11, pages021058, <10.1103/PhysRevX.11.021058> (year2021).

[33]
authorXu, H. et al.
journaltitleRealization of adiabatic and diabatic CZ gates in superconducting qubits coupled with a tunable coupler.
Chinese Physics B volume30, pages044212, <10.1088/1674-1056/abf03a> (year2021).

[34]
authorYan, F. et al.
journaltitleTunable coupling scheme for implementing high-fidelity two-qubit gates.
Phys. Rev. Appl. volume10, pages054062, <10.1103/PhysRevApplied.10.054062> (year2018).

[35]
authorChu, J. & authorYan, F.
journaltitleCoupler-assisted controlled-phase gate with enhanced adiabaticity.
Phys. Rev. Appl. volume16, pages054020, <10.1103/PhysRevApplied.16.054020> (year2021).

[36]
authorPapič, M., authorAuer, A. & authorde Vega, I.
titleError sources of quantum gates in superconducting qubits (year2023).

[37]
authorMarxer, F. et al.
journaltitleLong-distance transmon coupler with cz-gate fidelity above $99.8\%$.
PRX Quantum volume4, pages010314, <10.1103/PRXQuantum.4.010314> (year2023).

[38]
authorSete, E. A., authorChen, A. Q., authorManenti, R., authorKulshreshtha, S. & authorPoletto, S.
journaltitleFloating tunable coupler for scalable quantum computing architectures.
Phys. Rev. Appl. volume15, pages064063, <10.1103/PhysRevApplied.15.064063> (year2021).

[39]
authorGoto, H.
journaltitleDouble-transmon coupler: Fast two-qubit gate with no residual coupling for highly detuned superconducting qubits.
Phys. Rev. Appl. volume18, pages034038, <10.1103/PhysRevApplied.18.034038> (year2022).

[40]
authorStehlik, J. et al.
journaltitleTunable coupling architecture for fixed-frequency transmon superconducting qubits.
Phys. Rev. Lett. volume127, pages080505, <10.1103/PhysRevLett.127.080505> (year2021).

[41]
authorXu, Y. et al.
journaltitleHigh-fidelity, high-scalability two-qubit gate scheme for superconducting qubits.
Phys. Rev. Lett. volume125, pages240503, <10.1103/PhysRevLett.125.240503> (year2020).

[42]
authorJurcevic, P. & authorGovia, L. C. G.
journaltitleEffective qubit dephasing induced by spectator-qubit relaxation.
Quantum Science and Technology volume7, pages045033, <10.1088/2058-9565/ac8cad> (year2022).

[43]
authorMundada, P., authorZhang, G., authorHazard, T. & authorHouck, A.
journaltitleSuppression of qubit crosstalk in a tunable coupling superconducting circuit.
Phys. Rev. Appl. volume12, pages054023, <10.1103/PhysRevApplied.12.054023> (year2019).

[44]
authorWang, C. et al.
journaltitleTowards practical quantum computers: transmon qubit with a lifetime approaching 0.5 milliseconds.
npj Quantum Information volume8, pages3, <10.1038/s41534-021-00510-2> (year2022).

[45]
authorPlace, A. P. M. et al.
journaltitleNew material platform for superconducting transmon qubits with coherence times exceeding 0.3 milliseconds.
Nature Communications volume12, pages1779, <10.1038/s41467-021-22030-5> (year2021).

[46]
authorWillsch, D. et al.
titleObservation of josephson harmonics in tunnel junctions (year2023).

[47]
authorSerniak, K. et al.
journaltitleDirect dispersive monitoring of charge parity in offset-charge-sensitive transmons.
Physical Review Applied volume12, <10.1103/physrevapplied.12.014052> (year2019).

[48]
authorRol, M. et al.
journaltitleFast, high-fidelity conditional-phase gate exploiting leakage interference in weakly anharmonic superconducting qubits.
Physical Review Letters volume123, <10.1103/physrevlett.123.120502> (year2019).

[49]
authorCollodo, M. C. et al.
journaltitleImplementation of conditional phase gates based on tunable $zz$ interactions.
Phys. Rev. Lett. volume125, pages240502, <10.1103/PhysRevLett.125.240502> (year2020).

[50]
authorHeunisch, L., authorEichler, C. & authorHartmann, M. J.
titleTunable coupler to fully decouple superconducting qubits (year2023).

[51]
authorMcKay, D. C., authorWood, C. J., authorSheldon, S., authorChow, J. M. & authorGambetta, J. M.
journaltitleEfficient $z$ gates for quantum computing.
Phys. Rev. A volume96, pages022330, <10.1103/PhysRevA.96.022330> (year2017).

[52]
authorNielsen, M. A. & authorChuang, I. L.
titleQuantum Computation and Quantum Information (publisherCambridge University Press, year2000).

[53]
authorDiamond, S. et al.
journaltitleDistinguishing parity-switching mechanisms in a superconducting qubit.
PRX Quantum volume3, pages040304, <10.1103/PRXQuantum.3.040304> (year2022).

[54]
authorRistè, D. et al.
journaltitleMillisecond charge-parity fluctuations and induced decoherence in a superconducting transmon qubit.
Nature Communications volume4, <10.1038/ncomms2936> (year2013).

[55]
authorSerniak, K. et al.
journaltitleHot nonequilibrium quasiparticles in transmon qubits.
Phys. Rev. Lett. volume121, pages157701, <10.1103/PhysRevLett.121.157701> (year2018).

[56]
authorKurter, C. et al.
journaltitleQuasiparticle tunneling as a probe of josephson junction barrier and capacitor material in superconducting qubits.
npj Quantum Information volume8, pages31, <10.1038/s41534-022-00542-2> (year2022).

[57]
authorTennant, D. M. et al.
journaltitleLow-frequency correlated charge-noise measurements across multiple energy transitions in a tantalum transmon.
PRX Quantum volume3, pages030307, <10.1103/PRXQuantum.3.030307> (year2022).

[58]
authorThorbeck, T., authorEddins, A., authorLauer, I., authorMcClure, D. T. & authorCarroll, M.
journaltitleTwo-level-system dynamics in a superconducting qubit due to background ionizing radiation.
PRX Quantum volume4, pages020356, <10.1103/PRXQuantum.4.020356> (year2023).

[59]
authorMcEwen, M. et al.
journaltitleResolving catastrophic error bursts from cosmic rays in large arrays of superconducting qubits.
Nature Physics volume18, pages107–111, <10.1038/s41567-021-01432-8> (year2021).

[60]
authorCardani, L. et al.
journaltitleReducing the impact of radioactivity on quantum circuits in a deep-underground facility.
Nature Communications volume12, pages2733, <10.1038/s41467-021-23032-z> (year2021).

[61]
authorGrünhaupt, L. et al.
journaltitleLoss mechanisms and quasiparticle dynamics in superconducting microwave resonators made of thin-film granular aluminum.
Phys. Rev. Lett. volume121, pages117001, <10.1103/PhysRevLett.121.117001> (year2018).

[62]
authorWack, A. et al.
titleQuality, speed, and scale: three key attributes to measure the performance of near-term quantum computers (year2021).

[63]
authorBreuer, H., authorPetruccione, F. & authorPetruccione, S.
titleThe Theory of Open Quantum Systems (publisherOxford University Press, year2002).

[64]
authorEmerson, J., authorAlicki, R. & authorŻyczkowski, K.
journaltitleScalable noise estimation with random unitary operators.
Journal of Optics B: Quantum and Semiclassical Optics volume7, pagesS347, <10.1088/1464-4266/7/10/021> (year2005).

[65]
authorWood, C. J. & authorGambetta, J. M.
journaltitleQuantification and characterization of leakage errors.
Phys. Rev. A volume97, pages032306, <10.1103/PhysRevA.97.032306> (year2018).

[66]
authorHashim, A. et al.
journaltitleRandomized compiling for scalable quantum computing on a noisy superconducting quantum processor.
Phys. Rev. X volume11, pages041039, <10.1103/PhysRevX.11.041039> (year2021).

[67]
authorWallman, J. J. & authorFlammia, S. T.
journaltitleRandomized benchmarking with confidence.
New Journal of Physics volume16, pages103032, <10.1088/1367-2630/16/10/103032> (year2014).

[68]
authorSiddiqi, I.
journaltitleEngineering high-coherence superconducting qubits.
Nature Reviews Materials volume6, pages875–891, <10.1038/s41578-021-00370-4> (year2021).

[69]
authorMüller, C., authorCole, J. H. & authorLisenfeld, J.
journaltitleTowards understanding two-level-systems in amorphous solids: insights from quantum circuits.
Reports on Progress in Physics volume82, pages124501, <10.1088/1361-6633/ab3a7e> (year2019).

[70]
authorCarroll, M., authorRosenblatt, S., authorJurcevic, P., authorLauer, I. & authorKandala, A.
journaltitleDynamics of superconducting qubit relaxation times.
npj Quantum Information volume8, pages132, <10.1038/s41534-022-00643-y> (year2022).

[71]
authorPremkumar, A. et al.
journaltitleMicroscopic relaxation channels in materials for superconducting qubits.
Communications Materials volume2, pages72, <10.1038/s43246-021-00174-7> (year2021).

[72]
authorCho, Y. et al.
titleSimulating noise on a quantum processor: interactions between a qubit and resonant two-level system bath, <10.48550/ARXIV.2211.08535> (year2022).

[73]
authorLisenfeld, J. et al.
journaltitleElectric field spectroscopy of material defects in transmon qubits.
npj Quantum Information volume5, <10.1038/s41534-019-0224-1> (year2019).

[74]
authorAbdurakhimov, L. V. et al.
journaltitleIdentification of different types of high-frequency defects in superconducting qubits.
PRX Quantum volume3, pages040332, <10.1103/PRXQuantum.3.040332> (year2022).

[75]
authorRower, D. A. et al.
journaltitleEvolution of $1/f$ flux noise in superconducting qubits with weak magnetic fields.
Phys. Rev. Lett. volume130, pages220602, <10.1103/PhysRevLett.130.220602> (year2023).

[76]
authorBraumüller, J. et al.
journaltitleCharacterizing and optimizing qubit coherence based on squid geometry.
Phys. Rev. Appl. volume13, pages054079, <10.1103/PhysRevApplied.13.054079> (year2020).

[77]
authorYan, F. et al.
journaltitleThe flux qubit revisited to enhance coherence and reproducibility.
Nature Communications volume7, <10.1038/ncomms12964> (year2016).

[78]
authorMotzoi, F., authorGambetta, J. M., authorRebentrost, P. & authorWilhelm, F. K.
journaltitleSimple pulses for elimination of leakage in weakly nonlinear qubits.
Phys. Rev. Lett. volume103, pages110501, <10.1103/PhysRevLett.103.110501> (year2009).

[79]
authorBerke, C., authorVarvelis, E., authorTrebst, S., authorAltland, A. & authorDiVincenzo, D. P.
journaltitleTransmon platform for quantum computing challenged by chaotic fluctuations.
Nature Communications volume13, pages2495, <10.1038/s41467-022-29940-y> (year2022).

[80]
authorAlgaba, M. G. et al.
journaltitleCo-design quantum simulation of nanoscale nmr.
Phys. Rev. Res. volume4, pages043089, <10.1103/PhysRevResearch.4.043089> (year2022).

[81]
authorChristensen, B. G. et al.
journaltitleAnomalous charge noise in superconducting qubits.
Phys. Rev. B volume100, pages140503, <10.1103/PhysRevB.100.140503> (year2019).

[82]
authorWilen, C. D. et al.
journaltitleCorrelated charge noise and relaxation errors in superconducting qubits.
Nature volume594, pages369–373, <10.1038/s41586-021-03557-5> (year2021).

[83]
authorBravyi, S., authorDiVincenzo, D. P. & authorLoss, D.
journaltitleSchrieffer–wolff transformation for quantum many-body systems.
Annals of Physics volume326, pages2793–2826, <10.1016/j.aop.2011.06.004> (year2011).

[84]
authorAbad, T., authorFernández-Pendás, J., authorFrisk Kockum, A. & authorJohansson, G.
journaltitleUniversal fidelity reduction of quantum operations from weak dissipation.
Phys. Rev. Lett. volume129, pages150504, <10.1103/PhysRevLett.129.150504> (year2022).

[85]
authorBurnett, J. J. et al.
journaltitleDecoherence benchmarking of superconducting qubits.
npj Quantum Information volume5, <10.1038/s41534-019-0168-5> (year2019).

[86]
authorPapi čč, M. & authorde Vega, I.
journaltitleNeural-network-based qubit-environment characterization.
Phys. Rev. A volume105, pages022605, <10.1103/PhysRevA.105.022605> (year2022).

[87]
authorBergli, J., authorGalperin, Y. M. & authorAltshuler, B. L.
journaltitleDecoherence in qubits due to low-frequency noise.
New Journal of Physics volume11, pages025002, <10.1088/1367-2630/11/2/025002> (year2009).

[88]
authorPaladino, E., authorGalperin, Y. M., authorFalci, G. & authorAltshuler, B. L.
journaltitle$1/f$ noise: Implications for solid-state quantum information.
Rev. Mod. Phys. volume86, pages361–418, <10.1103/RevModPhys.86.361> (year2014).

§ ACKNOWLEDGMENTS

We would like to thank Gianluigi Catelani and all of the employees at IQM for fruitful discussions, especially Johannes Heinsoo, Attila Geresdi, Antti Vepsäläinen, Alejandro Gómez Frieiro, Alessandro Landra, Fabian Marxer, Frank Deppe, and Vladimir Milchakov. We would additionally like to acknowledge the support from the German Federal Ministry of Education and Research (BMBF) under Q-Exa (grant No. 13N16062), QSolid (grant No. 13N16161) and MUNIQC-SC (grant No. 13N16188).

§ AUTHOR CONTRIBUTIONS

M.P. performed the analysis of the parity-switch effects and optimal parameter regime and wrote the draft of the manuscript. J.T. performed the analytical analysis of the tunable coupler system. A.A. and I. de V. supervised the work. A.H. conceptualized and supervised the project. All authors read, revised and approved the final manuscript.

§ COMPETING INTERESTS

The authors declare the following competing interests: The transmon parameter design optimization procedure presented in this manuscript is part of the patent application number: EP 23193779.8.

§ DATA AVAILABILITY

The datasets used and/or analyzed during the current study available from the corresponding authors on reasonable request.