# Klein tunneling through triple barrier in AB bilayer graphene

Mouhamadou Hassane Saley<EMAIL_ADDRESS>Laboratory of Theoretical
Physics, Faculty of Sciences, Chouaïb Doukkali University, PO Box 20, 24000 El
Jadida, Morocco  Abderrahim El Mouhafid<EMAIL_ADDRESS>Laboratory of
Theoretical Physics, Faculty of Sciences, Chouaïb Doukkali University, PO Box
20, 24000 El Jadida, Morocco  Ahmed Jellal<EMAIL_ADDRESS>Laboratory of
Theoretical Physics, Faculty of Sciences, Chouaïb Doukkali University, PO Box
20, 24000 El Jadida, Morocco  Canadian Quantum Research Center, 204-3002 32
Ave Vernon, BC V1T 2L7, Canada  Ahmed Siari  Laboratory of Measurement and
Control Instrumentation, Faculty of Sciences, Chouaïb Doukkali University, PO
Box 20, 24000 El Jadida, Morocco

###### Abstract

We investigate the transport properties of charge carriers in AB bilayer
graphene through a triple electrostatic barrier. We calculate the transmission
and reflection using the continuity conditions at the interfaces of the triple
barrier together with the transfer matrix method. First, we consider the case
where the energy is less than the interlayer coupling $\gamma_{1}$ and show
that, at normal incidence, transmission is completely suppressed in the gap
for a large barrier width while it appears in the gap for a small barrier
width. For energies greater than $\gamma_{1}$, we show that in the absence of
an interlayer potential difference, transmission is less than that of a single
barrier, but in its presence, transmission in the gap region is suppressed, as
opposed to a double barrier. It is found that one, two, or three gaps can be
created depending on the number of interlayer potential differences applied.
Resonance in the $T_{-}^{+}$ transmission channel is observed that is not seen
in the single and double barrier cases. Finally, we compute the conductance
and show that the number of peaks is greater than the double barrier case.

###### pacs:

72.80.Vp, 73.21.Ac, 73.23.Ad
Keywords: Bilayer graphene, AB-Stacking, triple barrier, transmission, Klein
tunneling, conductance

##  I Introduction

Graphene graphene , a gapless semiconductor with just one layer of carbon
atoms, was discovered in 2004 by Geim and Novoselov and is now one of the most
active research topics in the field of condensed matter because of its
interesting physical properties. In fact, in addition to holding the record
for the highest thermal conductivity, graphene presents high electronic
quality with mobility of the order of 250 000 $\text{cm}^{2}/\text{Vs}$, 20
times higher than that of GaAs Morozov ; Geim . With a Young’s modulus of
$\sim$ 1 TPa and an intrinsic tensile strength of $\sim$ 130 GPa Changgu , it
is one of the strongest materials discovered. From the optical point of view,
graphene is transparent, it observes only 2.3% of visible light Nair . Due to
the symmetry of the graphene crystal structure, its charge carriers, which
behave like relativistic particles, are described by the Dirac-like
Hamiltonian for massless fermions Katsnelson_2006 ; Ramezani ; FONSECA20104359
. Graphene also has unusual properties such as Klein tunneling Katsnelson_2006
; Young_2009 ; Stander  and a non-conventional quantum Hall effect
Novoselov_2005 ; Zhang ; Gusynin . All of these outstanding features
demonstrate graphene’s relevance. However, the Klein tunneling is an obstacle
in terms of possible applications. A graphene-based transistor, for example,
cannot have an on/off switching state, according to Tudorovskiy_2012 ;
Yung_2013 ; Ramezani . It is then crucial to open a gap in the energy spectrum
of graphene for its application in semiconducting devices. Fortunately, there
are several methods to achieve this goal, such as coupling to a substrate or
doping with impurities Zhou ; Jung_2015 ; Filho ; San ; Zarenia .

In the case of graphene bilayer AB (Bernal Bernal ), a system of two layers of
carbon atoms placed on top of each other, it is easy to open and control a gap
by applying an external electric field McCanne ; Falko ; Barbier . In contrast
to single layer graphene, which has a linear energy band, bilayer graphene has
a parabolic energy spectrum with four bands. Furthermore, unlike single-layer
graphene, tunneling through a potential barrier at normal incidence is
characterized by a perfect reflection Katsnelson_2006 ; Tudorovskiy_2012 ,
anti-Klein tunneling, due to confined state cloaking GU . However, in cases of
double barrier Mouhafid ; Lu_2015 ; Bai  the anti-Klein tunneling is destroyed
and transmission resonances appear as a consequence of the bound electron
states present in the well. In addition, it is shown that in the presence of a
single magnetic barrier, the cloak effect is attenuated, so the anti-tunneling
is reversed GU ; Agrawal_Garg__2012 . Indeed, several studies on the transport
properties of graphene bilayers have been conducted Katsnelson_2006 ;
Tudorovskiy_2012 ; Barbier ; Ramezani ; Snymam ; Barbier_2010 , but only for
energies less than the interlayer coupling $\gamma_{1}\simeq 0.4$ eV McCanne .
In this case ($E<\gamma_{1}$) only one mode of propagation is available, and
then there is only one transmission channel. Recently, some studies Van ;
Mouhafid ; JELLAL2015149 ; Redouani ; BENLAKHOUY2021114835 , dealing with
tunneling through electrostatic and even magnetic barriers, have also taken
into account the case of energy higher than the interlayer coupling where two
propagation modes exist. As a result, four transmission channels are
accessible.

Based on previous works, especially Mouhafid , we investigate the transport
properties of Dirac fermions in bilayer graphene through a triple barrier
structure. Using the transfer matrix approach, we determine the transmission
and reflection probabilities as well as the associated conductance. We first
evaluate the transmission probability for energy less than $\gamma_{1}$ where
we find the dependence of the gap on the barrier widths. On the other hand,
four transmission and reflection channels are available for energy greater
than $\gamma_{1}$. We show that the transmission is suppressed in the gap
region and the number of gaps created depends strongly on the number of
interlayer potential differences applied in each region. In addition, we find
resonances that are not found in the case of single Van  and double barrier
Mouhafid , especially in the $T_{-}^{+}$ transmission channel. Moreover, the
conductance displays more peaks than that in Mouhafid .

This paper is organized as follows. In section II, we present our system’s
Hamiltonian and derive the corresponding energy eigenvalues and eigenvectors.
In section III, we use the transfer matrix approach and continuity conditions
to establish the transmission and reflection probabilities as well as the
related conductance. In section IV, we present and analyze the numerical
results obtained. Finally, we present a summary of the main findings in
section V.

##  II Theory and methods

Figure 1: A schematic representation of the triple barrier potential.

We consider AB-stacked bilayer graphene in the presence of triple barrier
potential along the $x$-direction as shown in Fig. 1. The Hamiltonian which
describes the Dirac fermions in bilayer graphene of AB-stacking is obtained by
a nearest-neighbor tight-binding approximation near the $K$ point Van ;
Mouhafid , which is

$H=\begin{pmatrix}V^{+}&v_{F}\pi^{\dagger}&-v_{4}\pi^{\dagger}&v_{3}\pi\\\
v_{F}\pi&V^{+}&\gamma_{1}&-v_{4}\pi^{\dagger}\\\
-v_{4}\pi&\gamma_{1}&V^{-}&v_{F}\pi^{\dagger}\\\
v_{3}\pi^{\dagger}&-v_{4}\pi&v_{F}\pi&V^{-}\end{pmatrix}$ (1)

where $\pi=p_{x}+ip_{y}$ denotes the in-plan momenta with
$p_{x,y}=-i\hbar\partial_{x,y}$ and the Fermi velocity is $v_{F}\approx
10^{6}$ m/s. The interlayer coupling term is
$v_{3,4}=\frac{v_{F}\gamma_{3,4}}{\gamma_{0}}$, which refers to the coupling
between the layers. As shown below, the potentials $V^{+}$ and $V^{-}$ on the
first and second layers are defined in each region denoted by $j$

$V^{\pm}_{j}=\begin{cases}0,&\qquad j=\text{1, 5}\\\
U_{j}\pm\zeta\delta_{j},&\qquad j=\text{2, 3, 4}\\\ \end{cases}$ (2)

where $U_{j}$ is the barrier strength, $\delta_{j}$ is the electrostatic
potential difference between layers, and $\zeta=+1$ for the first layer and
$\zeta=-1$ for the second layer. In the following calculations, we ignore the
skew parameters $\gamma_{3}=0.315$ eV and $\gamma_{4}=0.044$ eV because they
have no effect on the band structure at high energy or transmission at low
energy McCanne ; Van . Therefore, in the $j$-th region, the Hamiltonian (1)
can be written as

$H_{j}=\begin{pmatrix}U_{j}+\delta_{j}&v_{F}\pi^{\dagger}&0&0\\\
v_{F}\pi&U_{j}+\delta_{j}&\gamma_{1}&0\\\
0&\gamma_{1}&U_{j}-\delta_{j}&v_{F}\pi^{\dagger}\\\
0&0&v_{F}\pi&U_{j}-\delta_{j}\end{pmatrix}$ (3)

in the basis of four-component spinor of the form
$\psi^{j}(x,y)=[\psi^{j}_{A_{1}},\psi^{j}_{B_{1}},\psi^{j}_{A_{2}},\psi^{j}_{B_{2}}]^{\dagger}$,
the symbol $\dagger$ denotes the transpose row vector, with 1 and 2 indexing
the two layers. Since the momentum is a conserved quantity along the
$y$-direction, then the spinor can be separated as follows

$\psi^{j}(x,y)=e^{ik_{y}y}[\phi^{j}_{A_{1}},\phi^{j}_{B_{1}},\phi^{j}_{A_{2}},\phi^{j}_{B_{2}}]^{\dagger}.$
(4)

For simplicity, let us introduce a length scale defined by $l=\frac{\hbar
v_{F}}{\gamma_{1}}\approx$1.64 nm. As a result, we can define dimensionless
quantities that are $E_{j}\equiv\frac{E_{j}}{\gamma_{1}}$,
$U_{j}\equiv\frac{U_{j}}{\gamma_{1}}$,
$\delta_{j}\equiv\frac{\delta_{j}}{\gamma_{1}}$, $x\equiv\frac{x}{l}$, and
$k_{y}\equiv lk_{y}$.

A general solution of the spinors can be determined in regions with potential
barriers ($j$ = 2, 3, 4) by solving the eigenvalue equation

$\displaystyle H_{j}\psi_{j}=E_{j}\psi_{j}$ (5)

and then we just have to set the potentials $V^{\pm}_{j}=0$ to obtain the
solutions in both regions 1 and 5. By injecting Eqs. (3-4) into Eq. (5), we
end up with four coupled differential equations. They are given by

$\displaystyle-i(\partial_{x}+k_{y})\phi^{j}_{B_{1}}$
$\displaystyle=(\varepsilon_{j}-\delta_{j})\phi^{j}_{A_{1}}$ (6a)
$\displaystyle-i(\partial_{x}-k_{y})\phi^{j}_{A_{1}}$
$\displaystyle=(\varepsilon_{j}-\delta_{j})\phi^{j}_{B_{1}}-\phi^{j}_{A_{2}}$
(6b) $\displaystyle-i(\partial_{x}+k_{y})\phi^{j}_{B_{2}}$
$\displaystyle=(\varepsilon_{j}+\delta_{j})\phi^{j}_{A_{2}}-\phi^{j}_{B_{1}}$
(6c) $\displaystyle-i(\partial_{x}-k_{y})\phi^{j}_{A_{2}}$
$\displaystyle=(\varepsilon_{j}+\delta_{j})\phi^{j}_{B_{2}}$ (6d)

where we have set $\varepsilon_{j}=E_{j}-U_{j}$, $k_{y}$ is the wave vector
along $y$-direction. By eliminating progressively the unknowns in the above
equations, we obtain a differential equation of second degree. Choosing
$\phi^{j}_{B_{1}}$, for example, yields

