\section{Quantum dot interferometer}
In this section, we provide more details on the model for the quantum dot (QD) interferometer, which we discussed in the main text.
\subsection{Quantum dot Hamiltonian}
First, we introduce the model of an interacting QD tunnel-coupled to superconducting (SC) leads 
described by the Hamiltonian,
\begin{equation}
    H=H_{QD}+H_{l}+H_T\,.
    \label{Eq:QDmodel}
\end{equation}

The first contribution to the Hamiltonian, $H$, describes the spin-degenerate, interacting QD and is given by,
\begin{equation}
    H_{QD}=n_g\sum_{\sigma}n_{\sigma}+Un_\uparrow n_\downarrow\,,
    \label{H_QD}
\end{equation}
where $n_\sigma=d^{\dagger}_{\sigma}d_{\sigma}$ is the number operator on the QD with $d_{\sigma}$ annihilating a QD electron with spin $\sigma=\uparrow,\downarrow$. Moreover, $n_g$ denotes the offset charge and $U$ is the on-site Coulomb interaction strength.

The second contribution to the Hamiltonian, $H$ in Eq. \eqref{Eq:QDmodel}, describes the SC leads and reads,
\begin{equation}
    H_{l}=\sum_{k\nu\sigma}\xi_{k\nu\sigma}c^{\dagger}_{k\nu\sigma}c_{k\nu\sigma}+\sum_{k\nu}\Delta_{k\nu}c^{\dagger}_{k\nu\uparrow}c^{\dagger}_{k\nu\downarrow}+{\rm H.c.}\,,
    \label{H_leads}
\end{equation}
where $c_{k\nu\sigma}$ is the annihilation operator for an electron in lead $\nu=L,R$ with spin $\sigma$ and quasi-momentum $k$. In the normal state of the leads, the dispersion of such an electron is given by $\xi_{k\nu\sigma}$. In the SC state of the leads, $\Delta_{k\nu}$ is magnitude of the SC order parameter.

Finally, the tunneling Hamiltonian that connects the QD to the SC leads is given by, 
\begin{equation}
    H_{T}^{}=\sum_{\nu\sigma}\left|t_{k\nu}\right|e^{i\varphi_\nu}\Delta_\nu c^{\dagger}_{k\nu\sigma}d_{\sigma}+{\rm H.c.}\,,
    \label{H_leads}
\end{equation}
where $t_{k\nu}$ is denotes the tunneling amplitude and $\varphi_\nu$ is the phase of the SC order parameter on SC lead $\nu$. We define the phase difference between the SC leads as $\varphi=\varphi_L-\varphi_R$.

\begin{figure}[b!]\centering
\includegraphics[width=1\linewidth] {Fig_S1}
\caption{(a) Supercurrent, $I$, in a QD JJ as a function of the SC phase difference, $\varphi$, and the offset charge, $n_{g}$. (b) Critical current asymmetry, $\Delta I_{c}$, as a function of the external flux, $\Phi$, and the offset charge, $n_{g}$. For both panels the parameters are $(U,t_{L},t_{R},T)=(3,2/3,2/3,1/100)\Delta$. 
}\label{S4}
\end{figure}

