# Operating Envelopes under Probabilistic Electricity Demand and Solar
Generation Forecasts

Yu Yi, Gregor Verbič  The University of Sydney
NSW, Australia
{yu.yi<EMAIL_ADDRESS>
###### Abstract

The increasing penetration of distributed energy resources in low-voltage
networks is turning end-users from ‘consumers’ to ‘prosumers’. However, the
incomplete smart meter rollout and paucity of smart meter data due to the
regulatory separation between retail and network service provision makes
active distribution network management difficult. Furthermore, distribution
network operators oftentimes do not have access to real-time smart meter data,
which creates an additional challenge. For the lack of better solutions, they
use blanket rooftop solar export limits, leading to suboptimal outcomes. To
address this, we designed a conditional generative adversarial network
(CGAN)-based model to forecast household solar generation and electricity
demand, which serves as an input to chance-constrained optimal power flow used
to compute fair operating envelopes under uncertainty.

###### Index Terms:

Load forecasting, generation forecasting, operating envelopes, optimization
under uncertainty, chance-constrained optimal power flow, conditional
generative adversarial networks (CGAN).

##  I Introduction

With the increasing penetration of behind-the-meter distributed energy
resources (DER), including rooftop solar photovoltaic (PV) systems, home
battery storage, electric vehicles and flexible loads, the way of electricity
generation and consumption is changing. While DER can reduce energy
expenditure for the owners, they can also create network problems, including
reverse power flows and overvoltages. In countries with a high rooftop solar
PV penetration, low-voltage distribution networks are increasingly often at
the verge of operating limits. Yet, distribution network operators (DSO) are
ill-equipped to deal with the challenge. The issue is further exacerbated by
the regulatory separation between network service provision and electricity
retail and the partial smart meter rollout. In Australia, for example, only
retailers have access to real-time smart meter data. They can sell the data to
network companies, but only for the previous day. On top of that, the smart
meter rollout is only complete in Victoria, while in other states, most
customers still use accumulation meters. All these issues make active network
management difficult.

For the lack of better solutions, DSOs use blanket rooftop solar export
limits; that is, they limit solar exports conservatively to ensure that
network limits are not violated. However, this can lead to suboptimal
outcomes. Against this backdrop, the authors in [1, 2, 3] proposed to use
dynamic operating limits (also called operating envelopes) to limit PV
exports. However, they assume full real-time observability of end-user demand
and solar generation, which is unrealistic for the reasons explained above. To
address the above issues, we proposed a chance-constrained AC optimal power
flow (CC AC OPF) to compute operating envelopes under uncertainty [4].
However, we assumed that the probability distributions of uncertain
electricity demand and solar generation is known and given without going into
details of how they can be obtained in practice. We address the gap in this
paper by proposing a conditional generative adversarial network (CGAN)-based
model to forecast household solar generation and electricity demand.

The remainder of the paper is described as follows: Section II reviews the
literature of demand and solar generation forecasting. Section III formulates
the CC AC OPF problem for the computation of dynamic operating limits. The
forecasting methodology using the proposed conditional GAN is presented in
Section IV. Section V presents simulation results, and Section VI concludes
the paper.

##  II Literature Review

Forecasting information is the basis for the planning, management, and energy
trading in electricity networks. Both the electricity demand and solar
generation exhibit seasonal behavior. The prediction results at a certain time
slot are related to not only the historical data at the previous time slot but
also at the same time slot of the previous day or on the days in previous
weeks [5]. Moreover, short-term forecasting (seconds to minutes, hours, or
days ahead) is challenged by some unexpected issues like accidental human
behaviors or weather-related conditions [6]. The stepping variations of load
or energy output will weaken the operation economics and threaten system
stability. Hence, it is essential for network managers to use and integrate
the load and solar generation forecasts into power system optimal power flow
problems. Various forecasting approaches and research improve the techniques
to obtain the point forecast information. The approaches of which usually
contain input classification, modal regression formulation, and predictions
output[7]. A Long- and Short Time-series network (LSTNet) embedded with
Convolution Neural Network (CNN) and the Recurrent Neural Network (RNN) is
proposed in [8], which performs well on the points forecasts of both the load
and energy generation. The LSTNet model can either discover local dependency
patterns between multi-dimensional inputs or capture complex relevance in the
long period. A spatiotemporal graph dictionary learning (STGDL) model is used
in [9] to generate the behind-the-meter forecasts of load and PV power.
However, the deterministic forecasts are no longer sufficient in daily
operations due to the variability and uncertainty of the energy generation,
which increases the importance of the probabilistic forecasts.

The probabilistic forecasts are usually indicated in the form of intervals,
quantiles, or density functions. The forecasting interval is valuable for
network operators to correct the prediction errors or inconsistencies in the
day-ahead schedule [10]. The conditional forecast residuals are defined in
[11] to integrate both point forecasts and quantile regression in the
probabilistic forecast output. A constrained quantile regression averaging
(CQRA) method is studied in [12] to find the optimal weights for the
individual forecasting model to improve the performance of the combined
models. Authors in [13] propose a reinforcement learning-based framework to
choose the probabilistic forecasting model which has the highest performance
dynamically. Probabilistic forecasts take the probability distributions to
account for the variation information. Sometimes, in the process of solving
stochastic optimization models, it is necessary to generate random scenarios
with the given probability distribution of uncertainties. Then the GAN model
is introduced in the probabilistic forecasts. Authors in [14] use the GAN-
based model to generate long-term stochastic scenarios for renewable energy.
To accurately capture the short-term variations of loads and energy
generation, the conditional GAN (CGAN) model was formulated in [15], [16] to
output the probabilistic forecasts.

Inspired by the work above, we establish a CGAN-based framework to produce the
probabilistic load and solar generation forecasts across different time
periods. Then, the forecasting residuals are applied into the CC AC OPF model
to generate the operating envelops.

Figure 1: The proposed framework for the computation of operating envelopes
under uncertainty. Blocks T1-T3 produce probabilistic forecasts of electricity
demand and solar generation used in the chance-constrained optimal power flow
in block T4.

##  III Operating Envelopes under Uncertainty

###  III-A Framework Overview

Fig. 1 provides an overview of the computation of operating envelopes under
uncertainty across three time horizons. A point forecasting model is trained
first using a historical dataset (T1) and related factors such as day and
season type to capture features relevant for the long-term forecast. Next, the
historical dataset (T2) is used for testing and to compute forecasting
residuals. After that, the historical dataset (T2) and related factors (day
type, season type), point forecasts, forecasting residuals and random noise
are fed into the conditional forecasting residual model. Thirdly, a new
historical dataset (T3) is applied to the trained model in the last two steps.
The point forecasts and point forecasting residuals constitute the
probabilistic forecasts we need. Thus, the probabilistic forecasts (mean and
standard deviation) are used to input the uncertain variables in the CC AC OPF
model. Finally, we approximate the CC AC OPF model as a deterministic AC OPF
model and produce the operating envelopes for individual prosumers.

###  III-B Chance-constrained AC OPF Formulation

Since household load and solar generation are difficult to forecast
accurately, random forecast errors can adversely affect the computation of
operating envelopes. To this end, we proposed a CC AC OPF model in our
previous work [4] to compute nodal exporting limitations under uncertainty. In
our approach, the DSO computes nodal export limits using day-ahead
probabilistic forecasts of solar generation and customer load and submits them
to the customers as operating envelopes that serve as constraints in their
optimisation.

We consider a radial distribution network with a set of nodes
$i\in\mathcal{N}$ and a set of lines $(i,j)\in\mathcal{L}$. For the sake of
notational simplicity, we assume that one prosumer is connected to each node.
The computation of operating envelopes is formulated as a chance-constrained
OPF aiming to maximise network hosting capacity. We model the stochastic
household demand $\tilde{p}^{\mathrm{d}}_{i}$ and solar generation
$\tilde{p}^{\mathrm{pv}}_{i}$ as the sum of the point forecast part
$\widehat{p}_{i}^{\mathrm{d},\mu}$ and $\widehat{p}_{i}^{\mathrm{pv},\mu}$ and
forecasting residuals part $\varepsilon^{\mathrm{d}}_{i}$ and
$\varepsilon^{\mathrm{pv}}_{i}$, respectively. The methodology to obtain the
points forecasts and forecasting residuals is described in Section IV. The
chance constraints are introduced to guarantee that network constraints will
not be violated when the inputs are subject to uncertainty. To ensure fair
distribution of solar PV curtailment, we maximise the minimum export limit
across all prosumers in a time slot $t\in\mathcal{T}$111The max-min
formulation can be reformulated by introducing one auxiliary variable and a
set of linear constraints for each $t\in\mathcal{T}$.. We use a second-order
cone (SOC) branch-flow OPF model, which results in the following probabilistic
MISOCP OPF problem:

$\displaystyle\underset{x^{\text{p}}\in\mathcal{X}^{\text{p}},\>x^{\text{n}}\in\mathcal{X}^{\text{n}}}{\text{max}}\sum_{t\in\mathcal{T}}\sum_{i\in\mathcal{N}}\underset{i}{\text{min}}\
p^{\mathrm{exp}}_{i,t},$ (1a) s.t.
$\displaystyle\tilde{p}^{\text{pv}}_{i,t}\\!\\!-\\!\\!\tilde{p}^{\mathrm{d}}_{i,t}\\!\\!-\\!\\!\tilde{p}^{\mathrm{b}}_{i,t}\\!\\!=\\!\\!\sum_{jk\in\mathcal{L}}\\!\\!\left(\tilde{p}_{jk,t}\\!+\\!\tilde{\ell}_{jk,t}r_{jk}\right)\\!-\\!\\!\sum_{ij\in\mathcal{L}}\\!\\!\tilde{p}_{ij,t},\
i\\!\in\mathcal{N},$ (1b)
$\displaystyle\tilde{p}^{\text{pv}}_{i,t}-\tilde{p}^{\mathrm{d}}_{i,t}-\tilde{p}^{\mathrm{b}}_{i,t}\leq\overline{p}^{\mathrm{exp}}_{i,t},\
i\in\mathcal{N},$ (1c)
$\displaystyle\tilde{q}^{\mathrm{pv}}_{i,t}\\!-\tilde{q}^{\mathrm{d}}_{i,t}\\!=\\!\sum_{jk\in\mathcal{L}}\\!\left(\tilde{q}_{jk,t}\\!+\\!\tilde{\ell}_{jk,t}x_{jk}\\!\right)-\\!\sum_{ij\in\mathcal{L}}\\!\tilde{q}_{ij,t},\
i\in\\!\mathcal{N},$ (1d) $\displaystyle p^{\mathrm{exp}}_{i,t}\leq
p^{\mathrm{exp}},\ i\in\mathcal{N},$ (1e)
$\displaystyle\tilde{v}_{i,t}\\!-\\!\tilde{v}_{j,t}\\!=\\!2\left(r_{ij}\tilde{p}_{ij,t}\\!+\\!x_{ij}\tilde{q}_{ij,t}\right)\\!-\\!\left|z_{ij}\right|^{2}\\!{\tilde{\ell}}_{ij,t},\>(i,j)\in\mathcal{L},$
(1f) $\displaystyle\left\|\begin{array}[]{c}2\tilde{p}_{ij,t}\\\
2\tilde{q}_{ij,t}\\\
{\tilde{\ell}}_{ij,t}-\tilde{v}{i,t}\end{array}\right\|_{2}\leq{\tilde{\ell}}_{ij,t}+\tilde{v}_{i,t},\
(i,j)\in\mathcal{L},$ (1j)
$\displaystyle\tilde{p}_{ij,t}^{2}+\tilde{q}_{ij,t}^{2}\leq\overline{s}_{ij}^{2},\
(i,j)\in\mathcal{L},$ (1k) $\displaystyle
0\leq\tilde{p}^{\mathrm{pv}}_{t}\leq\overline{p}^{\mathrm{pv}}_{t},$ (1l)
$\displaystyle\underline{p}^{\mathrm{b}}\leq
b_{t}\tilde{p}^{\mathrm{b},+}_{t}+(1-b_{t})\tilde{p}^{\mathrm{b},-}_{t}\leq\overline{p}^{\mathrm{b}},$
(1m)
$\displaystyle{\tilde{e}^{\mathrm{b}}_{t}-\tilde{e}^{\mathrm{b}}_{t-1}}=\left(\eta
b_{t}\tilde{p}^{\mathrm{b},+}_{t}-\frac{(1-b_{t})}{\eta}\tilde{p}^{\mathrm{b},-}_{t}\right)\Delta{t},$
(1n)
$\displaystyle\underline{e}^{\mathrm{b}}\leq\tilde{e}^{\mathrm{b}}_{t}\leq\overline{e}^{\mathrm{b}},$
(1o) $\displaystyle b_{t}\in\left\\{0,1\right\\},$ (1p)
$\displaystyle\mathbb{P}\left(\widetilde{v}_{i}\leq\overline{v}_{i}\right)\geq
1-\xi^{\text{v}},$ (1q)
$\displaystyle\mathbb{P}\left(\widetilde{v}_{i}\geq\underline{v}_{i}\right)\geq
1-\xi^{\text{v}},$ (1r)
$\displaystyle\mathbb{P}\left(\widetilde{p}_{ij}^{2}+\widetilde{q}_{ij}^{2}\leq\overline{s}_{ij}^{2}\right)\geq
1-\xi^{\ell},$ (1s)
$\displaystyle\mathbb{P}\left(\widetilde{p}_{i,t}^{\mathrm{b}}\leq\overline{p}^{\mathrm{b}}\right)\geq
1-\xi^{\text{p}},$ (1t)
$\displaystyle\mathbb{P}\left(\widetilde{p}_{i,t}^{\mathrm{b}}\geq\underline{p}^{\mathrm{b}}\right)\geq
1-\xi^{\text{p}}.$ (1u)

where
$\mathcal{X}^{\text{n}}=\left\\{p^{\mathrm{exp}}_{i,t},\tilde{v}_{i,t},\tilde{u}_{i,t},\tilde{\ell}_{i,t},\tilde{p}_{ij,t},\tilde{q}_{ij,t}\right\\}$
for each $t\in\mathcal{T}$, $(i,j)\in\mathcal{L}$ and $i\in\mathcal{N}$ is the
vector of network decision variables;
$\mathcal{X}^{\text{p}}=\left\\{\tilde{p}^{\mathrm{pv}}_{i},\tilde{p}^{\mathrm{b}}_{i,t},\tilde{e}^{\mathrm{b}}_{t},b_{t}\right\\}$
for each $t\in\mathcal{T}$, $(i,j)\in\mathcal{L}$ and $i\in\mathcal{N}$ for
each $t\in\mathcal{T}$ is the vector of prosumer decision variables;
$\tilde{\bullet}$ represents the stochastic variables. The constraints include
active and reactive power balance for each node (1b)-(1d); maximum export
limit (1e); voltage drop (1f), where $r_{ij}$ and $x_{ij}$ are resistance and
reactance of the line, respectively, and $z_{ij}=r_{ij}+jx_{ij}$ is the line
impedance; (1j) is the SOC relaxation of the line power flow; (1k) is the line
power flow limit; (1l) - (1p) are prosumers solar PV curtailment and battery
constraints, respectively. For reactive powers $\tilde{q}^{\mathrm{pv}}$ and
$\tilde{q}^{\mathrm{d}}$ we assume a constant power factor
$\text{cos}\varphi=0.9$. (1q)-(1u) are stochastic constraints ensuring the
probability of a stochastic variable $\widetilde{x}$ remaining within limits
is $1-\xi^{\text{x}}$, where (1q)-(1r) are the voltage magnitude limits, (1s)
is the apparent power flow limit, and (1t)-(1u) ensures battery storage power
limits. The CC AC OPF model is then reformulated into the deterministic AC OPF
model referred to [17], which can be solved efficiently.

##  IV Demand and Generation Forecasting

Solar generation and load forecasting are affected by various factors such as
radiation, cloud coverage, temperature, user behavior, etc., and their
interaction is intricate. No forecasting approach can capture all the
influencing factors and generate absolute perfect results. Hence, the
forecasting residuals are introduced to model the characteristics not captured
in point forecasting. The point forecasts
$\widehat{y_{t}}=f\left(\mathbf{X}_{t}\right)$ and forecasting residuals
$\varepsilon_{t}$ are respectively considered as deterministic and uncertain
part to form the probabilistic forecasts $y_{t}$. The input factor
$\mathbf{X}_{t}$ contains the historical data, day type and season type. The
forecasting residuals denotes the errors between real and predicted values. In
this paper, we use the long short-term memory network (LSTN) model to train
and test for the point forecasts, and a CGAN model for forecasting the
residuals.

The whole forecasting process is conducted along three time horizons, with
each using a different dataset. In the first step, the LSTN
$f\left(\cdot\right)$ is trained using dataset T1. In the second step, we test
the model $f\left(\mathbf{X}_{T_{1}}\right)$ with dataset T2 and calculate the
forecast residuals
$\varepsilon_{T_{2}}=y_{T_{2}}-f\left(\mathbf{X}_{T_{2}}\right)$. The CGAN
model $g\left(\mathbf{C}_{t},z\right)$ is also trained, where $\mathbf{C}_{t}$
denotes condition factors (e.g. forecasts $y_{T_{2}}$, historical data T2, and
day type) and $z$ denotes the random noise. The generator of the CGAN model
can sample different residuals depending on the random noise $z$. In the third
step, the point forecasting model $f\left(\mathbf{X}_{T_{3}}\right)$ and the
residual forecasting model $g\left(\mathbf{C}_{T_{3}},z\right)$ use dataset T3
to generate scenario forecasts. The performance of the uncertainty and
variation of the forecasts are usually evaluated the mean square error (MSE)
or the mean absolute error (MAE). In this paper, the forecasts are evaluated
in terms of the continuous ranked probability (CRPS) and the pinball loss
(PL), which are more appropriate in our setting, as explained in Section IV-B.

Figure 2: Structure of the conditional generative adversarial network (CGAN).

###  IV-A Conditional Generative Adversarial Network (CGAN)

The general task of the Generative Adversarial Network (GAN) model is to train
two neural networks, called generator $G$ and discriminator $D$ in an
adversarial learning mechanism. The CGAN model is the improvement of GAN model
to capture variations comprehensively, where the condition $C$ is added as the
inputs. The structure of the CGAN is shown in Fig. 2, where the generator
$g\left(\mathbf{C}_{T_{3}},z\right)$ generates a large number of forecasting
residual scenarios $S_{g}$ with different random noise $z$ until they are
infinitely close to the real dataset $S_{r}$ (the real forecasting residuals).
The noise $z$ usually follows a pre-set distribution $z\sim\mathbb{P}_{z}$,
such as Gaussian or uniform distribution. The discriminator compares the
generated sample $S_{g}$ to the real ones to determine whether it should be
accepted or rejected. The required parameters for generator and discriminators
are $\theta_{G}$ and $\theta_{D}$, respectively. The sample vectors input to
discriminator are represented as $S_{g}=G(z,\theta_{G})$. The discriminator
$D(\cdot)$ distinguish whether samples are generated or real by estimating the
values of the probability
$p_{i}=D\left(G\left(z;\theta_{G}\right)\mid\mathbf{C},\theta_{D}\right)$.

The cross-entropy-based losses of generators and discriminators in the form of
expectation for sample vector $S_{g}$ is defined as:

$\displaystyle
L_{G}=\mathbb{E}_{Z}\left[\log\left(1-D\left(G\left(z;\theta_{G}\right)\mid\mathbf{C},\theta_{D}\right)\right)\right],$
(2a)
$\displaystyle\begin{split}L_{D}=-\mathbb{E}_{S_{r}}\left[\log\left(D\left(s_{r}\mid\mathbf{C};\theta_{D}\right)\right)\right]\\\
-\mathbb{E}_{Z}\left[\log\left(1-D\left(s_{g}\mid\mathbf{C};\theta_{D}\right)\right)\right],\end{split}$
(2b)

where $\mathbb{E}_{Z}$ and $\mathbb{E}_{S_{r}}$ represent generated sample
inputs and the expected values among real sample inputs, receptively. During
training, the generator and discriminator are optimized simultaneously. The
generated samples get more similar to the real ones with the iteration number
increasing.

The adversarial mechanism between two neural networks $G$ and $D$ is
formulated as a min-max optimization problem with the following objective:

$\begin{split}\min_{\theta_{G}}\max_{\theta_{D}}&\mathbb{E}_{S_{r}}\left[\log\left(D\left(s_{r}\mid\mathbf{C};\theta_{D}\right)\right)\right]+\\\
&\mathbb{E}_{Z}\left[\log\left(1-D\left(G\left(z;\theta_{G}\right)\mid\mathbf{C},\theta_{D}\right)\right)\right].\end{split}$
(3)

The authors in [16] improved the GAN model to the Wasserstein GAN (WGAN) by
adding a gradient penalty to accelerate the computation while stabilizing the
training process and ensuring that the convergence rates of the generator and
discriminator are synchronized. Compared to GAN, the forecasting model in WGAN
uses a linearized activation function in the last layer of the neutral network
$D$, so that the weights of network $D$ are limited into a fixed range after
each training. Finally, the logarithmic functions in the generator and
discriminator loss functions are replaced by linear ones since gradients of
logarithmic functions are negligible. Hence, the simplified loss functions can
be rewritten as:

$\displaystyle
L_{G}=-\mathbb{E}_{Z}\left[D\left(G\left(z;\theta_{G}\right)\mid\mathbf{C},\theta_{D}\right)\right],$
(4a) $\displaystyle
L_{D}=-\mathbb{E}_{S_{r}}\left[\left(D\left(s_{r}\mid\mathbf{C},\theta_{D}\right)\right)\right]+\mathbb{E}_{Z}\left[D\left(s_{g}\mid\mathbf{C},\theta_{D}\right)\right]$
(4b)

The gradient penalty ($GP$) is added to the discriminator loss function to
make sure the distinguishing ability of discriminator under simplification,
which is defined as:

$GP=\lambda\mathbb{E}_{S}\left[\|\nabla D(\hat{s})\|_{2}-1\right]^{2},$ (5)

where $\lambda$ represents the hyper-parameter that determines the weight of
$GP$; $\hat{s}=\rho S_{g}+\left(1-\rho\right)S_{r}$ represents the average of
the generated samples $S_{g}$ and real data samples $S_{r}$, where $\rho$ is a
random variable following uniform distribution $U(0,1)$; $\mathbb{E}_{S}$
represents the expected value among all $\hat{s}$. With the gradient penalty,
the discriminator loss function is reformulated as:

$L_{D}=-\mathbb{E}_{S_{r}}\left[\left(D\left(s_{r}\mid\mathbf{C},\theta_{D}\right)\right)\right]+\mathbb{E}_{Z}\left[D\left(s_{g}\mid\mathbf{C},\theta_{D}\right)\right]+GP.$
(6)

###  IV-B Evaluation Criteria

To evaluate stochastic forecasts, the continuous ranked probability (CRPS) and
pinball loss (PL) are introduced in the probabilistic domain. The calculation
of CRPS is formulated as:

$\operatorname{CRPS}\left(F,y_{t}\right)=\int_{-\infty}^{\infty}\left(F\left(y_{t}^{\prime}\right)-1\left(y_{t}^{\prime}-y_{t}\right)\right)^{2}dy_{t}^{\prime},$
(7)

where $F\left(y_{t}^{\prime}\right)$ denotes the cumulative distribution
function of the predicted variable at time $t$;
$1\left(y_{t}^{\prime}-y_{t}\right)$ is the term to compare the values between
$y_{t}^{\prime}$ and $y_{t}$, and the value of it equals to 1 if
$y_{t}^{\prime}\textgreater y_{t}$ and equals to 0 otherwise. The forecasts
are better with smaller CRPS.

The PL is capable of calculating both positive and negative errors for
different quantile compared to MAE, which is defined as:

$\operatorname{PL}\left(\hat{y}_{t,q},y_{t},q\right)=\begin{cases}(1-q)\left(\hat{y}_{t,q}-y_{t}\right)&y_{t}<\hat{y}_{t,q}\\\
q\left(y_{t}-\hat{y}_{t,q}\right)&y_{t}\geq\hat{y}_{t,q}\end{cases},$ (8)

where $q$ denotes the quantile; $\hat{y}_{t,q}$ denotes the predicted quantile
at time $t$.

##  V Case study

###  V-A Data Implementation

The condition inputs $\mathbf{C}$ include historical data of six days before
the forecasting day (1, 2, 3, 5, 7, 14 and 21 days before), the day type (day
of the week $\textit{W}_{d}$) and season type (day of the season
$\textit{S}_{d}$). Note the spring in the Southern Hemisphere begins in
September. Hence, the condition input vector is

$\mathbf{C}=\left[\mathbf{L}_{d-1},\mathbf{L}_{d-2},\mathbf{L}_{d-3},\mathbf{L}_{d-7},\mathbf{L}_{d-14},\mathbf{L}_{d-21},W_{d},S_{d}\right],$
(9)

where $\mathbf{L}_{d}-i$ represents the daily real data of $i$ days before the
forecasting day $d$. We assume 30-min dispatch, resulting in the dimension of
the condition input vector for one day to be 290 ($6\times 2\times 24+2=290$).

We use Ausgrid’s Solar Home Electricity dataset222www.ausgrid.com.au
consisting of solar generation and household electricity demand from
2010-2012. Dataset T1, between 1 July 2010 to 30 June 2011, is used to train
the LSTN model, while dataset T2, from 1 July 2011 to 30 June 2012 (except
7th, 14th and 28th of each month, 36 days in total) is used to test the point
forecasts and to train the CGAN model. The testing dataset T3, used to
evaluate the performance of the forecasting framework, consists of the days
that were excluded from the training dataset T2 (7th, 14th and 28th day of
each month, 36 days in total).

###  V-B Parameter Setup

We use a three-layer network for both generator and discriminator, containing
an input layer, a hidden layer, and an output layer. The number of neurons in
the input layer of generator is $802=290+512$. The hidden layer has 256
neurons while the output layer generate 48 scenarios. The dimension of random
noise $z$ is 512, which is larger than the condition input $\mathbf{C}$. For
the discriminator, there are 338 neurons in the input layer (48 for the number
of samples, 290 for the condition input vector). The number of neurons in the
hidden layer is 128 and the output layer has one neuron. The activation
functions of the hidden and the output layer are ReLU and Sigmoid,
respectively. The training strategy for the CGAN model is to train the
generator once and the discriminator five times, which balances the coverage
speed of two neurons network. The simulations were produced using the
optimisation modelling toolbox YALMIP with CPLEX solver in MATLAB® 2020 on a
64-bit laptop with a 2.60 GHz CPU and 16GB RAM.

###  V-C Neural Network training

Fig. 3 illustrates the learning process showing the loss functions of the
generator and discriminator over 20,000 iterations. Observe that the training
for both converges after about 7,000 iterations. The convergence value of the
discriminators is near 0, which means the sample scenarios generated by
generators is well approached to the real residuals. The CRPS are 27.38 and
48.62 for load forecasts and solar generation forecasts, respectively, and the
PL are 15.97 and 24.11 for load and solar generation forecasts, respectively,
which shows that the generated scenarios capture the variation characteristics
well for household load.

Figure 3: Losses of generator and discriminator

###  V-D Forecasting Results

Fig. 4 shows the real data, point forecasts and scenarios generated from
generators of uncertain variables (loads and solar generation) at the day of
7th, 14th and 28th on July and December, respectively. We select the days of
different to explore the regularity of predictions and forecasting residuals.
Observe that the variation of solar generation forecasting residuals is larger
than the one of the electricity demand. This is because of the natural
variability of solar insolation, which is affected by weather conditions, in
particular the intermittent cloud coverage. In contrast, the impact of
customer behavior on electricity demand forecast is less significant. Observe
also that the forecasting residuals are conditional on the time periods (e.g.
the variations of electricity demand forecasting residuals is more significant
at night while the solar generation forecasting residuals vary at daytime.)

Figure 4: Scenarios generated with CGAN model.

The Gaussian probability distributions of electricity demand and solar
generation forecasts at different time periods are shown in Figs. 5 and 6,
respectively. Observe that the electricity demand forecast is more volatile in
the afternoon and evening due to the increased customer activity when they
return from work, and less volatile in the morning and in the middle of the
day. By contrast, the solar generation forecast is more volatile during the
day due to the weather conditions, in particular the intermittent cloud
coverage.

Figure 5: Probability distributions of electricity demand forecasts at
different time periods.
Figure 6: Probability distribution of solar generation forecasts at different
time periods

###  V-E Operating Envelopes under Probabilistic Forecasts

We use a distribution network test system with 25 prosumers shown in Fig. 7
with HEMS, PV and battery storage systems (BSS) communicating with the DSO
through a smart meter. The Time-of-use (ToU) electricity tariff and feed-in-
tariff (FiT) are shown in Tab. I. The voltage limitation is $\pm 6\%$ of the
normalized values. Each prosumer is equipped with a
$6\text{\,}\mathrm{kW}$/$14\text{\,}\mathrm{kW}\text{\,}\mathrm{h}$ PV-battery
system. The battery charging/discharging limit is
$-3.5\text{\,}\mathrm{kW}$/$3.5\text{\,}\mathrm{kW}$, and the energy capacity
is limited between $4\text{\,}\mathrm{kW}\text{\,}\mathrm{h}$ and
$10\text{\,}\mathrm{kW}\text{\,}\mathrm{h}$. We fit a Gaussian distribution
using the generated scenarios to obtain standard deviations of the electricity
demand and solar generation in each time period. The parameters of the
distribution are used to approximate the stochastic constraints in (1).

Fig. 8 shows the weekly operating envelopes and mean nodal power flows with
different days of solar generation and demand forecasts. The nodal exported
power is limited with regard to the forecasting residuals, ensuring the
operations by DSO under uncertainties. The nodal exporting limitations are
obvious on the days with high solar generation, but the exporting limitations
relaxed on the days with low solar generation. The user-side solar PV
curtailment increases at the middle of the day as the PV system provide more
electricity.

TABLE I: Time-of-use (ToU) electricity tariff and feed-in-tariff (FiT). Time slot h) | 0-7 | 7-15 | 15-21 | 21-22 | 22-24
---|---|---|---|---|---
ToU(c/kWh) | 15.96 | 25.96 | 57.76 | 25.96 | 15.96
FiT(c/kWh) | 9 | 9 | 9 | 9 | 9

