\section{Introduction}
Semiconductor qubits controlled by gate electrodes provide a promising route to large scale quantum computation~\cite{Cha21}.
In these types of systems, individual electrons are confined in quantum dots (QDs) and qubits are formed, for example, by the spin or charge degrees of freedom of the electrons.
Controlling and manipulating the quantum state of the qubits is then achieved by applying external voltages via gate electrodes (gate voltages).
By carefully tuning these gate voltages, transitions between different states can be realized, which can then serve as the computational states for quantum computations.
The fabrication of large arrays (currently going up to 16 quantum dots~\cite{mills2019shuttling,volk2019loading, zajac2016scalable, borsoi202216}) creates a new challenge, as  the simultaneous manual tuning of many gate voltages quickly becomes infeasible for such large arrays. 

Changing the gate voltages and measuring the quantum dot occupations in the systems ground state leads to a so-called charge-stability diagram (CSD), mapping the high-dimensional voltage space to that of the number of electrons on each dot.
Constructing a CSD is typically done by performing many two-dimensional raster scans of pairs of gate voltages~\cite{oakes2020automatic, volk2019loading}. 
Based on these raster scans higher-precision line-scans are performed around areas where the QD occupations change to estimate the normal of the charge transition. 
Knowing the location and normal of the transitions allows one to define compensated gate voltages (i.e. linear combinations of voltages, see below).
While this approach is feasible for small arrays of quantum dots, mapping out the CSD for large arrays quickly becomes a very tedious task.
To date, hand-tuned loading and shuttling of electrons in QD arrays has been achieved with as large as 8 or 9 quantum dots~\cite{
mills2019shuttling,volk2019loading, zajac2016scalable}. 
Automating this process instead is highly desirable, freeing up valuable time towards more impactful experiments.
Moreover, since the number of control voltages grows linearly with the number of quantum dots, hand-tuning their values becomes more and more challenging due to cross-talk between the dots.
Only recently did we see the emergence of automatic tuning algorithms, often implemented using machine-learning \cite{durrer2020nn,chatterjee2021autonomous,krause2021estimation,oakes2020automatic,nguyen2021deep, zwolak2020ray, zwolak2021ray}. 
These approaches were used only for small arrays, and still lag behind the results achievable via manual tuning. 

In this article, we introduce an algorithm for the automated exploration of large quantum dot arrays, focusing on the identification of transitions between different charge states without the need to label transitions.
In particular, the algorithm is designed to reliably find single- and multi-electron charge-state transitions in large quantum dot arrays using simple measurement primitives.
We only rely on the use of line scans to detect the nearest transition in a chosen direction, and use a 
machine-learning model to estimate transitions from measurements.
Compared to previous works \cite{chatterjee2021autonomous,krause2021estimation}, we will also make use of explicit knowledge of the underlying physics to constrain the problem to make the algorithm more efficient.
We restrict ourselves to quantum dot devices described by the constant interaction model~\citep{schroer2007electrostatically}, appropriate for instance for arrays of gate-controlled spin qubits. In such quantum dot arrays, charge configurations are dictated by Coulomb energies (capacitive coupling), and not by relatively weak tunnel coupling. We further assume that in such a device, we can realize and measure transitions between charge-states reliably up to some precision $\delta$.

In the rest of this article we first introduce quantum dot arrays more formally, and proceed with describing the algorithmic solution.
We then analyze several relevant scenarios of device layouts and demonstrate the algorithm's performance.
The outlook then considers possible relaxations of the assumptions we make along the way.

\section{Quantum dot arrays and Coulomb diamonds}
We consider devices of $N$ quantum dots with significant charging energies ($E_c \gg k_BT$) that are weakly-tunnel coupled to each other and to $G$ external gate electrodes, such as the one schematically shown in Fig.~\ref{fig:device_sketch}. 
The gate voltages define a high-dimensional space that can be divided into regions where the device assumes a certain ground state  of QD occupations.
These extended regions are referred to as Coulomb diamonds, because it is the Coulomb blockade effect that prevents additional electrons from tunneling onto a quantum dot unless the potential is large enough~\cite{semiconductorsIhn}.
\begin{figure}
	\centering
	\includegraphics[width=0.42\textwidth]{Device_Sketch.pdf}
	\caption{Schematic of a $2\times 3$ array of capacitively coupled quantum dots, each occupied by $n_i$ electrons, capacitively coupled to external gate voltages $v_j$ described by a matrix $C^{DG}_{ij}$, and coupled to each other via capacitances $C^{DD}_{ik}$. For simplicity of visualisation, the couplings $C^{DG}_{ij}$ are only shown for $i=j$. In this device, the number of gate-voltages $G$ equals the number of dots $N=6$.}
	\label{fig:device_sketch}
\end{figure}
The Coulomb diamonds and transitions between them can be described well using the constant interaction model~\cite{schroer2007electrostatically}, which, given gate voltages $v = (v_1, \ldots, v_G) \in \mathbb{R}^G$ and quantum dot occupations $n = (n_1, \ldots, n_N) \in \mathbb{N}^N$, assigns the system a free energy given by:
\begin{equation}\label{eq:free_energy}
    F(n, v) = \frac 1 2 (\lvert e\rvert n - C^{DG} v)^T (C^{DD})^{-1} (\lvert e\rvert n -C^{DG} v)\enspace.
