[figure]style=plain,subcapbesideposition=top section[0em]  bx subsection[4em]
[.]1.5em

††thanks: Present address: School of Electrical Engineering and
Telecommunications, University of New South Wales, Sydney, New South Wales
2052, Australia.††thanks: Present address: Centre for Quantum Technologies,
National University of Singapore, 3 Science Drive 2, Singapore 117543,
Singapore.††thanks: Corresponding author<EMAIL_ADDRESS>
# Two-qubit sweet spots for capacitively coupled exchange-only spin qubits:
Supplementary Information

MengKe Feng Division of Physics and Applied Physics, School of Physical and
Mathematical Sciences, 21 Nanyang Link, Singapore 637371, Singapore. Lin Htoo
Zaw Division of Physics and Applied Physics, School of Physical and
Mathematical Sciences, 21 Nanyang Link, Singapore 637371, Singapore. Teck Seng
Koh Division of Physics and Applied Physics, School of Physical and
Mathematical Sciences, 21 Nanyang Link, Singapore 637371, Singapore.

## ontents

toc

## 0emupplementary Methods

###  1 Single Qubit Hamiltonian

We describe a single qubit system consisting of three electrons in a triple
quantum dot (TQD) with the Hubbard model [1, 2]:

$\displaystyle\mathcal{\hat{H}}=\hat{H}_{\varepsilon}+\hat{H}_{U}+\hat{H}_{t}$
(1)

where $\hat{H}_{\varepsilon}$ describes the detuning of the individual quantum
dots (QD), $\hat{H}_{t}$ describes tunneling between the dots, $\hat{H}_{U}$
describes the intra-dot and inter-dot Coulomb energies:

$\displaystyle\hat{H}_{\varepsilon}$
$\displaystyle=\sum_{i,\sigma}\varepsilon_{i}\hat{n}_{i\sigma}$ (2)
$\displaystyle\hat{H}_{t}$ $\displaystyle=\sum_{\langle
i,j\rangle,\sigma}t_{ij}\hat{c}^{\dagger}_{i,\sigma}\hat{c}_{j,\sigma}$ (3)
$\displaystyle\hat{H}_{U}$
$\displaystyle=\sum_{i}U_{i}\hat{n}_{i\uparrow}\hat{n}_{i\downarrow}+\frac{1}{2}\sum_{i}\sum_{j\neq
i}U_{ij}\hat{n}_{i}\hat{n}_{j}$ (4)

Here, $\hat{n}_{i,\sigma}=\hat{c}^{\dagger}_{i,\sigma}\hat{c}_{i,\sigma}$ is
the number operator and
$\hat{n}_{i}=\hat{n}_{i\uparrow}+\hat{n}_{i\downarrow}$, and $\varepsilon_{i}$
is the detuning of the $i$th-dot, $t_{ij}$ is the tunneling coupling between
neighboring dots, $U_{i}$ is the intra-dot Coulomb energy and $U_{ij}$ are
inter-dot Coulomb energy between all dots in the TQD. We assume identical QDs:
$U_{i}\equiv U$ for all $i$, $U_{12}=U_{23}=U_{45}=U_{56}\equiv U^{\prime}$
and $U_{13}=U_{46}\equiv U^{\prime\prime}$. We consider two identical qubits
(A and B) in a linear array. The basis states for qubit A are given by

$\displaystyle\ket{0_{A}}\equiv$
$\displaystyle\left(2\ket{\downarrow_{1}\uparrow_{2}\downarrow_{3}}-\ket{\downarrow_{1}\downarrow_{2}\uparrow_{3}}-\ket{\uparrow_{1}\downarrow_{2}\downarrow_{3}}\right)/\sqrt{6}$
(5) $\displaystyle\ket{1_{A}}\equiv$
$\displaystyle\left(\ket{\downarrow_{1}\downarrow_{2}\uparrow_{3}}-\ket{\uparrow_{1}\downarrow_{2}\downarrow_{3}}\right)/\sqrt{2}$
(6) $\displaystyle\ket{A_{1}}\equiv$
$\displaystyle\ket{\downarrow_{1}\uparrow_{3}\downarrow_{3}}$ (7)
$\displaystyle\ket{A_{2}}\equiv$
$\displaystyle\ket{\uparrow_{1}\downarrow_{1}\downarrow_{3}}$ (8)
$\displaystyle\ket{A_{3}}\equiv$
$\displaystyle\ket{\downarrow_{1}\uparrow_{2}\downarrow_{2}}$ (9)
$\displaystyle\ket{A_{4}}\equiv$
$\displaystyle\ket{\uparrow_{2}\downarrow_{2}\downarrow_{3}}.$ (10)

Those for qubit B follow accordingly. The doubly occupied states
$|A_{i}\rangle$ correspond to $(1,0,2),(2,0,1),(1,2,0)$ and $(0,2,1)$ charge
states for $i=1~{}\text{to}~{}4$. Note that these basis states are in the
$S_{z}=-1/2$ subspace, and the singlet and triplet states are symmetric about
the middle dot. Qubit A/B is the control/target qubit in CNOT and CPHASE
gates. The $(0,1,2)$ and $(2,1,0)$ charge configurations can be neglected
because they involve two tunneling processes and lead to higher order
admixtures. Dot detunings are defined by

$\displaystyle\varepsilon$
$\displaystyle=\frac{1}{2}(\varepsilon_{1}-\varepsilon_{3})$ (11)
$\displaystyle\varepsilon_{m}$
$\displaystyle=\varepsilon_{2}-\frac{1}{2}(\varepsilon_{1}+\varepsilon_{3})$
(12)

By shifting the energy of the system by the sum of the individual detunings,
the Hamiltonian can be expressed in terms of $\varepsilon$ and
$\varepsilon_{m}$. Performing a Schrieffer-Wolff transformation, we obtain the
effective Hamiltonian in the computational basis,

$\displaystyle\hat{H}_{\text{eff}}=-\frac{\hbar\omega}{2}\hat{\sigma}_{z}-\frac{\hbar
g}{2}\hat{\sigma}_{x}$ (13)

where $\hbar\omega=2\left(t_{l}^{2}b_{-}+t_{r}^{2}b_{+}\right)$ and $\hbar
g=2\sqrt{3}\left(t_{l}^{2}b_{-}-t_{r}^{2}b_{+}\right)$, with
$b_{\pm}=\frac{1}{\left(U-U^{\prime\prime}\pm\epsilon+\epsilon_{m}\right)}+\frac{1}{\left(U-2U^{\prime}+U^{\prime\prime}\mp\epsilon-\epsilon_{m}\right)}$.
The basis states for each qubit are now

$\displaystyle|{1}^{\prime}_{A}\rangle=$
$\displaystyle\frac{1}{\sqrt{N_{A}}}\left(|1_{A}\rangle+\sum_{i=1}^{4}\alpha_{i}|A_{i}\rangle\right),~{}~{}\text{and}~{}~{}|0_{A}\rangle,$
(14) $\displaystyle|{1}^{\prime}_{B}\rangle=$
$\displaystyle\frac{1}{\sqrt{N_{B}}}\left(|1_{B}\rangle+\sum_{i=1}^{4}\beta_{i}|B_{i}\rangle\right),~{}~{}\text{and}~{}~{}|0_{B}\rangle.$
(15)

where $\sqrt{N_{A/B}}$ are the normalization constants
($N_{A}=1+\alpha_{1}^{2}+\alpha_{2}^{2}+\alpha_{3}^{2}+\alpha_{4}^{2},N_{B}=1+\beta_{1}^{2}+\beta_{2}^{2}+\beta_{3}^{2}+\beta_{4}^{2}$)
and each of the $\alpha_{i},\beta_{i}$ terms are coefficients describing the
admixtures of the doubly occupied singlet states, given by

$\displaystyle\alpha_{1}$
$\displaystyle=\frac{t_{rA}}{\sqrt{2}(-U+U^{\prime\prime}+2\varepsilon_{A})},$
(16) $\displaystyle\alpha_{2}$
$\displaystyle=\frac{t_{lA}}{\sqrt{2}(-U+U^{\prime\prime}-2\varepsilon_{A})},$
(17) $\displaystyle\alpha_{3}$
$\displaystyle=\frac{t_{rA}}{\sqrt{2}(-U+2U^{\prime}-U^{\prime\prime}+\varepsilon_{A}+\varepsilon_{mA})},$
(18) $\displaystyle\alpha_{4}$
$\displaystyle=\frac{t_{lA}}{\sqrt{2}(-U+2U^{\prime}-U^{\prime\prime}-\varepsilon_{A}+\varepsilon_{mA})},$
(19) $\displaystyle\beta_{1}$
$\displaystyle=\frac{t_{rB}}{\sqrt{2}(-U+U^{\prime\prime}+2\varepsilon_{B})},$
(20) $\displaystyle\beta_{2}$
$\displaystyle=\frac{t_{lB}}{\sqrt{2}(-U+U^{\prime\prime}-2\varepsilon_{B})},$
(21) $\displaystyle\beta_{3}$
$\displaystyle=\frac{t_{rB}}{\sqrt{2}(-U+2U^{\prime}-U^{\prime\prime}+\varepsilon_{B}+\varepsilon_{mB})},$
(22) $\displaystyle\beta_{4}$
$\displaystyle=\frac{t_{lB}}{\sqrt{2}(-U+2U^{\prime}-U^{\prime\prime}-\varepsilon_{B}+\varepsilon_{mB}).}$
(23)

With the effective Hamiltonian (Eq. (13)) we are able to find the single-qubit
sweet spots, defined as flat points in the energy landscape satisfying [3]

$\displaystyle\frac{\partial\omega^{\prime}_{A/B}}{\partial\varepsilon_{A/B}}=$
$\displaystyle 0,$ (24)
$\displaystyle\frac{\partial\omega^{\prime}_{A/B}}{\partial\varepsilon_{mA/B}}=$
$\displaystyle 0,$ (25)

where $\omega^{\prime}_{A/B}=\sqrt{\omega_{A/B}^{2}+g_{A/B}^{2}}$ is the
eigenenergy difference between the ground and excited states for qubits $A$
and $B$. The double sweet spot where the AEON qubit is operated satisfies both
conditions. This exists at $\varepsilon_{A/B}=0$ meV and
$\varepsilon_{mA/B}=U^{\prime\prime}-U^{\prime}=-0.9~{}\mathrm{meV}$. On the
other hand, the RX qubit is generally defined to be operating in a region of
increased $\varepsilon_{mA/B}$ [3]. In our study, the RX qubit operates at the
same $\varepsilon_{A/B}=0~{}\mathrm{meV}$ single-qubit sweet spot, but at a
different middle detuning value, $\varepsilon_{mA/B}=-0.57~{}\mathrm{meV}$.
These choices of $\varepsilon_{mA/B}$ are so that we are able to obtain
comparable gate times for both the AEON and RX qubits, but generally, the RX
qubit can be operated at any value of $\varepsilon_{m}$ such that there is
strong admixture between the singly and doubly occupied states.

###  2 Triple Quantum Dot Potential

The confinement potential for a triple quantum dot system is modeled as a tri-
quadratic potential,

$\displaystyle
V_{\mathrm{pot}}(\vec{r})=\min{[v_{\mathrm{pot},1}(\vec{r}),v_{\mathrm{pot},2}(\vec{r}),v_{\mathrm{pot},3}(\vec{r})]},$
(26)

where $v_{\mathrm{pot},i}(\vec{r})$ refers to the potential well of each dot
centered at $\vec{R}_{i}$, given by

$\displaystyle
v_{\mathrm{pot},i}(\vec{r})=\frac{m\omega_{0}^{2}}{2}\left(|\vec{x}-\vec{R}_{i}|^{2}+\vec{y}^{2}\right)+\varepsilon_{i}.$
(27)

Here, $\omega_{0}$ is the confinement energy of the dot which depends on its
dimensions, and $\varepsilon_{i}$ is the detuning applied on the $i$-th dot.
The 2-dimensional character of this potential is a good approximation given
the tight confinement in the $z$-direction that is typical of
electrostatically gated quantum dots. The ground orbital of an electron in
such a potential is given by the Fock-Darwin wavefunction,

$\displaystyle\phi_{i}(\vec{r})=\left(\frac{1}{a_{B}\sqrt{\pi}}\right)\exp(-\frac{1}{2a_{B}^{2}}\left(|\vec{x}-\vec{R}_{i}|^{2}+\vec{y}^{2}\right)).$
(28)

This expression is normalized and $a_{B}=\sqrt{\hbar/m\omega_{0}}$ is the
radius of the quantum dot. We can now formulate dot-centered, orthonomalized
single-electron wavefunctions $\psi_{i}$ of the triple quantum dot potential
using the method of Löwdin orthogonalization [4, 5]. These are given by

$\displaystyle\Psi=\Phi M^{-1/2}$ (29)

where $\Psi=\\{\psi_{1},\psi_{2},\psi_{3}\\}$ is the set of dot-centered,
orthonomalized, single-electron wavefunctions,
$\Phi=\\{\phi_{1},\phi_{2},\phi_{3}\\}$ is the set of Fock-Darwin
wavefunctions, and $M$ is the overlap matrix given by

$\displaystyle M=\begin{pmatrix}1&s&s^{4}\\\ s&1&s\\\ s^{4}&s&1\end{pmatrix},$
(30)

where $s=\exp(-a^{2}/a_{B}^{2})$ and $2a$ is the distance between two
neighboring dots.

The single-electron wavefunctions are crucial in the calculation of the terms
in the interaction Hamiltonian, which is given by all pairwise direct Coulomb
interactions between electrons in each qubit. For three electrons in a triple
quantum dot, there are eighty-one pairwise direct Coulomb integrals. With the
quadratic potentials, each integral is analytically tractable and given
generally by

$\displaystyle\int\int\phi_{i}^{\star}(\vec{r}_{1})\phi_{j}^{\star}(\vec{r}_{2})\frac{\kappa}{|\vec{r}_{1}-\vec{r}_{2}|}\phi_{k}(\vec{r}_{1})\phi_{l}(\vec{r}_{2})~{}d\vec{r}_{1}d\vec{r}_{2}$
$\displaystyle=\kappa\sqrt{\frac{\pi}{2}}\frac{1}{a_{B}}\exp{-\frac{\left(R_{i}-R_{k}\right)^{2}-\left(R_{j}-R_{l}\right)^{2}}{4a_{B}^{2}}}\exp{-\frac{\left(R_{i}+R_{k}-R_{j}-R_{l}\right)^{2}}{16a_{B}^{2}}}I_{0}\left(\frac{\left(R_{i}+R_{k}-R_{j}-R_{l}\right)^{2}}{16a_{B}^{2}}\right)$
(31)

where $\kappa=q^{2}/4\pi\epsilon_{0}\epsilon_{r}$ and $I_{0}$ is the zeroth
order modified Bessel function of the first kind. The Bessel function appears
as the solution to a standard integral of the form
$\int_{0}^{1}\exp{-at^{2}}/\sqrt{1-t^{2}}~{}dt$. This result is similar to
that derived in Ref. [6]. In our calculations, we take relative permittivity
$\epsilon_{r}=11.68$ to be that for silicon.

###  3 Interaction Hamiltonian

The two-qubit capacitive interaction arises from the inter-dot Coulomb
interaction between the TQDs, given by

$\displaystyle\hat{H}_{\text{int}}=$
$\displaystyle\sum_{i=1}^{3}\sum_{j=4}^{6}\mathcal{V}_{ij}\hat{n}_{i}\hat{n}_{j},$
(32)

where $\mathcal{V}_{ij}$ is given by the direct Coulomb integral between the
orthonormalized electron densities of $i$-th and $j$-th QD,

$\displaystyle\mathcal{V}_{ij}=\int\int|\psi_{i}(\vec{r}_{1})|^{2}\frac{\kappa}{|\vec{r}_{1}-\vec{r}_{2}|}|\psi_{j}(\vec{r}_{2})|^{2}d\vec{r}_{1}d\vec{r}_{2}.$
(33)

Here, $i$-th dot-centered, normalized, single-electron wavefunctions
$\psi_{i}(\vec{r})$ are constructed using Löwdin’s orthogonalization method
from single QD wavefunctions as shown in Eq. (29). The three-electron
wavefunctions for qubits A and B follow from Eqs. (5)–(10), and the
corresponding hybridized basis states are as outlined in Eqs. (14)–(15).
Hence, the computational basis for the interaction Hamiltonian (Eq. (32))
comprises

$\displaystyle|0_{A}0_{B}\rangle$ , (34)
$\displaystyle|0_{A}1_{B}^{\prime}\rangle$
$\displaystyle=\frac{1}{\sqrt{N_{B}}}|0_{A}\rangle\left(|1_{B}\rangle+\sum_{n=1}^{4}\beta_{n}|B_{n}\rangle\right),$
(35) $\displaystyle|1_{A}^{\prime}0_{B}\rangle$
$\displaystyle=\frac{1}{\sqrt{N_{A}}}\left(|1_{A}\rangle+\sum_{m=1}^{4}\alpha_{m}|A_{m}\rangle\right)|0_{B}\rangle,$
(36) $\displaystyle|1_{A}^{\prime}1_{B}^{\prime}\rangle$
$\displaystyle=\frac{1}{\sqrt{N_{A}N_{B}}}\left(|1_{A}\rangle+\sum_{m=1}^{4}\alpha_{m}|A_{m}\rangle\right)\left(|1_{B}\rangle+\sum_{n=1}^{4}\beta_{n}|B_{n}\rangle\right).$
(37)