Figure 7: Distribution network topology.
Figure 8: Weekly operating envelopes and average nodal power flows calculated
by CC AF OPF model with probabilistic forecasts.

##  VI Conclusion

This paper has proposed a novel methodology for computing operating envelopes
under probabilistic electricity demand and solar generation forecasts. The
probability distributions of electricity demand and solar generation is
captured by a CGAN-based model, which is conditional on the historical
dataset, point forecasts and time horizons. The parameters of the Gaussian
distribution (mean and standard deviations) are calculated using the generated
residual scenarios, which is applied in the CC AC OPF problem to limit the
nodal operating envelopes. The effectiveness and performance of our proposed
framework is evaluated by the case studies.

## References

  * [1] K. Petrou, M. Z. Liu, A. T. Procopiou, L. F. Ochoa, J. Theunissen, and J. Harding, “Operating envelopes for prosumers in LV networks: A weighted proportional fairness approach,” in _IEEE PES Innovative Smart Grid Technologies Conference Europe_ , 2020.
  * [2] K. Petrou, A. T. Procopiou, L. Gutierrez-Lagos, M. Z. Liu, L. F. Ochoa, T. Langstaff, and J. Theunissen, “Ensuring distribution network integrity using dynamic operating limits for prosumers,” _IEEE Trans. Smart Grid._ , 2021.
  * [3] L. Blackhall, “On the calculation and use of dynamic operating envelopes, Evolve Project M4 Knowledge Sharing Report,” in _Technology report_ , 2020\.
  * [4] Y. Yi and G. Verbič, “Fair operating envelopes under uncertainty using chance constrained optimal power flow,” in _2022 Power Systems Computation Conference (PSCC)_ , 2022.
  * [5] H. S. Hippert, C. E. Pedreira, and R. C. Souza, “Neural networks for short-term load forecasting: A review and evaluation,” _IEEE Transactions on power systems_ , vol. 16, no. 1, pp. 44–55, 2001.
  * [6] D. Chaturvedi and I. Isha, “Solar power forecasting: A review,” _International Journal of Computer Applications_ , vol. 145, no. 6, pp. 28–50, 2016.
  * [7] M. Negnevitsky, P. Mandal, and A. K. Srivastava, “An overview of forecasting problems and techniques in power systems,” in _2009 IEEE Power & Energy Society General Meeting_. IEEE, 2009, pp. 1–4.
  * [8] G. Lai, W.-C. Chang, Y. Yang, and H. Liu, “Modeling long-and short-term temporal patterns with deep neural networks,” in _The 41st International ACM SIGIR Conference on Research & Development in Information Retrieval_, 2018, pp. 95–104.
  * [9] M. Khodayar, G. Liu, J. Wang, O. Kaynak, and M. E. Khodayar, “Spatiotemporal behind-the-meter load and pv power forecasting via deep graph dictionary learning,” _IEEE transactions on neural networks and learning systems_ , 2020\.
  * [10] T. Hong and S. Fan, “Probabilistic electric load forecasting: A tutorial review,” _International Journal of Forecasting_ , vol. 32, no. 3, pp. 914–938, 2016.
  * [11] Y. Wang, Q. Chen, N. Zhang, and Y. Wang, “Conditional residual modeling for probabilistic load forecasting,” _IEEE Transactions on Power Systems_ , vol. 33, no. 6, pp. 7327–7330, 2018.
  * [12] Y. Wang, N. Zhang, Y. Tan, T. Hong, D. S. Kirschen, and C. Kang, “Combining probabilistic load forecasts,” _IEEE Transactions on Smart Grid_ , vol. 10, no. 4, pp. 3664–3674, 2018.
  * [13] C. Feng, M. Sun, and J. Zhang, “Reinforced deterministic and probabilistic load forecasting via $q$-learning dynamic model selection,” _IEEE Transactions on Smart Grid_ , vol. 11, no. 2, pp. 1377–1386, 2019.
  * [14] C. Jiang, Y. Mao, Y. Chai, M. Yu, and S. Tao, “Scenario generation for wind power using improved generative adversarial networks,” _IEEE Access_ , vol. 6, pp. 62 193–62 203, 2018.
  * [15] Y. Chen, X. Wang, and B. Zhang, “An unsupervised deep learning approach for scenario forecasts,” in _2018 Power Systems Computation Conference (PSCC)_. IEEE, 2018, pp. 1–7.
  * [16] Y. Wang, G. Hug, Z. Liu, and N. Zhang, “Modeling load forecast uncertainty using generative adversarial networks,” _Electric Power Systems Research_ , vol. 189, p. 106732, 2020.
  * [17] M. Lubin, Y. Dvorkin, and L. Roald, “Chance constraints for improving the security of ac optimal power flow,” _IEEE Transactions on Power Systems_ , vol. 34, no. 3, pp. 1908–1917, 2019.