$\left[\partial_{x}^{2}+(k^{\tau}_{j})^{2}\right]\phi^{j}_{B_{1}}=0,$ (7)

where $\tau=\pm$ indicates the propagation modes and $k^{\tau}_{j}$ is the
corresponding longitudinal wave vector expressed as below

$k^{\tau}_{j}=\left[\varepsilon_{j}^{2}+\delta_{j}^{2}-k_{y}^{2}+\tau\sqrt{\varepsilon_{j}^{2}(1+4\delta_{j}^{2})-\delta_{j}^{2}}\right]^{1/2}$
(8)

and then the energy is given by

$\varepsilon^{\tau}_{\pm,j}=\pm\left[k^{2}+\delta_{j}^{2}+\frac{1}{2}+\tau\sqrt{k^{2}(1+4\delta_{j}^{2})+\frac{1}{4}}\right]^{1/2}$
(9)

where the wave vector is $k=[{(k^{\tau}_{j})^{2}+k_{y}^{2}}]^{\frac{1}{2}}$.
For regions 1 and 5, where there is no potential barrier, the corresponding
wave vector and energy spectrum can be derived from Eqs. (8-9) as follows

$\displaystyle
k^{\tau}_{0}=[\varepsilon_{j}^{2}+\tau\varepsilon_{j}-k^{2}_{y}]^{1/2}$ (10)
$\displaystyle
E^{\tau}_{\pm}=\pm\left[k^{2}+\frac{1}{2}+\tau\sqrt{k^{2}+\frac{1}{4}}\right]^{1/2}.$
(11)

Fig. 2 displays the energy spectrum of each region of the system presented in
Fig. 1. In the absence of an interlayer potential difference ($\delta_{j}=0$),
we observe two bands that touch at $k=0$ as well as two others that are
separated. As shown in Fig. 2(a), the result is that the upper band in the
conduction band of one region coincides with the lower band in the conduction
band of the adjacent region. When $\delta_{j}\neq 0$, the two touching bands
separate, resulting in a gap, as shown in Fig. 2(b). Hence, in the conduction
band, the upper band (dotted lines) of a region and the lower (dashed lines)
band of the next region are distinguished. If the wave vector in the first
region $k^{+}_{0}$ is real while $k^{-}_{0}$ is imaginary, then propagation is
only possible in the $k^{+}_{0}$ mode. when $E>\gamma_{1}$, $k^{-}_{0}$ is
real too, propagation is possible via both $k^{+}_{0}$ and $k^{-}_{0}$ modes.

Figure 2: (Color online) The energy as a function of wave vector $k_{y}$. (a):
Black lines represent the band’s outside barriers while blue dashed line
($U_{2}=1\gamma_{1}$), green dashed lines ($U_{3}=2\gamma_{1}$) and red dashed
line ($U_{4}=3\gamma_{1}$) show, respectively, the band’s inside regions 2, 3
and 4. (b): The same as in (a), but with
$\delta_{2}=\delta_{3}=\delta_{4}=0.15\gamma_{1}$, with the dotted lines
represent the upper band in the conduction bands in each region.

Now let us return to Eq. (7) and write its solution as a linear combination of
plane waves. This is

$\phi^{j}_{B_{1}}=c_{1}e^{ik^{+}_{j}x}+c_{2}e^{-ik^{+}_{j}x}+c_{3}e^{ik^{-}_{j}x}+c_{4}e^{-ik^{-}_{j}x}$
(12)

where $c_{n}$ denotes the normalization coefficients and $n=1,2,3,4$. One can
easily derive the other components of the spinor by injecting (12) into
(6a-6d). This process yields the following results

$\displaystyle\phi^{j}_{A_{1}}$
$\displaystyle=c_{1}f^{+}_{-}e^{ik^{+}_{j}x}-c_{2}f^{+}_{+}e^{-ik^{+}_{j}x}+c_{3}f^{-}_{-}e^{ik^{-}_{j}x}-c_{4}f^{-}_{+}e^{-ik^{-}_{j}x}$
(13) $\displaystyle\phi^{j}_{A_{2}}$
$\displaystyle=c_{1}h^{+}e^{ik^{+}_{j}x}+c_{2}h^{+}e^{-ik^{+}_{j}x}+c_{3}h^{-}e^{ik^{-}_{j}x}+c_{4}h^{-}e^{-ik^{-}_{j}x}$
(14) $\displaystyle\phi^{j}_{B_{2}}$
$\displaystyle=c_{1}h^{+}g^{+}_{+}e^{ik^{+}_{j}x}-c_{2}h^{+}g^{+}_{-}e^{-ik^{+}_{j}x}+c_{3}h^{-}g^{-}_{+}e^{ik^{-}_{j}x}-c_{4}h^{-}g^{-}_{-}e^{-ik^{-}_{j}x}$
(15)

and the quantities

$\displaystyle f^{\tau}_{\pm}=\frac{k^{\tau}_{j}\pm
ik_{y}}{\varepsilon_{j}-\delta_{j}}$ (16) $\displaystyle
h^{\tau}=\frac{(\varepsilon_{j}-\delta_{j})^{2}-(k^{\tau}_{j})^{2}-k^{2}_{y}}{\varepsilon_{j}-\delta_{j}}$
(17) $\displaystyle g^{\tau}_{\pm}=\frac{k^{\tau}_{j}\pm
ik_{y}}{\varepsilon_{j}+\delta_{j}}$ (18)

have been defined. It is more suitable to write the general solution in matrix
form in each region

$\psi^{j}(x,y)=G_{j}M_{j}(x)C_{j}e^{ik_{y}y}$ (19)

where $C_{j}$ represents the matrix form of coefficients $c_{n}$, and the
matrix $G_{j}$ is given by

$G_{j}=\begin{pmatrix}f^{+}_{-}&-f^{+}_{+}&f^{-}_{-}&-f^{-}_{+}\\\ 1&1&1&1\\\
h^{+}&h^{+}&h^{-}&h^{-}\\\
h^{+}g^{+}_{+}&-h^{+}g^{+}_{-}&h^{-}g^{-}_{+}&-h^{-}g^{-}_{-}\\\
\end{pmatrix}$ (20)

as well as $M_{j}(x)$

$M_{j}(x)=\begin{pmatrix}e^{ik^{+}_{j}x}&0&0&0\\\ 0&e^{-ik^{+}_{j}x}&0&0\\\
0&0&e^{ik^{-}_{j}x}&0\\\ 0&0&0&e^{-ik^{-}_{j}x}\\\ \end{pmatrix}.$ (21)

As is obvious, solutions for regions outside the barriers are obtained by
imposing $U_{j}=\delta_{j}=0$. Accordingly, in region 1, we have

$\displaystyle\phi^{\text{1}}_{A_{1}}$
$\displaystyle=\delta_{\tau,1}f^{+}_{-}e^{ik^{+}_{0}x}-r^{\tau}_{+}f^{+}_{+}e^{-ik^{+}_{0}x}+\delta_{\tau,-1}f^{-}_{-}e^{ik^{-}_{0}x}-r^{\tau}_{-}f^{-}_{+}e^{-ik^{-}_{0}x}$
(22) $\displaystyle\phi^{\text{1}}_{B_{1}}$
$\displaystyle=\delta_{\tau,1}e^{ik^{+}_{0}x}+r^{\tau}_{+}e^{-ik^{+}_{0}x}+\delta_{\tau,-1}e^{ik^{-}_{0}x}+r^{\tau}_{-}e^{-ik^{-}_{0}x}$
(23) $\displaystyle\phi^{\text{1}}_{A_{2}}$
$\displaystyle=\delta_{\tau,1}h^{+}e^{ik^{+}_{0}x}+r^{\tau}_{+}h^{+}e^{-ik^{+}_{0}x}+\delta_{\tau,-1}h^{-}e^{ik^{-}_{0}x}+r^{\tau}_{-}h^{-}e^{-ik^{-}_{0}x}$
(24) $\displaystyle\phi^{\text{1}}_{B_{2}}$
$\displaystyle=\delta_{\tau,1}h^{+}g^{+}_{+}e^{ik^{+}_{0}x}-r^{\tau}_{+}h^{+}g^{+}_{-}e^{-ik^{+}_{0}x}+\delta_{\tau,-1}h^{-}g^{-}_{+}e^{ik^{-}_{0}x}-r^{\tau}_{-}h^{-}g^{-}_{-}e^{-ik^{-}_{0}x}$
(25)

and we write in region 5

$\displaystyle\phi^{\text{5}}_{A_{1}}$
$\displaystyle=t^{\tau}_{+}f^{+}_{-}e^{ik^{+}_{0}x}+t^{\tau}_{-}f^{-}_{-}e^{ik^{-}_{0}x}$
(26) $\displaystyle\phi^{\text{5}}_{B_{1}}$
$\displaystyle=t^{\tau}_{+}e^{ik^{+}_{0}x}+t^{\tau}_{-}e^{ik^{-}_{0}x}$ (27)
$\displaystyle\phi^{\text{5}}_{A_{2}}$
$\displaystyle=t^{\tau}_{+}h^{+}e^{ik^{+}_{0}x}+t^{\tau}_{-}h^{-}e^{ik^{-}_{0}x}$
(28) $\displaystyle\phi^{\text{5}}_{B_{2}}$
$\displaystyle=t^{\tau}_{+}h^{+}g^{+}_{+}e^{ik^{+}_{0}x}+t^{\tau}_{-}h^{-}g^{-}_{+}e^{ik^{-}_{0}x}.$
(29)

The next section will be devoted to the calculation of the transmission and
reflection probabilities as well as the conductance. Given that there is no
potential in the incident ($j$=1) and the transmitted ($j$= 5) regions, we
will use the following relation

$G_{\text{1}}M_{\text{1}}(x)=G_{\text{5}}M_{\text{5}}(x).$ (30)

##  III Transport properties

In order to establish the transmission and reflection probabilities, it is
necessary to impose the continuity of the spinors at the different interfaces
of the triple barrier structure. In addition, by using the matrix transfer
method, it becomes easy to connect the coefficients of the incident region to
those of the transmitted one. These coefficients are given by

$C_{\text{1}}^{\tau}=\begin{pmatrix}\delta_{\tau,1}\\\ r_{+}^{\tau}\\\
\delta_{\tau,-1}\\\ r_{-}^{\tau}\end{pmatrix},\quad
C_{\text{5}}^{\tau}=\begin{pmatrix}t_{+}^{\tau}\\\ 0\\\ t_{-}^{\tau}\\\
0\end{pmatrix}$ (31)

where the Kronecker delta symbol is $\delta_{\tau,\pm 1}$. Using the
continuity at interfaces $x=0,a,b,c$, the following results are obtained

$\displaystyle
G_{\text{1}}M_{\text{1}}(0)C_{1}=G_{\text{2}}M_{\text{2}}(0)C_{\text{2}}$ (32)
$\displaystyle
G_{\text{2}}M_{\text{2}}(a)C_{\text{2}}=G_{\text{3}}M_{\text{3}}(a)C_{\text{3}}$
(33) $\displaystyle
G_{\text{3}}M_{\text{3}}(b)C_{\text{3}}=G_{\text{4}}M_{\text{4}}(b)C_{\text{4}}$
(34) $\displaystyle
G_{\text{4}}M_{\text{4}}(c)C_{\text{4}}=G_{\text{5}}M_{\text{5}}(c)C_{\text{5}}.$
(35)

As a result, using the transfer matrix $N$, the coefficients
$C_{\text{1}}^{\tau}$ and $C_{\text{5}}^{\tau}$ can be connected as follows

$C_{\text{1}}^{\tau}=NC_{\text{5}}^{\tau}.$ (36)

where $N$ is defined as

$N=\prod_{j=\text{1}}^{\text{4}}M_{j}^{-1}(x_{j})G_{j}^{-1}G_{j+1}M_{j+1}(x_{j}).$
(37)