\end{equation}
Here, $\lvert e\rvert$ is the charge of an electron, $C^{DG} \in \mathbb{R}^{N \times D}$ is a matrix whose entries $C^{DG}_{ij}$ store the capacitance between the $i$th dot and $j$th gate, and $C^{DD}$ is a matrix whose off-diagonal elements are inter-dot capacitances.
The diagonal elements of $C^{DD}$ are chosen such that $\sum_j^N C^{DD}_{ij} -\sum_k^G C^{DG}_{ik}=0$.

A Coulomb diamond is then a volume $P_n$ in the high-dimensional voltage space for which the system ground state assumes a set of occupations $n$:
\begin{align}\label{eq:const_interaction_trans}
   
    P_n &= \{v \mid \arg \min_{r} F(r, v)=n\}. \notag\\
    &=\{v \mid F(n,v) - F(r, v) \leq 0,\, \forall r \in \mathbb{N}^{N} \}. \notag\\
    &=\left\{ v \mid (r-n)^T \underbrace{(C^{DD})^{-1}C^{DG}}_A v + b_{nr} \leq 0,\, \forall r \in \mathbb{N}^{N} \right \}. \notag \\
    &\approx \left\{ v \mid t^TA v + b_{n,n+t} \leq 0,\, \forall t \in \{-1,0,1\}^{N}  \right \}
\end{align}
These inequalities describe high-dimensional half-spaces, whose boundaries form the facets of a convex polytope $P_n$. 
That is, for a given polytope $P$ with equations $\{ W_k v+ b_k \leq 0, k=1,\dots \}$, the facets are defined by the sets $\{v\in P |W_kv+b_k = 0\}$, $k=1,\dots$.
The $W_k$ define the normals of the facets, and the $b_k$ their offsets from the origin.
In Eq.~\ref{eq:const_interaction_trans}, the term $b_{nr}=\frac {\lvert e\rvert} 2 (n^T(C^{DD})^{-1}n-r^T(C^{DD})^{-1}r)$ does not depend on the gate voltages $v$.
The vector $t = r-n$ describes the \emph{changes} in occupation of the individual quantum dots when crossing a facet, and so the inequalities can be intuitively understood as  representing a transition $t$ from a state $n$ to a state $n+t$.
In the rest of this manuscript we limit the elements $t_i\in \{-1,0,1\}$, because multiple electrons entering or leaving a dot simultaneously is a process that can be neglected in practise.

Of special importance is the polytope $P_0$, the set of boundaries of the state with zero electrons on all dots. 
In this polytope, all transitions add a single electron on a quantum-dot, i.e. the $t$ are unit vectors, and the normals of the polytope facets are the (scaled) rows of the matrix $A$.
Knowing the (scaled) matrix $A$ allows one to compute compensated control voltages $u$, which are related to the gate voltages via a transformation matrix $v = U u$. The coordinate system of $u$ is chosen such, that $u_i$ only affects the potential of dot $i$. Knowledge of $u$ allows the device user to mitigate capacitive cross talk within the array, which simplifies the task of tuning the device. Geometrically, when measured in coordinates of $u$, the normal of the transition that adds one electron to the $i$th dot is parallel to the $i$th standard basis vector.

In the constant interaction model, the compensated control voltages can be computed by choosing $U$ such, that $A\cdot U$ is a diagonal matrix, for example by choosing $U$ as the (pseudo-)inverse of $A$. This can be verified by inserting the definition of $u$ into equation \eqref{eq:const_interaction_trans}, leading to
\[ 
    P_n = \left\{ v \mid t^T A U u  + b_{n,n+t} \leq 0,\, \forall t \in \{-1,0,1\}^{N}  \right \}\enspace.
\]
When $t$ is the $i$th standard basis vector, then $tAU$ is just a multiple of $t$. The compensated control voltages are sometimes referred to as virtual gates in the literature~\cite{hensgens2017quantum} and use of this technique was paramount in the recently achieved hand-tuned control over 8 and 9 qubits~\cite{volk2019loading,mills2019shuttling}.

\begin{figure}[!t]
	\centering
	\includegraphics[width=0.6\textwidth]{zero_state_alg.pdf}
	\caption{Graphical depiction of the algorithm for a double dot system. Starting with an initial dataset, the algorithm repeats steps a) and b) until termination. a) Algorithm fits model (black lines) to separate points inside the target state (blue) and outside (orange). b) Line-searches are performed from a common starting point (red dot) in the direction of a point of each facet of the model (red arrows). Line-searches produces new pairs of samples (orange/blue dots) with the true boundary (dotted lines) in-between. The distance between red and blue dots is smaller than $\delta$. See text for more details. }
	\label{fig:alg_high_level}
\end{figure}


\section{Method}\label{sec:method}
For our method, we assume that the device follows the constant interaction model and that it is equipped with non-invasive charge sensors that can detect the occurrence of a transition reliably without influencing the underlying capacitances. 
Further, we assume that the user can provide a line-search procedure that performs a raster scan along the line between two points of gate-voltages $v_{\text{start}}$ and $v_{\text{end}}$, and uses the charge sensors to detect the existence of a transition reliably along this line. We further assume, that the user can find two points $v^-, v^+$ along this raster scan that bound the position of the transition (i.e. $v^-$ and $v^+$ lie on opposite sides of a transition) and that the user can give a bound of the maximum distance $\lVert v^- -v^+ \rVert\leq \delta$.

