# Time dependent quantum graphs, the geometric phase and knots on closed
Lissajous figures

Uzy Smilansky1 and Gilad Sofer1,2 1 Department of Physics of Complex Systems,
Weizmann Institute of Science, Rehovot 7610001, Israel.
2 Faculty of Mathematics, Technion - Israel Institute of Technology, Haifa
32000, Israel.

###### Abstract

The purpose of the present paper is threefold. The first is to discuss the
time dependent Schrödinger equation on a metric graph with time-dependent
edges, and the proper way to pose the problem so that the corresponding time
evolution is unitary. The second is to compute the geometric phase acquired by
an eigenstate as it adiabatically evolves in time as the edge lengths change.
The third is to study the connection between the geometric phase and the
topology of the cycles in parameter space when they are taken to be Lissajous
knots. We show that the frequencies of the Lissajous knot can be deduced from
the geometric phase.

Dedicated to Michael Berry on his 81’th (+$\frac{2}{3}$) anniversary.

##  1 Introduction

In the present work we address three topics related to time evolution on
quantum graphs whose edge lengths change in time. When combined, these provide
a connection between the geometric phase acquired by an eigenfunction on the
graph, and the topology of the knot created by the closed path in parameter
space which describes the time evolution of the graph. In this section we
first explain the special features and problems each of these topics poses,
and finally state the main results.

### The time dependent Schrödinger equation on a metric graph

A stationary quantum graph consists of a set of vertices connected by edges of
prescribed lengths. On each edge, the Schrödinger operator is the one
dimensional Laplacian. The graph is supplemented with appropriate boundary
conditions at the vertices, which ensure that the operator is self-adjoint and
so the resulting time evolution is unitary.

To illustrate the difficulties which arise when the edge lengths vary in time,
we start by considering the simplest example. Namely, a graph that consists of
two vertices connected by an edge of length $L(t/T)$. Here, $T^{-1}$ measures
the rate of change of the edge lengths.

The evolution is dictated by the time-dependent Schrödinger equation:

$\displaystyle i$ $\displaystyle\hbar\frac{\partial\phi(x,t)}{\partial
t}=-\frac{\hbar^{2}}{2m}\frac{\partial^{2}\phi(x,t)}{\partial x^{2}}\ \ ,\ \ \
x\in\left[-\frac{1}{2}L(t/T),\frac{1}{2}L(t/T)\right].$ (1)

At the two vertices we impose the standard Neumann boundary condition:

$\displaystyle\ \ \ \frac{\partial\phi(x,t)}{\partial x}=0\ \ \ {\rm at\ }\
|\xi|=\frac{1}{2}L(t/T),\ \forall t\geq 0\ .$ (2)

The operator on the right hand side of (1) involves only derivatives with
respect to $x$. The time dependence is due to the boundary conditions, which
are imposed at the time dependent boundary. If one considers the time $t$ as a
parameter, the operator on the right is referred to as the instantaneous
Hamiltonian.

A straightforward computation gives the following expression for the time
derivative of the squared $L^{2}$ norm of $\phi(x,t)$:

$\frac{{\rm d}\ }{{\rm d}t}\|\phi\|^{2}=\frac{{\rm d}\ }{{\rm
d}t}\int_{-{\frac{1}{2}L(t/T)}}^{\frac{1}{2}L(t/T)}|\phi(x,t)|^{2}{\rm d}x\
=\frac{1}{T}\frac{{\rm d}L(t/T)}{{\rm d}t}\left|\phi(L(t/T),t)\right|^{2}\ .$
(3)

This rather unexpected change of the $L^{2}$ norm in time is due to the choice
of Neumann boundary condition. Recalling the standard interpretation of
$|\phi(x,t)|^{2}dx$ as the probability for finding the particle in the
interval $(x,x+dx)$ at time $t$, the above result implies that there exists a
probability flux which equals to the product of the density at the boundary
and its speed. Hence, the probability to find the particle within the
boundaries changes in time. This implies that the evolution dictated by (1) is
not unitary. In this context, it is worth noting:

1\. The instantaneous $L^{2}$ norm is constant in the adiabatic limit
$T\rightarrow\infty$.

2\. For the Dirichlet boundary condition ($\phi(x,t)=0$), the instantaneous
$L^{2}$ norm is constant. Equation (1) shows that the change in norm arises
also for the Robin boundary condition:

$\frac{\partial\phi}{\partial x}(x,t)=C\phi(x,t)\ \ C\in\mathbb{R}.$ (4)

In Section 2 we return to this problem, and derive the boundary conditions
which render the time evolution unitary, ensuring probability conservation. We
follow previous studies which treat time dependent domain boundaries in
different contexts. The time dependent Schrödinger equation on domains in
$\mathbb{R}^{n}$ with time dependent boundary was studied in [1, 2, 3]. The
propogation of electromagnetic waves in vibrating cavities is a closely
related problem studied in [4]. The present work can be considered as an
extension of these studies to quantum graphs. We shall note that the time
dependent Schrödinger equation on a star graph was studied in [5]. However,
the vertex boundary conditions considered there are not consistent with the
ones presented here, and do not address the problem discussed above.

### The geometric phase

In [6], Michael Berry considered quantum mechanical systems with a parameter
dependent potential which changes adiabatically in time. The system is
initially prepared at an eigenstate $|\phi_{n}(0)\rangle$ of the instantaneous
Hamiltonian at time $t=0$. It then evolves according to the Schrödinger
equation along a closed path in parameter space. Berry then shows that in the
adiabatic limit, the wave function at the end of a cycle $|\phi_{n}(T)\rangle$
returns to the initial function, up to a nontrivial phase:

$|\phi(T)\rangle=e^{-\frac{i}{\hbar}\int_{0}^{T}E_{n}(t^{\prime})dt^{\prime}}e^{-\int_{0}^{T}\langle\phi_{n}|\frac{\partial\phi_{n}}{\partial
t}\rangle dt^{\prime}}|\phi_{n}(0)\rangle.$ (5)

The first term in the phase is known as the dynamical phase, and is a direct
generalization of the phase acquired in the time independent Schrödinger
equation. The second term is known as the geometric phase. Unlike the
dynamical phase, the geometric phase depends only on the path transversed in
the space of parameters. The geometric phase has been the subject of many
studies since its discovery.

Unlike the work by Berry, where the time dependence enters through the
changing potential, in this work we are interested in metric graphs where the
edge lengths depend on time. In other words, it is not just the operator which
varies in time, but the Hilbert space itself. For such a system, the original
derivation of the geometric phase cannot be directly applied. Once the problem
is properly reformulated in Section 2, the effect of the moving boundary is
expressed in terms of a gauge field which must be included in order to render
the Schödinger operator self-adjoint. This is analogous to the derivation for
slowly varying time dependent domains presented in [3]. The geometric phase
can then be properly defined.

### The parameter space, Lissajous figures and knots

Figure 1: Three Lissajous knots with $\vec{\nu}=(2,3,5)$. Left - unknot with
$\vec{\alpha}=(0,\frac{\pi}{15},\frac{9\pi}{11})$; Middle - type 9/2 knot with
$\vec{\alpha}=(0,\frac{\pi}{15},\frac{6\pi}{11})$; Right - type 7/2 knot with
$\vec{\alpha}=(0,\frac{\pi}{15},\frac{5\pi}{11})$ (as seen in [7]).

In the present case, the space of parameters corresponds to the edge lengths
$L_{e}(t/T)$ for each edge $e$. Assuming that the graph has $E$ edges, we get
an $E$-dimensional space of parameters. For our model, we consider periodic
edge lengths described by simple trigonometric functions with coprime
frequencies $\\{\nu_{e}\\}_{e=1}^{E}$, which trace Lissajous figures in
parameter space:

$L_{e}(t/T)=\bar{L}_{e}+\rho_{e}\cos(2\pi\nu_{e}t/T+\alpha_{e})\ \ \ ;\ \ \
0\leq\rho_{e}<\bar{L}_{e}\ ,\ \nu_{e}\in{\mathbb{N}}\ ,\ \ {\rm for}\ \ 1\leq
e\leq E\ ,$ (6)

where the parameters $\alpha_{e}$ are phase shifts. To remove the trivial
freedom of choice of an initial point on the figure, one may assume
$\alpha_{1}=0$. By choosing $\rho_{e}=0$, one may naturally include time
independent edges. In the examples that are discussed in the following, we
usually focus on the case where $E=3$, but the formalism is valid for any
choice of $E$.

The curves (6) are simple and closed, and in the three dimensional case which
is studied here, correspond to a family of knots known as Lissajous knots.
These knots were studied by several authors, see e.g. [8, 7]. Nevertheless,
their classification in terms of the frequencies and phase-shift parameters is
still incomplete.

Due to their natural symmetries, not every knot can be described by a
Lissajous knot. For instance, all torus knots are not Lissajous, and in
particular the trefoil knot is not Lissajous.

The main quantity which determines the type of the Lissajous knot is the
vector of coprime frequencies $\vec{\nu}:=(\nu_{1},\nu_{2},\nu_{3})$. Up to
ambient isotopy, there are only finitely many Lissajous knots of any given
frequency vector. The class of the knot is then determined by the vector of
phases $\vec{\alpha}:=(\alpha_{1},\alpha_{2},\alpha_{3})$, as discussed in [8,
7]. Moreover, given $\vec{\nu}$, there exists always a set of phase shifts for
which the corresponding figure is an unknot.

The structure of the paper is as follows. In Section 2 we present the modified
version of the Schrödinger equation for a time dependent quantum graph which
guarantees self-adjointness of the operator and thus unitary evolution
(Theorem 1). In particular, we derive the secular equation whose zeros give
the instantaneous spectrum at any time $t$. In Section 3, the geometric phase
is derived by studying the time evolution in the adiabatic limit. The
resulting expression is discussed and it is shown that the frequency vector of
the Lissajous knot may be determined from the geometric phase (Theorems 2 and
3). Finally, in Section 4 we demonstrate the theory by studying the example of
a star graph. In particular, we improve the estimates from Theorems 2 and 3
and show that in some cases, the knot can be completely recovered from the
behavior of the geometric phase. The final section is devoted to further
discussion and summary.

##  2 Quantum graphs with time dependent edge lengths and Neumann boundary
conditions

This section consists of two parts, in which the time dependent Schrödinger
equation for quantum graphs with time dependent edges is formulated and
discussed. The first subsection addresses the problem in the simple case of a
single interval. It serves as a primer for the treatment of a general graph
which is presented in the second subsection.

###  2.1 A time dependent interval

Before proceeding further, it is worthwhile to rewrite the Schrödinger
equation (1) in dimensionless form. Observing that $\frac{\hbar}{2m}$ has the
dimension of $\frac{\text{length}^{2}}{\text{time}}$, then scaling $x$ by $L$
and $t$ by $T$ results in the dimensionless quantity
$\frac{L^{2}}{T}\cdot\frac{\hbar}{2m}$. Using the standard convention in which
$\frac{\hbar}{2m}$ takes the numerical value $1$, we remain with a
dimensionless $\frac{L^{2}}{T}$. This gives the following Schrödinger
equation:

$\displaystyle i$ $\displaystyle\frac{\partial\phi(x,t)}{\partial
t}=-\frac{\partial^{2}\phi(x,t)}{\partial x^{2}}\ \ ,\ \ \
x\in\left[-\frac{1}{2}L(t/T),\frac{1}{2}L(t/T)\right].$ (7)

Following [1, 2, 3], the Schrödinger equation is written in terms of the
scaled variables

$\tau=\frac{t}{T}\ ,\ \xi=\frac{x}{L(t/T)}\ .$ (8)

Using

$\displaystyle\frac{\partial\ }{\partial t}$ $\displaystyle=$
$\displaystyle\frac{1}{T}\left(\frac{\partial\ }{\partial\tau}-\xi\
\frac{\dot{L}}{L}\frac{\partial\ }{\partial\xi}\right)\ \ ;\ \ \frac{\partial\
}{\partial x}=\frac{1}{L}\frac{\partial\ }{\partial\xi}\ \ \ \ {\rm and}\ \ \
\ \ \dot{L}=\frac{\rm{d}L(\tau)}{\rm{d}\tau}\ ,$ (9)

