# Supercurrent in Bi4Te3 Topological Material-Based Three-Terminal Junctions

Jonas Kölzer  Peter Grünberg Institut (PGI-9), Forschungszentrum Jülich,
Wilhelm-Johnen-Straße, 52425 Jülich, Germany  JARA-Fundamentals of Future
Information Technology, Jülich-Aachen Research Alliance, Forschungszentrum
Jülich and RWTH Aachen University, 52425 Jülich, Germany  Abdur Rehman Jalil
Peter Grünberg Institut (PGI-9), Forschungszentrum Jülich, Wilhelm-Johnen-
Straße, 52425 Jülich, Germany  JARA-Fundamentals of Future Information
Technology, Jülich-Aachen Research Alliance, Forschungszentrum Jülich and RWTH
Aachen University, 52425 Jülich, Germany  Daniel Rosenbach  Peter Grünberg
Institut (PGI-9), Forschungszentrum Jülich, Wilhelm-Johnen-Straße, 52425
Jülich, Germany  JARA-Fundamentals of Future Information Technology, Jülich-
Aachen Research Alliance, Forschungszentrum Jülich and RWTH Aachen University,
52425 Jülich, Germany  Lisa Arndt  JARA Institute for Quantum Information,
RWTH Aachen University, Germany  Gregor Mussler  Peter Grünberg Institut
(PGI-9), Forschungszentrum Jülich, Wilhelm-Johnen-Straße, 52425 Jülich,
Germany  JARA-Fundamentals of Future Information Technology, Jülich-Aachen
Research Alliance, Forschungszentrum Jülich and RWTH Aachen University, 52425
Jülich, Germany  Peter Schüffelgen  Peter Grünberg Institut (PGI-9),
Forschungszentrum Jülich, Wilhelm-Johnen-Straße, 52425 Jülich, Germany  JARA-
Fundamentals of Future Information Technology, Jülich-Aachen Research
Alliance, Forschungszentrum Jülich and RWTH Aachen University, 52425 Jülich,
Germany  Detlev Grützmacher  Peter Grünberg Institut (PGI-9),
Forschungszentrum Jülich, Wilhelm-Johnen-Straße, 52425 Jülich, Germany  JARA-
Fundamentals of Future Information Technology, Jülich-Aachen Research
Alliance, Forschungszentrum Jülich and RWTH Aachen University, 52425 Jülich,
Germany  Hans Lüth  Peter Grünberg Institut (PGI-9), Forschungszentrum Jülich,
Wilhelm-Johnen-Straße, 52425 Jülich, Germany  JARA-Fundamentals of Future
Information Technology, Jülich-Aachen Research Alliance, Forschungszentrum
Jülich and RWTH Aachen University, 52425 Jülich, Germany  Thomas Schäpers
<EMAIL_ADDRESS>Peter Grünberg Institut (PGI-9), Forschungszentrum
Jülich, Wilhelm-Johnen-Straße, 52425 Jülich, Germany  JARA-Fundamentals of
Future Information Technology, Jülich-Aachen Research Alliance,
Forschungszentrum Jülich and RWTH Aachen University, 52425 Jülich, Germany

###### Abstract

In an in-situ prepared three-terminal Josephson junction based on the
topological insulator Bi4Te3 and the superconductor Nb the transport
properties are studied. The differential resistance maps as a function of two
bias currents reveal extended areas of Josephson supercurrent including
coupling effects between adjacent superconducting electrodes. The observed
dynamics for the coupling of the junctions is interpreted using a numerical
simulation of a similar geometry based on a resistively and capacitively
shunted Josephson junction model. The temperature dependency indicates that
the device behaves similar to prior experiments with single Josephson
junctions comprising topological insulators weak links. Irradiating radio
frequencies to the junction we find a spectrum of integer Shapiro steps and an
additional fractional step, which is interpreted by a skewed current-phase
relationship. In a perpendicular magnetic field we observe Fraunhofer-like
interference patterns of the switching currents.

##  I Introduction

Hybrid structures comprising three-dimensional topological insulator
nanoribbons combined with superconductors are a very promising platform for
realizing circuits for fault-tolerant topological quantum computing [1, 2, 3,
4]. For its operation Majorana bound states are employed, which are formed by
aligning an external magnetic field with a nanoribbon proximitized with an
s-type superconductor [5, 6, 7]. For the braiding of different pairs of
Majorana states for qubit operation multi-terminal structures are required [8,
2, 9]. Braiding can be performed by adjusting the superconducting phase of the
superconducting electrodes to each other.

Multi-terminal Josephson junctions are the backbone of Majorana braiding
mechanism in a topological qubit; where a three-terminal Josephson junction
acts as a basic building block [2]. Understanding the superconducting
transport in such a device holds a key importance for the realization of a
topological quantum system. Generally, the use of hybrid devices with multiple
connections leads to rich physics in terms of transport properties. Indeed,
theoretical studies have investigated singularities, such as Weyl nodes, in
the Andreev spectra of multi-terminal Josephson junctions [10, 11, 12].
Moreover, multi-terminal Josephson junctions with topologically trivial
superconducting leads may lead to realizations where the junction itself can
be regarded as an artificial topological material [13]. Furthermore, three-
terminal junctions also allow transport via the quartet mechanism and non-
local Andreev processes [14, 15, 16, 17].

On the experimental side, multi-terminal Josephson junctions were fabricated
with different materials for the weak link. In three-terminal Josephson
junctions with a Cu or InAs nanowire subgap states [18, 19] and half-integer
Shapiro steps [20] were observed, indicating transport via quartets of
entangled Cooper pairs. Supercurrent flow affected by dissipative currents in
an adjacent junction was studied on graphene-based junctions [21]. Moreover,
the higher-dimensional phase space was found to lead to fractional Shapiro
steps in this type of junctions due to the inverse AC Josephson effect [4]. By
combining a multi-terminal junction with a top gate, the effect of gate
voltage and magnetic field on the critical current contour has been studied
[3, 24, 25]. Recently, flakes of the topological insulator Bi2Se3 were also
used as a weak link in an interferometer structure, and evidence for a non-
sinusoidal current-phase relationship was observed [26]. In flux-controlled
three-terminal junctions based on Bi2Te3, the opening and closing of a minigap
was studied using normal probes [27].

Here, we report on the transport properties of a three-terminal Josephson
junction based on the Bi4Te3 material system as the weak link and Nb as the
superconductor. To fabricate the samples, we used selective-area growth for
the Bi4Te3 layer in combination with an in-situ bridge technology to define
the superconducting electrodes [2]. Bi4Te3 is a natural superlattice of
alternating Bi2 bilayers and Bi2Te3 quintuple layers. Initially, Bi4Te3 has
been reported to be a semimetal with zero band gap and a Dirac cone at the
$\Gamma$ point [29]. However, recent band structure calculations in
conjunction with scanning tunneling spectroscopy and angular photoemission
spectroscopy measurements suggest that the material is a semimetal with
topological surface states [30, 31, 32]. In particular, advanced $GW$-band
structure calculations have shown that a band gap of about
$0.2\text{\,}\mathrm{eV}$ opens at the $\Gamma$ point, which significantly
reduces the density of the bulk state in this energy range [32]. Bi4Te3 is
classified as a dual topological insulator, a strong topological insulator
with a non-zero mirror Chern number, i.e. a topological crystalline insulator
phase. Though Bi4Te3 does not exhibit the proposed Dirac semimetal phase, it
is still a very interesting material as it resides in close proximity to the
critical point of band crossing in the topological phase diagram of BixTey
alloys [33]. Such a transition is proposed by Yang et al. [34] where a
topological crystalline insulator (Bi2Te3) [35] can be topologically
transformed into a topological Dirac semimetal through alloying it with other
materials. On our multi-terminal junctions, we first investigated the DC
properties and related the results to simulations based on the resistively and
capacitively shunted Josephson junction (RCSJ) model. We then measured the
radio frequency (rf) response, finding evidence for coupling of adjacent
junctions. Finally, the behavior of our three-terminal junctions when an out-
of-plane magnetic field is applied is investigated.

##  II Experimental

Using the previously introduced technologies of topological insulator
selective-area growth and in-situ bridge technology we fabricated three-
terminal Josephson junctions, as illustrated in Fig. 1(a) [36, 2]. The
geometry of the nanoribbon T-shaped junction for selective-area growth is
defined by trenches in a SiO2/Si3N4 (5 nm/15 nm) layer on a highly-resistive
Si (111) substrate [37]. First, the 600-nm-wide nanotrenches are etched into
the top Si3N4 layer using a combination of electron beam lithography and
reactive ion etching. Subsequently, a second set of layers, i.e. a 100-nm-
thick SiO2 layer and a 300-nm-thick Si3N4 layer, is deposited on top to define
the stencil mask for the in-situ Nb deposition [2]. After patterning the
structures for the stencil mask into Si3N4, SiO2 is etched in hydrofluoric
acid (HF) forming the free-hanging bridge structures. Simultaneously, the
Si(111) surface in the selective-area growth trenches is released in the
bottom SiO2 layer defined by the Si3N4 layer on top. The Bi4Te3 layer is
selectively grown within these trenches, while the Si3N4 bridge structures are
employed to define the geometry of the in situ deposited superconducting
electrodes [2]. The Bi4Te3 layer is grown at a temperature of $310^{\circ}$C
using molecular beam epitaxy. Subsequently, the 50-nm-thick superconducting Nb
electrodes are deposited by electron beam evaporation followed by covering the
whole structure with a 5-nm-thick Al2O3 dielectric capping layer. Our
processing scheme ensured a high-quality crystalline topological insulator
material with clean superconductor interfaces [38, 2], as reported in previous
transmission electron microscopy studies. An electron microscopy image of the
investigated device is presented in Fig. 1(b).

The measurements of the three-terminal Josephson junction were carried out in
a dilution refrigerator with base temperature of $T=25\,\textnormal{mK}$.
containing a $1\text{ - }1\text{ - }6\,\textnormal{T}$ vector magnet. As
indicated in Fig. 1(b), the left, right, and bottom junction electrodes are
labeled as ”L”, ”R”, and ”B”, respectively. Two current sources supply
currents $I_{\mathrm{LB}}$ and $I_{\mathrm{RB}}$ from L and R to the bottom
electrode, respectively, with the according voltages $V_{\mathrm{LB}}$ and
$V_{\mathrm{RB}}$ measured. The differential resistances are measured by
adding an ac current of 10 nA to the DC current bias using a lock-in
amplifier. The rf-irradiation for the Shapiro step measurements was provided
by an antenna placed in close vicinity to the sample.