Our task is to estimate the facets of a polytope $P_n$ as best as we can, so that we can identify the voltage regions corresponding to fixed quantum dot occupations.
While $P_n$ has, in the worst case, $3^N-1$ facets, in practice only a small set of transitions is of interest. 
Our algorithm therefore starts with a pre-selected list of such transitions that are of relevance to an experiment.
For example, we may ask for all single-electron transitions from a specific charge state (e.g. $P_0$), or transitions that leave the number of electrons unchanged. In many applications, the set of all "one-electron transitions" is relevant, i.e., transitions in which exactly one electron moves from one dot to another, or from a reservoir to a dot (or vice versa).

We do not assume that the list is an exhaustive enumeration of all transitions present in the target polytope, nor that all transitions on the list exist. 
The algorithm will then use data gathered by the line-searches to select a subset of existing transitions from the list and produce a new set of candidate line-search directions in order to refine the dataset and fit the polytope again. 
This process is repeated until all transitions selected are well-supported by the measurements and the transitions which are not selected are similarly ruled out with high certainty.


The algorithm can be described by alternating two steps, depicted in Figure~\ref{fig:alg_high_level}. 
First, an initial dataset $\mathcal{D}$ of measurements is provided by the user. 
The dataset consists of pairs of points $(v^-_i,v^+_i)$, $i=1,\dots, \ell$, bounding a transition of interest. 
These pairs of points are acquired by the aforementioned line-search procedure.  
The initial dataset can contain data from line-searches acquired by picking random directions, but may also include educated guesses of informative search directions.

The algorithm then proceeds to estimate the set of transitions using a regularized maximum-likelihood approach (see section~\ref{sec:MLmodel} below).
This estimation can be erroneous and may represent a local optimum, or the initial dataset may have a flaw that makes estimation difficult, e.g., it does not feature enough data points on each transition to estimate it reliably. 
Because of this, we use the estimated model to acquire additional points. 
We sample a point on each estimated facet and then perform a line-search from a point inside $P_n$ through this point (see Figure~\ref{fig:alg_high_level}b and Appendix~\ref{app:sampling}). 
If the estimated model is correct, the points returned by the line-search will bound the estimated transition. 
Otherwise, they are an example of an error of the model and can be used to update the estimate. 
In either case, they are added to $\mathcal{D}$ and a new model is estimated. 
This loop is repeated until we are confident that all transitions found are either correct or too small to be estimated reliably. For the former, we count the number of pairs $(v_i^-,v_i^+)$ that are separated by a transition and for the latter, we use the radius of the largest $G-1$ dimensional hypersphere that can be inscribed into the facet. For details, we refer to Appendix~\ref{app:termination}.

In the following, we describe the model we have chosen to represent and learn a polytope and how we make it practical for searching relevant transitions of a given state. We will focus on the general idea of the model, while the detailed description are given in the appendices. Further, a reference implementation of the algorithm is available at \cite{code_repository}.

\subsection{Estimating the model}\label{sec:MLmodel}
\begin{figure}
	\centering
	\begin{subfigure}{0.3\textwidth}
        \centering
        \includegraphics[width=\textwidth]{model100}
        \caption{$\lVert W_k\rVert = 100$}
        \label{fig:model100}
    \end{subfigure}
    \begin{subfigure}{0.3\textwidth}
        \centering
        \includegraphics[width=\textwidth]{model1000}
        \caption{$\lVert W_k\rVert = 1000$}
        \label{fig:model1000}
    \end{subfigure}
    \begin{subfigure}{0.3\textwidth}
        \centering
        \includegraphics[width=\textwidth]{model_moved}
        \caption{removed facet}
        \label{fig:model_moved}
    \end{subfigure}
	\caption{Visualization of classifier $h(v)$ for several choices of model parameters. Black lines indicate the contour $h(v)=0$ while blue (red) pixels indicate $h(v)<0$ ($h(v) > 0$). Figures a) and b) use scaled parameters of $W_k$ and $b_k$ from a polytope of the (1,1) state of a simulated double dot. Parameters are scaled such that the norms $\lVert W_k\rVert$ have fixed values. For small norms in a), $h(v)=0$ has a rounded shape, while for larger norms in b) the shape $h(v)=0$ approximates the underlying convex polytope well. Figure c) demonstrates the removal of a facet by moving it outside the polytope (dashed line). }
\end{figure}

To learn a polytope $P_n$ with equations $W_k v+ b_k \leq 0$, $k=1,\dots,M$ we have to identify the correct set of transitions and then estimate their parameters.
We start with a dataset of point pairs of gate voltages, 
$\mathcal{D}=\{ v_i^+,v_i^-\in \mathbb{R}^G, i=1,\dots,\ell\}$ with
$\lVert v_i^+-v_i^-\rVert\leq \delta$, and assume that $v_i^-$ belongs to the state $n$, while $v_i^+$ belongs to a different state. 

We formalize learning the polytope from data as a binary classification task, where the goal is to find a classifier $h$ that separates all $v_i^-$ from $v_i^+$. 
We derive this model classifier from a multi-class logistic regression model, the construction of which is fully contained in Appendix~\ref{app:model}:
\begin{align}\label{eq:model}
    h(v)=\log \sum_{k=1}^M \exp(W_k^Tv+b_k)\enspace.
\end{align}
A point is deemed as inside the polytope if $h(v)<0$, and outside otherwise. 
This model is well suited for modeling convex polytopes as, when $\lVert W_k \rVert$ becomes large, the set $\{v\in \mathbb{R}^G|h(v) \leq 0\}$ approximates a convex polytope arbitrarily well. 
This is demonstrated in Figures~\ref{fig:model100}\&\ref{fig:model1000}. Figure~\ref{fig:model_moved} demonstrates that an unnecessary facet can be removed from the model by choosing its offset $b_k$ such, that it is moved away from the polytope.