the wave function $\phi(x,t)$ is mapped into
$\tilde{\omega}(\xi,\tau)=\phi(x(\xi,\tau),t(\xi,\tau))$, which satisfies the
equation

$\displaystyle\frac{i}{T}$
$\displaystyle\frac{\partial\tilde{\omega}(\xi,\tau)}{\partial\tau}=\frac{1}{L(\tau)^{2}}\left(-\frac{\partial^{2}\tilde{\omega}(\xi,\tau)}{\partial\xi^{2}}+i\xi\frac{L(\tau)\dot{L}(\tau)}{T}\frac{\partial\tilde{\omega}(\xi,\tau)}{\partial\xi}\right)\
\ ,\ -\frac{1}{2}\leq\xi\leq\frac{1}{2}\ .$ (10)

The adiabatic limit is achieved when $T\rightarrow\infty$, and in the modified
Schrödinger equation, the parameter $\frac{1}{T}$ replaces $\hbar$. In other
words, the adiabatic limit is the analogue of the semi-classical limit.

The Neumann boundary condition remains the same in the new coordinates:

$\frac{\partial\tilde{\omega}(\xi,\tau)}{\partial\xi}=0\ \ {\rm at\ }\
|\xi|=\frac{1}{2}\ .$ (11)

The resulting instantaneous operator is not self-adjoint, as can be checked by
direct computation. This is another manifestation of the non-unitary time
evolution of (1) observed in the introduction.

Substituting
$\tilde{\omega}(\xi,\tau)=(L(\tau))^{-\frac{1}{2}}\omega(\xi,\tau)$ in (10),
the equation for $\omega(\xi,\tau)$ reads

$\displaystyle\frac{i}{T}\frac{\partial\
}{\partial\tau}\omega(\xi,\tau)=\frac{1}{L(\tau)^{2}}\left[-\left(\frac{\partial\
}{\partial\xi}-i\xi\frac{\dot{L}(\tau)L(\tau)}{2T}\right)^{2}+\xi^{2}\left(\frac{\dot{L}(\tau)L(\tau)}{2T}\right)^{2}\right]\omega(\xi,\tau)\
.$ (12)

In the new Schrödinger equation, the Laplacian appears as a magnetic Laplacian
with $\xi\frac{\dot{L}L}{2T}$ playing the role of a vector potential. The
Schrödinger operator now also includes a harmonic potential.

To render the magnetic operator above self-adjoint, we replace the original
Neumann condition with the magnetic Neumann condition:

$\left(\frac{\partial\ }{\partial\xi}-i\xi\frac{\dot{L}L}{2T}\right)\
\omega(\xi,\tau)=0\ \ {\rm at\ }\ |\xi|=\frac{1}{2}.$ (13)

One may verify that under the magnetic boundary condition, our operator is now
self-adjoint, and the resulting time evolution is thus unitary.

The following alternative explanation for the magnetic boundary condition was
suggested by Michael Berry in [12]. Instead of moving to the scaled
coordinates $(\xi,\tau)$, one may account for the change in probability in (3)
by introducing a _modified_ Robin boundary condition:

$\frac{\partial\psi(x,t)}{\partial x}=C\psi(x,t)\ \ {\rm at\ }\
|x|=\frac{1}{2}L(t/T).$ (14)

While the usual Robin condition corresponds to $C\in\mathbb{R}$, a
straightforward computation shows that the requirement $\frac{{\rm d}\ }{{\rm
d}t}\|\phi\|^{2}=0$ gives the condition $Im(C)=\frac{\dot{L}}{4T}$. We may
thus solve the problem by suggesting the vertex condition

$\frac{\partial\psi(x,t)}{\partial x}=\frac{i\dot{L}}{4T}\psi(x,t)\ \ {\rm at\
}\ |x|=\frac{1}{2}L(t/T).$ (15)

And indeed, when translating condition (15) to $\omega(\xi,\tau)$, this is
exactly the magnetic Neumann condition (13).

Introducing the gauge transformation

$\omega(\xi,\tau)=e^{-\frac{i}{2}\theta(\tau)(\xi^{2}-\frac{1}{4})}g(\xi,\tau)\
,$ (16)

one may eliminate the magnetic Laplacian in Equation (12) if the phase is
chosen to be $\theta(\tau)=-\frac{\dot{L}L}{2T}$. Note that in the units used
in this work, $\theta$ is dimensionless.

The Schrödinger equation is now given by

$\frac{i}{T}\frac{\partial
g(\xi,\tau)}{\partial\tau}=\frac{1}{L^{2}}\left[-\frac{\partial^{2}g(\xi,\tau)}{\partial\xi^{2}}+\frac{\ddot{L}}{L}\left(\frac{{L^{2}}}{2T}\right)^{2}\
\xi^{2}g(\xi,\tau)\right]\ \ ,\ |\xi|\leq\frac{1}{2}\ .$ (17)

This equation was presented in [2, 5]. We derived it here as a prelude to the
treatment of the time dependent Schrödinger operator on a graph. Note that
also here the factor $\frac{1}{T}$ multiplies the $\tau$ derivative, showing
again that the adiabatic limit is the analogue of the semi-classical limit.

The instantaneous Schrödinger operator includes a one-dimensional quadratic
potential which is bounded in the interval $[-\frac{1}{2},\frac{1}{2}]$. Its
spectrum is discrete and the eigenfunctions can be written in terms of
parabolic cylindrical functions, as done in [5].

To distinguish between quantities which genuinely depend on two variables
$(x,t)$ or $(\xi,\tau)$ from those where $t$ or $\tau$ serve as parameters (as
in the instantaneous representation), we denote the latter by using semicolon
instead of comma and write them as $(x;t)$ and $(\xi;\tau)$. The instantaneous
eigenfunctions of (12) can thus be written as

$\omega^{(n)}(\xi;\tau)=e^{-\frac{i}{2}\theta(\tau)(\xi^{2}-\frac{1}{4})}g^{(n)}(\xi;\tau)\
.$ (18)

The phase $\theta(\tau)$ is determined solely by the time dependent length
$L(\tau)$. $\left\\{\omega^{(n)}(\xi;\tau)\right\\}_{n=1}^{\infty}$ are the
instantaneous eigenfunctions which enable us to approach the adiabatic limit,
by using the standard approach in [6], which is formulated for systems where
the domains supporting the eigenfunctions are time independent.

###  2.2 Compact, connected graphs with time dependent edge lengths

Consider now the metric graph $\mathcal{G}=(\mathcal{V,E,L})$ , where
$\mathcal{V}$ is the vertex set with $V:=|\mathcal{V}|$ and $\mathcal{E}$ is
the edge set with $E:=|\mathcal{E}|$ (both assumed finite). Moreover,
$\mathcal{L}:[0,T]\rightarrow\mathbb{R}^{E}$ is a family of edge lengths
parametrized as $\left\\{L_{e}(t/T)\right\\}_{e\in\mathcal{E}}$ with
$L_{e}(t/T)$ positive and sufficiently smooth. A coordinate $x_{e}$ with
$x_{e}\in\left[-\frac{1}{2}L_{e}(t/T),\frac{1}{2}L_{e}(t/T)\right]$ is
assigned to every edge in ${\mathcal{E}}$. We denote the set of edges
connected to the vertex $v$ by $S_{v}$ and the degree of the vertex by
$d_{v}:=|S_{v}|$.

The model which will be discussed in detail consists of edges which trace
Lissajous figures in the space of edge-lengths which was introduced in (6).
The frequency parameters $\nu_{e}$, which are relatively prime, and the phases
$\alpha_{e}\in[0,2\pi)$ determine the line traced in parameter space as $t$
varies.

Given an initial eigenstate $\psi_{n}(x,0)$, we assume that the time evolution
is such that at no time is the associated eigenvalue $\lambda_{n}(t)$
degenerate. This assumption is generic and satisfied almost surely for an
arbitrary choice of the length, frequency and phase parameters. This
assumption is necessary for the assumption of adiabatic evolution to hold,
which is the setting in which it is meaningful to discuss the geometric phase.

The time dependent Schrödinger operator for this model consists of the direct
sum of $E$ one-dimensional Laplacians attached to the edges
${\bf\Delta}=\bigoplus_{e\in{\mathcal{E}}}\frac{\partial^{2}\ }{\partial
x_{e}^{2}}$:

$i\frac{\partial\ }{\partial t}\Psi[{\bf x},t]=-{\bf\Delta}\Psi[{\bf x},t],\ \
\ \ \ {\bf x}=\\{x_{e}\\}_{e\in{\mathcal{E}}},\ $ (19)

with

$\Psi[{\bf x},t]=\left\\{\psi_{e}(x_{e},t)\right\\}_{e\in{\mathcal{E}}}\ ,\ \
\ 0\leq t\leq T.$ (20)

We have shown in the preceding section that for the graph which consists of a
single interval, the standard boundary conditions (with the exception of the
Dirichlet boundary condition) which are valid for the time independent case
should be modified when the length is time dependent. A similar computation
shows that the same is true for arbitrary graphs. In particular, one should
replace the commonly used Neumann-Kirchhoff vertex condition:

$\displaystyle\psi_{e}(x_{e}=v;t)=\psi_{e^{\prime}}(x_{e^{\prime}}=v;t)\ \ \ \
\forall e,e^{\prime}\in S_{v}\ \ ,$ (21) $\displaystyle\sum_{e\in
S_{v}}\frac{\partial\psi_{e}(v;t)}{\partial x_{e}}=0\ \ ({\rm derivatives\
directed\ away\ from\ the\ vertex}).$

The alternative formulation should provide a self-adjoint instantaneous
Schrödinger operator, which, in the limit of constant edge lengths, should
converge to the standard Neumann-Kirchhoff vertex condition.

Following the example of the interval, it is suggested that the proper way to
represent the problem is by introducing a new set of coordinates
$(\xi_{e}=x_{e}/L_{e}(t/T);\tau=t/T)$, in which the edge lengths are time
independent, and then scaling the corresponding eigenfunction in order to
obtain a self-adjoint magnetic operator. That is,

$\displaystyle\psi_{e}(x,t)=\tilde{\omega}_{e}(\xi_{e}(x,t),\tau(x,t))\ \ ;\ \
\tilde{\omega}_{e}(\xi_{e},\tau)=L_{e}(\tau)^{-1/2}{\omega}_{e}(\xi_{e},\tau).$
(22)

The new Schrödinger equation now reads

$\frac{i}{T}\frac{\partial\
}{\partial\tau}\omega_{e}(\xi_{e};\tau)=\frac{1}{L_{e}(\tau)^{2}}\left[-\left(\frac{\partial\
}{\partial\xi_{e}}-i\xi_{e}\frac{\dot{L}_{e}(\tau)L_{e}(\tau)}{2T}\right)^{2}+\xi_{e}^{2}\left(\frac{\dot{L}_{e}(\tau)L_{e}(\tau)}{2T}\right)^{2}\right]\omega_{e}(\xi_{e};\tau)\\\
.$ (23)

The vertex boundary conditions which correspond to the magnetic Schrödinger
operator above are given by

$\displaystyle
L_{e}^{-1/2}(\tau)\omega_{e}(v)=L_{e^{\prime}}^{-1/2}(\tau)\omega_{e^{\prime}}(v),\forall
e,e^{\prime}\in S_{v},\ \ $ (24) $\displaystyle\sum_{e\in
S_{v}}\frac{1}{L_{e}^{3/2}(\tau)}\left(\frac{\partial}{\partial\xi_{e}}-i\xi_{e}\frac{\dot{L}_{e}L_{e}}{2T}\right)\omega_{e}(\xi_{e};\tau)=0\
\ {\rm at\ }\ |\xi_{e}|=\frac{1}{2}.$ (25)

If we switch back to our original time dependent graph, the condition above
translates into the following magnetic vertex condition:

$\displaystyle\psi_{e}(v)=\psi_{e^{\prime}}(v),\forall e,e^{\prime}\in S_{v},\
$ (26) $\displaystyle\sum_{e\in S_{v}}\left(\frac{\partial}{\partial
x_{e}}-i\cdot\text{sgn}(x_{e})\frac{\dot{L}_{e}}{4T}\right)\psi_{e}(x,t)=0.$
(27)