# Readout using Resonant Tunneling in Silicon Spin Qubits

Tetsufumi Tanamoto  Department of Information and Electronic Engineering,
Teikyo University, Toyosatodai, Utsunomiya 320-8511, Japan
<EMAIL_ADDRESS>Keiji Ono  Advanced device laboratory, RIKEN,
Wako-shi, Saitama 351-0198, Japan

###### Abstract

Spin qubit systems are one of the promising candidates for quantum computing.
The quantum dot (QD) arrays are intensively investigated by many researchers.
Because the energy-difference between the up-spin and down-spin states is very
small, the detection of the qubit state is of prime importance in this field.
Moreover, many wires are required to control qubit systems. Therefore, the
integration of qubits and wires is also an important issue. In this study, the
measurement process of QD arrays is theoretically investigated using resonant
tunneling, controlled by a conventional transistor. It is shown that the
number of possible measurements during coherence time can exceed a hundred
under the backaction of the measurements owing to the nonlinear
characteristics of resonant tunneling. It is also discussed to read out the
measurement results by the conventional transistor.

##  I Introduction

Silicon quantum computers are among the most popular topics in the fields of
quantum physics and engineering. The legacy of integrated circuit technologies
is expected to be inherited by silicon qubit systems. Currently, one-
dimensional qubit arrays have been successfully fabricated Hensgens ; Mills ;
Fedele ; Veldhorst1 ; Ha ; Vandersypen . Silicon qubits use narrowly confined
carriers with their interactions governed by short-range spin–spin
interactions. As the distance between two qubits increases, exchange coupling
decreases exponentially. Accordingly, the control of distant qubits becomes
much more difficult than that of the nearest neighboring qubits. Moreover, the
number of wires required to manage a qubit is approximately greater than five
Yoneda . Therefore, a sufficiently small qubit is required. Because the wiring
structure at the bottom layer should have the same dimension as the qubit, the
tight structure is inevitable Veldhorst . In order to connect the distant
qubits, additional mechanism is necessary such as the shuttling qubits Noiri1
; Utsugi , or the mediated electrons Marcus ; Croot . Shuttling qubit systems
require overhead to precisely control electron shuttling. Therefore, an extra
area is required compared with mediated electrons.