Therefore, we can find the transmission and reflection coefficients from

$\begin{pmatrix}t_{+}^{\tau}\\\ r_{+}^{\tau}\\\ t_{-}^{\tau}\\\
r_{-}^{\tau}\end{pmatrix}=\begin{pmatrix}N_{11}&0&N_{13}&0\\\
N_{21}&-1&N_{23}&0\\\ N_{31}&0&N_{33}&0\\\
N_{41}&0&N_{43}&-1\end{pmatrix}^{-1}\begin{pmatrix}\delta_{\tau,1}\\\ 0\\\
\delta_{\tau,-1}\\\ 0\end{pmatrix}$ (38)

where $N_{ij}$ are the matrix elements of $N$ given in Eq. (37). Finally, we
end up with the transmission and reflection coefficients

$\displaystyle
t_{+}^{\tau}=\frac{N_{33}\delta_{\tau,1}-N_{13}\delta_{\tau,-1}}{N_{11}N_{33}-N_{13}N_{31}},\quad
t_{-}^{\tau}=\frac{N_{11}\delta_{\tau,-1}-N_{31}\delta_{\tau,1}}{N_{11}N_{33}-N_{13}N_{31}}$
(39) $\displaystyle r_{+}^{\tau}=N_{21}t_{+}^{\tau}+N_{23}t_{-}^{\tau},\quad
r_{-}^{\tau}=N_{41}t_{+}^{\tau}+N_{43}t_{-}^{\tau}.$ (40)

To obtain the transmission and reflection probabilities, we use the current
density. It is important to take into account the change in the wave velocity
since the electrons can be scattered into different modes of propagation. The
current density is given by

$\textbf{j}=v_{F}{\Psi}^{\dagger}\vec{\alpha}\Psi$ (41)

where $\vec{\alpha}$ is a $4\times 4$ matrix with two Pauli matrices
$\sigma_{x}$ on its diagonal and the other elements are zero. We can now
derive the incident $\textbf{j}_{\text{inc}}$, reflected
$\textbf{j}_{\text{ref}}$ and transmitted $\textbf{j}_{\text{tra}}$ current
densities from Eq. (41). As a result, the transmission and reflection
probabilities are given by

$\displaystyle
T_{\pm}^{\tau}=\frac{\left|\textbf{j}_{\text{tra}}\right|}{\left|\textbf{j}_{\text{inc}}\right|}=\frac{k_{0}^{\pm}}{k_{0}^{\tau}}\left|t_{\pm}^{\tau}\right|^{2}$
(42) $\displaystyle
R_{\pm}^{\tau}=\frac{\left|\textbf{j}_{\text{ref}}\right|}{\left|\textbf{j}_{\text{inc}}\right|}=\frac{k_{0}^{\pm}}{k_{0}^{\tau}}\left|r_{\pm}^{\tau}\right|^{2}.$
(43)

The conservation of the total probability results in

$\sum_{i,j}\left(T_{i}^{j}+R_{i}^{j}\right)=1$ (44)

where $i,j=\pm$. The index $i$ points to the arriving mode, while the index
$j$ points to the exiting mode. Consider the channel $k^{+}$, for example, and
we get $T_{+}^{+}+T_{+}^{-}+R_{+}^{+}+R_{+}^{-}=1.$

Since the transmission probabilities have already been determined, one can
calculate the conductance according to the formula of Landauer-Büttiker
BLANTER20001

$\mathrm{G}(E)=G_{0}\frac{L_{y}}{2\pi}\int_{-\infty}^{+\infty}\mathrm{d}k_{y}\sum_{\tau,n=\pm}T_{n}^{\tau}\left(E,k_{y}\right)$
(45)

with $L_{y}$ the width of the sample in the $y$-direction and
$G_{0}=4\frac{\mathrm{e}^{2}}{h}$, the factor 4 is the result of the valley
and spin degeneracy in graphene.

##  IV Results and discussions

###  IV.1 Two band Tunneling

Figure 3: (Color online): Transmission as a function of energy $E$, at normal
incidence ($k_{y}=0$), for barrier widths $b_{1}=b_{2}=b_{3}=15$ nm (blue
line), $25$ nm (red line), $60$ (green line) and $U_{2}=0.05\gamma_{1}$,
$U_{3}=0.06\gamma_{1}$, $U_{4}=0.07\gamma_{1}$. (a):
$\delta_{2}=\delta_{3}=\delta_{4}=0$. (b):
$\delta_{2}=\delta_{3}=\delta_{4}=0.03\gamma_{1}$. (c):
$\delta_{2}=\delta_{4}=0.03\gamma_{1}$ and $\delta_{3}=0$. (d):
$\delta_{2}=0.03\gamma_{1}$ and $\delta_{3}=\delta_{4}=0$.

As mentioned in section II, propagation is via a single mode in the low energy
regime ($E<\gamma_{1}$), so there is only one transmission channel. In Fig. 3,
we show the transmission as a function of energy $E$ at normal incidence
($k_{y}=0$) for various barrier widths: $b_{1}=b_{2}=b_{3}=15$ nm (blue line),
$25$ nm (red line), and $60$ nm (green line). Fig. 3(a) shows the transmission
for barriers of different heights $U_{2}=0.05\gamma_{1}$,
$U_{3}=0.06\gamma_{1}$, $U_{4}=0.07\gamma_{1}$ and
$\delta_{2}=\delta_{3}=\delta_{4}=0$. The transmission behaves analogously to
the one obtained for a single barrier ($d=100$ nm) in Van , even for small
barrier widths (blue and red lines). In fact, in the first energy range,
$E<U_{2}$, we observe the manifestation of anti-Klein tunneling, which means
that the transmission is zero. This is due to the chiral nature of charge
carriers in graphene Katsnelson_2006 . However, for energy greater than
$U_{2}$ transmission resonances appear. When an interlayer bias is introduced,
the transmission probability resembles that of a single barrier ($d=100$ nm)
in gapped bilayer graphene (see BENLAKHOUY2021114835 ), but without resonance
peaks in the energy range $E<U_{2}-\delta$. In addition, at the edge of the
induced gap, transmission appears for non-large barrier widths, as observed in
the case of double barrier Mouhafid . Nevertheless, when we significantly
increase the barrier width (green line), no transmission appears in the gap
region. We conclude that the gap size strongly depends on the barrier widths.
In Fig. 3(c), we set all parameters except $\delta_{3}$ to 0, as in Fig. 3(b).
For a large barrier width (green line), it is evident that the gap region
stays the same even though $\delta_{3}$ is zero, while for a small value, the
gap region narrows in favor of transmission. Once again, we see the gap can be
adjusted by the barrier widths. In Fig. 3 (d), with
$\delta_{2}=0.03\gamma_{1}$ and $\delta_{3}=\delta_{4}=0$, the transmission
presents a similar aspect to that obtained in Fig. 3(a). The bias introduced
in region 2 of the triple barrier system cannot open a gap on its own, even
for large values of the barrier widths.

###  IV.2 Four band Tunneling

Figure 4: (Color online): Density plot of transmission and reflection as a
function of energy $E$ and the wave vector $k_{y}$ for barrier heights
$U_{2}=1\gamma_{1}$, $U_{3}=2\gamma_{1}$, $U_{4}=3\gamma_{1}$ and
$\delta_{2}=\delta_{3}=\delta_{4}=0$. The barrier widths are set to be the
same, $b_{1}=b_{2}=b_{3}=20$ nm. The dashed black lines represent the band’s
outside barriers while the dashed white lines show the band’s inside barriers.

For higher energy $(E>\gamma_{1})$, the availability of two propagation modes
gives rise to four transmission and four reflection channels, which we plot as
a function of energy $E$ and transverse wave vector $k_{y}$ in Fig. 4. We
fixed the barrier heights at $U_{2}=\gamma_{1}$, $U_{3}=2\gamma_{1}$, and
$U_{4}=3\gamma_{1}$, as well as the barrier width at $b_{1}=b_{2}=b_{4}=20$ nm
and the interlayer potential difference at $\delta_{3}=\delta_{4}=0$. As shown
in Fig. 4, there are different energy regions delimited by dashed lines, which
represent the propagating modes inside and outside the barriers. The cloak
effect occurs in the region of energy beneath $U_{2}$ at nearly normal
incidence where the two modes $k^{+}$ outside and $k^{-}$ inside are
decoupled, resulting in no scattering between them. In the same region, for
non-normal incidence, the two modes $k^{+}$ outside and $k^{-}$ inside are
coupled, then transmission appears. However, unlike the case of unbiased
single barrier Van , no resonance occurs. For nearly normal incidence, the
cloak effect appears in the energy regions $U_{2}<E<U_{3}$ and
$U_{3}<E<U_{4}$. The asymmetry of the triple barrier system breaks the
equivalence observed between $T_{-}^{+}$ and $T_{+}^{-}$ in the case of
unbiased single barrier Van . It is clear that the cloak effect occurs
exclusively for nearly normal incidence in the $T_{-}^{+}$ channel and for
energy regions $U_{2}<E<U_{3}$ and $U_{3}<E<U_{4}$, as few transmission
resonances appear. On the other hand, in the $T_{+}^{-}$ channel, the cloaking
persists even for non-normal incidence so that there is no transmission in the
same energy regions. However, the reflection channels $R_{+}^{-}$ and
$R_{-}^{+}$ remain equivalent because the incident electrons return to the
electron state Van ; Mouhafid . Transmission is suppressed in the $T_{-}^{-}$
channel in the two energy ranges of $U_{2}<E<U_{3}$ and $U_{3}<E<U_{4}$, as in
single and double barrier cases Van ; BENLAKHOUY2021114835 ; Mouhafid .

Figure 5: (Color online) The same as in Fig. 4, but now for
$\delta_{2}=\delta_{3}=\delta_{4}=0.15$. The dotted black lines represent the
band’s outside barriers while the dashed and dotted white lines show the
band’s inside barriers.

To assess the impact of the interlayer potential difference, we plot the
transmission and reflection as a function of energy $E$ and wave vector
$k_{y}$ in Fig. 5 using the same parameters as in Fig. 4, but with
$\delta_{2}=\delta_{3}=\delta_{4}=0.15\gamma_{1}$. In the cases of double
barriers, without Mouhafid  and with a magnetic field, Redouani ,
transmissions have occurred in the gap as a result of the bound states in the
well. Because there is no well in the triple barrier system, transmission is
completely suppressed in the gap regions
$U_{j}-\delta_{j}<E<U_{j}+\delta_{j}$, $j$=2, 3, 4. However, despite the
absence of a well, for small barrier widths, transmission occurs. We notice
that the more we increase the width of the barriers, the more the transmission
is suppressed in the gap regions. The cloak effect is slightly modified in the
$T_{+}^{+}$ channel and for the energy range $E<U_{2}$, and thin transmission
resonances appear, but they are smaller than those observed in Van . The
$T_{-}^{+}$ channel exhibits more transmission resonances than the result
shown in Fig. 4. In the presence of bias, no resonance appears in the
$T_{-}^{+}$ channel in the case of single Van  and double barriers Mouhafid ,
in contrast to our case. In energy regions $U_{2}<E<U_{3}$ and
$U_{3}<E<U_{4}$, however, there is no transmission in $T_{+}^{-}$ and
$T_{-}^{-}$ channels.

Figure 6: (Color online) Density plot of transmission $T_{+}^{+}$ as a
function of energy $E$ and the wave vector $k_{y}$ for barrier heights
$U_{2}=1\gamma_{1}$, $U_{3}=2\gamma_{1}$, $U_{4}=3\gamma_{1}$. (a):
$\delta_{2}=0.15\gamma_{1}$, $\delta_{3}=\delta_{4}=0.$ (b):
$\delta_{3}=0.15\gamma_{1}$, $\delta_{2}=\delta_{4}=0.$ (c):
$\delta_{4}=0.15\gamma_{1}$, $\delta_{2}=\delta_{3}=0.$ (d): $\delta_{3}=0$,
$\delta_{2}=\delta_{4}=0.15\gamma_{1}.$ The barrier widths are set to be the
same, $b_{1}=b_{2}=b_{3}=20$ nm. The dashed black lines represent the band’s
outside barriers while the dashed white lines show the band’s inside barriers.