Because the basis states are eigenstates of the number operators, Eq. (32)
must be a diagonal matrix in the computational basis
$\\{\ket{0_{A}0_{B}},\ket{0_{A}1_{B}^{\prime}},\ket{1_{A}^{\prime}0_{B}},\ket{1_{A}^{\prime}1_{B}^{\prime}}\\}$,
i.e.

$\displaystyle\hat{H}_{\text{int}}=\hbar\times\text{diag}\\{V_{1},V_{2},V_{3},V_{4}\\}.$
(38)

The diagonal elements are explicitly calculated as follows.

$\displaystyle\hbar V_{1}\equiv$ $\displaystyle\langle
0_{A}0_{B}|\left(\sum_{i=1}^{3}\sum_{j=4}^{6}\mathcal{V}_{ij}\hat{n}_{i}\hat{n}_{j}\right)|0_{A}0_{B}\rangle$
$\displaystyle=$
$\displaystyle\mathcal{V}_{14}+\mathcal{V}_{24}+\mathcal{V}_{34}+\mathcal{V}_{15}+\mathcal{V}_{25}+\mathcal{V}_{35}+\mathcal{V}_{16}+\mathcal{V}_{26}+\mathcal{V}_{36}\equiv
C_{11},$ (39)

where we have defined $C_{11}$ as the sum of all the possible inter-qubit
capacitive interaction terms. The $V_{2}$ term is given as follows.

$\displaystyle\hbar V_{2}\equiv$ $\displaystyle\langle
0_{A}1_{B}^{\prime}|\left(\sum_{i=1}^{3}\sum_{j=4}^{6}\mathcal{V}_{ij}\hat{n}_{i}\hat{n}_{j}\right)|0_{A}1_{B}^{\prime}\rangle$
$\displaystyle=$ $\displaystyle\frac{1}{N_{B}}\left(\langle 0_{A}|\langle
1_{B}|+\sum_{m=1}^{4}\beta_{m}^{\star}\langle 0_{A}|\langle
B_{m}|\right)\left(\sum_{i=1}^{3}\sum_{j=4}^{6}\mathcal{V}_{ij}\hat{n}_{i}\hat{n}_{j}\right)\left(|0_{A}\rangle|1_{B}\rangle+\sum_{n=1}^{4}\beta_{n}|0_{A}\rangle|B_{n}\rangle\right)$
$\displaystyle=$
$\displaystyle\frac{1}{N_{B}}\left(\mathcal{V}_{14}+\mathcal{V}_{24}+\mathcal{V}_{34}+\mathcal{V}_{15}+\mathcal{V}_{25}+\mathcal{V}_{35}+\mathcal{V}_{16}+\mathcal{V}_{26}+\mathcal{V}_{36}\right)$
$\displaystyle+\frac{1}{N_{B}}\left(\beta_{1}^{2}\left(\mathcal{V}_{14}+\mathcal{V}_{24}+\mathcal{V}_{34}+2\mathcal{V}_{16}+2\mathcal{V}_{26}+2\mathcal{V}_{36}\right)\right)$
$\displaystyle+\frac{1}{N_{B}}\left(\beta_{2}^{2}\left(2\mathcal{V}_{14}+2\mathcal{V}_{24}+2\mathcal{V}_{34}+\mathcal{V}_{16}+\mathcal{V}_{26}+\mathcal{V}_{36}\right)\right)$
$\displaystyle+\frac{1}{N_{B}}\left(\beta_{3}^{2}\left(\mathcal{V}_{14}+\mathcal{V}_{24}+\mathcal{V}_{34}+2\mathcal{V}_{15}+2\mathcal{V}_{25}+2\mathcal{V}_{35}\right)\right)$
$\displaystyle+\frac{1}{N_{B}}\left(\beta_{4}^{2}\left(2\mathcal{V}_{15}+2\mathcal{V}_{25}+2\mathcal{V}_{35}+\mathcal{V}_{16}+\mathcal{V}_{26}+\mathcal{V}_{36}\right)\right)$
$\displaystyle\equiv$
$\displaystyle\frac{1}{N_{B}}\left(C_{11}+\beta_{1}^{2}C_{12}+\beta_{2}^{2}C_{13}+\beta_{3}^{2}C_{14}+\beta_{4}^{2}C_{15}\right)$
$\displaystyle\equiv$
$\displaystyle\frac{1}{N_{B}}\left(C_{11}+\vec{\nu}\cdot\vec{F}_{B}\right),$
(40)

where the exact terms can be derived from the basis states as defined in Eq.
(5)–(10) and we have defined

$\displaystyle C_{12}\equiv$
$\displaystyle\mathcal{V}_{14}+\mathcal{V}_{24}+\mathcal{V}_{34}+2\mathcal{V}_{16}+2\mathcal{V}_{26}+2\mathcal{V}_{36},$
(41) $\displaystyle C_{13}\equiv$ $\displaystyle
2\mathcal{V}_{14}+2\mathcal{V}_{24}+2\mathcal{V}_{34}+\mathcal{V}_{16}+\mathcal{V}_{26}+\mathcal{V}_{36},$
(42) $\displaystyle C_{14}\equiv$
$\displaystyle\mathcal{V}_{14}+\mathcal{V}_{24}+\mathcal{V}_{34}+2\mathcal{V}_{15}+2\mathcal{V}_{25}+2\mathcal{V}_{35},$
(43) $\displaystyle C_{15}\equiv$ $\displaystyle
2\mathcal{V}_{15}+2\mathcal{V}_{25}+2\mathcal{V}_{35}+\mathcal{V}_{16}+\mathcal{V}_{26}+\mathcal{V}_{36},$
(44) $\displaystyle\vec{F}_{B}\equiv$
$\displaystyle\\{C_{12},C_{13},C_{14},C_{15}\\}^{\text{T}},$ (45)
$\displaystyle\vec{\nu}\equiv$
$\displaystyle\\{\beta_{1}^{2},\beta_{2}^{2},\beta_{3}^{2},\beta_{4}^{2}\\}^{\text{T}}.$
(46)

Next,

$\displaystyle\hbar V_{3}\equiv$ $\displaystyle\langle
1_{A}^{\prime}0_{B}|\left(\sum_{i=1}^{3}\sum_{j=4}^{6}\mathcal{V}_{ij}\hat{n}_{i}\hat{n}_{j}\right)|1_{A}^{\prime}0_{B}\rangle$
$\displaystyle=$ $\displaystyle\frac{1}{N_{A}}\left(\langle 1_{A}|\langle
0_{B}|+\sum_{m=1}^{4}\alpha_{m}^{\star}\langle A_{m}|\langle
0_{B}|\right)\left(\sum_{i=1}^{3}\sum_{j=4}^{6}\mathcal{V}_{ij}\hat{n}_{i}\hat{n}_{j}\right)\left(|1_{A}\rangle|0_{B}\rangle+\sum_{n=1}^{4}\alpha_{n}|A_{m}\rangle|0_{B}\rangle\right)$
$\displaystyle=$
$\displaystyle\frac{1}{N_{A}}\left(\mathcal{V}_{14}+\mathcal{V}_{24}+\mathcal{V}_{34}+\mathcal{V}_{15}+\mathcal{V}_{25}+\mathcal{V}_{35}+\mathcal{V}_{16}+\mathcal{V}_{26}+\mathcal{V}_{36}\right)$
$\displaystyle+\frac{1}{N_{A}}\left(\alpha_{1}^{2}\left(\mathcal{V}_{14}+\mathcal{V}_{15}+\mathcal{V}_{16}+2\mathcal{V}_{34}+2\mathcal{V}_{35}+2\mathcal{V}_{36}\right)\right)$
$\displaystyle+\frac{1}{N_{A}}\left(\alpha_{2}^{2}\left(2\mathcal{V}_{14}+2\mathcal{V}_{15}+2\mathcal{V}_{16}+\mathcal{V}_{34}+\mathcal{V}_{35}+\mathcal{V}_{36}\right)\right)$
$\displaystyle+\frac{1}{N_{A}}\left(\alpha_{3}^{2}\left(\mathcal{V}_{14}+\mathcal{V}_{15}+\mathcal{V}_{16}+2\mathcal{V}_{24}+2\mathcal{V}_{25}+2\mathcal{V}_{26}\right)\right)$
$\displaystyle+\frac{1}{N_{A}}\left(\alpha_{4}^{2}\left(2\mathcal{V}_{24}+2\mathcal{V}_{25}+2\mathcal{V}_{26}+\mathcal{V}_{34}+\mathcal{V}_{35}+\mathcal{V}_{36}\right)\right)$
$\displaystyle\equiv$
$\displaystyle\frac{1}{N_{A}}\left(C_{11}+\alpha_{1}^{2}C_{21}+\alpha_{2}^{2}C_{31}+\alpha_{3}^{2}C_{41}+\alpha_{4}^{2}C_{51}\right)$
$\displaystyle\equiv$
$\displaystyle\frac{1}{N_{A}}\left(C_{11}+\vec{\xi}\cdot\vec{F}_{A}\right).$
(47)

where we have defined

$\displaystyle C_{21}\equiv$
$\displaystyle\mathcal{V}_{14}+\mathcal{V}_{15}+\mathcal{V}_{16}+2\mathcal{V}_{34}+2\mathcal{V}_{35}+2\mathcal{V}_{36},$
(48) $\displaystyle C_{31}\equiv$ $\displaystyle
2\mathcal{V}_{14}+2\mathcal{V}_{15}+2\mathcal{V}_{16}+\mathcal{V}_{34}+\mathcal{V}_{35}+\mathcal{V}_{36},$
(49) $\displaystyle C_{41}\equiv$
$\displaystyle\mathcal{V}_{14}+\mathcal{V}_{15}+\mathcal{V}_{16}+2\mathcal{V}_{24}+2\mathcal{V}_{25}+2\mathcal{V}_{26},$
(50) $\displaystyle C_{51}\equiv$ $\displaystyle
2\mathcal{V}_{24}+2\mathcal{V}_{25}+2\mathcal{V}_{26}+\mathcal{V}_{34}+\mathcal{V}_{35}+\mathcal{V}_{36},$
(51) $\displaystyle\vec{F}_{A}\equiv$
$\displaystyle\\{C_{21},C_{31},C_{41},C_{51}\\}^{\text{T}},$ (52)
$\displaystyle\vec{\xi}\equiv$
$\displaystyle\\{\alpha_{1}^{2},\alpha_{2}^{2},\alpha_{3}^{2},\alpha_{4}^{2}\\}^{\text{T}}$
(53)

The $V_{4}$ term is

$\displaystyle\hbar V_{4}\equiv$ $\displaystyle\langle
1^{\prime}_{A}1^{\prime}_{B}|\left(\sum_{i=1}^{3}\sum_{j=4}^{6}\mathcal{V}_{ij}\hat{n}_{i}\hat{n}_{j}\right)|1^{\prime}_{A}1^{\prime}_{B}\rangle$
$\displaystyle=$ $\displaystyle\frac{1}{N_{A}N_{B}}\left[\left(\langle
1_{A}|+\sum_{m=1}^{4}\alpha_{m}^{\star}\langle A_{m}|\right)\left(\langle
1_{B}|+\sum_{n=1}^{4}\beta_{n}^{\star}\langle
B_{n}|\right)\right]\left(\sum_{i=1}^{3}\sum_{j=4}^{6}\mathcal{V}_{ij}\hat{n}_{i}\hat{n}_{j}\right)\times$
$\displaystyle\left[\left(|1_{A}\rangle+\sum_{p=1}^{4}\alpha_{p}|A_{p}\rangle\right)\left(|1_{B}\rangle+\sum_{q=1}^{4}\beta_{q}|B_{q}\rangle\right)\right].$
(54)

Expanding the kets, we group terms according to the pre-factors,

$\displaystyle\left(|1_{A}\rangle+\sum_{p=1}^{4}\alpha_{p}|A_{p}\rangle\right)\left(|1_{B}\rangle+\sum_{q=1}^{4}\beta_{q}|B_{q}\rangle\right)$
$\displaystyle=|1_{A}\rangle|1_{B}\rangle+\sum_{p=1}^{4}\alpha_{p}|A_{p}\rangle|1_{B}\rangle+\sum_{q=1}^{4}\beta_{q}|1_{A}\rangle|B_{q}\rangle+\sum_{p=1}^{4}\sum_{q=1}^{4}\alpha_{p}\beta_{q}|A_{p}\rangle|B_{q}\rangle.$
(55)

There are no cross terms since all states are orthogonal, so the terms in
$V_{4}$ can be grouped into those without admixture factors, and terms with
admixture factors $\alpha_{p}^{2}$, $\beta_{q}^{2}$ and
$\alpha_{p}^{2}\beta_{q}^{2}$.

The first group of terms without admixtures are

$\displaystyle\langle 1_{A}|\langle
1_{B}|\left(\sum_{i=1}^{3}\sum_{j=4}^{6}\mathcal{V}_{ij}\hat{n}_{i}\hat{n}_{j}\right)|1_{A}\rangle|1_{B}\rangle=$
$\displaystyle\mathcal{V}_{14}+\mathcal{V}_{24}+\mathcal{V}_{34}+\mathcal{V}_{15}+\mathcal{V}_{25}+\mathcal{V}_{35}+\mathcal{V}_{16}+\mathcal{V}_{26}+\mathcal{V}_{36}\equiv
C_{11}.$ (56)

The second group of $\alpha_{p}^{2}$ terms are

$\displaystyle\left(\sum_{m=1}^{4}\alpha_{m}\langle A_{p}|\langle
1_{B}|\right)\left(\sum_{i=1}^{3}\sum_{j=4}^{6}\mathcal{V}_{ij}\hat{n}_{i}\hat{n}_{j}\right)\left(\sum_{p=1}^{4}\alpha_{p}|A_{p}\rangle|1_{B}\rangle\right)$
$\displaystyle=$
$\displaystyle\alpha_{1}^{2}\left(\mathcal{V}_{14}+\mathcal{V}_{15}+\mathcal{V}_{16}+2\mathcal{V}_{34}+2\mathcal{V}_{35}+2\mathcal{V}_{36}\right)$
$\displaystyle+\alpha_{2}^{2}\left(2\mathcal{V}_{14}+2\mathcal{V}_{15}+2\mathcal{V}_{16}+\mathcal{V}_{34}+\mathcal{V}_{35}+\mathcal{V}_{36}\right)$
$\displaystyle+\alpha_{3}^{2}\left(\mathcal{V}_{14}+\mathcal{V}_{15}+\mathcal{V}_{16}+2\mathcal{V}_{24}+2\mathcal{V}_{25}+2\mathcal{V}_{26}\right)$
$\displaystyle+\alpha_{4}^{2}\left(2\mathcal{V}_{24}+2\mathcal{V}_{25}+2\mathcal{V}_{26}+\mathcal{V}_{34}+\mathcal{V}_{35}+\mathcal{V}_{36}\right)$
$\displaystyle=$
$\displaystyle\alpha_{1}^{2}C_{21}+\alpha_{2}^{2}C_{31}+\alpha_{3}^{2}C_{41}+\alpha_{4}^{2}C_{51}=\vec{\xi}\cdot\vec{F}_{A}.$
(57)

The third group of $\beta_{q}^{2}$ terms are

$\displaystyle\left(\sum_{n=1}^{4}\beta_{n}\langle 1_{A}|\langle
B_{n}|\right)\left(\sum_{i=1}^{3}\sum_{j=4}^{6}\mathcal{V}_{ij}\hat{n}_{i}\hat{n}_{j}\right)\left(\sum_{q=1}^{4}\beta_{q}|1_{A}\rangle|B_{q}\rangle\right)$
$\displaystyle=$
$\displaystyle\beta_{1}^{2}\left(\mathcal{V}_{14}+\mathcal{V}_{24}+\mathcal{V}_{34}+2\mathcal{V}_{16}+2\mathcal{V}_{26}+2\mathcal{V}_{36}\right)$
$\displaystyle+\beta_{2}^{2}\left(2\mathcal{V}_{14}+2\mathcal{V}_{24}+2\mathcal{V}_{34}+\mathcal{V}_{16}+\mathcal{V}_{26}+\mathcal{V}_{36}\right)$
$\displaystyle+\beta_{3}^{2}\left(\mathcal{V}_{14}+\mathcal{V}_{24}+\mathcal{V}_{34}+2\mathcal{V}_{15}+2\mathcal{V}_{25}+2\mathcal{V}_{35}\right)$
$\displaystyle+\beta_{4}^{2}\left(2\mathcal{V}_{15}+2\mathcal{V}_{25}+2\mathcal{V}_{35}+\mathcal{V}_{16}+\mathcal{V}_{26}+\mathcal{V}_{36}\right)$
$\displaystyle=$
$\displaystyle\beta_{1}^{2}C_{12}+\beta_{2}^{2}C_{13}+\beta_{3}^{2}C_{14}+\beta_{4}^{2}C_{15}=\vec{\nu}\cdot\vec{F}_{B}.$
(58)

The fourth group of $\alpha_{p}^{2}\beta_{q}^{2}$ terms are