\subsection{Zero-bandwidth approximation}
There are different approaches for computing the supercurrent in a QD Josephson junction (JJ) such as perturbation theory in the tunneling amplitudes \cite{SpivakSM,Schrade2017SM,Rancic2019SM} or numerical renormalization group methods \cite{ChoiSM,Karrasch2008SM}. Here, we opt for an approach called the `zero-bandwith model'. This approximation, which has been shown to produce a good qualitative understanding of JJs with interacting QDs \cite{AllubSM,SaldanaSM}, approximates the SC leads by a single quasiparticle level. The Hamiltonian of the SC leads and the tunneling Hamiltonian thus take on the simplified form,
\begin{equation}
    \begin{split}
    H_{l}^{ZBM}&=\sum_{\nu}\Delta_\nu c^{\dagger}_{\nu\uparrow}c^{\dagger}_{\nu\downarrow}\,,
    \\
    H_{T}^{ZBM}&=\sum_{\nu\sigma}\left|t_\nu\right|e^{i\varphi_\nu}\Delta_\nu c^{\dagger}_{\nu\sigma}d_{\sigma}+{\rm H.c.}.
    \end{split}
\end{equation}
Here, zero energy corresponds to the center of the SC energy gap. 


\begin{figure}[b!]\centering
\includegraphics[width=1\linewidth] {Fig_S2}
\caption{Left panels: Supercurrent, $I$, as a function of the
SC phase difference, $\varphi$, and the offset charge, $n_g$, for 
various values of the interaction strength $U$. Right panels: Diode efficiency, $\eta$, as a function of the external flux, $\Phi$, and the offset charge, $n_g$ for the same interaction strength values, $U$, as in the left panels. System parameters are $(t_{L},t_{R},T)=(2/3,2/3,1/100)\Delta$.
}\label{S5}
\end{figure}

\subsection{Current-phase relation}
To compute the current-phase relation (CPR) of the QD JJ, we first perform a numerical exact diagonalization of the Hamiltonian, $H^{ZBM}$, given by, 
\begin{equation}
    H^{ZBM}=H_{QD}+H^{ZBM}_{l}+H^{ZBM}_T\,.
\end{equation}
From the energy levels of this Hamiltonian, denoted by $E_{j}(\varphi)$, we can compute the free energy at temperature $T$ as
\begin{equation}
    F(\varphi)=-k_B T\log\left(\sum_j e^{-E_j(\varphi)/k_BT}\right)\,.
\end{equation}
From this expression of the free energy, we obtain the CPR of the QD JJ as,
\begin{equation}
    I(\varphi)=\frac{2e}{\hbar}\frac{\partial F(\varphi)}{\partial \varphi}\,.\
    \label{current_F}
\end{equation}
The CPR for the QD JJ is shown in Fig. \ref{S4}(a) for the same parameters shown in Fig.\,3 of the main text. We see that there are two supercurrent reversal transitions as a function of the offset charge, $n_g$. Close to these transitions the CPR shows several sign changes within a period, resulting in a sizable higher harmonic content.


\subsection{Quantum dot interferometer}
In the main text, we considered a QD junction in parallel to a QPC junction. The diode efficiency, $\eta$, as a function of the offset charge, $n_{g}$, and the external flux, $\Phi$, 
reached a maximum of $\sim40\%$. 

\subsubsection{Critical current asymmetry}

We now want to complement the results of the main text with a numerical calculation of the critical current asymmetry, $\Delta I_{c}$, as a function of the the offset charge, $n_{g}$, and the external flux, $\Phi$. Our results are shown in Fig.\,\ref{S4}(b). We see that the critical current asymmetry is maximal in the vicinity of the supercurrent reversal transition. In this region, the CPR depicted in Fig.\,\ref{S4}(a) shows a strong non-sinusoidal behavior, characterized by a strong higher harmonic contribution, see Fig.\,2(a) in the main text.



\begin{figure}[b!] \centering
\includegraphics[width=1\linewidth] {Fig_S3}
\caption{(a) $n^{\text{th}}$ harmonic (for $n=1,\dots,4$) of a QD JJ as a function of the offset charge, $n_g$, when the Coulomb interaction strength on the QD vanishes, $U=0$. (b) Same as in (a) but for a finite Coulomb interaction strenght, $U=4\Delta$. The remaining parameters for both plots are $(t_{L},t_{R},T)=(2/3,2/3,1/100)\Delta$.
}\label{S6}
\end{figure}

\subsubsection{Effects of the quantum dot interaction strength}

We also want to provide further details on the on-site Coulomb interaction strength, $U$, modifies the JD effect. Our results are shown in Fig.\,\ref{S5} and Fig.\,\ref{S6}.

For the non-interacting case when $U=0$, 
we find that a sizable JD effect arises when the QD JJ is highly transmitting close to $n_{g}=0$, see Fig.\,\ref{S5}\,(a) and (b). In this regime, the CPR of the QD JJ exhibits sizable contributions from higher harmonics, as shown in Fig. \ref{S6}\,(a).

For the interacting case when $U>0$, we find that the effect of the interactions is two-fold: (1) The electron-hole symmetric point, defined by $I(\varphi,n_{g}-n^{*}_{g})=I(\varphi,-(n_{g}-n^{*}_{g}))$, shifts from $n^{*}_{g}=0$ to $n^{*}_{g}=-U/2$. (2) In the vicinity of the electron-hole symmetric point a supercurrent reversal transition emerges, see Fig.\,\ref{S6}\,(c),(e),(g). This transition amplifies the importance of higher harmonics in the CPR, as shown in Fig.\,\ref{S6}\,(a), and yields a JD effect with sizable efficiency, as shown in Fig.\,\ref{S6}\,(d),(f),(h). 




\section{Diode efficiency optimization}
\label{Sec::IdealDiode}

In this section, we provide more details on the diode efficiency optimization
discussed in the main text. 

\subsection{Conceptual idea for improving the diode efficiency}

In the main text, we argued that it is not possible to realize a $\eta=1$ SC diode with a finite number of JJs because the CPR of JJ lacks a phase-independent (`$0^{\text{th}}$ harmonic') term. We have, therefore, proposed 
a conceptual example of a CPR that asymptotically approaches $\eta=1$, 
\begin{equation}
    I(\varphi)=\sum^{m_{\text{max}}}_{m=1}\frac{{\rm cos}[m(\varphi-\pi/2)]}{\sqrt{m}}\,.
    \label{ideal_diode}