This proposed model is not practical however, as the number of parameters of $W$ is large and there is no easy way to infer which transition is modeled by a given row $W_k$. 
We can considerably improve on this by making use of physical knowledge of the transitions owing to the constant interaction model as follows. 
In the constant interaction model, following \eqref{eq:const_interaction_trans}, each transition has a normal $W_k \propto t_kA$. Thus, if $A$ were known, $W$ could be computed easily.

It is possible to estimate a row $A_k$ up to a factor by identifying a transition that adds an electron on the $k$th dot and measuring its position and normal (in voltage space). However, the length of the normal vector remains unknown. To represent this, we decompose $A = \Lambda \Gamma$, where $\Lambda$ is a real diagonal matrix and $\Gamma$ stores the estimated normal vectors.

The matrix $\Gamma$ can be learned efficiently, by measuring the transitions of the polytope $P_0$. 
This polytope has exactly $N$ transitions, one for each dot, and the transitions are $t_k=(0,\dots, 0,1,0,\dots)$, i.e. the unit vector with a 1 at index $k$. 
Consequently, $t_k^TA=A_k$, and hence identifying the transitions of $P_0$ gives a direct estimate of $\Gamma$.
The entries of $\Lambda$ can then be estimated by measuring transitions of electrons from one position in the array to the next. Moving an electron from position $i$ to position $j$ leads to a transition with normal $t^TA= A_j-A_i$. Using an estimate of the normal, we can find $\Lambda$ such, that $t^TA\propto t^T\Lambda \Gamma = \Lambda_{jj}\Gamma_j - \Lambda_{ii} \Gamma_i$. 


This logic needs to be translated into our fitting code and model Eq.~\eqref{eq:model}. To obtain an estimate of $\Gamma$, we set $W=\Gamma$, and train the model using regularized maximum likelihood using data $v_i^-, v_i^+$ obtained from the device initialized in state $n=0$.
This problem is non-convex and highly multi-modal, often resulting in getting trapped in bad fits, especially with only a small number of data points. 
For the minimization therefore, we include a regularizer $\Omega(\Gamma)$ that steers the optimizer away from the bad local optima using prior knowledge of the task.
Knowledge of the physics of the system comes in at this step, where we make use of two properties that are present in common device designs:
\begin{enumerate}
    \item For all transitions, we expect the absolute voltage value required to add an electron on the $k$th dot to be small 
    \item We assume that the capacitive coupling between a dot (i.e. its plunger gate) and its gate electrode dominates the influence of all other gate electrodes. That is, the angle between the normal of the transition adding an electron on the $i$th dot and the coordinate axis of the corresponding gate voltage is small.
\end{enumerate}
More details on how these are implemented are in Appendix~\ref{app:regularization}.

Having estimated $\Gamma$, we can now proceed with the estimation of $W$ (and $b$) of the target polytope $P_n$. 
We decompose $W=c T A= c T \Lambda \Gamma$, with $T$ a matrix whose rows are the transitions $t_k$ we are interested in, and $c$ a real diagonal matrix that enables sharpening the polytope (see Figure~\ref{fig:model1000}).
A single transition is thus described as 
\begin{equation}\label{eq:WPn}
W_k=c_{kk}\sum_{i=1}^N t_{ki} \Lambda_{ii} \Gamma_i\enspace.
\end{equation}

Together with a good estimate of $\Gamma$, we can now learn $P_n$ using regularized maximum likelihood given the pairs of data $v_i^-, v_i^+$ obtained from line-searches with the device initialized in state $n$. 
For simplicity, we assume that the rows $\Gamma_k$ are normalized s.t. $\lVert \Gamma_k \rVert=1$. 
With this, we keep $\Gamma$ fixed and only find the values of the diagonal entries of $c$, $\Lambda$ and the vector of offsets $b$. 
We again use regularized maximum likelihood on the probability of assigning the right class label to the points in a point-pair.
Here, too, we add a regularizer that we define in more detail in Appendix~\ref{app:regularization}. 

Due to the multi-modality of the solution, the optimisation performance varies a lot with the initialization for parameters $c,\Lambda$ and $b$. 
If we have no prior knowledge of a possible shape of the polytope (e.g., in the first iteration), we pick $c_{kk}=1/\delta$, $\Lambda=I_N$ and $b_k=-c_{kk}\max_i^\ell W_k^T v_i^{-}$. 
This choice of $b_k$ ensures that all points $v^{-}$ are classified correctly by the initial model and each facet is close to points in the dataset. If we already have a previous estimate, we choose $\Lambda$ as the same value from the previous iteration. 

