\section{Introduction}
This paper presents the ESPnet-ST group's English$\to$German speech translation (ST) system submitted to the IWSLT 2021 offline speech translation track.
ESPnet~\citep{watanabe2018espnet} has been widely used for many speech applications; automatic speech recognition (ASR), text-to-speech~\citep{hayashi2020espnet}, speech translation~\citep{inaguma-etal-2020-espnet}, machine translation (MT), and speech separation/enhancement~\citep{li2020espnet}.
The purpose of this submission is not only to show the recent progress on ST researches, but also to encourage future research by building strong systems along with the open-sourced project.

This year we focused on (1) sequence-level knowledge distillation (SeqKD)~\citep{kim-rush-2016-sequence}, (2) Conformer encoder~\citep{gulati2020}, (3) Multi-Decoder architecture~\citep{dalmia2021searchable}, (4) model ensembling, and (5) better segmentation with a neural network-based voice activity (VAD) system~\citep{bredin2020pyannote} and a novel algorithm to merge multiple short segments for long context modeling.
Our primary focus was E2E models, although we also compared them with cascade systems with our best effort.
All experiments were conducted with the ESPnet-ST toolkit~\citep{inaguma-etal-2020-espnet}, and the recipe is publicly available at \url{https://github.com/espnet/espnet/tree/master/egs/iwslt21}.

\section{Data preparation}
In this section, we describe data preparation for each task.
The corpus statistics are listed in Table~\ref{tab:corpus}.
We removed the off-limit talks following previous evaluation campaigns\footnote{\url{https://sites.google.com/view/iwslt-evaluation-2019/speech-translation/off-limit-ted-talks}}.
To fit the GPU memory, we excluded utterances having more than 3000 speech frames or more than 400 characters.
All sentences were tokenized with the {\tt tokenizer.perl} script in the Moses toolkit~\citep{koehn-etal-2007-moses}.

\begin{table}[t]
    \centering
    \small
    \begingroup
    \scalebox{0.95}{
    \begin{tabular}{lcc}\toprule
     & \textbf{\#Hour} & \textbf{\#Sentence}  \\ \midrule
     \textbf{ASR} & &  \\
      \ Must-C & 408 $\times$ 3 & \phantom{0}0.68M \\
      \ Must-C v2 & 458 $\times$ 3 & \phantom{0}0.74M \\
      \ ST-TED (cleaned) & 200 $\times$ 3 & \phantom{0}0.40M \\
      \ Librispeech & 960 & \phantom{0}0.28M \\
      \ TEDLIUM2 & 210 $\times$ 3 & \phantom{0}0.27M \\
      \midrule

     \textbf{E2E-ST} & &  \\
      \ Must-C & 408 $\times$ 3 & \phantom{0}0.68M \\
      \ Must-C v2 & 458 $\times$ 3 & \phantom{0}0.74M \\
      \ ST-TED (cleaned) & 200 $\times$ 3 & \phantom{0}0.40M \\
      \midrule

     \textbf{MT} & &  \\
      \ Must-C &  \multirow{9}{*}{-} & \phantom{0}0.68M \\
      \ Must-C v2 &  & \phantom{0}0.74M \\
      \ ST-TED (cleaned) &  & \phantom{0}0.40M \\
      \ Europarl &  & \phantom{0}1.82M \\
      \ Commoncrawl &  & \phantom{0}2.39M \\
      \ Paracrawl &  & 34.37M \\
      \ NewsCommentary &  & \phantom{0}0.37M \\
      \ WikiTitles &  & \phantom{0}1.38M \\
      \ RAPID &  & \phantom{0}1.63M \\
      \ WikiMatrix &  & \phantom{0}1.57M \\
      \bottomrule
    \end{tabular}
    }
    \endgroup
    \caption{Corpus statistics}
    \label{tab:corpus}
\end{table}

\subsection{ASR}
We used Must-C~\citep{di-gangi-etal-2019-must}, Must-C v2\footnote{\url{https://ict.fbk.eu/must-c-release-v2-0/}}, ST-TED~\citep{jan2018iwslt}, Librispeech~\citep{librispeech}, and TEDLIUM2~\citep{tedlium} corpora.
We used the cleaned version of ST-TED following~\citep{inaguma19asru}.
The speech data was augmented by three-fold speed perturbation~\citep{speed_perturbation} with speed ratios of 0.9, 1.0, and 1.1 except for Librispeech.
We removed case information and punctuation marks except for apostrophes from the transcripts.
The 5k unit vocabulary was constructed based on the byte pair encoding (BPE) algorithm~\citep{sennrich-etal-2016-neural} with the \texttt{sentencepiece} toolkit\footnote{\url{https://github.com/google/sentencepiece}} using the English transcripts only.

\subsection{E2E-ST}
We used Must-C, Must-C v2, and ST-TED only.
The shared source and target vocabulary of BPE16k units was constructed using cased and punctuated transcripts and translations.

\subsection{MT}
We used available bitext for WMT20\footnote{Europarl, Commoncrawl, Paracrawl, NewsCommentary, WikiTitles, RAPID and WikiMatrix} in addition to the in-domain TED data used for E2E-ST systems.
We first performed perplexity-based filtering with an in-domain n-gram language model (LM)~\citep{moore-lewis-2010-intelligent}.
We controlled the WMT data size by thresholding and obtained three data pools: 5M, 10M, and 20M sentences.
Next, we removed non-printing characters and performed language identification with the \texttt{langid.py} toolkit~\citep{lui-baldwin-2012-langid}\footnote{\url{https://github.com/saffsd/langid.py}} and kept sentences whose language IDs were identified correctly on both English and German sides.
We also removed sentences having more than 250 tokens in either language or a source-target length ratio of more than 1.5 with the \texttt{clean-corpus-n.perl} script in Moses.
Finally, we removed sentences having CJK and other unrelated characters in either language with the built-in \texttt{regex} module in Python.
The resulting data size is shown in Table~\ref{tab:mt_bitext_filtering}.
We found that our filtering strategy removed 22-37\% of data.
Note that the above filtering process was performed over the WMT data only.
For each data size, the joint source and target vocabulary of BPE32k units was constructed using cased and punctuated sentences after the filtering.
We did not use additional monolingual data.

\begin{table}[t]
    \centering
    \small
    \begingroup
    \scalebox{0.95}{
    \begin{tabular}{lccc}\toprule
      \multirow{2}{*}{Filtering method} & \multicolumn{3}{c}{\textbf{\#Sentence}} \\ \cmidrule{2-4}
      & \texttt{WMT5M} & \texttt{WMT10M} & \texttt{WMT20M} \\ \midrule
      In-domain LM & 5.00M & 10.00M & 20.00M \\
      \ + \texttt{langid} & 3.42M & \phantom{0}7.90M & 15.33M \\
      \ \ + length/character & 3.15M & \phantom{0}7.77M & 15.01M \\
      \bottomrule
    \end{tabular}
    }
    \endgroup
    \caption{MT bitext filtering}
    \label{tab:mt_bitext_filtering}
\end{table}

\begin{figure}[t]
    \centering
    \includegraphics[width=0.37\linewidth]{Conformer_block.pdf}
    \caption{Block diagram of Conformer architecture}
    \label{fig:conformer}
\end{figure}

\section{System}
\subsection{Conformer encoder}
Conformer encoder~\citep{gulati2020} is a stacked multi-block architecture and has shown consistent improvement over a wide range of E2E speech processing applications~\citep{guo2020recent}.
The architecture of each block is depicted in Figure~\ref{fig:conformer}.
It includes a multi-head self-attention module, a convolution module, and a pair of position-wise feed-forward modules in the Macaron-Net style.
While the self-attention module learns the long-range global context, the convolution module aims to model the local feature patterns synchronously.
Recent studies have shown improvements by introducing Conformer in the E2E-ST task~\citep{guo2020recent,inaguma2021source}, which motivated us to adopt this architecture as our system.

\subsection{SeqKD}
Sequence-level knowledge distillation (SeqKD)~\citep{kim-rush-2016-sequence} is an effective method to transfer knowledge in a teacher model to a student model via discrete symbols.
Our recent studies~\citep{inaguma2020orthros,inaguma2021source} showed a large improvement in ST performance with this technique.
Unlike the previous studies, however, we used more training data than bitext in ST training data to train teacher MT models.
We translated source transcripts in the ST training data by the teacher MT models with a beam width of 5 and then replaced the original ground-truth translation with the generated translation.
We used cased and punctuated transcripts as inputs to the MT teachers.
We also combined both the original and pseudo translations as data augmentation (\textit{multi-referenced training})~\citep{gordon2019explaining}.

\begin{figure}
    \centering
    \includegraphics[width=0.99\linewidth]{md-fig.pdf}
    \caption{The Multi-Decoder (MD) architecture decomposes the overall ST task with ASR and MT sub-nets while maintaining E2E differentiability.}
    \label{fig:md-fig}
\end{figure}

\subsection{Multi-Decoder architecture}
The Multi-Decoder is an E2E-ST model using Searchable Hidden Intermediates to decompose the overall ST task into ASR and MT sub-tasks~\citep{dalmia2021searchable}.
As shown in Figure~\ref{fig:md-fig}, the Multi-Decoder consists of two encoder-decoder models, an ASR sub-net and a subsequent MT sub-net, where the hidden representations of the ASR decoder are passed as inputs to the encoder of the MT sub-net.
During inference, the best ASR decoder hidden representations are retrieved using beam search decoding at this intermediate stage.

Since this framework decomposes the overall ST task, it brings several advantages of cascaded approaches into the E2E setting.
For instance, the Multi-Decoder allows for greater search capabilities and separation of speech and text encoding.
However, one trade-off is a greater risk of error propagation from the ASR sub-net to the downstream MT sub-net.
To alleviate this issue, we condition the decoder of the MT sub-net on the ASR encoder hidden representations in addition to the MT encoder hidden representations using multi-source cross-attention.
This improved variant of the architecture is called the Multi-Decoder with Speech Attention.

\subsection{Model ensembling}
\label{sec:ensembling}
We use posterior probability combination to ensemble models trained with different data and architectures.
During inference, we perform a posterior combination at each step of beam search decoding by first computing the softmax normalized posterior probabilities for each model in the ensemble and then taking the mean value.
In this ensembling approach, a single unified beam search operates over the combined posteriors of the models to find the most likely decoded sequence.

\subsection{Segmentation}
How to segment audio during inference significantly impacts ST performances~\citep{Gaido2020,Pham2020,potapczyk-przybysz-2020-srpols,gaido2021beyond}.
This is because the ST systems are usually trained with utterances segmented based on punctuation marks~\citep{di-gangi-etal-2019-must} while the audio segmentation by voice activity detection (VAD) at test time does not access such meta information.
Since VAD splits a long speech recording into chunks by silence regions, it would prevent models from extracting semantically coherent contextual information.
Therefore, it is very important to seek a better segmentation strategy in order to minimize this gap in training and test conditions and evaluate models correctly.
In fact, the last year's winner obtained huge improvements by using their own segmentation strategy.

Motivated by this fact, we investigated two VAD systems apart from the provided segmentation.
Specifically, we used WebRTC\footnote{\url{https://github.com/wiseman/py-webrtcvad}} and pyannote.audio~\citep{bredin2020pyannote}\footnote{\url{https://github.com/pyannote/pyannote-audio}} toolkits.
For WebRTC, we set the frame duration, padding duration, and aggressive mode to 10ms, 150ms, and 3, respectively.
For pyannote.audio, we used a publicly available model pre-trained on the DIHARD corpus~\citep{Ryant2019}.

However, we observed that VAD systems are more likely to generate short segments because they do not take contextual information into account.
Therefore, we propose a novel algorithm to merge multiple short segments into a single chunk to enable long context modeling by self-attention in both encoder and decoder modules.
The proposed algorithm is shown in Algorithm~\ref{algo:merge_segment}.
We first perform VAD and obtain multiple segments.
Then, we check the segments in a greedy way from left to right and merge adjacent segments if (1) the total utterance duration is below a threshold $M_{\rm dur}$ [10ms] and (2) the time interval of the two segments is below a threshold $M_{\rm int}$ [10ms].
This process continues until no segment is merged in an iteration.
Although recent studies proposed similar methods~\citep{potapczyk-przybysz-2020-srpols,gaido2021beyond}, our algorithm is a bottom-up approach while theirs are top-down.

\begin{algorithm}[t]
\caption{Merge short segments after VAD for long context modeling}
\footnotesize
\begin{algorithmic}[1]
\Function{MergeSegment}{${\bm x}, M_{\rm dur}, M_{\rm int}$}
    \State $Q \gets VAD({\bm x})$  {\color{blue} \Comment{$\{(s_{1},e_{1}),\cdots,(s_{M},e_{M})\}$}}
    \While {True}
        \State $N_{\rm merge} \gets 0$
        \State $Q_{\rm next} \gets \{\}$ {\color{blue} \Comment{Queue}}
        \State $S, T \gets s_{1}, e_{1}$ {\color{blue} \Comment{Start/End time}}
        \For {$(s_{m},e_{m}) \in Q$}
            \If {$e_{m} - S < M_{\rm dur}$ and $s_{m} - E < M_{\rm int}$}
                \State $N_{\rm merge} \gets N_{\rm merge} + 1$ {\color{blue} \Comment{Merge segments}}
            \Else
                \State $Q_{\rm next}.enqueue((S,E))$
                \State $S \gets s_{m}$  {\color{blue} \Comment{Reset}}
            \EndIf
            \State $E \gets e_{m}$
        \EndFor
        \State $Q \gets Q_{\rm next}$

        \If {$N_{\rm merge} = 0$}
            \State \textbf{break}
        \EndIf
    \EndWhile

    \State \Return $Q$
\EndFunction
\end{algorithmic}\label{algo:merge_segment}
\end{algorithm}

\section{Experimental setting}
In this section, we describe the experimental setting for each task.
The detailed configurations for each task are summarized in Table~\ref{tab:training_config}.

\subsection{Feature extraction}
We extracted 80-channel log-mel filterbank coefficients computed with 25-ms window size and shifted every 10-ms with 3-dimensional pitch features using the Kaldi toolkit~\citep{kaldi}.
The features were normalized by the mean and the standard deviation calculated on the entire training set.
We applied SpecAugment~\citep{specaugment} with mask parameters $(m_{T}, m_{F}, T, F)=(2, 2, 40, 30)$ and time-warping for both ASR and E2E-ST tasks.

\begin{table}[t]
    \centering
    \begingroup
    \scalebox{0.72}{
    \begin{tabular}{lcccc}\toprule
    \multicolumn{1}{c}{\multirow{2}{*}{Configuration}} & \multicolumn{1}{c}{\multirow{2}{*}{ASR}} & \multicolumn{2}{c}{E2E-ST} & \multirow{2}{*}{MT} \\ \cmidrule{3-4}
    \multicolumn{1}{c}{} & \multicolumn{1}{c}{} & non-MD & MD &  \\
    \midrule
     Warmup step & 25k & 25k & 25k & 8k \\
     Learning rate factor & 10.0 & 2.5 & 12.5 & 1.0 \\
     Batch size & 200 utt & 128 utt & 120 utt & 65k tok \\
     Epoch & 30 & 30 & 30 & 40 \\
     Validation metric & Accuracy & BLEU & BLEU & BLEU \\
     Model average & 5 & 5 & 5 & 5 \\
     Beam width & 10 & 4 & 16, 10 & 4 \\
     \bottomrule
    \end{tabular}
    }
    \endgroup
    \caption{Summary of training configuration}\label{tab:training_config}
\end{table}

\subsection{ASR}
We used both Transformer and Conformer architectures.
The encoder had two CNN blocks followed by 12 Transformer/Conformer blocks following~\citep{karita2019comparative,guo2020recent}.
Each CNN block consisted of a channel size of 256 and a kernel size of 3 with a stride of 2 $\times$ 2, which resulted in time reduction by a factor of 4.
Both architectures had six Transformer blocks in the decoder.
In both encoder and decoder blocks, the dimensions of the self-attention layer $d_{\rm model}$ and feed-forward network $d_{\rm ff}$ were set to 512 and 2048, respectively.
The number of attention heads $H$ was set to 8.
The kernel size of depthwise separable convolution in Conformer blocks was set to 31.
We optimized the model with the joint CTC/attention objective~\citep{hybrid_ctc_attention} with a CTC weight of 0.3.
We also used CTC scores during decoding but did not use any external LM for simplicity.
We adopted the best model configuration from the Librispeech ASR recipe in ESPnet.

\subsection{MT}
We used the Transformer-Base and -Big configurations in~\citep{vaswani2017attention}.

\subsection{E2E-ST}
We used the same Conformer architecture as ASR except for the vocabulary.
We initialized the encoder parameters with those of the Conformer ASR.
On the decoder side, we initialized parameters like BERT~\citep{devlin-etal-2019-bert}, where weight parameters were sampled from ${\mathcal N}(0, 0.02)$, biases were set to zero, and layer normalization parameters were set to $\beta=0$, $\gamma=1$.
This technique led to better translation performance and faster convergence.

\begin{table}[t]
    \centering
    \begingroup
    \scalebox{0.72}{
    \begin{tabular}{lccc}\toprule
    \multirow{3}{*}{Model} & \multicolumn{3}{c}{\bf{WER} ($\downarrow$)} \\
    \cmidrule{2-4}
       & Librispeech & TEDLIUM2 & Must-C \\
       & \texttt{test-other} & \texttt{test} & \texttt{tst-COMMON} \\
       \midrule
      Transformer & 9.4 & 6.4 & 7.0	 \\
      Conformer & \bf{7.1} & \bf{6.2} & \bf{5.6} \\
      \bottomrule
    \end{tabular}
    }
    \endgroup
    \caption{Word error rate (WER) of ASR systems}
    \label{tab:result_asr}
\end{table}

\begin{table}[t]
    \centering
    \begingroup
    \scalebox{0.63}{
    \begin{tabular}{lccccccc}\toprule
    \multirow{2}{*}{VAD} & \multirow{2}{*}{$M_{\rm dur}$} & \multirow{2}{*}{$M_{\rm int}$} & \multicolumn{5}{c}{\bf{WER} ($\downarrow$)} \\
    \cmidrule(lr){4-7} \cmidrule(lr){8-8}
      &  &  & tst2010 & tst2015 & tst2018 & tst2019 & Avg. \\ \midrule
      \multirow{4}{*}{Provided}
      & -- & -- & 18.2 & 32.1 & 23.5 & 20.8 & 23.65 \\
      & 1500 & 200 & 14.4 & 29.3 & 18.4 & 15.5 & 19.40 \\
      & 2000 & 200 & 12.7 & 27.7 & 16.4 & 11.5 & 17.08 \\
      & 2500 & 200 & 14.5 & 29.9 & 15.1 & 12.2 & 17.93 \\ \midrule

      \multirow{4}{*}{WebRTC} & -- & -- & 35.3 & 35.1 & 44.0 & 22.7 & 34.28 \\
      & 1500 & 200 & 19.4 & 26.7 & 27.7 & 13.8 & 21.90 \\
      & 2000 & 200 & 19.8 & 27.7 & 27.1 & 11.9 & 21.63  \\
      & 2500 & 200 & 22.9 & 29.5 & 27.1 & 11.6 & 22.78 \\  \midrule

      \multirow{7}{*}{pyannote} & -- & -- & \phantom{0}9.5 & 24.0 & 15.5 & \phantom{0}7.3 & 14.08 \\
      & 1500 & 200 & \phantom{0}8.0 & 23.0 & 12.4 & \phantom{0}7.3 & 12.68 \\
      & 1500 & 100 & \phantom{0}7.5 & 22.2 & 12.4 & \phantom{0}6.5 & 12.15 \\
      & 2000 & 200 & 10.3 & 22.5 & 12.2 & \phantom{0}6.5 & 12.88 \\
      & 2000 & 150 & \phantom{0}9.6 & 21.8 & 12.3 & \phantom{0}6.1 & 12.45 \\
      & 2000 & 100 & \phantom{0}8.1 & \bf{21.5} & \bf{12.0} & \phantom{0}5.8 & 11.90 \\ %
      & 2000 & \phantom{0}50 & \phantom{0}\bf{7.3} & 21.9 & 12.4 & \phantom{0}\bf{5.9} & \bf{11.88} \\
      \bottomrule
    \end{tabular}
    }
    \endgroup
    \caption{Impact of audio segmentation for ASR}
    \label{tab:result_asr_segmenation}
\end{table}

\begin{table*}[t]
    \centering
    \small
    \begingroup
    \scalebox{0.96}{
    \begin{tabular}{lcccccccc}\toprule
      \multirow{3}{*}{Model} & \multicolumn{8}{c}{\bf{BLEU} ($\uparrow$)} \\
     \cmidrule(lr){2-3} \cmidrule(lr){4-4} \cmidrule(lr){5-8} \cmidrule(lr){9-9}
       & \multicolumn{2}{c}{Must-C} & Must-C v2 & \multirow{2}{*}{tst2010} & \multirow{2}{*}{tst2015} & \multirow{2}{*}{tst2018} & \multirow{2}{*}{tst2019} & Must-C \\
       & \texttt{dev} & \texttt{tst-COMMON} & \texttt{tst-COMMON} &  & & & & \texttt{Train} \\
      \midrule
      Base (Must-C only) & -- & 30.02 & 29.86 & 27.28 & 24.92 & 21.13 & 20.37 & \\
      \midrule
      Base (\texttt{WMT5M}) & 31.31 & 34.13 & 33.85 & 31.61 & 32.44 & 28.30 & 28.28 & 45.68 \\
      \ + Big & 27.32 & 29.11 & 28.85 & 27.61 & 28.44 & 24.42 & 23.92 & -- \\
      Base (\texttt{WMT10M}) & \bf{33.28} & \bf{35.09} & 34.80 & \bf{33.58} & 33.26 & 29.24 & 28.87 & 38.31 \\
      \ + In-domain finetune & 30.67 & \bf{35.50} & \bf{35.30} & 30.79 & 31.43 & 25.35 & 26.10 & -- \\
      Base (\texttt{WMT20M}) & 33.15 & 35.06 & \bf{34.87} & 33.26 & \bf{33.56} & \bf{29.94} & \bf{29.08} & 33.60 \\

      \bottomrule
    \end{tabular}
    }
    \endgroup
    \caption{BLEU scores of text-based MT systems}
    \label{tab:result_mt}
\end{table*}

\section{Results}
\subsection{ASR}
\subsubsection{Architecture}
We compared Transformer and Conformer ASR architectures in Table~\ref{tab:result_asr}.
We observed that Conformer significantly outperformed Transformer.
Therefore, we use the Conformer encoder in the following experiments.

\subsubsection{Segmentation}\label{sssec:segmentation_asr}
Next, we investigated the VAD systems and the proposed segment merging algorithm for long context modeling in Table~\ref{tab:result_asr_segmenation}.
We used the same decoding hyperparameters tuned on Must-C.
We firstly observed that merging short segments was very effective probably because it alleviated frame classification errors in the VAD systems.
Among three audio segmentation methods, we confirmed that pyannote.audio significantly reduced the WER while WebRTC had negative impacts compared to the provided segmentation.
Specifically, we found that the \texttt{dihard} option in pyannote.audio worked very well while the rest options did not.
The optimal maximum duration $M_{\rm dur}$ was around 2000 frames (i.e., 20 seconds).
In the last experiments, we tuned the maximum interval $M_{\rm int}$ among \{50, 100, 150, 200\} and found 50 and 100 (i.e., 0.5 and 1 second) was best on average.
Compared to the provided segmentation, we obtained a 49.6\% improvement on average.

\begin{table*}[t]
    \centering
    \begingroup
    \scalebox{0.66}{
    \begin{tabular}{llccccccccc}\toprule
    \multirow{3}{*}{ID} & \multirow{3}{*}{Model} & \multicolumn{8}{c}{\bf{BLEU} ($\uparrow$)} \\
    \cmidrule(lr){3-5} \cmidrule(lr){6-6} \cmidrule(lr){7-10}
       & & \multicolumn{3}{c}{Must-C} & Must-C v2 & \multirow{2}{*}{tst2010} & \multirow{2}{*}{tst2015} & \multirow{2}{*}{tst2018} & \multirow{2}{*}{tst2019} \\
       & & \texttt{dev} & \texttt{tst-COMMON} & \texttt{tst-HE} & \texttt{tst-COMMON} \\
       \midrule

       \multirow{6}{*}{-} & Bidir SeqKD (E2E)~\citep{inaguma2021source} & 25.67 & 27.01 & 25.36 & -- & -- & -- & -- & -- \\
       & Multi-Decoder (E2E)~\citep{dalmia2021searchable} & -- & 26.4\phantom{0} & -- & -- & -- & -- & -- & -- \\
       & RWTH (Cascade)~\citep{bahar2021tight} & -- & 26.50 & 26.80 & -- & -- & 28.4\phantom{0} & -- & -- \\
       & KIT (E2E)~\citep{Pham2020} & -- & 30.60 & -- & -- & 24.27 & 21.82 & -- & -- \\
       & KIT (Cascade)~\citep{Pham2020} & -- & -- & -- & -- & 26.68 & 24.95 & -- & -- \\
       & SRPOL (E2E)~\citep{potapczyk-przybysz-2020-srpols} & -- & -- & -- & -- & 29.44 & 24.6\phantom{0} & -- & 23.96 \\
      \midrule

      \texttt{A1} & Baseline (\texttt{X}) & 25.14 & \bf{35.63} & 22.63 & \bf{36.07} & 21.40 & 18.18 & 16.69 & 17.39 \\
      \texttt{A2} & \ + SeqKD (\texttt{Y}) & 26.31 & 29.29 & 26.33 & 29.50 & 23.34 & 21.24 & 21.09 & 22.25 \\
      \texttt{A3} & \ + 2ref SeqKD (\texttt{X}+\texttt{Y}) & 26.50 & 30.59 & 26.21 & 30.92 & 23.00 & 22.18 & 20.38 & 21.59 \\
      \texttt{A4} & \ + 3ref SeqKD (\texttt{X}+\texttt{Y}+\texttt{Z}) & \bf{27.66} & 30.90 & \bf{27.44} & 31.07 & \bf{24.97} & \bf{22.66} & \bf{22.20} & \bf{23.41} \\ \midrule

      \texttt{B1} & MD + 2ref SeqKD & -- & 30.78 & -- & -- & -- & -- & -- & \textbf{23.78} \\ \midrule

      \texttt{C1} & Conformer ASR $\to$ Base MT (\texttt{WMT10M}) & 27.01 & 29.42 & 26.13 & 29.75 & \bf{25.04} & \bf{23.17} & \bf{23.05} & 23.19 \\

      \bottomrule
    \end{tabular}
    }
    \endgroup
    \caption{BLEU scores of ST systems. \texttt{X}: original, \texttt{Y}: \texttt{WMT5M}, \texttt{Z}: \texttt{WMT10M}. For unsegmented test sets, we used \underline{pyannote.audio} with $M_{\rm dur}=2000$ and $M_{\rm int}=100$.}
    \label{tab:result_st}
\end{table*}

\subsection{MT}
In this section, we show the results of our MT systems used for cascade systems and pseudo labeling in SeqKD.
We report case-sensitive detokenized BLEU scores~\citep{papineni-etal-2002-bleu} with the \texttt{multi-bleu-detok.perl} script in Moses.
We carefully investigated the effective amount of WMT training data to improve the performance of the TED domain.
The results are shown in Table~\ref{tab:result_mt}.
We confirmed that adding the WMT data improved the performance by more than 4 BLEU.
Regarding the WMT data size, using up to 10M sentences was helpful, but 20M did not show clear improvements, probably because of the undersampling of the TED data.
Oversampling as in multilingual NMT~\citep{arivazhagan2019massively} could alleviate this problem, but this is beyond our scope.

After training with a mix of the WMT and TED data, we also tried to finetune the model with the TED data only, but this did not lead to clear improvement, especially for the IWSLT test sets.
Increasing the model capacity was not helpful, although the conclusion might change by adding more training data and evaluating the model in other domains such as news.
Because our primary focus to use MT systems was pseudo labeling for SeqKD, we decided to use the Base configuration to speed up decoding.

Finally, we checked the BLEU scores on the Must-C training data used for SeqKD.
We observed that adding more WMT data decreased the BLEU score, from which we can conclude that using more WMT data gradually changed the MT output from the TED style.
Therefore, we decided to use the models trained on \texttt{WMT5M} and \texttt{WMT10M} as teachers for SeqKD.

\subsection{Speech translation}

\subsubsection{E2E-ST}

\paragraph{SeqKD} The results are shown in Table~\ref{tab:result_st}.
We first observed the baseline Conformer model (\texttt{A1}) achieved 35.63 BLEU on the Must-C \texttt{tst-COMMON} set, and it is the new state-of-the-art record to the best of our knowledge.
Surprisingly, it even outperformed text-based MT systems in Table~\ref{tab:result_mt}.
On the other hand, unlike our observations in~\citep{inaguma2020orthros,inaguma2021source}, SeqKD (\texttt{A2-4}) degraded the performance on the Must-C \texttt{tst-COMMON} set.
However, the results on the Must-C \texttt{dev} and \texttt{tst-HE} sets showed completely different trends, where we observed better BLEU scores by SeqKD in proportion to the WMT data used for training the teachers.
Therefore, after tuning audio segmentation, we also evaluated the models on the unsegmented IWSLT test sets.
Here, we used the pyannote.audio based segmentation with $(M_{\rm dur},M_{\rm int})=(2000,100)$ as described in~\Sref{sssec:segmentation_asr}.
Then, we confirmed large improvements with SeqKD by 2-6 BLEU, and therefore we decided to determine the best model based on the IWSLT test sets.
Multi-referenced training consistently improved the BLEU scores on the IWSLT sets.
For example, \texttt{A4} outperformed \texttt{A1} by 6.02 BLEU on tst2019 although the tst2019 set was well-segmented (WER: 6.0\%).
Given these observations, we recommend evaluating ST models on multiple test sets for future research.

\paragraph{Multi-Decoder architecture}
We combined the SeqKD and Multi-Decoder techniques in our \texttt{B1} system.
\texttt{B1}, which used a conformer ASR encoder and \texttt{2ref} SeqKD, showed an improvement of 2.19 BLEU on tst2019 over \texttt{A3}, the encoder-decoder which also used \texttt{2ref} SeqKD.
\texttt{B1} also achieved a slightly higher result on tst2019 compared to \texttt{A4} which used \texttt{3ref} SeqKD.
These results suggest that the Multi-Decoder architecture is indeed compatible with SeqKD.

\begin{table}[t]
    \centering
    \small
    \begin{tabular}{llc}\toprule
      ID & Ensembled Models & tst2019 \\
      \midrule
      \phantom{0}- & \texttt{B1} & 21.06 \\
      \texttt{E1} & \texttt{B1}, \texttt{A4} & 22.51 \\
      \texttt{E2} & \texttt{B1}, \texttt{A4}, \texttt{A1} & 22.83 \\
      \texttt{E3} & \texttt{B1}, \texttt{A4}, \texttt{A1}, \texttt{A3} & 23.36 \\
      \texttt{E4} & \texttt{B1}, \texttt{A4}, \texttt{A1}, \texttt{A3}, \texttt{A2} & \textbf{23.61} \\
      \bottomrule
    \end{tabular}
    \caption{BLEU ($\uparrow$) scores of ensembled E2E-ST systems on tst2019, using the \underline{provided} segmentation with $M_{\rm dur}=2000$ and $M_{\rm int}=100$}
    \label{tab:result_ensembling}
\end{table}

\begin{table}[t]
    \centering
    \begingroup
    \scalebox{0.63}{
    \begin{tabular}{lccccccc}\toprule
    \multirow{2}{*}{VAD} & \multirow{2}{*}{$M_{\rm dur}$} & \multirow{2}{*}{$M_{\rm int}$} & \multicolumn{5}{c}{\bf{BLEU} ($\uparrow$)} \\
    \cmidrule(lr){4-7} \cmidrule(lr){8-8}
       &  &  & tst2010 & tst2015 & tst2018 & tst2019 & Avg. \\ \midrule
       Provided${}^\dagger$ & -- & -- & -- & -- & -- & 20.1\phantom{0} & -- \\ \midrule

      \multirow{4}{*}{\shortstack{Provided\\(E2E)}}
      & -- & -- & 21.99 & 19.94 & 19.29 & 19.70 & 20.23 \\
       & 1000 & 200 & 22.62 & 20.54 & 19.80 & 20.54 & 20.88 \\
       & 1500 & 200 & 23.00 & 21.66 & 20.14 & 21.50 & 21.58 \\
       & 2000 & 200 & 22.95 & 21.58 & 20.03 & 21.34 & 21.48 \\
      \midrule

       \multirow{4}{*}{\shortstack{WebRTC\\(E2E)}}
       & -- & -- & 13.13 & 12.97 & 11.07 & 13.32 & 12.62 \\
       & 1000 & 200 & 20.95 & 20.66 & 17.09 & 20.87 & 19.89 \\
       & 1500 & 200 & 21.00 & 20.99 & 17.67 & 21.05 & 20.18 \\
       & 2000 & 200 & 20.25 & 21.81 & 17.08 & 20.71 & 19.96 \\
       \midrule

       \multirow{7}{*}{\shortstack{pyannote\\(E2E)}}
       & -- & -- & 22.26 & 16.84 & 17.78 & 19.98 & 19.22 \\
       & 1500 & 200 & 25.00 & 22.22 & 21.97 & 22.67 & 22.97 \\
       & 1500 & 100 & \bf{25.92} & \bf{22.81} & \bf{22.51} & 22.88 & \bf{23.53} \\
       & 2000 & 200 & 24.10 & 21.98 & 21.00 & 22.71 & 22.45 \\
       & 2000 & 150 & 24.25 & 22.26 & 21.41 & 22.99 & 22.73 \\
       & 2000 & 100 & 24.97 & 22.66 & 22.20 & \bf{23.41} & 23.31 \\
       & 2000 & \phantom{0}50 & 24.50 & 20.67 & 22.14 & 22.89 & 22.55 \\ \midrule \midrule

       \multirow{6}{*}{\shortstack{pyannote\\(Cascade)}}
       & 1500 & 200 & 25.06 & 22.65 & 23.01 & 22.51 & 23.31 \\
       & 1500 & 100 & \bf{25.56} & 22.85 & 23.03 & 22.82 & 23.57 \\
       & 2000 & 200 & 24.41 & 22.76 & 22.15 & 22.08 & 22.85 \\
       & 2000 & 150 & 24.50 & 23.03 & \bf{23.12} & 23.11 & 23.44 \\
       & 2000 & 100 & 25.04 & \bf{23.17} & 23.05 & \bf{23.19} & \bf{23.61} \\
       & 2000 & \phantom{0}50 & 24.33 & 20.79 & \bf{23.12} & 23.11 & 22.84 \\

      \bottomrule
    \end{tabular}
    }
    \endgroup
    \caption{Impact of audio segmentation for ST. \texttt{A4} was used for the E2E model. ${}^\dagger$~\citep{potapczyk-przybysz-2020-srpols}}
    \label{tab:result_st_segmenation}
\end{table}

\paragraph{Model ensemble}
As shown in Table~\ref{tab:result_ensembling}, ensembling our various ST systems using the posterior combination method described in~\Sref{sec:ensembling} showed improvements over the best single model, \texttt{B1}.
We found that an ensemble of all of our models, \texttt{A1-4} and \texttt{B1}, achieved the best result of 23.61 BLEU on tst2019 and outperformed \texttt{B1} by 2.55 BLEU.
Although \texttt{A1} as a single system performs worse on tst2019 than the other single systems as shown in Table~\ref{tab:result_st}, including it in an ensemble with the two best single systems, \texttt{B1} and \texttt{A4}, still yielded a slight gain of 0.32 BLEU (\texttt{E2}).
Therefore, we can conclude that weak models are still beneficial for ensembling.

\begin{table*}[t]
    \centering
    \small
    \begin{tabular}{lcccccccc}\toprule
      \multirow{3}{*}{System} & \multirow{3}{*}{Segmentation} & \multirow{3}{*}{\shortstack{Segment\\merging}} & \multirow{3}{*}{$M_{\rm int}$} &  \multicolumn{5}{c}{\bf{BLEU} ($\uparrow$)} \\
    \cmidrule(lr){5-6} \cmidrule(lr){7-9}
      &  &  & & \multirow{1}{*}{tst2019} & \multirow{1}{*}{tst2020} & \multicolumn{3}{c}{tst2021} \\
       &  & & & & & ref1 & ref2 & both \\
      \midrule
      \multirow{2}{*}{IWSLT'20 winner${}^{\clubsuit}$} & given & -- & -- & 20.1\phantom{0} & 21.5 & -- & -- & -- \\
      & own & -- & -- & 23.96 & 25.3 & -- & -- & -- \\
      \midrule

      \texttt{E4} (primary) & pyannote & \cmark & 200 & \bf{24.14} & \bf{25.6} & \bf{19.3} & \bf{21.2} & \bf{31.4} \\ %
      \midrule
      \texttt{E4}+* & pyannote & \cmark & 200 & 24.41 & 25.5 & \bf{19.7} & 20.6 & 30.8 \\  %
      \texttt{E4}+* & pyannote & \cmark & 100 & \bf{24.87} & \bf{26.0} & 19.5 & 21.1 & 31.3 \\  %

      \texttt{E4}+* & given & \cmark & 100 & 23.72 & 25.1 & 19.4 & \bf{21.4} & \bf{31.5} \\ %
      \texttt{E4}+* & given & \xmark & -- & 21.10 & 22.3 & 17.4 & 18.4 & 27.7 \\ %
      \texttt{B1} & pyannote & \cmark & 100 & 23.78 & 25.0 & 18.9 & 20.9 & 31.1 \\ %

      \bottomrule
    \end{tabular}
    \caption{BLEU scores of submitted systems on tst2020 and tst2021. ${}^{\clubsuit}$~\citep{potapczyk-przybysz-2020-srpols}. $M_{\rm dur} =2000$ was used for the segment merging algorithm. *Late submission (not official). \texttt{E4}+ denotes \texttt{E4} trained for more steps.}
    \label{tab:result_final}
\end{table*}

\subsubsection{Segmentation}
\label{sec:segmentation_results}
Similar to~\Sref{sssec:segmentation_asr}, we also investigated the impact of audio segmentation for E2E-ST models.
To this end, we used the \texttt{A4} model.
Note that we used the same decoding hyperparameters tuned on Must-C.
The results are shown in Table~\ref{tab:result_st_segmenation}.
We confirmed a similar trend to ASR.
Although $(M_{\rm dur},M_{\rm int})=(1500,100)$ showed the best performance on average, we decided to use $(M_{\rm dur},M_{\rm int})=(2000,100)$ for submission considering the best performance on the latest IWSLT test, tst2019.

\subsubsection{Cascade system}
We also evaluated the cascade system with the Conformer ASR and the Transformer-Base MT trained on the \texttt{WMT10M} data (\texttt{C1}).
The MT model was trained by feeding source sentences without case information and punctuation marks.
The results in Table~\ref{tab:result_st_segmenation} showed that the BLEU scores correlated to the WER in Table~,\ref{tab:result_asr_segmenation} and the performance was comparable with that of \texttt{A4}.
Although there is some room for improving the performance of the cascade system further by using in-domain English LM, it is difficult to conclude which modeling (cascade or E2E) is effective because the cascade system had more model parameters in the ASR decoder and MT encoder.
This means that the E2E model could also be enhanced by using a similar amount of parameters.

\subsubsection{Final system}
Our final system was the best ensemble system \texttt{E4}, using the pyannote.audio based segmentation with $(M_{\rm dur},M_{\rm int})=(2000,200)$\footnote{Because of time limitation, we submitted the systems before completing tuning segmentation hyperparameters.}.
This system, which was our primary submission, scored 24.14 BLEU on tst2019 as shown in Table~\ref{tab:result_final}.
Compared to the result in Table~\ref{tab:result_ensembling}, it improved by 0.53 BLEU thanks to better audio segmentation.
It was also slightly higher than the IWSLT20 winner's submission by SPROL~\citep{potapczyk-przybysz-2020-srpols}.

We also present the results on tst2020 and tst2021 in Table~\ref{tab:result_final}.
Our primary submission \texttt{E4} outperformed the result of last year's winner system on tst2020.

\section{Conclusion}
In this paper, we have presented the ESPnet-ST group's offline systems on the IWSLT 2021 submission.
We significantly improved the baseline Conformer performance with multi-referenced SeqKD, Multi-Decoder architecture, segment merging algorithm, and model ensembling.
Our future work includes scaling training data and careful analysis of the performance gap in different test sets.

\section{Acknowledgement}
This work was partly supported by ASAPP and JHU HLTCOE.
This work used the Extreme Science and Engineering Discovery Environment (XSEDE)~\citep{towns2014xsede}, which is supported by National Science Foundation grant number ACI-1548562.
Specifically, it used the Bridges system~\citep{nystrom2015bridges}, which is supported by NSF award number ACI-1445606, at the Pittsburgh Supercomputing Center (PSC).

\newpage
\bibliographystyle{acl_natbib}


\section{Device}

\begin{figure*}[htbp]
\center{\includegraphics[width=1\linewidth]{figures/Figure_1.eps}}
\caption{\textbf{a.} Scanning electron microscope images of the silicon quantum processor showing the quantum dot gate pattern and the micromagnet on top. \textbf{b.} Control paths for determining the symmetry operation point in the charge stability diagram. ($M$,$N$) represent the number of electrons in the dots underneath the tip of LP and RP respectively. \textbf{c.} Pulse sequence schematic of a decoupled controlled-phase operation interleaved in a Ramsey interference sequence on Q1. \textbf{d.} Spin-up probability of Q1 after the Ramsey sequence in \textbf{c}, as a function of the detuning in the double dot potential and the total duration of the barrier voltage pulses.}
\label{fig:device}
\end{figure*}

We use a gate-defined double quantum dot in an isotopically enriched $^{28}$Si/SiGe heterostructure (Fig.~\ref{fig:device}a)~\cite{xue2020cmos}, with each dot occupied by one single electron (see Methods). The spin states of the electrons serve as qubits. The spin states are measured with the help of a sensing quantum dot (SQD), which is capacitively coupled to the qubit dots~\cite{elzerman2004single}. A micromagnet on top of the device provides a magnetic field gradient enabling electric-dipole spin resonance (EDSR)~\cite{pioro2008electrically}, and separates the resonance frequencies of the qubits in the presence of an external magnetic field ($\sim$\unit[320]{mT}) to \unit[11.993]{GHz} (Q1) and \unit[11.890]{GHz} (Q2). Single-qubit $X$ and $Y$ gates are implemented by frequency-multiplexed microwave signals applied to gate MW, and virtual $Z$ gates are implemented by a phase update of the reference frame~\cite{vandersypen2005nmr}. The plunger gates (LP and RP) control the chemical potentials of the quantum dots.\par

The native two-qubit gate for spin qubits utilizes the exchange interaction~\cite{loss1998quantum, petta2005coherent}, originating from the wave-function overlap of electrons in neighbouring dots. This selectively shifts the energy of the anti-parallel spin states and thus allows for an electrically pulsed adiabatic \textsc{cphase} gate~\cite{meunier2011efficient,veldhorst2015two,watson2018programmable}. The barrier gate (B) controls the tunnel coupling between the dots, allowing to precisely tune the exchange coupling from < \unit[100]{kHz} to \unit[20]{MHz}. In order to minimize the sensitivity to charge noise, we activate the exchange coupling while avoiding a tilt in the double dot potential~\cite{martins2016noise,reed2016reduced,zhang_2017_tilt} (Fig.~\ref{fig:device}a).
This symmetric condition can be determined accurately by decoupled adiabatic exchange pulses inside a Ramsey sequence (Fig.~\ref{fig:device}c-d).
The tunnel barrier is controlled by simultaneously pulsing gate B and compensating on LP and RP to avoid shifts of the electrochemical potentials~\cite{martins2016noise}. The detuning between quantum dots is controlled by additional offsets on the LP and RP pulses in opposite directions. As the decoupling pulses remove additional single-qubit phase accumulation from electron movement in the magnetic field gradient, the spin-up probability of Q1 results in a symmetric Chevron pattern, with the symmetry point at the center (Fig.~\ref{fig:device}e).\par

\section{Conclusion}

The two-qubit gate with fidelity above 99.5\% and single-qubit gate fidelities in the two-qubit gate space above 99\% on average, place semiconductor spin qubit logic at the error threshold of the surface code. Recently, a two-qubit operation between nuclear spin qubits in silicon, mediated by an electron spin qubit, has been demonstrated to surpass 99\% fidelity as well, further highlighting that semiconductor spin qubits offer precise two-qubit logic~\cite{madzik2021precision}. Independent studies have shown spin qubit readout with a fidelity above 98\% in only a few $\mu$s~\cite{zheng2019rapid}, with further improvements underway~\cite{schaal2020fast}. With a modest effort in reducing crosstalk errors and in extending the device designs, we are optimistic that the individually demonstrated advantages of semiconductor spin qubits can be combined into a fault-tolerant and highly-integrated quantum computer. The same advances will allow us to implement more sophisticated algorithms in the NISQ era, such as solving energies involving excited states of more complex molecules.

\textbf{Acknowledgements}
We acknowledge fruitful discussions with P. Cerfontaine, C. Bureau-Oxton, M. Madzik, A. Morello, J. Helsen, B. Terhal, M. Veldhorst, and all the members of the spin qubit team, and technical assistance by O. Benningshof, M. Sarsby, R. Schouten and R. Vermeulen. This research was funded by the Dutch Ministry for Economic Affairs through the allowance for Topconsortia for Knowledge and Innovation (TKI) and the Army Research Office (ARO) under grant numbers W911NF-17-1-0274. The views and conclusions contained in this document are those of the authors and should not be interpreted as representing the official policies, either expressed or implied, of the ARO or the US Government. The US Government is authorized to reproduce and distribute reprints for government purposes notwithstanding any copyright notation herein.

\textbf{Author contributions} X.X. performed the experiment with help from N.S. and B.U., M.R developed the theory model and analyzed the data with X.X., N.S. fabricated the quantum dot device, A.S. and G.S. designed and grew the Si/SiGe heterostructure, X.X. and L.M.K.V. conceived the project, L.M.K.V. supervised the project, X.X., M.R. and L.M.K.V. wrote the manuscript with input from all authors.

\textbf{Competing interests}
The authors declare no competing interests.

\textbf{Data availability}
Data supporting this work are available at zenodo, https://doi.org/10.5281/zenodo.5044450.

\textbf{Code availability}
The codes used for data acquisition and processing are from the open source python packages QCoDeS, which is available at https://github.com/QCoDeS/Qcodes, QTT, which is available at https://github.com/QuTech-Delft/qtt, and PycQED, which is available at https://github.com/DiCarloLab-Delft/PycQED{\textunderscore}py3. The codes used for the design and analysis of the gate set tomography experiment are from pyGSTi, which is available at https://github.com/pyGSTio/pyGSTi. And the codes used for the design and analysis of the variational quantum eigensolver experiment are from OpenFermion, which is available at https://github.com/quantumlib/OpenFermion.

\section{Fitting qubit frequencies and dephasing times}
The transition energies Eqs.~\eqref{eq:resonance1}-\eqref{eq:resonance4} are fitted simultaneously to the measured results from the Ramsey experiment (see Fig.~\ref{fig:pulse_engineer}a). For the fitting we use the \texttt{NonLinearModelFit} function from the software Mathematica with the least square method. The best fits yield the following parameters $\alpha=\unit[12.1\pm 0.05]{V^{-1}}$, $\beta_1=\unit[-2.91\pm 0.11]{MHz^\gamma/V^\gamma}$, $\beta_2=\unit[67.2\pm 0.63]{MHz^\gamma/V^\gamma}$, and $\gamma=\unit[1.20\pm 0.01]{}$, and $J_\text{res}=\unit[58.8\pm 1.8]{kHz}$.\par

The dephasing times Eqs.~\eqref{eq:T2star_1}-\eqref{eq:T2star_4} are fitted simultaneously to the measured results from the Ramsey experiment (see Fig.~\ref{fig:pulse_engineer}c) using the same method. The best fits yield the following parameters $\delta v_B=\unit[0.40\pm 0.01]{m V}$, $\delta f_{Q1}=\unit[11\pm 0.1]{kHz}$, and $\delta f_{Q2}=\unit[24\pm 0.7]{kHz}$.
\section{Gate calibration}

In the gate set used in this work, $\{I$, $X_{Q1}$, $Y_{Q1}$, $X_{Q2}$, $Y_{Q2}$, $\textsc{cphase} \}$, the duration of the $I$ gate and the \textsc{cphase} are set to \SI{100}{\nano\second}, and we calibrate and keep the amplitudes of the single-qubit drives fixed and in the linear response regime. The envelope of the single-qubit gates are shaped following a ``tukey'' window, as it allows adiabatic single-qubit gates with relatively small amplitudes, thus avoids distortion caused by nonlinear response.
The general tukey window of length $t_p$ is given by
\begin{align}
    W(t,r) = \begin{cases}
      \frac{1}{2}\left[1-\cos\left(\frac{2\pi t }{r t_p}\right)\right] & 0 \leq t \leq \frac{r t_p}{2} \\
      1 & \frac{r t_p}{2} < t\leq t_p -\frac{r t_p}{2} \\
       \frac{1}{2}\left[1-\cos\left(\frac{2\pi (t_p-t) }{r t_p}\right)\right] & t_p -\frac{r t_p}{2}\leq t \leq t_p,
   \end{cases}
   \label{eq:tukey_window}
\end{align}
where $r=0.5$ for our pulses.
Apart from these fixed parameters, there are 11 free parameters that must be calibrated: single-qubit frequencies $f_{Q1}$ and $f_{Q2}$, burst lengths for single-qubit gates $t_{XY1}$ and $t_{XY2}$, phase shifts caused by single-qubit gates on the addressed qubit itself $\phi_{11}$ and $\phi_{22}$, phase shifts caused by single-qubit gates on the unaddressed ``victim qubit'' $\phi_{12}$ and $\phi_{21}$ ($\phi_{12}$ is the phase shift on Q1 induced by a gate on Q2 and similar for $\phi_{21}$), the peak amplitude of the \textsc{cphase} gate $A_{v_B}$, and phase shifts caused by the gate voltage pulses used for \textsc{cphase} gate on the qubits $\theta_1$ and $\theta_2$ (in addition, we absorb into $\theta_1$ and $\theta_2$ the 90 degree phase shifts needed to transform $\text{diag}(1,i,i,1)$ into $\text{diag}(1,1,1,-1)$).\par
For single-qubit gates, $f_{Q1}$ and $f_{Q2}$ are calibrated by standard Ramsey sequences, which are automatically executed every two hours. The EDSR burst times $t_{XY1}$ and $t_{XY2}$ are initially calibrated by an AllXY calibration protocol~\cite{reed2013entanglement}. The phases $\phi_{11}$, $\phi_{12}$, $\phi_{21}$, and $\phi_{22}$ are initially calibrated by measuring the phase shift of the victim qubit (Q1 for $\phi_{11}$ and $\phi_{21}$; Q2 for $\phi_{22}$ and $\phi_{12}$) in a Ramsey sequence interleaved by a pair of [$X_{Qi}$, $-X_{Qi}$] gates on the addressed qubit (Q1 for $\phi_{11}$ and $\phi_{12}$; Q2 for $\phi_{22}$ and $\phi_{21}$) (Extended Data Fig.~\ref{fig:Manual_calibration}). \par
The optimal pulse design presented in Fig.~\ref{fig:pulse_engineer} gives a rough guidance of the pulse amplitude $A_{v_B}$. In a more precise calibration of the \textsc{cphase} gate, an optional $\pi$-rotation is applied to the control qubit (e.g. Q1) to prepare it into the $\ket{0}$ or $\ket{1}$ state, followed by a Ramsey sequence on the target qubit (Q2) interleaved by an exchange pulse. The amplitude is precisely tuned to bring Q2 completely out of phase (by 180 degree) between the two measurements (Extended Data Fig.~\ref{fig:Manual_calibration} d-e). The phase $\theta_2$ is determined such that the phase of Q2 changes by zero ($\pi$) when Q1 is in the state $\ket{0}$ ($\ket{1}$),  corresponding to the $\textsc{cphase}=\text{diag}(1,1,1,-1)$ in the standard basis. The same measurement is then performed again with the Q2 as the control qubit and Q1 as the target qubit to determine $\theta_1$~\cite{watson2018programmable}.\par
In such a ``conventional'' calibration procedure of the \textsc{cphase} gate, we notice that the two qubits experience different conditional phases (Extended Data Fig.~\ref{fig:Manual_calibration}). We believe that this effect is caused by off-resonant driving from the optional $\pi$-rotation on the control qubit. Similar effects can also affect the calibration of the phase crosstalk from single-qubit gates.\par

This motivates us to use the results from GST as feedback to adjust the gate parameters. The error generators not only describe the total errors of the gates, but also distinguish Hamiltonian errors (coherent errors) from stochastic errors (incoherent errors). We use the information on 7 different Hamiltonian errors ($IX$, $IY$, $XI$, $YI$, $ZI$, $IZ$ and $ZZ$) of each gate, to correct all 11 gate parameters, except $f_{Q1}$ and $f_{Q2}$, for which calibrations using standard Ramsey sequences are sufficient.
For single-qubit gates, $t_{XY1}$ and $t_{XY2}$ are adjusted according to the $IX$, $IY$, $XI$ and $YI$ errors. The phases $\phi_{11}$, $\phi_{12}$, $\phi_{21}$, and $\phi_{22}$ are adjusted according to the $ZI$ and $IZ$ errors. For the \textsc{cphase} gate, $\theta_1$ and $\theta_2$ are adjusted according to the $ZI$ and $IZ$ errors, and $A_{v_B}$ is adjusted according to the $ZZ$ error. The adjusted gates are then used in a new GST loop.
\section{Measurement setup}

The measurement setup and device are similar to the one used in Ref.~\cite{xue2020cmos}. We summarize a few key points and all the differences here. The gates LP, RP, and B are connected to arbitrary waveform generators (AWGs, Tektronix 5014C) via coaxial cables. The position in the charge stability diagram of the quantum dots is controlled by voltage pulses applied on LP and RP. Linear combinations of the voltage pulses on B, LP and RP are used to control the exchange coupling between the two qubits at the symmetry point. The compensation coefficients are: $v_{LP}/v_B = -0.081, v_{RP}/v_B = -0.104$. A vector signal generator (VSG, Keysight E8267D) is connected to gate MW and sends frequency-multiplexed microwave bursts (not necessarily time-multiplexed) to implement electric-dipole spin resonance (EDSR). The VSG has two I/Q input channels, receiving I/Q modulation pulses from two channels of an AWG. I/Q modulation is used to control the frequency, phase, and length of the microwave bursts. The current signal of the sensing quantum dot is converted to a voltage signal and recorded by a digitizer card (Spectrum M4i.44), and then converted into 0 or 1 by comparing it to a threshold value.\par

Two differences between the present setup and the one in Ref.~\cite{xue2020cmos} are that 1)  the programmable mechanical switch is configured such that gate MW is always connected to the VSG, and not to the cryo-CMOS control chip; 2) a second AWG of the same model is connected to gate B with its clock synchronized with the first AWG.\par

\section{Numerical simulations}
\label{sec:num_sim}

For all numerical simulations,  we solve the time-dependent Schr{\"o}dinger equation
\begin{align}
    i\hbar \ket{\psi(t)} = H \ket{\psi(t)}
\end{align}
and iteratively compute the unitary propagator according to
\begin{align}
    U(t+\Delta t) = e^{-\frac{i}{\hbar} H(t+\Delta t)}U(t).
\end{align}
Here, $H(t+\Delta t$ is discretized into $N$ segments of length $\Delta t$ such that $H(t)$ is constant in the time-interval $\left[t,t+\Delta t\right]$.
All simulations are performed in the rotating frame of the external magnetic field $(B_{z,1}+B_{z,2})/2$ and neglecting the counter-rotating terms, making the so-called rotating wave approximation (RWA). This allows us to chose $\Delta t=\unit[10]{ps}$ as a sufficiently small time step.

For the noise simulations, we included classical fluctuations of $f_{Q1}\rightarrow f_{Q1} + \delta f_{Q1}$, $f_{Q2}\rightarrow f_{Q2} + \delta f_{Q2}$, and $v_B\rightarrow v_B + \delta v_B$. We assume the noise coupling to the resonance frequencies $\delta f_{Q1}$ and $\delta f_{Q2}$ to be quasistatic and assume $1/f$ noise for $v_B$ which we describe by its spectral density $S_{\delta v_B}(\omega)= \delta v_B/\omega$. To compute time traces of the fluctuation we use the approach introduced in Refs.~\cite{yang_achieving_2019,koski_strong_2020} to generate time-correlated time traces. The fluctuations are discretized into $N$ segments with time $\Delta t$ such that $\delta v_B(t)$ is constant in the time interval $\left[t,t+\Delta t\right)$, with the same $\Delta t$ as above. Consequently, fluctuations which are faster than $f_\text{max}=\frac{1}{\Delta t}$ are truncated.
\section{\textsc{cphase} gate}
We realize a universal $\textsc{cphase}=\text{diag}(1,1,1,-1)$ gate by adiabatically pulsing the exchange interaction using a carefully designed pulse shape. Starting from Eq.~\eqref{eq:ham_matrix}, the full dynamics can be projected on the odd-parity space spanned by $\ket{01}$ and $\ket{10}$. The entangling exchange gate is reduced in this subspace to a global phase shift thus the goal is to minimize any dynamics inside the subspace. Introducing a new set of Pauli operators in this subspace
$\sigma_x=\ket{01}\bra{10}+\ket{10}\bra{01}$, $\sigma_y=-i \ket{01}\bra{10}+i\ket{10}\bra{01}$, and $\sigma_z=\ket{01}\bra{01}-\ket{10}\bra{10}$, we find
\begin{align}
H_{\text{sub}}(t)=\frac{1}{2}\big( - J\left(v_B(t)\right) + \Delta E_z \,\sigma_z+J\left(v_B(t)\right)\, \sigma_x\big).
\label{eq:HheisSub}
\end{align}
In order to investigate the adiabatic behaviour, it is convenient to switch into the adiabatic frame defined by $U_\text{ad}=e^{-\frac{i}{2} \tan^{-1}\left(\frac{J\left(v_B(t)\right)}{\Delta E_z}\right)\sigma_y}$. The Hamiltonian accordingly transforms as
\begin{align}
H_{\text{ad}}&=U_\text{ad}^\dagger(t) H_{\text{sub}}(t) U_\text{ad}(t) -i \hbar U_\text{ad}^\dagger(t)\dot{U}_\text{ad}(t)\label{eq:Ham_ramp_part1}\\
&\approx \frac{1}{2}\bigg(- J\left(v_B(t)\right) +  \Delta E_z\,\sigma_z
-\frac{\hbar\dot{J}}{\Delta E_z}\sigma_y \bigg),
\label{eq:Ham_ramp_part2}
\end{align}
where the first term is unaffected and describes the global phase accumulation due to the exchange interaction, the second term describes the single-qubit phase accumulations, and the last term $f(t)=\hbar\dot{J}/(2\Delta E_z)$ describes the diabatic deviation proportional to the derivative of the exchange pulse. From Eq.~\eqref{eq:Ham_ramp_part1} to Eq.~\eqref{eq:Ham_ramp_part2} we assumed a constant $\Delta E_z(t)\approx \Delta E_z$, and $J(t)\ll \Delta E_z$. The transition probability from state $\ket{\uparrow\downarrow}$ to $\ket{\downarrow\uparrow}$ using a pulse of length $t_p$ is then given by~\cite{martinis_fast_2014}
\begin{align}
P_{\ket{\uparrow\downarrow}\rightarrow\ket{\downarrow\uparrow}} &\approx \left|\int_{0}^{t_p} f(t) e^{-\frac{i}{\hbar} \Delta E_z t }\, dt \right|^2 \label{eq:ESD_part1} \\
&= S\big(f(t)\big).
\label{eq:ESD}
\end{align}
From the first to the second line we replaced the integral by the (short time-scale) Fourier transform, allowing us to describe the spin-flip error probability by the energy spectral density (ESD) of the input signal $f(t)$. Minimizing such errors is therefore identical to minimizing the ESD of a pulse, a well-known and solved problem from classical signal processing and statistics. Optimal shapes are commonly referred to as window functions $W(t)$ due to their property to restrict the spectral resolution of signals.
A high-fidelity exchange pulse is consequently given by $J(0)=J(t_p)$ and
\begin{align}
    \int_0^{t_p}dt J\left(v_B(t)\right)/\hbar=\pi/2,
    \label{eq:cphase_condition}
\end{align}
while setting $J(t)= A_{v_B} W(t) J_\text{res}$~\cite{martinis_fast_2014}, with a scaling factor $A_{v_B}$ that is to be determined. In this work, we have chosen the cosine window
\begin{align}
    W(t)=\frac{1}{2}\left[1-\cos\left(\frac{2\pi t }{t_p}\right)\right],
\end{align}
from signal processing which has a high spectral resolution. The amplitude $A_{v_B}$ follows from condition Eq.~\eqref{eq:cphase_condition}. For a pulse length of $t_p=\unit[100]{ns}$ and a cosine pulse shape we find $A_{v_B} J_\text{res}=\unit[10.06]{MHz}$. Our numerical simulations predict an average gate infidelity $1-F_\text{gate}<10^{-6}$ without noise and $1-F=0.22\times 10^{-3}$ with the inclusion of noise through the fluctuations $\delta f_{Q1}$, $\delta f_{Q2}$, and $\delta v_B$ discussed in section~\ref{sec:num_sim}.
As explained in the main text, due to the exponential voltage-exchange relation the target pulse shape for J(t) must be converted to a barrier gate pulse, following~\cite{russ_soon_to_be_appear}
\begin{align}
    v_B(t)=\frac{1}{2\alpha}\text{log}\left(A_{v_B}\, W(t)\right).
\end{align}
The numerical simulations with the fitted noise parameters in the simplified model from Section~\ref{sec:num_sim} predict a gate fidelity above 99.97\%. The measured PTMs reveal significantly higher rates of incoherent errors, which we attribute to drifts in the barrier voltage on a timescale much longer than the timescale on which $\delta f_{Q1}$, $\delta f_{Q2}$, and $\delta v_B$ were determined.

\section{Theoretical model}
In this section we describe the theoretical model used for the fitting, the pulse optimization, and the numerical simulations.
The dynamics of two electron spins in the $(1,1)$ charge configuration can be well-described by an extended Heisenberg model~\cite{loss1998quantum}
\begin{align}
H &=J\,(\vv{S}_1\cdot \vv{S}_2-\frac{1}{4}) +  g\mu_B\vv{B}_1\cdot\vv{S}_1+  g\mu_B\vv{B}_2\cdot\vv{S}_2,
\label{eq:ham_matrix}
\end{align}
with $\vv{S}_j = \hbar(X_j,Y_j,Z_j)^T/2$, where $X_j,Y_j,Z_j$ are the single-qubit Pauli-matrices acting on spin $j=1,2$, $\mu_B$ the Bohr's magneton, $g\approx 2$ the g-factor in silicon, and $\hbar=h/(2\pi)$ the reduced Planck constant. The first and second term describe the interaction of the electron spin in dot 1 and dot 2 with the magnetic fields $\vv{B}_j=(B_{x,j},0,B_{z,j})^T$ originating from the externally applied field and the micromagnet. The transverse components $B_{x,j}$ induce spin-flips, thus, single-qubit gates if modulated resonantly via EDSR.
For later convenience we define the resonance frequencies $hf_{Q1}=g\mu_B B_{z,1}$ and $hf_{Q2}=g\mu_B B_{z,2}$, and the energy difference between the qubits $\Delta E_z= g\mu_B (B_{z,2}- B_{z,1})$. The last term in the Hamiltonian of Eq.~\eqref{eq:ham_matrix} describes the exchange interaction $J$ between the spins in neighboring dots. The exchange interaction originates from the overlap of the wave-functions through virtual tunneling events and is in general a non-linear function of the applied barrier voltage $v_B$.
We note that $v_B$ determines the compensation pulses on LP and RP for virtual barrier control.
We model $J$ as an exponential function~\cite{cerfontaine2020high,pan_resonant_2020}
\begin{align}
J(v_B) = J_\text{res} e^{2\alpha v_B},
\label{eq:J_sim_expression}
\end{align}
where $J_\text{res}\approx\unit[20-100]{kHz}$ is the residual exchange interaction during idle and single-qubit operations and $\alpha$ the lever arm.
In general the magnetic fields $\vv{B}_j$ depend on the exact position of the electron. We include this in our model $B_{z,j}\rightarrow B_{z,j}(v_B)= B_{z,j}(0) + \beta_j v_B^\gamma$, where $\beta_j$ amounts for the impact of the barrier voltage on the resonance frequency of qubit $j$. The transition energies described in the main text are now given by diagonalizing Hamiltonian from Eq.~\eqref{eq:ham_matrix} and computing the energy difference between the eigenstates corresponding to the computational basis states $\lbrace\ket{00},\ket{01},\ket{10},\ket{11}\rbrace$~\cite{russ_high-fidelity_2018}. We have
\begin{align}
    hf_\text{Q1 (Q2$=\ket{0}$)} &= \mathcal{E}(\ket{10}) - \mathcal{E}(\ket{00}),\label{eq:resonance1} \\
    hf_\text{Q1 (Q2$=\ket{1}$)} &= \mathcal{E}(\ket{11}) - \mathcal{E}(\ket{01}),\label{eq:resonance2} \\
    hf_\text{Q2 (Q1$=\ket{0}$)} &= \mathcal{E}(\ket{01}) - \mathcal{E}(\ket{00}),\label{eq:resonance3} \\
    hf_\text{Q2 (Q1$=\ket{1}$)} &= \mathcal{E}(\ket{11}) - \mathcal{E}(\ket{10}),\label{eq:resonance4}
\end{align}
where $\mathcal{E}(\ket{\xi})$ denotes the eigenenergy of eigenstate $\ket{\xi}$ and $\ket{0}=\ket{\downarrow}$ is defined by the magnetic field direction.\par

In the presence of noise, qubits start to loose information. In silicon, charge noise and nuclear noise are the dominating noise sources. In the absence of two-qubit coupling and correlated charge noise, both qubits decohere largely independently of each other, giving rise to a decoherence time set by the interaction with the nuclear spins and charge noise coupling to the qubit via intrinsic and artificial (via the inhomogeneous magnetic field) spin-orbit interaction. We describe this effect by $f_{Q1}\rightarrow f_{Q1} + \delta f_{Q1}$ and $f_{Q2}\rightarrow f_{Q2} + \delta f_{Q2}$, where $\delta f_{Q1}$ and $\delta f_{Q1}$ are the single-qubit frequency fluctuations. Charge noise additionally can affect both qubits via correlated frequency shifts and the exchange interaction through the barrier voltage, which we model as $v_B\rightarrow v_B + \delta v_B$. In the presence of finite exchange coupling one can define four distinct pure dephasing times, each corresponding to the dephasing of a single qubit with the other qubit in a specific basis state. In a quasistatic approximation the four dephasing times are then given by
\begin{widetext}

\begin{align}
    T_2^\star(\text{Q1 (Q2$=\ket{0}$)}) &= \frac{1}{\sqrt{2}\pi\sqrt{\left[\frac{d\left(hf_\text{Q1 (Q2$=\ket{0}$)})\right)}{d v_B}\right]^2\delta v_B^2
    + \left[\frac{d\left(hf_\text{Q1 (Q2$=\ket{0}$)}\right)}{d hf_{Q1}}\right]^2 \delta f_{Q1}^2
    +\left[\frac{d\left(hf_\text{Q1 (Q2$=\ket{0}$)}\right)}{d hf_{Q2}}\right]^2\delta f_{Q2}^2 }},\label{eq:T2star_1} \\
    T_2^\star(\text{Q1 (Q2$=\ket{1}$)}) &= \frac{1}{\sqrt{2}\pi\sqrt{\left[\frac{d\left(hf_\text{Q1 (Q2$=\ket{1}$)})\right)}{d v_B}\right]^2\delta v_B^2
    + \left[\frac{d\left(hf_\text{Q1 (Q2$=\ket{1}$)}\right)}{d hf_{Q1}}\right]^2 \delta f_{Q1}^2
    +\left[\frac{d\left(hf_\text{Q1 (Q2$=\ket{1}$)}\right)}{d hf_{Q2}}\right]^2\delta f_{Q2}^2 }}, \label{eq:T2star_2} \\
    T_2^\star(\text{Q2 (Q1$=\ket{0}$)}) &= \frac{1}{\sqrt{2}\pi\sqrt{\left[\frac{d\left(hf_\text{Q2 (Q1$=\ket{0}$)})\right)}{d v_B}\right]^2\delta v_B^2
    + \left[\frac{d\left(hf_\text{Q2 (Q1$=\ket{0}$)}\right)}{d hf_{Q1}}\right]^2 \delta f_{Q1}^2
    +\left[\frac{d\left(hf_\text{Q2 (Q1$=\ket{0}$)}\right)}{d hf_{Q2}}\right]^2\delta f_{Q2}^2 }}, \label{eq:T2star_3} \\
    T_2^\star(\text{Q2 (Q1$=\ket{1}$)}) &= \frac{1}{\sqrt{2}\pi\sqrt{\left[\frac{d\left(hf_\text{Q2 (Q1$=\ket{1}$)})\right)}{d v_B}\right]^2\delta v_B^2
    + \left[\frac{d\left(hf_\text{Q2 (Q1$=\ket{1}$)}\right)}{d hf_{Q1}}\right]^2 \delta f_{Q1}^2
    +\left[\frac{d\left(hf_\text{Q2 (Q1$=\ket{1}$)}\right)}{d hf_{Q2}}\right]^2\delta f_{Q2}^2 }}.
    \label{eq:T2star_4}
\end{align}
\end{widetext}
\section{VQE}

We follow the approach of Ref.~\cite{hempel2018quantum} to using the VQE algorithm to compute the ground state of molecular hydrogen, after mapping this state onto the state of two qubits. We include this information here for completeness.
The Hamiltonian of a molecular system in atomic units ($\hbar = 1$) reads
\begin{align}
    H = &-\sum_{i} \frac{\nabla_{R_i}^2}{2M_i} -\sum_{j} \frac{\nabla_{r_j}^2}{2} - \sum_{i,j} \frac{Q_i}{\abs{R_i - r_j}}\nonumber\\
    &+\sum_{i,j>i} \frac{Q_iQ_j}{\abs{R_i - R_j}} + \sum_{i,j>i} \frac{1}{\abs{r_i - r_j}},
\end{align}
where $R_i$, $M_i$ and $Q_i$ are the position, mass and charge of the $i$-th nuclei, and $r_j$ is the position of the $j$-th electron. The first two sums describe the kinetic energies of the nuclei and electrons, respectively. The last three sums describe the Coulomb repulsion between nuclei and electrons, nuclei and nuclei, and electrons and electrons, respectively. As we are primarily interested in the electronic structure of the molecule, and nuclear masses are a few orders of magnitude larger than the electron masses, the nuclei are treated as static point charges under the Born-Oppenheimer approximation. Consequentially, the electronic Hamiltonian can be simplified to
\begin{align}
    H_e = &-\sum_{i} \frac{\nabla_{r_i}^2}{2} - \sum_{i,j} \frac{Q_i}{\abs{R_i - r_j}} + \sum_{i,j>i} \frac{1}{\abs{r_i - r_j}}.
\end{align}
Switching into the second-quantization representation, described by fermionic creation and annihilation operators, $a_p^{\dagger}$ and $a_q$, acting on a finite basis, the Hamiltonian becomes
\begin{align}
    H_e = \sum_{pq}h_{pq}a_p^{\dagger}a_q + \sum_{pqrs}h_{pqrs}a_p^{\dagger}a_q^{\dagger}a_ra_s.
\end{align}
The anti-symmetry under exchange is retained through the anti-commutation relation of the operators.
The weights of the two sums are given by the integrals
\begin{align}
    &h_{pq} = \int d\sigma\psi_p^*(\sigma)(\frac{\nabla_{r_i}^2}{2} - \sum_{i} \frac{Q_i}{\abs{R_i - r}})\psi_q^*(\sigma), \\
    &h_{pqrs} = \int d\sigma_1d\sigma_2\frac{\psi_p^*(\sigma_1)\psi_q^*(\sigma_2)\psi_s(\sigma_1)\psi_r(\sigma_2)}{\abs{r_1 - r_2}}.
\end{align}
Such a second-quantized molecular Hamiltonian can be mapped onto qubits using the Jordan-Wigner (JW) or the Bravyi-Kitaev (BK) transformation~\cite{mcardle2020quantum}. The JW transformation directly encodes the occupation number (0 or 1) of the $i$-th spin-orbital into the state ($\ket{0}$ or $\ket{1}$) of the $i$-th qubit. The number of qubits required after JW transformation is thus the same as the number of spin-orbitals that are of interest. The BK transformation, on the other hand, encodes the information in both the occupation number and parities -- whether there is an even or odd occupation in a subset of spin-orbitals.\par

Taking molecular hydrogen in the Hartree-Fock basis as an example, we are interested in investigating the bonding ($\ket{O_1\uparrow}$, $\ket{O_1\downarrow}$) and the anti-bonding orbital state ($\ket{O_2\uparrow}$, $\ket{O_2\downarrow}$). The initial guess of the solution is the Hartree-Fock (HF) state in which both electrons occupy the $\ket{O_1}$ orbital. The JW transformation encodes the HF initial state as $\ket{0011}$, representing $\ket{N_{O_2\downarrow} N_{O_2\uparrow} N_{O_1\downarrow} N_{O_1\uparrow}}$ from left to right, where $N_{O_iS}$ is the occupation of the $O_iS$ spin-orbital with $S=\uparrow,\downarrow$. The BK transformation encodes the HF initial state as $\ket{0001}$, where the first and the third qubit (counting from the right) encode the occupation number of the first and third spin-orbital ($N_{O_1\uparrow} = 1$ and $N_{O_2\uparrow} = 0$), the second qubit encodes the parity of the first two spin-orbitals ($(N_{O_1\uparrow} + N_{O_1\downarrow})$ mod 2 = 0), and the fourth qubit encodes the parity of all four spin-orbitals ($(N_{O_1\uparrow} + N_{O_1\downarrow} + N_{O_2\uparrow} + N_{O_2\downarrow})$ mod 2 = 0). With the standard transformation rules for fermionic creation and annihilation operators, the system Hamiltonian becomes a four-qubit Hamiltonian
\begin{align}
    H_{JW} = &~g_0I + g_1Z_1 + g_2Z_2 + g_3Z_3 + g_4Z_4\nonumber\\
             &~+ g_5Z_1Z_2 + g_6Z_1Z_3 + g_7Z_1Z_4\nonumber\\
             &~+ g_8Z_2Z_3 + g_9Z_2Z_4 + g_{10}Z_3Z_4\nonumber\\
             &~+ g_{11}Y_1X_2X_3X_4 + g_{12}Y_1Y_2X_3X_4\nonumber\\
             &~+ g_{13}X_1X_2Y_3Y_4 + g_{14}X_1Y_2Y_3Y_4,\\
    H_{BK} = &~g_0I + g_1Z_1 + g_2Z_2 + g_3Z_3\nonumber\\
             &~+ g_4Z_1Z_2 + g_5Z_1Z_3 + g_6Z_2Z_4\nonumber\\
             &~+ g_7Z_1Z_2Z_3 + g_8Z1Z_3Z_4 + g_{9}Z_2Z_3Z_4\nonumber\\
             &~+ g_{10}Z_1Z_2Z_3Z_4 + g_{11}X_1Z_2X_3\nonumber\\
             &~+ g_{12}Y_1Z_2Y_3 + g_{13}Y_1Z_2Y_3Z_4.
\end{align}
We see that due to the symmetry of the represented system in $H_{BK}$, qubit 2 and qubit 4 are never flipped, allowing us to reduce the dimension of the Hamiltonian to
\begin{align}
    H_{BK}^\text{reduced} = &~h_0I + h_1Z_1 + h_2Z_2 + h_3Z_1Z_2 \nonumber\\
                        &~+ h_4X_0X_1 + h_5Y_0Y_1.
\end{align}
This reduced representation requires only two qubits to simulate the hydrogen molecule. The HF initial state thus becomes $\ket{01}$. We emphasize that such a reduction of the BK Hamiltonian is not a special case for $H_2$ molecule but is connected to symmetry considerations to reduce the complexity of systems, in a scalable way.\par

The variational quantum eigensolver (VQE) is a method to compute the ground state energy of the Hamiltonian.
The total energy can be directly calculated by measuring the expectation value of each Hamiltonian term. This can be done easily by partial quantum state tomography. All the expectation values are then added up with a set of weights ($h_0$ through $h_5$). The weights are only functions of the internuclear separation ($R$) and can be computed efficiently by a classical computer. Here we use the \texttt{OpenFermion} python package to compute these weights~\cite{mcclean2020openfermion}.\par
The main task of the quantum processor is then to encode the molecular spin-orbital state into the qubits. The starting point is the HF initial state, which is believed to largely overlap with the actual ground state. In order to find the actual ground state, the initial state needs to be ``parameterized'' into an ansatz to explore a subspace of all possible states. We apply the unitary coupled cluster (UCC) theory to the parameterized ansatz state, which is widely believed to be a powerful approach and cannot be efficiently executed on a classical computer~\cite{taube2006new}. The UCC operator has a format
\begin{align}
    U_{UCC}(\vv{\theta}) = e^{\sum_n (T_n(\vv{\theta}) - T_n^{\dagger}(\vv{\theta}))},
\end{align}
with
\begin{align}
    &T_1(\vv{\theta}) = \sum_{m,i} \vv{\theta}_i^m a_m^{\dagger}a_i,\\
    &T_2(\vv{\theta}) = \sum_{m,n,i,j} \vv{\theta}_{i,j}^{m,n} a_m^{\dagger}a_n^{\dagger}a_ia_j
\end{align}
representing single-excitation and double-excitation of the electrons. The indices $i, j$ label the occupied spin-orbitals and $m, n$ are the labels of the unoccupied spin-orbitals. The vector $\vv{\theta}$ is the set of all parameters to optimize. In the case of a $H_2$ molecule, the UCC operator is transformed into a qubit operator as
\begin{align}
    U_{UCC}^{BK}(\vv{\theta}) = e^{-i\theta XY},
\end{align}
where $\theta$ is a single parameter to variationally optimize.
\section{Gate set tomography analysis}
We designed a customized gate set tomography (GST) experiment using the gate set $\{ I$, $X_{Q1}$, $Y_{Q1}$, $X_{Q2}$, $Y_{Q2}$, $\textsc{cphase}\}$, where $I$ is a $\unit[100]{ns}$ idle gate, $X_{Q1}$ ($Y_{Q1}$) and $X_{Q2}$ ($Y_{Q2}$) are single qubit $\pi/2$ gates on qubit 1 and 2 with rotation axis $\hat{x}$ ($\hat{y}$), and $\textsc{cphase}=\text{diag}(1,1,1,-1)$. A classic two-qubit GST experiment consists of a set of germs designed to amplify any error rate in the sequence when repeated,  and a set of 36 fiducials composed by the 11 elementary operations $\{\textit{null}$, $X_{Q1}$, $X_{Q1}X_{Q1}$, $X_{Q1}X_{Q1}X_{Q1}$, $Y_{Q1}$, $Y_{Q1}Y_{Q1}Y_{Q1}$, $X_{Q2}$, $X_{Q2}X_{Q2}$, $X_{Q2}X_{Q2}X_{Q2}$, $Y_{Q2}$, $Y_{Q2}Y_{Q2}Y_{Q2}\}$ required to do quantum process tomography of the germs~\cite{greenbaum2015introduction}. Note, that the $\textit{null}$ gate is the instruction for doing nothing in zero time, different from the idle gate. The germs and fiducials are then compiled into GST sequences such that each sequence consists of two fiducials interleaved by a single germ or power of germs (as illustrated in Fig.~\ref{fig:1QGST}a of the main text)~\cite{nielsen2020gate}. The GST sequences are classified by their germ powers into lengths $L=1,2,4,8,16\cdots$, where a sequence of length $n$ consists of $n$ gates plus the fiducial gates. After the execution of all sequences a maximum-likelihood estimation (MLE) is performed to estimate the process matrices of each gate in the gate set and the SPAM probabilities. We use the open source \texttt{pyGSTi} python package~\cite{nielsen_python_2019,nielsen_probing_2020} to perform the MLE, as well as to design a reduced GST experiment by eliminating redundant circuits, and to provide statistical error bars by computing all involved Hessians. The circuit optimization allows us to perform GST with a maximum sequence length $L_\text{max}=16$ using 1685 different sequences in total.
The \texttt{pyGSTi} package quantifies the Markovian-model violation of the experimental data counting the number of standard deviations exceeding their expectation values under the $\chi^2$ hypothesis~\cite{nielsen_probing_2020}. This model violation is internally translated into a more accessible goodness ratio from $0-5$ with 5 being the best~\cite{nielsen_python_2019}, where we get a 4 out of 5 rating indicating remarkably small deviations from expected results.

From the gate set tomography experiment, we have extracted the Pauli transfer matrix (PTM) $\mathcal{M}_\text{exp}$ describing each gate in our gate set $\{I$,  $X_{Q1}$, $Y_{Q1}$, $X_{Q2}$, $Y_{Q2}$, $\textsc{cphase}\}$. The PTM is isomorphically related to the conventionally used $\chi$-matrix describing a quantum process. A completely positive trace-preserving (CPTP) two-qubit PTM has 240 parameters describing the process. To get insight in the errors of the gates in the experiment, we first compute the error in the PTM given by $E=\mathcal{M}_\text{exp}\mathcal{M}_\text{ideal}^{-1}$, where we have adapted the convention to add the error after the ideal gate. The average gate fidelity is then conveniently given by
\begin{align}
F_\text{gate} = \frac{\text{Tr}(\mathcal{M}_\text{exp}^{-1}\mathcal{M}_\text{ideal})+d}{d(d+1)}.
\label{eq:averageGate}
\end{align}
It is related to the entanglement fidelity via $F_\text{ent}=\frac{d+1}{d}F_\text{gate}$~\cite{white_measuring_2007}, where $d$ is the dimension of the two-qubit Hilbert space. While the PTM $\mathcal{M}$ perfectly describes the errors, it is more intuitive to analyze the corresponding error generator $\mathcal{L}=\log(E)$ of the process~\cite{blume2021taxonomy}. The error generator relates $\mathcal{L}$ to the error PTM $E$ in a similar way as a Hamiltonian relates $H$ to a unitary operation $U=e^{-i H}$. The error generator can be separated into several blocks. A full discussion about the error generator can be found in Ref.~\cite{blume2021taxonomy}. In this work, we have used the error generator to distinguish the dynamics originating from coherent Hamiltonian errors, which can be corrected by adjusting gate parameters (see Extended Data Figure~\ref{fig:GST_optimization}, from noisy/stochastic dynamics which cannot be corrected easily. The coherent errors can be extracted by projecting $\mathcal{L}$ onto the $4\times 4$-dimensional Hamiltonian space $H$. In the Hilbert-Schmidt space, the Hamiltonian projection is given by~\cite{blume2021taxonomy}
\begin{align}
    H_{mn} = -\frac{i}{d^2}\text{Tr}\left[\left(P_m^T\otimes P_n^T\otimes \boldsymbol{1}_d - \boldsymbol{1}_d\otimes P_m\otimes P_n \right)\mathcal{L}_\text{sup}\right],
\end{align}
where $\mathcal{L}_\text{sup}$ is the error generator in Liouville superoperator form, $P_m\in{I,X,Y,Z}$ are the extended Pauli matrices with $m,n=0,1,2,3$, $\boldsymbol{1}_d$ is the $d$-dimensional Identity matrix, and $d=4$ is the dimension of the two-qubit Hilbert space. To improve the calibration of our gate set, we use the information of 7 different Hamiltonian errors ($IX$, $IY$, $XI$, $YI$, $ZI$, $IZ$ and $ZZ$).  To estimate coherent Hamiltonian errors and incoherent stochastic errors, two new metrics are considered~\cite{blume2021taxonomy}; the Jamiolkowski probability
\begin{align}
    \epsilon_J(\mathcal{L}) = -\text{Tr}(\rho_J(\mathcal{L})\ket{\Psi}\bra{\Psi})),
    \label{eq:Jam_prob}
\end{align}
which describes the amount of incoherent error in the process, and the Jamiolkowski amplitude
\begin{align}
    \theta_J(\mathcal{L}) = ||(1-\ket{\Psi}\bra{\Psi})\rho_J(\mathcal{L})\ket{\Psi}||_2,
    \label{eq:Jam_amp}
\end{align} which approximately describes the amount of coherent Hamiltonian errors (Extended Data Table.~\ref{tab:GST_results}). Here, $\rho_J(\mathcal{L}) = (\mathcal{L}\otimes\boldsymbol{1}_{d^2})[\ket{\Psi}\bra{\Psi}]$ is the Jamiolkowski state and $\ket{\Psi}$ is a maximally entangling four-qubit state which originates from the relation of quantum processes to states in a Hilbert space twice the dimension via the Choi-Jamiolkowski isomorphism~\cite{jamiolkowski_linear_1972}. For small errors, the average gate infidelity can be approximated by~\cite{blume2021taxonomy}
\begin{align}
    1-F_\text{gate}=\frac{d}{d+1}\left[\epsilon_J(\mathcal{L})+\theta_J(\mathcal{L})^2\right].
\end{align}
For a comparison of the performance of the single-qubit gates with previous experiments reporting single-qubit gate fidelities, we compute the fidelities projected to the single-qubit space from the PTMs or the error generators. In Fig.~\ref{fig:1QGST} and Extended Data Fig.~\ref{fig:Subspace_PTM}, single-qubit gate fidelities are estimated by projecting the PTMs onto corresponding subspace. Let $\mathcal{P}_{j}$ be the projector on the subspace of qubit $j$ then the fidelity is given by
\begin{align}
F_\text{sub} = \frac{\text{Tr}(\mathcal{P}_{j}\mathcal{M}_\text{exp}^{-1}\mathcal{P}_{j}\mathcal{M}_\text{ideal})+(d/2)}{(d/2)((d/2)+1)}.
\label{eq:averageGate_sub}
\end{align}
Error bars for the fidelity projected to the subspace are computed using standard error propagation of the confidence intervals of $\mathcal{M}_\text{exp}$ provided by the \texttt{pyGSTi} package.
A more optimistic estimation for the fidelities in the single-qubit subspace is given by projecting the error generators instead of the PTMs.

\section*{Methods}
\input{Measurement_setup}
\input{Gate_calibration}
\input{Theoretical_model}
\input{Fitting}
\input{Numerical_simulations}
\input{Optimal_pulse_shapes}
\input{gate_set_tomography}
\input{VQE}
\vfill

\input{Extended_figures}
\clearpage