\end{equation}
To further clarify the behavior of the CPR, we represent this CPR for 
values of $m_{\text{max}}$ in Fig. \ref{S2}\,(a). We see that upon increasing $m_{\text{max}}$, the CPR develops a sharp peak at $\varphi=\pi/2$, since
\begin{equation}
    I(\varphi=\pi/2)=\sum^{m_{\text{max}}}_{m=1}\frac{1}{\sqrt{m}}\,
\end{equation}
diverges upon increasing $m_{\text{max}}$. At the same time, the CPR approaches zero upon increasing $m_{\text{max}}$ at $\varphi=3\pi/2$ because 
\begin{equation}
    I(\varphi=3\pi/2)=\sum^{m_{\text{max}}}_{m=1}\frac{(-1)^{m}}{\sqrt{m}}\,
\end{equation}
converges to zero upon increasing $m_{\text{max}}$. 

\subsection{Concatenation of interference loops}
In the main text, we showed that the concatenation of multiple interference loops 
can significantly improve the diode efficiency. In Fig.\,3(b) of the main text, we have plotted the CPRs that yield the maximum diode efficiency for an array with $N_{J}$ JJs. 
The CPR of such an array is of the form,
\begin{equation}
    I(\varphi)=\sum_{\ell=1}^{N_J}I_{\ell}(\varphi-2\pi(\Phi^{(\ell)}-\Phi^{(\ell-1)})/\Phi_0)\,,
\end{equation}
where $\Phi^{(\ell)}$ is the flux piercing through the $\ell^{\text{th}}$ interferometer loop and $\Phi^{(0)}=0$. In Table \ref{table1} we provide the parameters the best JD (maximal $\eta$) found, where we fixed the transmission and the phase of one of the junctions to improve the convergence of our Monte Carlo maximization code. 

\begin{widetext}

\begin{table}[!t]
\centering
\caption{Parameters, $\{\tau_{\ell},2\pi(\Phi^{(\ell)}-\Phi^{(\ell-1)})/\Phi_0\}$, used to generate the CPRs in Fig.\,3(b) of the main text for $N_{J}$ JJs that are concatenated in an interferometr array. Here, $\Phi^{(\ell)}$ is the external flux that pierces through the $\ell^{\text{th}}$ interferometer loop.}
\begin{tabular}[t]{|l|c|c|c|c|c|c|c|c|c|c|c|c|c|c|cccc}
\hline
$N_J$ & 2 & 3 & 4 & 5 & 6 & 7 & 8 & 9 & 10 & 11\\
\hline
1 & $1.0,0.0$ & $1.0,0.0$& $1.0,0.0$& $1.0,0.0$ & $1.0,0.0$ & $1.0,0.0$ & $1.0,0.0$ & $1.0,0.0$ & $1.0,0.0$ & $1.0,0.0$\\
2 & $0.746,4.22$ & $0.833,4.90$& $0.9,5.21$& $0.891,0.55$& $0.935,4.50$ & $0.953,5.43$ & $0.894,0.48$ & $0.935,5.75$ & $0.951,5.64$ & $0.947,5.74$\\
3 & - & $0.828,3.02$& $0.888,2.37$& $0.887,1.83$& $0.914,0.66$ & $0.943,1.55$ & $0.884,4.53$ & $0.924,1.54$ & $0.893,1.50$ & $0.935,1.46$\\
4 & - & -& $0.823,3.78$& $0.879,4.14$& $0.865,3.35$ & $0.888,2.56$ & $0.780,3.82$ & $0.845,2.49$ & $0.815,2.37$ & $0.835,4.85$\\
5 & -  & - & - &  $0.877,3.13$ & $0.814,1.81$ & $0.687,4.62$ & $0.697,2.60$ & $0.811,4.68$ & $0.799,4.58$ & $0.797,2.38$\\
6 & -  & - & - & - & $0.635,2.43$ & $0.681,3.91$ & $0.677,1.30$ & $0.700,3.35$ & $0.724,3.43$ & $0.748,4.05$\\
7 & -  & - & - & - & - & $0.652,3.62$ & $0.651,1.76$ & $0.543,3.26$ & $0.544,4.51$ & $0.615,3.02$\\
8 & -  & - & - & - & - & - & $0.592,3.10$ & $0.541,4.54$ & $0.336,3.00$ & $0.507,3.08$\\
9 & -  & - & - & - & - & - & - & $0.106,4.45$ & $0.255,3.33$ & $0.481,5.02$\\
10 & -  & - & - & - & - & - & - & - & $0.097,2.86$ & $0.356,2.73$\\
11 & -  & - & - & - & - & - & - & - & - & $0.143,2.74$\\
\hline
$\eta$ & 0.374  & 0.538 & 0.631 & 0.681 & 0.708 & 0.722 & 0.737 & 0.744 & 0.750 & 0.761\\
\hline
\end{tabular}
\label{table1}
\end{table}%
\end{widetext}