The transmission probability $T_{+}^{+}$ is presented in Fig. 6 with the same
parameters as in Fig. 5 but for variable interlayer potential in each barrier.
We chose $\delta_{2}=0.15\gamma_{1}$ and $\delta_{3}=\delta_{4}=0.$ in Fig.
6(a), $\delta_{3}=0.15\gamma_{1}$ and $\delta_{2}=\delta_{4}=0.$ in Fig. 6(b),
$\delta_{4}=0.15\gamma_{1}$ and $\delta_{2}=\delta_{3}=0.$ in Fig. 6(c),
finally $\delta_{3}=0$ and $\delta_{2}=\delta_{4}=0.15\gamma_{1}.$ in Fig.
6(d). Introducing bias in one of the triple barrier regions, in contrast to
the low energy regime ($E<\gamma_{1}$), opens a gap. Furthermore, we see that
the cloak effect is more adjusted, so the number of thin resonances in the
energy zone under $U_{2}$ has increased in comparison to the $T_{+}^{+}$
channel in Fig. 5. It is also worth noting that Figs. 6(a,d) have fewer
resonances than Figs. 6(b,c), which is due to the gap induced in the region
$U_{2}-\delta_{2}<E<U_{2}+\delta_{2}$.

###  IV.3 Conductance

Figure 7: (Color online) Conductance of the triple barrier structure as a
function of energy $E$ for $U_{2}=1\gamma_{1}$, $U_{3}=2\gamma_{1}$,
$U_{4}=3\gamma_{1}$ and $b_{1}=b_{2}=b_{3}=20$ nm. (a):
$\delta_{2}=\delta_{3}=\delta_{4}=0$. (b):
$\delta_{2}=\delta_{3}=\delta_{4}=0.15\gamma_{1}$. (c): variable interlayer
potential difference.

In Fig. 7, we investigate the energy dependence of the conductance for barrier
heights $U_{2}=1\gamma_{1}$, $U_{3}=2\gamma_{1}$, $U_{4}=3\gamma_{1}$ and
barrier widths $b_{1}=b_{2}=b_{3}=20$ nm. Fig. 7(a) shows the conductance in
absence of interlayer potential difference, i.e.
$\delta_{2}=\delta_{3}=\delta_{4}=0$. The resonances observed in the
transmission probabilities (Fig. 4) result in peaks in the triple barrier
conductance. For energies less than $U_{2}$, the $T_{+}^{+}$ channel is the
only one contributing to total conductance, i.e. $G_{Tot}=G_{+}^{+}$. In the
energy region smaller than $U_{2}$, it is clear that, as a consequence of the
contribution of $T_{-}^{+}$ shoulders appear in the total conductance, which
increased. This is because resonances in transmission $T_{-}^{+}$ exist in
this energy region. Also because of the cloaking, the contributions of
$T_{+}^{-}$ and $T_{-}^{-}$ remain zero. However, for energies larger than
$U_{4}$, $G_{+}^{-}$ starts conducting while $G_{-}^{-}$ stays close to zero.
The contributions of $T_{-}^{+}$ and $T_{+}^{-}$ are equal in the case of
simple barrier Van ; BENLAKHOUY2021114835 , implying that
$G_{-}^{+}=G_{+}^{-}$. In our case, the asymmetric aspect of the triple
barrier breaks this equality, i.e. $G_{-}^{+}\neq G_{+}^{-}$. The conductance
for non-null values of interlayer potential difference,
$\delta_{2}=\delta_{3}=\delta_{4}=0.15\gamma_{1}$, is plotted in Fig. 7(b). We
see suppression of conductance in the gap regions
$U_{j}-\delta_{j}<E<U_{j}+\delta_{j}$ ($j$ = 2, 3, 4), resulting in a
reduction in the number of shoulders. Fig. 7(c) shows the conductances
$G_{+}^{+}$ associated with transmission $T_{+}^{+}$ plotted in Fig. 6. When
$\delta_{2}=0.15\gamma_{1}$ (black and blue dashed lines) and
$E<U_{2}+\delta_{2}$ are used, $G_{+}^{+}$ has fewer peaks than when
$\delta_{2}=0$ (green dashed and red dashed lines). This is a consequence of
the gap opened by $\delta_{2}=0.15\gamma_{1}$ around $E=U_{2}$. All of the
conductances $G_{+}^{+}$ are roughly equivalent in the energy range
$U_{2}+\delta_{2}$. It is worth noting that the conductance nears zero around
$E=U_{3}$, even when $\delta_{3}=0$. To sum up, the main differences between
the different cases are found around $E=U_{2}$ and $E=U_{4}$.

##  V Conclusion

We have investigated the transport properties of AB bilayer graphene through a
triple barrier system. Starting from the four band Hamiltonian, we have
determined the energy spectrum and the eigenvectors of each region.
Thereafter, based on the transfer matrix method and using the continuity
conditions at the different interfaces, we have calculated the transmission
and reflection probabilities. When the barrier width is reduced in the low
energy regime ($E<\gamma_{1}$), at normal incidence ($k_{y}=0$), and in the
presence of bias, it was seen that the number of transmission’s peaks increase
in the gap region. We have also observed that the application of bias in only
one region is not enough to create a gap, even when the barrier width is
increased. However, the application of bias in two regions resulted in an
opening gap for high barrier widths.

Subsequently, we have evaluated the transmission for an energy regime greater
than $\gamma_{1}$. As a result, we have shown the appearance of extra
resonances, especially in the $T_{-}^{+}$ transmission channel. Consequently,
we have reached the following conclusions. Firstly, the application of the
interlayer potential difference controls the number of gaps created. Secondly,
the number of peaks in the conduction is increased compared to the case of the
double barrier Mouhafid  because of the resonances present in $T_{-}^{+}$.

## References

  * (1) K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, Y. Zhang, S. V. Dubonos, I. V. Grigorieva, and A. A. Firsov, Science 306, 666 (2004).
  * (2) S. V. Morozov, K. S. Novoselov, M.,I. Katsnelson, F. Schedin, D. C. Elias, J. A. Jaszczak, A. K. Geim, Phys. Rev. Lett. 100, 016602 (2008).
  * (3) A. K. Geim, Science 324, 1530 (2009).
  * (4) C. Lee, X. Wei, J. W. Kysar, and J. Hone, Science 321, 385 (2008).
  * (5) R. R. Nair, P. Blake, A. N. Grigorenko, K. S. Novoselov, T. J. Booth, T. Stauber, N. M. R. Peres, and A. K. Geim, Science 320, 1308 (2008).
  * (6) M. I. Katsnelson, K. S. Novoselov, and A. K. Geim, Nat. Phys. 2, 620 (2006).
  * (7) M. Ramezani Masir, P. Vasilopoulos, and F. M. Peeters, Phys. Rev. B 79, 035409 (2009).
  * (8) J. M. Fonseca, W. A. Moura-Melo, and A. R. Pereira, Phys. Lett. A 374, 4359 (2010).
  * (9) A. F. Young and P. Kim, Nat. Phys. 5, 222 (2009).
  * (10) N. Stander, B. Huard, and D. Goldhaber-Gordon, Phys. Rev. Lett. 102, 026807 (2009).
  * (11) K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, M. I. Katsnelson, I. V. Grigorieva, S. V. Dubonos, and A. A. Firsov, Nature 438, 197 (2005).
  * (12) Y. B. Zhang, Y. W. Tan, H. L. Stormer, and P. Kim, Nature 438, 201 (2005).
  * (13) V. P. Gusynin and S. G. Sharapov, Phys. Rev. Lett. 95, 146801 (2005).
  * (14) T. Tudorovskiy, K. J. A. Reijnders, and M. I. Katsnelson, Phys. Scr. T 146, 014010 (2012).
  * (15) K.C. Yung, W.M. Wu, M.P. Pierpoint, and F.V. Kusmartsev, Contemp. Phys. 54, 233 (2013).
  * (16) S. Zhou, G. Gweon, A. Fedorov, P. First, W. Heer, D. H. Lee, F. Guinea, A. C. Neto, and A. Lanzara, Nat. Mater. 6, 770 (2007).
  * (17) J. Jung, A. M. DaSilva, A. H. MacDonald, and S. Adam, Nat. Commun. 6, 6308 (2015).
  * (18) R. N. Costa Filho, G. A. Farias, and F. M. Peeters, Phys. Rev. B 76, 193409 (2007).
  * (19) P. San-Jose, A. Gutiérrez-Rubio, M. Sturla, and F. Guinea, Phys. Rev. B 90, 075428 (2014).
  * (20) M. Zarenia, O. Leenaerts, B. Partoens, and F. M. Peeters, Phys. Rev. B 86, 085451 (2012).
  * (21) J. D. Bernal, Proc. R. Soc. A 106, 749 (1924).
  * (22) E. McCann, D.S.L. Abergel, and V.I. Fal’ko, Solid State Commun. 143, 110 (2007).
  * (23) E. McCann, D. S. L. Abergel, and V. I. Fal’ko, Eur. Phys. J. Spec. Top. 148, 91 (2007)
  * (24) M. Barbier, P. Vasilopoulos, F. M. Peeters, and J. Milton Pereira, Phys. Rev. B 79, 155402 (2009).
  * (25) N. Gu, M. Rudner, and L. Levitov, Phys. Rev. Lett. 107, 156603 (2011).
  * (26) H. M. Abdullah, E. Mouhafid, H. Bahlouli, and A. Jellal, Mater. Res. Express 4, 025009 (2017).
  * (27) W. Lu, W. Li, C. Xu, and C. Ye, J. Phys. D 48, 285102 (2015).
  * (28) C. Bai and X. Zhang, Phys. Rev. B 76, 075430 (2007).
  * (29) N. Agrawal, S. Grover, S. Ghosh, and M. Sharma, J. Condens. Matter Phys. 24, 175003 (2012).
  * (30) I. Snyman and C. Beenakker, Phys. Rev. B 75, 045322 (2007).
  * (31) M. Barbier, P. Vasilopoulos, and F. M. Peeters, Phys. Rev. B 82, 235408 (2010).
  * (32) B. Van Duppen and F. M. Peeters, Phys. Rev. B 87, 205427 (2013).
  * (33) A. Jellal, I. Redouani, and H. Bahlouli, Physica E 72, 149 (2015).
  * (34) I. Redouani, A. Jellal, and H. Bahlouli, J. Low Tem. Phys. 181, 197 (2015).
  * (35) N. Benlakhouy, A. El Mouhafid, and A. Jellal, Physica E 134, 114835 (2021).
  * (36) Ya. M. Blanter and M. B$\ddot{\text{u}}$ttiker, Phys. Rep. 336, 1 (2000).


††thanks: These authors contributed equally.††thanks: These authors
contributed equally.

# All rf-based tuning algorithm for quantum devices using machine learning

Barnaby van Straaten  Department of Materials, University of Oxford, Oxford
OX1 3PH, United Kingdom  Federico Fedele  Department of Materials, University
of Oxford, Oxford OX1 3PH, United Kingdom  Florian Vigneau  Department of
Materials, University of Oxford, Oxford OX1 3PH, United Kingdom  Joseph Hickie
Department of Materials, University of Oxford, Oxford OX1 3PH, United Kingdom
Daniel Jirovec  Institute of Science and Technology Austria, Am Campus 1, 3400
Klosterneuburg, Austria  Andrea Ballabio  L-NESS, Physics Department,
Politecnico di Milano, via Anzani 42, 22100, Como, Italy  Daniel Chrastina
L-NESS, Physics Department, Politecnico di Milano, via Anzani 42, 22100, Como,
Italy  Giovanni Isella  L-NESS, Physics Department, Politecnico di Milano, via
Anzani 42, 22100, Como, Italy  Georgios Katsaros  Institute of Science and
Technology Austria, Am Campus 1, 3400 Klosterneuburg, Austria  Natalia Ares
<EMAIL_ADDRESS>Department of Engineering Science, University of
Oxford, Oxford OX1 3PJ, United Kingdom