$\displaystyle\left(\sum_{m=1}^{4}\sum_{n=1}^{4}\alpha_{m}^{\star}\beta_{n}^{\star}\langle
A_{m}|\langle
B_{n}|\right)\left(\sum_{i=1}^{3}\sum_{j=4}^{6}\mathcal{V}_{ij}\hat{n}_{i}\hat{n}_{j}\right)\left(\sum_{p=1}^{4}\sum_{q=1}^{4}\alpha_{p}\beta_{q}|A_{p}\rangle|B_{q}\rangle\right)$
$\displaystyle=$
$\displaystyle\alpha_{1}^{2}\left(\beta_{1}^{2}C_{22}+\beta_{2}^{2}C_{23}+\beta_{3}^{2}C_{24}+\beta_{4}^{2}C_{25}\right)$
$\displaystyle+\alpha_{2}^{2}\left(\beta_{1}^{2}C_{32}+\beta_{2}^{2}C_{33}+\beta_{3}^{2}C_{34}+\beta_{4}^{2}C_{35}\right)$
$\displaystyle+\alpha_{3}^{2}\left(\beta_{1}^{2}C_{42}+\beta_{2}^{2}C_{43}+\beta_{3}^{2}C_{44}+\beta_{4}^{2}C_{45}\right)$
$\displaystyle+\alpha_{4}^{2}\left(\beta_{1}^{2}C_{52}+\beta_{2}^{2}C_{53}+\beta_{3}^{2}C_{54}+\beta_{4}^{2}C_{55}\right)$
$\displaystyle=$
$\displaystyle\vec{\xi}~{}^{\text{T}}\overleftrightarrow{F}\vec{\nu}.$ (59)

In the last line, we have defined the matrix $\overleftrightarrow{F}\equiv
C_{ij}$ for $i,j=2,3,4,5$, and column vectors $\vec{\xi}$, $\vec{\nu}$. The
definitions of $C_{ij}$ above are

$\displaystyle C_{22}\equiv$
$\displaystyle\mathcal{V}_{14}+2\mathcal{V}_{16}+2\mathcal{V}_{34}+4\mathcal{V}_{36},$
(60) $\displaystyle C_{23}\equiv$ $\displaystyle
2\mathcal{V}_{14}+\mathcal{V}_{16}+4\mathcal{V}_{34}+2\mathcal{V}_{36},$ (61)
$\displaystyle C_{24}\equiv$
$\displaystyle\mathcal{V}_{14}+2\mathcal{V}_{15}+2\mathcal{V}_{34}+4\mathcal{V}_{35},$
(62) $\displaystyle C_{25}\equiv$ $\displaystyle
2\mathcal{V}_{15}+\mathcal{V}_{16}+4\mathcal{V}_{35}+2\mathcal{V}_{36},$ (63)
$\displaystyle C_{32}\equiv$ $\displaystyle
2\mathcal{V}_{14}+4\mathcal{V}_{16}+\mathcal{V}_{34}+2\mathcal{V}_{36},$ (64)
$\displaystyle C_{33}\equiv$ $\displaystyle
4\mathcal{V}_{14}+2\mathcal{V}_{16}+2\mathcal{V}_{34}+\mathcal{V}_{36},$ (65)
$\displaystyle C_{34}\equiv$ $\displaystyle
2\mathcal{V}_{14}+4\mathcal{V}_{15}+\mathcal{V}_{34}+2\mathcal{V}_{35},$ (66)
$\displaystyle C_{35}\equiv$ $\displaystyle
4\mathcal{V}_{15}+2\mathcal{V}_{16}+2\mathcal{V}_{35}+\mathcal{V}_{36},$ (67)
$\displaystyle C_{42}\equiv$
$\displaystyle\mathcal{V}_{14}+2\mathcal{V}_{16}+2\mathcal{V}_{24}+4\mathcal{V}_{26},$
(68) $\displaystyle C_{43}\equiv$ $\displaystyle
2\mathcal{V}_{14}+\mathcal{V}_{16}+4\mathcal{V}_{24}+2\mathcal{V}_{26},$ (69)
$\displaystyle C_{44}\equiv$
$\displaystyle\mathcal{V}_{14}+2\mathcal{V}_{15}+2\mathcal{V}_{24}+4\mathcal{V}_{25},$
(70) $\displaystyle C_{45}\equiv$ $\displaystyle
2\mathcal{V}_{15}+\mathcal{V}_{16}+4\mathcal{V}_{25}+2\mathcal{V}_{26},$ (71)
$\displaystyle C_{52}\equiv$ $\displaystyle
2\mathcal{V}_{24}+4\mathcal{V}_{26}+\mathcal{V}_{34}+2\mathcal{V}_{36},$ (72)
$\displaystyle C_{53}\equiv$ $\displaystyle
4\mathcal{V}_{24}+2\mathcal{V}_{26}+2\mathcal{V}_{34}+\mathcal{V}_{36},$ (73)
$\displaystyle C_{54}\equiv$ $\displaystyle
2\mathcal{V}_{24}+4\mathcal{V}_{25}+\mathcal{V}_{34}+2\mathcal{V}_{35},$ (74)
$\displaystyle C_{55}\equiv$ $\displaystyle
4\mathcal{V}_{25}+2\mathcal{V}_{26}+2\mathcal{V}_{35}+\mathcal{V}_{36}.$ (75)

Therefore,

$\displaystyle\hbar V_{4}=$
$\displaystyle\frac{1}{N_{A}N_{B}}\left(C_{11}+\vec{\xi}\cdot\vec{F}_{A}+\vec{\nu}\cdot\vec{F}_{B}+\vec{\xi}~{}^{\text{T}}~{}\overleftrightarrow{F}\vec{\nu}\right).$
(76)

In summary, the interaction Hamiltonian in the computational basis is

$\displaystyle\hat{H}_{\text{int}}=\hbar\times\text{diag}\\{V_{1},V_{2},V_{3},V_{4}\\},$
(77)

where the matrix elements are

$\displaystyle\hbar V_{1}=$ $\displaystyle C_{11},$ (78) $\displaystyle\hbar
V_{2}=$
$\displaystyle\frac{1}{N_{B}}\left(C_{11}+\vec{\nu}\cdot\vec{F}_{B}\right),$
(79) $\displaystyle\hbar V_{3}=$
$\displaystyle\frac{1}{N_{A}}\left(C_{11}+\vec{\xi}\cdot\vec{F}_{A}\right),$
(80) $\displaystyle\hbar V_{4}=$
$\displaystyle\frac{1}{N_{A}N_{B}}\left(C_{11}+\vec{\xi}\cdot\vec{F}_{A}+\vec{\nu}\cdot\vec{F}_{B}+\vec{\xi}~{}^{\text{T}}~{}\overleftrightarrow{F}\vec{\nu}\right).$
(81)

###  4 Noise Simulations and Fidelity Calculations

Supplementary Figure 1: $1/f$ noise spectrum. (a) Typical numerically
generated $1/f$ noise spectrum ($\delta\varepsilon$, $\delta\varepsilon_{m}$
or $\delta t_{l(r)}$) with $15\times 10^{5}$ samples. (b) Power density
spectrum of a typical noise spectrum with red dashed line showing the $1/f$
dependence. The spectrum corresponds to that which is described in Eq. (83),
with higher and lower frequency cutoffs of 50 GHz and 66.7 kHz respectively.

In this section, we will explain how $1/f$ noise is generated and incorporated
into the Hamiltonian, and therefore forming the noisy gate, of which we can
then calculate the fidelity. The typical noise spectrum, Fig. 1(a) is
generated in Python, modified from code written by Felix Patzelt [7, 8] using
an algorithm in which we generate Gaussian distributed noise that is inversely
proportional to the number of samples. The power spectrum density, Fig. 1(b),
is also plotted to demonstrate that the noise has $1/f$ characteristics.

The noise is generated in the following way:

  1. 1.

First, based on the number of samples (i.e. number of data points) that is to
be generated, we calculate the frequency spectrum. In our simulations, we use
a total of $15\times 10^{5}$ samples for a time interval of
$[0,\tau_{\mathrm{expt}}]$, where $\tau_{\mathrm{expt}}$ is the total length
of time of the experiment. In terms of frequency, the sampling frequency will
be defined as $f_{\mathrm{samp}}=1/\Delta t$, where $\Delta t$ is the size of
the time steps as defined by the time evolution of the Hamiltonian. This means
that $\tau_{\mathrm{expt}}=N\Delta t$ can change depending on the size of the
time evolution time step. The higher frequency cutoff will consequentially be
$f_{h}=f_{\mathrm{samp}}/2$ and the lower frequency cutoff will be defined as
$f_{l}=f_{\mathrm{samp}}/N$ for all simulations except for results where the
value of $f_{l}$ is independently adjusted to examine its effect on fidelity.

  2. 2.

Second, calculate the scaling factors corresponding to each frequency. The
lower frequency cutoff is applied here with all lower values of frequency set
to the cutoff value. Following which, we construct the scaling factors, $s$,
in this way:

$\displaystyle s=f^{-\alpha/2}$ (82)

such that we account for the frequency characteristic of the spectrum as given
by $\alpha$ in the exponent of the frequency. Note that the exponent is
divided by $1/2$ because the noise spectrum is proportional to the square root
of the power spectrum density, $\sqrt{S(2\pi f)}$.

  3. 3.

Then, we calculate the ideal variance of the spectrum assuming a Gaussian
distribution. This determines the amplitude of the output noise spectrum and
will be proportional to the scaling factors calculated earlier.

  4. 4.

Next, based on the scaling, we generate random normal distributed values for
both the power and the phase of the noise spectrum.

  5. 5.

Finally, a Fourier transformation is performed on these randomly generated
values to obtain a real time series based on the previously generated noise
spectrum, which is then scaled to unit variance. The use of the scaling factor
ensures that the generated spectrum follows the intended frequency power law,
which is also shown in Fig. 1(b). The final noise spectrum can be described as
follows:

$\displaystyle S(2\pi f)=\begin{cases}\Delta^{2}_{n_{i}}/f_{l}&f\leq f_{l}\\\
\Delta^{2}_{n_{i}}/f&f_{l}\leq f\leq f_{h}\\\
0&\mathrm{,otherwise}\end{cases}$ (83)

In the simulations of time evolution, time steps of 0.01 ns is used,
indicating a sampling frequency, $f_{\mathrm{samp}}=100$ GHz. Correspondingly,
the higher frequency cutoff, $f_{h}$ is 50 GHz, and the lower frequency
cutoff, $f_{l}$ is 66.7 kHz.

These noise spectra will be generated over many iterations of noisy time
evolution, and we will now describe here how noise enters the Hamiltonian. In
each iteration of noisy time evolution, we generate two separate noise spectra
for four different noisy parameters:

$\displaystyle\\{\delta\varepsilon,\delta\varepsilon_{m},\delta t_{l},\delta
t_{r}\\}$ (84)

We need only two separate noise spectra for four parameters because we are
considering uncorrelated noise with independent spectra, so we are only
effectively comparing the effects of only two noisy parameters at the same
time which would then mean that only two independent noise spectra are needed
at the same time.

Depending on the parameter being considered, we consider a fixed magnitude of
fluctuations in the parameter, with differing magnitudes for the detuning and
tunneling parameters since they have different magnitudes. The output spectrum
generated from the procedure is scaled by a value of variance that is
calculated from the scaling factors which may not be of the amplitude that we
want. In addition, the mean of the generated spectrum is not exactly at zero.
Therefore, we want to first scale the noise spectra to unit variance and zero
mean:

$\displaystyle y_{i}=\frac{x_{i}-\bar{x}}{\bar{\sigma_{x}}}$ (85)

where $x$ is the generated noise spectrum and $y$ is the normalized noise
spectrum. After this normalization, we can now adjust freely the amplitude and
standard deviation of the noise spectrum. The spectra is then multiplied by
the desired noise amplitude (which changes effectively the standard deviation
of the noise) as shown in Eq. (84). We define the noise levels in these
parameters to vary from $10^{-5}$ to $10^{-3}$ for detuning parameters and
$10^{-6}$ to $10^{-4}$ for tunneling parameters. These noisy parameters are
now added to the corresponding parameters such that

$\displaystyle\\{\varepsilon\rightarrow\varepsilon+\delta\varepsilon,~{}\varepsilon_{m}\rightarrow\varepsilon_{m}+\delta\varepsilon_{m},~{}t_{l}\rightarrow
t_{l}+\delta t_{l},~{}t_{r}\rightarrow t_{r}+\delta t_{r}\\}$ (86)

Fluctuations in the detuning and tunnel coupling parameters will directly
result in fluctuations in the admixtures ($\alpha_{i}$ and $\beta_{i}$) of the
doubly occupied singlet states in the dressed singlet excited state
($\ket{1^{\prime}_{A}}$ and $\ket{1^{\prime}_{B}}$).

$\displaystyle\alpha^{\prime}_{1}$ $\displaystyle=\frac{t_{rA}+\delta
t_{rA}}{\sqrt{2}(-U+U^{\prime\prime}+2(\varepsilon_{A}+\delta\varepsilon_{A}))}$
(87) $\displaystyle\alpha^{\prime}_{2}$ $\displaystyle=\frac{t_{lA}+\delta
t_{lA}}{\sqrt{2}(-U+U^{\prime\prime}-2(\varepsilon_{A}+\delta\varepsilon_{A}))}$
(88) $\displaystyle\alpha^{\prime}_{3}$ $\displaystyle=\frac{t_{rA}+\delta
t_{rA}}{\sqrt{2}(-U+2U^{\prime}-U^{\prime\prime}+\varepsilon_{A}+\delta\varepsilon_{A}+\varepsilon_{mA}+\delta\varepsilon_{mA})}$
(89) $\displaystyle\alpha^{\prime}_{4}$ $\displaystyle=\frac{t_{lA}+\delta
t_{lA}}{\sqrt{2}(-U+2U^{\prime}-U^{\prime\prime}-\varepsilon_{A}-\delta\varepsilon_{A}+\varepsilon_{mA}+\delta\varepsilon_{mA})}$
(90) $\displaystyle\beta^{\prime}_{1}$ $\displaystyle=\frac{t_{rB}+\delta
t_{rB}}{\sqrt{2}(-U+U^{\prime\prime}+2(\varepsilon_{B}+\delta\varepsilon_{B}))}$
(91) $\displaystyle\beta^{\prime}_{2}$ $\displaystyle=\frac{t_{lB}+\delta
t_{lB}}{\sqrt{2}(-U+U^{\prime\prime}-2(\varepsilon_{B}+\delta\varepsilon_{B}))}$
(92) $\displaystyle\beta^{\prime}_{3}$ $\displaystyle=\frac{t_{rB}+\delta
t_{rB}}{\sqrt{2}(-U+2U^{\prime}-U^{\prime\prime}+\varepsilon_{B}+\delta\varepsilon_{B}+\varepsilon_{mB}+\delta\varepsilon_{mB})}$
(93) $\displaystyle\beta^{\prime}_{4}$ $\displaystyle=\frac{t_{lB}+\delta
t_{lB}}{\sqrt{2}(-U+2U^{\prime}-U^{\prime\prime}-\varepsilon_{B}-\delta\varepsilon_{B}+\varepsilon_{mB}+\delta\varepsilon_{mB})}$
(94)

This will lead to fluctuations in the Hamiltonian terms ($V_{1}$, $V_{2}$,
$V_{3}$, and $V_{4}$), such that the interaction Hamiltonian as outlined in
Eq. (77) is now $\mathrm{diag}\\{V_{1},V_{2}+\delta V_{2},V_{3}+\delta
V_{3},V_{4}+\delta V_{4}\\}$.

## 0emupplementary Notes

###  1 Comparison to the Ising Coupling

Since the interaction Hamiltonian is diagonal, it includes only the
$\Bqty{\hat{\mathbb{1}}_{A},\hat{\sigma}_{z,A}}\otimes\Bqty{\hat{\mathbb{1}}_{B},\hat{\sigma}_{z,B}}$
terms. We can recast the Hamiltonian in the Pauli basis by first writing the
Hamiltonian in terms of the two-qubit basis states, then rewriting the state
projectors with the Pauli operators