Figure 1: Rendering of a selective-area grown three-terminal Josephson
junction and false color scanning electron micrograph with circuit: $(a)$ The
three-terminal junction is composed of the silicon substrate (gray bottom
layer), the first hard mask composed out of a silicon oxide (white)/ silicon
nitride (blue) layer (as indicated by the labels). On top of this another hard
mask layer composed of silicon oxide (white) and silicon nitride (blue) is
deposited and patterned as a shadow mask. The topological insulator (red) is
grown selectively into the first hard mask trench and the shadow mask is used
for the definition of the junction in the metal deposition (silver) step.
$(b)$ False-color scanning electron micrograph of the in-situ prepared three-
terminal junction device. Niobium contacts (cyan) are deposited on top of the
TI (red). The measurement configuration is also indicated.

##  III Results and Discussion

### DC characteristics

Information about the basic junction characteristics is obtained by measuring
the differential resistances $R_{\mathrm{LB}}=\Delta V_{\mathrm{LB}}/\Delta
I_{\mathrm{LB}}$ and $R_{\mathrm{RB}}=\Delta V_{\mathrm{RB}}/\Delta
I_{\mathrm{RB}}$ as a function of the bias currents $I_{\mathrm{LB}}$ and
$I_{\mathrm{RB}}$, respectively. Starting with the left junction we find that
$R_{\mathrm{LB}}$ shown in Figs 2(a) and (b) contains a superconducting region
in the center when $I_{\mathrm{LB}}$ and $I_{\mathrm{RB}}$ are varied. The
observed critical current contour is similar to what has been observed in
induced superconducting nano junctions made of high mobility materials such as
InAs/Al [3, 24] or graphene [21].

Figure 2: Differential resistance maps: (a) shows $R_{\mathrm{LB}}$ as a
function of the bias currents $I_{\mathrm{LB}}$ and $I_{\mathrm{RB}}$ at 25 mK
with corresponding line cuts given in (b). In (c) the differential resistance
map of $R_{\mathrm{RB}}$ is depicted with a selection of line cuts given in
(d). The dashed lines in a and (c) indicate the superconducting regions of
compensating bias currents. The differential resistances was measured by using
lock-in technique, i.e. $R_{\mathrm{LB}}=\Delta V_{\mathrm{LB}}/\Delta
I_{\mathrm{LB}}$ and $R_{\mathrm{RB}}=\Delta V_{\mathrm{RB}}/\Delta
I_{\mathrm{RB}}$.

The superconducting region extends along an inclined line indicated by the
dashed line in Fig. 2(a). The switching to the superconducting state can be
seen in the line cuts at fix values $I_{\mathrm{RB}}=0$ and $\pm 0.7\mu$A
provided in Fig. 2(b). The extension of the superconducting state originates
from a part of $I_{\mathrm{RB}}$ which flows via R to L through the junction
between L and B compensating the current $I_{\mathrm{LR}}$ partly and by that
reducing the total current. For our three-terminal device no reduced
differential resistance is observed along the line
$I_{\mathrm{LB}}=I_{\mathrm{RB}}$, which would indicate the presence of a
Josephson supercurrent between the junction formed between electrodes L and R
[3, 25]. We attribute this to the fact that the distance between these
electrodes is slightly larger than for the other junctions so that no
Josephson supercurrent is obtained. However, the junction between L and R acts
as a shunt resistor taking care that the switching to the superconducting
state is non-hysteretic. The differential resistance $R_{\mathrm{RB}}$
measured between R and B electrodes, depicted in Figs 2(c) and (d), shows a
similar behaviour as $R_{\mathrm{LB}}$, i.e. featuring also an extended
superconducting range due a compensation provided by part of
$I_{\mathrm{LR}}$. The tilt of the superconducting range indicated by the
dashed line in Figure 2(c) is lower compared to Fig. 2(a) since now
$I_{\mathrm{LR}}$ is the compensating current.

### Simulations

The experimental results are modeled by assuming a network of two resistively
and capacitively shunted Josephson (RCSJ) junctions coupled by a resistor
$R_{C}$, as illustrated in Fig. 3(a). Solving the related system of
differential equations numerically, in analogy to what was presented in
previous works [3, 4], we simulate the behaviour of the experimental system
(information about the procedure see Supplementary Material). The results of
the simulations are shown in Figs. 3(b) to (e), where the differential
resistance $R_{\mathrm{LB}}$ is given as a function of the bias currents
$I_{\mathrm{LB}}$ and $I_{\mathrm{RB}}$.

Figure 3: Numerical simulation of different coupling scenarios: (a) The three-
terminal circuit is modeled by two RCSJ shunted Josephson junctions
$J_{\mathrm{LB}}$ and $J_{\mathrm{RB}}$ (green and blue dashed line boxes),
which are each modeled by a resistor $R_{N}$, a capacitor c, and an ideal
Josephson junction $J$. Currents $I_{\mathrm{LB}}$ and $I_{\mathrm{RB}}$ are
supplied via current sources while the voltage drops $V_{\mathrm{LB}}$ and
$V_{\mathrm{RB}}$ across the junctions are measured. Both junctions are
coupled via a coupling resistance $R_{C}$. (b) Differential resistance
$R_{\mathrm{LB}}$ as a function of current biases for a realistic scenario for
$R_{C}$ close to the one extracted in the experiment:
$R_{N}=$40\text{\,}\mathrm{\SIUnitSymbolOhm}$$,
$R_{C}=$160\text{\,}\mathrm{\SIUnitSymbolOhm}$$,
$I_{c}=$538\text{\,}\mathrm{n}\mathrm{A}$$,
$\beta_{c}=(2e/\hbar)I_{c}R_{N}^{2}C=0.1$. The zero resistance range is
observed as a tilted line due to a compensation by a part of
$I_{\mathrm{RB}}$. Additionally, the influence of the second junction is
observed as a similar line close to horizontal orientation. The corresponding
line cuts indicated in (b) are presented in (c). The scenario for a very small
coupling resistance ($R_{C}\rightarrow 0$) is shown as a color map of
$R_{\mathrm{LB}}$ and selected line cuts in (d) and (e).