###### Abstract

Radio-frequency measurements could satisfy DiVincenzo’s readout criterion in
future large-scale solid-state quantum processors, as they allow for high
bandwidths and frequency multiplexing. However, the scalability potential of
this readout technique can only be leveraged if quantum device tuning is
performed using exclusively radio-frequency measurements i.e. without
resorting to current measurements. We demonstrate an algorithm that
automatically tunes double quantum dots using only radio-frequency
reflectometry. Exploiting the high bandwidth of radio-frequency measurements,
the tuning was completed within a few minutes without prior knowledge about
the device architecture. Our results show that it is possible to eliminate the
need for transport measurements for quantum dot tuning, paving the way for
more scalable device architectures.

##  I Introduction

Radio-frequency (rf) reflectometry allows for high-bandwidth measurements of
quantum devices Vigneau  _et al._ (2022). Due to their potential for
scalability, rf techniques play an increasingly important role in developing
quantum device circuits. A recent achievement has been the high-fidelity
single-shot readout of spins Hogg  _et al._ (2022); Madzik  _et al._ (2022).
However, until now, the much slower measurement of current through the device
was still necessary to automatically tune quantum devices to their operating
condition Moon  _et al._ (2020); Baart  _et al._ (2016); Mills  _et al._
(2019); Zwolak  _et al._ (2020). This was a limiting factor for scalability
because such measurements are incompatible with scalable device architectures
Cai  _et al._ (2019); Li  _et al._ (2018); Veldhorst  _et al._ (2017);
Gonzalez-Zalba  _et al._ (2021) and are too slow for error correction Laucht
_et al._ (2021).

We demonstrate an algorithm that tunes a quantum device using rf measurements
exclusively, eliminating the need for current measurements through the device.
This opens up the possibility of new, more scalable device architectures
designed only to use rf measurements. The algorithm aimed to define double
quantum dots, a cornerstone of the solid-state qubit effort. These devices can
encode either a singlet-triple qubit or a pair of Loss-Divincenzo spin qubits.
Chatterjee  _et al._ (2021a); Burkard  _et al._ (2021).

Our algorithm exploits the bandwidth of rf measurements to acquire high-
resolution two-dimensional charge stability diagrams within milliseconds
Stehlik  _et al._ (2015); Schupp  _et al._ (2020). The device’s gate voltage
space is efficiently explored using Gaussian processes, principal component
analysis to perform blind signal separation Jolliffe and Cadima (2016); Pal
_et al._ (2013); com (2010), and Kolmogorov-Smirnov statistical tests Hodges
(1958), along with our fast scans. Quantum dot features are identified using a
score function based on the Fourier transform of the charge stability diagram.

We demonstrate our tuning algorithm in a hole-based quantum dot array hosted
in a Ge/SiGe heterostructure. Ge hole spin qubits are currently enjoying
intense research interest owing to their ease of operation and compatibility
with existing Si technology Scappucci  _et al._ (2021). From 2018, and within
just four years, a Loss-DiVincenzo qubit Watzinger  _et al._ (2018), a single-
triplet qubit Jirovec  _et al._ (2021), two-qubit devices Hendrickx  _et al._
(2020) and a four-qubit Ge quantum processor Hendrickx  _et al._ (2021) have
been realised, demonstrating the potential of Ge for quantum information.
Furthermore, Ge has also hosted record-breaking ultra-fast control
frequencies, with Rabi frequencies exceeding $540\text{\,}\mathrm{MHz}$ Wang
_et al._ (2022).

Our algorithm contributes significantly towards scalable readout, an essential
step in realising the scalability potential of Si-compatible quantum circuits.

Figure 1: (a) A scanning electron micrograph of the Ge/SiGe quantum dot array,
with an L-matching circuit connected to the device ohmic. Bias-tees are used
on gates shaded in red to apply both ac and dc voltages. Devices A and B
demonstrate how two different double dots can be defined in the array. (b) A
diagram showing the waveforms sent to gates $V_{2}$ ($V_{4}$) and $V_{4}$
($V_{6}$) to perform a fast two-dimensional (2D) scan of $10\times 10$ pixels
and $100\times 100$ $\text{\,}\mathrm{mV}$ on device A (B). (c) An
illustrative diagram showing an early iteration of the algorithm. For clarity,
we restrict the schematic to the plane defined by the two fast plunger gates
of a device. Firstly, the algorithm performs an exploratory ramp along a
randomly chosen direction, where a Gaussian process chooses the starting gate
voltage coordinate to lie as close to, but still outside of, the rf
hypervolume as its confidence allows. Along the length of the exploratory
ramp, fast 2D scans are performed using gates $V_{2}$ ($V_{4}$) and $V_{4}$
($V_{6}$) for device A (B). The rf classifier evaluates each of these scans to
determine whether they show rf features or just noise. The first instance of
the rf classifier classifying features marks the point where the exploratory
ramp enters the rf hypervolume and the rf becomes sensitive. At this point,
the exploratory ramp is terminated, and then the algorithm explores the local
gate voltage region by random walk, again while taking 2D scans. The score
function evaluates these scans to determine the quality of the double-dot
features the algorithm has found.

##  II Methods

###  II.1 Device and rf readout

We demonstrate our algorithm in a Ge/SiGe depletion mode quantum dot array
operated as two double dot devices, labelled A and B, where five voltage gates
shaped the electrostatic potential (Fig. 1(a)). The fabrication process of
this device is similar to that described in Jirovec  _et al._ (2021). Ac and
dc signals were applied to the plunger gates for each quantum dot,
$V_{2},V_{4}$ and $V_{6}$, via bias-tees. A constant voltage of
$2\text{\,}\mathrm{V}$ was set to the splitter gates $V_{\text{sl}}$ and
$V_{\text{sr}}$ to isolate the triple dot array from a charge sensor, which
was not used for this work. Measurements were performed at
$50\text{\,}\mathrm{mK}$.

One device ohmic was connected to an L-matching network for rf readout. This
L-matching network was formed by a $92\text{\,}\mathrm{pF}$ decoupling
capacitor, a $2.2\text{\,}\mathrm{\SIUnitSymbolMicro H}$ inductor, and the
parasitic capacitance to the ground of the printed circuit board (PCB) and the
device, $C_{P}$ (Fig. 1(a))Ares  _et al._ (2016). These values were chosen to
approximately match the impedance of the double quantum dot device Vigneau
_et al._ (2022). The power of the rf drive at the device was fixed at
$-90\text{\,}\mathrm{d}\mathrm{B}\mathrm{m}$. The rf signal’s in-phase and
out-of-phase components were acquired with homodyne demodulation.

The device measurements consisted of fast two-dimensional gate voltage scans.
Plunger gates were swept in a raster pattern using an arbitrary waveform
generator (Fig. 1(b)). These scans had an amplitude and resolution of 100 by
$100\text{\,}\mathrm{mV}$ and $100$ by $100$ pixels, respectively. An
integration time of $1\text{\,}\mathrm{\SIUnitSymbolMicro s}$ per pixel was
used, enabling an entire scan to be completed in $10\text{\,}\mathrm{ms}$. It
was necessary to precompensate the waveforms responsible for the raster
pattern for the distortion introduced by the bias tees. See Appendix A for a
complete description of how these scans were performed and precompensated.

###  II.2 Algorithm

Our rf tuning algorithm is built from two components: a strategy to navigate
through the device’s gate voltage space and a score function to quantify the
quality of any double dot features found.

Previous strategies to navigate gate voltage spaces cannot be trivially
applied to rf measurements, as matching networks are only sensitive over a
small window of device impedances. These strategies used current measurement
to search for the current hypersurface, which marks the transition between the
device being ‘open’ and ‘closed’ Moon  _et al._ (2020); Severin  _et al._
(2021). Quantum dot features are found near this current hypersurface.

To circumvent this limitation, we formulate an rf hypervolume to navigate by.
This volume is defined as where the rf matching network is sensitive to the
changes in quantum device impedance. Outside the rf hypervolume, the matching
network is not sensitive, and the signal in the demodulated in-phase and out-
phase components falls below the noise floor. The current hypersurface and the
rf hypervolume do not necessarily intersect, and we can control the position
in the gate voltage space of the rf hypervolume by choice of the matching
network components. With an appropriate choice of components, quantum dot
features should be found near the rf hypervolume. In general, a range of
occupations can be accessed depending on the rf circuit, including the few-
carrier regime Petersson  _et al._ (2010). In our case, these features
correspond to carrier occupancies on the order of 10.

Our algorithm navigates to the rf hypervolume using exploratory gate voltage
ramps. These ramps slowly (at approximately
$1\text{\,}\mathrm{V}\text{\,}{\mathrm{s}}^{-1}$) move gate voltages towards
$0\text{\,}\mathrm{V}$ along a randomly-chosen radial direction in the space
defined by voltage gates $V_{1}-V_{5}$ for device A, or $V_{3}-V_{7}$ for
device B (Fig. 1(c)). The first of these ramps starts from a gate voltage
coordinate on the largest sphere bounded by the gate voltage upper limits,
which are 4V for all gates. This upper bound was chosen to span the full
operational range of the gate electrodes constrained by leakage currents.
Along the length of a ramp, the algorithm performs fast 2D gate voltage scans
(described in II.1) every $50\text{\,}\mathrm{mV}$. The algorithm evaluates
each of these scans using a classifier to determine whether the scan shows
only noise. When the ramp crosses into the rf hypervolume, the classifier will
cease to detect noise, and the ramp is terminated.

The algorithm then searches for double-dot features near the gate voltages
where the rf hypervolume was found using a random walk. Double dot features
are identified and quantified using a score function based on the discrete-
time two-dimensional Fourier transform (described in II.2.2). The random walk
is performed by perturbing a randomly-chosen gate voltage by an amount sampled
from a Gaussian distribution with a standard deviation of
$50\text{\,}\mathrm{mV}$.

The algorithm efficiently explores the gate voltage space by repeatedly
undertaking these exploratory ramps. With each ramp, it has the possibility of
finding new double-dot features which might score better than any dot features
found previously. Over time, with repeated exploratory ramps, the algorithm
explores an increasing fraction of gate voltage space; therefore, the score
associated with the highest-scoring double-dot features should improve over
the algorithm’s run time. Concurrently, each ramp’s termination gate voltage
coordinate is used to update a Gaussian process, as in Ref. Moon  _et al._
(2020). These observations of the hypervolume location inform the starting
gate voltage coordinate of subsequent ramps, allowing them to start closer to
the hypervolume and waste less time measuring regions of gate voltage space
where only noise can be observed. Therefore, as the algorithm improves its
model for the location of the hypervolume, it becomes more efficient at
exploring the gate voltage space.

The following sections describe the rf classifier and the score function used
by the algorithm.

####  II.2.1 The rf Classifier

Figure 2: (a) The demodulated in-phase $(I)$ and out-of-phase $(Q)$
quadratures of a 2D scan, showing double quantum dot features on the edge of
the rf hypervolume. (b) A scatter graph of the demodulated $I$ and $Q$
quadratures of the 2D scans in (a). We observe two principal directions of
variance. The information-bearing component can be found along $e_{1}$, while
the less informative component is along $e_{2}$. (c) A histogram of the $I$
and $Q$ data once projected onto $e_{1}$, shown in green. For comparison, we
show in grey an histogram of pre-measured noise projected similarly. The noise
has a Gaussian distribution, while the double dot features produce a
distinctly different distribution, providing a means to distinguish them.