$\displaystyle\hat{H}_{\text{int}}$ $\displaystyle=\ket{0_{A}0_{B}}\hbar
V_{1}\bra{0_{A}0_{B}}+\ket{0_{A}1^{\prime}_{B}}\hbar
V_{2}\bra{0_{A}1^{\prime}_{B}}+\ket{1^{\prime}_{A}0_{B}}\hbar
V_{3}\bra{1^{\prime}_{A}0_{B}}+\ket{1^{\prime}_{A}1^{\prime}_{B}}\hbar
V_{4}\bra{1^{\prime}_{A}1^{\prime}_{B}}$ (95) $\displaystyle=\hbar
V_{1}\pqty\Big{\tfrac{1}{2}\pqty{\hat{\mathbb{1}}_{A}+\hat{\sigma}_{z,A}}\otimes\tfrac{1}{2}\pqty{\hat{\mathbb{1}}_{B}+\hat{\sigma}_{z,B}}}+\hbar
V_{2}\pqty\Big{\tfrac{1}{2}\pqty{\hat{\mathbb{1}}_{A}+\hat{\sigma}_{z,A}}\otimes\tfrac{1}{2}\pqty{\hat{\mathbb{1}}_{B}-\hat{\sigma}_{z,B}}}$
$\displaystyle\qquad+\hbar
V_{3}\pqty\Big{\tfrac{1}{2}\pqty{\hat{\mathbb{1}}_{A}-\hat{\sigma}_{z,A}}\otimes\tfrac{1}{2}\pqty{\hat{\mathbb{1}}_{B}+\hat{\sigma}_{z,B}}}+\hbar
V_{4}\pqty\Big{\tfrac{1}{2}\pqty{\hat{\mathbb{1}}_{A}-\hat{\sigma}_{z,A}}\otimes\tfrac{1}{2}\pqty{\hat{\mathbb{1}}_{B}-\hat{\sigma}_{z,B}}}$
$\displaystyle=\tfrac{1}{4}\pqty{\hbar V_{1}+\hbar V_{2}+\hbar V_{3}+\hbar
V_{4}}~{}\hat{\mathbb{1}}_{A}\otimes\hat{\mathbb{1}}_{B}+\tfrac{1}{4}\pqty{\hbar
V_{1}-\hbar V_{2}+\hbar V_{3}-\hbar
V_{4}}~{}\hat{\sigma}_{z,A}\otimes\hat{\mathbb{1}}_{B}$
$\displaystyle\qquad+\tfrac{1}{4}\pqty{\hbar V_{1}+\hbar V_{2}-\hbar
V_{3}-\hbar
V_{4}}~{}\hat{\mathbb{1}}_{A}\otimes\hat{\sigma}_{z,B}+\tfrac{1}{4}\pqty{\hbar
V_{1}-\hbar V_{2}-\hbar V_{3}+\hbar
V_{4}}~{}\hat{\sigma}_{z,A}\otimes\hat{\sigma}_{z,B}$

Upon the recognition that the $\hat{\sigma}_{z,A}\otimes\hat{\mathbb{1}}_{B}$
and $\hat{\mathbb{1}}_{A}\otimes\hat{\sigma}_{z,B}$ terms are single-qubit
operations that can be grouped into the single-qubit Hamiltonians as
adjustments to the single-qubit energies, and ignoring the global phase from
$\hat{\mathbb{1}}_{A}\otimes\hat{\mathbb{1}}_{B}$, it leaves the $Z$-$Z$
coupling term

$\displaystyle\tfrac{1}{4}\pqty{\hbar V_{1}-\hbar V_{2}-\hbar V_{3}+\hbar
V_{4}}~{}\hat{\sigma}_{z,A}\otimes\hat{\sigma}_{z,B}$ (96)
$\displaystyle\qquad=\frac{1}{4}\pqty{\pqty{1-\tfrac{1}{N_{A}}}\pqty{1-\tfrac{1}{N_{B}}}C_{11}-\pqty{1-\tfrac{1}{N_{A}}}\pqty{\tfrac{\vec{\nu}}{N_{B}}}\cdot\vec{F}_{B}-\pqty{1-\tfrac{1}{N_{B}}}\pqty{\tfrac{\vec{\xi}}{N_{A}}}\cdot\vec{F}_{A}+\pqty{\tfrac{\vec{\xi}}{N_{A}}}^{T}\overleftrightarrow{F}\pqty{\tfrac{\vec{\nu}}{N_{B}}}}~{}\hat{\sigma}_{z,A}\otimes\hat{\sigma}_{z,B}.$

This agrees with previous work done on capacitive coupling that approximates a
capacitively-coupled two-qubit interaction as an Ising Hamiltonian of the form
$\hat{H}^{(\text{Ising})}=-J\pqty{\hat{\sigma}_{z}\otimes\hat{\sigma}_{z}}$[9,
10].

###  2 Noise Autocorrelation Function

In this paper, the control parameters are taken to have $1/f$ noise. This has
the spectral density

$S_{n_{i}}(\omega)=\begin{cases}\Delta_{n_{i}}^{2}\frac{2\pi}{\omega_{l}},&\text{for
}\absolutevalue{\omega}\leq\omega_{l}\\\
\Delta_{n_{i}}^{2}\frac{2\pi}{\absolutevalue{\omega}},&\text{for
}\omega_{l}\leq\absolutevalue{\omega}\leq\omega_{h}\\\
0,&\text{otherwise}\end{cases},$ (97)

where $\Delta_{n_{i}}^{2}$ is the noise amplitude. The exact form of the auto-
correlation function $C_{n_{i}}(t)$ can be worked out using the
Wiener–Khinchin theorem [11] to give

$\displaystyle C_{n_{i}}(t)$
$\displaystyle=\int_{-\infty}^{\infty}\frac{\differential{\omega}}{2\pi}e^{i\omega
t}S_{n_{i}}(\omega)$ (98)
$\displaystyle=\Delta_{n_{i}}^{2}\pqty{\int_{-\omega_{l}}^{\omega_{l}}\differential{\omega}\frac{e^{i\omega
t}}{\omega_{l}}+\int_{\omega_{l}}^{\omega_{h}}\differential{\omega}\frac{e^{i\omega
t}}{\omega}+\int_{-\omega_{h}}^{-\omega_{l}}\differential{\omega}\frac{e^{i\omega
t}}{-\omega}}$
$\displaystyle=\Delta_{n_{i}}^{2}\pqty{\int_{0}^{\omega_{l}}\differential{\omega}\frac{e^{i\omega
t}+e^{-i\omega
t}}{\omega_{l}}+\int_{\omega_{l}}^{\omega_{h}}\differential{\omega}\frac{e^{i\omega
t}+e^{-i\omega t}}{\omega}}$
$\displaystyle=2\Delta_{n_{i}}^{2}\pqty{\int_{0}^{\omega_{l}}\differential{\omega}\frac{\cos(\omega
t)}{\omega_{l}}+\int_{\omega_{l}}^{\omega_{h}}\differential{\omega}\frac{\cos(\omega
t)}{\omega}}$
$\displaystyle=2\Delta_{n_{i}}^{2}\pqty{\frac{\sin(\omega_{l}t)}{\omega_{l}t}+\operatorname{Ci}\pqty{\omega_{h}t}-\operatorname{Ci}\pqty{\omega_{l}t}}.$

Meanwhile, the noise amplitude can be found by relating the auto-correlation
to the variance, $\lim_{t\to 0}C_{n_{i}}(t)=\sigma_{n_{i}}^{2}$, given by

$\displaystyle\sigma_{n_{i}}^{2}$ $\displaystyle=\lim_{t\to 0}C_{n_{i}}(t)$
(99) $\displaystyle=2\Delta_{n_{i}}^{2}\lim_{t\to
0}\pqty{\frac{\sin(\omega_{l}t)}{\omega_{l}t}+\operatorname{Ci}\pqty{\omega_{h}t}-\operatorname{Ci}\pqty{\omega_{l}t}}$
$\displaystyle=2\Delta_{n_{i}}^{2}\lim_{t\to
0}\pqty{1+\sum_{k=1}^{\infty}\pqty{-1}^{k}\frac{\pqty{\omega_{l}t}^{2k}}{\pqty{2k+1}!}+\pqty{\gamma+\ln(\omega_{h}t)+\sum_{k=1}^{\infty}\frac{\pqty{-\omega_{h}^{2}t^{2}}^{k}}{2k\pqty{2k}!}}-\pqty{\gamma+\ln(\omega_{l}t)+\sum_{k=1}^{\infty}\frac{\pqty{-\omega_{l}^{2}t^{2}}^{k}}{2k\pqty{2k}!}}}$
$\displaystyle=2\Delta_{n_{i}}^{2}\Bigg{(}1+\ln(\omega_{h}/\omega_{l})+\lim_{t\to
0}\sum_{k=1}^{\infty}\underbrace{\pqty{\pqty{-1}^{k}\frac{\pqty{\omega_{l}t}^{2k}}{\pqty{2k+1}!}+\frac{\pqty{-\omega_{h}^{2}t^{2}}^{k}}{2k\pqty{2k}!}-\frac{\pqty{-\omega_{l}^{2}t^{2}}^{k}}{2k\pqty{2k}!}}}_{\text{goes
to }0\text{ when }t\rightarrow 0}\Bigg{)}$
$\displaystyle=2\Delta_{n_{i}}^{2}\pqty{1+\ln(\omega_{h}/\omega_{l})}.$

Hence, the auto-correlation function is

$C_{n_{i}}(t)=\sigma_{n_{i}}^{2}\times\overline{C}(t)=\sigma_{n_{i}}^{2}\times\frac{1}{1+\ln(\omega_{h}/\omega_{l})}\pqty{\frac{\sin(\omega_{l}t)}{\omega_{l}t}+\operatorname{Ci}\pqty{\omega_{h}t}-\operatorname{Ci}\pqty{\omega_{l}t}},$
(100)

where $\overline{C}(t)$ is the normalized auto-correlation function for $1/f$
noise.

A common appearance in the analytical formula for fidelity is the double
integral
$\varsigma(t)=\int_{0}^{t}\differential{t^{\prime}}\int_{0}^{t^{\prime}}\differential{t^{\prime\prime}}\overline{C}(t)$,
which is

$\displaystyle\varsigma(t)$
$\displaystyle=\frac{1}{2\pqty{1+\ln(\omega_{h}/\omega_{l})}}\bigg{(}\operatorname{Ci}(\omega_{h}t)t^{2}-\operatorname{Ci}(\omega_{l}t)t^{2}-\frac{\sin(\omega_{h}t)t}{\omega_{h}}+\frac{\sin(\omega_{l}t)t}{\omega_{l}}+\frac{2\operatorname{Si}(\omega_{l}t)t}{\omega_{l}}+\frac{\cos(\omega_{h}t)-1}{\omega_{h}^{2}}+\frac{\cos(\omega_{l}t)-1}{\omega_{l}^{2}}\bigg{)}$
(101)

for $1/f$ noise. This can be easily differentiated twice to verify the result.

###  3 Analytic Expression for Makhlin Invariants

We use Makhlin invariants to determine if the interaction Hamiltonian is able
to form a two-qubit CPHASE or CNOT gate [12]. Given the interaction
Hamiltonian in Eq. (32), we find the two Makhlin invariants to be:

$\displaystyle G_{1}$
$\displaystyle=\frac{1}{4}\exp(-i(V_{1}+V_{2}+V_{3}+V_{4})t_{0})\times\left[\exp(i(V_{2}+V_{3})t_{0})+\exp(i(V_{1}+V_{4})t_{0})\right]^{2}$
$\displaystyle G_{2}$
$\displaystyle=2+\frac{1}{2}\exp(-i(V_{1}+V_{2}+V_{3}+V_{4})t_{0})\times\left[\exp(2i(V_{2}+V_{3})t_{0})+\exp(2i(V_{1}+V_{4})t_{0})\right]$

and there exists a value of gate time $t_{0}$ where $G_{1}=0$ and $G_{2}=1$,
which indicates that the two-qubit gate of interest is indeed a CNOT or CPHASE
gate up to a sequence of local gates. Expressing $G_{1}$ and $G_{2}$ in their
trigonometric identities, we find that the two-qubit gate time $t_{0}$ must
fulfill:

$\displaystyle t_{0}=\frac{k\pi}{(V_{1}-V_{2}-V_{3}+V_{4})}$ (102)

where $k$ is any odd integer value.

###  4 Equivalence of Fidelity Definitions

In this paper, we have used the operator fidelity for both the numerical
simulations and the analytical calculations. However, as unitary operators
cannot be directly measured in experiments, this figure of merit might not be
very appropriate for use with experimental data. To enable a more direct
comparison, and to ensure that our analysis is still applicable with
experimental considerations, we show in this section that the operator
fidelity is equivalent to a figure of merit more commonly used and directly
measureable in experiments.

#### The Pauli Basis

Every $n$-qubit operator $\hat{A}$ can be expanded into the Pauli basis

$\hat{A}=\sum_{i=0}^{d^{2}-1}a_{i}\hat{\sigma}_{i}^{(n)},$ (103)

where $d=2^{n}$ is the dimension of the $n$-qubit Hilbert space. The Pauli
basis operators for $n=1$ are the standard Pauli operators

$\displaystyle\hat{\sigma}_{0}^{(1)}$ $\displaystyle=\hat{\mathbb{1}}$
$\displaystyle\hat{\sigma}_{1}^{(1)}$ $\displaystyle=\hat{\sigma}_{x}$
$\displaystyle\hat{\sigma}_{2}^{(1)}$ $\displaystyle=\hat{\sigma}_{y}$
$\displaystyle\hat{\sigma}_{3}^{(1)}$ $\displaystyle=\hat{\sigma}_{z},$ (104)

while the basis operators for $n$ qubits are formed via the tensor product of
the single-qubit Pauli operators

$\hat{\sigma}_{i}^{(n)}\equiv\hat{\sigma}_{i_{1}}^{(1)}\otimes\hat{\sigma}_{i_{2}}^{(1)}\otimes\dots\otimes\hat{\sigma}_{i_{n}}^{(1)}.$
(105)

These Pauli operators form an orthogonal basis

$\tr(\hat{\sigma}_{i}^{(n)}\hat{\sigma}_{j}^{(n)})=d\delta_{ij},$ (106)

which, together with the linearity of the trace, can be used to find the
coefficients $a_{i}$ from the operator $\hat{A}$

$\displaystyle\frac{1}{d}\tr(\hat{\sigma}_{i}^{(n)}\hat{A})$
$\displaystyle=\frac{1}{d}\tr(\hat{\sigma}_{i}^{(n)}\sum_{j}a_{j}\hat{\sigma}_{j}^{(n)})$
(107)
$\displaystyle=\sum_{j}\frac{a_{j}}{d}\underbrace{\tr(\hat{\sigma}_{i}^{(n)}\hat{\sigma}_{j}^{(n)})}_{d\delta_{ij}}$
$\displaystyle=a_{i}.$

For the rest of the section, the superscripts of
$\hat{\sigma}_{\square}^{(n)}$ will be dropped with the understanding that the
appropriate Pauli operators will be used based on the number of qubits of the
system, the ideal unitary operator will be written as
$\hat{U}=\sum_{i}u_{i}\hat{\sigma}_{i}$, and the noisy unitary operator as
$\widetilde{U}=\sum_{i}\tilde{u}_{i}\hat{\sigma}_{i}$.

#### Operator Fidelity

The operator fidelity is defined [13] as

$\mathcal{F}\equiv\frac{1}{d^{2}}\absolutevalue{\tr(\hat{U}^{\dagger}\widetilde{U})}^{2},$
(108)

where $d$ is the dimension of the Hilbert space, $\hat{U}$ is the ideal
operator, and $\widetilde{U}$ is its noisy implementation. $\mathcal{F}=1$
only in the noiseless case where
$\widetilde{U}=\hat{U}\implies\hat{U}^{\dagger}\widetilde{U}=\hat{\mathbb{1}}$.
In the presence of noise, $\mathcal{F}<1$, and the fidelity decreases with an
increase in noise.

Following from the definition, the operator fidelity in terms of the Pauli
basis coefficients is

$\displaystyle\mathcal{F}$
$\displaystyle=\frac{1}{d^{2}}\absolutevalue{\tr(\hat{U}^{\dagger}\widetilde{U})}^{2}$
$\displaystyle=\frac{1}{d^{2}}\absolutevalue{\tr(\sum_{i}u_{i}^{*}\hat{\sigma}_{i}\sum_{j}\tilde{u}_{j}\hat{\sigma}_{j})}^{2}$
$\displaystyle=\frac{1}{d^{2}}\absolutevalue{\sum_{ij}u_{i}^{*}\tilde{u}_{j}\operatorname{tr}(\hat{\sigma}_{i}\hat{\sigma}_{j})}^{2}$
(109)
$\displaystyle=\frac{1}{d^{2}}\absolutevalue{\sum_{ij}u_{i}^{*}\tilde{u}_{j}d\delta_{ij}}^{2}$
$\displaystyle=\absolutevalue{\sum_{i}u_{i}^{*}\tilde{u}_{i}}^{2}$

#### Process Fidelity

In [14], the process fidelity is defined as

$\mathcal{F}^{(p)}=\tr(\chi\widetilde{\chi}),$ (110)

where the $\chi$ is the process matrix, whose elements $\chi_{ij}$ are the
coefficients that appear when expanding a quantum process in the Pauli basis

$\mathcal{E}\pqty{\hat{\rho}}=\sum_{ij}\chi_{ij}\,\hat{\sigma}_{i}\hat{\rho}\hat{\sigma}_{j}.$
(111)

For an ideal unitary map, this is

$\displaystyle\mathcal{U}\pqty{\hat{\rho}}$
$\displaystyle=\hat{U}\hat{\rho}\hat{U}^{\dagger}$ (112)
$\displaystyle=\pqty{\sum_{i}u_{i}\hat{\sigma}_{i}}\hat{\rho}\pqty{\sum_{j}u_{j}^{*}\hat{\sigma}_{j}}$
$\displaystyle=\sum_{ij}\underbrace{u_{i}u_{j}^{*}}_{\chi_{ij}}~{}\hat{\sigma}_{i}\hat{\rho}\hat{\sigma}_{j}.$

Similarly, the noisy unitary map has the process matrix elements
$\tilde{\chi}_{ij}=\tilde{u}_{i}\tilde{u}_{j}^{*}$. Placing this back into Eq.
(110),