The model describes the experiment well by reproducing the Josephson
supercurrent along the inclined lines originating from compensating currents
from both electrodes with a superconducting region at the center. The
inclination is determined by the coupling resistance $R_{C}$. In Figs. 3(b)
and (c), the coupling resistance was taken as $R_{C}=4\cdot R_{\mathrm{LB}}$,
with $R_{\mathrm{LB}}=$40\text{\,}\mathrm{\SIUnitSymbolOhm}$$ which results in
the same tilt as observed experimentally. Taking these values into account the
normal state resistance is given by $R_{N}=6/5\cdot
R_{\mathrm{LB}}=$48\text{\,}\mathrm{\SIUnitSymbolOhm}$$. In our simulations
for the critical current and for the Steward-McCumber parameter we assumed
$I_{c}=$538\text{\,}\mathrm{n}\mathrm{A}$$ and
$\beta_{c}=(2e/\hbar)I_{c}R_{N}^{2}C=0.1$, respectively, with c the junction
capacitance. We found that the superconducting state in the junction between R
and B leads to some weak feature as a similar line inclined towards horizontal
orientation. Note, that for this line $R_{\mathrm{LB}}$ is non-zero, as the
supercurrent in the other junction only partly reduces the current in the
junction between L and B and hence only partially reduces the voltage drop. A
noticeable difference between experiment and simulation is that in the
measurements the extension of the superconducting state observed along the
inclined line (cf. Fig. 2(a) is decreased compared to the simulation depicted
in Fig. 3(b). As discussed by Draelos et al. [21], this effect can be
explained by dissipation in the neighboring junction being in the normal state
resulting in an effective heating, in particular for junctions with small
dimensions. In our simulation the direct coupling between the different
junctions was neglected. As shown by Arnault et al. [4], including coupling
results in a more complex contour of the critical current area. If the
coupling resistance becomes very small, i.e. $R_{C}\rightarrow 0$, the
observed lines in the differential resistance shift towards the diagonal (cf.
Figs. 3(d) and (e). Thus, both junctions are maximally correlated to both
current biases $I_{\mathrm{LB}}$ and $I_{\mathrm{RB}}$.

### Temperature dependence

In Figs. 4(a) to (f) the differential resistance maps are shown for
$R_{\mathrm{LB}}$ and $R_{\mathrm{RB}}$ measured at temperatures of 100 mK,
200 mK, and 800 mK.

Figure 4: Differential resistance maps at various temperatures: Left column
(a), (c), (e) shows the differential resistance $R_{\mathrm{LB}}$, right
column (b), (d), ($\mathbf{f}$) $R_{\mathrm{RB}}$, accordingly. The
temperatures displayed in the rows from up to down are
$100\text{\,}\mathrm{mK}$, $200\text{\,}\mathrm{mK}$, and
$800\text{\,}\mathrm{mK}$, respectively.

One finds that with increasing temperature the area of the central
superconducting region shrinks. This is in accordance with the temperature
dependence of the critical current of a single Nb/Bi4Te3/Nb reference
junction, as shown in the Supplementary Material. It is noteworthy that the
superconducting feature along the inclined lines basically does not change
with increasing temperature. This can be explained by the fact, that the
dissipation in the neighboring junction already leads to an increased
temperature being larger than the substrate temperature.

### rf characteristics

Next the radio frequency response of the system is investigated in order to
confirm that the experiment is described well by Josephson junction physics
and to analyze the rf response of the Josephson current. This is done by first
choosing a frequency and an amplitude for the rf irradiation so that both
junctions show a large rf response in the differential resistance.
Subsequently the same DC bias sweeps are performed as in the prior
experiments. Figures 5(a) and (b) show Shapiro step measurements of the
differential resistances $R_{\mathrm{LB}}$ and $R_{\mathrm{RB}}$,
respectively, as a function of bias currents $I_{\mathrm{LB}}$ and
$I_{\mathrm{RB}}$. The differential resistances are calculated by numerical
differentiation. Differential resistances obtained by lock-in amplifier
measurements can be found in the Supplementary Material.

Figure 5: Shapiro step measurements at $5.8\text{\,}\mathrm{GHz}$: (a)
Numerically determined differential resistance $R_{\mathrm{LB}}$ as a function
of $I_{\mathrm{LB}}$ and $I_{\mathrm{RB}}$ at $5.8\text{\,}\mathrm{GHz}$ and
rf power of 0 dBm. (b) Corresponding map of the differential resistance
$R_{\mathrm{RB}}$.

The rf frequency $f_{\mathrm{rf}}$ and the according power was set to
$5.8\text{\,}\mathrm{GHz}$ and 0 dBm, respectively. The differential
resistances show clear intercrossing stripe-like patterns which can be
attributed to the presence of Shapiro steps confirming the presence of a
Josephson supercurrent. The intercrossing parallel stripes indicating a
coupling of both junctions. By calculating the related voltage drop we find
that for both junctions the Shapiro steps are located at integer multiples,
$n=1,\,2,\,3\dots$, of the characteristic voltage $V_{0}=hf_{\mathrm{rf}}/2e$.

In Figs. 6(a) and (b) the differential resistance maps of $R_{\mathrm{LB}}$
and $R_{\mathrm{LB}}$, now taken at $8.5\text{\,}\mathrm{GHz}$ at 0 dBm, are
depicted, respectively.

Figure 6: Shapiro step measurements at $8.5\text{\,}\mathrm{GHz}$: (a)
Numerically determined differential resistance $R_{\mathrm{LB}}$ as a function
of the normalized voltage drops $V_{\mathrm{LB}}/V_{0}$ and
$V_{\mathrm{RB}}/V_{0}$ at $8.5\text{\,}\mathrm{GHz}$ and rf power of 0 dBm,
with $V_{0}=hf_{\mathrm{rf}}/2e$. The blue curves represent the averaged
signal along $V_{\mathrm{LB}}/V_{0}$ and $V_{\mathrm{RB}}/V_{0}$,
respectively. The dashed lines indicate the half-integer steps. (b)
Corresponding map of the differential resistance $R_{\mathrm{RB}}$.

Here, the color maps are plotted as a function of the normalized voltages
$V_{\mathrm{LB}}/V_{0}$ and $V_{\mathrm{RB}}/V_{0}$. On first sight one finds
that the Shapiro step pattern is more pronounced in $R_{\mathrm{LB}}$. We
attribute this to a stronger coupling of the rf signal compared to the
neighbouring junction due to spatial variations of the rf field. As for the
measurements at $5.8\text{\,}\mathrm{GHz}$ a coupling of both junctions,
although weaker, is observed. Our experimental results of Shapiro step
measurements are supported by comparison to simulations based on the
previously introduced RCSJ model. In Supplementary Figures 4(a) and (b) maps
of the simulated values of $R_{\mathrm{LB}}$ and $R_{\mathrm{LB}}$ as a
function of the normalized bias voltages are shown. There, one finds that the
coupling by $R_{C}$ results in a weak cross coupling of the Shapiro signal
resulting in intercrossing stripe-like patterns of different contrast.

A closer inspection of the resistance map presented in Fig. 6(a) reveals that
apart from the integer Shapiro steps also half-integer Shapiro steps, e.g. at
$n=1/2$, are observed. The half integer steps are also clearly resolved in the
averaged value of $R_{\mathrm{LB}}$ along $V_{\mathrm{LB}}/V_{0}$ shown in
Fig. 6(a). In single Josephson junctions such fractional steps are interpreted
by assuming a skewed current-phase relationship [39, 40, 41] (a simulation for
this case using our model is provided in the Supplementary Material). More
specifically for multi-terminal junctions the rf response of superconductivity
induced into normal metal has been studied previously by Duvauchelle et al.
[20]. Here, half-integer steps have been found and interpreted as a feature
due to the presence of coherent quartet states. However, in Fig. 2 we did not
find indications of quartet states, which would be visible by a feature in the
differential resistance at opposite voltage drops on the left and right
terminal [18]. Other experimental observations of such fractional steps in
multi-terminal junctions are interpreted on the basis of highly connected
nonlinear networks of Josephson junctions, where (due to the higher phase
space) different transitions of the phase particle in the washboard potential
are possible [4]. However, since fractional Shapiro steps were observed in
single junctions made with similar materials [42], we favor the explanation
based on a skewed current-phase relationship, which can be attributed to
contributions of quasi-ballistic transport. In our measurements under rf
radiation we did not find indications of missing odd Shapiro steps, as
predicted when Majorana bound states are present in topological junctions [43,
2]. Probably, for our samples the narrow width of the Bi4Te3 ribbons prevents
the formation of these states, since due to the finite Berry phase a magnetic
field along the junctions is required to gain a gap closure for the coherent
surface states around the nanoribbon cross section [36].

### Magnetic field response

The junction characteristics were also analyzed in a perpendicularly oriented
magnetic field $B_{\perp}$. In Fig. 7(a) the magnetic field dependence
$R_{\mathrm{LB}}$ is plotted as a function of $B_{\perp}$ and
$I_{\mathrm{LB}}$, while $I_{\mathrm{RB}}$ is kept at zero. One clearly
observes a Fraunhofer-like interference pattern of the switching current, i.e.
the boundary between the red superconducting areas and the areas with finite
resistance. The blue line in Fig. 7(a) indicates the according fitting based
on the Fraunhofer interference relation. The close resemblance of the
experimental data to an ideal Fraunhofer pattern points towards a relatively
homogeneous distribution of the supercurrent density. From the fit we extract
a period of about $\Delta B=$$14\text{\,}\mathrm{mT}$, which corresponds to a
junction area of $152\times 10^{3}\,\mathrm{nm}^{2}$. Relating these values to
the dimensions of the left junction $J_{\mathrm{LB}}$ one finds that the
period is about a factor of ten smaller than expected. Based on the actual
junction size of $200\times 72\,\mathrm{nm}^{2}$ a period of
$144\text{\,}\mathrm{mT}$ is expected for a $h/2e$ flux periodicity. We
attribute the discrepancy to the experimental period to a pronounced flux
focusing effect, where the magnetic field is expelled from the edge regions of
the superconducting electrodes and bundled in the junction area. As a matter
of fact, a comparably large flux focusing effect was previously observed in
similar planar Josephson junctions based on topological insulators and Nb
superconducting electrodes [36].

Figure 7: Differential resistances under perpendicular magnetic field sweep:
(a) shows a map of $R_{\mathrm{LB}}$ as a function of $B_{\perp}$ and
$I_{\mathrm{LB}}$ for $I_{\mathrm{RB}}=0$. (b) represents the corresponding
map of $R_{\mathrm{RB}}$ as a function of $B_{\perp}$ and $I_{\mathrm{RB}}$
for $I_{\mathrm{LB}}=0$. In (c) $R_{\mathrm{LB}}$ is plotted with the sweep
current chosen to be $I_{\mathrm{LB}}=I_{\mathrm{RB}}$, which corresponds a
sweep along the diagonal in the current plane. In all cases a standard
Fraunhofer pattern is fitted indicated as blue lines.

In Fig. 7(b) the magnetic field dependence $R_{\mathrm{RB}}$ is shown as a
function of $B_{\perp}$ and $I_{\mathrm{RB}}$ at $I_{\mathrm{LB}}=0$. Once
again, a Fraunhofer-like interference is observed, although with a smaller
period, i.e. an larger effective area where the magnetic flux is picked up.
The reason for the difference compare to the measurements shown in Fig. 7(a)
might be some inhomogeneity in the supercurrent density in the junction.
Finally, the $R_{\mathrm{LB}}$ maps are scanned diagonally, i.e.
$I_{\mathrm{LB}}=I_{\mathrm{RB}}$, as shown in Fig. 7(c). Here, once again a
regular Fraunhofer pattern is observed, which is almost identical to the
pattern shown in Fig. 7(d), indication, that the current $I_{\mathrm{RB}}$
through the neighboring junction basically has not effect.

##  IV Conclusion

We have succeeded in extending the previously developed in situ fabrication
technology for Josephson junctions to a working more complex design of a
three-terminal junction. Analysis of the transport experiments shows that our
system indeed behaves like a coupled network of Josephson junctions in DC
transport, rf response, as well as magnetic field response. This is the first
report on the topological multi-terminal devices where an interaction between
the individual Josephson junctions is observed. Moreover, the observation of
fractional steps in the rf response opens a window that provides a first
insight into the novel physics of this type of device. On a more technical
level, our results demonstrate the realization of more complex devices
required for network structures in topological quantum circuits.

Further investigations and detailed understanding of such a system are crucial
for the realization of complex topological quantum systems. In future, similar
experiments with more intricate circuit designs and superconducting phase
controlled measurements will be performed. The complexities in the junction
characteristics arose from the selected weak-link material Bi4Te3. In future
experiments we plan to incorporate conventional three-dimensional topological
insulators, e.g. Bi2Te3, Sb2Te3, Bi2Se3, and the topological Dirac semimetal
exhibited by the correctly tuned BixTey stoichiometric alloy.

## Acknowledgement

We thank Herbert Kertz for technical assistance as well as Kristof Moors and
Roman Riwar for fruitful discussion. This work was partly funded by the
Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under
Germany’s Excellence Strategy—Cluster of Excellence Matter and Light for
Quantum Computing (ML4Q) EXC 2004/1—390534769, the German Federal Ministry of
Education and Research (BMBF) via the Quantum Futur project “MajoranaChips”
(Grant No. 13N15264) within the funding program Photonic Research Germany, as
well as the Bavarian Ministry of Economic Affairs, Regional Development and
Energy within Bavaria’s High-Tech Agenda Project “Bausteine für das
Quantencomputing auf Basis topologischer Materialien mit experimentellen und
theoretischen Ansätzen“ (grant allocation no. 07 02/686 58/1/21 1/22 2/23).

## References

  * Kitaev [2003] A. Kitaev, Fault-tolerant quantum computation by anyons, Annals of Physics 303, 2 (2003).
  * Hyart  _et al._ [2013] T. Hyart, B. van Heck, I. C. Fulga, M. Burrello, A. R. Akhmerov, and C. W. J. Beenakker, Flux-controlled quantum computation with Majorana fermions, Phys. Rev. B 88, 035121 (2013).
  * Aasen  _et al._ [2016] D. Aasen, M. Hell, R. V. Mishmash, A. Higginbotham, J. Danon, M. Leijnse, T. S. Jespersen, J. A. Folk, C. M. Marcus, K. Flensberg, and J. Alicea, Milestones toward Majorana-based quantum computing, Phys. Rev. X 6, 031016 (2016).
  * Manousakis  _et al._ [2017] J. Manousakis, A. Altland, D. Bagrets, R. Egger, and Y. Ando, Majorana qubits in a topological insulator nanoribbon architecture, Phys. Rev. B 95, 165424 (2017).
  * Cook and Franz [2011] A. Cook and M. Franz, Majorana fermions in a topological-insulator nanowire proximity-coupled to an $s$-wave superconductor, Phys. Rev. B 84, 201105 (2011).
  * Cook  _et al._ [2012] A. M. Cook, M. M. Vazifeh, and M. Franz, Stability of Majorana fermions in proximity-coupled topological insulator nanowires, Phys. Rev. B 86, 155431 (2012).
  * Legg  _et al._ [2021] H. F. Legg, D. Loss, and J. Klinovaja, Majorana bound states in topological insulators without a vortex, Phys. Rev. B 104, 165405 (2021).
  * Fu and Kane [2008] L. Fu and C. L. Kane, Superconducting proximity effect and Majorana fermions at the surface of a topological insulator, Phys. Rev. Lett. 100, 096407 (2008).
  * Stenger  _et al._ [2019] J. P. T. Stenger, M. Hatridge, S. M. Frolov, and D. Pekker, Braiding quantum circuit based on the $4\pi$ Josephson effect, Phys. Rev. B 99, 035307 (2019).
  * Yokoyama and Nazarov [2015] T. Yokoyama and Y. V. Nazarov, Singularities in the Andreev spectrum of a multiterminal Josephson junction, Phys. Rev. B 92, 155437 (2015).
  * Eriksson  _et al._ [2017] E. Eriksson, R.-P. Riwar, M. Houzet, J. S. Meyer, and Y. V. Nazarov, Topological transconductance quantization in a four-terminal Josephson junction, Phys. Rev. B 95, 075417 (2017).
  * Xie  _et al._ [2018] H.-Y. Xie, M. G. Vavilov, and A. Levchenko, Weyl nodes in Andreev spectra of multiterminal Josephson junctions: Chern numbers, conductances, and supercurrents, Phys. Rev. B 97, 035443 (2018).
  * Riwar  _et al._ [2016] R.-P. Riwar, M. Houzet, J. S. Meyer, and Y. V. Nazarov, Multi-terminal Josephson junctions as topological matter, Nature Comm. 7, 1 (2016).
  * Houzet and Samuelsson [2010] M. Houzet and P. Samuelsson, Multiple Andreev reflections in hybrid multiterminal junctions, Phys. Rev. B 82, 060517 (2010).
  * Freyn  _et al._ [2011] A. Freyn, B. Douçot, D. Feinberg, and R. Mélin, Production of nonlocal quartets and phase-sensitive entanglement in a superconducting beam splitter, Phys. Rev. Lett. 106, 257005 (2011).
  * Nowak  _et al._ [2019] M. P. Nowak, M. Wimmer, and A. R. Akhmerov, Supercurrent carried by nonequilibrium quasiparticles in a multiterminal Josephson junction, Phys. Rev. B 99, 075416 (2019).
  * Melo  _et al._ [2022] A. Melo, V. Fatemi, and A. R. Akhmerov, Multiplet supercurrent in Josephson tunneling circuits, SciPost Phys. 12, 17 (2022).
  * Pfeffer  _et al._ [2014] A. H. Pfeffer, J. E. Duvauchelle, H. Courtois, R. Mélin, D. Feinberg, and F. Lefloch, Subgap structure in the conductance of a three-terminal Josephson junction, Phys. Rev. B 90, 075401 (2014).
  * Cohen  _et al._ [2018] Y. Cohen, Y. Ronen, J.-H. Kang, M. Heiblum, D. Feinberg, R. Mélin, and H. Shtrikman, Nonlocal supercurrent of quartets in a three-terminal Josephson junction, Proceedings of the National Academy of Sciences 115, 6991 (2018).
  * Duvauchelle  _et al._ [2016] J. E. Duvauchelle, A. H. Pfeffer, H. Courtois, and F. Lefloch, Quantum coherence of the quartet scheme observed by Shapiro resonance under radio-frequency irradiation in three terminal Josephson junctions, IEEE Transactions on Applied Superconductivity 26, 1 (2016).
  * Draelos  _et al._ [2019] A. W. Draelos, M.-T. Wei, A. Seredinski, H. Li, Y. Mehta, K. Watanabe, T. Taniguchi, I. V. Borzenets, F. Amet, and G. Finkelstein, Supercurrent flow in multiterminal graphene Josephson junctions, Nano Letters 19, 1039–1043 (2019).
  * Arnault  _et al._ [2021] E. G. Arnault, T. F. Q. Larson, A. Seredinski, L. Zhao, S. Idris, A. McConnell, K. Watanabe, T. Taniguchi, I. Borzenets, F. Amet, and G. Finkelstein, Multiterminal inverse ac Josephson effect, Nano Letters 21, 9668 (2021).
  * Graziano  _et al._ [2020] G. V. Graziano, J. S. Lee, M. Pendharkar, C. J. Palmstrøm, and V. S. Pribiag, Transport studies in a gate-tunable three-terminal Josephson junction, Phys. Rev. B 101, 054510 (2020).
  * Pankratova  _et al._ [2020] N. Pankratova, H. Lee, R. Kuzmin, K. Wickramasinghe, W. Mayer, J. Yuan, M. G. Vavilov, J. Shabani, and V. E. Manucharyan, Multiterminal Josephson effect, Phys. Rev. X 10, 031051 (2020).
  * Graziano  _et al._ [2022] G. V. Graziano, M. Gupta, M. Pendharkar, J. T. Dong, C. P. Dempsey, C. Palmstrøm, and V. S. Pribiag, Selective control of conductance modes in multi-terminal Josephson junctions, arXiv preprint arXiv:2201.01373 (2022).
  * Kurter  _et al._ [2015] C. Kurter, A. Finck, Y. S. Hor, and D. J. Van Harlingen, Evidence for an anomalous current–phase relation in topological insulator Josephson junctions, Nature Comm. 6, 1 (2015).
  * Yang  _et al._ [2019] G. Yang, Z. Lyu, J. Wang, J. Ying, X. Zhang, J. Shen, G. Liu, J. Fan, Z. Ji, X. Jing, F. Qu, and L. Lu, Protected gap closing in Josephson trijunctions constructed on Bi2Te3 , Phys. Rev. B 100, 180501 (2019).
  * Schüffelgen  _et al._ [2019] P. Schüffelgen, D. Rosenbach, C. Li, T. W. Schmitt, M. Schleenvoigt, A. R. Jalil, S. Schmitt, J. Kölzer, M. Wang, B. Bennemann, U. Parlak, L. Kibkalo, S. Trellenkamp, T. Grap, D. Meertens, M. Luysberg, G. Mussler, E. Berenschot, N. Tas, A. A. Golubov, A. Brinkman, T. Schäpers, and D. Grützmacher, Selective area growth and stencil lithography for in situ fabricated quantum devices, Nature Nanotechnology 14, 825 (2019).
  * Saito  _et al._ [2017] Y. Saito, P. Fons, K. Makino, K. V. Mitrofanov, F. Uesugi, M. Takeguchi, A. V. Kolobov, and J. Tominaga, Compositional tuning in sputter-grown highly-oriented Bi–Te films and their optical and electronic structures, Nanoscale 9, 15115 (2017).
  * Chagas  _et al._ [2020] T. Chagas, G. A. S. Ribeiro, P. H. R. Gonçalves, L. Calil, W. S. Silva, Â. Malachias, M. S. C. Mazzoni, and R. Magalhães-Paniago, Bi2:Bi2Te3 stacking influence on the surface electronic response of the topological insulator Bi4Te3, Electronic Structure 2, 015002 (2020).
  * Chagas  _et al._ [2022] T. Chagas, O. A. Ashour, G. A. S. Ribeiro, W. S. Silva, Z. Li, S. G. Louie, R. Magalhães Paniago, and Y. Petroff, Multiple strong topological gaps and hexagonal warping in Bi4Te3, Phys. Rev. B 105, L081409 (2022).
  * Nabok  _et al._ [2022] D. Nabok, M. Tas, S. Kusaka, E. Durgun, C. Friedrich, G. Bihlmayer, S. Blügel, T. Hirahara, and I. Aguilera, Bulk and surface electronic structure of Bi4Te3 from GW calculations and photoemission experiments, Phys. Rev. Materials 6, 034204 (2022).
  * Jalil [2022] A. R. Jalil,  _Engineering topological superlattices and their epitaxial integration in selectively grown hybrid nanostructures via MBE_ , Ph.D. thesis, RWTH Aachen University (2022).
  * Yang and Nagaosa [2014] B.-J. Yang and N. Nagaosa, Classification of stable three-dimensional dirac semimetals with nontrivial topology, Nature Comm. 5, 1 (2014).
  * Rauch  _et al._ [2014] T. c. v. Rauch, M. Flieger, J. Henk, I. Mertig, and A. Ernst, Dual topological character of chalcogenides: Theory for Bi2Te3, Phys. Rev. Lett. 112, 016802 (2014).
  * Rosenbach  _et al._ [2021] D. Rosenbach, T. W. Schmitt, P. Schüffelgen, M. P. Stehno, C. Li, M. Schleenvoigt, A. R. Jalil, G. Mussler, E. Neumann, S. Trellenkamp, A. A. Golubov, A. Brinkman, D. Grützmacher, and T. Schäpers, Reappearance of first Shapiro step in narrow topological Josephson junctions, Science Advances 7, eabf1854 (2021).
  * Schmitt  _et al._ [2022] T. W. Schmitt, M. R. Connolly, M. Schleenvoigt, C. Liu, O. Kennedy, J. M. Chávez-Garcia, A. R. Jalil, B. Bennemann, S. Trellenkamp, F. Lentz, _et al._ , Integration of topological insulator Josephson junctions in superconducting qubit circuits, Nano Letters 22, 2595 (2022).
  * Kölzer  _et al._ [2021] J. Kölzer, K. Moors, A. R. Jalil, E. Zimmermann, D. Rosenbach, L. Kibkalo, P. Schüffelgen, G. Mussler, D. Grützmacher, T. L. Schmidt, H. Lüth, and T. Schäpers, In-plane magnetic field-driven symmetry breaking in topological insulator-based three-terminal junctions, Comm. Materials 2, 116 (2021).
  * Snyder  _et al._ [2018] R. A. Snyder, C. J. Trimble, C. C. Rong, P. A. Folkes, P. J. Taylor, and J. R. Williams, Weak-link Josephson junctions made from topological crystalline insulators, Phys. Rev. Lett. 121, 097701 (2018).
  * Panghotra  _et al._ [2020] R. Panghotra, B. Raes, C. C. de Souza Silva, I. Cools, W. Keijers, J. E. Scheerder, V. V. Moshchalkov, and J. Van de Vondel, Giant fractional Shapiro steps in anisotropic Josephson junction arrays, Comm. Physics 3, 53 (2020).
  * Raes  _et al._ [2020] B. Raes, N. Tubsrinuan, R. Sreedhar, D. S. Guala, R. Panghotra, H. Dausy, C. C. de Souza Silva, and J. Van de Vondel, Fractional Shapiro steps in resistively shunted Josephson junctions as a fingerprint of a skewed current-phase relationship, Phys. Rev. B 102, 054507 (2020).
  * Rosenbach [2021] D. Rosenbach,  _Quantum transport and induced superconductivity in selectively deposited topological insulator devices_ , Dissertation, RWTH Aachen University, Aachen (2021).
  * Domínguez  _et al._ [2017] F. Domínguez, O. Kashuba, E. Bocquillon, J. Wiedenmann, R. S. Deacon, T. M. Klapwijk, G. Platero, L. W. Molenkamp, B. Trauzettel, and E. M. Hankiewicz, Josephson junction dynamics in the presence of $2\pi$\- and $4\pi$-periodic supercurrents, Phys. Rev. B 95, 195430 (2017).

Supplementary Material: Supercurrent in Bi4Te3 Topological Material-Based
Three-Terminal Junctions

##  SI Single junction measurements

As a reference a single Nb/Bi4Te3/Nb junction was measured. The junction has a
length of $140\text{\,}\mathrm{nm}$ and a width of $500\text{\,}\mathrm{nm}$.
In Supplementary Figure S1(a) the current voltage characteristics is shown at
temperatures in the range from $30\text{\,}\mathrm{m}\mathrm{K}$ to
$0.77\text{\,}\mathrm{K}$. At lowest temperature a critical current of
$750\text{\,}\mathrm{nA}$ is obtained. In contrast to the three terminal
junction, here, a hysteretic behaviour is observed, which can be explained by
the missing shunt for the single Josephson junction. We attribute the
hysteresis to heating resulting in a lower return current $I_{r}$ compared to
$I_{c}$ [1]. The critical current monotonously decreases with temperature with
some kink around $0.4\text{\,}\mathrm{K}$. The latter might be attributed to a
switching from a more diffusive to a more ballistic transport in the weak link
[2].

Supplementary Figure S1: Current-voltage characteristics of a single
Nb/Bi4Te3/Nb junction: (a) Current-voltage characteristics at temperatures
ranging from $30\text{\,}\mathrm{m}\mathrm{K}$ to $0.77\text{\,}\mathrm{K}$.
(b) Critical current $I_{c}$ as well as return current $I_{r}$ as a function
of temperature.

##  SII RCSJ model for a three-terminal junction

The characteristics of our three-terminal junctions is simulated by employing
a two-dimensional resistively and capacitively shunted Josephson junction
(RCSJ) Ansatz in analogy to what was presented in previous works [3, 4]. In
Fig. 3(a) in the main text the corresponding network is depicted including two
resistively and capacitively shunted Josephson junctions with the normal state
resistance $R_{N}$ and the capacitance $C$. We assume two identical junctions
each having a critical current of $I_{c}$. The junctions are connected by a
coupling resistor $R_{C}$ representing the non-superconducting junction
between electrodes L and R. Following the RCSJ Ansatz the characteristics of
the three-terminal junction can be described by a set of coupled differential
equations of the form:

$\frac{I_{\mathrm{LB}}}{I_{c}}=\sin(\varphi_{\mathrm{LB}})+\frac{d\varphi_{\mathrm{LB}}}{d\widetilde{\tau}}+\beta_{c}\frac{d^{2}\varphi_{\mathrm{LB}}}{d\widetilde{\tau}^{2}}+\frac{R_{N}}{R_{C}}\left(\frac{d\varphi_{\mathrm{RB}}}{d\widetilde{\tau}}-\frac{d\varphi_{\mathrm{LB}}}{d\widetilde{\tau}}\right)\,,$
(S1)
$\frac{I_{\mathrm{RB}}}{I_{c}}=\sin(\varphi_{\mathrm{RB}})+\frac{d\varphi_{\mathrm{RB}}}{d\widetilde{\tau}}+\beta_{c}\frac{d^{2}\varphi_{\mathrm{RB}}}{d\widetilde{\tau}^{2}}-\frac{R_{N}}{R_{C}}\left(\frac{d\varphi_{\mathrm{RB}}}{d\widetilde{\tau}}-\frac{d\varphi_{\mathrm{LB}}}{d\widetilde{\tau}}\right)\,,$
(S2)

with $\varphi_{\mathrm{LB}}$ and $\varphi_{\mathrm{RB}}$ the phase differences
between junctions $J_{\mathrm{LB}}$ and $J_{\mathrm{RB}}$, respectively,
$\widetilde{\tau}=t/\tau_{J}$ the normalized time, $\tau_{J}=\Phi_{0}/(2\pi
I_{c}R_{N})$, with $\Phi_{0}=h/2e$ the magnetic flux quantum, and
$\beta_{c}=(2e/\hbar)I_{c}R_{N}^{2}C$ the Stewart-McCumber parameter [5]. The
equations are similar to the standard RCSJ model for a single junction, except
of the last term, which introduces the current through the resistor, coupling
the two junctions. This current is a result of the voltage difference between
the two junctions and the coupling resistance. For $R_{C}\to\infty$ the
coupling term goes to zero, leading to two individual junctions (decoupled
system) and for $R_{C}\to 0$ the system is dominated by the coupling term.

##  SIII Shapiro Steps in Three-Terminal Junction Experiments

The differential resistances $R_{\mathrm{LB}}$ and $R_{\mathrm{RB}}$ exposed
to an rf radiation with a frequency of $5.8\text{\,}\mathrm{GHz}$ at 0 dBm
recorded as a function of the applied DC currents are presented in
Supplementary Figures S2(a) and (b). In contrast to the corresponding figure,
which was gained by numerical differentiation, here, the resistance is
directly taken using a lock-in amplifier. In Supplementary Figures S3 the
corresponding measurements at a frequency of $5.8\text{\,}\mathrm{GHz}$ at 0
dBm are shown.

Supplementary Figure S2: Shapiro Step response at 5.8 GHz: (a) shows the
measured differential resistance across the first junction $R_{\mathrm{LB}}$
as a function of the direct current $I_{\mathrm{LB}}$ and $I_{\mathrm{RB}}$
across the junction. (b) shows $R_{\mathrm{RB}}$ for the same current
constellation. Supplementary Figure S3: Shapiro Step response at 8.5 GHz:
$\mathbf{a}$ shows the measured differential resistance across the first
junction $R_{\mathrm{LB}}$ as a function of the direct current
$I_{\mathrm{LB}}$ and $I_{\mathrm{RB}}$ across the junction. $\mathbf{b}$
shows $R_{\mathrm{RB}}$ for the same current constellation.

##  SIV Shapiro Steps in Three-Terminal Junction Simulation

Using the model described in Supplementary Information SII the Shapiro
response was simulated by adding an oscillation contribution
$i_{\mathrm{j,rf}}\sin(2\pi f_{\mathrm{rf}}t)$, $j=\mathrm{LB,RB}$, to the dc
bias currents. The simulated differential resistances $R_{\mathrm{LB}}$ and
$R_{\mathrm{RB}}$ as a function of the normalized voltage drops at a frequency
of $8.5\text{\,}\mathrm{GHz}$ are presented in Supplementary Figures S4(a) and
(b). One finds that the Shapiro response is strong in the corresponding
junctions, while the coupling from the neighboring junction is weak.

Supplementary Figure S4: Shapiro step simulations at 8.5 GHz: (a) Numerically
determined differential resistance $R_{\mathrm{LB}}$ as a function of the
normalized voltage drops $V_{\mathrm{LB}}/V_{0}$ and $V_{\mathrm{RB}}/V_{0}$
at $8.5\text{\,}\mathrm{GHz}$. The blue curves represent the averaged signal
along $V_{\mathrm{LB}}/V_{0}$ and $V_{\mathrm{RB}}/V_{0}$, respectively. (b)
Corresponding map of the differential resistance $R_{\mathrm{RB}}$ with the
blue curves representing the averaged differential resistance along
$V_{\mathrm{LB}}/V_{0}$ and $V_{\mathrm{RB}}/V_{0}$, respectively.

In order to simulate the appearance of the fractional Shapiro steps a non-
sinusoidal current-phase relationship was assumed for the Josephson junction
by including a $\sin(2\varphi)$ contribution. In Supplementary Figures S5(a)
and (b) the respective simulation outcomes $R_{\mathrm{LB}}$ and
$R_{\mathrm{RB}}$ for junctions $J_{\mathrm{LB}}$ and $J_{\mathrm{RB}}$ are
shown as a function of bias currents. One finds that by increasing the $\sin
2\varphi$ contribution fractional steps appear.

Supplementary Figure S5: Simulation of Shapiro maps at
$8.5\text{\,}\mathrm{GHz}$ with $2\phi$ term: (a) differential resistance of
the first junction Shapiro steps as a function of
$I_{\mathrm{LB}},I_{\mathrm{RB}}$ with an equal contribution of a
$\sin{2\phi}$-term in the system, (b) the same for the second junction. (c)
and (d) show the same after doubling the $\sin{2\phi}$ contribution in the
system and (e) and (f) show the same after doubling the contribution of (c)
and (d).

## References

  * Courtois  _et al._ [2008] H. Courtois, M. Meschke, J. T. Peltonen, and J. P. Pekola, Origin of hysteresis in a proximity Josephson junction, Physical Review Letters 101, 067002 (2008).
  * Schüffelgen  _et al._ [2019] P. Schüffelgen, D. Rosenbach, C. Li, T. W. Schmitt, M. Schleenvoigt, A. R. Jalil, S. Schmitt, J. Kölzer, M. Wang, B. Bennemann, U. Parlak, L. Kibkalo, S. Trellenkamp, T. Grap, D. Meertens, M. Luysberg, G. Mussler, E. Berenschot, N. Tas, A. A. Golubov, A. Brinkman, T. Schäpers, and D. Grützmacher, Selective area growth and stencil lithography for in situ fabricated quantum devices, Nature Nanotechnology 14, 825 (2019).
  * Graziano  _et al._ [2020] G. V. Graziano, J. S. Lee, M. Pendharkar, C. J. Palmstrøm, and V. S. Pribiag, Transport studies in a gate-tunable three-terminal Josephson junction, Phys. Rev. B 101, 054510 (2020).
  * Arnault  _et al._ [2021] E. G. Arnault, T. F. Q. Larson, A. Seredinski, L. Zhao, S. Idris, A. McConnell, K. Watanabe, T. Taniguchi, I. Borzenets, F. Amet, and G. Finkelstein, Multiterminal inverse ac Josephson effect, Nano Letters 21, 9668 (2021).
  * Tinkham [2004] M. Tinkham,  _Introduction to Superconductivity_ (Dover Publications, New York, 2004).


#  Topology-free optimal power dispatch for distribution network considering
security constraints and flexible building thermal inertia

Ge Chen, Hongcai Zhang, Ningyi Dai, and Yonghua Song  This work is supported
in part by the Science and Technology Development Fund, Macau SAR (File no.
0137/2019/A3). G. Chen, H. Zhang, N. Dai, and Y. Song are with the State Key
Laboratory of Internet of Things for Smart City and Department of Electrical
and Computer Engineering, University of Macau, Macao, 999078 China (email:
hczhang@um.edu.mo).

###### Abstract

With the increasing integration of distributed PV generation, the distribution
network requires more and more flexibility to achieve the security-constrained
optimal power dispatch. However, the conventional flexibility sources usually
require additional investment cost for equipment. Moreover, involving the
security constraints is very challenging due to the requirements of accurate
network model that may be unavailable in practice. This paper addresses the
aforementioned challenge by proposing a topology-free optimal power dispatch
framework for distribution networks. It utilizes building thermal inertia to
provide flexibility to avoid additional investment. To guarantee the operation
safety, a multi-layer perception (MLP) is trained based on historical
operational data and then reformulated as mixed-integer constraints to replace
the inexplicit original security constraints. Numerical results confirm that
the proposed framework can derive a feasible and optimal strategy without any
topology information.

###### Index Terms:

Demand response, building thermal inertia, optimization, security constraints,
neural networks.

##  I Introduction

In recent years, the investment in PV generation grows rapidly in response to
the challenges of energy crisis and $\text{CO}_{2}$ emission. In 2017, the
cumulative installed PV generation capacity reaches to 130.25 GW in China [1].
Because that PV generation is intermittent, stochastic and often
uncontrollable, its increasing capacity may deteriorate operational efficiency
and security of power distribution networks, e.g. leading to voltage
instability and line current violation. In order to enhance operational
efficiency and security in distribution networks with high penetration of PV
generation, more flexibility sources and smarter power scheduling are required
[2].

Many scholars focus on coordinating distribution networks with distributed
energy storage to improve system efficiency and security. References [3, 4]
combined the battery with grid to reduce the PV curtailment. The security
constraints, such as the line current and bus voltage limitations, were
introduced based on the power flow equations. In [5], the grid was coordinated
with a pumped hydro system to store the excess renewable generation. Reference
[6] treated electric vehicles as distributed energy storage systems and
shifted their charging period to effectively utilize PV power. In [7], an
electricity-gas integrated energy system was proposed to enhance the
utilization of renewable energy and reduce operational cost.

Recently, an increasing number of works utilized building thermal inertia,
such as controlling heating, ventilation, and air conditioning, to shift the
peak loads and provide flexibility for distribution network. Compared with
energy storage systems, thermal inertia does not require any additional
devices, so the investment cost can be reduced. Reference [8] utilized thermal
inertia of water in pipelines to reduce the operation cost of integrated
electricity and district heating systems. In [9], the thermal inertia of
pipeline network was treated as thermal storage systems to realize the optimal
power dispatch for district cooling systems. Reference [10] considered the
flexibility of heat sources, pipelines, and loads to avoid unnecessary
renewable power curtailment. Reference [11] developed an approach to quantify
the demand response flexibility provided by the building thermal inertia. In
references [12, 13], aggregation models for heating, ventilation and air
conditioning systems were established for providing regulation service.

Although the above papers have made great progress on the power dispatch of
distribution network with high-penetration renewable energy, there are still
significant gaps to fill. References [8, 11, 9, 12, 13] utilized the demand-
side sources to avoid the additional investment, but the power flow security
constraints are ignored. References [3, 4, 5, 6, 7] built the security
constraints via the power flow equations, but they require accurate
distribution network model that may be unknown in many practical cases.
Recently, deep reinforcement learning was employed to optimize the power
dispatch [14]. However, unlike the white-box models which can be solved by
off-the-shelf solvers efficiently, those learning-based models are time-
consuming and can not guarantee global optimality.

To overcome the aforementioned challenges, we propose a topology-free optimal
power dispatch framework for distribution networks. This method advances the
published literature in the following aspects:

  1. 1.

We propose a topology-free optimal power dispatch method for distribution
networks. It trains a multi-layer perceptron (MLP) neural network based on
historical data to describe the security constraints of the distribution
network with unknown topology and parameters.

  2. 2.

We utilize buildings’ thermal inertia as distributed energy storage systems to
provide flexibility. This can significantly enhance the power dispatch
performance and reduce PV curtailment.

  3. 3.

A white-box optimal model is proposed by equivalently reformulating the
intractable black-box MLP as mixed-integer linear constraints so that the
global optimal strategy can be found efficiently by off-the-shelf solvers.

The modeling accuracy and dispatch optimality of the proposed framework are
validated by numerical experiments.

The remaining parts are organized as follows. Section II formulates the
optimal power dispatch model. Section III presents the construction and
reformulation of MLP. Section IV conducts simulations and Section V concludes
this paper.

##  II Problem formulation

###  II-A Modeling of building thermal dynamics

A distribution network contains multiple buses and transmission lines. Each
bus can provide electricity for multiple buildings. This paper considers the
flexibility of building thermal inertia in the power dispatch. We aggregate
all buildings in one bus as a big thermal zone to simplify the problem.111Note
that this will not significantly affect the power dispatch performance and
considering multiple buildings under one bus is trivial. We use $i$ and $t$ to
index buses/zones ($i\in\mathcal{I}=\\{1,2,3,...,I\\}$) and time slots
($t\in\mathcal{T}=\\{1,2,3,...,T\\}$). The thermal dynamics of each aggregated
zone can be expressed as follows [15]:

$\displaystyle\begin{split}C_{i}\frac{d\theta_{i,t}^{\text{in}}}{dt}=&\frac{\theta_{t}^{\text{out}}-\theta_{i,t}^{\text{in}}}{R_{i}}+q^{\text{h}}_{i,t}-q^{\text{c}}_{i,t},\quad\forall
i\in\mathcal{I},\quad\forall t\in\mathcal{T},\end{split}$ (1)

where $\theta_{i,t}^{in}$ and $\theta_{t}^{\text{out}}$ are the temperatures
of indoor environment and ambience, in ℃, respectively. Parameter $R_{i}$
represents the thermal resistances from the $i$-th zone to the surroundings,
in ℃/kW. Symbol $q^{\text{h}}_{i,t}$ and $q^{\text{c}}_{i,t}$ are heat load
and cooling supply of the $i$-th space, in kW, respectively. Eq. (1) can be
converted into a linear form by using the finite difference method:

$\displaystyle\begin{split}\theta_{i,t}^{\text{in}}=&\alpha_{i}\theta_{i,t-1}^{\text{in}}+\beta_{i}(q^{\text{h}}_{i,t-1}-q^{\text{c}}_{i,t-1})\\\
&+\gamma_{i}\theta_{t-1}^{\text{out}},\quad\forall
i\in\mathcal{I},\quad\forall t\in\mathcal{T},\end{split}$ (2)

where $\alpha_{i}=1-\Delta t/(R_{i}C_{i})$, $\beta_{i}=\Delta t/C_{i}$ and
$\gamma_{i}=\Delta t/(R_{i}C_{i})$. The indoor temperature should keep in a
proper range to ensure the thermal comforts, as follows:

$\displaystyle\underline{\theta}\leq\theta_{i,t}\leq\overline{\theta},\quad\forall
i\in\mathcal{I},\quad\forall t\in\mathcal{T},$ (3)

where $\underline{\theta}$ and $\overline{\theta}$ denote the lower and upper
temperature bounds of thermal comfortable region.

Cooling systems need to consume power to supply cooling power to each zone.
The corresponding active power consumption $P_{i,t}^{\text{cs}}$ in the $i$-th
space in time slot $t$ is calculated by:

$\displaystyle P_{i,t}^{\text{cs}}=q^{\text{c}}_{i,t}/COP_{i},\quad\forall
i\in\mathcal{I},\quad\forall t\in\mathcal{T},$ (4)

where the parameter $COP$ represents the coefficient of performance of cooling
systems. The reactive power demands of cooling systems are set as zeros
because modern cooling systems usually contain correction devices to maintain
their power factor close to one [16].

###  II-B Modeling of power flow

To guarantee the safe operation, all buses and lines have to operate at an
appropriate working condition. For example, the bus voltages must keep in
[0.9p.u., 1.1p.u.], and the line currents shall not violate allowable upper
bounds.

In order to build the security constraints, we first construct an operation
variable vector $\bm{x}_{t}\in\mathbb{R}^{3I}$ as:

$\displaystyle(\bm{x}_{t})^{\intercal}=[(\bm{P}_{t})^{\intercal},(\bm{Q}_{t})^{\intercal},(\bm{G}_{t}^{\text{PV}})^{\intercal}],\quad\forall
t\in\mathcal{T},$ (5)

where $\bm{P}_{t}$, $\bm{Q}_{t}$ and $\bm{G}_{t}$ denote the vector form of
active, reactive demands and actual used PV generation at all buses in time
slot $t$. The active power demands $\bm{P}_{t}$ can be expressed as:

$\displaystyle\bm{P}_{t}=\bm{P}_{t}^{\text{cs}}+\bm{P}_{t}^{\text{e}},\quad\forall
t\in\mathcal{T},$ (6)

where $\bm{P}_{t}^{\text{e}}$ represents the base load, i.e., the other active
power load in the building except that of the cooling system in $t$.222Note
that $\bm{Q}_{t}$ is the reactive base load because that we assume the cooling
system consume or generate no reactive power.

The actually used PV generation should be less than the available PV
generation $\bm{G}_{\text{av}}^{\text{PV}}$:

$\displaystyle\bm{G}_{t}^{\text{PV}}\leq\bm{G}_{\text{av}}^{\text{PV}},\quad\forall
t\in\mathcal{T}.$ (7)

When the operation variable vector $\bm{x}_{t}$ is determined, the operation
state variables, e.g., currents and voltages, can be calculated based on power
flow equations. Hence, the security constraints of power flows can be
expressed as:

$\displaystyle\bm{s}_{t}(\bm{x}_{t})\in\mathcal{S},\quad\forall
t\in\mathcal{T},$ (8)

where $\bm{s}_{t}$ is the vector containing all operation state variables at
time $t$. Symbol $\mathcal{S}$ denotes the feasible set for safe operation,
which is constructed based on security constraints.

The operation variable vector $\bm{x}_{t}$ also determines the power loss.
Thus, the power loss $l_{t}$ is a function of $\bm{x}_{t}$:

$\displaystyle l_{t}=f(\bm{x}_{t}),\quad\forall t\in\mathcal{T}.$ (9)

The electricity supply from the main grid $G_{t}^{\text{grid}}$ can be
calculated via power balance:

$\displaystyle
G_{t}^{\text{grid}}=G_{t}^{\text{buy}}-G_{t}^{\text{sell}}=\bm{1}^{\intercal}\bm{P}_{t}+l_{t}-\bm{1}^{\intercal}\bm{G}_{t}^{\text{PV}},\quad\forall
t\in\mathcal{T},$ (10) $\displaystyle G_{t}^{\text{buy}}\geq
0,G_{t}^{\text{sell}}\geq 0,\quad\forall t\in\mathcal{T}.$ (11)

If $G_{t}^{\text{grid}}=G_{t}^{\text{buy}}\geq 0$, then the distribution
network purchases power from the main grid; otherwise
$G_{t}^{\text{grid}}=-G_{t}^{\text{sell}}\leq 0$, it sells electricity. The
energy cost $EC_{t}$ can be expressed as:

$\displaystyle
EC_{t}=c^{\text{buy}}G_{t}^{\text{buy}}-c^{\text{sell}}G_{t}^{\text{sell}},\quad\forall
t\in\mathcal{T},$ (12)

where $c^{\text{buy}}$ and $c^{\text{sell}}$ denotes the prices for purchasing
and selling electricity.

###  II-C Formulation of optimal power dispatch problem

We propose a security-constrained optimal power dispatch strategy for
distribution networks. The building thermal inertia is treated as energy
storage systems to provide flexibility. The proposed model can be summarized
as follows:

$\displaystyle\min_{\\{q_{t}^{\text{t}},\forall i\in\mathcal{I},\forall
t\in\mathcal{T}\\}}\quad\sum_{t=1}^{T}EC_{t},\quad\text{s.t.:}\quad\text{Eqs.
(\ref{eqn_T_in})-(\ref{eqn_EC})}.$

Note that the decision variables in the dispatch model include cooling
supplies $q_{t}^{\text{t}}$ and actual PV generation $\bm{G}_{t}^{\text{PV}}$;
The other parts in the operation vector $\bm{x}_{t}$ are treated as given
parameters.

######  Remark 1.

Conventional methods require accurate power flow models to describe function
$\bm{s}_{t}(\bm{x}_{t})$ in Eq. (8) and $f(\bm{x}_{t})$ in Eq. (9). However,
in many practical cases, an accurate power flow model may be unavailable
because of unclear topology information or unknown network parameters.

##  III Topology-free optimal power dispatch

To realize the topology-free power dispatch framework, a multi-layer
perception (MLP) neural network is trained to classify whether the state
variables $\bm{s}_{t}$ stay in the feasible set $\mathcal{S}$ with any given
operation vector $\bm{x}_{t}$. Then, the forward propagation of MLP is exactly
reformulated into mixed-integer constraints that are used to substitute
constraint (8) in problem $\mathcal{P}_{1}$. A linear regression (LR) is also
trained to approximate the power loss function (9). As a result, the optimal
dispatch strategy becomes a mixed integer linear program, which can be
efficiently solved by off-the-shelf solvers.

###  III-A Construction of MLP

Figure 1: A MLP with 3 hidden layers, where $h_{n}^{k}$ denotes the output of
$n$-th neuron in the $k$-th hidden layer.

A MLP is trained to judge the safety of current operation variable
$\bm{x}_{t}$. Fig. 1 shows the typical structure of MLP. It consists of one
input layer, $K$ hidden layers and one output layer. Each layer contains
multiple neurons. Every neuron is composed of a linear mapping and a non-
linear activation function. We use $k$ and $n_{k}$ to index the hidden layers
($k\in\mathcal{K}={1,\cdots,K}$) and the neurons in the $k$-th hidden layer
($n_{k}\in\mathcal{N}_{k}={1,\cdots,N_{k}}$). ReLU is employed as the non-
linear activation function. Then, the forward propagation of MLP is expressed
as

$\displaystyle\bm{x}_{t}=\bm{h}^{0}_{t},\quad\forall t\in\mathcal{T},$ (13)
$\displaystyle\bm{z}^{k}_{t}=\bm{W}^{k}\bm{h}^{k-1}_{t}+\bm{b}^{k},\quad\forall
k\in\mathcal{K},\quad\forall t\in\mathcal{T},$ (14)
$\displaystyle\bm{h}^{k}_{t}=\max(\bm{z}^{k}_{t},0),\quad\forall
k\in\mathcal{K},\quad\forall t\in\mathcal{T},$ (15)

where $h^{k}_{t}\in\mathbb{R}^{N_{k}}$ represents the outputs of layer $k$.
Parameters $\bm{W}^{k}\in\mathbb{R}^{N_{k}\times N_{k-1}}$ and
$\bm{b}^{k}\in\mathbb{R}^{N_{k}}$ denote the weight matrix and bias of layer
$k$, which are parameters to be trained. The outputs of this MLP, i.e.,
$\bm{y}_{t}\in\mathbb{R}^{2}$, is given by:

$\displaystyle\bm{y}_{t}=\bm{W}^{K+1}\bm{h}^{K}_{t}+\bm{b}^{K+1},\quad\forall
t\in\mathcal{T}.$ (16)

The final classification result is determined by comparing the two elements of
$\bm{y}_{t}$. The sample is labeled as “unsafe” if $y_{1,t}>y_{2,t}$;
otherwise, this sample is labeled as “safe”.

###  III-B Construction of LR model

In distribution networks, the effects of power loss on the total energy cost
is not significant because it is usually small compared with the total active
demands. It is not necessary to train a time-consuming model for predicting
power loss with high accuracy. Thus, a simple LR is used to build the mapping
between operation variable $\bm{x}_{t}$ and power loss $l_{t}$, as follows:

$\displaystyle
l_{t}=\bm{W}^{\text{loss}}\bm{x}_{t}+b^{\text{loss}},\quad\forall
t\in\mathcal{T},$ (17)

where $W^{\text{loss}}\in\mathbb{R}^{3I}$ and $b^{\text{loss}}\in\mathbb{R}$
are parameters to be trained.

###  III-C Reformation of topology-free power dispatch model

The constraints (14)-(15) are intractable because of the maximum operator.
However, they can be reformulated into mixed-integer forms by introducing
auxiliary continuous variables $\bm{r}^{k}_{t}$ and binary variable
$\bm{\mu}^{k}_{t}$ based on the big-M method [17]:

$\displaystyle\bm{h}^{k}_{t}-\bm{r}^{k}_{t}=\bm{W}^{k}\bm{h}^{k-1}_{t}+\bm{b}^{k},\quad\forall
k\in\mathcal{K},\quad\forall t\in\mathcal{T},$ (18) $\displaystyle
0\leq\bm{h}^{k}_{t}\leq M\cdot\bm{\mu}^{k}_{t},\quad\forall
k\in\mathcal{K},\quad\forall t\in\mathcal{T},$ (19) $\displaystyle
0\leq\bm{r}^{k}_{t}\leq M\cdot(1-\bm{\mu}^{k}_{t}),\quad\forall
k\in\mathcal{K},\quad\forall t\in\mathcal{T},$ (20)
$\displaystyle\bm{\mu}^{k}_{t}\in\\{0,1\\}^{N_{k}},\quad\forall
k\in\mathcal{K},\quad\forall t\in\mathcal{T},$ (21)

in which $M$ is a big enough constant.

By replacing the security constraint (8) and the power loss function (9) with
the reformulated MLP (18)-(21) and LR model (17), respectively, problem
$\mathcal{P}_{1}$ can be reformulated as:

$\displaystyle\min_{\\{q_{t}^{\text{t}},\forall i\in\mathcal{I},\forall
t\in\mathcal{T}\\}}\quad\sum_{t=1}^{T}EC_{t}$ ($\mathcal{P}_{2}$)
$\displaystyle\begin{array}[]{r@{\quad}r@{}l@{\quad}l}\text{s.t.:}&&\text{Eqs.
(\ref{eqn_T_in})-(\ref{eqn_P}), (\ref{eqn_balance})-(\ref{eqn_input}),
(\ref{eqn_lr})-(\ref{eqn_r4})},\\\ &&y_{t,1}\leq y_{t,2},\quad\forall
t\in\mathcal{T}.\end{array}$ (24)

In $\mathcal{P}_{2}$, constraints (2)-(4) describe the thermal dynamics of
buildings. Constraints (5)-(6) give the relation between original decision
variables and inputs of MLP and LR. Constraints (10)-(12) define the power
balance and calculate the energy cost in objective function. Constraint (13)
defines the inputs of MLP. Constraint (17) calculates the power loss.
Constraints (18)-(21) represent the forward propagation of MLP. The last
constraint $y_{t,1}\leq y_{t,2}$ ensures that the obtained solution satisfies
the security constraints of power flow.

######  Remark 2.

The proposed framework only needs historical data to train the MLP and LR.
Thus, the security constraints and power loss can be introduced to our
optimization without any topology information. Besides, $\mathcal{P}_{2}$ is a
mixed-integer linear program. It can be efficiently solved by the off-the-
shelf solvers. The global optimality can be also guaranteed.

##  IV Case study

We implement two case studies (i.e. heavy and light load cases) based on the
IEEE 33-bus system [18] (shown in Fig. 2) to demonstrate the benefits of the
proposed framework. The heavy-load case is used to validate the optimality and
reliability of the proposed MLP-based framework, while the light-load case is
implemented to show the benefits of treating building thermal inertia as
flexibility source. We assume PV power stations are installed at 5 different
buses (i.e. buses 6, 9, 12, 18 and 30). Note that the information of topology
and parameters of the IEEE 33-bus system are only used to generate training
data and assumed to be unknown in the view of the proposed framework. The
total heat load $\bm{q}_{t}^{h}$, ambient temperature
$\theta_{t}^{\text{out}}$, active based load $\bm{P}_{t}^{e}$, reactive base
load $\bm{Q}_{t}$ and available PV generation $\bm{G}_{\text{av}}^{\text{PV}}$
are illustrated in Fig. 3. All loads in the heavy-load case are the same with
those in Fig.3, while the loads are reduced by half in the light-load case.
The optimization time horizon is 24 hours. The values of other parameters are
summarized in Table I.

Figure 2: Structure of the 33-bus system. PV power stations are installed at
green buses.

Figure 3: (a) Total heat loads, ambient temperature, and (b) active based load, reactive base load and available PV generation. TABLE I: Parameters in case study Parameters | Value | Parameters | Value
---|---|---|---
$C_{i}$ | 1 MWh/℃ | $\underline{\theta}$ | 24℃
$R_{i}$ | 50 ℃/MW | $\overline{\theta}$ | 28℃
$COP_{i}$ | 3.6 | $\eta^{\text{buy}}$ | 0.1122 $/kWh
$\Delta t$ | 1h | $\eta^{\text{sell}}$ | 0.056 $/kWh

To create the labeled dataset for MLP and LR, we randomly generate the
operation vector $\bm{x}_{t}$ as features and input them to Pandapower, a
power system modeling tool in Python, to calculate state variables, i.g.,
nodal voltages and distribution line currents, to obtain the corresponding
labels. If the vector $\bm{x}_{t}$ satisfies all security constraints, this
sample’s label is set as “safe”; otherwise, its label is “unsafe”. In this
paper, the voltage limitation (i.e. the bus voltages are limited in [0.9 p.u.,
1.1 p.u.]) and current limitation (i.e. the line currents should be less than
0.249kA) are considered as our security constraints. Finally, 10,000 samples
are generated, in which 40% of them are “safe”; 60% of them are “unsafe”.
Then, 70% of samples are randomly selected to form the training set and the
remainings are used as the testing set.

A model without considering the security constraints is implemented as
benchmark 1 to show the reliability of the proposed $\mathcal{P}_{2}$.
Benchmark 1 has the same objective function with $\mathcal{P}_{2}$ but its
constraints only contain Eqs. (2-7), (10-12) and (17). To demonstrate the
optimality of the proposed framework, a conventional power-flow-based model is
also presented as benchmark 2. This benchmark is a second-order cone program
and requires exact topology parameters to build the power flow equations in
order to describe security constraints and calculate power losses.

All simulations are performed on an Intel(R) Core(TM) 8700 3.20GHz CPU with 16
GB memory. The MLP and LR are implemented by Pytorch. The corresponding
optimization problem is built by CVXPY and solved by GUROBI.

###  IV-A Heavy load case

####  IV-A1 Reliability

We use a MLP with 2 hidden layer with 8 neurons in each hidden layer to judge
whether the current operation variables $\bm{x}_{t}$ stratifies the security
constraint (8). The accuracy on test set and loss function value during
training are shown in Fig. 4. The model accuracy in the testing set can reach
to 99.09% with a low loss function value, which indicates the high accuracy of
the MLP.

Figure 4: (a) Testing accuracy, and (b) value of loss function during
training.

To demonstrate the reliability of the proposed framework, we calculate the
corresponding true currents and voltages with Pandapower by treating the
solutions of $\mathcal{P}_{2}$ and benchmark 1 as inputs. The corresponding
maximum violations of security constraints among all buses and lines are shown
in Fig. 5. The proposed framework can obtain relatively safe strategy, while
the solution derived by benchmark 1 is much worse. For example, the maximum
voltage violation of the proposed framework almost keeps at zero. The current
violation only appears at $t=8$ and the violation value is 9.99A. The maximum
violation of benchmark 1 are 5.53V ($t=19$) and 66.19A ($t=19$), which are
much larger than those in the proposed framework. The violations in benchmark
1 also occurs much more frequently. These results confirm the proposed
framework can drive a high reliable strategy without any topology information.

Figure 5: Maximum violations of (a) bus voltages, and (b) line currents.

####  IV-A2 Optimality

Fig. 6 compares the hourly costs obtained by $\mathcal{P}_{2}$ and benchmark 2
to demonstrate the optimality of the proposed framework. In most hours, the
hourly cost obtained by the proposed framework is almost the same with that in
benchmark 2. During 12:00-16:00, the pre-cooling in the proposed model is not
so sufficient like that in benchmark 2, so it has to pay more in the next few
hours. Nevertheless, the total cost of the proposed framework is $9120.8,
which is only 1.6% larger than that in the benchmark 2. These results indicate
the high optimality of the proposed framework.

Figure 6: Comparison between the hourly costs obtained by the proposed
framework and benchmark 2.

###  IV-B Light load case

We employ the proposed framework to solve the light-load case to demonstrate
the benefits of employing building thermal inertia to provide flexibility.

####  IV-B1 Indoor temperature

Fig. 7 illustrates the average indoor temperatures obtained by the proposed
$\mathcal{P}_{2}$. In the hours with low PV generation (i.e. 00:00am-11:00am),
all indoor temperature keep at the upper bound of the comfortable region to
minimize the heating effects from outdoor to indoor environments. The
temperature drop occurs at the time with high PV generation (i.e.
11:00am-18:00pm). We call this temperature drop as pre-cooling. With the
decrease of the indoor temperature, more and more power can be stored in the
indoor environments by pre-cooling to compensate the cooling demands in the
next few hours.

Figure 7: Results of indoor temperatures with building thermal inertia.

####  IV-B2 PV curtailments

Fig. 8 compares the results with and without the utilization of building
thermal inertia. The scenario without thermal inertia is implemented by
restricting the indoor temperature at its maximum value $\overline{\theta}$.
Because the indoor environments can serve as thermal storage system, the
excess PV generation can be stored by pre-cooling at the time with high
available PV power. Then, the PV curtailment can be significantly reduced, as
shown in Fig. 8(a). Meanwhile, because of the thermal inertia, the stored
cooling power can satisfy parts of cooling demands in the next few hours
without sunshine (i.e. 18:00pm-24:00pm). Thus, the electricity consumed by
cooling systems in those hours will also decrease. Besides, lower power
consumption means smaller power flow, which can also reduce power loss.

Figure 8: Comparison of (a) PV curtailments and (b) total cost and loss with
and without the utilization of building thermal inertia.

##  V Conclusions

This paper proposes a topology-free framework to realize the security-
constrained optimal power dispatch for distribution networks. The building
thermal inertia is employed as thermal storage systems to store the extra PV
generation. Unlike the conventional power-flow-based approach which requires
accurate topology information, the proposed framework trains a MLP to replace
the security constraints. To achieve the global optima efficiently, the MLP is
reformulated as new constraints with mixed-integer form. Simulation results
show that the classification accuracy of MLP can achieve 99.09%. Then, two
case studies are conducted to demonstrate the feasibility, optimality of the
proposed framework and the benefits of utilizing building thermal inertia. The
results show that the proposed framework can obtain a proper operation
strategy with fewer violations of security constraints. Moreover, the relative
optimality gap is only 1.6%. The results also indicate that employing thermal
inertia can effectively reduce PV curtailment and power losses, and save total
energy cost.

## References

  * [1] N. Tang, Y. Zhang, Y. Niu, and X. Du, “Solar energy curtailment in china: Status quo, reasons and solutions,” Renew. Sust. Energ. Rev., vol. 97, pp. 509–528, 2018.
  * [2] B. Mohandes, M. S. El Moursi, N. Hatziargyriou, and S. El Khatib, “A review of power system flexibility with high penetration of renewables,” IEEE Trans. Power Syst., vol. 34, no. 4, pp. 3140–3155, 2019.
  * [3] Y. Guo, K. Baker, E. Dall’Anese, Z. Hu, and T. H. Summers, “Data-based distributionally robust stochastic optimal power flow—part i: Methodologies,” IEEE Trans. Power Syst., vol. 34, no. 2, pp. 1483–1492, 2018.
  * [4] Y. Guo, K. Baker, E. Dall’Anese, Z. Hu, and T. H. Summers, “Data-based distributionally robust stochastic optimal power flow—part ii: Case studies,” IEEE Trans. Power Syst., vol. 34, no. 2, pp. 1493–1503, 2018\.
  * [5] P. P. Biswas, P. Suganthan, and G. A. Amaratunga, “Optimal power flow solutions incorporating stochastic wind and solar power,” Energy Convers. Manag., vol. 148, pp. 1194–1207, 2017.
  * [6] H. Kikusato, Y. Fujimoto, S.-i. Hanada, D. Isogawa, S. Yoshizawa, H. Ohashi, and Y. Hayashi, “Electric vehicle charging management using auction mechanism for reducing pv curtailment in distribution systems,” IEEE Trans. Sustain. Energy, vol. 11, no. 3, pp. 1394–1403, 2019.
  * [7] Z. Zhang, J. Zhou, Z. Zong, Q. Chen, P. Zhang, and K. Wu, “Development and modelling of a novel electricity-hydrogen energy system based on reversible solid oxide cells and power to gas technology,” Int. J. Hydrog. Energy, vol. 44, no. 52, pp. 28305–28315, 2019.
  * [8] X. Li, W. Li, R. Zhang, T. Jiang, H. Chen, and G. Li, “Collaborative scheduling and flexibility assessment of integrated electricity and district heating systems utilizing thermal inertia of district heating network and aggregated buildings,” Appl. Energ., vol. 258, p. 114021, 2020.
  * [9] G. Chen, B. Yan, H. Zhang, and Y. Song, “Optimal power dispatch for district cooling system considering cooling water transport delay,” in 2020 12th IEEE PES Asia-Pacific Power and Energy Engineering Conference (APPEEC), pp. 1–5, IEEE, 2020.
  * [10] Y. Xue, Z. Li, C. Lin, Q. Guo, and H. Sun, “Coordinated dispatch of integrated electric and district heating systems using heterogeneous decomposition,” IEEE Trans. Sustain. Energy, 2019.
  * [11] L. Zhang, N. Good, and P. Mancarella, “Building-to-grid flexibility: Modelling and assessment metrics for residential demand response from heat pump aggregations,” Appl. Energ., vol. 233, pp. 709–723, 2019.
  * [12] Q. Shi, F. Li, G. Liu, D. Shi, Z. Yi, and Z. Wang, “Thermostatic load control for system frequency regulation considering daily demand profile and progressive recovery,” IEEE Trans. Smart Grid, vol. 10, no. 6, pp. 6259–6270, 2019.
  * [13] R. Adhikari, M. Pipattanasomporn, and S. Rahman, “Heuristic algorithms for aggregated hvac control via smart thermostats for regulation service,” IEEE Trans. Smart Grid, vol. 11, no. 3, pp. 2023–2032, 2019.
  * [14] R. Lu and S. H. Hong, “Incentive-based demand response for smart grid with reinforcement learning and deep neural network,” Appl. Energ., vol. 236, pp. 937–949, 2019.
  * [15] H. Hao, C. D. Corbin, K. Kalsi, and R. G. Pratt, “Transactive control of commercial buildings for demand response,” IEEE Trans. Power Syst., vol. 32, no. 1, pp. 774–783, 2016.
  * [16] S. M. R. Sanhueza, F. L. Tofoli, F. L. de Albuquerque, J. C. de Oliveira, and G. C. Guimaraes, “Analysis and evaluation of residential air conditioners for power system studies,” IEEE Trans. Power Syst., vol. 22, no. 2, pp. 706–716, 2007.
  * [17] M. Fischetti and J. Jo, “Deep neural networks and mixed integer linear optimization,” Constraints, vol. 23, no. 3, pp. 296–309, 2018.
  * [18] M. E. Baran and F. F. Wu, “Network reconfiguration in distribution systems for loss reduction and load balancing,” IEEE Power Eng. Rev., vol. 9, no. 4, pp. 101–102, 1989.