In Ref.TanaAIP , we have proposed a common gate structure using mediated
electrons which are directly embedded between conventional transistor
channels. However, since most investigations thus far have been conducted in
QD systems Hensgens ; Mills ; Fedele ; Veldhorst1 ; Ha , a mediated electron
system in QD array should be considered presently. Not all QDs in QD systems
have to be qubits; some of them can act as measurement devices. In Ref.
TanaPRB , we have theoretically investigated the effect of the measurement of
QD system by using the Green’s function method Jauho , and we proved that the
energy-difference of two QDs, which corresponds to the qubit energy-difference
under magnetic fields, becomes detectable for appropriate coupling between the
QDs and the electrodes. However, in Ref. TanaPRB , we only considered a small-
bias region with a fixed-value resistor, while Zeeman splitting in channel QDs
is ignored. In this study, we investigate a side-QD system including Zeeman
splitting of all QDs. Because there are many literature regarding the quantum
operations in several QDs Fei ; Medford ; Burkard ; Kandel ; Noiri , our
emphasis is on the measurement properties.

Generally, QD with electrodes exhibits resonant tunneling in a wide bias
region. The specific features of resonant tunneling include a nonlinear
current with coherent tunneling Capasso ; Fujita ; Medford ; Kim2 , and the
transport properties have been discussed in literature Reed ; Jiang . In this
study, we theoretically show that, by focusing on the nonlinear
characteristics of resonant current, efficient measurements against the
backaction of the measurements can be obtained. Specifically, we propose a
system to measure the energy difference between two qubits using the nonlinear
change in current depending on the spin direction. We illustrate that the
number of the detections of the spin states during decoherence exceeds a
hundred, which can lead to the possibility of surface coding in the
measurement phase Fowler .

As aforementioned, an advantage of silicon qubits is that they can be directly
integrated into conventional complementary metal-oxide-semiconductor (CMOS)
circuits on the same chip. Therefore, some parts on the chip should connect
qubits and CMOS transistors using local wires. In this study, we investigate
the transport properties of the aforementioned qubit system directly connected
to the conventional transistors BSIM . A transistor-connected qubit system
offers several advantages to control the qubit system. In general, the
variations in the size of QDs are inevitable, resulting in the variations of
energy levels of the QDs. The transistor can adjust the difference of the
energy levels of the qubits. Moreover, as the switching speed of a
conventional CMOS is in the order of ps Jena , the spin qubits can be quickly
controlled by the CMOS TanaCMOS , compared with the previous methods of
separating the measurement apparatus from the cryogenic refrigerator using
long wires ISSCC1 ; ISSCC2 . As semiconductor technologies continue to advance
daily and the era of 3D circuits is on the rise Chung ; Tsutsui ; Kim ;
Ryckaert ; Tanaka ; Lee ; Lue , active usage of the advances in state-of-the-
art technologies is necessary. Hence, the implementation of the proposed qubit
system has been discussed as a near-future structure.

The rest of this study is organized as follows. In Section II, our proposed
model is explained in detail with the formalism using the standard Green
function method. In Section III, the numerical results of our method are
presented. In Section IV, discussions regarding our results are provided.
Section V summarizes and concludes this study. In appendix, additional
explanations are presented including a possible qubit structure using the
advanced transistor system.

Figure 1: Quantum dot (QD) system considered in this study. (a) Array of
coupled qubit system connected to transistors. The electric potentials of the
qubits (yellow circles) are controlled by a common gate. The measurement
current flows from the transistors to the channel-QD (green ellipse). (b)
Measurement unit of two qubits using three QDs. $\Gamma_{\alpha s}$ denotes
the tunneling couplings between the electrodes ($\alpha=L,R$, $s=\pm 1/2$),
and $W_{ij}$s are the coupling strengths between the QDs. In the measurement
mode, the current interacts with the upper energy levels of the qubit-QD
(right-side figure). $\Delta$ and $U$ denote Zeeman splitting energy and on-
site Coulomb energy, respectively. $E_{i}$ denotes an energy-level of the
down-spin ($i=1,2,3$). (c) Depending on the position of the energy level of
the center-QD, two modes are available: interaction ($V_{D}=0$) and
measurement ($V_{D}\neq 0$) modes. The black and blue circles denote the
electrons of center-QD and qubits, respectively. The right figure represents
one of the four measurement mode types (Fig.2). (d) Typical current-voltage
characteristics of resonant tunneling. Owing to a finite magnetic field, the
bottom of the Fermi energy of the up-spin current (down-spin current) is
shifted downwards (upward). With an increase in voltage, the current increases
(point “A”), and when the bottom of the Fermi level of the drain is above the
resonant level, the current begins to decrease (point “B”).

##  II Model

Figure 1(a) shows the proposed qubit array connected to CMOS transistors.
Qubit-QDs and channel QDs are coupled by tunneling oxides such as SiO2. The
channel QDs are connected to the transistors. Figure 1(b) illustrates the
concrete model estimated in this study, which consists of two QDs designated
qubits and a channel QD. In the qubit-QD, an electron exists at the lowest
energy level and works as a spin qubit, as proposed in Ref. TanaAIP . The
next-higher energy level of the qubit-QD is placed over the on-site Coulomb
interaction ($U$) to form a singlet state ($S$), as shown in Fig.1(c).
Electrons of the center-QD enter the qubit-QDs only when aligned energy levels
of the qubit-QD exist. If $U$ is estimated as the charging energy, then we
have $U\propto 46.4$ meV for a QD of 10 nm TanaAIP , which is significantly
larger than the operating region $V_{D}<$ 3meV considered here. We follow many
literature regarding the qubit-qubit operation in three QD system and focus on
the measurement mode here. Depending on the up and down-spin of the qubits,
singlet energy levels are represented by $E_{1}$ and $E_{1}+\Delta$,
respectively. ($\Delta$ denotes the Zeeman splitting energy). The energy
levels $E_{1}$ and $E_{1}+\Delta$ allow only down and up-spin of electrons,
respectively, as shown in the inset of Fig.1(b).

Measurements are performed when the lowest energy level of the center-QD is
increased to the upper energy levels for the second electron in the qubit QDs,
such that the current through the center-QD can detect $E_{1}$ or
$E_{1}+\Delta$, which depends on the spin state of the qubits. The qubit
states (spin state of the electron in the lowest energy level of the qubit-QD)
are indirectly affected by the measurements.

The measurement ($V_{D}\neq 0$) and interaction ($V_{D}=0$) modes are
separated by the potential energies of the center-QD, as shown in Fig. 1(c).
The energy levels of the qubits are controlled by the common gate, while that
of the center-QD are controlled by the source-to-drain voltage. The energy
levels of the qubits are assumed to change independently from that of the
center-QD. In the interaction mode (left side of Fig. 1(c)), the lowest energy
level of the center-QD is reduced such that an exchange interaction occurs
between the electrons. In the measurement mode (right side of Fig. 1(c)), the
lowest energy level of the center-QD is raised such that the qubits do not
interact directly with each other through electrons in the center-QD.

When a magnetic field is applied, the energy levels of the QDs and electrodes
are separated (Fig 1(d)) Sanchez . The measurement currents of up
($\uparrow$-current) and down ($\downarrow$-current) spins are assumed to be
treated independently by neglecting the spin flips between transports. Then,
the $\uparrow$-current and $\downarrow$-current can be treated to have
different Fermi energies $E_{F\pm}=E_{F}\pm\Delta/2$. A negative differential
conductance appears when the bottom of the Fermi energy of the source exceeds
the discrete energy level of the channel QD, as shown in Fig.1(d). As both the
QDs and electrodes are affected by the Zeeman splitting energy, the end bias
of the negative differential conductance coincides with $\uparrow$-current and
$\downarrow$-current

In the measurement mode (Fig. 2), four cases of energy levels exist in the
qubit states on both sides of the center-QD. The current of the channel-QDs
changes reflecting the different energy levels of the qubit-QDs. When two
qubit-QDs do not exhibit the same spin state (Fig. 2[i] and [ii]), the current
reflects one of the two qubit states with $E_{1}$ or $E_{1}+\Delta$. When the
two qubit-QDs have the same spin state (Fig. 2[iii] and [vi]), the current
reflects both the qubit states at the same energy level. For example, in Fig.
2[i], the left and right-qubits have the down and up-spin states,
respectively. The higher energy level of the left-qubit-QD is represented by
$E_{1}+\Delta$, while that of the right-qubit-QD is denoted as
$E_{3}(=E_{1})$. Here, the $\uparrow$-current resonates with the left-qubit,
but no resonance energy levels of the $\downarrow$-current between the center-
QD and the left-qubit exist. To identify the four cases of the currents
([i]$\sim$[vi]), we set the reference current as that of $W_{12}=W_{23}=0$,
and estimate the current deviations $\Delta I$ between the four cases of the
currents ([i]$\sim$[vi]) and the reference current. The reference current is
obtained when the energy level of the center-QD exists between the lowest and
highest energy levels of the qubit-QDs.