$\displaystyle\mathcal{F}^{(p)}$ $\displaystyle=\tr(\chi\widetilde{\chi})$
$\displaystyle=\sum_{ij}\chi_{ij}\widetilde{\chi}_{ji}$
$\displaystyle=\sum_{ij}u_{i}u_{j}^{*}\tilde{u}_{j}\tilde{u}_{i}^{*}$ (113)
$\displaystyle=\pqty{\sum_{j}u_{j}^{*}\tilde{u}_{j}}\pqty{\sum_{i}u_{i}\tilde{u}_{i}^{*}}$
$\displaystyle=\absolutevalue{\sum_{i}u_{i}^{*}\tilde{u}_{i}}^{2}$

Therefore, we show that the operator fidelity and process fidelity are
equivalent.

#### Experimentally Measured Fidelity and the Ensemble Average

In actual setups, the fidelity is obtained by first measuring the process
matrix elements through an experimental procedure [15, 16, 17]. As the
experimentally-determined process matrix requires the averaging of repeated
measurements of the matrix elements to obtain an expected value of
$\expectationvalue{\widetilde{\chi}_{ij}}$, what can actually be measured in
an experiment is in fact

$\mathcal{F}^{(\text{exp})}=\tr(\chi\\!\expectationvalue{\widetilde{\chi}})=\expectationvalue{\tr(\chi\widetilde{\chi})}=\expectationvalue{\mathcal{F}},$
(114)

which is the average fidelity. The cumulant average can be brought out of the
trace as the trace is a linear operation.

Since the two definitions of fidelity are equivalent, we report the operator
fidelity Eq. (108), averaged over noise realizations, to emulate the result of
experimentally-determined fidelities.

###  5 Analytical Derivation of Fidelity

#### Time Evolution in the Presence of Noise

Consider the control parameters
$\mathbf{n}_{A/B}\equiv(t_{lA/B},~{}t_{rA/B},~{}\varepsilon_{A/B},~{}\varepsilon_{mA/B})^{T}$.
Note that this notation is slightly different from the main text, where the
parameters are instead combined into a single vector
$\vec{n}=(\mathbf{n}_{A},~{}\mathbf{n}_{B})^{T}$. In the presence of noise,
the noisy parameters
$\tilde{\mathbf{n}}_{A/B}(t)=\mathbf{n}_{A/B}+\delta\mathbf{n}_{A/B}(t)$
introduce fluctuations into the interaction Hamiltonian. We take these
fluctuations to have zero mean as we are assuming that upon many repetitions
of the experiment, the noisy parameters average out to the target (ideal)
parameters $\expectationvalue{\tilde{\mathbf{n}}_{A/B}(t)}=\mathbf{n}_{A/B}$.

The resulting interaction due to the noise can be written as the ideal
Hamiltonian with a noise term

$\widetilde{H}_{\text{int}}=\hat{H}_{\text{int}}+\delta\hat{H}_{\text{int}}.$
(115)

As we do not expect the magnitude of the noise to be very large in comparison
to the control parameters, $\delta\hat{H}_{\text{int}}$ can be found by taking
the first term of the series expansion of $\hat{H}$. This means that

$\delta\hat{H}_{\text{int}}\approx\gradient_{\\!A}\hat{H}_{\text{int}}~{}\delta\mathbf{n}_{A}+\gradient_{\\!B}\hat{H}_{\text{int}}~{}\delta\mathbf{n}_{B}$
(116)

where
$\gradient_{A/B}=(\partial_{t_{lA/B}},~{}\partial_{t_{rA/B}},~{}\partial_{\varepsilon_{A/B}},~{}\partial_{\varepsilon_{mA/B}})$.
In addition, only the noise terms with coefficients to the linear order of
$t_{l/r,A/B}/(U+\dots)$ are kept. This is because the tunneling parameters are
much smaller in magnitude than the inter/intra-dot energies when operating a
triple quantum dot qubit.

Writing the normalization factor as $N_{A}=1+\vec{\xi}\cdot\vec{1}$ and
$N_{B}=1+\vec{\nu}\cdot\vec{1}$ where $\vec{1}=\\{1,1,1,1\\}^{T}$, the last
diagonal term in $\delta\hat{H}$ is

$\displaystyle\hbar\delta V_{4}$
$\displaystyle\approx\gradient_{\\!A}\bqty{\frac{C_{11}+\vec{\xi}\cdot\vec{F}_{A}+\vec{\nu}\cdot\vec{F}_{B}+{\vec{\xi}}^{T}\overleftrightarrow{F}\vec{\nu}}{N_{A}N_{B}}}\delta\mathbf{n}_{A}+\gradient_{\\!B}\bqty{\frac{C_{11}+\vec{\xi}\cdot\vec{F}_{A}+\vec{\nu}\cdot\vec{F}_{B}+{\vec{\xi}}^{T}\overleftrightarrow{F}\vec{\nu}}{N_{A}N_{B}}}\delta\mathbf{n}_{B}$
$\displaystyle=\frac{1}{N_{A}N_{B}}\Big{(}{\vec{F}_{A}}^{T}\gradient_{\\!A}\vec{\xi}~{}\delta\mathbf{n}_{A}+{\vec{F}_{B}}^{T}\gradient_{\\!B}\vec{\nu}~{}\delta\mathbf{n}_{B}+{\vec{\xi}}^{T}\overleftrightarrow{F}\gradient_{\\!B}\vec{\nu}~{}\delta\mathbf{n}_{B}+{\vec{\nu}}^{T}\overleftrightarrow{F}^{T}\gradient_{\\!A}\vec{\xi}~{}\delta\mathbf{n}_{A}$
$\displaystyle\qquad-\pqty{N_{B}\gradient_{\\!A}N_{A}~{}\delta\mathbf{n}_{A}+N_{A}\gradient_{\\!B}N_{B}~{}\delta\mathbf{n}_{B}}~{}\hbar
V_{4}\Big{)}$
$\displaystyle=\frac{1}{N_{A}N_{B}}\bigg{(}{\vec{F}_{A}}^{T}\gradient_{\\!A}\vec{\xi}~{}\delta\mathbf{n}_{A}+{\vec{F}_{B}}^{T}\gradient_{\\!B}\vec{\nu}~{}\delta\mathbf{n}_{B}+{\vec{\xi}}^{T}\overleftrightarrow{F}\gradient_{\\!B}\vec{\nu}~{}\delta\mathbf{n}_{B}+{\vec{\nu}}^{T}\overleftrightarrow{F}^{T}\gradient_{\\!A}\vec{\xi}~{}\delta\mathbf{n}_{A}$
$\displaystyle\qquad-\pqty{\pqty{1+\vec{\nu}\cdot\vec{1}}\pqty{\vec{1}^{T}\gradient_{\\!A}\vec{\xi}}~{}\delta\mathbf{n}_{A}+\pqty{1+\vec{\xi}\cdot\vec{1}}\pqty{\vec{1}^{T}\gradient_{\\!B}\vec{\nu}}~{}\delta\mathbf{n}_{B}}~{}\hbar
V_{4}\bigg{)}$
$\displaystyle=\frac{1}{N_{A}N_{B}}\bigg{(}{\vec{F}_{A}}^{T}\gradient_{\\!A}\vec{\xi}~{}\delta\mathbf{n}_{A}+{\vec{F}_{B}}^{T}\gradient_{\\!B}\vec{\nu}~{}\delta\mathbf{n}_{B}-\pqty{{\vec{1}^{T}\gradient_{\\!A}\vec{\xi}}~{}\delta\mathbf{n}_{A}+{\vec{1}^{T}\gradient_{\\!B}\vec{\nu}}~{}\delta\mathbf{n}_{B}}~{}\hbar
V_{4}$
$\displaystyle\qquad+{\vec{\xi}}^{T}\overleftrightarrow{F}\gradient_{\\!B}\vec{\nu}\;\delta\mathbf{n}_{B}+{\vec{\nu}}^{T}\overleftrightarrow{F}^{T}\gradient_{\\!A}\vec{\xi}\;\delta\mathbf{n}_{A}-\pqty{\pqty{\vec{\nu}\cdot\vec{1}}\pqty{{\vec{1}}^{T}\gradient_{\\!A}\vec{\xi}}\delta\mathbf{n}_{A}+\pqty{\vec{\xi}\cdot\vec{1}}\pqty{{\vec{1}}^{T}\gradient_{\\!B}\vec{\nu}}\delta\mathbf{n}_{B}}~{}\hbar
V_{4}\bigg{)}.$ (117)

The second row of the final expression clearly has terms larger than
$\mathcal{O}(\xi,\nu)=\mathcal{O}(\alpha^{2},\beta^{2})$, so they will be
dropped. Meanwhile, since the first row already includes linear-ordered terms,
the normalisation also has to be expanded to
$1/N_{A/B}=1-\mathcal{O}(\alpha^{2},\beta^{2})$. Finally, the only term that
survives from $\hbar V_{4}$ on the right-hand side of the equation is the
constant $C_{11}$, as it is multiplied by a factor that is already of the
order $\mathcal{O}(\alpha,\beta)$. Therefore, $\hbar\delta V_{4}$ is

$\hbar\delta
V_{4}=\pqty{\vec{F}_{A}-C_{11}\vec{1}}^{T}\gradient_{\\!A}\vec{\xi}~{}\delta\mathbf{n}_{A}+\pqty{\vec{F}_{B}-C_{11}\vec{1}}^{T}\gradient_{\\!B}\vec{\nu}~{}\delta\mathbf{n}_{B}+\mathcal{O}(\alpha^{2},\beta^{2})$
(118)

Applying a similar process to the remaining matrix elements,

$\displaystyle\hbar\delta V_{1}$ $\displaystyle\approx 0$ (119)
$\displaystyle\hbar\delta V_{2}$
$\displaystyle\approx\pqty{\vec{F}_{B}-C_{11}\vec{1}}^{T}\gradient_{\\!B}\vec{\nu}~{}\delta\mathbf{n}_{B}$
(120) $\displaystyle\hbar\delta V_{3}$
$\displaystyle\approx\pqty{\vec{F}_{A}-C_{11}\vec{1}}^{T}\gradient_{\\!A}\vec{\xi}~{}\delta\mathbf{n}_{A}$
(121) $\displaystyle\hbar\delta V_{4}$
$\displaystyle\approx\pqty{\vec{F}_{A}-C_{11}\vec{1}}^{T}\gradient_{\\!A}\vec{\xi}~{}\delta\mathbf{n}_{A}+\pqty{\vec{F}_{B}-C_{11}\vec{1}}^{T}\gradient_{\\!B}\vec{\nu}~{}\delta\mathbf{n}_{B}.$
(122)

The time evolution of the system in the presence of noise can be then found
from the noisy Hamiltonian.

$\displaystyle\widetilde{U}(t)$
$\displaystyle=\mathcal{T}_{+}\exp\Bqty{-\frac{i}{\hbar}\int_{0}^{t}\differential{t^{\prime}}\widetilde{H}_{\text{int}}(t^{\prime})}$
(123)
$\displaystyle=\underbrace{\exp\Bqty{-\frac{it}{\hbar}\hat{H}_{\text{int}}}}_{\equiv\hat{U}(t)}\underbrace{\exp\Bqty{-\frac{i}{\hbar}\int_{0}^{t}\differential{t^{\prime}}\delta\hat{H}_{\text{int}}(t^{\prime})}}_{\equiv\hat{U}_{\delta}(t)}$

Here, $\hat{U}(t)$ is the ideal time evolution of the system in the absence of
noise, while $\hat{U}_{\delta}(t)$ consists of just the contributions due to
noise. Note that the second step can be arrived from the first step only
because

  1. 1.

The interaction Hamiltonian, its ideal term, and its noise term are all
diagonal. So, the time-ordering operator is not required as they commute among
themselves, and for the same reason, the sum in the exponent can be seperated
into a product of exponentials

  2. 2.

The control parameters are constant over the operation of the two-qubit gate,
so the only time-dependent terms are their fluctuations

#### Cumulant Expansion of Fidelity

Considering the operator fidelity, Eq. (108) in the context of Eq. (123), the
fidelity of the noisy two-qubit interaction is

$\displaystyle\mathcal{F}$
$\displaystyle=\frac{1}{4^{2}}\absolutevalue{\tr(\hat{U}^{\dagger}\widetilde{U})}^{2}$
(124)
$\displaystyle=\frac{1}{16}\absolutevalue{\tr(\hat{U}^{\dagger}\pqty{\hat{U}\hat{U}_{\delta}})}^{2}$
$\displaystyle=\frac{1}{16}\absolutevalue{\tr(\hat{U}_{\delta})}^{2}.$

In an actual setup, this fidelity will be obtained through the repetition of
many experiments — or in other words, over many noise realizations. This means
that the final value obtained will be a cumulant average

$\displaystyle\expectationvalue{\mathcal{F}}$
$\displaystyle=\frac{1}{16}\expectationvalue{\absolutevalue{\tr(\exp\Bqty{-\frac{i}{\hbar}\int_{0}^{t}\differential{t^{\prime}}\delta\hat{H}_{\text{int}}(t^{\prime})})}^{2}}$
(125)
$\displaystyle=\frac{1}{16}\expectationvalue{\absolutevalue{\sum_{i}\exp\Bqty{-i\int_{0}^{t}\differential{t^{\prime}}\delta
V_{i}(t^{\prime})}}^{2}}$
$\displaystyle=\frac{1}{16}\sum_{i,j}\expectationvalue{\exp\Bqty{-i\int_{0}^{t}\differential{t^{\prime}}\pqty{\delta
V_{i}(t^{\prime})-\delta V_{j}(t^{\prime})}}}.$

This averaged operator fidelity is the same as the experimentally-measurable
fidelity — this equivalence is laid out in Sec. 4.

From the cumulant expansion [18], to the first order,

$\displaystyle\expectationvalue{\mathcal{F}}$
$\displaystyle\approx\frac{1}{16}\sum_{i,j}\exp\Bqty{-i\int_{0}^{t}\differential{t^{\prime}}\underbrace{\expectationvalue{\delta
V_{i}(t^{\prime})-\delta
V_{j}(t^{\prime})}}_{0}-\frac{1}{2}\int_{0}^{t}\differential{t^{\prime}}\int_{0}^{t^{\prime}}\differential{t^{\prime\prime}}\expectationvalue{\pqty{\delta
V_{i}(0)-\delta V_{j}(0)}\pqty{\delta V_{i}(t^{\prime\prime})-\delta
V_{j}(t^{\prime\prime})}}}$ (126) $\displaystyle=\frac{1}{16}\sum_{i\neq
j}\exp\Bqty{-\frac{1}{2}\int_{0}^{t}\differential{t^{\prime}}\int_{0}^{t^{\prime}}\differential{t^{\prime\prime}}\expectationvalue{\delta
V_{i}(0)\delta V_{i}(t^{\prime\prime})}-2\expectationvalue{\delta
V_{i}(0)\delta V_{j}(t^{\prime\prime})}+\expectationvalue{\delta
V_{j}(0)\delta V_{j}(t^{\prime\prime})}},$

where the first term vanishes as the noise has zero mean. In addition, it is
assumed that the correlation is only defined by the time difference between
the two factors, so $\expectationvalue{\delta V_{i}(0)\delta
V_{j}(t^{\prime\prime})}=\expectationvalue{\delta
V_{i}(t^{\prime\prime})\delta V_{j}(0)}$.

This expression can be further simplified by substituting $\delta V_{1}=0$ and
$\delta V_{4}=\delta V_{2}+\delta V_{3}$ from Eq. (119) to Eq. (122)

$\displaystyle\expectationvalue{\mathcal{F}}$
$\displaystyle=\frac{1}{4}+\frac{1}{4}\exp\Bqty{-\frac{1}{2}\int_{0}^{t}\differential{t^{\prime}}\int_{0}^{t^{\prime}}\differential{t^{\prime\prime}}\expectationvalue{\delta
V_{2}(0)\delta
V_{2}(t^{\prime\prime})}}+\frac{1}{4}\exp\Bqty{-\frac{1}{2}\int_{0}^{t}\differential{t^{\prime}}\int_{0}^{t^{\prime}}\differential{t^{\prime\prime}}\expectationvalue{\delta
V_{3}(0)\delta V_{3}(t^{\prime\prime})}}$ (127)
$\displaystyle\qquad~{}+\frac{1}{8}\exp\Bqty{-\frac{1}{2}\int_{0}^{t}\differential{t^{\prime}}\int_{0}^{t^{\prime}}\differential{t^{\prime\prime}}\expectationvalue{\delta
V_{2}(0)\delta V_{2}(t^{\prime\prime})}-2\expectationvalue{\delta
V_{2}(0)\delta V_{3}(t^{\prime\prime})}+\expectationvalue{\delta
V_{3}(0)\delta V_{3}(t^{\prime\prime})}}$
$\displaystyle\qquad~{}+\frac{1}{8}\exp\Bqty{-\frac{1}{2}\int_{0}^{t}\differential{t^{\prime}}\int_{0}^{t^{\prime}}\differential{t^{\prime\prime}}\expectationvalue{\delta
V_{2}(0)\delta V_{2}(t^{\prime\prime})}+2\expectationvalue{\delta
V_{2}(0)\delta V_{3}(t^{\prime\prime})}+\expectationvalue{\delta
V_{3}(0)\delta V_{3}(t^{\prime\prime})}}.$