One can readily verify that our operator is now self-adjoint and that the
corresponding time evolution is unitary. We thus conclude that for the problem
of the time dependent graph to be well posed, the usual Neumann-Kirchhoff
vertex condition should be replaced with the magnetic vertex condition (26),
(27). Moreover, note that in the limit $\dot{L}(\tau)\rightarrow 0$, one
recovers the standard Neumann-Kirchhoff vertex condition of a stationary
quantum graph. We state this as the main theorem for this section:

######  Theorem 1.

Let $\mathcal{G}=(\mathcal{V,E,L})$ be a time dependent quantum graph. Denote
by $\mathcal{G}_{0}$ the stationary metric graph obtained by taking $L_{e}=1$
for all $e\in\mathcal{E}$. Consider the time dependent Schrödinger equation
(23) on $\mathcal{G}_{0}$ under the magnetic vertex condition (24), and denote
its solution on each edge $e$ by $\omega_{e}(\xi_{e},\tau)$.

Then the solution to the time dependent Schrödinger equation (7) on
$\mathcal{G}$ with the magnetic Neumann-Kirchhoff vertex condition (26), (27)
is given by

$\psi_{e}(x,t)=L_{e}^{-1/2}(t/T)\omega_{e}\left(\frac{x_{e}}{L_{e}(t/T)},\frac{t}{T}\right).$
(28)

The resulting time evolution is unitary.

We now go back to our representation of the problem using the stationary graph
(24), (25). It turns out that using this representation, the problem can be
solved more easily.

Note that the total magnetic flux induced by the magnetic potential of the
instantaneous Hamiltonian through every cycle of the graph is equal to zero.
Thus, by following the method presented in [9], one can eliminate the magnetic
term by applying the gauge transformation

$\omega_{e}(\xi_{e};\tau)=e^{-\frac{i}{2}\theta_{e}(\tau)(\xi_{e}^{2}-\frac{1}{4})}g_{e}(\xi_{e};\tau)\
,\ \ \theta_{e}=-\frac{{\dot{L}_{e}}L_{e}}{2T}.$ (29)

The Schr̈odinger equation then reads:

$\frac{i}{T}\frac{\partial
g_{e}(\xi_{e};\tau)}{\partial\tau}=\frac{1}{L_{e}^{2}(\tau)}\left[-\frac{\partial^{2}g_{e}(\xi_{e};\tau)}{\partial\xi_{e}^{2}}+\frac{\ddot{L_{e}}}{L_{e}}\left(\frac{{L_{e}^{2}}}{2T}\right)^{2}\
\xi_{e}^{2}g_{e}(\xi_{e};\tau)\right]\ \ ,\ |\xi_{e}|\leq\frac{1}{2}\ ,\
\forall e\in{\mathcal{E}},$ (30)

where the vertex conditions are now

$\displaystyle
L_{e}^{-1/2}(\tau)g_{e}(v;\tau)=L_{e^{\prime}}^{-1/2}(\tau)g_{e^{\prime}}(v;\tau),\forall
e,e^{\prime}\in S_{v},\ \ $ (31) $\displaystyle\sum_{e\in
S_{v}}\frac{1}{L_{e}^{3/2}(\tau)}\frac{\partial
g_{e}(v;\tau)}{\partial\xi_{e}}=0.\ \ $ (32)

In the adiabatic limit $T\rightarrow\infty$, the harmonic potential in (30)
(which is bounded for $|\xi_{e}|\leq\frac{1}{2}$) is negligible, and we remain
with a weighted free Laplacian. The instantaneous equation for the
eigenfunctions $\\{g_{e}(\xi_{e};\tau)\\}_{e=1}^{E}$ can thus be solved via
the standard scattering approach presented in [10]. If on each edge we write
our eigenfunctions as

$g_{e}(\xi_{e};\tau)=a_{e}(\tau)e^{ikL_{e}(\tau)\xi_{e}}+a_{\hat{e}}(\tau)e^{ikL_{e}(\tau)(1-\xi_{e})},$
(33)

then a straightforward computation shows that the following linear relation
holds due to the vertex conditions:

$a_{e}=\frac{2L_{e}^{1/2}}{d_{v}}\sum_{j\in
S_{v}}L_{j}^{-1/2}e^{ikL_{j}}a_{\hat{j}}-a_{\hat{e}}.$ (34)

We can equivalently write this in terms of the bond scattering matrix; We say
that the directed edge $j^{\prime}$ follows the directed edge $j$ if the
starting vertex of $j^{\prime}$ is the end vertex of $j$. Then by the above,
we see that $k^{2}>0$ is an eigenvalue of the instantaneous problem if and
only if it is a root of the secular equation:

$\zeta(k;\tau):=\det\left(I-S(\tau)D\left(k;\tau\right)\right)=0.$ (35)

Here, $D(k;\tau)$ is the diagonal matrix of size $2E$:

$D\left(k;\tau\right)_{j,j^{\prime}}=e^{ikL_{j}}\delta_{jj^{\prime}}\ ,$ (36)

and $S(\tau)$ is the instantaneous bond scattering matrix:

$\displaystyle
S(\tau)_{j^{\prime},j}=\cases{\frac{2}{d_{v}}\left(\frac{L_{j^{\prime}}}{L_{j}}\right)^{1/2}-1&$j^{\prime}=\hat{j}$\cr\frac{2}{d_{v}}\left(\frac{L_{j^{\prime}}}{L_{j}}\right)^{1/2}&$j^{\prime}$\text{
follows $j$, $j^{\prime}\neq\hat{j}$}\cr 0&Otherwise\cr}$ (37)

Thus, one may obtain the instantaneous spectrum in the adiabatic limit by
solving the instantaneous secular equation for each $\tau$. Note that unlike
many common cases, here the bond scattering matrix $S(\tau)$ is not unitary,
since the edges of the graph were scaled, giving a Laplacian with different
(length dependent) weights on each edge.

##  3 The geometric phase

We follow the approach suggested by Berry in [6]. As before, we denote the
eigenfunctions of the instantaneous Hamiltonian by $\omega^{(n)}(\xi;\tau)$
with eigenvalues $\lambda^{(n)}(\tau)$. Given our initial eigenstate, we
suggest the following solution for the state at time $\tau$:

$\omega(\xi,\tau)=\sum_{k\in\mathbb{N}}a_{k}(\tau)\omega^{(k)}(\xi;\tau).$
(38)

Plugging this into the Schrödinger equation (30) we get

$i\sum_{k\in\mathbb{N}}\left(\dot{a}_{k}(\tau)\omega^{(k)}(\xi;\tau)+a_{k}(\tau)\dot{\omega}^{(k)}(\xi;\tau)\right)=\sum_{k\in\mathbb{N}}\lambda^{(k)}(\tau)a_{k}(\tau)\omega^{(k)}(\xi;\tau).$
(39)

Recall that we assume that the system evolves adiabatically, and that the
eigenstates are non-degenerate for all $\tau$. In first order perturbation
theory, we get the following expression for the coefficients:

$i\dot{a}_{n}(\tau)\approx\lambda^{(n)}(\tau)a_{n}(\tau)-i\langle\omega^{(n)}|\dot{\omega}^{(n)}\rangle
a_{n}(\tau).$ (40)

We thus conclude that just as in the original derivation by Berry, the
expression for the rate of change of the geometric phase is given by

$\dot{\gamma}^{(n)}(\tau)=\langle\omega^{(n)}|\dot{\omega}^{(n)}\rangle.$ (41)

We now return to the example of a single interval and compute the
corresponding geometric phase. We then move to the general graph setting.

###  3.1 The geometric phase for a time dependent interval

Using (29) and (40), we get the following formula for the rate of change of
the geometric phase for an interval:

$\dot{\gamma}^{(n)}(\tau)=-\frac{1}{2}\dot{\theta}(\tau)\int_{0}^{1}\left(\xi^{2}-\frac{1}{4}\right)|g^{(n)}(\xi;\tau)|^{2}{\rm
d}\xi\ ,$ (42)

where we have used the fact that $g^{(n)}(\xi;\tau)$ are real and normalized
for all $\tau$.

Plugging in the expression for $\theta(\tau)$ from (29) and using the
definition of the expectation value one obtains:

$\dot{\gamma}^{(n)}(\tau)=\frac{1}{8T}\frac{{\rm d}^{2}L^{2}(\tau)}{{\rm
d}{\tau}^{2}}\ \ \left\langle\xi^{2}-\frac{1}{4}\right\rangle_{\tau}^{(n)}.$
(43)

Integrating the rate over $\tau$ one then gets:

$\Gamma^{(n)}:=T[\gamma^{(n)}(1)-\gamma^{(n)}(0)]=\frac{1}{8}\int_{0}^{1}{\rm
d}\tau\frac{{\rm d}^{2}L^{2}(\tau)}{{\rm d}{\tau}^{2}}\ \
\left\langle\xi^{2}-\frac{1}{4}\right\rangle_{\tau}^{(n)}\ .$ (44)

Consider the limit $T\rightarrow\infty$: In this limit, the instantaneous
quadratic potential which appears in (17) can be neglected since
$|\xi|<\frac{1}{2}$. Hence, the spectrum and normalized eigenfunction do not
depend on $T$, and therefore $\Gamma^{(n)}$ defined above depends only on the
closed path in parameter space $L(\tau)$. In this limit, the instantaneous
eigenfunctions are

$g^{(n)}(\xi;\tau)=\frac{1}{\sqrt{2}}\cos(\pi n\xi),\\\
\left\langle\xi^{2}-\frac{1}{4}\right\rangle_{\tau}^{(n)}=\frac{1}{4\sqrt{2}}\left(\frac{\left(-1\right)^{n}}{\pi^{2}n^{2}}-\frac{1}{3}\right),$
(45)

and the corresponding geometric phase is given by

$\Gamma^{(n)}=\frac{1}{32\sqrt{2}}\left(\frac{\left(-1\right)^{n}}{\pi^{2}n^{2}}-\frac{1}{3}\right)\left.\frac{{\rm
d}L^{2}(\tau)}{{\rm d}\tau}\right|_{\tau=0}^{\tau=1}=0\ .$ (46)

That is, in a one-dimensional parameter space, the geometric phase is trivial.

Figure (2) displays the accumulated phase $\Gamma(\tau)$ for $\tau\in[0,1]$
for the interval. As $\tau\rightarrow 1$ (which corresponds to completing a
cycle in parameter phase), the phase goes to zero, as expected.

Figure 2: The accumulated phase $\Gamma(\tau)$ of the first eigenstate for an
interval.

###  3.2 The geometric phase for a time dependent quantum graph

Continuing just as in the case of an interval and taking the limit
$T\rightarrow\infty$, one obtains the following analogous expression for the
geometric phase:

$\displaystyle\Gamma^{(n)}$ $\displaystyle=$
$\displaystyle\frac{1}{8}\sum_{e=1}^{E}\int_{0}^{1}{\rm d}\tau\frac{{\rm
d}^{2}L_{e}^{2}(\tau)}{{\rm d}{\tau}^{2}}\ \
\left\langle\left(\xi^{2}-\frac{1}{4}\right)_{e}\right\rangle_{\tau}^{(n)}$
$\displaystyle=$ $\displaystyle\frac{1}{4}\sum_{e=1}^{E}\int_{0}^{1}{\rm
d}\tau\left[L_{e}(\tau)\frac{{\rm d}^{2}L_{e}(\tau)}{{\rm
d}{\tau}^{2}}+\left(\frac{{\rm d}L_{e}(\tau)}{{\rm
d}{\tau}}\right)^{2}\right]\
\left\langle\left(\xi^{2}-\frac{1}{4}\right)_{e}\right\rangle_{\tau}^{(n)}$

From the expression above we see that unlike the case of an interval, the
geometric phase may be nontrivial, and its value depends on the relations
between the expectation values
$\left\langle\left(\xi^{2}-\frac{1}{4}\right)_{e}\right\rangle^{(n)}_{\tau}$
which multiply the components of the vector $\frac{{\rm
d}^{2}L_{e}^{2}(\tau)}{{\rm d}{\tau}^{2}}$ in parameter space.

Integrating by parts, we get the alternative form