As aforementioned, the energy levels of the qubit-QDs are controlled by the
common gate and that of the center-QDs are controlled by the source and drain
voltage. Here, the common gate voltage is assumed to be fixed. Thus,
$E_{1}(=E_{3})$ is fixed and $E_{2}$ changes such that
$E_{2}=E_{2}^{(0)}+V_{D}/2$. Subsequently, cases [i]([iii]) and [ii]([iv])
occur when $E_{2}=E_{1}+\Delta$ ($E_{2}=E_{1}-\Delta$). In addition,
$E_{1}(=E_{3})>E_{2}^{(0)}$ is set such that $E_{2}$ crosses $E_{1}(=E_{3})$
as $V_{D}$ increases from $V_{D}=0$.

Figure 2:  Four current cases in the measurement mode reflect the four qubit
states. The black and blue circles show the electrons of the center-QD and
qubits-QDs, respectively. The white circle shows the vacant energy level in
which other electrons can enter. $E_{S}$ and $E_{T}$ show the singlet and
triplet energy-levels, respectively.

The Hamiltonian of the QD system depicted in Fig. 1(b) is given by

$\displaystyle H$ $\displaystyle=$
$\displaystyle\sum_{s}\sum_{i=1}^{3}E_{is}d_{is}^{\dagger}d_{is}\\!+\\!\\!\sum_{\alpha=L,R}\sum_{k_{\alpha},s}E_{k_{\alpha}s}c_{k_{\alpha}s}^{\dagger}c_{k_{\alpha}s}$
(1) $\displaystyle+$
$\displaystyle\\!\sum_{\alpha=L,R}\sum_{k_{\alpha},s}[V_{k_{\alpha},s,2}c_{k_{\alpha},s}^{\dagger}d_{2s}+V_{k_{\alpha}s}^{*}d_{2s}^{\dagger}c_{k_{\alpha}s}]$
$\displaystyle+$
$\displaystyle\sum_{i=1}^{2}\sum_{s}W_{i,i+1}(d_{is}^{\dagger}d_{i+1,s}+h.c.),$

where $c_{ks}^{\dagger}$ ($c_{ks}$) creates (annihilates) an electron of
momentum $k$ and spin $s(=\pm 1/2)$ in the electrodes, while
$d^{\dagger}_{is}$ ($d_{is}$) creates (annihilates) an electron in the QDs
($i=1,2,3$). $E_{k_{\alpha}s}=E_{k_{\alpha}}+sg\mu_{B}B$ is the energy level
of the electrode($\alpha=L,R$), and $E_{is}=E_{i0}+sg\mu_{B}B$ is the energy
level for three QDs ($i=1,2,3$).

The Zeeman splitting energy is expressed as $\Delta\equiv g\mu_{B}B$ in the
following calculations. One energy-level in each QD is assumed. As the Zeeman
splitting energy for the channel-QD and electrodes is considered in this
study, the Hamiltonian is slightly different from that of Ref. TanaPRB . It is
assumed that the current of the up-spin state is independent of that of the
down-spin state. Therefore, the formulation proposed in Ref. TanaPRB  can be
independently applied to $\uparrow$-current and $\downarrow$-current. The
coupling coefficients of the electrodes to the QDs are given by

$\Gamma_{\alpha,s}(\omega)=2\pi\sum_{k_{\alpha}}|V_{k_{\alpha},s}|^{2}\delta(\omega-
E_{k_{\alpha}s}),$ (2)

($\alpha=L,R$). Strictly speaking, $\Gamma_{\alpha,s}(\omega)$ depends on the
spin direction through the density of states. However, for the sake of
simplicity, we take
$\Gamma_{s}\equiv\Gamma_{\alpha,\uparrow}(\omega)=\Gamma_{\alpha,\downarrow}(\omega)(\alpha=L,R)$.
The current can be expressed by the transmission probability
$\mathcal{T}_{2s}(\omega)$, such that

$I_{D}=\sum_{s}\frac{e^{2}}{h}\int
d\omega\mathcal{T}_{2s}(\omega)[f_{S}(\omega)-f_{D}(\omega)],$ (3)

where $f_{\alpha}(\omega)=[\exp[(\omega-\mu_{\alpha})/(k_{B}T)]+1]^{-1}$
($k_{B}$, $\mu_{\alpha}$, and $T$) denote the Boltzmann constant, chemical
potential of the $\alpha$-electrode, and temperature, respectively (see
Appendix A)). The magnitude of $\Gamma_{\alpha,s}$ is determined such that the
calculated current is less than a few hundred pA Elzerman ; Shaji ; Lu . In
the present setup, the qubit states do not directly interact with the channel
electrons. Thus, it is considered that the effect of the channel current on
the coherence of the qubit is weak. A detailed calculation of the decoherence
time $t_{\rm dec}$ is formulated below (Eq.(7)). As $W_{ij}$ increases, the
current reflects the qubit state more strongly. In the following calculations,
$W_{ij}$ is selected such that the number of measurements exceeds a hundred.

The core model for $I_{D}$-$V_{D}$ of the fin field-effect transistor (FinFET)
provided by BSIM  is considered, given by

$I_{\rm D}=\beta(V_{\rm g}-V_{\rm th}-V_{\rm ds}/2)V_{\rm ds},$ (4)

where $\beta\equiv(\mu{\mathcal{W}}/L)(\epsilon/EOT)$ ($L=1~{}\mu$m,
${\mathcal{W}}=80$nm, $\mu=1000$cm2V-1s-1, $EOT$=1 nm, and $\epsilon=3.9\times
8.854\times 10^{-12}$ F/m represent the gate length, gate width, mobility,
oxide thickness, and dielectric constant, respectively BSIM ). The transistor
characteristics at low temperatures are under an intensive study by the CMOS
Society Pahwa . Therefore, the simple form of Eq.(4) is considered.

Numerical methods are applied to solve the equations of the QD device and
transistor, given by

$V_{D}=V_{\rm out}+V_{\rm ds}.$ (5)

For a given $V_{D}$, the current of the channel-QD, expressed by Eq.(3),
should be equal to that of the transistor, expressed by Eq.(4). Thus, this
equation is solved numerically using Newton’s method.

The purpose of the measurement is to distinguish between the different spin
states based on different magnitudes of the measured currents. Thus, the
current difference $\Delta I$ should be larger than the average noise level.
The possible number of measurements is estimated by comparing the decoherence
time $t_{\rm dec}$ with the measurement time. Following Ref. Schon , the
measurement time $t_{\rm meas}$ is defined using $\Delta I$ and shot noise$S$,
given by

$t_{\rm meas}^{-1}\equiv\frac{(\Delta I)^{2}}{4S},$ (6)

Because there are four cases of currents (Fig.2), four cases of the
measurement times are considered. The quantum shot noise $S$ is originally
expressed by numerous terms, and it has a simple form at $T=0$ Kobayashi ;
Lopez . The classical shot noise of $S=2e\langle I\rangle$ is found to be
larger than that in the quantum case. As shot noise increases, the measurement
time $t_{\rm meas}$ becomes longer and undesirable. Here, the classical form
is used as the worst case.

The decoherence time caused by the interactions is described using the golden
rule Schoelkopf . The last term of the Hamiltonian Eq. (1) is treated as the
perturbation term, given by $H_{\rm
int}(t)\equiv\sum_{i=1}^{2}\sum_{s}W_{i,i+1}(d_{i}^{\dagger}(t)d_{i+1}(t)+h.c.)$,
while the non-perturbated terms are those of the current lines ${\rm S}-{\rm
QD}_{2}-{\rm D}$. Then, $t_{\rm dec}$ can be defined as TanaPRB

$\displaystyle t_{\rm dec}^{-1}$ $\displaystyle\approx$
$\displaystyle\frac{1}{\hbar^{2}}\int_{-\infty}^{\infty}d\tau
e^{-i\omega_{01}\tau}\langle H_{\rm int}(\tau)H_{\rm int}(0)\rangle,$ (7)

(see Appendix A). Similar to $t_{\rm meas}$, four cases of $t_{\rm dec}$
corresponding to the four currents are considered.

Figure 3:  Current characteristics as a function of the applied voltage
$V_{D}$ in Fig.1(b) without the transistor. (a)Individual $\uparrow$-current
and $\downarrow$-current are shown separately. (b)Total currents are shown
where the black solid, dashed, and dotted lines denote the case of $W=0$
(reference current). $\Gamma_{0}=0.02u_{0}$ and $W=2u_{0}$ ($u_{0}=10^{-4}$
eV) is considered for $\Delta=0.16$ meV. The other solid lines represent the
four cases of measurement currents shown in Fig.2. $E_{1}=E_{3}=E_{F}+2u_{0}$,
$E_{2}^{(0)}=E_{F}+u_{0}$, $E_{F}$=1 meV, and $T=100$ mK.  Figure 4:  Current
characteristics as a function of the applied voltage $V_{D}$ in Fig.1(b) with
the transistor. (a) $W=2u_{0}$ and (b) $W=0.5u_{0}$ for $\Delta=0.16$ meV.
Other parameters are the same as those considered for Fig. 3. The black solid
lines show the case of $W=0$ (reference current), while the other solid lines
show the four cases of measurement currents shown in Fig.2.  Figure 5:  Zeeman
splitting energy $\Delta$ dependence of $t_{\rm dec}/t_{\rm meas}$ in
$V_{D}<0.5$meV region. $t_{\rm dec}/t_{\rm meas}$ is considered to correspond
to the approximate number of possible measurements before the measurement
backaction breaks the coherence. (a) $W=2u_{0}$ and (b) $W=0.5u_{0}$. Other
parameters are the same as those considered for Fig. 4.  Figure 6:  (a)(b)
Zeeman splitting energy $\Delta$ dependence of $t_{\rm dec}/t_{\rm meas}$ at
$V_{D}=2.225$ meV. The solid lines present the sum of currents of different
spin states. (c)(d) Decoherence time $t_{\rm dec}$ calculated using Eq.(7),
(a)(c) for the cases of [i] and [ii] and (b)(d) for the cases of [iii] and
[iv] shown in Fig.2. The parameters are the same as those of Fig. 4.

##  III Numerical Results

Figure 3 shows the $I_{D}$–$V_{D}$ characteristics of the system shown in
Fig.1(b) for different spin states without the transistor. The currents
increase from $V_{D}=0$ and decrease from approximately $V_{D}\sim 2.5$ meV.
Thus, a single resonant peak is observed in the present calculations. In Fig.
3(a), the $\uparrow$-current and $\downarrow$-current are presented
independently, where the dotted line (dashed line) represents the
$\uparrow$-current ($\downarrow$-current) for $W_{12}=W_{23}=0$ (reference
current). In Fig. 3(b), the four cases of channel currents, shown in Fig.2,
with the reference current have been illustrated. As an example, the case [i]
current consists of the case [i] $\uparrow$-current in Fig. 3(a), and the
$\downarrow$-spin reference current (dashed line in Fig. 3(a)). The range of
$V_{D}$ is chosen to cover the peak structure of the resonant tunneling. Cases
[i] and [iii], with the $\uparrow$-currents reflecting the qubit states, show
similar behaviors, while cases [ii] and [iv], with the $\downarrow$-currents
reflecting the qubit states, show similar behaviors.