Finally, $\expectationvalue{\delta V_{2}(0)\delta V_{3}(t^{\prime\prime})}$
can be taken to be zero if there is no correlation between the noise in the
two qubits, as $\delta V_{2}$ only includes terms in qubit B and $\delta
V_{3}$ only includes terms in qubit A

$\displaystyle\expectationvalue{\mathcal{F}}$
$\displaystyle=\frac{1}{4}+\frac{1}{4}\exp\Bqty{-\frac{1}{2}\int_{0}^{t}\differential{t^{\prime}}\int_{0}^{t^{\prime}}\differential{t^{\prime\prime}}\expectationvalue{\delta
V_{2}(0)\delta
V_{2}(t^{\prime\prime})}}+\frac{1}{4}\exp\Bqty{-\frac{1}{2}\int_{0}^{t}\differential{t^{\prime}}\int_{0}^{t^{\prime}}\differential{t^{\prime\prime}}\expectationvalue{\delta
V_{3}(0)\delta V_{3}(t^{\prime\prime})}}$ (128)
$\displaystyle\qquad~{}+\frac{1}{4}\exp\Bqty{-\frac{1}{2}\int_{0}^{t}\differential{t^{\prime}}\int_{0}^{t^{\prime}}\differential{t^{\prime\prime}}\expectationvalue{\delta
V_{2}(0)\delta V_{2}(t^{\prime\prime})}+\expectationvalue{\delta
V_{3}(0)\delta V_{3}(t^{\prime\prime})}}$
$\displaystyle=\frac{1}{4}\pqty{1+\exp\Bqty{-\frac{1}{2}\int_{0}^{t}\differential{t^{\prime}}\int_{0}^{t^{\prime}}\differential{t^{\prime\prime}}\expectationvalue{\delta
V_{2}(0)\delta
V_{2}(t^{\prime\prime})}}}\pqty{1+\exp\Bqty{-\frac{1}{2}\int_{0}^{t}\differential{t^{\prime}}\int_{0}^{t^{\prime}}\differential{t^{\prime\prime}}\expectationvalue{\delta
V_{3}(0)\delta V_{3}(t^{\prime\prime})}}}$

As the expression for $\delta V_{2}$ and $\delta V_{3}$ are identical, it
suffices to work out just one of the expressions and apply it to the other by
switching $A\leftrightarrow B$ and $\vec{\xi}\leftrightarrow\vec{\nu}$.
However, note that this similarity is purely symbolic, in the sense that
$\vec{F}_{A}$ and $\vec{F}_{B}$ do not share the same numerical values in
general, so care must be taken to ensure that $\vec{F}_{A}$ always appears
alongside $\vec{\xi}$ in $\delta V_{3}$ and $\vec{F}_{B}$ alongside
$\vec{\nu}$ in $\delta V_{2}$. With that in mind, the qubit subscripts will be
neglected in the intermediate steps.

$\displaystyle\expectationvalue{\delta V(0)\delta V(t^{\prime\prime})}$
$\displaystyle=\expectationvalue{\pqty{\vec{F}-C_{11}\vec{1}}^{T}\gradient\vec{\xi}~{}\delta\mathbf{n}(0)\pqty{\vec{F}-C_{11}\vec{1}}^{T}\gradient\vec{\xi}~{}\delta\mathbf{n}(t^{\prime\prime})}$
(129)
$\displaystyle=\sum_{i,j}\expectationvalue{\pqty{\vec{F}-C_{11}\vec{1}}^{T}\partial_{i}\vec{\xi}~{}\delta
n_{i}(0)\;\pqty{\vec{F}-C_{11}\vec{1}}^{T}\partial_{j}\vec{\xi}~{}\delta
n_{j}(t^{\prime\prime})}$
$\displaystyle=\sum_{i,j}\pqty{\pqty{\vec{F}-C_{11}\vec{1}}^{T}\partial_{i}\vec{\xi}}\pqty{\pqty{\vec{F}-C_{11}\vec{1}}^{T}\partial_{j}\vec{\xi}}\expectationvalue{\delta
n_{i}(0)\delta n_{j}(t^{\prime\prime})}$

Here, two further simplifying assumptions are made:

  1. 1.

The noise in the control parameters are independent of each other. Hence,

$\expectationvalue{\delta n_{i}(0)\delta
n_{j}(t)}=\delta_{ij}\expectationvalue{\delta n_{i}(0)\delta n_{i}(t)}$ (130)

where $\delta_{ij}$ is the Kronecker delta

  2. 2.

Every noise statistics follow the same behaviour, in the sense that apart from
the scale, the noise distributions are identical. So,

$\expectationvalue{\delta n_{i}(0)\delta
n_{i}(t)}=C_{n_{i}}(t)=\sigma_{n_{i}}^{2}\overline{C}(t)$ (131)

where $\sigma_{n_{i}}$ is the standard deviation of the noise in control
parameter $n_{i}$, and $\overline{C}(t)$ is the autocorrelation function
normalized to $\overline{C}(0)=1$, so that $C_{n_{i}}(0)=\sigma_{n_{i}}^{2}$

If necessary, these assumptions can be relaxed by assigning an auto-
correlation function $\overline{C}_{n_{i},n_{j}}(t)$ to the (auto-)correlation
of parameter(s). For the purposes of our analysis, however, the
simplifications above suffices to investigate the effects of noise on the
fidelity of our system.

With that,

$\displaystyle\expectationvalue{\delta V(0)\delta V(t^{\prime\prime})}$
$\displaystyle=\sum_{i,j}\pqty{\pqty{\vec{F}-C_{11}\vec{1}}^{T}\partial_{i}\vec{\xi}}\pqty{\pqty{\vec{F}-C_{11}\vec{1}}^{T}\partial_{j}\vec{\xi}}\delta_{ij}\sigma_{n_{i}}^{2}C(t)$
(132)
$\displaystyle=\sum_{i}\pqty{\pqty{\vec{F}-C_{11}\vec{1}}^{T}\partial_{i}\vec{\xi}}^{2}\sigma_{n_{i}}^{2}C(t)$

Therefore, the fidelity is

$\displaystyle\expectationvalue{\mathcal{F}}$
$\displaystyle=\frac{1}{4}\pqty{1+\exp\Bqty{-\frac{1}{2}\sum_{i}\pqty{\pqty{\vec{F}_{A}-C_{11}\vec{1}}^{T}\partial_{i}\vec{\xi}}^{2}\sigma_{n_{iA}}^{2}\varsigma(t)}}\pqty{1+\exp\Bqty{-\frac{1}{2}\sum_{i}\pqty{\pqty{\vec{F}_{B}-C_{11}\vec{1}}^{T}\partial_{i}\vec{\nu}}^{2}\sigma_{n_{iB}}^{2}\varsigma(t)}}$
(133)

which leaves the time dependence of the fidelity purely in
$\varsigma(t)\equiv\int_{0}^{t}\differential{t^{\prime}}\int_{0}^{t^{\prime}}\differential{t^{\prime\prime}}C(t^{\prime\prime})$,
the double integral of the autocorrelation function.

This final form is instructive in approaching the optimization of the control
parameters to maximize the fidelity of the two-qubit gate. As $\sigma_{n_{i}}$
and $C(t)$ arise from the noise and cannot be controlled, its coefficients
${\vec{F}_{A}}^{T}\partial_{i}\vec{\xi}$ and
${\vec{F}_{B}}^{T}\partial_{i}\vec{\nu}$ are the parameters to be minimized
instead.

These coefficients for qubit A are, explicitly,

$\displaystyle\pqty{\vec{F}_{A}-C_{11}\vec{1}}^{T}\partial_{t_{lA}}\vec{\xi}$
$\displaystyle=\frac{t_{lA}\pqty{C_{31}-C_{11}}}{\pqty{U-U^{\prime\prime}+2\varepsilon_{A}}^{2}}+\frac{t_{lA}\pqty{C_{51}-C_{11}}}{\pqty{U-2U^{\prime}+U^{\prime\prime}-\varepsilon_{mA}+\varepsilon_{A}}^{2}}$
(134a)
$\displaystyle\pqty{\vec{F}_{A}-C_{11}\vec{1}}^{T}\partial_{t_{rA}}\vec{\xi}$
$\displaystyle=\frac{t_{rA}\pqty{C_{21}-C_{11}}}{\pqty{U-U^{\prime\prime}-2\varepsilon_{A}}^{2}}+\frac{t_{rA}\pqty{C_{41}-C_{11}}}{\pqty{U-2U^{\prime}+U^{\prime\prime}-\varepsilon_{mA}-\varepsilon_{A}}^{2}}$
(134b)
$\displaystyle\pqty{\vec{F}_{A}-C_{11}\vec{1}}^{T}\partial_{\varepsilon_{A}}\vec{\xi}$
$\displaystyle=\frac{2t_{rA}^{2}\pqty{C_{21}-C_{11}}}{\pqty{U-U^{\prime\prime}-2\varepsilon_{A}}^{3}}-\frac{2t_{lA}^{2}\pqty{C_{31}-C_{11}}}{\pqty{U-U^{\prime\prime}+2\varepsilon_{A}}^{3}}$
$\displaystyle\qquad\quad+\frac{t_{rA}^{2}\pqty{C_{41}-C_{11}}}{\pqty{U-2U^{\prime}+U^{\prime\prime}-\varepsilon_{mA}-\varepsilon_{A}}^{3}}-\frac{t_{lA}^{2}\pqty{C_{51}-C_{11}}}{\pqty{U-2U^{\prime}+U^{\prime\prime}-\varepsilon_{mA}+\varepsilon_{A}}^{3}}$
(134c)
$\displaystyle\pqty{\vec{F}_{A}-C_{11}\vec{1}}^{T}\partial_{\varepsilon_{mA}}\vec{\xi}$
$\displaystyle=\frac{t_{rA}^{2}\pqty{C_{41}-C_{11}}}{\pqty{U-2U^{\prime}+U^{\prime\prime}-\varepsilon_{mA}-\varepsilon_{A}}^{3}}+\frac{t_{lA}^{2}\pqty{C_{51}-C_{11}}}{\pqty{U-2U^{\prime}+U^{\prime\prime}-\varepsilon_{mA}+\varepsilon_{A}}^{3}}$
(134d)

and for qubit B,

$\displaystyle\pqty{\vec{F}_{B}-C_{11}\vec{1}}^{T}\partial_{t_{lB}}\vec{\xi}$
$\displaystyle=\frac{t_{lB}\pqty{C_{13}-C_{11}}}{\pqty{U-U^{\prime\prime}+2\varepsilon_{B}}^{2}}+\frac{t_{lB}\pqty{C_{15}-C_{11}}}{\pqty{U-2U^{\prime}+U^{\prime\prime}-\varepsilon_{mB}+\varepsilon_{B}}^{2}}$
(135a)
$\displaystyle\pqty{\vec{F}_{B}-C_{11}\vec{1}}^{T}\partial_{t_{rB}}\vec{\xi}$
$\displaystyle=\frac{t_{rB}\pqty{C_{12}-C_{11}}}{\pqty{U-U^{\prime\prime}-2\varepsilon_{B}}^{2}}+\frac{t_{rB}\pqty{C_{14}-C_{11}}}{\pqty{U-2U^{\prime}+U^{\prime\prime}-\varepsilon_{mB}-\varepsilon_{B}}^{2}}$
(135b)
$\displaystyle\pqty{\vec{F}_{B}-C_{11}\vec{1}}^{T}\partial_{\varepsilon_{B}}\vec{\xi}$
$\displaystyle=\frac{2t_{rB}^{2}\pqty{C_{12}-C_{11}}}{\pqty{U-U^{\prime\prime}-2\varepsilon_{B}}^{3}}-\frac{2t_{lB}^{2}\pqty{C_{13}-C_{11}}}{\pqty{U-U^{\prime\prime}+2\varepsilon_{B}}^{3}}$
$\displaystyle\qquad\quad+\frac{t_{rB}^{2}\pqty{C_{14}-C_{11}}}{\pqty{U-2U^{\prime}+U^{\prime\prime}-\varepsilon_{mB}-\varepsilon_{B}}^{3}}-\frac{t_{lB}^{2}\pqty{C_{15}-C_{11}}}{\pqty{U-2U^{\prime}+U^{\prime\prime}-\varepsilon_{mB}+\varepsilon_{B}}^{3}}$
(135c)
$\displaystyle\pqty{\vec{F}_{B}-C_{11}\vec{1}}^{T}\partial_{\varepsilon_{mB}}\vec{\xi}$
$\displaystyle=\frac{t_{rB}^{2}\pqty{C_{14}-C_{11}}}{\pqty{U-2U^{\prime}+U^{\prime\prime}-\varepsilon_{mB}-\varepsilon_{B}}^{3}}+\frac{t_{lB}^{2}\pqty{C_{15}-C_{11}}}{\pqty{U-2U^{\prime}+U^{\prime\prime}-\varepsilon_{mB}+\varepsilon_{B}}^{3}}$
(135d)

###  6 Two-Qubit Sweet Spots

As defined in the main text, two-qubit sweet spots (2QSS) are points in the
parameter space where the first-order derivatives (Eqs. (119)–(122)) vanish.
While the 2QSS do not arise from the maximization of the fidelity, the same
first-order derivatives also appear as coefficients in the exponent of the
fidelity in Eqs. (134) & (135). Hence, when operating at the sweet spot, the
fidelity is also maximized as a consequence.

#### Tunneling Sweet Spots

Before attempting to find the sweet spots, we shall take a closer look at the
coefficients $\vec{F}_{A/B}-C_{11}\vec{1}$

$\displaystyle\vec{F}_{A}-C_{11}\vec{1}$ $\displaystyle=\pmqty{undef}$ (136)
$\displaystyle=\sum_{j=4}^{6}\pmqty{undef}$
$\displaystyle=\sum_{j=4}^{6}\pmqty{undef}$ $\implies\begin{aligned}
C_{31}-C_{11}&=-\pqty{C_{51}-C_{11}},\\\
C_{41}-C_{11}&=-\pqty{C_{21}-C_{11}}.\end{aligned}$ (137)

Following the same steps with $\vec{F}_{B}$ leads us to

$\displaystyle C_{12}-C_{11}$ $\displaystyle=-\pqty{C_{14}-C_{11}},$ (138)
$\displaystyle C_{15}-C_{11}$ $\displaystyle=-\pqty{C_{13}-C_{11}}.$

With this in hand, Eqs. (134a) & (134b) give

$\displaystyle\pqty{\vec{F}_{A}-C_{11}\vec{1}}^{T}\partial_{t_{lA}}\vec{\xi}$
$\displaystyle=-t_{lA}\pqty{C_{51}-C_{11}}\pqty{\frac{1}{\pqty{U-U^{\prime\prime}+2\varepsilon_{A}}^{2}}-\frac{1}{\pqty{U-2U^{\prime}+U^{\prime\prime}-\varepsilon_{mA}+\varepsilon_{A}}^{2}}}$
(139)
$\displaystyle\pqty{\vec{F}_{A}-C_{11}\vec{1}}^{T}\partial_{t_{rA}}\vec{\xi}$
$\displaystyle=t_{rA}\pqty{C_{21}-C_{11}}\pqty{\frac{1}{\pqty{U-U^{\prime\prime}-2\varepsilon_{A}}^{2}}-\frac{1}{\pqty{U-2U^{\prime}+U^{\prime\prime}-\varepsilon_{mA}-\varepsilon_{A}}^{2}}}.$
(140)

When the terms in the larger brackets go to zero, we get the sweet spots for
the tunneling parameters. This occurs when the denominators are equal

The equations with $-\varepsilon_{A}$ and $+3\varepsilon_{A}$ terms are the
$t_{lA}$ two-qubit sweet spots ($t_{lA}$ SS), while the $+\varepsilon_{A}$ and
$-3\varepsilon_{A}$ terms are the $t_{rA}$ SS. The corresponding equations for
qubit B can be found by swapping $t_{lA}\to t_{rB}$, $t_{rA}\to t_{lB}$, and
$\varepsilon_{mA}\to-\varepsilon_{mB}$.

The sweet spot given by Eqs. (141b) & (142b) rarely appears (see Fig. 6),
while the condition for Eq. (141a) to exist is constrained by the requirement
for the qubit to remain in the (1,1,1) region

$\begin{gathered}-(U-U^{\prime})<\pqty{\varepsilon_{mA}+U^{\prime}-U^{\prime\prime}}\pm\varepsilon_{A}<U-U^{\prime}\\\
-(U-U^{\prime})<\pqty{\mp\varepsilon_{A}-2\pqty{U^{\prime}-U^{\prime\prime}}}+\pqty{U^{\prime}-U^{\prime\prime}}\pm\varepsilon_{A}<U-U^{\prime}\\\
\implies U^{\prime}-U^{\prime\prime}<U-U^{\prime},\end{gathered}$ (143)

where only one side of the inequality was needed as
$U>U^{\prime}>U^{\prime\prime}$.

It seems that we can find a double sweet spot for both tunneling noises by
equating both forms of Eq. (141a) to find
$(\varepsilon_{A},\varepsilon_{mA})=(0,-2(U^{\prime}-U^{\prime\prime}))$.
Unfortunately, the Makhlin invariants, as laid out in Sec. 3, are not
satisfied with these operating parameters. This causes the two-qubit gate
times to blow up to infinity as discussed in the main text, and hence cannot
be actually implemented.