The role of our rf classifier is to evaluate whether one of our fast 2D scans
shows noise or quantum dot features. If a scan at a given gate voltage
coordinate shows quantum dot features, then that coordinate must lie within
the rf hypervolume. Both in-phase ($I$) and out-of-phase ($Q$) components of
the rf signal could contain useful information about these features (Fig.
2(a)). This can be regarded as a blind signal separation (BSS) problem since
we wish to extract the information-bearing signal shared between the two
quadratures. We use a popular method to deal with this type of problem known
as principal component analysis (PCA). PCA considers the signal’s covariance
matrix to optimally find the direction in the $I-Q$ plane with the largest
variance and projects the data onto it Tharwat . Our algorithm uses PCA to
optimally project a pair of 2D scans measured in both quadratures into a main
one (Fig. 2(b)). This optimally-projected data is then classified using a
hypothesis test to compare the pixel value distribution of the projected scan
against a reference noise measurement. The noise is premeasured by taking 10
2D scans with all gate voltages set to their maximum value, far from the rf
hypervolume. As a way to illustrate this test, a histogram of the PCA-
projected pixel values of the scan in Fig. 2(a) is overlaid on the measured
reference noise distribution (Fig. 2(c)). Quantum dot features can be observed
in Fig. 2(a), and a long-tailed distribution captures this compared to the
noise in Fig. 2(c). This 2D scan would thus be classified as containing
features and terminate an exploratory gate voltage ramp.

The hypothesis test is a two-sample Kolmogorov-Smirnov test Hodges (1958);
Virtanen  _et al._ (2020). The null hypothesis is that the two distributions
are identical; therefore, the 2D scan contains just noise. The alternative
hypothesis is that the distributions differ, so the scan shows features
instead. Our critical p-value is $10^{-3}$, so p-values below this critical
value lead to the rejection of the null hypothesis in favour of the
alternative. This choice of p-value means that, on average, the null
hypothesis is falsely rejected only once in every $10^{3}$ scans.

####  II.2.2 Score Function

The role of the score function is to identify scans showing double dot
features by scoring them highly whilst scoring other features, such as those
from single dots, poorly. Score functions used in the literature to evaluate
current measurements have been based on neural networks Kalantre  _et al._
(2019); Zwolak  _et al._ (2020); Ziegler  _et al._ (2022), Hough transforms
Lapointe-Major  _et al._ (2020); Mills  _et al._ (2019), custom fitting models
Baart  _et al._ (2016), handcrafted heuristics Moon  _et al._ (2020), and even
ray-based classification Zwolak  _et al._ (2021). These score functions are
too noise-sensitive for fast measurements or cannot generalise to complex
signals. We designed a score function based on the discrete-time Fourier
transform of the 2D scans. This score function is naturally phase-independent
and thus can evaluate both in-phase and out-of-phase scans. Owing to the
properties of the Fourier transform, it is also robust against noise, allowing
us to benefit from the high bandwidth of rf measurements. See Appendix C for a
mathematical description of how the score function was evaluated and a
discussion of the noise tolerance inherited from the Fourier transform.

The discrete-time Fourier transform of the 2D scans quantifies any
periodicities present in these scans. Each periodicity results in a prominent
maximum in the Fourier transform, and the more periodic the features, the
higher the value of these maxima. Our score function searches the Fourier
transform of 2D scans for the two maxima resulting from double quantum dot
features. We constrain the search for these maxima so that we only consider
periodicities consistent with two quantum dots, each coupled enough to their
respective plunger gates, showing a few charge transitions over a
$100\text{\,}\mathrm{mV}$ plunger gate sweep. In practice, if the plunger
gates are $V_{x}$ and $V_{y}$, we search Fourier frequencies
$\nu_{x},\nu_{y}\in[0,12]$ a.u. and only considering maxima on opposite sides
of $\nu_{x}=\nu_{y}$. The score function then assigns the score as the
strength of the weaker of the two maxima. A scan showing single dot features
should result in just one maximum and other features in none, and thus low
scores will be assigned to such scans. A comparison of both double dot and
single dot features with respective Fourier transforms is displayed in Fig.
3(a-d).

Figure 3: (a-b) A 2D scan showing single dot features and corresponding time-
discrete Fourier transform of the standardised $I$ and $Q$ data. The single
dot features show one distinct periodicity, resulting in a single strong
maximum in the Fourier transform (with fringes resulting from the finite
window size). (c-d) A 2D scan showing double dot features and the
corresponding time-discrete Fourier transform of the standardised $I$ and $Q$
data. The double dot features show two distinct periodicities (and associated
harmonics), resulting in two strong maxima in the Fourier transform on either
side of the diagonal $\nu_{x}=\nu_{y}$ marked as a dashed red line.

##  III Results

To evaluate the performance of our algorithm both in devices A and B, we
performed 24 runs for a fixed time of an hour. This resulted in the algorithm
attempting to tune devices A and B 13 and 11 times, respectively, which
produced $\sim 2\times 10^{4}$ scans in total, each with an associated score.
To demonstrate the success of these runs, we need first to establish that the
score function is a good metric for the identification of double quantum dot
regimes and then show that the score increases as a function of tuning time
and that it does so faster than more straightforward search strategies.

To quantify how correlated our score function was with experts’ perception of
double dot quality, we created three categories: non-double, imperfect double,
and satisfactory double dot. Two human experts classified the scans according
to such labels without knowledge of the associated score. Fig. 4(a) shows a
stacked histogram of the scores, where each bar has been divided into colours
based on the human experts’ labels. We observe that non-double dot
(satisfactory quantum dot) features reliably receive low (high) scores.

Based on Fig. 3(a), we define two score thresholds to be 0.04 and 0.08. These
values are chosen to produce confusion matrices reflecting the score
function’s performance and have no bearing on the double quantum dot regimes
that the algorithm finds. Tables 1 and 2 show the confusion matrices for
devices A and B, respectively. Our score function produces very few false
positives even when faced with an overwhelming number of scans in which no
double quantum dot features are present whilst retaining the ability to
separate satisfactory from imperfect double quantum dots.

Table 1: Confusion matrix for device A. Device A (%) | Expert Classification
---|---
| Non-double
---
dot
| Imperfect
---
double dot
| Satisfactory
---
double dot
score $<0.04$ | 96.96 | 0.93 | 0.00
$0.04\leq$ score $<0.08$ | 0.19 | 1.43 | 0.06
score $\geq 0.08$ | 0.00 | 0.11 | 0.30
Table 2: Confusion matrix for device B. Device B (%) | Expert Classification
---|---
| Non-double
---
dot
| Imperfect
---
double dot
| Satisfactory
---
double dot
score $<0.04$ | 96.93 | 0.89 | 0.00
$0.04\leq$ score $<0.08$ | 0.30 | 1.41 | 0.09
score $\geq 0.08$ | 0.01 | 0.10 | 0.27

We now show the improvement of the score with time. Over its run time, the
algorithm explores a greater fraction of the rf hypervolume and discovers more
dot features. Each newly-discovered dot feature has a chance of scoring higher
than anything found before; therefore, on average, the score associated with
the highest-scoring features should improve over time. The purple curves in
Fig. 4(b-c) show the evolution of the best score found by the algorithm as a
function of run time for devices A and B, respectively. For a baseline
comparison, we ran a reduced version of our algorithm, which uses the
exploratory ramps but did not benefit from any of the Gaussian processes, the
rf classifier (which terminated the gate voltage ramps close to the rf
hypervolume), nor the subsequent random walk, shown in orange. We also ran a
random search algorithm that randomly chose a coordinate in gate voltage
space, set the gate voltages, and then performed a scan, shown in black.

We find that our algorithm outperforms these more simplistic versions. This is
likely because our algorithm spends a more significant fraction of its time
exploring regions that have the potential to correspond to double dot regimes.

Overall, the full algorithm found and identified 16.5 and 16.8 scans an hour
on average, corresponding to double dots (satisfactory and non-ideal) in
devices A and B, respectively. If we restrict our attention to satisfactory
double dots, the rates are 2.4 and 2.7 scans an hour. These tuning times are
the fastest recorded in laterally-defined quantum dots Moon  _et al._ (2020);
Severin  _et al._ (2021). See Appendix D for plots of satisfactory double dots
tuned by the algorithm in each hour-long run.

Figure 4: (a) A histogram of the double dot feature scores plotted on a linear
scale, where the colour coding of each bin indicates relative proportions of
the scans judged by human labelling to be non-double (blue), imperfect double
(green) and satisfactory double dot (red). The red double dot region lies
above the blue non-double dot region, whilst the green region corresponding to
satisfactory double dots lies above both the orange and the blue. For clarity,
we plotted only scores above $0.25$ on the linear scale due to the
overwhelming number of scores below this threshold; the inset shows a
histogram of all the scores on a log scale. (b-c) The highest score found for
increasing tuning times for algorithms which randomly search (grey), a reduced
version of the algorithm (orange) and the full algorithm (pink) in devices A
and B. To account for the statistical variation between runs, the curves
presented are the mean of many runs. Our algorithm explores the gate voltage
space more efficiently than a random search or random rays; therefore, its
best score improves fastest.

##  IV Conclusion

In conclusion, we have demonstrated the first algorithm for tuning double dots
exclusively utilising rf measurements. Such a result is an essential first
step toward the scalability of quantum devices. Our algorithm reliably tuned
satisfactory double quantum dots dot in approximately 15 minutes, the fastest
recorded time in laterally-defined quantum dots. In creating this algorithm,
we developed a robust method for distinguishing noise from rf features and
proposed a score function based on the Fourier transform with significant
noise tolerance.

Faster device tuning times could be achieved using a more informative Gaussian
prior for the Gaussian process, particularly one which encodes knowledge of
the gate architecture. Such a prior could use an electrostatic model of the
device Craig  _et al._ (2021). Furthermore, an additional Gaussian process and
Bayesian optimisation techniques could be used to choose the ramp directions.

To probe the charge states of large quantum dot arrays, the fast 2D scans
could be generalised to radial rays, as performed slowly in Chatterjee  _et
al._ (2021b); Zwolak  _et al._ (2021). Along with frequency multiplexing
Hornibrook  _et al._ (2014), our algorithm could allow for scalable tuning
and, thus, fully-automatic tuning of multiple small quantum dot arrays.

###### Acknowledgements.

We thank Nicholas Sim for providing help with the experiment and David Craig
for his feedback on the manuscript. This work was supported by the Royal
Society (URF-R1-191150), the EPSRC National Quantum Technology Hub in
Networked Quantum Information Technology (EP/M013243/1), Quantum Technology
Capital (EP/N014995/1), EPSRC Platform Grant (EP/R029229/1), the European
Research Council (Grant agreement 948932), the Scientific Service Units of IST
Austria through resources provided by the nanofabrication facility, the FWF-P
30207 and FWF-I 05060 projects, and grant number FQXi-IAF19-01 from the
Foundational Questions Institute Fund, a donor advised fund of Silicon Valley
Community Foundation.

##  Appendix A High Resolution Scans

In this section, we describe how fast two-dimensional scans were performed.
For this task, we used a Tektronix AWG (model 5014c) to feed a slow (red) and
a fast (blue) sawtooth waveforms to the device gates $V_{2}$, $V_{4}$ and
$V_{6}$, as shown in Fig. 1(b), such that their combination would produce a
raster scan in the gate voltage plane defined by $V_{2}$ and $V_{4}$ ($V_{4}$
and $V_{6}$) for device A (B), centred around the gates dc-voltage. Since the
amplitude of the ramp determines the scan size, this is ultimately limited by
the amount of attenuation in the cryostat lines. We chose our scans to be
$100\times 100$ $\text{\,}\mathrm{mV}$ in size.