Energy levels of $E_{1}=E_{3}=E_{F}+2u_{0}$ and $E_{2}^{(0)}=E_{F}+u_{0}$ are
selected such that these energy levels are within the range of the nonlinear
region around $V_{D}\approx 0$. The resonant peak ends when $E_{2}$ reaches
$V_{D}$ such that $E_{2}(=E_{2}^{(0)}+V_{D}/2)=V_{D}$ as shown in Fig.1(d),
implying that $V_{D}\approx 2E_{2}^{(0)}$. For $E_{F}=10u_{0}$ and
$E_{2}^{(0)}=E_{F}+u_{0}$, the peak structure ends at approximately
$V_{D}=22u_{0}=2.2$ meV. The end of the peaks of the $\uparrow$-current is the
same as that of the $\downarrow$-current (Fig.1(d)). On the other hand, the
current flow begins when $E_{2}$ is given by
$E_{2}(=E_{2}^{(0)}+V_{D}/2)=V_{D}+E_{F}$ which leads to
$V_{D}=2(E_{2}^{(0)}-E_{F})$. For the $\uparrow$-current
($\downarrow$-current) in which the effective Fermi surface is given by
$E_{F}+\Delta/2$ ($E_{F}-\Delta/2$), the bias for the beginning of the current
is given by $V_{D\uparrow}=2(E_{2}^{(0)}-E_{F}-\Delta/2)$
($V_{D\downarrow}=2(E_{2}^{(0)}-E_{F}+\Delta/2)$). Thus, the difference of the
beginning of the current is given by $V_{D\uparrow}-V_{D\downarrow}=2\Delta
E_{F}=2\Delta$, as seen in Fig. 3. The current difference $\Delta I$ in Eq.
(6) is calculated as the difference between currents [i]$\sim$[iv] and the
reference currents. As shown in Fig. 3, the difference $\Delta I$ is prominent
around $V_{D}\sim 0$ and $V_{D}>2.0$ meV.

Figure 4 illustrates the $I_{D}$-$V_{D}$ characteristics for the system shown
in Fig.1(b), with the transistor of gate length $L$=1 $\mu$m and $V_{\rm
g}-V_{\rm th}$=0.1V. By comparing Figs. 4 with Figs. 3, the $I_{D}$-$V_{D}$
characteristics shown in Fig. 4 are similar to that shown in Fig. 3. This is
because the resistance of the transistor of $L$=1 $\mu$m is considerably
smaller than that of the QD system. The effect of coupling $W$ for $W=0.05$
meV is smaller than that for the $W=0.2$ meV case. For a smaller $V_{D}$
region (left of the resonant peak), different effects are observed on the
$\uparrow$-current (Fig.2[i] and [iii]) and $\downarrow$-current(Fig.2[ii] and
[iv]). Conversely, for a larger $V_{D}$ region (right of the resonant peak),
the change from the reference current looks small for $W=0.05$ meV.

A large $\Delta I$ is desirable for a short measurement time $t_{\rm meas}$.
As can be seen in Fig. 4, the differences between the four currents
[i]$\sim$[iv] are prominent at $V_{D}<0.5$meV and 2.2$<V_{D}<$ 2.5 meV. We
consider that $t_{\rm dec}/t_{\rm meas}$ corresponds to the number of possible
measurements. Figure 5 shows $t_{\rm dec}/t_{\rm meas}$ in $V_{D}<0.5$meV
region. $t_{\rm dec}/t_{\rm meas}$ is found to exceed 100 for the cases [i]
and [iii]. However, those of cases [ii] and [iv] are less than 100. Thus, the
region $V_{D}<0.5$meV is preferable only for the cases [i] and [iii]. Figures
6 (a) and (b) show $t_{\rm dec}/t_{\rm meas}$ at $V_{D}=2.225$ meV. More than
100 measurements are expected for all cases [i]$\sim$[iv]. Figures 6 (c) and
(d) show the corresponding $t_{\rm dec}$ values as a function of $\Delta$.
$t_{\rm dec}$ increases for smaller $\Delta$, where the maximum value of
$t_{\rm dec}$ is taken as 1$\mu$s, reflecting the recent improvement in the
decoherence time of spin qubits Veldhorst ; Kobayashi2 . Thus, when the limit
of the decoherence time is increased, a greater number of measurements is
expected. Conversely, when the limit of the coherence time is decreased, a
smaller number of measurements is expected. Such a limit of the decoherence
time is considered to arise from the other origin of the qubit system, which
must be determined experimentally. The recessed patterns shown in Fig. 6(b)
are obtained from the abrupt structure of $t_{\rm dec}$, as provided in Ref.
TanaPRB .

$t_{\rm dec}/t_{\rm meas}$ varies depending on $E_{1}$ and $E_{2}$. Guidelines
for achieving the optimal $E_{2}^{(0)}$ and $E_{1}$ have not been established
other than that $E_{2}^{(0)}>E_{1}$ and $E_{2}^{(0)}$ is close to $E_{1}$ from
the results of the repeated simulations. $t_{\rm dec}/t_{\rm meas}>100$ for
$\Delta=0.25$ meV is feasible for $B_{z}=2.16$ T because 1 meV$\sim$8.64 T. If
we regard the possibility of more than hundred measurements as the less than
1% measurement error, the side-QD structure is considered to lead to the
future surface code.

Figure 7(a) shows $V_{\rm out}$ when the gate length of the transistor is
$L=1~{}\mu$m. The change of $V_{\rm out}$ is less than 50 $\mu$V because the
resistance of the channel-QD part ($\sim$ 1 meV /1nA$\sim 10^{6}\Omega$) is
much larger than that of the transistor of $L=1~{}\mu$m. Figure 7(b) shows the
case of $L=10~{}\mu$m, where the change in $V_{\rm out}$ exceeds 500 $\mu$V.
This result shows that long channel transistor is required to obtain
sufficiently large $V_{\rm out}$. Depending on $V_{\rm g}-V_{\rm th}$,
$I_{D}$-$V_{D}$ characteristics changes. Generally, the variations in the size
of QDs are inevitable, and appropriate bias $V_{D}$ changes for different
channel-QDs. The direct connection of the transistor to the qubit-system has
the advantage that variations of $I_{D}$-$V_{D}$ characteristics can be
adjusted by controlling the connected transistor. A larger $t_{\rm dec}/t_{\rm
meas}$ is observed for a larger $V_{\rm g}-V_{\rm th}$, but with a smaller
change in $V_{\rm out}$. Thus, a trade-off exists between the changes in
$V_{\rm out}$ and $V_{\rm g}-V_{\rm th}$. Detailed relationships are a topic
for future research. The $I_{D}$–$V_{D}$ characteristics and $t_{\rm
dec}/t_{meas}$ for $L=10~{}\mu$m are provided in Appendix B.

Figure 7:  Output $V_{\rm out}$ of the system shown in Fig.1(b), when the gate
length of the transistor $L$ changes. $V_{\rm g}-V_{\rm th}=0.1V$. (a) $L$=1
$\mu$m, and (b) $L$=10 $\mu$m. Other parameters are the same as those
considered for Fig. 6.

##  IV Discussion

In the present measurement, spin directions are not determined by a single
shot. Thus, the ratio $t_{\rm meas}/t_{\rm dec}$ does not exactly represent
the required number of measurements. The qubit state can be determined by
taking an average of repeated measurements. In any case, a large number of
possible measurements will lead to the future surface code. In previous
experiments ISSCC1 ; ISSCC2 , several amplification circuits have been used to
detect small differences between the up and down-spin states. Direct
connection of the transistor to the qubit is observed to increase the
resolution of the current differences (see Appendix C).

The effect of unexpected trap sites in the transistor should be considered
TanaAIP . The local electric field produced by static traps could be adjusted
by changing the transistor voltages. The effect of dynamic traps such as
random telegraph noise (RTN) would be significant because the voltage
fluctuation is in the order of mV. Owing to the long interval of RTN
($\sim\mu$s), repeated measurements might be desirable to identify the effect
of RTN.

Although Eq. (7) shows the first-order effect of tunneling between the qubit-
QD and center-QD, the decoherence time is a complicated function of the energy
levels of QDs as discussed in Ref. TanaPRB . The ratio $W/\Gamma_{0}$ is
larger than that in Ref. TanaPRB , indicating that the electrons have to spend
some time in the QD region to obtain information about the qubits. When
higher-order tunneling, such as cotunneling effects, is included, the formula
presented in Eq. (7) becomes significantly more complicated, and is a topic of
future research.

Note that, in the cases of [i] and [ii], which of QD1 and QD3 has larger
energy is not determined. What is known by the current detection for the cases
of [i] and [ii] is that the two qubits have different spin states. The
absolute energy of each qubit-QD is determined by detecting the edge qubits
because the corresponding current reflects only one qubit. The edge qubit is
described such that either $W_{12}$ or $W_{23}$ is zero. In this study, we
focused on the case for $W_{12}\neq 0$ and $W_{23}\neq 0$, because the edge
qubits behave similarly.

In this study, we did not mention quantum operations. Regarding the quantum
operations between the qubits, several previous proposals Fei ; Medford ;
Burkard ; Kandel ; Noiri  would be applied. In these quantum operations, a
magnetic-field gradient is required, and it is better that the micromagnets
are placed on the structure of Fig.1(a).

As an extension of Fig.1(a), the stacked structure can be further considered
based on a nanosheet which is one the near-future three-dimensional (3D)
semiconductor systems Chung ; Tsutsui ; Kim ; Ryckaert ; Tanaka ; Lee ; Lue
(see Appendix D). However, the fabrication of nanosheets itself requires
advanced fabrication technologies. This is also a future issue.

In the present setup, the qubit state is indirectly affected by the channel
current. To protect the qubit state more efficiently, addition of ancilla
qubit between the channel-QD and qubit-QD can be considered (Appendix E). A
detailed analysis will be conducted in the future.

##  V Conclusions

To summarize, the measurement process for several QD arrays is theoretically
investigated using resonant tunneling. Moreover, direct connections of the
qubit system to the conventional transistor are considered for reading
measurements. By using the nonlinear region of resonant tunneling phenomenon,
the small energy difference between the qubits is enhanced, resulting in the
possibility of more than hundreds of measurements during the decoherence time.
By connecting a conventional transistor to the qubits, the measurement results
are directly transmitted to the conventional CMOS circuits without long wires.
Although the proposed structure is expected to be an extension of state-of-
the-art semiconductor technology, a concrete fabrication process will be
investigated in the future. Since the proposed system has several parameters,
the determination of optimal parameters will depend on future experiments.

### DATA AVAILABILITY

The data that supports the findings of this study are available within the
article.

###### Acknowledgements.

We are grateful to S. Takagi, T. Mori and H. Fuketa for fruitful discussions.
This work was partly supported by MEXT Quantum Leap Flagship Program (MEXT
Q-LEAP), Grant Number JPMXS0118069228, Japan. This work was also supported by
JSPS KAKENHI Grant Number JP22K03497.

### AUTHOR CONTRIBUTIONS

T.T. conceived and designed the theoretical calculations. K.O. discussed the
results from viewpoint of experimentalist.

### CONFLICT OF INTEREST

The authors have no conflicts to disclose.

##  Appendix A Green’s Function Method

###  A.1 $I_{D}$-$V_{D}$ Characteristics

Green’s functions are derived using the equation-of-motion method Jauho .
Here, the formulation has been briefly explained. The time-dependent behavior
of the operator $d_{i}(t)$ is derived from
$i\hbar\frac{d~{}d_{i}(t)}{dt}=[H,d_{i}(t)]$, and we have

$\displaystyle\omega d_{i}(\omega)=[H,d_{i}(\omega)].$ (8)

Generally, the current through the QD system is described by the transport
properties of the left and right electrodes. Hereafter, the notations $L(=S)$
and $R(=D)$ have been used. The current $I_{Ls}$ for spins $s$ between the
central and left electrodes is given by

$\displaystyle I_{Ls}(t)$ $\displaystyle=$ $\displaystyle\frac{e}{\hbar}\int
dE{\rm
Re}\left\\{\sum_{k_{L}}V_{k_{L}s,i}G^{<}_{d_{is},c_{ik_{L}s}}(E)\right\\},$
(9)