#### Detuning Sweet Spots

To find the middle detuning sweet spot, we bring together Eqs. (134d) & (137),
and require that
${[{(\vec{F}_{A}-C_{11}\vec{1})}^{T}\partial_{\varepsilon_{mA}}\vec{\xi}]}_{\varepsilon_{mA}\text{SS}}=0$

$\begin{gathered}0=-\frac{t_{rA}^{2}\pqty{C_{21}-C_{11}}}{\pqty{U-2U^{\prime}+U^{\prime\prime}-\varepsilon_{mA}-\varepsilon_{A}}^{3}}+\frac{t_{lA}^{2}\pqty{C_{51}-C_{11}}}{\pqty{U-2U^{\prime}+U^{\prime\prime}-\varepsilon_{mA}+\varepsilon_{A}}^{3}}\\\
\therefore\frac{t_{lA}}{t_{rA}}=\sqrt{\frac{C_{51}-C_{11}}{C_{21}-C_{11}}}\pqty{\frac{U-2U^{\prime}+U^{\prime\prime}-\varepsilon_{mA}-\varepsilon_{A}}{U-2U^{\prime}+U^{\prime\prime}-\varepsilon_{mA}+\varepsilon_{A}}}^{\frac{3}{2}}.\end{gathered}$
(144)

Hence, for the $\varepsilon_{mA}$ sweet spot, we need to adjust the ratios of
the tunneling parameters to the one given in Eq. (144). This ratio has to be
positive for the sweet spot to exist, and the term in $(\cdots)^{\frac{3}{2}}$
is always positive in the (1,1,1) region, so the requirement is

$\operatorname{sgn}\pqty{C_{51}-C_{11}}=\operatorname{sgn}\pqty{C_{21}-C_{11}},$
(145)

which is purely determined by the geometry of the setup. In the linear array,
as we have used for our numerical investigations,
$C_{51}-C_{11}=\sum_{j}\mathcal{V}_{2j}-\mathcal{V}_{1j}$ and
$C_{21}-C_{11}=\sum_{j}\mathcal{V}_{3j}-\mathcal{V}_{2j}$ are both positive,
since the first (second) dot is farther from qubit B than the second (third)
dot. However, this might not be true in a different geometry, in which case an
$\varepsilon_{mA}$ sweet spot will not exist.

For $\varepsilon_{mB}$, we find

$\frac{t_{lB}}{t_{rB}}=\sqrt{\frac{C_{13}-C_{11}}{C_{14}-C_{11}}}\pqty{\frac{U-2U^{\prime}+U^{\prime\prime}-\varepsilon_{mB}-\varepsilon_{B}}{U-2U^{\prime}+U^{\prime\prime}-\varepsilon_{mB}+\varepsilon_{B}}}^{\frac{3}{2}}$
(146)

with the condition

$\operatorname{sgn}\pqty{C_{13}-C_{11}}=\operatorname{sgn}\pqty{C_{14}-C_{11}},$
(147)

for $C_{13}-C_{11}=\sum_{j}\mathcal{V}_{j4}-\mathcal{V}_{j5}$ and
$C_{14}-C_{11}=\sum_{j}\mathcal{V}_{j5}-\mathcal{V}_{j6}$, which are again
positive for the linear array because the fifth (sixth) dot is farther from
qubit A than the fourth (fifth) dot.

The strategy for finding the $\varepsilon_{A/B}$ sweet spots are identical,
and we find

$\displaystyle\frac{t_{lA}}{t_{rA}}$
$\displaystyle=\sqrt{\pqty{-\frac{C_{51}-C_{11}}{C_{21}-C_{11}}}\pqty{\frac{\frac{2}{\pqty{U-U^{\prime\prime}+2\varepsilon_{A}}^{3}}-\frac{1}{\pqty{U-2U^{\prime}+U^{\prime\prime}-\varepsilon_{mA}+\varepsilon_{A}}^{3}}}{\frac{2}{\pqty{U-U^{\prime\prime}-2\varepsilon_{A}}^{3}}-\frac{1}{\pqty{U-2U^{\prime}+U^{\prime\prime}-\varepsilon_{mA}-\varepsilon_{A}}^{3}}}}}$
(148a) $\displaystyle\frac{t_{lB}}{t_{rB}}$
$\displaystyle=\sqrt{\pqty{-\frac{C_{13}-C_{11}}{C_{14}-C_{11}}}\pqty{\frac{\frac{2}{\pqty{U-U^{\prime\prime}+2\varepsilon_{B}}^{3}}-\frac{1}{\pqty{U-2U^{\prime}+U^{\prime\prime}-\varepsilon_{mB}+\varepsilon_{B}}^{3}}}{\frac{2}{\pqty{U-U^{\prime\prime}-2\varepsilon_{B}}^{3}}-\frac{1}{\pqty{U-2U^{\prime}+U^{\prime\prime}-\varepsilon_{mB}-\varepsilon_{B}}^{3}}}}}.$
(148b)

Unlike the previous scenario, there are no simple conditions for the existence
of the $\varepsilon_{A/B}$ sweet spots apart from ensuring that the terms
inside the square root in Eqs. (148) are not negative. For the linear array
chosen for our paper, no $\varepsilon_{A/B}$ sweet spots exist.

###  7 Parameters of the Hubbard Model

In this study, we used the Hubbard model to describe both intra-TQD and inter-
TQD interactions. Intra-TQD interactions comprise QD detunings, tunnel
couplings as well as intra- and inter-dot Coulomb energies. Inter-TQD
interaction comprise inter-dot Coulomb energies only, when tunnel coupling
between the TQDs (QD 3 and 4) is zero.

The quadratic QD confinement potential was used to calculate inter-TQD Coulomb
interactions. This approach enables the analytical derivation of the
functional dependence of capacitive coupling on a minimal set of QD
geometrical parameters, Bohr radius $a_{B}$ and interdot distance $a$.
Importantly, these dependencies helped to identify the existence and
conditions of 2QSS on simple experimental parameters. A key feature of the
quadratic model is that simple analytical formulae connect geometrical
parameters with experimental parameters of tunnel coupling and Coulomb
energies, quantities utilized in the Hubbard model. This points to the
compatibility between both models.

Ideally, intra-qubit tunnel couplings and Coulomb energies should also be
derived with the same quadratic model, for consistency. Unfortunately, the
quadratic model does not correspond perfectly with experiment – tunnel
couplings are tuned in experiment by application of gate voltages that
modulate the barrier potential whereas in the quadratic model, tunnel coupling
can only be varied by changing interdot distance. In addition, changing
interdot distances modifies Coulomb energies in the quadratic model, which is
undesirable. Even though terms corresponding to tunnel barriers could be
added, or alternative Gaussian confinement models used, the drawback is that
simple analytical formula relating the inputs of these models with
experimental parameters do not exist. Because the goal was to utilize
theoretical models that would yield desired insights into 2QSS, we chose to
use the Hubbard model for its compatibility with the quadratic confinement
model and its excellent theory-experiment agreement [1, 2]:, in order to
describe intra-qubit interactions. Parameters for the Hubbard model were
estimated with reference to experiment [19] and theory [2] for silicon QDs,
and checked against those derived from the quadratic model to be within the
same order of magnitude.

## 0emupplementary Discussion

###  1 Effect of Lower Cutoff Frequency

For $1/f$ noise to have finite power, the spectrum has to be limited within a
higher and lower cutoff frequency. Increasing the lower cutoff shifts the
frequency range of the noise spectrum to the higher frequencies, and as the
power spectral density is inversely proportional to the frequency, this should
result in an improved fidelity, as can be seen in Fig. 2. The effect on
fidelity arises from the term $\varsigma(t)$ in the exponent of the fidelity
formula (Eq. (133)),

$\displaystyle\varsigma(t)=$
$\displaystyle\frac{1/2}{1+\ln(\omega_{h}/\omega_{l})}\bigg{(}\operatorname{Ci}(\omega_{h}t)t^{2}-\operatorname{Ci}(\omega_{l}t)t^{2}-\frac{\sin(\omega_{h}t)t}{\omega_{h}}+\frac{\sin(\omega_{l}t)t}{\omega_{l}}+\frac{2\operatorname{Si}(\omega_{l}t)t}{\omega_{l}}+\frac{\cos(\omega_{h}t)-1}{\omega_{h}^{2}}+\frac{\cos(\omega_{l}t)-1}{\omega_{l}^{2}}\bigg{)}$
(149)

where $\text{Ci}(x)$ and $\text{Si}(x)$ are the cosine and sine integrals
respectively. When $\omega_{l}t_{0}\ll 1$, increasing lower cutoff frequency
$\omega_{l}$ reduces the magnitude of $\varsigma(t_{0})$ and increases the
fidelity.

Supplementary Figure 2:  Effect of lower cutoff frequency $\omega_{l}/2\pi$ on
infidelity $1-\left\langle\mathcal{F}\right\rangle$. The standard deviation of
noise in detunings was fixed at
$\sigma_{\varepsilon_{A}}=\sigma_{\varepsilon_{mA}}=10^{-4}~{}\textrm{meV}$
for AEON (blue) and RX (red), and the standard deviation of noise in tunnel
couplings was fixed at
$\sigma_{t_{lA}}=\sigma_{t_{rA}}=10^{-5}~{}\textrm{meV}$ for AEON (magenta)
and RX (yellow), with both qubits sitting at the $\varepsilon_{mA}$ 2QSS.
Circles are numerical simulations (Eq. (13), main text), averaged over 500
noise realizations. Lines are analytical calculations (Eq. (14), main text).
As expected, infidelity decreases as the lower cutoff frequency increases,
while keeping $\omega_{h}$ constant.

###  2 Noise in Individual Qubits

In the main text, we only reported the cases where noise is present in the
control qubit. This is as there is not much more to be learnt from the cases
where noise is present in the target qubit, or where noise is present in both
qubits.

From Eq. (133), we can write the total fidelity as a product of two fidelities
$\expectationvalue{\mathcal{F}}=\expectationvalue{\mathcal{F}}_{A}\times\expectationvalue{\mathcal{F}}_{B}$
where $\expectationvalue{\mathcal{F}}_{A}$
($\expectationvalue{\mathcal{F}}_{B}$) is the fidelity of the system when
there is only noise in qubit A (B). Explicitly,

$\displaystyle\expectationvalue{\mathcal{F}}_{A}$
$\displaystyle=\left.\expectationvalue{\mathcal{F}}\right\rvert_{\sigma_{\mathbf{n}_{B}}=0}$
(150)
$\displaystyle=\frac{1}{4}\pqty{1+\exp\Bqty{-\frac{1}{2}\sum_{i}\pqty{\pqty{\vec{F}_{A}-C_{11}\vec{1}}^{T}\partial_{i}\vec{\xi}}^{2}\sigma_{n_{iA}}^{2}\varsigma(t_{0})}}\pqty{1+\exp\Bqty{0}}$
$\displaystyle=\frac{1}{2}\pqty{1+\exp\Bqty{-\frac{1}{2}\sum_{i}\pqty{\pqty{\vec{F}_{A}-C_{11}\vec{1}}^{T}\partial_{i}\vec{\xi}}^{2}\sigma_{n_{iA}}^{2}\varsigma(t_{0})}},$

and $\expectationvalue{\mathcal{F}}_{B}$ is obtained by swapping all
appearances of $\xi\to\nu$ and $A\to B$.

If both qubits are operated symmetrically —
$\varepsilon_{A}=-\varepsilon_{B}$, $\varepsilon_{mA}=\varepsilon_{mB}$,
$t_{l/rA}=t_{r/lB}$ — the noise coefficients
${\vec{F}_{A}}^{T}\partial_{i}\vec{\xi}={\vec{F}_{B}}^{T}\partial_{i}\vec{\nu}$
are equal for both qubits. Since $\varsigma(t_{0})$ is the same for both
qubits, the behaviour of the individual fidelities for both the control and
target qubits will be affected by noise in the same manner, as evident in Fig.
3.

Supplementary Figure 3: Noise in control qubit versus noise in target qubit.
The introduction of noise in either qubit results in the same effect on the
fidelity (note that $t_{lA}\leftrightarrow t_{rB}$ mirrors
$t_{rA}\leftrightarrow t_{lB}$ )

In addition, the individual fidelities combine in a predictable way. Comparing
the case where noise is present in either qubit individually, against that
where noise is present in both qubits concurrently, we will find that the
multiplication of the two individual fidelities in the former case gives the
same fidelity in the latter case, as shown through our numerical simulations
in Fig. 4. The result that
$\expectationvalue{\mathcal{F}}_{A}\times\expectationvalue{\mathcal{F}}_{B}=\expectationvalue{\mathcal{F}}$
holds true even when the qubit operation is not symmetric.

Supplementary Figure 4: Comparison of noise individually affecting either
qubits, versus the concurrent presence of noise in both qubits. We see that
the individual fidelities multiply to give the concurrent case.

###  3 Direct vs Exchange Coulomb Energies

Supplementary Figure 5: Dependence of direct (blue) and exchange (red) Coulomb
interactions on the distance between qubits, $R$. The ratio of direct to
exchange Coulomb interactions is plotted as a dashed purple line. The
different markers are used to represent the values for these parameters at the
value of $R$ chosen ($R=160$ nm), with the circle and triangle markers used
for the direct and exchange Coulomb interactions respectively and the square
marker fused for the ratio between them. At the $R$ value of choice, we
achieve an exchange Coulomb interaction ($\simeq 6\times 10^{-5}$ meV) about
four orders smaller than the direct Coulomb interaction ($\simeq 1$ meV).

One of the main motivations for this study was the advantage of using only a
single two-qubit gate to achieve a CPHASE or CNOT gate. This is in stark
contrast to complicated gate sequences that are required with exchange
interaction [20, 21, 22]. A key assumptions in our two-qubit system is that
the exchange interaction is weak enough such that it can be neglected,
allowing us to describe the two-qubit interaction as a purely diagonal
Hamiltonian. To achieve this, we could simply move the qubits away from each
other and at large enough distances, the direct Coulomb interaction dominates
as compared to the exchange Coulomb interaction, as shown in Fig. 5. The
exchange energies are calculated similarly to the direct Coulomb energies,

$\displaystyle\mathcal{X}_{ij}$
$\displaystyle=\int\int\psi_{i}^{\star}(\vec{r}_{1})\psi_{i}^{\star}(\vec{r}_{2})\frac{\kappa}{|\vec{r}_{1}-\vec{r}_{2}|}\psi_{j}(\vec{r}_{1})\psi_{j}(\vec{r}_{2})d\vec{r}_{1}d\vec{r}_{2},$
(151)

where $\kappa=q^{2}/4\pi\epsilon_{0}\epsilon_{r}$, as before. In our system,
we considered the quantum dot radius, $a_{B}=25$ nm, the inter-dot distance,
$2a=100$ nm, and the inter-qubit distance, $2R=320$ nm as measured from the
center dot of each qubit. This ensures dots 3 and 4 are well separated with a
distance of $120$ nm between them. At this distance, we compare the magnitudes
of the exchange and direct Coulomb interaction between the two closest dots in
the two-qubit system, namely dots 3 and 4 as labeled in Fig. 1(a) in the main
text. The interaction between these two dots are expected to dominate as
compared to the other inter-qubit interactions since they are the two closest
dots between the two qubits. As evidenced in Fig. 5, the exchange energy
($\mathcal{X}_{34}$ in red) is about four orders of magnitude smaller than the
corresponding direct Coulomb energy ($\mathcal{V}_{34}$ in blue) at this
distance. As seen in Fig. 6 the trade-off is longer gate times as inter-qubit
distance is increased.

Supplementary Figure 6: Plot of gate time $t_{0}$ with distance between qubits
$R$, at the single qubit sweet spot ($\varepsilon_{A/B}=0$). The shortest
possible gate times correspond to $k=1$. When $k=7$ for the RX qubit, its gate
time (448 ns) is comparable with the fastest AEON qubit gate time (450 ns). At
the value of $R$ chosen for the simulations ($R=160$nm) as shown by the
markers, we observe that $t_{0}$ remains reasonable at less than 0.5 $\mu$s.