Measurements were performed by rf reflectometry, using a Zurich lock-in
amplifier (UHFLI) to feed a frequency tone of 137.5 MHz to the L-matching
network connected to the device ohmic and demodulate the reflected component
of the signal. The resulting $I$ and $Q$ quadratures were finally recorded
using a fast digitiser card (Alazar tech ATS9440) using an AWG marker to
trigger the card acquisition such that it would be synchronised with the
voltage ramps. We used an integration time per pixel of $\tau_{\text{int}}=$
$1\text{\,}\mathrm{\SIUnitSymbolMicro s}$ and a resolution of 100 by 100
pixels, allowing us to perform a full 2D scan within $10\text{\,}\mathrm{ms}$.

The exact shape of the sawtooth waveform needed to be adjusted to avoid
artefacts introduced by high-pass filters in the setup. The action of a first-
order high-pass filter with a time constant $\tau$ on a signal $s(t)$ is
$s_{\text{filtered}}(t)=\mathcal{F}^{-1}[\mathcal{F}[s(t)](\omega)\cdot
H_{HP}(\omega)]$ where $H_{HP}(\omega)=i\omega\tau/(1+i\omega\tau)$, and
$\mathcal{F}[\cdot]$ denotes the time to frequency domain Fourier transform.
To correct for this effect, we pre-compensated the waveforms according to
$s_{\text{compensated}}(t)=\mathcal{F}^{-1}[\mathcal{F}[s(t)](\omega)\cdot
H_{HP}^{-1}(\omega)]$, such that the final shape reaching the device gates
would be exactly $s(t)$. Precompensating the desired waveform for this filter
requires only the knowledge of the filter’s time constant. At cryogenic
temperatures, we measured the time constants to be $1.043,1.051,1.025$
$\text{\,}\mathrm{ms}$ for gates $V_{2}$, $V_{4}$ and $V_{6}$, respectively.
Given that the slow sawtooth waveform period was $10\text{\,}\mathrm{ms}$,
precompensation was necessary.

##  Appendix B Projection of $I$ and $Q$ data with principal component
analysis

When performing the fast 2D scans, we measured both the in-phase and out-of-
phase quadratures. In general, both components could carry information. At the
price of recording both quadratures, PCA extracts the information-carrying
signal from the noise com (2010) and eliminates the need of phase adjustment.
Let us consider a 2D scan with in-phase data, $\vec{I}$, and the out-of-phase
data, $\vec{Q}$, where each index in these vectors corresponds to a pixel in
the 2D scan. To perform blind signal separation via principal component
analysis on this data, we obtained the eigenvector, $\vec{e}_{1}$,
corresponding to the largest eigenvalue of the covariance matrix,

$\left[\begin{array}[]{cc}\operatorname*{cov}(\vec{I},\vec{I})&\operatorname*{cov}(\vec{I},\vec{Q})\\\
\operatorname*{cov}(\vec{Q},\vec{I})&\operatorname*{cov}(\vec{Q},\vec{Q})\end{array}\right].$
(1)

We then obtained the extract the information bearing-signal,
$\vec{Y}_{\text{PCA}}$, by projecting the in-phase and out-of-phase data onto
this eigenvector, such that

$\vec{Y}_{\>\text{PCA}}=\vec{e}_{1}^{\;\intercal}\mathbf{X}\quad\text{where}\quad\mathbf{X}=\begin{bmatrix}(\vec{I}-\bar{I})\rightarrow\\\
(\vec{Q}-\bar{Q})\rightarrow\end{bmatrix},$ (2)

where $\bar{I}$ and $\bar{Q}$ denote the mean of $\vec{I}$ and $\vec{Q}$
respectively.

##  Appendix C The discrete time Fourier transform score function

This section describes how to evaluate the DTFT score function on the in-
phase, $\vec{I}$, and out-of-phase, $\vec{Q}$, data produced by a 2D scan of
resolution $N\times M$ pixels. Firstly, we combine the in-phase and out-of-
phase components into one complex dataset, to obtain
$\vec{Z}:=\vec{I}+i\vec{Q}\in\mathbb{C}^{N\times M}$. Then, we standardise
this complex data, $\vec{Z}$ according to
$\vec{\mathcal{Z}}:=(\vec{Z}-\bar{Z})/{\sigma_{Z}}$, where
$\bar{Z}:=\bar{I}+i\bar{Q}$ and
$\sigma^{2}_{Z}:=\operatorname*{var}[I]+\operatorname*{var}[Q]$. The reason
for standardising the data is so that the score function reflects the quality
of the double dot features and not how well-matched they are to the L-matching
circuit. Next, we compute the discrete-time Fourier transform of
$\vec{\mathcal{Z}}$ according to

$\text{DTFT}[\vec{\mathcal{Z}}]_{\delta\gamma}=\frac{1}{NM}\sum_{a=1}^{N}\sum_{b=1}^{M}e^{-2\pi
i(\nu^{x}_{\delta}a+\nu^{y}_{\gamma}b)}\mathcal{Z}_{ab},$ (3)

where $\vec{\nu}^{x}$ and $\vec{\nu}^{y}$ are a high-resolution linspace of
frequencies between 0 and 12, to obtain
$\text{DTFT}[\vec{\mathcal{Z}}]\in\mathbb{C}^{100\times 100}$. From here the
score, $s$, is computed according to

$s:=\min\left(\max|\text{DTFT}[\vec{\mathcal{Z}}]|_{\delta<\gamma},\max|\text{DTFT}[\vec{\mathcal{Z}}]|_{\delta>\gamma}\right).$
(4)

To paraphrase this formula, the score function finds the strongest periodicity
presented above and below the line $\nu_{x}=\nu_{y}$ and then assigns the
score to be the weaker of the two. In this manner, the score function favours
double dots where each dot is coupled most strongly to its plunger gate. If
there is just a single dot present, then there could be strong periodicity and
even perhaps harmonics on one side of the line $\nu_{x}=\nu_{y}$. However,
single-dot features cannot score highly without a dot coupled to the opposite
plunger, as the two periodicities necessary for a high score would not be
present.

As discussed in the following section, this score function inherits
considerable noise tolerance from the Fourier transform, such that the noise
in the score function is suppressed compared to the noise in 2D scans by a
factor of the square root of the number of pixels, $NM$. In our work, we used
$N$ and $M=100$, meaning that the noise in the elements of the DTFT was
suppressed by a factor of 100.

###  C.1 The discrete-time Fourier transform’s noise characteristics

If $\vec{Z}\in\mathbb{C}^{N\times M}$ is a complex dataset of $NM$ elements
and $\vec{N}\in\mathbb{C}^{N\times M}$ is the corresponding gaussian noise,
where each element is complex-normally distributed such that
$N_{ij}\sim\mathcal{\mathcal{C}N}(0,\sigma_{N})$ (meaning that the real and
imaginary parts of elements of $\vec{N}$ are independently normally
distributed with a mean of zero and variance of $\sigma_{N}^{2}/2$). Then

$\text{DTFT}[\vec{Z}+\vec{N}]_{\delta\gamma}\sim\mathcal{\mathcal{C}N}\left(\text{DTFT}[\vec{Z}]_{\delta\gamma},\frac{\sigma_{N}}{\sqrt{NM}}\right),$
(5)

where $\text{DTFT}[\cdot]$ is the discrete-time Fourier transform. Therefore,
the noise is suppressed by a factor of $1/\sqrt{NM}$ in the discrete-time
Fourier transform. This can be shown by using the additive nature of the
discrete-time Fourier transform. We can write
$\text{DTFT}[\vec{Z}+\vec{N}]_{\delta\gamma}=\text{DTFT}[\vec{Z}]_{\delta\gamma}+\text{DTFT}[\vec{N}]_{\delta\gamma}$,
where the time-discrete Fourier transform of the noise $\vec{N}$ is

$\text{DTFT}[\vec{N}]_{\delta\gamma}=\frac{1}{NM}\sum_{a=1}^{N}\sum_{b=1}^{M}e^{-2\pi
i(\nu^{x}_{\delta}a+\nu^{y}_{\gamma}b)}N_{ab}.$ (6)

As the complex-normal distribution is circularly-symmetric the multiplication
by the complex phase factor, $\exp(-2\pi
i(\nu^{x}_{\delta}a+\nu^{y}_{\gamma}b))$, leaves distribution unchanged, thus
$\exp(-2\pi
i(\nu^{x}_{\delta}a+\nu^{y}_{\gamma}b))N_{ab}\sim\mathcal{CN}(0,\sigma_{N})$,
which means that $\text{DTFT}[\vec{N}]_{\delta\gamma}$ can be regarded as the
mean of $NM$ samples from a complex-normal distribution. The central limit
theorem predicts that this mean will be distributed according to
$\text{DTFT}[\vec{N}]_{\delta\gamma}\sim\mathcal{\mathcal{C}N}\left(0,\sigma_{N}/\sqrt{NM}\right)$.
As there is no uncertainty in the Fourier transform of the noiseless signal,
$\vec{N}]_{\delta\gamma}$ the sum of the Fourier transform of the signal and
noise will be distributed according to equation 5.

##  Appendix D Outcomes from the hour-long tuning sessions

The highest scoring of the double dot features found by the algorithm in each
hour-long tuning session (Fig. 5).

Figure 5: The highest-scored double dot features found by the algorithm by the
end of each of the hour-long tuning sessions performed on devices A or B,
labelled in the top corner. For scans taken in device A (B), the horizontal
and vertical axis are gates $V_{2}$ ($V_{4}$) and $V_{4}$ ($V_{6})$
respectively, both of these gate voltages are swept over a
$100\text{\,}\mathrm{mV}$ range. The pixel values correspond to the PCA
projection described in appendix B.