where $G^{<}_{d_{is},c_{ik_{\alpha}s}}(t,t^{\prime})\equiv i\langle
c_{ik_{\alpha}s}^{\dagger}(t^{\prime})d_{is}(t)\rangle$, and
$G^{<}_{c_{ik_{\alpha}s},d_{is}}(t,t)=-[G^{<}_{d_{is},c_{ik_{\alpha}s}}(t,t)]^{*}$.
The total current is assumed to be conserved between the source (left
electrode) and drain (right electrode). Thus, the left electrode currents
$I_{L}$ and the right electrode currents $I_{R}$ satisfy the relation
$I_{L}=-I_{R}$. Then, the current through the device can be expressed as
$I=(I_{L}+I_{L})/2=(I_{L}-I_{R})/2$. The spin-flip process is neglected.
Although $\uparrow$-current and $\downarrow$-current have to be considered, in
this section, the suffix $s$ is omitted for the sake of convenience. The
current formula is expressed by the element free Green’s functions of QDs and
electrodes. The Green’s function of free QDs is given by

$\displaystyle g_{di}^{<}(\omega)$ $\displaystyle=$ $\displaystyle 2\pi
if_{i}(\omega)\delta(\omega-E_{i}),$ (10) $\displaystyle g_{di}^{>}(\omega)$
$\displaystyle=$ $\displaystyle 2\pi i[f_{i}(\omega)-1]\delta(\omega-E_{i}),$
(11) $\displaystyle g_{di}^{r}$ $\displaystyle=$
$\displaystyle\frac{1}{\omega-E_{i}+i\delta},$ (12) $\displaystyle
g_{d2}^{r}(\omega)$ $\displaystyle=$ $\displaystyle\frac{1}{\omega-
E_{2}+i\Gamma_{s}/2},$ (13) $\displaystyle\Gamma_{s}$ $\displaystyle=$
$\displaystyle\Gamma_{Ls}+\Gamma_{Rs},$ (14)

($i=1,3$).

The current $I_{Ls}$ provided in Ref. TanaPRB  is expressed as

$\displaystyle I_{Ls}$ $\displaystyle=$ $\displaystyle\frac{e}{h}\int
d\omega\Bigl{\\{}\frac{\Gamma_{Rs}\Gamma_{Ls}}{|\Delta_{c}^{r}|^{2}D_{2}(\omega)}[f_{2L}(\omega)-f_{2R}(\omega)]\Bigr{\\}},$
(15)

where

$\displaystyle D_{2s}(\omega)$ $\displaystyle=$ $\displaystyle(\omega-
E_{2})^{2}+(\Gamma_{Ls}+\Gamma_{Rs})^{2}/4,$ (16) $\displaystyle C_{ij}$
$\displaystyle=$ $\displaystyle|W_{ij}|^{2}g_{di}g_{dj},$ (17)
$\displaystyle\Delta_{c}$ $\displaystyle=$ $\displaystyle 1-C_{12}-C_{32}.$
(18)

Thus, the transmission probability is given by

$\displaystyle\mathcal{T}_{2s}(\omega)$ $\displaystyle=$
$\displaystyle\frac{\Gamma_{Rs}\Gamma_{Ls}}{|\Delta_{c}^{r}|^{2}D_{2s}(\omega)}.$
(19)

###  A.2 Derivation of the Decoherence Time

The decoherence time $t_{\rm dec}$ of the three QDs (Fig.1(b)) are the same as
that provided in Ref. TanaPRB , and is expressed as

$\displaystyle t_{\rm
dec}^{-1}\approx\frac{|W_{21}|^{2}}{\hbar^{2}}\Big{\\{}f(E_{1})\left(\frac{F_{2}^{>}(E_{1}-\omega_{01})}{[E_{1}-\omega_{01}-E_{2}]^{2}+\Gamma_{s}^{2}/4}\right)+[1-f(E_{1})]\left(\frac{F_{2}^{<}(E_{1}+\omega_{01})}{[E_{1}+\omega_{01}-E_{2}]^{2}+\Gamma_{s}^{2}/4}\right)\Big{\\}}$
(20) $\displaystyle+$
$\displaystyle\frac{|W_{32}|^{2}}{\hbar^{2}}\Big{\\{}f(E_{3})\left(\frac{F_{2}^{>}(E_{2}-\omega_{01})}{[E_{2}-\omega_{01}-E_{2}]^{2}+\Gamma_{s}^{2}/4}\right)+[1-f(E_{3})]\left(\frac{F_{2}^{<}(E_{2}+\omega_{01})}{[E_{2}+\omega_{01}-E_{2}]^{2}+\Gamma_{s}^{2}/4}\right)\Big{\\}},$

where $\omega_{01}=|E_{1}-E_{2}|$ and,

$\displaystyle F_{2}^{<}(\omega)$ $\displaystyle=$
$\displaystyle\Gamma_{Ls}f_{L}(\omega)+\Gamma_{Rs}f_{R}(\omega),$ (21)
$\displaystyle F_{2}^{>}(\omega)$ $\displaystyle=$
$\displaystyle\Gamma_{Ls}[1-f_{L}(\omega)]+\Gamma_{Rs}[1-f_{R}(\omega)].$ (22)

Here, $f(E_{i})$ and $1-f(E_{i})$ indicate the existence of an electron at the
$E_{i}$ level and no electron at the $E_{i}$ level ($i=1,3$), respectively.
For the average case, $f(E_{i})=1/2$ is considered in the following
calculations.

##  Appendix B Temperature Dependence and the Case of $L=10~{}\mu$m

Figure 8:  (a) $I_{D}$-$V_{D}$ characteristics and (b) $t_{\rm dec}/t_{\rm
meas}$ at $T=200$ mK, $W=0.2$ meV, and $V_{\rm g}-V_{\rm th}=0.1$ V. Other
parameters are the same as those considered for Fig. 6.

Figure 8 illustrates the number of possible measurements at $T=200$ mK, which
is compared to the results shown in Fig. 6 in the main text. As the
temperature increases, the resonant peak blurs and the sensitivity weakens.
Thus, the operating temperature is important while performing the
measurements.

Figure 9 shows the $I_{D}$–$V_{D}$ characteristics and $t_{\rm meas}/t_{\rm
dec}$, where the gate length of the transistor is given by $L=10\mu$m with
$W=0.2$ meV and $V_{\rm g}-V_{\rm th}=0.1$ eV. On comparing Fig. 4 and Fig. 6,
the $I_{D}$–$V_{D}$ characteristics are observed to have changed owing to the
large resistance of the transistor.

Figure 9:  (a) $I_{D}$-$V_{D}$ characteristics and (b) $t_{\rm dec}/t_{\rm
meas}$ of $L=10~{}\mu$m at $T=100$ mK. $\Gamma_{0}=2\mu$eV, and $V_{\rm
g}-V_{\rm th}=0.1$ eV. Other parameters are the same as those considered for
Fig. 6.

##  Appendix C Two-dimensional Qubit Array with Common Gate Structure

Figure 10:  (a) A two-dimensional qubit array as an extension of Fig.1 in the
main text. (b) Qubit-qubit operations are conducted in both x and
y-directions. The channel-QDs work as the controlling line to mediated
interactions between the qubits and for the measurement of two adjacent
qubits. (c) Different phase pulses used to avoid undesirable operations.

A two-dimensional qubit array consisting of the units shown in Fig.1 of the
main text is considered. Figure 10(a) shows the schematic of a qubit array. A
single qubit operation is assumed to have been performed using the magnetic-
field gradient produced by the micromagnets around the array Noiri1 ; Noiri .
The channel-QDs (green circles) are shared by two qubits (yellow and white
circles), such that the changes in the electric potential of the channel-QDs
affect both sides of the qubits (crosstalk), as shown in Fig. 10(b). In other
word, each spin direction is changed by the two surrounding nearest channel-
QDs. For example, by changing the pulse voltage on the nearest channel-QDs,
the qubit sandwiched between the two channel-QDs experiences an oscillating
field, which results in a spin direction change under the magnetic-field
gradient. Because each channel-QD is shared with two qubits, when a pulse
voltage is applied to the channel-QD, the voltage change affects the qubits on
both sides of the channel-QDs. Thus, the random-access process is not suitable
in this structure. However, the cross-talk can be avoided, as shown in Fig.
10(c). By applying different phase pulses to the nearest channel-QDs, the
effect of the applied pulse between the center qubits can be negated.

When the channel-QDs are used as the control line to change the spin
direction, as shown in Fig. 10(a) and (b), the resistance between the QD and
the source and drain is greater than 50 $\Omega$ owing to the tunneling
resistance between the channel-QD and electrodes. Thus, a high-frequency
signal does not pass through the channel-QD and is reflected. It produces an
additional electric field to the qubits, as if the channel-QD component
functions as the gate electrode. Detailed discussion on the operations are
beyond the scope of this study.

Figure 11:  Cross-section of the stacked qubits. (a) Qubit system embedded
into nanosheet as a modification of the system discussed in Ref. TanaAIP . (b)
As an extension of (a), toy system of two-dimensional qubit array. The system
shown in Figure 10 is stacked on the substrate using a common gate, which
controls the potential of $E_{2}$, shown in Fig.1 in the main text. The qubit-
QDs and the channel-QDs are placed side by side. The qubits A and B have the
same structure.  Figure 12:  Bird’s eye view of the 3D stacked qubit system.
(a) Qubits are surrounded by channel-QDs. (b) Sources and drains of different
channel-QDs are independently accessed by attaching different bias lines to
the sources and drains.

##  Appendix D Common Gate Stacked Qubit Structure

Here, we would like to propose a stacked qubit structure as a typical
embodiment of Fig. 10. Particularly, the common gate concept is applied to the
nanosheet to reduce the number of wires required in the system. The nanosheet
is one of the most advanced semiconductor structures available, where the
channel parts are stacked on the substrate. Qubits are represented by QDs
embedded between the channels, as described in Ref. TanaAIP . Figure 11 shows
a cross-section of the proposed stacked qubits as an extension of the
nanosheet. The gate electrode shared by the channel-QDs controls the electric
potential of the QD system. Figure 12(a) shows a bird’s eye view of the
stacked qubit system.

Figure 11(a) shows a simple extension of the system discussed in Ref. TanaAIP
, which is rotated by 90 degree of the previous proposal using a FinFET
structure. The qubits are embedded into the oxide layer between the two
channels. Here, the magnetic gradient method is applied and the two nearest
metallic wires are used to generate an oscillating electric field gradient on
the electrons in the qubit. In this structure, the distance between qubits in
the lateral direction is still larger than that in the vertical direction;
therefore, the lateral coupling is weaker than the vertical coupling. Figure
11(b) shows the advanced structure in which the coupling magnitude of the
lateral qubits is equal to that of the vertical coupling. The channel-QDs play
both the roles of detection and interaction between the nearest qubits, as
discussed in Fig. 10.

Figure 12 (a) shows the bird’s eye view of the 3D stacked qubit system. The
source and drain of each channel-QDs are separated to independently control
the Fermi energies of the QDs. The interaction of the qubits is switched on
and off according to the position of the Fermi energy. The size of the contact
to each source and drain should be sufficiently small. Such a structure is
different from that of a typical nanosheet but its feasibility will be high
when fabrication technologies for 3D NAND flash memories are applied. Although
the case in which the channel consists of a QD has been discussed, nanowires
can replace the center-QD.

Figure 13:  Extended version of the qubit system in which additional ancilla
qubits (QD2 and QD4) are placed between the channel-QD(QD3) and target qubits
(QD1 and QD5) to reduce the effect of backaction. If the target qubit size is
smaller than the ancilla qubits, the tunneling of electrons between the target
qubits and the electrodes is prohibited. In addition, because the target
qubits are set close to the electrodes, their energy levels can be controlled
by the potential of the electrode relative to the common gate. Due to the
different sizes of the QDs, the energy levels are different.

##  Appendix E Five QD Operation