\section{Results}\label{sec:experiments}
We consider five different scenarios to test and demonstrate our algorithm in.
Each scenario is characterized by a different 2D device geometry (c.f. Figure~\ref{fig:device_sketch}) and a different set of target states or transitions.
In all examples we assume for simplicity that the number of gates equals the number of dots and thus $G=N$.
For each scenario we generate devices with different capacitance matrices $C^{DD}$ and $C^{DG}$, using the scheme given in Appendix~\ref{app:generation}. The scheme depends on a parameter $\rho$ that changes the interaction strength (or cross-talk) between gates and dots. We generate $20$ devices with weak ($\rho=1$) and strong ($\rho=3$) interaction strength, and estimate their ground truth polytopes (see Appendix~\ref{app:baseline}). The interaction strengths are chosen such that the polytopes obtained from $\rho=1$ devices are similar to many practical devices, while for $\rho=3$, the resulting polytopes have large distortions, which are likely larger than encountered in practice and more difficult to estimate, especially for the task of estimating $\Gamma$. The generated device parameters are chosen such that the gate-voltages are measured in V and the size of the polytopes are around 0.2 (i.e., a range of 200mV) along each axis.
We consider the following scenarios:
\begin{description}
    \item[$S_1$] As an example for a sequence of operations involving several states, we select a 3x2 device and aim to find all transitions needed to exchange a pair of electrons. 
    This exchange is done in such a way that spin can be preserved and thus electrons can not be positioned on the same dot, due to the Pauli exclusion principle. 
    We choose a sequence of transitions that requires computing polytopes for the following 6 states (the corresponding transitions $t$ are just the differences between consecutive states). 
    Colors indicate the intended electron positions after an operation:
    \[
        \begin{pmatrix}
            \color{red}{1} & 0\\
            0 & 0\\
            \color{blue}{1} & 0        
        \end{pmatrix}
        \rightarrow
        \begin{pmatrix}
            \color{red}{1} & 0\\
            \color{blue}{1} & 0\\
            0 & 0        
        \end{pmatrix}
        \rightarrow
        \begin{pmatrix}
            \color{red}{1} & 0\\
            0 & \color{blue}{1}\\
            0 & 0        
        \end{pmatrix}
        \rightarrow
        \begin{pmatrix}
            0 & 0\\
            \color{red}{1} & \color{blue}{1}\\
            0 & 0        
        \end{pmatrix}
        \rightarrow
        \begin{pmatrix}
            0 & 0\\
            0 & \color{blue}{1}\\
            \color{red}{1} & 0        
        \end{pmatrix}
         \rightarrow
        \begin{pmatrix}
            0 & 0\\
            \color{blue}{1} & 0\\
            \color{red}{1} & 0        
        \end{pmatrix}
         \rightarrow
        \begin{pmatrix}
            \color{blue}{1} & 0\\
            0 & 0\\
            \color{red}{1} & 0        
        \end{pmatrix}
    \]
    Here, for each computed polytope, $T$ contains all transitions that add or remove an electron as well as all transitions of a single electron to a different location.
    \item[$S_2$] As an example for a more complex transition involving multiple electrons, we consider a 3x2 device with 3 electrons in a zig-zag configuration. The goal is to find the transition that lets all electrons simultaneously change side, i.e. we need to find the transition
    \[
        \begin{pmatrix}
            1 & 0\\
            0 & 1\\
            1 & 0        
        \end{pmatrix}
        \rightarrow
        \begin{pmatrix}
            0 & 1\\
            1 & 0\\
            0 & 1        
        \end{pmatrix},\;
        t= \begin{pmatrix}
            -1 & 1\\
            1 & -1\\
            -1 & 1        
        \end{pmatrix}\enspace.
    \]
    To find this facet reliably, we need to include most facets that are near it or intersect with it. However, we don't want to include all transitions, since on larger devices, the total set of all possible transitions would be too large. As a result, we choose $T$ to include the following transitions:
    \begin{itemize}
        \item $t$ and transitions that can be obtained from $t$ by setting one entry to 0
        \item Transitions that affect up to two electrons (and 4 locations in the array) simultaneously. 
    \end{itemize}
    The second group includes transitions in which up to two electrons move from one dot to another, or from a reservoir to a dot (or vice versa).
    
    \item[$S_3$] In a 3x3 device with one electron on each dot, find all transitions that add, remove or move a single electron. $T$ is chosen using the same principle as in $S_1$.
    \item[$S_4$] The same scenario as $S_3$ with a 4x4 device.
    \item[$S_5$] We consider a 4x4 device with one electron on each dot. 
    However, in this case we assume that we can completely detach the reservoir and thus the device can not exchange electrons with it. 
    Technically, this means that the computed polytopes no longer represent ground states of the device.
    
    This way only transitions are possible that keep the amount of electrons constant, i.e., transitions $t$ that move electrons within the array, thus $t^T\mathbbm{1}_N=0$ (here, $\mathbbm{1}_N$ denotes the all-one vector with $N$ entries). 
    These polytopes are special as they are unbounded independent of the state the device is in, which we show in the following. 
    Assuming $t^T\mathbbm{1}_N=0$, we can compute the normal of a transition $W_k=t_k^TA$. 
    If we take the direction $q=A^{-1}\mathbbm{1}_N$, we obtain 
    \[
        W_k^Tq = t_k^TA A^{-1} \mathbbm{1}_N = t_k^T\mathbbm{1}_N=0 \enspace.
    \]
    Thus, the normals of all facets of the polytope are orthogonal to the direction $q$. 
    From this follows, that for any point $v \in P_n$, $v+tq \in P_n$, $\forall t \in \mathbb{R}$.
    
    To solve this issue and to obtain bounded polytopes, after estimating $\Gamma$, we choose $N-1$ independent vectors orthogonal to $\hat{q}=\Gamma^{-1}\mathbbm{1}$ and find the projected polytope in the resulting $N-1$ dimensional subspace. 
    This projection can be computed using a householder reflection matrix that maps the direction $q$ to the first basis vector that can then be discarded. 
    While the shape of the resulting polytope depends on the chosen $\Gamma$, the polytope will always represent a finite polytope, as long as $q^T\hat{q} \neq 0$.
    
    For $T$, we consider all single electron inter-dot transitions in this example. Due to the fact that the large transitions that include exchanges with the reservoir cannot take place in this scenario, all $ N\cdot (N-1)=240$ such transitions are present on the ground truth polytope.