$\displaystyle\Gamma^{(n)}$ $\displaystyle=$
$\displaystyle-\frac{1}{4}\sum_{e,e^{\prime}=1}^{E}\int_{0}^{1}{\rm
d}\tau\frac{{\rm d}L_{e}(\tau)}{{\rm d}{\tau}}\left\\{\ L_{e}(\tau)\
\frac{\partial\ \ }{\partial
L_{e^{\prime}}}\left\langle\left(\xi^{2}-\frac{1}{4}\right)_{e}\right\rangle_{\tau}^{(n)}\right\\}\frac{{\rm
d}L_{e^{\prime}}(\tau)}{{\rm d}{\tau}}\ $ $\displaystyle=$
$\displaystyle-\pi^{2}\sum_{e,e^{\prime}=1}^{E}\rho_{e}\rho_{e^{\prime}}\nu_{e}\nu_{e^{\prime}}\int_{0}^{1}{\rm
d}\tau\ \sin(2\pi\nu_{e}\tau+\alpha_{e})\ A_{(n);e,e^{\prime}}(\tau)\
\sin(2\pi\nu_{e^{\prime}}\tau+\alpha_{e^{\prime}})$

where,

$A_{e,e^{\prime}}((n);\tau)\ =\ L_{e}(\tau)\ \frac{\partial\ \ }{\partial
L_{e^{\prime}}}\left\langle\left(\xi^{2}-\frac{1}{4}\right)_{e}\right\rangle_{\tau}^{(n)}.$
(49)

In leading order, the last line in (3.2) implies that the geometric phase
depends at least quadratically on the parameters $\rho_{e}$ which determine
the size of the Lissajous figure in parameter space. If we denote the term in
$A_{e,e^{\prime}}$ which is of zero order in these parameters (and therefore
independent of $\tau$) by $A_{e,e^{\prime}}^{(0)}(n)$, we obtain the
following:

######  Theorem 2.

Expanding the integrand in (3.2) with respect to the amplitudes $\rho_{e}$,
the geometric phase is given in leading order by

$\Gamma^{(n)}=-\frac{\pi^{2}}{2}\sum_{e=1}^{E}\rho_{e}^{2}\ \nu_{e}^{2}\
A^{(0)}(n)_{e,e}.$ (50)

This expression depends exclusively on quantities which depend on the graph
with constant lengths $\\{\bar{L}_{e}\\}$. Thus, one can deduce the Lissajous
frequencies $\vec{\nu}$ from the value of the geometric phase $\Gamma^{(n)}$
in the limit $\|\vec{\rho}\|\rightarrow 0$.

To demonstrate the theorem, we perform the explicit computation for the ground
state of some compact, connected metric graph. Note that for our instantaneous
operator, the ground state is given by $k_{0}(\tau)=0$ with a piecewise
constant eigenfunction:

$g_{e}=\frac{L_{e}^{\frac{1}{2}}}{\left(\sum_{j\in\mathcal{E}}L_{j}\right)^{\frac{1}{2}}},\
\ {\rm and}\ \
\left\langle\xi_{e}^{2}-\frac{1}{4}\right\rangle^{(0)}=-\frac{L_{e}}{6\sum_{j\in\mathcal{E}}L_{j}}.$
(51)

The corresponding geometric phase is given by

$\Gamma^{(0)}=-\frac{1}{48}\sum_{e\in\mathcal{E}}\int_{0}^{1}\left(\frac{L_{e}(\tau)}{\sum_{j\in\mathcal{E}}L_{j}(\tau)}\frac{{\rm
d}^{2}}{{d}\tau^{2}}L_{e}(\tau)^{2}\right)\ {\rm d}\tau\ .$

Denote the sum of the mean edge lengths by
$\bar{L}:=\sum_{e\in\mathcal{E}}\bar{L}_{e}$. Then plugging in the expression
for $L_{e}(\tau)$ into (3.2) and expanding with respect to the amplitudes
$\rho_{e}$, we obtain the following:

######  Theorem 3.

The geometric phase of the ground state is given in leading order by

$\Gamma^{(0)}=\frac{\pi^{2}}{12\bar{L}^{2}}\sum_{e\in\mathcal{E}}\nu_{e}^{2}\bar{L}_{e}(\bar{L}-\bar{L}_{e})\rho_{e}^{2}.$
(52)

The phases $\vec{\alpha}$ can be deduced from the next to leading order terms.
The general expression for these terms is quite long and not simple to
compute. In the next section we show that for the ground state of the 3-star
graph, explicit expressions can be obtained.

Figure 3: The 3-star graph.

##  4 Example - A star graph with three edges

In the present section we illustrate the general theory by applying it to a
simple graph - the star graph with three edges. After writing down the
explicit expressions for the spectral parameters and the resulting geometric
phase, we show how perturbation theory can be used to sharpen Theorems 2 and 3
in this particular case. We also provide numerical demonstrations for the
results.

The 3-star graph consists of a central vertex (denoted by $v=0$) and three
edges (denoted by $e=1,2,3$) with lengths $L_{e}(\tau)$ connecting the central
vertex to the external ones (denoted by $v=1,2,3$). The modified Neumann-
Kirchhoff condition (31), (32) is imposed at all vertices. The scaled
coordinates $\xi_{e}$ start at the central vertex with the value
$\xi_{e}=-\frac{1}{2}$ and terminate at the exterior vertices where
$\xi_{e}=+\frac{1}{2}$ .

Since we are interested in the instantaneous representation, we suppress the
$\tau$ dependence of the eigenfunctions when it is irrelevant. The
instantaneous Schrödinger equation at the edges is given by

$-\frac{\partial^{2}g_{e}(\xi_{e})}{\partial\xi_{e}^{2}}=(kL_{e}(\tau))^{2}g_{e}(\xi_{e})\
,$ (53)

with the following vertex conditions:

$\displaystyle\frac{\partial
g_{e}}{\partial\xi_{e}}\left(\frac{1}{2}\right)=0\ {\rm for}\ \
e\in\\{1,2,3\\},\ \ $ (54) $\displaystyle
L_{e}^{-1/2}g_{e}\left(-\frac{1}{2}\right)=L_{e^{\prime}}^{-\frac{1}{2}}g_{e^{\prime}}\left(-\frac{1}{2}\right)\
\ \forall e,e^{\prime},\ \ $ (55)
$\displaystyle\sum_{e=1}^{3}\frac{1}{L_{e}^{3/2}}\frac{\partial
g_{e}}{\partial\xi_{e}}\left(-\frac{1}{2}\right)=0.\ \ $ (56)

The soultions are of the form:

$g_{e}(\xi_{e})=L_{e}^{\frac{1}{2}}a_{e}\cos\left(kL_{e}\left(\xi_{e}-\frac{1}{2}\right)\right),$
(57)

where $k=k_{n}(\tau)$ is some zero of the secular equation:

$\zeta(k;\tau):=\sum_{e=1}^{3}\tan(kL_{e})=0.$ (58)

A typical $\tau$ dependent spectrum of the star graph for two knots is shown
in Figure 4.

Figure 4: The $\tau$ dependent spectrum for the first 60 levels of the star
graph for the unknot (left) and 9/2 knot (right), with the same parameters as
in Figure 1.

For $k^{2}>0$ in the spectrum, the amplitudes $a_{e}$ are given by

$a_{1}=\frac{\cos(kL_{2})\cos(kL_{3})}{N},\ \ \
a_{2}=\frac{\cos(kL_{3})\cos(kL_{1})}{N},\ \ \
a_{3}=\frac{\cos(kL_{1})\cos(kL_{2})}{N}\ ,$ (59)

where the normalization factor $N$ satisfies

$N^{2}=\frac{1}{2}\sum_{e=1}^{3}L_{e}a_{e}^{2}\left(1+\frac{\sin
2kL_{e}}{2kL_{e}}\right).$ (60)

A straightforward computation then gives

$\left\langle\xi_{e}^{2}-\frac{1}{4}\right\rangle^{(n)}_{\tau}=\frac{a_{e}^{2}}{24N^{2}(kL_{e})^{3}}\left(3\sin\left(2kL_{e}\right)-3kL_{e}\cos\left(2kL_{e}\right)-3kL_{e}+2(kL_{e})^{3}\right).$
(61)

The above is a rather complicated expression of $\tau$, since $L_{e}$ itself
is a periodic function of $\tau$, and both $k_{n}$, $a_{e}$ and $N^{2}$ depend
on all three edge lengths. Substituting this in (3.2), one should still
compute the integral over $\tau$ to obtain the geometric phase. This task can
be carried out numerically, where as a preliminary step, the spectrum
$k_{n}(\tau)$ is computed numerically from the secular equation.

Attempts to obtain approximate expressions for the geometric phase rely on the
fact that the time dependence appears always through the terms
$\rho_{e}\cos(2\pi\nu_{e}\tau+\alpha_{e})$. This calls for a perturbative
expansion with respect to the coefficients $\rho_{e}$. This was also the
approach used in Theorems 2 and 3. However, going beyond the leading
approximation is technically hard. The ground state with $k=0$ seems to be the
only case for which the analytic computation of the geometric phase to higher
order is possible.

Using Theorem 3, we know that the geometric phase of the ground state is given
in leading order by (52). Plugging our expressions for the spectral parameters
of the ground state for the 3-star into (52), one may expand it in powers of
the $\rho_{e}$. Doing so, and writing
$\cos(2\pi\nu_{e}\tau+\alpha_{e})=\frac{1}{2}(e^{i(2\pi\nu_{e}\tau+\alpha_{e})}+e^{-i(2\pi\nu_{e}\tau+\alpha_{e})})$
and then expanding each of the powers and integrating, one finds that a
typical term of order $q$ is of the form