## Supplementary References

  * [1] Yang, S., Wang, X. & Das Sarma, S. Generic Hubbard model description of semiconductor quantum-dot spin qubits. Phys. Rev. B 83, 161301 (2011). https://doi.org/10.1103/PhysRevB.83.161301
  * [2] Das Sarma, S., Wang, X. & Yang, S. Hubbard model description of silicon spin qubits: Charge stability diagram and tunnel coupling in Si double quantum dots. Phys. Rev. B 83, 235314 (2011). https://doi.org/10.1103/PhysRevB.83.235314
  * [3] Russ, M. & Burkard, G. Three-electron spin qubits. J. Phys. Condens. Matter 29, 393001 (2017). https://doi.org/10.1088/1361-648X/aa761f
  * [4] Zhang, C., Yang, XC. & Wang, X. Leakage and sweet spots in triple-quantum-dot spin qubits: A molecular-orbital study. Phys. Rev. A 97, 042326 (2018). https://doi.org/10.1103/PhysRevA.97.042326
  * [5] Annavarapu, R. N. Singular Value Decomposition and the Centrality of Löwdin Orthogonalizations. Am. J. Comput. Appl. Math. 3, 33–35 (2013). https://doi.org/10.5923/j.ajcam.20130301.06
  * [6] Calderon-Vargas, F. A. & Kestner, J. P. Directly accessible entangling gates for capacitively coupled singlet-triplet qubits. Phys. Rev. B 91, 035301 (2015). https://doi.org/10.1103/PhysRevB.91.035301
  * [7] Patzelt, F. Python package to generate gaussian (1/f)**beta noise. GitHub (2017). Available at: https://github.com/felixpatzelt/colorednoise. (Accessed: 15th August 2019)
  * [8] Timmer, J. & Koenig, M. On generating power law noise. Astron. Astrophys. 300, 707–710 (1995).
  * [9] Watson, T., Philips, S., Kawakami, E. et al. A programmable two-qubit quantum processor in silicon. Nature 555, 633–637 (2018). https://doi.org/10.1038/nature25766
  * [10] Frees, A., Mehl, S., Gamble, J.K. et al. Adiabatic two-qubit gates in capacitively coupled quantum dot hybrid qubits. npj Quantum Inf. 5, 73 (2019). https://doi.org/10.1038/s41534-019-0190-7
  * [11] Chatfield, C. The Analysis of Time Series. (Chapman and Hall/CRC, 2003).
  * [12] Makhlin, Y. Nonlocal Properties of Two-Qubit Gates and Mixed States, and the Optimization of Quantum Computations. Quantum Inf. Process. 1, 243–252 (2002). https://doi.org/10.1023/A:1022144002391
  * [13] Green, T. J., Sastrawan, J., Uys, H. & Biercuk, M. J. Arbitrary quantum control of qubits in the presence of universal noise. New J. Phys. 15, 095004 (2013). https://doi.org/10.1088/1367-2630/15/9/095004
  * [14] Yang, YC., Coppersmith, S.N. & Friesen, M. Achieving high-fidelity single-qubit gates in a strongly driven charge qubit with 1/f charge noise. npj Quantum Inf. 5, 12 (2019). https://doi.org/10.1038/s41534-019-0127-1
  * [15] James, D. F., Kwiat, P. G., Munro, W. J. & White, A. G. Measurement of qubits. Phys. Rev. A 64, 052312 (2001). https://doi.org/10.1103/PhysRevA.64.052312
  * [16] O’Brien, J. L. et al. Quantum Process Tomography of a Controlled-NOT Gate. Phys. Rev. Lett. 93, 080502 (2004). https://doi.org/10.1103/PhysRevLett.93.080502
  * [17] Yamamoto, T. et al. Quantum process tomography of two-qubit controlled-Z and controlled-NOT gates using superconducting phase qubits. Phys. Rev. B 82, 184515 (2010). https://doi.org/10.1103/PhysRevB.82.184515
  * [18] Kubo, R. Generalized Cumulant Expansion Method. J. Phys. Soc. Jpn. 17, 1100–1120 (1962). https://doi.org/10.1143/JPSJ.17.1100
  * [19] Neyens, S. F. et al. Measurements of Capacitive Coupling Within a Quadruple-Quantum-Dot Array. Phys. Rev. Appl. 12, 064049 (2019). https://doi.org/10.1103/PhysRevApplied.12.064049
  * [20] Fong, B. H. & Wandzura, S. M. Universal quantum computation and leakage reduction in the 3-qubit decoherence free subsystem. Quantum Inf. Comput. 11, 1003–1018 (2011). https://doi.org/10.5555/2230956.2230965
  * [21] Shi, Z. et al. Fast Hybrid Silicon Double-Quantum-Dot Qubit. Phys. Rev. Lett. 108, 140503 (2012). https://doi.org/10.1103/PhysRevLett.108.140503
  * [22] DiVincenzo, D., Bacon, D., Kempe, J. et al. Universal quantum computation with the exchange interaction. Nature 408, 339–342 (2000). https://doi.org/10.1038/35042541


# The CubeSpec space mission: Asteroseismology of massive stars from time-
series optical spectroscopy

D. M. Bowman1 B. Vandenbussche1 H. Sana1 A. Tkachenko1
G. Raskin1 T. Delabie2 B. Vandoren2 P. Royer1 S. Garcia1
T. Van Reeth1 and the CubeSpec Collaboration  1Institute of Astronomy, KU
Leuven, Celestijnenlaan 200D, 3001 Leuven, Belgium
email<EMAIL_ADDRESS>
2Arcsec NV, Blijde Inkomststraat 22, 3000 Leuven, Belgium

(2022)

###### Abstract

The ESA/KU Leuven CubeSpec mission is specifically designed to provide low-
cost space-based high-resolution optical spectroscopy. Here we highlight the
science requirements and capabilities of CubeSpec. The primary science goal is
to perform pulsation mode identification from spectroscopic line profile
variability and empower asteroseismology of massive stars.

###### keywords:

instrumentation: spectrographs, stars: early-type, stars: oscillations

††volume: 361††journal: Massive Stars Near and Far††editors: N. St-Louis, J.
S. Vink & J. Mackey, eds.

##  1 Scientific Context

There is a need to constrain the interior rotation, mixing and angular
momentum transport mechanisms of massive stars, as they strongly influence
stellar evolution. Asteroseismology – the study of stellar structure from
pulsations – is an excellent method for probing the interiors of massive stars
[Bowman (2020), (Bowman 2020)]. The ongoing NASA TESS mission is providing
high-precision light curves of thousands of massive stars [Ricker et al.
(2015), (Ricker et al. 2015)]. Complementary time-series spectroscopy is
highly advantageous in fully characterising pulsations: spectral line profile
variability reveals the geometry of pulsation modes in terms of spherical
harmonics [Aerts et al. (1994), Briquet2003, (e.g. Aerts et al. 1994; Briquet
et al. 2009)]. The combination of TESS photometry and CubeSpec spectroscopy
will unlock the physics of massive star interiors through forward
asteroseismic modelling.

##  2 CubeSpec Specifications

CubeSpec is a 6U cubesat on schedule to be launched in 2024. The platform
includes a Cassegrain telescope with a rectangular primary mirror of $9\times
19$ cm2, and a compact high-resolution echelle optical spectrograph with a
resolving power of $R=55000$ [Raskin et al. (2018), (Raskin et al. 2018)].
Schematics of the platform and optical design are shown in Fig. 1, and on the
CubeSpec project website: https://fys.kuleuven.be/ster/instruments/cubespec.

To achieve the primary science goal to assemble time-series spectroscopy for
massive-star asteroseismology, CubeSpec has the following science
requirements: (i) High resolving power ($R>50000$) to resolve perturbations to
narrow spectral lines in slowly rotating stars; (ii) High signal-to-noise
ratio (SNR $>200$) spectra to detect small perturbations relative to the
median line profile; (iii) Observing cadence of order 90 min to effectively
sample the pulsation periods of order several hours; (iv) Short exposure times
of order minutes to avoid pulsational smearing of spectral lines; (v) Mission
duration exceeding three months to provide a high frequency precision; and
(vi) Optical wavelength range to include the silicon triplet at
$\lambda\lambda$ 4552, 4567, and 4574 Å, which is sensitive to pulsations.

Figure 1: Platform layout (left) and optical design (right) of the CubeSpec
space mission.

##  3 Prioritised Target List

The compact design and science requirements of CubeSpec mean that bright (i.e.
$V<4$ mag) pulsating early-type stars are optimal targets to achieve its
asteroseismic science case (see Fig. 2). CubeSpec targets $\beta$ Cep stars,
which have pulsation periods between $2<P<8$ hr and masses in the approximate
range of $6<M<30$ M⊙ [Bowman (2020), (Bowman 2020)]. A search for pulsating
stars with $V<4$ mag and known spectral types between O8 and B3 (including
dwarfs, giants and supergiants) was performed a priori using previous
literature studies of photometric and spectroscopic variability and new TESS
light curves [Bowman et al. (2022), (Bowman et al. 2022)]. In total, 23
$\beta$ Cep stars were identified from a sample of 90 stars, with the highest
priority targets being slow rotators with large amplitude and long period
pulsations. Their sky location is shown in Fig. 2. Therefore, ample pulsating
massive stars exist to be targeted by CubeSpec and empower asteroseismology of
massive stars.

Figure 2: Left: SNR curves for simulated B0 V stars of different $V$ mag
brightnesses for conservative instrument performance allocation and best case
performance scenario using exposure times of 120 and 900 s. The SNR is
estimated at a wavelength of 5000 Å assuming $R=55000$. Right: Sky
distribution of high-priority CubeSpec targets from [Bowman et al. (2022),
Bowman et al. (2022)].

## References

  * [Aerts et al. (1994)] Aerts et al. 1994, A&A, 286, 109
  * [Briquet et al. (2009)] Briquet et al. 2009, A&A, 506, 269–280
  * [Bowman et al. (2022)] Bowman et al. 2022, A&A, 658, A96
  * [Bowman (2020)] Bowman 2020, Frontiers Astron. Space Sci., 7, 70
  * [Raskin et al. (2018)] Raskin et al. 2018, SPIE proceedings, 10698, id. 106985R
  * [Ricker et al. (2015)] Ricker at al. 2015, Astron. Tel. Inst. Sys., 1, 014001


* Dusling et al. (2016) Dusling, K.; Li, W.; Schenke, B.  Novel collective phenomena in high-energy proton–proton and proton–nucleus collisions.  Int. J. Mod. Phys. E 2016, 25, 1630002.  https://doi.org/10.1142/S0218301316300022.
  * Zhang et al. (2019) Zhang, C.; Marquet, C.; Qin, G.Y.; Wei, S.Y.; Xiao, B.W.  Elliptic Flow of Heavy Quarkonia in $pA$ Collisions.  Phys. Rev. Lett. 2019, 122, 172302.  https://doi.org/10.1103/PhysRevLett.122.172302.
  * Zhang et al. (2020) Zhang, C.; Marquet, C.; Qin, G.Y.; Shi, Y.; Wang, L.; Wei, S.Y.; Xiao, B.W.  Collectivity of heavy mesons in proton-nucleus collisions.  Phys. Rev. D 2020, 102, 034010.  https://doi.org/10.1103/PhysRevD.102.034010.
  * Adam et al. (2016a) Adam, J.; et al. (ALICE Collaboration)  Elliptic flow of muons from heavy-flavour hadron decays at forward rapidity in Pb–Pb collisions at $\sqrt{s_{\rm NN}}=2.76$ TeV.  Phys. Lett. 2016, B753, 41–56.  https://doi.org/10.1016/j.physletb.2015.11.059.
  * Adam et al. (2016b) Adam, J.; et al. (ALICE Collaboration)  Elliptic flow of electrons from heavy-flavour hadron decays at mid-rapidity in Pb-Pb collisions at $\sqrt{{\mathrm{s}}_{\mathrm{NN}}}=2.76$ TeV.  J. High Energy Phys. 2016, 09, 28.  https://doi.org/10.1007/JHEP09(2016)028.
  * Acharya et al. (2017) Acharya, S.; et al. (ALICE Collaboration)  D-meson azimuthal anisotropy in mid-central Pb-Pb collisions at $\mathbf{\sqrt{s_{\rm NN}}=5.02}$ TeV. _arXiv_ 2017, arXiv:1707.01005.
  * Abelev et al. (2014) Abelev, B.B.; et al. (ALICE Collaboration)  Azimuthal anisotropy of D meson production in Pb-Pb collisions at $\sqrt{s_{\rm NN}}=2.76$ TeV.  Phys. Rev. C 2014, 90, 034904.  https://doi.org/10.1103/PhysRevC.90.034904.
  * Acharya et al. (2017) Acharya, S.; et al. (ALICE Collaboration)  J/$\psi$ elliptic flow in Pb-Pb collisions at $\sqrt{s_{\mathrm{NN}}}=5.02$ TeV.  Phys. Rev. Lett. 2017, 119, 242301.  https://doi.org/10.1103/PhysRevLett.119.242301.
  * CMS (2022) Measurements of azimuthal anisotropy of nonprompt D0 mesons in PbPb collisions at $\sqrt{s_{\mathrm{NN}}}$ = 5.02 TeV. _arXiv_ 2022, arXiv:2212.01636.
  * Acharya et al. (2023) Acharya, S.; et al. (ALICE Collaboration)  Measurement of Non-prompt $\rm D^{0}$-meson Elliptic Flow in Pb-Pb Collisions at $\sqrt{s_{\rm NN}}$ = 5.02 TeV. _arXiv_ 2023, arXiv:2307.14084.
  * Li et al. (2019) Li, H.; Lin, Z.W.; Wang, F.  Charm quarks are more hydrodynamic than light quarks in final-state elliptic flow.  Phys. Rev. C 2019, 99, 044911.  https://doi.org/10.1103/PhysRevC.99.044911.
  * Acharya et al. (2018) Acharya, S.; et al. (ALICE Collaboration)  Search for collectivity with azimuthal J/$\psi$-hadron correlations in high multiplicity p-Pb collisions at $\sqrt{s_{\rm NN}}$ = 5.02 and 8.16 TeV.  Phys. Lett. B 2018, 780, 7–20.  https://doi.org/10.1016/j.physletb.2018.02.039.
  * Sirunyan et al. (2019) Sirunyan, A.M.; et al. (CMS Collaboration)  Observation of prompt J/$\psi$ meson elliptic flow in high-multiplicity pPb collisions at $\sqrt{s_{\mathrm{NN}}}=$ 8.16 TeV.  Phys. Lett. B 2019, 791, 172–194.  https://doi.org/10.1016/j.physletb.2019.02.018.
  * Acharya et al. (2023) Acharya, S.; et al. (ALICE Collaboration)  Measurements of azimuthal anisotropies at forward and backward rapidity with muons in high-multiplicity p–Pb collisions at $\sqrt{s_{\mathrm{NN}}}=$8.16 TeV.  Phys. Lett. B 2023, 846, 137782.  https://doi.org/10.1016/j.physletb.2023.137782.
  * Aad et al. (2020) Aad, G.; et al. (ATLAS Collaboration)  Measurement of azimuthal anisotropy of muons from charm and bottom hadrons in Pb+Pb collisions at $\sqrt{s_{\mathrm{NN}}}=$ 5.02 TeV with the ATLAS detector.  Phys. Lett. B 2020, 807, 135595.  https://doi.org/10.1016/j.physletb.2020.135595.
  * Nahrgang et al. (2014) Nahrgang, M.; Aichelin, J.; Gossiaux, P.B.; Werner, K.  Influence of hadronic bound states above $T_{c}$ on heavy-quark observables in Pb + Pb collisions at at the CERN Large Hadron Collider.  Phys. Rev. C 2014, 89, 014905.  https://doi.org/10.1103/PhysRevC.89.014905.
  * Ke et al. (2019) Ke, W.; Xu, Y.; Bass, S.A.  Modified Boltzmann approach for modeling the splitting vertices induced by the hot QCD medium in the deep Landau-Pomeranchuk-Migdal region.  Phys. Rev. C 2019, 100, 064911.  https://doi.org/10.1103/PhysRevC.100.064911.
  * Katz et al. (2020) Katz, R.; Prado, C.A.G.; Noronha-Hostler, J.; Noronha, J.; Suaide, A.A.P.  Sensitivity study with a D and B mesons modular simulation code of heavy flavor RAA and azimuthal anisotropies based on beam energy, initial conditions, hadronization, and suppression mechanisms.  Phys. Rev. C 2020, 102, 024906.  https://doi.org/10.1103/PhysRevC.102.024906.
  * Lin et al. (2005) Lin, Z.W.; Ko, C.M.; Li, B.A.; Zhang, B.; Pal, S.  A Multi-phase transport model for relativistic heavy ion collisions.  Phys. Rev. C 2005, 72, 064901.  https://doi.org/10.1103/PhysRevC.72.064901.
  * Lin and Zheng (2021) Lin, Z.W.; Zheng, L.  Further developments of a multi-phase transport model for relativistic nuclear collisions.  Nucl. Sci. Tech. 2021, 32, 113.  https://doi.org/10.1007/s41365-021-00944-5.
  * Zhang (1998) Zhang, B.  ZPC 1.0.1: A Parton cascade for ultrarelativistic heavy ion collisions.  Comput. Phys. Commun. 1998, 109, 193–206.  https://doi.org/10.1016/S0010-4655(98)00010-1.
  * He et al. (2016) He, L.; Edmonds, T.; Lin, Z.W.; Liu, F.; Molnar, D.; Wang, F.  Anisotropic parton escape is the dominant source of azimuthal anisotropy in transport models.  Phys. Lett. B 2016, 753, 506–510.  https://doi.org/10.1016/j.physletb.2015.12.051.
  * Du and Rapp (2015) Du, X.; Rapp, R.  Sequential Regeneration of Charmonia in Heavy-Ion Collisions.  Nucl. Phys. A 2015, 943, 147–158.  https://doi.org/10.1016/j.nuclphysa.2015.09.006.
  * Khachatryan et al. (2017) Khachatryan, V.; et al. (CMS Collaboration)  Evidence for collectivity in pp collisions at the LHC.  Phys. Lett. B 2017, 765, 193–220.  https://doi.org/10.1016/j.physletb.2016.12.009.