\end{description}

In all scenarios, we use our algorithm to estimate $\Gamma$ and the relevant polytopes, and use two different choices, $\delta \in \{0.001, 0.002\}$, for the line-search precision.
To prevent endless runs we terminate either when, for computing $\Gamma$, more than 4000 line-searches were conducted, or if for the individual polytopes in a scenario more than 15000 line-searches were conducted. 
As lower bound for voltages considered by line-searches, we take $v_i \geq -2$. 
As a lower bound for the size of a facet to be considered we require the radius $r$ of the largest hypersphere that can be inscribed into it to  be $r \geq r_{\min}=2\delta$.

As initial datasets, for estimating $\Gamma$ in a device with $N$ dots, we sample $4N(N+5)$ initial points, a number which worked well empirically in our experiments. 
As starting point, we choose $v=-2\mathbbm{1}_N$, the lower bound in all directions and we sample in direction $p=\exp(2y)$, where $y\in \mathbb{R}^G$ is a standard normal variable. 
This prevents sampling in directions that can violate the lower bound. 
This strategy is designed to give enough informative measurements that allow the algorithm to terminate quickly.

For the target polytopes $P_n$, we sample $N^2$ initial points. 
We first create a starting point for line-searches by sampling a point close to the boundary. 
This is uninformative and thus reflects the fact that we often have no good estimate for the center of the polytope.
Then we sample direction vectors uniformly on the unit sphere.
This setup only gives little information of the small facets of the polytope, and thus finding the small facets purely relies on our ability to sample informative candidate points during the run of the algorithm.

\begin{figure}
    \centering
    \begin{subfigure}{0.7\textwidth}
        \centering
        \includegraphics[width=\textwidth]{zero_angle_err}
        \caption{Error of normals $\Gamma$}
        \label{fig:error_angle}
    \end{subfigure}\\
    \begin{subfigure}{0.7\textwidth}
        \centering
        \includegraphics[width=\textwidth]{zero_vangle_err}
        \caption{Error of compensated control voltages}
        \label{fig:error_vangle}
    \end{subfigure}
    \caption{Accuracy of the algorithm for estimating normal vectors $\Gamma$. We show both the angular deviations of the estimated rows of $\Gamma$ to the true values as well as the misalignment of the compensated coordinate system.}
    \label{fig:sim_error}
\end{figure}

Next, we present our evaluation of the results. We will focus here on measurements for precision and accuracy of the estimated polytopes. Our results on running time of the algorithm and required number of line-searches can be found in Appendix~\ref{app:runtime}.
\subsection{Evaluation of \texorpdfstring{$\Gamma$}{Gamma}}
We first evaluate the results of our algorithm for computing $\Gamma$. For this, we evaluate the values of $\Gamma$ computed in scenarios $S_1$ (3x2), $S_3$ (3x3) and $S_4$ (4x4). We first evaluated in how many cases the algorithm succeeded. This means that the algorithm terminated and for each row of $\Gamma$, we found at least $N+3$ pairs of samples separated by it. 

In all our experiments, we only observed failed trials on the 4x4 device with setting $\rho=3$, where the solver could get trapped in a local optimum. This happened in approximately 10\% of trials.  
This also means that for these trials no other polytopes of the scenario could be computed and we exclude them from our experiments. Still, all our results reported for $S_4$ and $S_5$ had at least 17 successful trials for $\rho=3$.

For the successful trials, we computed the angular differences between $A_k$ and $\Gamma_k$ for all dots and trials and report their 5\%,50\% and 95\% quantiles of the angular errors in all scenarios in Figure~\ref{fig:error_angle}. 
In all cases, the error is below $0.1$ degrees. 
Further, we measured the correctness by computing the angle between the transitions of $P_0$ in the coordinate system of the compensated control voltages (i.e. the angle between the rows of $A\Gamma^{-1}$) and the standard basis. This is the most important metric, since it measures how perpendicular the transitions of $P_0$ are in this coordinate system and thus, how well the compensation works.
Again, we report their 5\%, 50\% and 95\% quantiles in all scenarios in Figure~\ref{fig:error_vangle}. The error is slightly higher but still below $0.15$ degrees.


\begin{figure}[ht]
    \centering
    \begin{subfigure}{0.7\textwidth}
        \centering
        \includegraphics[width=\textwidth]{poly_num_found}
        \caption{Fraction of target transitions found}
        \label{fig:poly_num_found}
    \end{subfigure}\\
    \begin{subfigure}{0.7\textwidth}
        \centering
        \includegraphics[width=\textwidth]{poly_radii_missed}
        \caption{Relative size of missed transitions}
        \label{fig:poly_error_radii}
    \end{subfigure}
    \caption{Success rate of the algorithm for computing $P_n$. Top: fraction of facets found. We take only facets of the polytope into account that are included in $T$ (e.g., "all one-electron transitions") and would be selected by the stopping criterion. Bottom: for missed facets, we measure the relative size compared to our stopping criterion. See the text for more details.}
    \label{fig:poly_error}
\end{figure}

\subsection{Evaluation of \texorpdfstring{$P_n$}{Pn}}
Next, we evaluate the estimated polytope transitions in the five scenarios. Our evaluation focuses on three different error conditions:
\begin{description}
    \item[False Positive] The algorithm returns a transition that does not exist on the device.
    \item[False Negative] The algorithm misses to learn a facet that does exist on the device, is included in $T$ and has largest inscribed hypersphere radius $r \geq r_{\min}$.
    \item[Unusable transition] An existing transition is found, but the estimates are so poor that using it results in a different transition.