Yoneda et al.  Yoneda  proposed that ancilla QD enhances the coherence of
target qubits. Hence, we can also consider the ancilla QD between the channel-
QD and target qubits. Figure 13 shows the schematic of the 5-QD case, where
the ancilla QDs (QD2 and QD4) are placed between the center-QD (QD3) and the
target qubits (QD1 and QD5). The target qubits, which are placed close to the
electrodes, are smaller than the ancilla QDs. This leads to stronger
capacitive coupling between the target qubits and electrodes where no
tunneling between the qubits and electrodes is assumed. As QD1 and QD5 are
close to the electrodes, their potentials can change significantly by the
applied bias, as compared with those of the ancilla QD2 and QD4. The electrons
are introduced into QD1 and QD5 by lowering the potentials of the electrodes
to a value lower than that of the channel-QD controlled by the common gate.
When a large bias is applied, to avoid current flow through the channel-QD,
the S/D voltages are set to be equal. Detailed discussion on the operations
are beyond the scope of this study.

## References

  * (1) T. Hensgens, T. Fujita, L. Janssen, Xiao Li, C. J. Van Diepen, C. Reichl, W. Wegscheider, S. Das Sarma, and L. M. K. Vandersypen, Nature 548, 70 (2017).
  * (2) A. R. Mills, M. M. Feldman, C. Monical, P. J. Lewis, K. W. Larson, A. M. Mounce, and J. R. Petta, Appl. Phys. Lett. 115, 113501 (2019).
  * (3) F. Fedele, A. Chatterjee, S. Fallahi, G. C. Gardner, M. J. Manfra, and F. Kuemmeth, PRX Quantum 2, 040306 (2021).
  * (4) M. Veldhorst, H. G. J. Eenink, C. H. Yang, and A. S. Dzurak, Nat. Commun. 8, 1766 (2017).
  * (5) W. Ha, S. D. Ha, M. D. Choi, Y. Tang, A. E. Schmitz, M. P. Levendorf, K. Lee, J. M. Chappell, T. S. Adams, D. R. Hulbert, E. Acuna, R. S. Noah, J. W. Matten, M. P. Jura, J. A. Wright, M. T. Rakher, and M. G. Borselli, Nano Lett. 22, 1443 (2022).
  * (6) S. G. J. Philips, M. T. Madzik, S. V. Amitonov, S. L. De Snoo, M. Russ, N. Kalhor, C. Volk, W. I. L. Lawrie, D. Brousse, L. Tryputen, B. P. Wuetz, A. Sammak, M. Veldhorst, G. Scappucci, and L. M. K. Vandersypen, Nature 609, 919 (2022).
  * (7) J. Yoneda, K. Takeda, A. Noiri, T. Nakajima, S. Li, J. Kamioka, T. Kodera, and S. Tarucha, Nat. Commun. 11, 1144 (2020).
  * (8) M. Veldhorst, J. C. C. Hwang, C. H. Yang, a. W. Leenstra, B. de Ronde, J. P. Dehollain, J. T. Muhonen, F. E. Hudson, K. M. Itoh, A. Morello, and A. S. Dzurak, Nat. Nanotechnol. 9, 981 (2014).
  * (9) A. Noiri, K. Takeda, T. Nakajima, T. Kobayashi, A. Sammak, G. Scappucci, and S. Tarucha, Nat. Communications 13, 5740 (2022).
  * (10) T. Utsugi, N. Lee, R. Tsuchiya, T. Mine, R. Mizokuchi, J. Yoneda, T. Kodera, S. Saito, D. Hisamoto, and H. Mizuno, Jpn. J. Appl. Phys. 62, SC1020 (2023).
  * (11) F. K. Malinowski, F. Martins, T. B. Smith, S. D. Bartlett, A. C. Doherty, P. D. Nissen, S. Fallahi, G. C. Gardner, M. J. Manfra, C. M. Marcus, and F. Kuemmeth, Nat Commun 10, 1196 (2019).
  * (12) X. G. Croot, S. J. Pauka, J. D. Watson, G. C. Gardner, S. Fallahi, M. J. Manfra, and D. J. Reilly Phys. Rev. Appl. 10, 044058 (2018).
  * (13) T. Tanamoto and K. Ono, AIP Advances 11, 045004 (2021).
  * (14) T. Tanamoto and T. Aono, Phys. Rev. B 106, 125401 (2022).
  * (15) A. P. Jauho, N. S. Wingreen, and Y. Meir, Phys. Rev. B 50, 5528 (1994).
  * (16) J. Fei, D. Zhou, Y. -P. Shim, S. Oh, X. Hu, and M. Friesen, Phys. Rev. A 86, 062328 (2012).
  * (17) M. Russ, J. R. Petta, and G. Burkard, Phys. Rev. Lett. 121, 177701 (2018).
  * (18) Y. P. Kandel, H. Qiao, and J. M. Nichol Appl. Phys. Lett. 119, 030501 (2021).
  * (19) A. Noiri, T. Nakajima, J. Yoneda, M. R. Delbecq, P. Stano, T. Otsuka, K. Takeda, S. Amaha, G. Allison, K. Kawasaki, Y. Kojima, A. Ludwig, A. D. Wieck, D. Loss, and S. Tarucha, Nat. Communications 9,5066 (2018).
  * (20) J. Medford, J. Beil, J. M. Taylor, E. I. Rashba, H. Lu, A. C. Gossard, and C. M. Marcus Phys. Rev. Lett. 111, 050501 (2013).
  * (21) F. Capasso, K. Mohammed, and A. Y. Cho, Appl. Phys. Lett. 48, 478 (1986).
  * (22) K. Fujita, T. Edamura, S. Furuta, and M. Yamanishi, Appl. Phys. Lett. 96, 241107 (2010).
  * (23) D. Kim, D. R. Ward, C. B. Simmons, D. E. Savage, M. G. Lagally, M. Friesen, S. N. Coppersmith, and M. A. Eriksson, npj Quantum Inf 1, 15004 (2015).
  * (24) ”Nanostructures and Mesoscopic Systems”, edited by W. P. Kirk and M. A. Reed (Academic Press, San Diego, 1992).
  * (25) Proceedings of the 21st International Conference on the Physics of Semiconductors, edited by P. Jiang and H-Z Zheng (World Scientific, Singapore, 1992).
  * (26) A. G. Fowler, A. C. Whiteside, A. L. McInnes, and A. Rabbani, Phys. Rev. X 2, 041003 (2012).
  * (27) Y. S. Chauhan, D. D. Lu, V. Sriramkumar, S. Khandelwal, J. P. Duarte, N. Payvadosi, A. Niknejad, and C. Hu, FinFET Modeling for IC Simulation and Design: Using the BSIM-CMG Standard. (Academic Press, London, 2015).
  * (28) B. Jena, S. Dash, and G. P. Mishra, IEEE Trans. Electron Devices, 65, 2422, June (2018).
  * (29) T. Tanamoto and K. Ono, Appl. Phys. Lett. 119, 174002 (2021).
  * (30) Y. Peng, A. Ruffino, T. -Y. Yang, J. Michniewicz, M. F. Gonzalez-Zalba and E. Charbon, IEEE J. Solid-State Circuits, 57, 2374 (2022).
  * (31) J. -S. Park, S. Subramanian, L. Lampert, T. Mladenov, I. Klotchkov, D. Kurian, E. Juárez-Hernández, B. P. Esparza, S. R. Kale, K. T. A. Beevi, S. P. Premaratne, T. Watson, S. Suzuki, M. Rahman, J. Timbadiya, S. Soni, and S. Pellerano 2021 IEEE International Solid- State Circuits Conference (ISSCC), 208 (2021).
  * (32) Y. -Y. Chung, B. -J. Chou, C. -F. Hsu, W. -S. Yun, M. -Y. Li, S. -K. Su, Y. -T. Liao, M. -C. Lee, G .-W. Huang, S. -L. Liew, Y. -Y. Shen, W. -H. Chang, C. -W. Chen, C. -C. Kei, H. Wang, H. -S. P. Wong, T. Y. Lee, C. -H. Chien, C. -C. Cheng, and I. P. Radu, 2022 International Electron Devices Meeting (IEDM), 34.5 (2022).
  * (33) G. Tsutsui, S. Song, J. Strane, R. Xie, L. Qin, C. Zhang, D. Schmidt, S. Fan, B. Hong, Y. Jung, C. -W. Sohn, I. Hwang, J. Yim, G. H. Son, G. Jo, K-I. Kim, M. Sankarapandian, S. Mochizuki, I. Seshadri, E. Miller, J. Li, J. Demarest, C. Waskiewicz, R. G. Southwick, H. Zhou, R. N. Pujari, P. Nieves, M. Wang, H. Jagannathan, B. Anderson, D. Guo, R. Divakaruni, T. Wu, K-I. Seo, and H. Bu, 2022 International Electron Devices Meeting (IEDM), 34.4 (2022).
  * (34) S. D. Kim, M. Guillorn, I. Lauer, P. Oldiges, T. Hook, and M.-H. Na, Microelectron. Technol. Unified Conf. (S3S), 1, (2015).
  * (35) J. Ryckaert, M. H. Na, P. Weckx, D. Jang, P. Schuddinck, B. Chehab, S. Patli, S. Sarkar, O. Zografos, R. Baert, and D. Verkest, 2019 IEEE International Electron Devices Meeting (IEDM), 29.4 (2019).
  * (36) T. Tanaka et al., 2016 IEEE International Solid-State Circuits Conference (ISSCC), 142 (2016).
  * (37) J. Lee, J. Jang, J. Lim, Y. G. Shin, K. Lee, and E. Jung, 2016 IEEE International Electron Devices Meeting (IEDM), 11.2 (2016).
  * (38) H. -T. Lue et al., 2020 IEEE International Electron Devices Meeting (IEDM), 6.4 (2020).
  * (39) D. Sánchez, C. Gould, G. Schmidt, and L. W. Molenkamp, IEEE Trans. Electron Devices 54, 984 (2007).
  * (40) J. M. Elzerman, R. Hanson, L. H. Willems van Beveren, B. Witkamp, L. M. K. Vandersypen, and L. P. Kouwenhoven, Nature 430, 431 (2004).
  * (41) N. Shaji, C. B. Simmons, Madhu Thalakulam, Levente J. Klein, Hua Qin, H. Luo, D. E. Savage, M. G. Lagally, A. J. Rimberg, R. Joynt, M. Friesen, R. H. Blick, S. N. Coppersmith, and M. A. Eriksson, Nature Phys 4, 540 (2008).
  * (42) W. Lu, Z. Ji, L. Pfeiffer, K. W. West, and A. J. Rimberg, Nature 423, 422 (2003).
  * (43) G. Pahwa, P. Kushwaha, A. Dasgupta, S. Salahuddin, and C. Hu, IEEE Trans. Electron Devices, 68, 4223 (2021).
  * (44) Y. Makhlin, G. Schön, and A. Shnirman, Rev. Mod. Phys. 73, 357 (2001).
  * (45) K. Kobayashi and M. Hashisaka, J. Phys. Soc. Jpn. 90, 102001 (2021).
  * (46) R. López, R. Aguado, and G. Platero, Phys. Rev. B 69, 235305 (2004). Physics Reports 336, 1 (2000).
  * (47) R. J. Schoelkopf, A. A. Clerk, S. M. Girvin, K. W. Lehnert, M. H. Devoret, Proceedings, Conference, Subtitle or Series: Quantum Noise in Mesoscopic Physics 175-203 (2003).
  * (48) T. Kobayashi, J. Salfi, C. Chua, J. van der Heijden, M. G. House, D. Culcer, W. D. Hutchison, B. C. Johnson, J. C. McCallum, H. Riemann, N. V. Abrosimov, P. Becker, H. -J. Pohl, M. Y. Simmons and S. Rogge, Nat. Mater. 20, 38 (2021).