$\displaystyle\prod_{e=1}^{3}\rho_{e}^{q_{e}}\int_{0}^{1}{\rm
d}\tau\cos^{q_{e}}(2\pi\nu_{e}\tau+\alpha_{e})$ (62) $\displaystyle=$
$\displaystyle\frac{\prod\rho_{e}^{q_{e}}}{2^{q}}\int_{0}^{1}{\rm
d}\tau\sum_{j_{1},j_{2},j_{3}=0}^{q_{1},q_{2},q_{3}}\prod_{e=1}^{3}\left(\begin{array}[]{c}q_{e}\\\
j_{e}\end{array}\right)\ \delta[(\sum_{e=1}^{3}n_{e}(2j_{e}-q_{e})],$ (65)

where $q=\sum_{e=1}^{3}q_{e}$ (In some terms one of the cosine factors may be
replaced by a sine which does not change the argument) . The terms which do
not cancel after integration are selected by the Kronecker delta, which gives
the demand that there exist triplets of integers $\\{j_{e}\\}_{e=1}^{3}$ with
$0\leq j_{e}\leq q_{e}$ such that the frequency vector of the Lissajous curve
$\vec{\nu}$ satisfies

$2\ \sum_{e=1}^{3}\nu_{e}j_{e}=\sum_{e=1}^{3}\nu_{e}q_{e}\ .$ (66)

Condition (66) on the frequency vector is quite restrictive, where the most
obvious restriction is that $\sum_{e=1}^{3}\nu_{e}q_{e}$ must be even. Another
consequence of (66) is that if the three $\nu_{i}$ are odd, then at least one
of the $q_{i}$ must be even. Thus, in this case, the leading order $q$ cannot
be less than four. However, for $\vec{\nu}=(2,3,5)$, the power $q=3$ is
allowed since (66) is satisfied for $q_{1}=q_{2}=q_{3}=1$ and
$j_{1}=j_{2}=1,j_{3}=0$. This is demonstrated by the blue line in (6). For a
knot with frequency vector $(2,3,7)$ or $(3,5,7)$, the lowest leading order is
$q=4$, which is shown by the green and red lines in (6).

It is instructive to study in detail the knot with $\vec{\nu}=(2,3,5)$, $q=3$.
One can expand the geometric phase in this case to third order and get:

$\displaystyle\Gamma^{(0)}$ $\displaystyle=$ $\displaystyle\
\frac{\pi^{2}}{12\bar{L}^{2}}\sum_{e=1}^{3}\nu_{e}^{2}\
\bar{L}_{e}(\bar{L}-\bar{L}_{e})\rho_{e}^{2}$ $\displaystyle+$
$\displaystyle\sum_{cyc}\
\frac{\pi^{2}\rho_{1}\rho_{2}\rho_{3}\left(\bar{L}_{1}^{2}\nu_{1}^{2}+\bar{L}_{2}^{2}\nu_{2}^{2}+\bar{L}_{3}^{2}\nu_{3}^{2}(\nu_{1}+\nu_{2})^{2}\right)}{12\bar{L}^{3}}\cos(\alpha_{1}+\alpha_{2}-\alpha_{3})\delta_{\nu_{1}+\nu_{2}-\nu_{3}}\
$ $\displaystyle+$ $\displaystyle{\mathcal{O}}(\|\vec{\rho}\|^{4})\ ,$ (68)

where $cyc$ refers to the cyclic summation over all cyclic permutations of the
indices $1,2,3$.

Figure 5: Ground state geometric phases for Lissajous knot of frequencies
(2,3,5), (2,3,7),(3,5,7). Upper line: $\Gamma_{0}/\rho^{2}$ as a function of
$Log[\rho]$ . Lower line: $\Gamma_{0}/\rho^{2}$ as a function of the phase
shifts difference $\phi=\alpha_{3}-\alpha_{2},\ ({\rm with}\ \alpha_{1}=0)$ .

The computation above demonstrates the fact that the phases $\vec{\alpha}$ can
be determined by the leading order for certain combinations of $\vec{\nu}$,
and thus recover the knot.

For the more general case, one can continue to write even higher order terms
which will depend on the phases, but the resulting expressions are lengthy and
not very instructive. The sensitivity of the geometric phase to the phase
shifts $\vec{\alpha}$ is demonstrated in the lower part of Figure 5

Figure 6: Ground state data for a Lissajous knot of frequencies
(2,3,5)-(blue), (2,3,7)-(green),(3,5,7)-(red), showing
$Log[\delta\Gamma_{0}/\rho^{2}]$ as a function of $Log[\rho]$

The data presented in Figures (5, 6) results from numerical integration of
(3.2) for Lissajoux knots with frequency vectors
$\vec{\nu}=(2,3,5),(2,3,7),(3,5,7)$, phase
$\vec{\alpha}=(0,\frac{\pi}{3},\frac{\pi}{5})$, and identical amplitudes
$\rho_{e}:=\rho$. The upper part of Figure (5) shows
$\Gamma_{0}(\rho)/\rho^{2}$ as a function of $\log(\rho)$. The continuous dark
line is the leading expression for $\Gamma_{0}/\rho^{2}$ as derived in (4).
Figure (6) shows the difference ($\delta\Gamma_{0}$) between
$\Gamma_{0}/\rho^{2}$ and its leading order as a function of $\rho$ on a
double Log plot. The expected power laws are demonstrated.

Figure (7) shows the geometric phase for the first $60$ levels of the star
graph for the two different knots for which the instantaneous spectra are
shown in Figure (4). We cannot offer an analytical expression for the
numerical results. However, we can heuristically explain their main features.
The values of $\Gamma(n)$ are shown in three colors which distinguish between
$n$ such that red stands for $n=_{{\rm mod\ 3}}1$, green for $n=_{{\rm mod\
3}}2$ and blue for $n=_{{\rm mod\ 3}}3$ This is done in order to draw the
attention to the incipient braided nature of $\Gamma(n)$ which looks like the
union of three oscillatory functions. This behavior is most clearly observed
for lower values of $\rho$, where the spectral curves $k_{n}(\tau)$ are
bunched into triplets where the two lower levels are closer to one another and
are distant from the third. This is typical for the 3-star graph. The
oscillatory behavior follows from the fact that the time dependence appears as
a trigonometric function of a trigonometric function of $\tau$. This suggest
towards a possible Bessel-like oscillation property of the geometric phase as
a function of the energy level, when one considers the integral representation
of the Bessel functions (Hansen-Bessel formula). While the formula for the
geometric phase is not identical to that of the Bessel functions, it seems
that at least qualitively, similar oscillations do occur.

The second feature is that the three branches of the phase seem to converge to
a constant. This constant is increasing with respect to $\rho$.

Figure 7: The geometric phase for the first 60 levels of the star graph for
the unknot (left) and 9/2 knot (right), with the same parameters as in Figure
1.

##  5 Further discussion and summary

The present work lays the ground for the study of evolution on time dependent
quantum graphs and for the computation of the associated geometric phase. It
turns out that for most vertex conditions (including the standard Neumann-
Kirchhoff condition, but excluding Dirichlet), the movement of the graph edges
generates a probability flux through the boundary of the graph, which means
that the associated instantaneous operator is not self–adjoint. Theorem 1
treats this problem by providing the modified vertex conditions on the graph
which render the time evolution unitary, and allow us to consider the
geometric phase presented in Formula (3.2) in the adiabatic limit.

It is worthwhile to note that the resulting Schrödinger operator (23) depends
on the changing edge lengths $L_{e}(t/T)$, which results in a different wave
velocity on each edge. This is analogous to the Schrödinger operator which
corresponds to a quantum graph made of a composite materials, where each edge
has different propogation speed (see [11]). Thus, the time dependent quantum
graph may, in a sense, be treated as a network of composite conductors whose
wave velocity depends on the dynamics of the graph evolution. It is
interesting to think whether in such systems, a similar geometric phase may
arise.

Theorems 2 and 3 show that the leading order behavior of the geometric phase
can be used to deduce the frequencies which define the Lissajous knot. To
completely determine the knot type, one has to also know the higher order
behavior, in order to extract the phase shifts of the knot. While the example
in Section 4 shows that this can be done for the 3–star graph, the computation
for a general graph is difficult to perform, even with perturbation theory.
This difficulty stems from the fact that the geometric phase is given by an
integral of a periodic function, whose argument is itself a periodic function
(the edge lengths).

Formula (52) presented in Theorem 3 suggests that in leading order, the
geometric phase of the ground state depends only on the given knot in
parameter space, and in fact, not on the graph itself (only its edge lengths).
For instance, up to leading order, the ground state of the 3–star graph should
have the same geometric phase as that of a triangle graph, although the two
graph have different geometric structures. For instance, the 3–star graph is
simply connected while the triangle graph is not. It would be interesting to
understand if when moving to higher orders, the geometric phase can be used to
distinguish between more subtle geometric features of the two different graphs
(like their Betti numbers). In other words - if one can ‘hear’ additional
geometric features of the graph from the geometric phase, rather than just the
edge lengths.

## Acknowledgements

We are obliged to Professor Michael Berry for explaining the physical origin
of the magnetic phase. Thanks are due to Professor Kirill Cherednichenko for
his continuous interest, and for pointing out the similarity between the
boundary conditions in the present problem and the corresponding conditions
which hold for graphs with variable properties which appear in the study of
composite conductors. We also thank Lior Alon for some interesting discussions
about the behavior of the phase near degenerate states. One of us (US) thanks
the University of Bath for the hospitality while holding a David Parkin
Professorship (September 2021- March 2022).

## References

## References

  * [1] S. W. Doescher, and M. H. Rice Infinite Square-Well Potential with a Moving Wall, American Journal of Physics 37  , 1246 (1969)
  * [2] A. Munier, J. R. Burgan, M. Feix, E. Fijalkow Schrödinger Equation with Time-Dependent Boundary Conditions J. Math. Phys. 22 6, (1981)
  * [3] A. Duca, R. Joly, Schrödinger Equation in Moving Domains Ann. Henri Poincare 22 , 2029–2063 (2021)
  * [4] Gerald T. Moore, Quantum Theory of the Electromagnetic Field in a Variable-Length One-Dimensional Cavity J. Math. Phys. , 11, 2679 (1970).
  * [5] D. U. Matrasulov, J. R. Yusupov, K. K. Sabirov,Z. A. Sobirov Time-dependent Quantum Graph, Nanosystems: Physics, Chemistry, Mathematics, 2015,6 (2), P. 173–181 .
  * [6] M. V. Berry Quantal Phase Factors Accompanying Adiabatic Changes, Proc. R. Soc. Lond. A 392, 45-57 (1984)
  * [7] A. Boocher, J. Daigle, J. Hoste, W. Zheng, Sampling Lissajous and Fourier Knots, Experimental Mathematics 18 , 481 (2009).
  * [8] M. G. V. Bogle, J. E. Hearst, V. F. R. Jones, L. Stoilov Lissajous knots, Jour. Knot Theory and its Ramifications, 3 (1994) 121
  * [9] G. Berkolaiko, P. Kuchment, Introduction to Quantum Graphs (2013)
  * [10] T. Kottos, U. Smilansky, Chaotic Scattering on Graphs (2000)
  * [11] Kirill Cherednichenko, Private Correspondence (2022)
  * [12] M.V. Berry, Private Correspondence (2022)


# Passivation and characterization of charge defects in ambipolar silicon
quantum dots

Paul C. Spruijtenburg  NanoElectronics Group, MESA+ Institute for
Nanotechnology, University of Twente, P.O. Box 217, 7500 AE Enschede The
Netherlands  Sergey V. Amitonov  NanoElectronics Group, MESA+ Institute for
Nanotechnology, University of Twente, P.O. Box 217, 7500 AE Enschede The
Netherlands  Filipp Mueller  NanoElectronics Group, MESA+ Institute for
Nanotechnology, University of Twente, P.O. Box 217, 7500 AE Enschede The
Netherlands  Wilfred G. van der Wiel  NanoElectronics Group, MESA+ Institute
for Nanotechnology, University of Twente, P.O. Box 217, 7500 AE Enschede The
Netherlands  Floris A. Zwanenburg  NanoElectronics Group, MESA+ Institute for
Nanotechnology, University of Twente, P.O. Box 217, 7500 AE Enschede The
Netherlands<EMAIL_ADDRESS>
###### Abstract

In this Report we show the role of charge defects in the context of the
formation of electrostatically defined quantum dots. We introduce a barrier
array structure to probe defects at multiple locations in a single device. We
measure samples both before and after an annealing process which uses an
$\textrm{Al}_{2}\textrm{O}_{3}$ overlayer, grown by atomic layer deposition.
After passivation of the majority of charge defects with annealing we can
electrostatically define hole quantum dots up to 180 nm in length. Our
ambipolar structures reveal amphoteric charge defects that remain after
annealing with charging energies of $\simeq$10 meV in both the positive and
negative charge state.

###### keywords:

atomic layer deposition; Pb centers; dangling bonds; quantum dots; holes;
charge defects

In order to perform sufficient coherent operations in a proposed quantum
computer [1], the quantum states of the corresponding qubits are required to
be long-lived. In the scheme proposed by Loss and DiVincenzo[2], quantum logic
gates perform operations on coupled spin states of single electrons in
neighboring quantum dots (QDs). Most experiments have focused on quantum dots
formed in III-V semiconductors, especially GaAs [3, 4]. However, electron spin
coherence in those materials is limited by hyperfine interactions with nuclear
spins and spin-orbit coupling. Group IV materials have been shown to have long
spin lifetimes because of weak spin-orbit interactions and the predominance of
spin-zero nuclei. This prospect has stimulated significant experimental effort
to isolate single charges in carbon nanotubes [5, 6], Si/SiGe heterostructures
[7, 8], Si nanowires [9], planar Si MOS structures [10], and dopants in Si
[11, 12, 13, 14, 15]. Silicon not only holds promise for very long coherence
times [16, 17, 18], but also for bringing scalability of quantum devices one
step closer, and has thus attracted much attention for quantum computing
purposes [19, 20]. For the quantum states to be well-defined, it is necessary
to define singular quantum dots. Disorder in the confinement potential can
lead to unintentional quantum dots[21, 22], resonant tunneling features[23,
24], and Anderson localization[25]. In silicon planar quantum dots, nanometer-
size Coulomb islands are electrostatically defined in silicon in a gated
MOSFET-type heterostructure. Apart from the quality of the crystalline silicon
the quality and disorder in the heterostructure is directly linked to the
disorder of the confinement potential.

In this Report, we focus on disorder caused by defects at the
Si/$\textrm{Si}\textrm{O}_{2}$ interface and investigate their effect on the
formation of electrostatically defined quantum dots. We do this by creating an
array of barrier electrodes, each capable of locally controlling the potential
for quantum dot formation and probing the respective crystalline environment.
We show that an annealing process with $\textrm{Al}_{2}\textrm{O}_{3}$, grown
by atomic layer deposition (ALD), reduces the disorder such that quantum dots
180 nm in length can be formed. Furthermore, we observe that after annealing a
small percentage of barriers show charge transitions with an amphoteric
character.

The role of defects at the Si/$\textrm{Si}\textrm{O}_{2}$ interface in this
heterostructure has been extensively studied in the context of MOSFET
technology. One of the most studied defects at this interface is the
paramagnetic $\textrm{P}_{\textrm{b}}$ center[26, 27, 28] which is linked to
surface charges, decreased mobility, and the negative-bias temperature
instability (NBTI) effect[29, 30]. For quantum computation,
$\textrm{P}_{\textrm{b}}$ centers have been used in coherent manipulation of
spin-dependent charge-carrier recombination of phosphorus donors.[31] Defects
at the interface have also been associated with random telegraph signals at
low temperatures and characterized by means of single-electron spin resonance
[32, 33] and magnetic-field dependent measurements[34].

The $\textrm{P}_{\textrm{b}}$ center is in essence a Si atom with a dangling
bond that has formed due to the incommensurability of the crystalline Si and
the amorphous $\textrm{Si}\textrm{O}_{2}$. Two types of
$\textrm{P}_{\textrm{b}}$ center can be further distinguished. The
$\textrm{P}_{\textrm{b0}}$ center is backbonded to two Si atoms and one oxygen
atom, while the $\textrm{P}_{\textrm{b1}}$ center is backbonded to three Si
atoms. For the $\textrm{P}_{\textrm{b0}}$ center, the density of interface
states D${}_{\text{it}}$ in the gap is amphoteric in nature with its maxima at
$0.25$ and $0.85$ eV above the valence band [35]. Passivating
$\textrm{P}_{\textrm{b}}$ centers is commonly achieved by introducing hydrogen
to the dangling bonds in an annealing process. This process reduces the
dangling bond spatial density from order $10^{13}\leavevmode\nobreak\
\text{cm}^{-2}$ to $10^{10}\leavevmode\nobreak\ \text{cm}^{-2}$ and improves
numerous parameters related to MOSFET-like transistor operation and has also
been shown to improve parameters related to quantum dot formation [37, 38].

Addressing only a single quantum dot defined by electrostatic gating can be
hampered by a charge defect. If a charge defect is located beneath an
electrode intended to form a tunnel barrier, this electrode acts a gate on
this defect. Instead of forming a single tunnel barrier to the QD, the levels
of the defect remain available for transport resulting in additional single-
charge tunneling events. This situation is sketched for an unpassivated charge
defect in Figure 1a for the right tunnel barrier.

While the valence band of the right tunnel barrier is below the Fermi level,
the level of the charge defect is above the Fermi level and is thus available
for hole transport. If the defect is in close enough proximity to the QD, an
ill-defined double quantum dot can form. If the defect is passivated by
hydrogen, the potential is smooth and behaves as a normal tunnel barrier for
the QD. The same arguments hold for an impurity level just below the
conduction band when measuring an electron QD.

Figure 1b shows an atomic force microscope image and Figure 1c a schematic
cross section of our device structure, made with a combination of optical and
electron-beam lithography (EBL), based on the recipe as described by Angus _et
al._[38]. We further distinguish two device types — A and B — specified below.
Near-intrinsic (100) silicon ($\rho\geq 10\leavevmode\nobreak\
\text{k}\Omega\cdot\text{cm}$) is used as the substrate. Source and drain
regions are implanted with boron/phosphorus dopant atoms, which are activated
by rapid thermal annealing, and serve as hole/electron reservoirs. Ohmic
contacts to these regions are made by sputtering Al-Si alloy (99:1) contact
pads. A 10 nm thick high-quality $\textrm{Si}\textrm{O}_{2}$ oxide window is
thermally grown at $900^{\circ}\text{C}$ and serves as an insulating barrier
between the substrate and the aluminum gates. Contact pads for gates are
defined using optical lithography followed by development, evaporation of
Ti/Pd, and subsequent lift-off. At this point, for device type B, a $5$ nm
$\textrm{Al}_{2}\textrm{O}_{3}$ “underlayer” using ALD at a process
temperature of $250^{\circ}\text{C}$ is grown on top of the
$\textrm{Si}\textrm{O}_{2}$. Subsequently, EBL followed by evaporation and
liftoff is used to define the sub-micron aluminum gates. These gates will
electrostatically control charge accumulation. In the first step, the barrier
gates are created and thermally oxidized. In the second step, the lead gate is
created. Atomic force microscopy images show barrier gates with a typical
width of 35 nm, and a minimum separation of 180 nm, while the lead gate has a
width of approximately 80 nm, see Figure 1b. Following the creation of the
aluminum gate structures, we perform a first set of measurements. Then we grow
a final “overlayer” of 5 nm $\textrm{Al}_{2}\textrm{O}_{3}$ at
$100^{\circ}\text{C}$, and the sample is heated to $300^{\circ}\text{C}$ for
45 minutes in pure Ar ambient of 10 mbar. As a final step, the sample is
cleaned in UV ozone for 2 minutes to remove any moisture present.

(a)

(b)

(c)

Figure 1: Si barrier array device structure.
LABEL:sub@fig:schematic_potentialprofile Electrochemical potential profile of
a hole quantum dot with and without a passivated and unpassivated charge
defect in the left and right tunnel barrier respectively. In case of an
unpassivated charge defect (“CD”), an extra level in the bandgap (green) is
available for transport. LABEL:sub@fig:AFM Atomic force microscopy image of
the device, showing the lead gate L horizontally across the image. The
barriers B1 through B10 allow local probing of the 2DHG/2DEG at many
locations. LABEL:sub@fig:schematic_layout A schematic cross-sectional image of
the device. The hydrogen (H) present in the $\textrm{Al}_{2}\textrm{O}_{3}$
can passivate charge defects during annealing. The applied voltage on the
aluminum gates creates a 2DHG/2DEG, indicated by the dashed lines.

To study both defects and their influence on QD formation, we use a twofold
approach: by using an array of electrodes (seen in an AFM image in Figure 1b),
each acting as a single barrier, we can locally probe the 2D carrier gas and
observe single-charge tunneling events. At the same time, using an ambipolar
device structure[39, 40] it is possible to infer if accessible levels are
located near the conduction or valence band edge by either inducing a 2D hole
gas (2DHG) or 2D electron gas (2DEG). Thus, our ambipolar structures allow
investigation of both the hole and electron nature of these defects.

Recent work has shown that $\textrm{Al}_{2}\textrm{O}_{3}$ grown with ALD
contains hydrogen and can be used for passivation[41]. We use this property of
ALD-grown $\textrm{Al}_{2}\textrm{O}_{3}$ in a separate anneal step after
$\textrm{Al}_{2}\textrm{O}_{3}$ deposition. We further postulate that the
mechanical stability of our gate structure with pure Al is improved by
encapsulation in this $\textrm{Al}_{2}\textrm{O}_{3}$ oxide. It is
energetically more favorable to form $\textrm{Al}_{2}\textrm{O}_{3}$ than
$\textrm{Si}\textrm{O}_{2}$ because the Gibbs free energy change of
$\textrm{Al}_{2}\textrm{O}_{3}$ formation is lower than that of
$\textrm{Si}\textrm{O}_{2}$. Al is known to reduce $\textrm{Si}\textrm{O}_{2}$
as confirmed by transmission electron microscopy of our samples[42]. It has
also been shown that Al can spike through $\textrm{Si}\textrm{O}_{2}$ causing
leakages and shorts to the channel[43]. Diffusivity of Al at temperatures
around $300^{\circ}\text{C}$ can cause dewetting[44], or void formation can
occur due to different thermal expansion coefficients[45].

We have found that the overlayer of 5 nm of $\textrm{Al}_{2}\textrm{O}_{3}$
alone in device type A stabilizes the pure aluminum electrodes. Without this
overlayer void formation and dewetting-like behaviour was observed in all
devices after the $300^{\circ}\text{C}$ annealing step (see Supplementary
Information). The presence of an $\textrm{Al}_{2}\textrm{O}_{3}$ underlayer in
device type B, between the $\textrm{Si}\textrm{O}_{2}$ and the Al gate
structure could prevent the reduction reaction that would otherwise occur at
the Al/$\textrm{Si}\textrm{O}_{2}$ interface. The under- and overlayer of
$\textrm{Al}_{2}\textrm{O}_{3}$ respectively change the surface energy and
suppress surface diffusivity of the Al film.

(a)

(b)

(c)

(d)

Figure 2: Transport measurements of device A before and after annealing.
LABEL:sub@fig:traces Current-voltage relation at $T\simeq
4.2\leavevmode\nobreak\ \mathrm{K}$ before (red) and after (green) annealing.
Individual barriers are swept to 0 V while all other gate voltages are
maintained above the threshold voltage with
$V_{\mathrm{SD}}=1\leavevmode\nobreak\ \text{mV}$ and
$V_{\mathrm{L}}=-4.6\leavevmode\nobreak\ \text{V}$ before, and
$V_{\mathrm{L}}=-2.4\leavevmode\nobreak\ \text{V}$ after anneal.
LABEL:sub@fig:linear Linear transport measurements of current $I$ as a
function of barrier gate voltages $V_{\mathrm{B10}}$ and $V_{\mathrm{B9}}$
before (left) and after LABEL:sub@fig:linear (right) annealing.
LABEL:sub@fig:bias_pre The differential conductance $dI/dV_{\mathrm{SD}}$, as
a function of source-drain bias $V_{SD}$ and barrier gate voltage $V_{B10}$
before annealing, taken at $V_{\mathrm{B9}}=-950\leavevmode\nobreak\
\text{mV}$. Irregular Coulomb-diamond like features are visible with charging
energies up to 10 meV. LABEL:sub@fig:bias_post The same measurement after
annealing with $V_{\mathrm{B10}}=-350\leavevmode\nobreak\ \text{mV}$ and
$T\simeq 10\leavevmode\nobreak\ \text{mK}$. Well-defined Coulomb diamonds with
charging energies of 2-3 meV are visible, corresponding to an
electrostatically defined Coulomb island 180 nm in length and 80 nm in width.

We have measured various devices and report here characteristic behaviour of
two types of samples A and B. First, we will focus on device A and the
influence of charge defects on the formation of electrostatically defined QDs.
Device A (shown in Figures 1b and 1c) does not contain a layer of
$\textrm{Al}_{2}\textrm{O}_{3}$ between the $\textrm{Si}\textrm{O}_{2}$ and
the Al gate structures. For its transport characteristics we will focus on
hole transport. Device A is measured once before
$\textrm{Al}_{2}\textrm{O}_{3}$ overlayer growth and subsequent annealing
step. After this passivation procedure the device is measured again. The
current $I$ is measured at the drain while the drain is kept at ground.
Additionally, all voltages are given with respect to ground.

Figure 2a shows the current-voltage relation of the individual barriers before
and after annealing. For these measurements all electrodes except the barrier
electrode in question are at sufficiently high voltage to allow conduction,
while one of the barrier electrodes is swept to such a voltage that the
barrier will become opaque and thus “cut off” the channel locally. Cut-off
voltages decrease from a spread of -2 V to -1 V to a more narrow region of -1
V to -0.5 V. This is expected to be due to the reduction of a charge layer
formed by charge defects present before annealing. The decreased number of
resonances in the current $I$ indicates disorder in the device has decreased
substantially after annealing. The threshold voltage $V_{\mathrm{Th}}$, here
defined as the voltage at which the entire channel becomes conducting when
applying the same voltage to all gates, is also lowered after annealing from
$V_{\mathrm{Th}}\simeq-4$ V to $V_{\mathrm{Th}}\simeq-2.3$ V (not shown).

We have observed that nearly all resonances associated with the barriers
disappear after $\textrm{Al}_{2}\textrm{O}_{3}$ growth and the subsequent
annealing procedure. We have confirmed this in separate bias spectroscopy
measurements where Coulomb diamonds of order $E_{C}\simeq
10\leavevmode\nobreak\ \text{meV}$ coupled exclusively to one barrier gate are
no longer present after annealing.

Before annealing, Figure 2b (left) contains at least two sets of parallel
lines of Coulombs peaks with distinct coupling to B9 and B10. One set of
nearly horizontal lines runs from -2.4 V to -2.3 V parallel to the $V_{B9}$
axis. We ascribe these horizontal lines to a charge defect underneath B10[21].
Another feature is seen in the increased tunnel current running vertically
around $V_{\mathrm{B9}}=-900\leavevmode\nobreak\ \text{mV}$, indicating
another feature coupled strongly to B9 with different periodicity from the one
underneath B10. A set of irregular diagonal lines is also visible. Thus, this
feature is coupled to both B9 and B10, signifying a quantum dot centered
between B9 and B10. However, these lines are interrupted and not continuous
when intersecting the features associated with the defects, which implies a QD
in series with at least 2 defects. After annealing the same device shows very
ordered linear transport data in Figure 2b (right) ($T\simeq
10\leavevmode\nobreak\ \mathrm{mK}$) at
$V_{\mathrm{SD}}=0.3\leavevmode\nobreak\ \text{mV}$ with uninterrupted
diagonal lines equally coupled to B9, B10. We extract
$\sfrac{C_{\mathrm{B10}}}{C_{\mathrm{B9}}}\simeq 1$, where
$C_{\textrm{B}_{i}}$ is the capacitance of barrier gate i to the dot. The
absence of a charge sensor does not allow us to probe how many electrons are
on the dot but the amount of visible charge transitions ($\simeq$ 50) and the
lack of interruptions in the charge transition lines indicates the dot has a
very controllable charge occupation over a large range of voltages.

We next perform bias spectroscopy in order to more accurately determine the
ability to form quantum dots. A bias spectroscopy before annealing reveals
several overlapping sets of Coulomb diamonds, see Figure 2c. This is
indicative of charge centers in series with the QD. Coulomb diamonds with a
charging energy of $E_{C}\simeq 10\leavevmode\nobreak\ \text{meV}$ are visible
as well as another set of diamonds with a charging energy 3-5 meV. After
annealing, Figure 2d shows very regular Coulomb diamonds with a well-defined
charging energy of $E_{C}\simeq 2\leavevmode\nobreak\ \text{meV}$. Below
$V_{\mathrm{B9}}=400\leavevmode\nobreak\ \text{meV}$ the diamonds are closing.
Around $V_{\mathrm{B9}}=-400\leavevmode\nobreak\ \text{mV}$ a band of Coulomb
blockade, with $|V_{\mathrm{SD}}|<0.3\leavevmode\nobreak\ \text{mV}$, causes
the diamonds not to close. This additional blockade is not due to the dot
itself but rather from blockade present further along the 3 micron long
channel formed by the lead gate which we were unable to mitigate. The many
charge transitions and regular spacing of the diamonds indicates the low level
of disorder in the device.

(a)

(b)

Figure 3: Intentional quantum dots at both ends of device A. Linear transport
measurements of current $I$ as a function of adjacent barrier gate voltages.
An intentional quantum dot created LABEL:sub@fig:fig3_left between B2 and B3
with $V_{\mathrm{L}}=-2.9\leavevmode\nobreak\ \text{V}$,
$V_{\mathrm{SD}}=0.5\leavevmode\nobreak\ \text{mV}$ and
LABEL:sub@fig:fig3_right between B8 and B9 with
$V_{\mathrm{L}}=-3.1\leavevmode\nobreak\ \text{V}$,
$V_{\mathrm{SD}}=1\leavevmode\nobreak\ \text{mV}$. Diagonal current peaks
indicate charge transitions coupled almost equally to both B2, B3
($\sfrac{C_{\mathrm{B2}}}{C_{\mathrm{B3}}}\simeq 0.75$) and B8, B9
($\sfrac{C_{\mathrm{B8}}}{C_{\mathrm{B9}}}\simeq 0.95$).

Figure 3 shows additional linear transport measurements taken after annealing
at two different ends of device A. Quantum dots QD2 and QD3, indicated in
Figure 1 are about 2 micron apart. Both Figure 3a and Figure 3b have the same
regularly spaced charge transitions as in Figure 2b (right). The charge
transitions have equal capacitive coupling to the barrier gates defining the
quantum dot and there is no interference from unwanted charge defects. We
conclude from Figures 2 and 3 that our annealing process passivates the
majority of defects, and allows formation of hole quantum dots of 180 nm long.

(a)

(b)

Figure 4: Amphoteric nature of the charge center in device
B.LABEL:sub@fig:amphoteric_hole the hole operation regime
($V_{\mathrm{L}}=-1.6\leavevmode\nobreak\ \text{V}$), and the
LABEL:sub@fig:amphoteric_electron electron operation regime
($V_{\mathrm{L}}=1.7\leavevmode\nobreak\ \text{V}$) of the same device. Both
the hole and electron operation regime show clear charge transitions, with
regions of charge stability indicated with h and e for hole- and electron-
occupancy, respectively. The charging energy of the first charge state of both
electron and hole are found to be $E_{C}\simeq 10\leavevmode\nobreak\
\text{meV}$.

Now, we will investigate the nature of a single electrically active defect,
left after the annealing procedure. Device B, in contrast to device A,
contains a layer of $\textrm{Al}_{2}\textrm{O}_{3}$ between the
$\textrm{Si}\textrm{O}_{2}$ and the gate electrodes. We have observed that
this further increases the yield of our devices. The specific effects of this
interlayer on transport characteristics before and after annealing, on the
same device, were not studied. After annealing we observe, in all devices with
this interlayer, comparable threshold voltages and the same low probability of
finding charge traps underneath individual barriers.

In our ambipolar device structure, the same heterostructure and its associated
imperfections can be probed in both the hole and electron operation regime.
Figure 4 shows a bias spectroscopy of a defect still present after annealing
underneath a single barrier. The Coulomb diamonds near the valence band and
the conduction band are very similar, and several charge states can be clearly
observed. After the last charge transition ($1\textrm{h}\rightarrow
0\textrm{h}$ and $1\textrm{e}\rightarrow 0\textrm{e}$) the diamonds opens up
and no charge transitions are observed up to a very high bias, for both the
hole and electron operation regime. This indicates that the dot has been
emptied completely, although without additional verification from e.g. a
charge sensing experiment this cannot be said with absolute certainty. The
charging energy on both the hole and electron side is found to be $E_{C}\simeq
10\leavevmode\nobreak\ \text{meV}$.

We will now discuss 4 different explanations of these remaining charge
transitions, namely dopant atoms, $\textrm{E}^{\prime}$ centers, unintentional
quantum dots, and $\textrm{P}_{\textrm{b}}$ centers.

The chance that we are dealing with a dopant atom is small, since our
background doping is $\leq 10^{12}\textrm{cm}^{-3}$. Additionally, these would
not show amphoteric behaviour and charging energies below 20 meV.

$\textrm{E}^{\prime}$ centers, another type of $\textrm{Si}\textrm{O}_{2}$
defect, are described as non-amphoteric deep hole traps and responsible for
fixed positive charge.[48] Primarily induced by radiation,
$\textrm{E}^{\prime}$ centers are often characterized as border traps[49],
located further from the conduction channel. Whether or not this defect could
cause Coulomb blockade would depend on its distance from conduction channel.
Radiation induced damage caused by electron beam evaporation[50], used in our
sample fabrication, could lead to higher incidence of this type of defect.

One other explanation for the amphoteric nature of the observed charge
transition is the formation of a so-called unintentional quantum dot, created
by e.g. disorder or strain[47]. In case of the latter the conduction or
valence band is modulated by strain, induced by the difference in thermal
expansion coefficients between the silicon and the aluminum gates in the
heterostructure. This can create tunnel barriers on each side of the barrier
gate, resulting in an unintentional quantum dot. The fact that we observe
quantum dot behaviour in a small percentage of barrier gates could be
explained by variations in grain size of our constituent Al, resulting in non-
uniform strain and tunnel barriers which vary from barrier to barrier.

A system very comparable to the $\textrm{P}_{\textrm{b1}}$ center is the
single dangling bond on hydrogen-terminated intrinsic silicon, i.e. without an
oxide. Density functional theory (DFT) calculations of this single dangling
bond ($\textrm{D}^{\textrm{0}}$) have shown that the associated bound state
has an energy level approximately 13 meV above the valence band maximum[46].
This agrees very well with our extracted value for the charging energy.

The multiple charge states visible in both the hole and electron operation
regime preclude us from identifying the simple case of a single
$\textrm{P}_{\textrm{b}}$ center with only three charge states:
$\textrm{D}^{-}$, $\textrm{D}^{0}$, $\textrm{D}^{+}$. The origin of the
Coulomb diamonds and associated charge states can be explained by the non-
singular nature of the defect: Two $\textrm{P}_{\textrm{b}}$ centers in
parallel would also give rise to overlapping Coulomb diamonds, while the
different addition energies of the neighbouring Coulomb diamonds would not be
expected. It is also known that very proximate dangling bonds can form
molecular type bonding, which could account for the observation of more than
three charge states[46].

The ability to resolve single $\textrm{P}_{\textrm{b}}$ centers could be
attributed to the equilibrium density and the corresponding possibility of
finding a small but definite number of unpassivated $\textrm{P}_{\textrm{b}}$
centers underneath a barrier. This corresponds to our findings that a finite
number of barriers still shows an electrically active defect. To confirm this
type of model, the device structure proposed in this article combined with
statistical analysis on the incidence of charge defects with varying annealing
processes would be needed. Electrical detection of spin resonance (EDSR) [33]
would further provide strong experimental evidence to confirm if these
remaining charge defects are caused by unpassivated $\textrm{P}_{\textrm{b}}$
centers.

In conclusion, we have demonstrated that an annealing process using ALD-grown
$\textrm{Al}_{2}\textrm{O}_{3}$ is able to passivate the majority of the
electrically active defects present at the Si/$\textrm{Si}\textrm{O}_{2}$
interface. We attribute this to the passivating properties of hydrogen present
in the $\textrm{Al}_{2}\textrm{O}_{3}$ layer. After annealing the ability to
form electrostatically defined quantum dots markedly improves with dot lengths
of at least 180 nm. These quantum dots show many charge transitions,
indicating the low level of disorder in the devices. We believe this is most
likely due to the elimination of the majority of $\textrm{P}_{\textrm{b}}$
centers. The ambipolar device architecture then allows us to reveal the
amphoteric behaviour of a remaining charge defect underneath a single barrier
with energy levels $\simeq$10 meV above the valence band and below the
conduction band.

The ability to form singular quantum dots in this heterostructure is strongly
dependent on the density of charge defects at the interface. The control of
the spatial density of charge defects in Si planar quantum dots is thus
necessary in order to succesfully and reproducibly scale up to many devices in
a proposed quantum computer. For $\textrm{P}_{\textrm{b}}$ centers it is known
that the passivation/depassivation reaction leads to an equilibrium, which
implies there will always be a finite defect density at the
Si/$\textrm{Si}\textrm{O}_{2}$ interface. In order to eliminate defects at
this interface, another crystalline material to supplant the amorphous
$\textrm{Si}\textrm{O}_{2}$ would be ideal. Candidates for this could either
be crystalline dielectrics or other crystalline materials such as SiGe[8].

The architecture presented here successfully avoids dewetting of Al in the
annealing step by use of $\textrm{Al}_{2}\textrm{O}_{3}$, which simultaneously
provides hydrogen for the passivation reaction. The ability of the ambipolar
architecture to investigate the amphoteric nature of defects could be a very
useful tool in the arsenal to further investigate charge defects in all manner
of materials.

## References

  * [1] Ladd, T. D. _et al._ Quantum computers.  _Nature_ 464, 45–53.  DOI: 10.1038/nature08812 (2010).
  * [2] Loss, D. & DiVincenzo, D. P.  Quantum computation with quantum dots.  _Physical Review A_ 57, 120–126.  DOI: 10.1103/PhysRevA.57.120 (1998).
  * [3] van der Wiel, W. G. _et al._ Electron transport through double quantum dots.  _Rev. Mod. Phys._ 75, 1–22.  DOI: 10.1103/RevModPhys.75.1 (2002).
  * [4] Hanson, R., Kouwenhoven, L. P., Petta, J. R., Tarucha, S. & Vandersypen, L. M. K.  Spins in few-electron quantum dots.  _Reviews of Modern Physics_ 79, 1217–1265.  DOI: 10.1103/RevModPhys.79.1217 (2007).
  * [5] Jarillo-Herrero, P., Sapmaz, S., Dekker, C., Kouwenhoven, L. P. & Van Der Zant, H. S. J.  Electron-hole symmetry in a semiconducting carbon nanotube quantum dot.  _Nature_ 429, 389–392.  DOI: 10.1038/nature02568 (2004).
  * [6] Higginbotham, A. P., Kuemmeth, F., Hanson, M. P., Gossard, A. C. & Marcus, C. M.  Coherent Operations and Screening in Multielectron Spin Qubits.  _Physical Review Letters_ 112, 026801.  DOI: 10.1103/PhysRevLett.112.026801 (2014).
  * [7] Simmons, C. B. _et al._ Charge Sensing and Controllable Tunnel Coupling in a Si/SiGe Double Quantum Dot.  _Nano Letters_ 9, 3234–3238.  DOI: 10.1021/nl9014974 (2009).
  * [8] Borselli, M. G. _et al._ Undoped accumulation-mode Si/SiGe quantum dots.  _Nanotechnology_ 26, 375202.  DOI: 10.1088/0957-4484/26/37/375202 (2015).
  * [9] Zwanenburg, F. A., van Rijmenam, C. E. W. M., Fang, Y., Lieber, C. M. & Kouwenhoven, L. P.  Spin States of the First Four Holes in a Silicon Nanowire Quantum Dot.  _Nano Letters_ 9, 1071–1079.  DOI: 10.1021/nl803440s (2009).
  * [10] Lim, W. H. _et al._ Electrostatically defined few-electron double quantum dot in silicon.  _Applied Physics Letters_ 94, 173502.  DOI: 10.1063/1.3124242 (2009).
  * [11] Lansbergen, G. P. _et al._ Gate-induced quantum-confinement transition of a single dopant atom in a silicon FinFET.  _Nature Physics_ 4, 656–661.  DOI: 10.1038/nphys994 (2008).
  * [12] van der Heijden, J. _et al._ Probing the Spin States of a Single Acceptor Atom.  _Nano Letters_ 14, 1492–1496.  DOI: 10.1021/nl4047015 (2014).
  * [13] Fuechsle, M. _et al._ A single-atom transistor.  _Nature Nanotechnology_ 7, 242–246.  DOI: 10.1038/nnano.2012.21 (2012).
  * [14] Prati, E. _et al._ Few electron limit of n-type metal oxide semiconductor single electron transistors.  _Nanotechnology_ 23, 215204.  DOI: 10.1088/0957-4484/23/21/215204 (2012).
  * [15] Pierre, M. _et al._ Single-donor ionization energies in a nanoscale CMOS channel.  _Nature Nanotechnology_ 5, 133–137.  DOI: 10.1038/nnano.2009.373 (2010).
  * [16] Steger, M. _et al._ Quantum Information Storage for over 180 s Using Donor Spins in a 28Si ”Semiconductor Vacuum”.  _Science_ 336, 1280–1283.  DOI: 10.1126/science.1217635 (2012).
  * [17] Veldhorst, M. _et al._ An addressable quantum dot qubit with fault-tolerant control-fidelity.  _Nature Nanotechnology_ 9, 981–985.  DOI: 10.1038/nnano.2014.216 (2014).
  * [18] Veldhorst, M. _et al._ A two-qubit logic gate in silicon.  _Nature_ 526, 410–414.  DOI: 10.1038/nature15263 (2015).
  * [19] Morton, J. J. L., McCamey, D. R., Eriksson, M. A. & Lyon, S. A.  Embracing the quantum limit in silicon computing.  _Nature_ 479, 345–353.  DOI: 10.1038/nature10681 (2011).
  * [20] Zwanenburg, F. A. _et al._ Silicon quantum electronics.  _Reviews of Modern Physics_ 85, 961–1019.  DOI: 10.1103/RevModPhys.85.961 (2013).
  * [21] Spruijtenburg, P. C. _et al._ Single-hole tunneling through a two-dimensional hole gas in intrinsic silicon.  _Applied Physics Letters_ 102, 192105.  DOI: 10.1063/1.4804555 (2013).
  * [22] Li, R., Hudson, F. E., Dzurak, A. S. & Hamilton, A. R.  Single hole transport in a silicon metal-oxide-semiconductor quantum dot.  _Applied Physics Letters_ 103, 163508.  DOI: 10.1063/1.4826183 (2013).
  * [23] Jovanovic, D., Leburton, J. P., Chang, H., Grundbacher, R. & Adesida, I.  Disorder-induced resonant tunneling in planar quantum-dot nanostructures.  _Physical Review B_ 50, 5412–5419.  DOI: 10.1103/PhysRevB.50.5412 (1994).
  * [24] Hofheinz, M. _et al._ Individual charge traps in silicon nanowires.  _The European Physical Journal B_ 54, 299–307 (2006).
  * [25] Prati, E., Hori, M., Guagliardo, F., Ferrari, G. & Shinada, T.  Anderson – Mott transition in arrays of a few dopant atoms in a silicon transistor.  _Nature Nanotechnology_ 7, 443–447.  DOI: 10.1038/nnano.2012.94 (2012).
  * [26] Nishi, Y., Tanaka, K. & Ohwada, A.  Study of Silicon-Silicon Dioxide Structure by Electron Spin Resonance II.  _Japanese Journal of Applied Physics_ 11, 85–91.  DOI: 10.1143/JJAP.11.85 (1972).
  * [27] Thoan, N. H., Keunen, K., Afanas’ev, V. V. & Stesmans, A.  Interface state energy distribution and Pb defects at Si(110)/SiO2 interfaces: Comparison to (111) and (100) silicon orientations.  _Journal of Applied Physics_ 109, 013710.  DOI: 10.1063/1.3527909 (2011).
  * [28] Mishima, T. & Lenahan, P.  A spin-dependent recombination study of radiation-induced Pb1 centers at the [001] Si/SiO2 interface.  _IEEE Transactions on Nuclear Science_ 47, 2249–2255.  DOI: 10.1109/23.903761 (2000).
  * [29] Schroder, D. K. & Babcock, J. A.  Negative bias temperature instability: Road to cross in deep submicron silicon semiconductor manufacturing.  _Journal of Applied Physics_ 94, 1.  DOI: 10.1063/1.1567461 (2003).
  * [30] Campbell, J. P., Lenahan, P. M., Cochrane, C. J., Krishnan, A. T. & Krishnan, S.  Atomic-Scale Defects Involved in the Negative-Bias Temperature Instability.  _IEEE Transactions on Device and Materials Reliability_ 7, 540–557.  DOI: 10.1109/TDMR.2007.911379 (2007).
  * [31] Stegner, A. R. _et al._ Electrical detection of coherent 31P spin quantum states.  _Nature Physics_ 2, 835–838.  DOI: 10.1038/nphys465 (2006).
  * [32] Xiao, M., Martin, I. & Jiang, H. W.  Probing the spin state of a single electron trap by random telegraph signal.  _Physical review letters_ 91, 078301.  DOI: 10.1103/PhysRevLett.91.078301 (2003).
  * [33] Xiao, M., Martin, I., Yablonovitch, E. & Jiang, H. W.  Electrical detection of the spin resonance of a single electron in a silicon field-effect transistor.  _Nature_ 430, 435–439.  DOI: 10.1038/nature02727 (2004).
  * [34] Prati, E., Fanciulli, M., Ferrari, G. & Sampietro, M.  Effect of the triplet state on the random telegraph signal in Si-n-MOSFETs.  _Physical Review B_ 74, 033309.  DOI: 10.1103/PhysRevB.74.033309 (2006).
  * [35] Ragnarsson, L.-Å. & Lundgren, P.  Electrical characterization of Pb centers in (100)Si–SiO2 structures: The influence of surface potential on passivation during post metallization anneal.  _Journal of Applied Physics_ 88, 938.  DOI: 10.1063/1.373759 (2000).
  * [36] Lannoo, M.  The role of dangling bonds in the properties of surfaces and interfaces of semiconductors.  _Revue de Physique Appliquée_ 25, 887–894.  DOI: 10.1051/rphysap:01990002509088700 (1990).
  * [37] Cartier, E., Stathis, J. H. & Buchanan, D. a.  Passivation and depassivation of silicon dangling bonds at the Si/SiO2 interface by atomic hydrogen.  _Applied Physics Letters_ 63, 1510.  DOI: 10.1063/1.110758 (1993).
  * [38] Angus, S. J., Ferguson, A. J., Dzurak, A. S. & Clark, R. G.  Gate-Defined Quantum Dots in Intrinsic Silicon.  _Nano Letters_ 7, 2051–2055.  DOI: 10.1021/nl070949k (2007).
  * [39] Mueller, F., Konstantaras, G., Spruijtenburg, P. C., van der Wiel, W. G. & Zwanenburg, F. A.  Electron–Hole Confinement Symmetry in Silicon Quantum Dots.  _Nano Letters_ 15, 5336–5341.  DOI: 10.1021/acs.nanolett.5b01706 (2015).
  * [40] Betz, A. C., Gonzalez-Zalba, M. F., Podd, G. & Ferguson, A. J.  Ambipolar quantum dots in intrinsic silicon.  _Applied Physics Letters_ 105, 153113.  DOI: 10.1063/1.4898704 (2014).
  * [41] Dingemans, G., Beyer, W., van de Sanden, M. C. M. & Kessels, W. M. M.  Hydrogen induced passivation of Si interfaces by Al2O3 films and SiO2/Al2]O3.  _Applied Physics Letters_ 97, 152106.  DOI: 10.1063/1.3497014 (2010).
  * [42] Müller, F.  _Single-charge tunneling in ambipolar silicon quantum dots_.  Ph.D. thesis, University of Twente, Enschede, The Netherlands.  DOI: 10.3990/1.9789036539036 (2015).
  * [43] Bierhals, A., Aberle, A. G. & Hezel, R.  Improved understanding of thermally activated structural changes in Al/SiOx/p-Si tunnel diodes by means of infrared spectroscopy.  _Journal of Applied Physics_ 83, 1371.  DOI: 10.1063/1.366839 (1998).
  * [44] Thompson, C. V.  Solid-State Dewetting of Thin Films.  _Annual Review of Materials Research_ 42, 399–434.  DOI: 10.1146/annurev-matsci-070511-155048 (2012).
  * [45] Hinode, K., Asano, I. & Homma, Y.  Void formation mechanism in VLSI aluminum metallization.  _IEEE Transactions on Electron Devices_ 36, 1050–1055.  DOI: 10.1109/16.24347 (1989).
  * [46] Schofield, S. R. _et al._ Quantum engineering at the silicon surface using dangling bonds.  _Nature Communications_ 4, 1649.  DOI: 10.1038/ncomms2679 (2013).
  * [47] Thorbeck, T. & Zimmerman, N. M.  Formation of strain-induced quantum dots in gated semiconductor nanostructures.  _AIP Advances_ 5, 087107\.  DOI: 10.1063/1.4928320 (2015).
  * [48] Kim, Y. Y. & Lenahan, P. M.  Electron-spin-resonance study of radiation-induced paramagnetic defects in oxides grown on (100) silicon substrates.  _Journal of Applied Physics_ 64, 3551.  DOI: 10.1063/1.341494 (1988).
  * [49] Fleetwood, D. M., Pantelides, S. T. & Schrimpf, R. D. (eds.) _Defects in Microelectronic Materials and Devices_ (CRC Press, 2008).
  * [50] Schiele, Y., Hahn, G. & Terheiden, B.  Investigation of Radiation Damage to the Al2O3 / Si Wafer Interface during Electron Beam Evaporation by Means of C-V and Lifetime Measurements.  in _EU PVSEC Proceedings_ , 1068–1072 (2011).  DOI: 10.4229/26thEUPVSEC2011-2CO.14.5.

## Acknowledgements

The authors would like to thank J.W. Mertens and A.A.I. Aarnink for the growth
of the $\textrm{Al}_{2}\textrm{O}_{3}$, and J. Schmitz for valuable discussion
and input. We acknowledge financial support through the EC FP7-ICT initiative
under Project SiAM No 610637. W.G. van der Wiel and F.A. Zwanenburg
acknowledge support from the Foundation for Fundamental Research on Matter
(FOM), which is part of the Netherlands Organization for Scientific Research
(NWO).

## Author contributions statement

P.C.S. conducted the experiment(s) and wrote the manuscript text, P.C.S.,
S.V.A., and F.M. fabricated the samples, P.C.S, S.V.A, and F.A.Z. analysed the
results. All authors reviewed the manuscript.

## Additional information

All authors declare that there are no potential competing financial interests.