\end{description}
For false negatives, the condition on size might sound odd at first. After all, we would hope to find all transitions in $T$ that exist in the polytope. However, our line-searches operate with limited accuracy, making it impossible to detect facets that are smaller than the accuracy of the line-search.

This is reflected in the design of the algorithm, as it only returns facets with largest inscribed hypersphere radius $\geq r_{\min}=2\delta$ (See Appendix~\ref{app:termination}). Consequently, it only makes sense to use the same threshold on the facets of the ground truth polytope to check whether they should be even considered for comparison. As a result, changing $\delta$ effectively changes the task, as not only the precision of the estimated polytope is changed, but also the ground truth we compare it to. For example in $S_4$, with $\delta=0.002$ we only expect to find approximately 80 facets on average, while at $\delta=0.001$, this number rises to 100.

To evaluate whether a facet is usable, we take the center of the computed largest inscribed hypersphere on the facets of the estimates polytope and perform a line-search through this point from inside the polytope. We then evaluate whether the facet hit in the ground truth polytope has the same label as the label assigned to the transition by the model. This mimics how a practitioner would use these facets in order to change the ground-state of the device.

In all converged trials, we did not find a single instance of a false positive or of an unusable transition. This leaves evaluation of false negatives. The fraction of facets found in all scenarios is shown in Figure~\ref{fig:poly_num_found}. In all scenarios and for all settings, more than 96\% of the transitions have been found and for $S_5$ all facets have been found in all trials.

For the missed facets with radius $r\geq r_{\min}$, we computed the relative size-difference $\frac {r-r_{\min}} {r_{\min}}$ of the radius of the inscribed hypersphere to the stopping criterion. The reported 50\% and 95\% quantiles are shown in Figure~\ref{fig:poly_error_radii}. One can see that for the scenarios $S_3$, $S_4$ and $S_5$ the missing facets have a difference in radius consistently smaller than $1.1r_{\min}$, while for $S_1$ we observed a few facets being missed with radius $r\approx 2r_{\min}$ in the case of $\rho=1$ and $\delta=0.002$. However, this number is based on a total of 7 missed facets over all estimated polytopes.

For scenarios $S_1$ and $S_2$ we were targeting specific (sequences of) transitions. We evaluated whether the tasks were solved by checking whether all required transitions were found. For $S_2$ the task was to find a single transition. This was achieved in all but 2 trials in the the setting $\rho=1,\delta=0.002$. The remaining 78 trials succeeded. For $S_1$ we had to check whether all 6 transitions were found. For this we made use of the fact that each transition appears as candidate in two polytopes: either as part of the polytope belonging to the initial state, or as transition in the polytope from the target state back to the initial state. If any of them succeeded in finding the transition, we counted the transition as found. With this strategy, we managed to find the sequence of transitions in all cases but one for $\rho=1, \delta=0.002$. In this case, we found that one of the required facets in the ground truth polytope was smaller than our cut-off radius and thus the task was not solvable for the algorithm with our chosen stopping criterion.

\section{Discussion}

In this manuscript we introduced an algorithm for determining Coulomb diamonds (polytopes) in large quantum dot arrays. 
We used regularized maximum likelihood optimization to estimate a model (Eq.~\ref{eq:model}), and tested it on several scenarios. 
In almost all of the cases, the resulting polytopes and transitions were found with high accuracy. Most importantly, our algorithm did neither return any transition that did not exist on the device, nor mislabeled any of the transitions it found. Moreover, we showed that we can construct a point on each facet such, that a linear ramp from a point inside the polytope through it results in the desired ground state transition.

Additionally, our method for proposing new line-search directions is efficient. Our algorithm required less than 15000 measurements even in devices with 16 dots with up to 100 facets of interest. This places our algorithm within a factor 10 of the minimum amount ($G\cdot\lVert T\rVert$) of measurements required to verify the learned polytope. This is in contrast to the bounds obtained for the number of line-search directions for uninformed sampling strategies (see Theorem 2 in \cite{weber2022theoretical}), which are exponential in the dimension of the number of dots and which assume perfect line-search accuracy.

To achieve this, we made several assumptions. 
First, we use the constant interaction model to describe the array of quantum dots. We also assume that the devices can realize, and measure, transitions between ground-state transitions reliably. 
This means that we can measure whether an electron moved inside the array and that the device is manufactured such, that existing electron transitions occur when the device control parameters are chosen accordingly.
For example, an existing transition adding an electron to a quantum dot requires that this dot can easily exchange electrons with a reservoir and does not need to rely on co-tunneling events. However, this is not a limitation of the algorithm, but of the device. if a theoretically existing ground state transition does not occur during our measurements, we can as well treat the transition as non-existing and we have shown that our algorithm handles missing transitions gracefully.


Our assumptions on the measurements of transitions are weaker. 
The only assumption we make here is that when linearly ramping the gate voltages between $v_\text{start}$ and $v_\text{end}$, we can reliably measure whether a transition occurred and that we can determine the transition voltages with within a precision $\delta$.
We do not require that we know \textit{which} transition occurred nor how the measurement is conducted. 
We treat $\delta$ as an external parameter: if the sensor signal is noisy and uncertain, our algorithm will miss facets of the Coulomb diamonds that are small since they can not be detected reliably anymore.