## References

  * Vigneau  _et al._ (2022) Florian Vigneau, Federico Fedele, Anasua Chatterjee, David Reilly, Ferdinand Kuemmeth, Fernando Gonzalez-Zalba, Edward Laird, and Natalia Ares, “Probing quantum devices with radio-frequency reflectometry,” arXiv:2202.10516 (2022).
  * Hogg  _et al._ (2022) Mark R Hogg, Prasanna Pakkiam, Samuel K Gorman, Andrey V Timofeev, Yousun Chung, Gurpreet K Gulati, Matthew G House, and Michelle Y Simmons, “Single-shot readout of multiple donor electron spins with a gate-based sensor,” arXiv:2203.09248 (2022).
  * Madzik  _et al._ (2022) Mateusz T. Madzik, Serwan Asaad, Akram Youssry, Benjamin Joecker, Kenneth M. Rudinger, Erik Nielsen, Kevin C. Young, Timothy J. Proctor, Andrew D. Baczewski, Arne Laucht, Vivien Schmitt, Fay E. Hudson, Kohei M. Itoh, Alexander M. Jakob, Brett C. Johnson, David N. Jamieson, Andrew S. Dzurak, Christopher Ferrie, Robin Blume-Kohout, and Andrea Morello, “Precision tomography of a three-qubit donor quantum processor in silicon,” Nature 601, 348–353 (2022).
  * Moon  _et al._ (2020) H. Moon, D. T. Lennon, J. Kirkpatrick, N. M. van Esbroeck, L. C. Camenzind, Liuqi Yu, F. Vigneau, D. M. Zumbühl, G. A. D. Briggs, M. A. Osborne, D. Sejdinovic, E. A. Laird, and N. Ares, “Machine learning enables completely automatic tuning of a quantum device faster than human experts,” Nature Communications 11, 4161 (2020).
  * Baart  _et al._ (2016) T. A. Baart, P. T. Eendebak, C. Reichl, W. Wegscheider, and L. M. K. Vandersypen, “Computer-automated tuning of semiconductor double quantum dots into the single-electron regime,” Applied Physics Letters 108, 213104 (2016).
  * Mills  _et al._ (2019) A. R. Mills, M. M. Feldman, C. Monical, P. J. Lewis, K. W. Larson, A. M. Mounce, and J. R. Petta, “Computer-automated tuning procedures for semiconductor quantum dot arrays,” Applied Physics Letters 115, 113501 (2019), https://doi.org/10.1063/1.5121444 .
  * Zwolak  _et al._ (2020) Justyna P. Zwolak, Thomas McJunkin, Sandesh S. Kalantre, J.P. Dodson, E.R. MacQuarrie, D.E. Savage, M.G. Lagally, S.N. Coppersmith, Mark A. Eriksson, and Jacob M. Taylor, “Autotuning of Double-Dot Devices In Situ with Machine Learning,” Physical Review Applied 13 (2020).
  * Cai  _et al._ (2019) Zhenyu Cai, Michael A Fogarty, Simon Schaal, Sofia Patomäki, Simon C Benjamin, and John JL Morton, “A silicon surface code architecture resilient against leakage errors,” Quantum 3, 212 (2019).
  * Li  _et al._ (2018) Ruoyu Li, Luca Petit, David P. Franke, Juan Pablo Dehollain, Jonas Helsen, Mark Steudtner, Nicole K. Thomas, Zachary R. Yoscovits, Kanwal J. Singh, Stephanie Wehner, Lieven M. K. Vandersypen, James S. Clarke, and Menno Veldhorst, “A crossbar network for silicon quantum dot qubits,” Science Advances 4, 3960 (2018), https://www.science.org/doi/pdf/10.1126/sciadv.aar3960 .
  * Veldhorst  _et al._ (2017) M. Veldhorst, H. G. J. Eenink, C. H. Yang, and A. S. Dzurak, “Silicon CMOS architecture for a spin-based quantum computer,” Nat Commun 8, 1766 (2017).
  * Gonzalez-Zalba  _et al._ (2021) M. F. Gonzalez-Zalba, S. de Franceschi, E. Charbon, T. Meunier, M. Vinet, and A. S. Dzurak, “Scaling silicon-based quantum computing using CMOS technology,” Nat Electron 4, 872–884 (2021).
  * Laucht  _et al._ (2021) Arne Laucht, Frank Hohls, Niels Ubbelohde, M Fernando Gonzalez-Zalba, David J Reilly, Søren Stobbe, Tim Schröder, Pasquale Scarlino, Jonne V Koski, Andrew Dzurak, Chih-Hwan Yang, Jun Yoneda, Ferdinand Kuemmeth, Hendrik Bluhm, Jarryd Pla, Charles Hill, Joe Salfi, Akira Oiwa, Juha T Muhonen, Ewold Verhagen, M D LaHaye, Hyun Ho Kim, Adam W Tsen, Dimitrie Culcer, Attila Geresdi, Jan A Mol, Varun Mohan, Prashant K Jain, and Jonathan Baugh, “Roadmap on quantum nanotechnologies,” Nanotechnology 32, 162003 (2021).
  * Chatterjee  _et al._ (2021a) Anasua Chatterjee, Paul Stevenson, Silvano De Franceschi, Andrea Morello, Nathalie P. de Leon, and Ferdinand Kuemmeth, “Semiconductor qubits in practice,” Nat Rev Phys 3, 157–177 (2021a).
  * Burkard  _et al._ (2021) Guido Burkard, Thaddeus D Ladd, John M Nichol, Andrew Pan, and Jason R Petta, “Semiconductor spin qubits,” arXiv:2112.08863 (2021).
  * Stehlik  _et al._ (2015) J. Stehlik, Y.-Y. Liu, C. M. Quintana, C. Eichler, T. R. Hartke, and J. R. Petta, “Fast charge sensing of a cavity-coupled double quantum dot using a Josephson parametric amplifier,” Physical Review Applied 4, 014018 (2015).
  * Schupp  _et al._ (2020) F. J. Schupp, F. Vigneau, Y. Wen, A. Mavalankar, J. Griffiths, G. A. C. Jones, I. Farrer, D. A. Ritchie, C. G. Smith, L. C. Camenzind, L. Yu, D. M. Zumbühl, G. A. D. Briggs, N. Ares, and E. A. Laird, “Sensitive radiofrequency readout of quantum dots using an ultra-low-noise SQUID amplifier,” Journal of Applied Physics 127, 244503 (2020).
  * Jolliffe and Cadima (2016) Ian T Jolliffe and Jorge Cadima, “Principal component analysis: a review and recent developments,” Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 374, 20150202 (2016).
  * Pal  _et al._ (2013) Madhab Pal, Rajib Roy, Joyanta Basu, and Milton S. Bepari, “Blind source separation: A review and analysis,” in  _2013 International Conference Oriental COCOSDA held jointly with 2013 Conference on Asian Spoken Language Research and Evaluation (O-COCOSDA/CASLRE)_ (2013) pp. 1–5.
  * com (2010) “Contributors,” in  _Handbook of Blind Source Separation_, edited by P. Comon and C. Jutten (Academic Press, Oxford, 2010) pp. xxiii–xxiv.
  * Hodges (1958) J. L. Hodges, “The significance probability of the smirnov two-sample test,” Ark. Mat. 3, 469–486 (1958).
  * Scappucci  _et al._ (2021) Giordano Scappucci, Christoph Kloeffel, Floris A. Zwanenburg, Daniel Loss, Maksym Myronov, Jian-Jun Zhang, Silvano De Franceschi, Georgios Katsaros, and Menno Veldhorst, “The germanium quantum information route,” Nat Rev Mater 6, 926–943 (2021).
  * Watzinger  _et al._ (2018) Hannes Watzinger, Josip Kukučka, Lada Vukušić, Fei Gao, Ting Wang, Friedrich Schäffler, Jian-Jun Zhang, and Georgios Katsaros, “A germanium hole spin qubit,” Nat Commun 9, 3902 (2018).
  * Jirovec  _et al._ (2021) Daniel Jirovec, Andrea Hofmann, Andrea Ballabio, Philipp M. Mutter, Giulio Tavani, Marc Botifoll, Alessandro Crippa, Josip Kukucka, Oliver Sagi, Frederico Martins, Jaime Saez-Mollejo, Ivan Prieto, Maksim Borovkov, Jordi Arbiol, Daniel Chrastina, Giovanni Isella, and Georgios Katsaros, “A singlet-triplet hole spin qubit in planar Ge,” Nature Materials 20, 1106–1112 (2021).
  * Hendrickx  _et al._ (2020) N. W. Hendrickx, D. P. Franke, A. Sammak, G. Scappucci, and M. Veldhorst, “Fast two-qubit logic with holes in germanium,” Nature 577, 487–491 (2020).
  * Hendrickx  _et al._ (2021) Nico W. Hendrickx, William I. L. Lawrie, Maximilian Russ, Floor van Riggelen, Sander L. de Snoo, Raymond N. Schouten, Amir Sammak, Giordano Scappucci, and Menno Veldhorst, “A four-qubit germanium quantum processor,” Nature 591, 580–585 (2021).
  * Wang  _et al._ (2022) Ke Wang, Gang Xu, Fei Gao, He Liu, Rong-Long Ma, Xin Zhang, Zhanning Wang, Gang Cao, Ting Wang, Jian-Jun Zhang, Dimitrie Culcer, Xuedong Hu, Hong-Wen Jiang, Hai-Ou Li, Guang-Can Guo, and Guo-Ping Guo, “Ultrafast coherent control of a hole spin qubit in a germanium quantum dot,” Nat Commun 13, 206 (2022).
  * Ares  _et al._ (2016) N Ares, FJ Schupp, A Mavalankar, G Rogers, J Griffiths, GAC Jones, I Farrer, DA Ritchie, CG Smith, A Cottet, _et al._ , “Sensitive radio-frequency measurements of a quantum dot by tuning to perfect impedance matching,” Physical Review Applied 5, 034011 (2016).
  * Severin  _et al._ (2021) B Severin, Dominic T Lennon, Leon C Camenzind, Florian Vigneau, F Fedele, D Jirovec, A Ballabio, D Chrastina, G Isella, M de Kruijf, _et al._ , “Cross-architecture tuning of silicon and sige-based quantum devices using machine learning,” arXiv:2107.12975 (2021).
  * Petersson  _et al._ (2010) K. D. Petersson, C. G. Smith, D. Anderson, P. Atkinson, G. A. C. Jones, and D. A. Ritchie, “Charge and Spin State Readout of a Double Quantum Dot Coupled to a Resonator,” Nano Lett. 10, 2789–2793 (2010).
  * (30) Alaa Tharwat, “Independent component analysis: An introduction,”  17, 222–249, publisher: Emerald Publishing Limited.
  * Virtanen  _et al._ (2020) Pauli Virtanen, Ralf Gommers, Travis E. Oliphant, Matt Haberland, Tyler Reddy, David Cournapeau, Evgeni Burovski, Pearu Peterson, Warren Weckesser, Jonathan Bright, Stéfan J. van der Walt, Matthew Brett, Joshua Wilson, K. Jarrod Millman, Nikolay Mayorov, Andrew R. J. Nelson, Eric Jones, Robert Kern, Eric Larson, C J Carey, İlhan Polat, Yu Feng, Eric W. Moore, Jake VanderPlas, Denis Laxalde, Josef Perktold, Robert Cimrman, Ian Henriksen, E. A. Quintero, Charles R. Harris, Anne M. Archibald, Antônio H. Ribeiro, Fabian Pedregosa, Paul van Mulbregt, and SciPy 1.0 Contributors, “SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python,” Nature Methods 17, 261–272 (2020).
  * Kalantre  _et al._ (2019) Sandesh S Kalantre, Justyna P Zwolak, Stephen Ragole, Xingyao Wu, Neil M Zimmerman, Michael D Stewart, and Jacob M Taylor, “Machine learning techniques for state recognition and auto-tuning in quantum dots,” npj Quantum Information 5, 1–10 (2019).
  * Ziegler  _et al._ (2022) Joshua Ziegler, Thomas McJunkin, E. S. Joseph, Sandesh S. Kalantre, Benjamin Harpt, D. E. Savage, M. G. Lagally, M. A. Eriksson, Jacob M. Taylor, and Justyna P. Zwolak, “Toward Robust Autotuning of Noisy Quantum Dot Devices,” Physical Review Applied 17, 024069 (2022).
  * Lapointe-Major  _et al._ (2020) Maxime Lapointe-Major, Olivier Germain, J Camirand Lemyre, Dany Lachance-Quirion, Sophie Rochette, F Camirand Lemyre, and Michel Pioro-Ladrière, “Algorithm for automated tuning of a quantum dot into the single-electron regime,” Physical Review B 102, 085301 (2020).
  * Zwolak  _et al._ (2021) Justyna P. Zwolak, Thomas McJunkin, Sandesh S. Kalantre, Samuel F. Neyens, E.R. MacQuarrie, Mark A. Eriksson, and Jacob M. Taylor, “Ray-based framework for state identification in quantum dot devices,” PRX Quantum 2 (2021).
  * Craig  _et al._ (2021) David L Craig, H Moon, Federico Fedele, Dominic T Lennon, B Van Straaten, Florian Vigneau, Leon C Camenzind, Dominik M Zumbühl, G Andrew D Briggs, Michael A Osborne, _et al._ , “Bridging the reality gap in quantum devices with physics-aware machine learning,” arXiv:2111.11285 (2021).
  * Chatterjee  _et al._ (2021b) Anasua Chatterjee, Fabio Ansaloni, Torbjørn Rasmussen, Bertram Brovang, Federico Fedele, Heorhii Bohuslavskyi, Oswin Krause, and Ferdinand Kuemmeth, “Autonomous estimation of high-dimensional coulomb diamonds from sparse measurements,” arXiv:2108.10656 (2021b).
  * Hornibrook  _et al._ (2014) J. M. Hornibrook, J. I. Colless, A. C. Mahoney, X. G. Croot, S. Blanvillain, H. Lu, A. C. Gossard, and D. J. Reilly, “Frequency multiplexing for readout of spin qubits,” Appl. Phys. Lett. 104, 103108 (2014).