\begin{figure} \centering
\includegraphics[width=1\linewidth] {Fig_S4}
\caption{(a) Plot of the CPR given by Eq.\,\eqref{ideal_diode} that asymptotically approaches $\eta=1$. (b) Example CPRs used for the construction of the
bound, $\eta^{*}$. 
}\label{S2}
\end{figure}

\subsection{Diode efficiency bound}
In the main text, we argued that $\eta^{*}=(N_{J}-1)/N_{J}$ can provide an upper bound to the diode efficiency of $N_{J}$ JJs that are not diodes themselves. Here, we present more details on how we arrived at this bound with the help of a geometrical argument: 

Given the $N_{J}$ JJs, we assume that the the first Josephson junction (JJ) hosts a flat CPR with $I_{c}^{+}=I_{c}^{-}\equiv I_c$, black line in Fig. \ref{S2} (b). To achieve a high efficiency, the supercurrent maximum and minimum of a second JJ should align with maximum of the first one,  see Fig. \ref{S2} (b). The maximal efficiency is achieved for a critical current of the second JJ of $2I_c$. The sum of the described CPRs leads to $I_{c}^{+}=3I_{c}^{-}=3I_c$, leading to $\eta=0.5$ when $N_{J}=2$. It is possible to continue adding more peaked Josephson elements to increase the efficiency. The requirement is to align the peaks in one direction, while they should not overlap in the other one, dashed blue line in Fig. Fig. \ref{S2}. The critical currents in both directions scales with the number of JJs as
\begin{eqnarray}
    I_{c}^{+}=(2N_J-1)I_c\,,\qquad I_{c}^{-}=I_c\,,
\end{eqnarray}
leading to an upper bound for the diode efficiency
\begin{eqnarray}
    \eta^{*}=\frac{N_J-1}{N_J}\,.
\end{eqnarray}
We note that we have checked that $\eta^{*}$ also provides a bound for the CPR 
in Eq. \eqref{ideal_diode} when summing up a finite number of harmonics $m_{\rm max}$. 



\begin{figure}[!t] \centering
\includegraphics[width=1\linewidth] {Fig_S5}
\caption{Differential resistance, $\text{d}V/\text{d}i_{\text{dc}}$, versus $i_{\rm ac}$ and $i_{\rm dc}$ for $\Phi/\Phi_{0}=1/2$ and $(\tau_{1},\tau_{2})=(0.99,0.96)$ in panel (a) and $(\tau_{1},\tau_{2})=(0.99,0.99)$ in panel (b), showing the appearance
of additional, fractional Shapiro steps. 
}\label{S7}
\end{figure}


\section{Fractional Shapiro steps}
In the main text, we argued that in a supercurrent 
interferometer comprised of QPC JJs that contain higher harmonic content in their CPR, the first harmonic can be cancelled when
the transmissions of the two JJs are balanced by local gates ($\tau_{1}=\tau_{2}$) and a half-flux quantum ($\Phi=\Phi_{0}/2$) pierces through the interferometer loop.
The dominant contribution to the supercurrent is then 
due to the coherent tunneling of pairs of Cooper-pairs. 

To reveal this unusual tunneling effect of pairs of Cooper-pairs, we suggested in the main text the use of an AC drive that should give rise to a doubling of the Shapiro steps. 
While we showed this doubling of the Shapiro steps with the
help of the $i_{dc}-V$ curve in the main text, it can also be revealed by plotting the differential resistance, $\text{d}V/\text{d}i_{\text{dc}}$, versus $i_{\rm ac}$ and $i_{\rm dc}$. Here, we provide this additional information in Fig.\,\ref{S7}. When approaching the  balanced situation, the the differential resistance peaks split. In the balanced configuration, ($\tau_1=\tau_2$ and $\Phi=\Phi_{0}/2$) the peaks have doubled, showing the same height. 



Finally, we would like to point out that is is possible to engineer junctions where the leading contribution to the supercurrent arises from the tunneling of $n>2$ Cooper pairs. The criterion is to use at least $n$ JJs with tunable transmission and magnetic flux. The interferometer CPR in this case is given by
\begin{equation}
    I(\varphi)=\sum_{m=0}^{n-1}I_m(\varphi-2m\pi\Phi/n\Phi_0)\,,
\end{equation}
where $I$ is given by Eq. (4) of the main text and the transmission of every CPR has to be equal.