While the devices we consider are idealized, our work is an important stepping stone towards fully autonomous tuning of these transitions in spin qubit arrays. 
If we fail to find an algorithm under the idealized conditions described in this work, we can not expect to find an algorithm that performs well on a real device. 
On the other hand, if we find an algorithm that works reliably, there is hope that this algorithm can be generalized to less favourable conditions in the future, or that improvements in material design, manufacturing and sensing progressively narrow the gap. In this case, we have shown that our algorithm provides the essential tools to program a device reliably.

Still, our algorithm is not perfect. For estimating $P_0$, which is important for computing compensated control voltages, we have observed that our algorithm does not find the solution in all cases for large devices. We expect that this problem can be solved simply by re-running the algorithm with a new set of measurements. Failures here are expected, as the problem of estimating convex polytopes using a binary classification dataset is known to be a NP-hard machine-learning problem \cite{gottlieb2018learning}. Our results also showed that our strategy for obtaining additional measurements does not help to solve the issue: In most trials the algorithm managed to find the solution with the chosen initial set of randomly sampled points, or not at all. We suspect that this problem can be solved more efficiently by obtaining measurements not only of the position of the transition but also of its normal, as has been utilized by experimentalists in the literature \cite{volk2019loading}. However, adapting this approach requires a more careful modeling of the sensor signal, which is outside the scope of our work.

A limitation of our work arises due to the charge sensors themselves. In our algorithm, we do not model the effects of charge sensing on the constant interaction model, thus adding the implicit assumption that charge sensors do not interact with it. 
This is an assumption that may not be valid in practice, because charge sensors are capacitively coupled to the array and therefore change the capacitance matrices. Still, our results carry over to weakly coupled charge sensors, when their effect on the constant interaction model can be approximated by adapting \eqref{eq:const_interaction_trans} to
\[
    P_n \approx \left\{ v\in \mathbb{R}^G \mid t^TA v + t^TA^S v^S + b_{n,n+t} \leq 0,\, \forall t \in \{-1,0,1\}^{N}  \right \},
\]
where $v^S$ is the vector of gate voltages of the sensor gate electrodes and $A^S$ is the effect of the gate electrodes on the dots of the device due to capacitive coupling. In practice, gate-voltages $v^S$ are chosen to compensate for changes of the control voltages $v$ to keep the potential of the sensor dots constant. This compensation is done by choosing $v^S$ as affine linear functions of $v$, i.e., $v^S=U^Sv+b^S$. The resulting polytope as a function of $v$ using sensor compensation is described as
\[
    P^\text{Comp}_n \approx \Bigg\{ v \mid t^T\underbrace{(A+A^SU^S)}_{\tilde{A}} v + \underbrace{b_{n,n+t}+t^TA^Sb^S}_{\tilde{b}_{n,n+t}} \leq 0,\, \forall t \in \{-1,0,1\}^{N}  \Bigg \}\enspace.
\]
The polytope still abides to the general form of \eqref{eq:const_interaction_trans}, especially the connection of transition $t$ and its normal, which allows to compute the polytope $P^\text{Comp}_0$ to obtain an estimate of $\tilde{A}$ that can be used to learn $P^\text{Comp}_n$ with the adapted offset $\tilde{b}_{n,n+t}$.
Thus, our proposed algorithm can still be used to learn the transitions of the \emph{compensated} gate voltages, assuming that the sensor compensation function is the same for all polytopes of interest and the number of electrons on the sensor dot remains constant. We leave more complex interaction of sensors for future work.

More generally, future work must incorporate more general types of deviations of the device from the constant interaction model, which requires adaptations of our model and fitting process. While the work presented here is able to adapt to some deviations due to the freedom of picking parameters $b$ and $\Lambda$, deviations in the normals of the learned facets are currently not modeled well. 
This holds especially for deviations that cause the facets of transitions $t$ and $-t$ (e.g., adding and removing an electron at a dot location) to not be parallel. These deviations can be described by generalisations of the constant interaction model, in which the capacitance matrices $C^{DD}$ and $C^{DG}$ may depend on the occupation numbers $n$ of the array. These changes reflect that dots with large $n_i$ are typically somewhat larger and thus have typically larger capacitances, than a dot with small $n_i$.
In this case, using an estimate of $\Gamma$ obtained in the empty device ($P_0$) might not be a good description of the normals of the coulomb diamond of a state with large electron occupation. However, if we still assume that deviations from parallelity in the polytope $P_n$ are small, we can try to adapt $\Gamma$, while learning the polytope. This requires a careful choice of regularisation and adaptation of the stopping criteria to expect larger errors. For larger deviations, additional terms must be incorporated into the model.

Further, our work does not cover the tuning of additional gate electrodes, for example "barrier electrodes" that tune the tunnel barriers between the dots. We currently assume that practitioners can fix those parameters before using our algorithm for estimating $\Gamma$ and other polytopes $P_n$. However, these electrodes are fine-tuned after certain states have been found (e.g., the state $n=(1,1,\dots, 1)$ or are changed during device operation, and thus affect the learned polytopes due to capacitive interactions. A full incorporation of these barrier electrodes in our algorithm is therefore a challenge and a desirable next step to allow automatic tuning of quantum dot arrays.

\section*{Acknowledgements}
This project has received funding from the European Union’s Horizon 2020 research and innovation program under the Marie Sklodowksa-Curie grant agreement No. 895439 `ConQuER'.

\section*{Code and Data Repository}
The code (and link to the associated data) can be found at \url{https://github.com/Ulfgard/quantum_polytopes}.

