# Quantum graph models for transport in filamentary switching

Alison A. Silva Programa de Pós-Graduação em Ciências/Física, Universidade
Estadual de Ponta Grossa, 84030-900 Ponta Grossa, Paraná, Brazil  Fabiano M.
Andrade Programa de Pós-Graduação em Ciências/Física, Universidade Estadual de
Ponta Grossa, 84030-900 Ponta Grossa, Paraná, Brazil  Departamento de
Matemática e Estatística, Universidade Estadual de Ponta Grossa, 84030-900
Ponta Grossa, Paraná, Brazil  Francesco Caravelli Theoretical Division (T4),
Los Alamos National Laboratory, Los Alamos, New Mexico 87545, USA

###### Abstract

The formation of metallic nanofilaments bridging two electrodes across an
insulator is a mechanism for resistive switching. Examples of such phenomena
include atomic synapses, which constitute a distinct class of memristive
devices whose behavior is closely tied to the properties of the filament.
Until recently, experimental investigation of the low-temperature regime and
quantum transport effects has been limited. However, with growing interest in
understanding the true impacts of the filament on device conductance,
comprehending quantum effects has become crucial for quantum neuromorphic
hardware.

We discuss quantum transport resulting from filamentary switching in a narrow
region where the continuous approximation of the contact is not valid, and
only a few atoms are involved. In this scenario, the filament can be
represented by a graph depicting the adjacency of atoms and the overlap
between atomic orbitals. Using the theory of quantum graphs with node
scattering which is locally diffusive, we calculate the scattering amplitude
of charge carriers on this graph and explore the interplay between filamentary
formation and quantum transport effects.

##  I Introduction

There has been a huge interest in resistive switching and memristive behavior,
as the complexity of many nanoscale devices requires new techniques to
understand their behavior.

For instance, memristors offer the possibility of harnessing both nonlinear
behavior and non-trivial memory in electronic circuits. Specialized circuits
composed of large numbers of such devices promise a new generation of
computational hardware operating orders of magnitude faster, and at far lower
power, than traditional digital circuitry [1, 2, 3, 4, 5]. Since memristive
behavior is associated to a 1-port device which is current- or voltage-
controlled, these components support volumetric memory scaling for dense
storage capacity, with the possibility of being simply operated on a crossbar
array [6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21]. Moreover,
recent results suggest that networks of memristive devices [22, 24, 23] can be
used for various tasks in analog computing, and exhibit chaos and system
switching phenomena [25, 26].

One advantage of these nanoscale devices is that with a displacement of a few
atoms, we have possibly access to resistive states at ultimately low energy
costs. However, device engineering lacks the complete quantum theoretical
characterization of resistive switching when this is driven by a conductance
change due to a filament formation, also considering that the nature of the
filament formation is inherently stochastic. Although some analytical methods
exist, these are restricted to the semiclassical domain [27].

Filament formation, e.g. the nonequilibrium process of chemically forming a
filament in a narrow region [28, 29, 30, 31, 32], is a mechanism for resistive
devices. However, with the increasing capability of scaling down to a few
atoms, and the ability to measure quantum effects in these filaments at a few
Kelvins, new techniques are necessary to capture the phenomenology of these
devices.

Quantum transport, on the other hand, has been extremely successful in the
study of nanoscale devices [72, 73, 74], whether using the Kubo [34] or
Landauer approach [35, 36, 37, 38, 39, 40, 41, 42, 43]. For the case of non-
filamentary switching, nanoimaging techniques have shown recently that charge
density waves emerge when two carriers of opposite charge are allowed in the
narrow region, for instance in the case of perovskite nickelates doped with
Hydrogen [33]. The present paper focuses in particular on nanoscale
filamentary switching. Recently, state-of-the-art experimental results
managed, via Andreev reflection processes [44, 45, 46, 47, 48], to measure
quantum effects in niobium pentoxide [49, 54, 55, 56, 57, 58], emerging at the
filament terminals when superconducting electrodes are utilized [60, 61, 62,
63, 64, 65, 59]. In particular, in certain cases, the quantum PIN code, i.e.
the transmission probabilities of each conduction channel contributing to the
conductance of the nanojunctions can be obtained [66] using Andreev
reflections. These studies suggest that atomic-size metallic filaments are
indeed the source of high conductance in niobium pentoxide (the so-called ON
state). These recent studies motivate the theoretical investigation of
transport in heterogeneous filamentary structures [28, 69, 70, 71].

However, immediately one faces an obstruction. Previous studies focusing on
filament formation, directly simulating filament formation, have used either
continuum models and/or well-defined lattices, or classical Kirchhoff and
Ohm’s laws to calculate the resistance. On the other hand, quantum mesoscopic
transport techniques are based, typically, on the assumption of continuity,
e.g. the conductor acts as a waveguide for the charge carriers [72, 73, 74].
As such, the study of quantum transport on this mesoscopic network is a
challenging task already at the modeling level, even before any calculation
can be attempted.

Figure 1: Sketch of filamentary switching for a few atoms. The structure is
far from continuous, and the charge carriers hop across the structure as
induced by the voltage. The formalism we use is the one by Landauer-Büttiker,
with two reservoirs (leads).

For this reason, the general approach we consider in this paper is represented
in Fig. 1. We consider a left and right lead and the dynamical formation of a
filament in the region. We assume that the filament formation is classical,
and on a timescale much longer than the effective speed of the electrons
across the filament. The filament is composed of an assembly of $N$ atoms,
with overlapping bands. However, the scattering process of the electron across
the filament is quantum mechanical, and thus one has to consider the
scattering amplitude $T_{ij}$ from the lead $i$ to the lead $j$ quantum
mechanically.

However, over the past three decades, the study of electrical conductors has
evolved considerably from a macroscopic level of description to the nanoscale.
Assuming that the characteristic length of the quantum correlations is more or
less of a few atoms (e.g. dephasing or Fermi lengths are small), we can assume
coherent quantum transport of electrons. The setup we study is the typical one
of mesoscopic transport [34], e.g. our system has two leads, denoted as left
and right, and a voltage applied to them. Such a picture has been extremely
successful in the study of molecular and atomic point contacts [75, 76, 77,
78, 79, 80], carbon nanotubes [81, 82] or graphene [83, 84, 85, 86, 87]. More
recently, atomic point contacts have been also investigated for resistive
switching [88, 89] and there has been an interest in generic quantum wires
[90] or Josephson junction networks [91].

The idea that the discreteness and yet heterogeneity of the structure might
play a role has been investigated as early as the 30’s by Pauling [92, 93, 94,
95]. Thankfully, the study of scattering on discrete structures has obtained
some considerable attention over the past decades (the so-called quantum
graphs) [96, 97, 98, 99, 100, 101, 102, 103, 104, 105, 106, 107, 108, 109],
driven by the study of quantum chaotic systems. Mathematically, the overlap
between the atoms can be described as a graph [110, 111], with the atoms
represented as vertices and the overlap as edges of such graph.

The assumption is that electrons on the left or right leads are asymptotically
free and incoming and outgoing states [112, 113, 114, 115, 116]. The electrons
can propagate freely or in a constant potential $u_{ij}$ along the edge of the
graph (a one dimensional structure), whose distance is assumed to be constant
$\ell$. Then, the problem of transport can be described as the problem of
calculating the scattering matrix which transforms the incoming into outgoing
free states on these “quantum” graphs [117, 118, 119, 120, 121, 122, 123, 124,
125, 126, 127]. Intuitively, Kirchhoff’s currents conservation for all the
nodes of a circuit, and which is classically written in terms of the directed
incidence matrix of a graph as $B^{t}=\vec{i}=0$ [110, 24], has to be replaced
by a “quantum” Kirchhoff law in which a unitary matrix has to be inserted in
the sum over currents at a node. In quantum graphs, these unitary operators
can be obtained by the scattering process a local one-body Hamiltonian.

For an entire network connecting two leads, such a unitary matrix contains the
information about the transmission amplitudes; at this point, we can use the
transmission amplitudes within the context of a Landauer approach. In fact,
the probably best-known result in the theory of quantum transport obtained
using the scattering matrix approach is the Landauer formula [35], who was the
first to use scattering matrices to describe transport problems. For multi-
leads, we use a generalization of this approach, also known as the Landauer-
Büttiker formula. Specifically, we will treat our system as a particular case
of a multi-lead system, where however some of the leads are identical. In the
case of several channels, the expression for the conductance

$G=\frac{2e^{2}}{h}\sum_{n,m}T_{nm}$ (1)

contains the sum of transmission probabilities $T_{nm}$ from one mode
(channel) to another. The expression above will have to be adapted to the case
of quantum graphs which we discuss below. In the more specific case of a
filament, we will see that we can have $K$ channels attached to the left and
$K^{\prime}$ channels attached to the right lead.

The approach we use here is different from the case of the non-equilibrium
Green’s function [73, 74] in various ways. First, the scattering process is on
a discrete structure rather than a continuous one, as it is commonly assumed
in analytical calculations [129, 130, 131, 132]. Moreover, we make the
simplifying assumption that the electrons are faster than the typical
timescale of the filament formation. As such, we consider the equilibrium
Green’s functions on a particular structure. This is motivated by the fact
that the typical formation of a filament in these materials is of the order of
a few milliseconds, in which this annealed approximation makes sense.

The paper is organized as follows. We first review the standard Landauer-
Büttiker approach for the multichannel scattering. We then introduce the setup
of this paper and the construction of the scattering matrix based on quantum
graphs. We then discuss the molecular dynamics simulations that we employed
for the stochastic filament formation, and give the assumptions on the
dielectric. Conclusions follow.

##  II Standard Landauer approach and graph scattering matrix

The approach we consider here is the one introduced in [39].

###  II.1 Landauer-Büttiker formalism

To describe a quasi-one-dimensional coherent conductor, we first consider a
purely one-dimensional problem [73, 74]. We consider for instance electrons
with energies up to a certain scale $\mu$ in a scattering state, moving (for
instance) left to right, in a Lippmann-Schwinger scattering states [133, 134]
of the form

$\Psi_{L,E}(x)=\begin{cases}e^{ikx}+r(E)e^{-ikx},&x\to-\infty,\\\
t(E)e^{ikx},&x\to\infty,\end{cases}$ (2)

where we introduced $k=\sqrt{2mE}/\hbar$, and we assume the state above to be
normalized. The quantities $r(E)$ and $t(E)$ are the reflection and
transmission amplitudes respectively. For right-moving states, we have an
analogous state to the one of eqn. (2). The density of states in one dimension
is given by the expression $\nu(E)=\frac{dk}{dE}=\frac{m}{\hbar^{2}k}$.
Because of the electron statistics, the probability of having an electron at a
certain energy $E$, is given by

$f(E)=\frac{1}{e^{(E-\mu)/\kappa T}+1}.$ (3)

where $\mu$ is a chemical potential, which in the Landauer approach is assumed
to be dependent on the voltage potential at the lead. If we have two leads,
then, the left and right leads depend on two chemical potentials $\mu_{L}$ and
$\mu_{R}$, which we assume to be dependent on the voltage at the two
electrodes, anodes and cathodes. The total current is the sum of the left and
right contribution, given by the expression [112, 113, 114, 115, 116]

$\displaystyle I(V)$ $\displaystyle=$ $\displaystyle I_{L}(V)+I_{R}(V)$ (4)
$\displaystyle=$
$\displaystyle\frac{2e}{h}\int_{-\infty}^{\infty}dE(f_{L}(E)-f_{R}(E))T(E),$

where $\frac{2e}{h}=G_{0}/e$, where $G_{0}$ is the conductance quantum. The
transmission probability from left to right, $T=|t|^{2}$, is equal to the
transmission probability from right to left, $T=T^{\prime}=|t^{\prime}|^{2}$
because of unitarity, and then in this sense, all we need to know is the
scattering amplitudes for the system given the left and right scattering.

The quantity $T(E)$ is the sum over all channels,
$T(E)=\text{Tr}(\tau^{\dagger}(E)\tau(E))$, where $\tau(E)$ is the
transmission matrix of the system.

We see that within this approach, we only need to solve the scattering process
across the structure, whether this is a continuum or discrete. As we have said
above, in this paper we consider only local effective one-dimensional
movements of the electrons in between the underlying atoms forming the
filaments. For this purpose, then, we need to introduce the graph scattering
process associated with our generalized Lippman-Schwinger states.

In general, if the system has many leads, the total transition probability can
be written as

$\displaystyle I_{k}(V)$ $\displaystyle=$
$\displaystyle\frac{G_{0}}{e}\int_{-\infty}^{\infty}dE(f_{L}(E)-f_{R}(E))\sum_{j}T_{kj}(E).$
(5)

where $T_{kj}$ is the transition probability from the lead $j$ to the lead
$k$.

If the leads $k$ are then attached to the same contact, we can write

$\displaystyle
I_{tot}(V)=\frac{G_{0}}{e}\int_{-\infty}^{\infty}dE(f_{L}(E)-f_{R}(E))\sum_{jk}T_{kj}(E).$
(6)

Since in quantum graphs we do obtain the quantities $T_{kj}$, we will use this
expression. At small temperatures and for $\mu_{L}-\mu_{R}\approx eV$, this
reduces to

$\displaystyle
I_{tot}(V)=\frac{G_{0}}{e}e\sum_{jk}T_{kj}(E_{f})V=G_{0}(\sum_{jk}T_{kj}(E_{f}))V,$
(7)

where $E_{f}$ is the Fermi energy for the system and $V$ the applied voltage
$G_{0}\approx$ (12900 Ohms)-1, the quantum of conductance. This equation
cannot be directly applied in the theory of quantum graphs in a sense we will
see in the next subsection, but a simple modification will be sufficient to
allow us to use it.

###  II.2 Adapting the LB formula to quantum graphs

Before we discuss the theory of quantum graphs, it is worth discussing what
one can calculate and how the results should be interpreted. Quantum graphs
are a generalization of the one-dimensional Schrödinger equation, in which
single particles enter along a one-dimensional path between the scattering
nodes, and the electron can move across a heterostructure. In this sense,
particles do not have any transverse modes, and the scattering process is
purely one-dimensional. The starting point is a certain graph representing the
atomic structure and the overlapping orbitals. In such a graph, the nodes
represent where the particles scatter and the edges are the one-dimensional
paths in between scattering events.

What one can calculate with this mathematical construct is the transmission
probability $T_{ij}(E)$, e.g. the chances that the particle being coherent
along the mesoscopic region enters the region through a certain node $i$ and
leaves from node $j$. Now let us assume that our scattering region,
approximated as a graph, is attached via $N$ nodes on the left lead and $M$
nodes on the right lead.

If the particle incomes from the left lead, and want to calculate the right-
moving current $I^{L}$, then we assume that will have classical probability
$p^{L}_{i}(E)$ to enter at node $i$, satisfying
$\sum_{i=1^{N}}p^{L}_{i}(E)=1$. Then, the total current, following a similar
argument to the one of Landauer, is given by

$\displaystyle I^{L}(V)=G_{0}/e\int
dEf_{L}(E,V)\sum_{i=1}^{N}\sum_{j=1}^{M}T_{ij}(E)p_{i}^{L}(E).$ (8)

A similar argument can be made for the left moving particles incoming from the
right lead, obtaining

$\displaystyle I^{R}(V)=-G_{0}/e\int
dEf_{R}(E,V)\sum_{i=1}^{N}\sum_{j=1}^{M}T_{ij}(E)p_{j}^{R}(E),$ (9)

from which we get the total current expression

$\displaystyle I(V)$ $\displaystyle=$ $\displaystyle I_{L}(V)+I_{R}(V)$
$\displaystyle=$ $\displaystyle G_{0}/e\int
dE\sum_{i=1}^{N}\sum_{j=1}^{M}\Big{(}f_{L}p_{i}^{L}-f_{R}p_{j}^{R}\Big{)}T_{ij}$

where we have suppressed the functional dependence for simplicity. We note
that this expression is similar to the one of the Landauer-Büttiker formula,
except for the correcting factors $p_{i}^{L}$ and $p_{j}^{R}$. The method
above is, however, time-consuming, as it requires the evaluation of $ML$
channels’ probabilities. We use instead the following trick, commonly used in
the literature [135, 136, 137] of quantum dragons to study the transmissions
to a single channel. Instead of repeating the process $ML$ times, we attach
two extra nodes to the graph attached to the boundary nodes (connected to the
leads). Using this method, we reduce our problem to a single Lippman-Schwinger
asymptotic in and out states. This is shown in Fig. 2.

Figure 2: Reduction of the lead connections to a single channel and incoming-
outgoing Lippman-Schwinger asymptotic states.

###  II.3 Quantum graphs: quantum scattering on graphs

In the Landauer-Büttiker approach all one needs to know is the scattering
amplitudes on the underlying structure. As mentioned earlier, we take the
point of view in which the scattering occurs along one-dimensional paths
connecting the underlying atoms. This is a departure from previous works in
quantum transport, but it can be considered as a model approximation to
calculate the scattering amplitudes as a function of the incoming energy. It
is thus worth formalizing the process, as it is not a standard approach in
transport. The goal of this and the next few sections is to introduce the
formalism that we use.

A graph $\mathcal{G}(V,E)$ is the tuple given by the set of vertices vertices
$V(\mathcal{G})=\\{1,\ldots,n\\}$ and the set of edges between nodes,
$E(\mathcal{G})=\\{e_{1},\ldots,e_{l}\\}$, where each edge is a pair of
vertices [110].

Figure 3:  Locally, any graph looks like a star graph.

In the near vicinity of any node, a quantum graph always looks like a star
graph of the type shown in Fig. 3.

The graph connectivity is encoded in the $n\times n$ adjacency matrix
$A(\mathcal{G})$, $A_{ij}$ is one if the vertices $i$ and $j$ are connected
and zero otherwise. The set $E_{i}=\left\\{j:\\{i,j\\}\in
E(\mathcal{G})\right\\}$ is the neighborhood of the vertex $i\in
V(\mathcal{G})$ (e.g. at distance one). We denote with
$E_{i}^{k}=E_{i}\setminus\\{k\\}$ the set of neighbors of the vertex $i$
except for $k$. The degree of the vertex $i$ is defined as
$d_{i}=|E_{i}|=\sum_{j=1}^{n}A_{ij}(\mathcal{G})$. These definitions are
standard in graph theory and are called simple or discrete graphs.

However, since our graph is embedded in space, it is useful to also consider
metric graphs, in which the notion of distance is encoded. To discuss quantum
graphs, it is necessary to equip the graphs with a metric. A metric graph
$\Gamma(V,E)$ is a graph in which every edge is assigned a certain length
$\ell_{e_{s}}\in(0,+\infty)$, defining
$\boldsymbol{\ell}=\\{\ell_{e_{1}},\ldots,\ell_{e_{l}}\\}$.

Given this definition, we call leads those edges which have a semi-infinite
length $e_{s}$ ($\ell_{e_{s}}=+\infty$); these edges are associated with the
Lippman-Schwinger states (in the literature, these are also called “open
quantum graphs").

####  II.3.1 Basic theoretical framework

Quantum graphs are metric graphs in which the Schrödinger operator is defined
along the edges of the graph. At every vertex, we have boundary conditions
(BCs). The triple $\\{\Gamma(V,E),\mathbb{H},\mathbb{bc}\\}$ with $H$ a
differential operator and $\mathbb{bc}$ a set of BCs defines a quantum graph.
For instance, along a certain edge $\\{i,j\\}$, the free Schrödinger operator
$H_{ij}=-(\hbar^{2}/2m)d^{2}/d{x_{e_{ij}}}^{2}$ is associated to the
eigenvalue equation on the edge $\\{i,j\\}$ is given by

$-\psi_{\\{i,j\\}}^{\prime\prime}(x)=k^{2}\psi_{\\{i,j\\}}(x),$ (11)

where $k=\sqrt{2mE/\hbar^{2}}$, $m$ is the mass of the particle, $E$ is its
energy, and $\psi_{\\{i,j\\}}$ is the wave function. We call naked quantum
graphs the triple in which $H$ is the free evolution, e.g. there is no
potential, and dressed quantum graphs those in which edge potentials are
present.

In general, formally, the total wave function of the system is a collection of
local wavefunctions assigned to the edges. We have

$\displaystyle\Psi=\begin{pmatrix}\psi_{e_{1}}(x_{e_{a_{1}b_{1}}})\\\
\vdots\\\ \psi_{e_{m}}(x_{e_{a_{m}b_{m}}})\\\ \end{pmatrix}$ (12)

where the pairs $(a_{m}b_{m})$ represent an edge, e.g. a pair of nodes, and
$m$ is the total number of edges. In general, $x_{e_{ij}}$ is a local
parameter such that for $x=0$ we have the wave function on the vertex $i$,
while $x=\ell$ we have the wavefunction on the vertex $j$. We will often write
$e_{k}$ to denote a certain edge, dropping the two vertices notation unless
necessary. The Hamiltonian operator

$\displaystyle
H_{e_{k}}(x_{e_{j}})=-\frac{\hbar^{2}}{2m}\frac{d^{2}}{dx_{e_{j}}^{2}}+V(x_{e_{j}})$
(13)

is the standard 1-d particle Hamiltonian on the compact domain $[0,\ell_{j}]$.

In the theory of quantum graphs, vertices are interpreted as zero-range
potentials, where boundary conditions need to be supplied. The collection of
boundary conditions is, for every node, the continuity and differentiability
of the wave functions. If $\\{e_{1},\cdots,e_{k}\\}$ are the edges adjacent to
a certain vertex $n$, then we must have

$\displaystyle(1)$ $\displaystyle\ \ \ \ \
\psi_{e_{1}}|_{n}=\psi_{e_{2}}|_{n}=\cdots=\psi_{e_{k}}|_{n}$ (14)
$\displaystyle(2)$ $\displaystyle\ \ \ \ \
\sum_{j=1}^{k}\partial_{x_{e_{j}}}\psi_{e_{j}}|_{n}=0$ (either at
$x_{e_{j}}=0$ or $x_{e_{j}}=\ell_{j}\\}$ depending on the definition

In the case of zero potential for the edge, the wavefunction is generically of
the form

$\displaystyle-$
$\displaystyle\frac{d^{2}\psi_{e_{k}}}{dx_{e_{k}}^{2}}=k^{2}\psi_{e_{k}}$ (16)
$\displaystyle\ \ \
\rightarrow\psi_{e_{j}}(x_{e_{j}})=c^{+}_{j}e^{ikx_{e_{j}}}+c^{-}_{j}e^{-ikx_{e_{k}}}$

A Lippman-Schwinger state is a particular choice of an edge as in the incoming
“left" state, and accordingly a right state [133]. For instance, a “left"
(incoming) state on the edge $i$ and a “right" (outgoing) state on the edge
$f$ is given by the pair

$\psi_{LS,\\{e_{i},e_{f}\\}}=\begin{cases}e^{ikx_{i}}+r(E)e^{-ikx_{i}},&x_{i}\to-\infty,\\\
t(E)e^{ikx_{f}},&x_{f}\to\infty.\end{cases}$ (17)

We write now the left and right in quotation simply because our system is not
simply one dimensional any more, but the notation stands.

Eqs give a solution of the quantum graph equations. (14)-(16) plus the
external boundary conditions given by the Lippman-Schwinger asymptotic states
of Eq. (17); these are called “open” quantum graphs, and are represented in
Fig. 4. In this state, the quantum graph formalism is equivalent to the
scattering approach necessary for the Landauer-Büttiker approach, or Eq. (4).
While Green’s function approach is unitarily equivalent to the Schrödinger
approach for quantum graphs, the former can be automatized by using the
adjacency matrix, which is the technique we use here to study transport
numerically and we will review below.

Figure 4: Graph with two leads added turning it into an open quantum graph.

####  II.3.2 Green’s function approach

An important quantity in the theory of quantum scattering in general, and in
particular in the theory of quantum graphs, is the Green’s function approach
(GFA). The Green’s function for a certain edge is defined as

$\displaystyle[E-H(x_{f})]G(x_{f},x_{i};E)=\delta(x_{f}-x_{i}),$ (18)

which is the Green’s function of a 1-dimensional quantum particle. Formally,
the solution of the equation above can be written exactly, and we invite you
to read the review [101] on the subject. We only wish to add that there are a
few methods to calculate Green’s function exactly, including the exact
solution of the differential equation, using a spectral representation of the
exact solution, or alternatively using a Feynmann path integral
representation.

The approach that we emphasize here is the path integral approach. It is in
fact remarkable that in fact for quantum graphs, the semiclassical approach
introduced by Gutzwiller, or also Gutzwiller trace formula, is exact [101]. As
such, Green’s function can be written exactly as the sum over all classical
periodic orbits. Specifically, the Green’s function can be written as [106,
105] (see also [101])

$\displaystyle G(x_{f},x_{i};E)=\frac{m}{i\hbar^{2}k}\sum_{\rm sp}W_{\rm
sp}e^{\frac{i}{\hbar}S_{sp}(x_{f},x_{i};k)},$ (19)

where $\rm sp$ stands for scattering paths: these are all the paths that begin
in $x_{i}$ and end in $x_{f}$. Clearly, the sum above is possibly over an
infinite number of paths, given by the fact that a particle can be in
principle reflected in a countable many ways on an edge before leaving the
graph. This sum can be performed exactly. The quantity $W_{\rm sp}$ is the
quantum amplitude of the path: for each time a particle is reflected, it is
multiplied by a factor $r_{j}(E)$ on the specific edge if it is reflected, and
if transmitted it is multiplied by $t_{j}(E)$. The quantity $S_{\rm sp}$ is
instead the classical action of the particle over that specific path: for each
edge $j$, the particle has traversed along the path, we need to add
$S_{j}=\ell_{j}k$. Let us make a quick example which shows why this is useful.

Let us use now the following notation: for a vertex $j$, $r_{j}$, and $t_{j}$
are the reflection and transmission probabilities associated with the
wavefunction parametrization. As an example, let us look at Fig. 5. We have
two incoming and outgoing Lippman-Schwinger states, and one internal edge
bordered by the vertices $A$ and $B$, whose distance we assume to be simply
$l$. It is not hard to see that any scattering path can be written in the
forms of $sp_{1},\cdots,sp_{3},\cdots$: if a particle is reflected at $A$,
then because the edge $i$ is infinite it will never return back, and it cannot
be part of a scattering path.

The amplitude $W_{1}$ associated with $sp_{1}$ of Fig. 5 (top) is simply
$t_{A}t_{B}$, while the classical action is given by
$S_{sp1}=kl+kx_{i}+kx_{f}$. For $sp_{2}$, we have
$W_{2}=t_{A}r_{B}r_{A}t_{B}$, and $S_{sp1}=3kl+kx_{i}+kx_{f}$. For $sp_{3}$,
we have $W_{2}=t_{A}(r_{B}r_{A})^{2}t_{B}$, and $S_{sp1}=5kl+kx_{i}+kx_{f}$.
Any path can be written in this form, and it is not hard to see that the
Green’s function can be written as a geometrical series. We obtain then

$\displaystyle
G_{fi}(x_{f},x_{i};k)=\frac{m}{i\hbar^{2}k}\frac{t_{A}t_{B}e^{ikl}}{1-r_{A}r_{B}e^{2ikl}}e^{ik(x_{f}+x_{i})}.$
(20)

What is interesting about this expression is that we can read out the
transmission amplitude for the whole system, $T_{fi}$, given by

$\displaystyle T_{fi}=\frac{t_{A}t_{B}e^{ikl}}{1-r_{A}r_{B}e^{2ikl}}.$ (21)

At this point, we see that this total scattering amplitude has the exact same
form as the as the scattering of a quantum particle in one dimension in a
potential $V(x)$ of the form of Fig. 5 (bottom). This example shows that the
quantum graph technique is a legitimate option to calculate the scattering
amplitude in a mesoscale system under the single-body approximation, although
it also provides a generalization to more complex structures. In general, the
Green’s function can then be written as

$\displaystyle
G_{fi}(x_{f},x_{i};k)=\frac{m}{i\hbar^{2}k}T_{\Gamma}e^{ik(x_{f}+x_{i})}.$
(22)

where $T_{\Gamma}$ is the transmission, which is what we aim to calculate.

In what follows, we will also use another technique to evaluate the total
transmission probability which has been obtained recently and is based on the
adjacency matrix of the graph.

Figure 5: Top: The simplest graph which has a nontrivial transmission
probability. The quantities $t_{A}$ and $t_{B}$ are the transmission
amplitudes for the incoming wave function, while $r_{A}$ and $r_{B}$ are the
reflection amplitudes. The two “open" edges, on the left and right, are the
semi-infinite leads associated with the asymptotic Lippman-Schwinger states.
Calculating the Green’s function requires the enumeration of all possible
paths connecting the two asymptotic states. Bottom: The equivalent shape of
the potential associated with the top graph.

####  II.3.3 Adjacency matrix approach

One of the main ingredients of the GFA is the local origin of the total
scattering amplitude and the fact that the scattering at each vertex can be
defined locally.

For each vertex, we generalize the local reflection and transmission
amplitudes into a tensor, $\boldsymbol{\sigma}_{j}$, depending on each vertex
$j$ label. Locally, the scattering amplitudes are in fact those of a star
graph, a property which has been very fruitful in the study of quantum graphs.
More details on the connection between the boundary conditions at each star,
quantum probability flux conservation, and self-adjointness can be found in
[106]. We will now take a more practical approach, introduced in [127]. For a
star graph, let the local wavefunction be given by
$\Psi(j)=\left(\psi_{\\{j,1\\}}(j),\ldots,\psi_{\\{j,n\\}}(j)\right)^{T}$. The
most general BCs that are consistent with the self-adjoint condition [101] are
totally defined by two $d_{j}\times d_{j}$ matrices. At each node, in fact,
given two matrices $\mathcal{A}_{j}$ and $\mathcal{B}_{j}$, the boundary
conditions we defined earlier can be generalized to

$\mathcal{A}_{j}\Psi(j)+\mathcal{B}_{j}\Psi^{\prime}(j)=0.$ (23)

Note that $\mathcal{A}_{j}$ is not the adjacency matrix of the graph, but
represents the nodal boundary condition.

The only requirement is that the matrix $\mathcal{A}\mathcal{B}^{*}$ is self-
adjoint, while the matrix $\left(\mathcal{A}_{j},\mathcal{B}_{j}\right)$ has
the maximal rank $d_{j}$ [127]. We now want to see that defining a local
scattering matrix is equivalent to the boundary conditions. The boundary
conditions of eqns. (23) can be in fact determined by considering a plane wave
on the edge $e_{i,j}$ incoming into the vertex $j$, which we assume to be of
degree $d_{j}$.

Figure 6: Transmission and reflection probabilities for an incoming wave from
$i\rightarrow j$.

Let us introduce the matrix $\boldsymbol{\sigma}_{j}$, which is related to the
reflection and scattering amplitudes as
$\sigma_{j}^{[\\{j,i\\},\\{i,j\\}]}(k)=r_{j}^{[\\{j,i\\},\\{i,j\\}]}(k)$ and
$\sigma_{j}^{[\\{j,p\\},\\{i,j\\}]}(k)=t_{j}^{[\\{j,p\\},\\{i,j\\}]}(k)$,
where $r_{j}^{[\\{j,i\\},\\{i,j\\}]}$ is the reflection in vertex $j$ for the
wave incoming from $i$, while $t_{j}^{[\\{j,p\\},\\{i,j\\}]}(k)$ the
transmission through vertex $j$, for the wave incoming from $i$ and scattering
to $p$. Unitarity typically requires that these are not independent
quantities.

Given the definition above, the “local" Lippman-Schwinger states at the node
$j$ are given by

$\displaystyle\psi_{\\{i,j\\}}(x)={}$ $\displaystyle
e^{-ikx}+\sigma_{j}^{[\\{j,i\\},\\{i,j\\}]}(k)e^{ikx},$
$\displaystyle\psi_{\\{j,l\\}}(x)={}$
$\displaystyle\sigma_{j}^{[\\{j,l\\},\\{i,j\\}]}(k)e^{ikx}.$ (24)

A direct application of eqn. (23) gives

$\boldsymbol{\sigma}_{j}(k)=-(\mathcal{A}_{j}+ik\mathcal{B}_{j})^{-1}(\mathcal{A}_{j}-ik\mathcal{B}_{j}).$
(25)

A few comments are in order. First, one can either set the matrix
$\boldsymbol{\sigma}_{j}$ or the boundary conditions [106]. Second, the
scattering matrix depends on $k$, and thus on the energy
$E_{k}=\frac{\hbar^{2}k^{2}}{2m}$, in a non-trivial manner (although certain
boundary conditions are $k$-independent [104]).

The conservation of probability is encoded in the unitarity of the matrix
$\boldsymbol{\sigma}_{j}(k)$, e.g.
$\boldsymbol{\sigma}_{j}(k)\boldsymbol{\sigma}_{j}^{\dagger}(k)=\mathbbm{1}$,
with $\boldsymbol{\sigma}_{j}(k)=\boldsymbol{\sigma}_{j}^{\dagger}(-k)$. These
conditions are equivalent, in matricial form, to

$\displaystyle\sigma_{j}^{[\\{j,l\\},\\{i,j\\}]}(k)=\left[\sigma_{j}^{[\\{i,j\\},\\{j,l\\}]}(-k)\right]^{*},$
$\displaystyle\sum_{i\in
E_{j}}\sigma_{j}^{[\\{j,l\\},\\{i,j\\}]}(k)\left[\sigma_{j}^{[\\{j,m\\},\\{i,j\\}]}(k)\right]^{*}=\delta_{lm},$
(26) $\displaystyle\sum_{i\in
E_{j}}\sigma_{j}^{[\\{i,j\\},\\{j,l\\}]}(k)\left[\sigma_{j}^{[\\{i,j\\},\\{j,m\\}]}(k)\right]^{*}=\delta_{lm},$

which can be seen as a generalization of one-dimensional scattering amplitudes
[107]. For instance, for a node of degree 2, we can write

$\displaystyle\boldsymbol{\sigma}(k)=\begin{pmatrix}r^{+}(k)&t^{-}(k)\\\
t^{+}(k)&r^{-}(k),\end{pmatrix}$ (27)

and the conditions
$\boldsymbol{\sigma}(k)\boldsymbol{\sigma}(k)^{\dagger}=\boldsymbol{\sigma}(k)^{\dagger}\boldsymbol{\sigma}(k)=\boldsymbol{1}$,
and $\boldsymbol{\sigma}(-k)=\boldsymbol{\sigma}(k)^{\dagger}$, reduce to

$\displaystyle r^{+}(k)r^{-}(k)+t^{+}(k)r^{-}(k)$ $\displaystyle=$
$\displaystyle 0$ (28)

Referring to Fig. 6, we must have at node $j$ and for the incoming wave
$i\rightarrow j$, that

$\displaystyle|\sigma_{j}^{\\{j,i\\},\\{i,j\\}}|^{2}+\sum_{p=1}^{d_{j}-1}|\sigma_{j}^{\\{j,p\\},\\{i,j\\}}|^{2}=1.$
(30)

Let us now introduce the adjacency matrix technique for the evaluation of the
Green’s function, eqn. (22). The goal of this section is to state an
alternative formulation based on the adjacency matrix $A_{ij}$ (not to be
confused with $\mathcal{A}_{j}$), instead of defining it in terms of the
scattering paths. Such a formulation is more compact and can be written,
instead of having to enumerate all possible scattering paths, by the solution
of a linear problem.

The construction of the total scattering amplitude $T_{\Gamma}$ is equivalent
to the following construction, based on the adjacency matrix. Given the
adjacency matrix $A(\mathcal{G})$, the Green’s function between the incoming
edge $e_{i}$ to the outgoing edge $e_{f}$, can be written as [127]

$G_{\Gamma,fi}=\frac{m}{i\hbar^{2}k}T_{\Gamma}e^{ik(x_{i}+x_{f})},$ (31)

where, given $\mathcal{N}_{i}$ the neighborhood of the node $i$ (the nodes at
distance one) we have that the total transmission probability is given by

$T_{\Gamma}=\sum_{j\in\mathcal{N}_{i}}\sigma_{i}^{[\\{i,j\\},e_{i}]}A_{ij}p_{ij},$
(32)

where $i$ is the entry node and where the definition of $p_{ij}$ will be given
in a moment. The quantity $\sigma_{i}^{[\\{i,j\\},e_{i}]}$ is the transmission
probability incoming from the “open" edge $e_{i}$, incoming to node $i$, and
scattering into edge $j$. Similarly, we define
$\sigma_{n}^{[e_{f},\\{i,n\\}]}$ as the transmission probability at the node
$n$ from the edge $\\{i,n\\}$ into the ”open" edge $e_{f}$.

Let us note that in the expression above, the sum over the scattering paths
has disappeared, and we are left with an expression that depends only on the
transmission and reflection probabilities, which are equivalent to the
boundary conditions, as we have seen. We will assume these as given for the
time being. The equivalence between the Green’s function via the Gutzwiller
formula relies upon the following assumptions [127]:

  1. 1.

At every vertex $j$ of the graph we define a scattering matrix
$\boldsymbol{\sigma}_{j}(k)$ associated with the boundary condition used at
the vertex $j$;

  2. 2.

a particle that propagates along the edge ($i,j$) (if present), contributes a
factor $z_{ij}=z_{ji}=e^{ik\ell_{ij}}$, where $\ell_{ij}$ is assumed to be the
metric distance between the nodes $i$ and $j$;

  3. 3.

at each existing edge, let it be ($i,j$), we introduce a factor $p_{ij}$,
representing the probability flowing between $i$ and $j$, and $p_{ji}$,
flowing in the opposite direction. The quantity $p_{ij}$ is interpreted as the
weight associated with the flow of the particle on that edge. For each node
$j$, and $\mathcal{N}_{j}^{n}$ as the neighborhood of the node $j$ minus the
node $n$. For the unitary equivalence between the scattering paths Green’s
function of eqn. (22) and the Green’s function based on the adjacency matrix
of eqn. (31), the functions $p_{ij}$ must satisfy the equation:

$p_{ij}=\sum_{l\in{\mathcal{N}_{j}^{n}}}z_{ij}\sigma_{j}^{[\\{j,l\\},\\{i,j\\}]}A_{jl}p_{jl}+\delta_{jn}z_{in}\sigma_{n}^{[e_{f},\\{i,n\\}]},$
(33)

and similarly for $p_{ji}$, provided that we simply swap the indices $ij$. At
each vertex $i$ we associate one $p_{ij}$ for every $j\in\mathcal{N}_{i}$. The
last term in (33) is, as said above, the transmission amplitude at the vertex
$n$ from the edge $\\{i,n\\}$ to the lead $e_{f}$. By construction, the graph
is simple, and the particle cannot hop on the node, thus $p_{ii}=0$ for all
nodes.

Given the definition above of the $p_{ij}$’s, eqn. (32) make sense: the total
transmission probability is simply equal to the total flow from the incoming
node, which must be equal to the flow to the exit node.

The remarkable property of the construction above is that infinite sum over
the scattering paths is replaced by a finite matrix inverse. This approach is
powerful, as constructing the Green’s function is equivalent to knowing
everything about our quantum system [109]. For instance, bound-state energies
can be obtained from the poles of the Green function, and the wave functions
from the associated residues [101]. In fact, $|T_{\Gamma}|^{2}$ represents the
global transmission probability from the lead $e_{i}$ to the lead $e_{f}$ and
it is constructed from the individual quantum amplitudes. A similar
construction was employed, before [127], in [118, 119].

An example at this point might be useful to clarify how to use this
construction. Consider for instance the graph in Fig. 5, and the associated
$p_{AB}$ and $p_{BA}$, as in Fig. 7. Because there is only one entry node of
degree $2$, then we have $T_{\mathcal{G}}=t_{A}p_{AB}$, from eqn. (32). The
equations for the $p$’s are then given by eqns. (33):

$\displaystyle p_{AB}$ $\displaystyle=$ $\displaystyle
z_{AB}(r_{B}p_{BA}+t_{B})$ (34) $\displaystyle p_{BA}$ $\displaystyle=$
$\displaystyle z_{BA}r_{A}p_{AB},$ (35)

from which we obtain

$p_{AB}=\frac{{z_{AB}}t_{B}}{1-z_{AB}z_{BA}r_{A}r_{B}},$

and finally, using $z_{AB}=z_{BA}=e^{ikl}$

$\displaystyle
T_{\Gamma}=t_{A}p_{AB}=\frac{t_{A}t_{B}{e^{ikl}}}{1-r_{A}r_{B}e^{{2}ikl}},$
(36)

which is the same expression as in eqn. (21), and that we had obtained earlier
from the sum over the infinite scattering path.

The construction presented so far is for general quantum graphs. But we need
to supply the definitions for the transmission and reflection probabilities,
which is where the Physics of the scattering is contained.

Clearly, we know from basic quantum mechanics that the $r$’s and the $t$’s
cannot be independent, and so the entries of the matrix
$\boldsymbol{\sigma}_{j}$ cannot either. For instance, in a one-dimensional
scattering, it is known that we must have $|r|^{2}+|t|^{2}=1$, following from
the boundary conditions of eqns. (14) and (II.3.1). In the case of a star
graph, if the wave is incoming from node $i$ to node $n$, the particle can be
either reflected, or it can be transmitted in multiple directions, then we
must have

$\displaystyle|\sigma_{n}^{{i,n},{n,i}}|^{2}+\sum_{s\in\mathcal{N}^{i}_{n}}|\sigma_{n}^{{i,n},{n,s}}|^{2}=1,$
(37)

or, written in terms of the elements of $\sigma^{\prime}s$

$\displaystyle|r_{n}^{i,n}|^{2}+\sum_{s\in\mathcal{N}^{i}_{n}}|t_{n}^{i,s}|^{2}=1.$
(38)

We can write one of these for each scattering process where the particle
incomes from one of the nodes and scatters into the other neighbors, at the
given node $n$, and these represent the constraints of our local scattering
parameters. The equations above clearly generalize the standard one-
dimensional scattering on a heterogenous structure, like the one of a filament
with multiple branches.

For naked graphs, in which there are no potentials on the edges, and the node
is in a certain discrete sense isotropic, a set of quantum parameters
conditions are given by

$\displaystyle r_{n}$ $\displaystyle=$ $\displaystyle\frac{2}{d_{j}}-1$ (39)
$\displaystyle t_{j}$ $\displaystyle=$ $\displaystyle\frac{2}{d_{j}},$ (40)

independently from the directionality. It is easy to see that these satisfy
the conditions above, and physically represent the particle diffusion
isotropically at the node in all possible directions, with a probability only
dependent on the local degree. For dressed graphs, e.g. a non-zero potential,
this has to be generalized.

Figure 7: The graph above is the same scattering problem as in Fig. 5, but the
quantities $p_{AB}$ and $p_{BA}$ contain the information the graph’s total
Green’s function. Their definition is in eqn. (33). The transmission and
reflection probabilities are the same as those of Fig. 5.

####  II.3.4 Dressed graphs

Insofar we have assumed that the propagation on the edge of the graph is free.
We wish now to introduce a potential on the edges, e.g. if we have a local
Hamiltonian of the form

$-\psi_{\\{i,j\\}}^{\prime\prime}(x)=k^{2}\psi_{\\{i,j\\}}(x)+V_{\\{i,j\\}}(x)\psi_{\\{i,j\\}}(x),$
(41)

with a constant potential $V_{\\{i,j\\}}=u_{ij}\equiv u$ along the edges and
across the graph, the formalism of the previous section still applies,
provided that we renormalize the edge length. This effective energy $u_{ij}$
represents an effective barrier in the propagation between sites. First, let
us briefly explain why a constant potential can act both as an approximation
and why it can be physically motivated. First, a constant approximation of the
potential is not uncommon in transport. In fact, the Kronig-Penney potential
[138, 139, 140] is a common approximation in one-dimensional transport on a
lattice. Since here our lattice is replaced by a generic graph, the assumption
that the one-dimensional potential along each is the analog of a Kronig-Penney
on a quantum graph, and can be seen merely as an approximation to estimate the
dependence of the transport on such potential. The second comment is that a
constant approximation can be solved exactly locally within the framework we
introduced, and thus it is not only physically motivated but also of
convenience. Nonetheless, we will then consider the general case first, and
reduce to the constant case after.

As we have discussed earlier, in an arbitrary graph we can always treat a
vertex $j$ with its edges as a star graph. A star graph on $n$ vertices,
$S_{n}$, is a graph where one central vertex has degree $n-1$ and all other
vertices have degree $1$. Consider a star graph as the one depicted in Fig. 3
and let $\Psi=\left(\psi_{\\{j,1\\}},\ldots,\psi_{\\{j,n\\}}\right)^{T}$. The
scattering solutions are given by

$\displaystyle\psi_{\\{i,j\\}}(x)={}$ $\displaystyle
e^{-ik_{ij}x}+r_{j}^{[\\{j,i\\},\\{i,j\\}]}e^{ik_{ij}x},$ (42)
$\displaystyle\psi_{\\{j,l\\}}(x)={}$
$\displaystyle\sqrt{\frac{k_{ij}}{k_{jl}}}\,t_{j}^{[\\{j,l\\},\\{i,j\\}]}e^{ik_{jl}x},$
(43)

where

$k_{ij}=k_{ji}=\sqrt{\frac{2m(E-u_{ij})}{\hbar^{2}}}=k\sqrt{1-\frac{u_{ij}}{E}},$
(44)

and $u_{ij}$ is the constant potential along the edge $\\{i,j\\}$. Note that
if $u_{ij}>E$, then $k_{ij}$ becomes complex and thus modes become evanescent.
Thus, $u_{ij}$ acts as an effective cutoff for the energy of the particles
that are transmitted on the graph.

Assuming the Neumann boundary condition

$\displaystyle\psi_{\\{j,l\\}}={}$ $\displaystyle\varphi_{j},\qquad\forall
l\in E_{j},$ (45) $\displaystyle\sum_{l\in E_{j}}\psi_{\\{j,l\\}}^{\prime}={}$
$\displaystyle 0,$ (46)

we obtain the transmission and reflection amplitudes

$\displaystyle r_{j}^{[\\{j,i\\},\\{i,j\\}]}={}$
$\displaystyle\frac{k_{ij}-\sum_{l\in E_{j}^{i}}k_{jl}}{K},$ $\displaystyle
t_{j}^{[\\{j,l\\},\\{i,j\\}]}={}$
$\displaystyle\frac{2\sqrt{k_{ij}k_{jl}}}{K},$ (47)

where $K=\sum_{l\in E_{j}}k_{jl}$, $E_{j}$ is the neighborhood of the vertex
$j$, and $E_{j}^{i}$ is the neighborhood of the vertex $j$ but with the vertex
$i$ excluded.

For a vertex with degree two

$r_{j}^{[\\{j,i\\},\\{i,j\\}]}=\frac{k_{ij}-k_{jl}}{k_{ij}+k_{jl}},$ (48)
$t_{j}^{[\\{j,l\\},\\{i,j\\}]}=\frac{2\sqrt{k_{ij}k_{jl}}}{k_{ij}+k_{jl}},$
(49)

which are the quantum amplitudes for a one-dimensional step potential.

The exact scattering Green’s function for a dressed quantum graph with
adjacency matrix $A(\Gamma)$ can be written as

$G_{\Gamma}=\frac{m}{i\hbar^{2}k}T_{\Gamma}e^{ik_{i}x_{i}+ik_{f}x_{f}},$ (50)

where $T_{\Gamma}=\sum_{j\in
E_{i}}\sigma_{i}^{[\\{i,j\\},e_{i}]}A_{ij}p_{ij}$, and

$p_{ij}=\sum_{l\in{E_{j}^{n}}}z_{ij}\sigma_{j}^{[\\{j,l\\},\\{i,j\\}]}A_{jl}p_{jl}+\delta_{jn}z_{in}\sigma_{n}^{[e_{f},\\{i,n\\}]},$
(51)

with

$z_{ij}=e^{ik_{ij}\ell_{ij}}=e^{ik\sqrt{(1-u_{ij}/E_{k})}\ell_{ij}}.$ (52)

Let us now assume that $u_{ij}=u$ and $l_{ij}=l$. Then, we can replace
$l\equiv l(k)=l\sqrt{1-\frac{u}{E_{k}}}$ on every edge of the graph.

If all the wave vectors are equal, i.e, $k_{ij}=k$, we have

$\displaystyle r_{j}$ $\displaystyle=$ $\displaystyle\frac{2}{d_{j}}-1,$ (53)
$\displaystyle t_{j}$ $\displaystyle=$ $\displaystyle\frac{2}{d_{j}},$ (54)

which are the well-known quantum amplitudes for the Neumann boundary condition
for “naked” quantum graphs we discussed earlier. As such, we can use the same
framework as before, provided that we renormalize the length at each
particular value of $k$, which explains why these types of graphs are
“dressed". Thus, in the dressed graph case, we will have an extra parameter to
analyze, which is $u$.

##  III Classical model of filament growth and effective quantum graph

At this point, we need to discuss how the filaments are formed in our model.
Below, we provide the simplest non-trivial model which incorporates the
voltage for the filament growth.

First, we consider a certain narrow atomic region, of size $W\times W\times
H$, where $H$ is the horizontal line dimension of the region, and $H$ is its
vertical dimension, in atomic units. Let us now assume that the top (anode)
and bottom (cathode) of the region $z=H$ and $z=1$, are connected to a certain
battery applying a voltage $V$. We assume then that the anode and cathode have
a certain voltage $v_{a}=\frac{V(t)}{2}$ and $v_{c}=-\frac{V(t)}{2}$, and with
an opposite sign if the battery’s voltage changes sign. These are the simplest
boundary conditions. We follow the following algorithm for every time $t$

  * •

We assume that at the anode and cathode, we have a certain time-dependent
voltage $V(t)$ applied, driven by the battery. Inside the region, the vector
field is given by $\vec{E}(x)=-\vec{\nabla}V(p)$, where $V(p)$ is the
potential inside the region. To calculate the potential $V(p)$, we solve for
the equation

$\displaystyle\triangle V_{t}(\vec{x})$ $\displaystyle=$ $\displaystyle 0$
(55) $\displaystyle V|_{\mathcal{B}_{a}(t)}$ $\displaystyle=$
$\displaystyle\frac{V(t)}{2}$ (56) $\displaystyle V|_{\mathcal{B}_{c}(t)}$
$\displaystyle=$ $\displaystyle-\frac{V(t)}{2}$ (57)

where ${\mathcal{B}_{a}}$ is the subset of points associated with the anode,
while ${\mathcal{B}_{c}}$ is associated with the cathode and
$\triangle=\vec{\nabla}\cdot\vec{\nabla}$. The equation above can be solved
numerically using an iterative method. Once we have obtained the solution
$V(\vec{x})$ everywhere, $\vec{E}(x)$ can be obtained.

  * •

Insofar, however, we have not discussed how the filament is formed
dynamically. We assume that if $V>0$, a certain density per unit of time and
unit of voltage $\rho_{0}$ is associated to the creation of particles of mass
$m$, formed at the center of the cathode, and if $V<0$ it is formed at the
center of the anode, with zero momentum. Then, we solve numerically for the
equation of motion for each particle $i$, which are assumed to be all
identical ions of charge $e_{c}$ abd $m_{p}$,

$\displaystyle\frac{d^{2}\vec{x}_{i}}{dt^{2}}=-\frac{e_{c}}{m_{p}}\vec{E}(\vec{x})+W(t).$
(58)

It is easy to see from eqn. (58) that we are using a molecular dynamics
approach to the construction of the filament between the two leads. If a
particle reaches the anode, or cathode, or touches another particle, a
filament is formed. If the filament touches the cathode, then the number of
particles is associated with $\mathcal{B}_{c}(t)$, and if at the anode they
are associated with $\mathcal{B}_{a}(t)$. Following this prescription, the
electric field $\vec{E}_{t}$ is also dynamic, but the time label is only to
capture the electric field at a particular time $t$ in the simulation. Because
the speed of the particles is assumed to be much lower than the speed of light
($V$ small enough), $p_{i}/m\ll c$, we are still solving for electrostatic
equations. The quantity $W(t)$ is a Wiener noise associated with the random
formation of the filament. Because of this, filament formation is a stochastic
process. Particles can also tick to the partly formed filament: if a certain
particle is directly in the vicinity of the filament, the particle can stick
and grow the filament. Thus, in principle, the filament can grow both in
horizontal and vertical directions.

  * •

If the electrodes are connected by a sequence of particles, we construct the
connected component of the effective graph $\mathcal{G}_{(}t)$, which encodes
the adjacency of the atoms. We construct the set of vertices $\mathcal{N}_{c}$
and $\mathcal{N}_{a}$ which are the vertices connected to the cathode and
anode respectively via a depth-first-search algorithm. These are the vertices
which associated with the Lippman-Schwinger states.

  * •

If the voltage flips sign, then the graph $\mathcal{G}(t)$ can decay. We
associate a certain probability $p_{d}$ per unit of time such that the
vertices (the atoms) can be removed. In parallel, the growth process can
restart from the other lead, and particles appear at random at the other end.
The filament formation then restarts. As such, the device we are simulating is
apolar.

An example of the filament formed via this process is shown in Fig. 8, with
the corresponding current and voltage cycle, the I-V diagram and the
transmission probability as a function of the energy for the final junction.

Figure 8: In the center images, we show a filament formed during an
aggregation process, due to the ion migration, induced by the electric field,
and equivalent to a quantum graph. We see the transmission spectrum on the top
right panel. The voltage is shown on the left panels. In the central figure,
the size of the atoms has been reduced to a point to show the connectivity
graph. For this figure, we have chosen an insulator with constant conductance.
At the bottom of the graph in the center panel, all the points on the surface
are connected to a single node as in Fig. 2. Five cycles of filament formation
(color-coded). In the bottom panels, we see the final cycle filament formed.
In the corresponding I-V curves, the color-coded hysteresis curves. The
corresponding currents are shown on the left of each I-V curve panel. The
difference between the left and right is the choice of the insulator. In the
case of the right, we have chosen a Poole-Frenkel insulator, with exponential
leakage. The Schottky barrier is approximately $V_{0}=0.3V$. This implies an
exponential envelope, usually observed in experiments [66]. In the left panel,
we have chosen an insulator constant with respect to the voltage. The
transmission probabilities as a function of the energy are shown in the bottom
right of each filament.

An important comment is that at this point, the device’s effective
conductivity can be calculated from the Landauer formula, using
$\sigma_{r}=\frac{I_{t}}{V_{t}}$, where $I_{t}$ is the effective. However, in
experimental realizations, it is hardly possible to isolate the filament
uniquely. As a result, the effective conductivity of the device is typically
the one of the filament plus a parallel conductive device, whether this is a
dielectric or an actual resistor. It is often argued in fact that even at a
few Kelvins, it is not possible to isolate the linear resistive switching
effect, but one has to consider other types of tunneling effects [66]. For
this reason, we add to the resistivity an effective conductivity, mimicking a
classical or semiclassical resistor. The effective current is then given by

$\displaystyle I(V)$ $\displaystyle=$ $\displaystyle I_{L}(V)+I_{R}(V)$
$\displaystyle=$ $\displaystyle
G_{0}\int_{-\infty}^{\infty}dE(f_{L}(E)/N-f_{R}(E)/M)\cdot$
$\displaystyle\hskip 56.9055pt\cdot(\frac{h}{2e^{2}}\sigma_{0}(V)+T(E)).$

For this paper, we use both a classical resistor and thus $\sigma_{0}$ can
either be a constant or a voltage-dependent quantity. In the case it is a
constant, then we will observe at low temperatures just a resistive switching,
without nonlinearity, but even if the device’s filament is not present, a
current. If $\sigma_{0}(V)$ is nonlinear, then the observed device’s current
can be nonlinear. In this paper we also consider a tunneling effect [142,
143], with $\sigma_{0}(V)\sim\sigma_{0}e^{V/V_{c}}$, with $V_{c}\approx 0.15$
volts, with $\frac{h}{2e^{2}}\sigma_{0}\sim 10^{3}\Omega^{-1}$. This is the
number of the same order of magnitude as the one typically observed in
tunneling experiments, and such that at voltages of the order $V\sim 1$ volts,
we obtain currents of the order of the microamperes.

Figure 9:  Current-Voltage characteristic within the context of the model
developed in this paper averaged over 100 Monte Carlo samples, and as a
function of the dressing parameter $u$, for $u=0$, $u=0.01$ and $u=0.05$ and
$x=0$.

##  IV Numerical simulations

We now discuss the results of our numerical experiments. The goal is to
provide some qualitative $I-V$ curves and the dependence on some key
underlying parameters.

###  IV.1 Parameters choice

In our simulations, the parameters have been chosen as follows. First, the
depth of the region is such that approximately $10$ atoms can be in a straight
line, a number compatible with the approximations we wished to make from the
start (e.g. a discrete rather than continuous filament). In particular, we
have chosen an arbitrary time unit fixed, however, by the following
prescription. In the molecular dynamics of eqn. (58) the parameters have been
chosen such that within $10-50$ time units, we observe a filament formation at
a certain target frequency $\omega^{*}$, assuming
$V(t)=V_{0}\cos(\omega^{*}t)$, with $V_{0}=0.35$ V. In terms of the voltage
drive, the time units can be expressed such that under the sinusoidal drive at
$V\approx 0.3$ we do observe a switching, which is a physically plausible
number for the switching.

Also, for simplicity, we have assumed that the Fermi chemical energy $\mu_{F}$
is much smaller than the voltage applied, $eV$. In doing so, we have
essentially constrained our parameters to be at the typical frequency at
which, depending on the ions that form the filament, the observation of the
$I-V$ hysteresis first shows a memristive effect. At this point, the free
parameters are then the strength of the stochastic force in eqn. (58) and the
potential $u$ in the quantum graph. The strength of the stochastic force has
been chosen to be $\langle W(t)^{2}\rangle_{t}=\sqrt{0.005}$, which was the
smallest value such that at every single simulation we could obtain a
different filament. This is dictated by the necessity that we wish to perform
a Monte Carlo analysis of the $I-V$ curves. The molecular dynamics simulations
were instead performed integrating the Newtonian equations of motion for the
particles in the electric field in the narrow region.

###  IV.2 Analysis of the results

Provided that our model of filament growth is limited in many ways, let us
first comment that we had first written a 2-dimensional code. In spirit, the
results are similar to those that we have obtained here, but some comments are
in order. In the three-dimensional case we discussed here, the filament
attracts particles from three dimensions, meaning that it is much more likely
for a shortcut to be found. This is why we found that there are many more ways
to have significant jumps in current in three dimensions rather than two.
Moreover, we have found that a good model for the growth of the filament in
two dimensions, as an alternative for the voltage-induced molecular dynamics,
was a diffusion-limited aggregation process [144], without having to solve for
the Laplace equation. These results will be discussed elsewhere, but from the
point of view of the present paper, it is experimentally more relevant to
study the three-dimensional case.

While this paper aims to provide a qualitative agreement between previously
obtained experimental results and theoretical ones, some comments are in
order. For that purpose, let us look at Fig. 8, where several voltage cycles
were applied to the device with $u=0$. We see the filament formed on the
negative voltage side at the center.111The video of the filament formation for
Panel 8 is shown at https://youtu.be/mhgoE0MWGk0.

We have found that in the vast majority of our numerical experiments, only one
filament was formed. The reason is that as the filament is formed, it attracts
any free particle in the vicinity as a result of eqn. (57). This implies that
although the filament has a finite thickness, the charge propagation is quasi-
one-dimensional. However, gaps in the filaments can be filled by floating
particles, leading to discrete jumps in the resistivity. This is the jump we
see for instance in Fig. 8 (bottom left IV). The (middle right) is the voltage
profile applied to the device. The vertical jumps, both on the positive and
negative side, due to the resistive switching, e.g. the filament reaching the
cathode/anode. The curves in Fig. 8 represent several realizations of the
filament. The averaged curves are shown in Fig. 9 for various sizes of the
junction and values of the parameter $u$.

We have also performed a slightly more refined statistical analysis of the
conductance distributions. In particular, in Fig. 9 we plot the averaged
hysteresis curve over 100 realizations of the filament, and for three values
of $u=0,0.01,0.05$ and junction size of $5,10,15,20$ atoms. The shadow region
represents the error across various realizations of the filament. We see that
for larger values of the gap, the memory effect is more noticeable, due to the
intuitive fact that there are more possibilities for the filament to be
formed. In Fig. 8 (bottom) we show different cycles of hysteresis for two
types of insulators.222A video for the Frenkel-Poole insulator is
https://youtu.be/uUrPZnSrIpQ . For the insulator model which is constant in
the voltage, a video is available at https://youtu.be/TfrhEuasAD4 .

We then note that, for $u=0.01$, the switching to the ON state is to a lower
resistive state. In particular, because of this, the tunneling phenomenon is
such that the switching to the ON is almost hidden in the nonlinearity caused
by the tunneling. This feature was also observed in the measurements of
niobium oxides [66]. Thus, this qualitative behavior of the pinched hysteresis
obtained by this model is realistic. The distribution of the conductance is
plotted in Fig. 10. Overall, this provides a hint about the relationship
between the size of the gap and the conductance (see for instance [67], Fig.
19), expected from an Anderson model with finite localization length [68]. An
important point to make is that the parameter $u$, in particular for larger
values, can be interpreted in terms of an (effective) Schottky barrier to
electron propagation.

Figure 10: Conductance distribution as a function of the gap size (in atoms,
Dim) and the parameter $u$. For smaller regions, we see that the probability
of having a conductance near one is non-zero. For larger gaps, the probability
decreases and peaks near smaller values. The parameter $u$ broadens the
distribution. Each figure is averaged over 100 realizations of the filament.

##  V Conclusions

The present paper introduced a framework to study quantum transport within the
context of filamentary switching. The method we developed is not based on the
quantum scattering in a continuum or a lattice approximation but assumes that
the filament is finite-dimensional and discrete, e.g. charge carriers can move
on a discrete graph. This paper is based on the Landauer-Büttiker formula
applied to quantum graphs, which provides the theory to study one body
scattering on discrete structures such as metric graphs. In addition, we have
used a model for the filament growth which is based on a molecular dynamics
simulation. The molecular dynamics we used is the simplest model that supports
filament growth and is based on the solution of the Laplace equation and the
derivation of the electrical vector field in the narrow region.

Yet, although our results are not compatible with previous numerical studies
for a continuous saddle for silver ions [70], and the experimental results for
niobium pentoxide [66], in which filamentary switching occurs, they are more
similar to the junction model by Milano and collaborators [50, 51, 27], for
which stochastic models are currently being developed [52, 53]. In addition,
the distribution of conductances in Fig. 10 qualitatively is consistent with
the conductance distributions arising from random matrix theory [67].

As far as we understand this is the first time that a non-many body quantum
approach has been used for a realistic dynamical filament formation and where
we model the discrete nature of the filament, rather than it being assumed
continuous. In this sense, we believe our paper provides a new methodology to
analyze quantum transport on these structures. The advantage is in the
scalability of our numerical codes: our codes scale as $\mathcal{O}(E^{3})$
(which comes from the matrix inversion for evaluating $p_{ij}$), where $E$ is
the number of active channels, or edges of the graph (between atoms/molecules)
in the narrow region. For the case with ballistic motion, the code scales as
$\mathcal{O}(d_{max}NE^{3})$, but still subexponentially, where $N$ is the
number of atoms and $d_{max}$ the maximum graph degree. As we discussed in the
text, this extra factor comes from the construction of the scattering
amplitudes $\sigma_{j}$’s. This paper did not consider a particular type of
material, but we focused on generic features of the I-V curves which are
phenomenological in nature. While this is just a first step in this direction,
various improvements can be made. First of all, a more precise model of the
filament formation can be introduced, which carefully models the evolution of
the filament over time. Moreover, it allows the introduction of a variety of
modeling assumptions about the scattering and channel properties of the
filament and atoms. One important simplifying assumption used here was that
all the graph edges were of equal length $l$, the unit of distance between
atoms. Relaxing such assumption allows also to have localization phenomena and
drastic reductions of the conductance (induced by small $T_{fi}$ between the
leads) due to chaotic scattering in the graph [145].

One drawback of our methodology is that the evaluation of the scattering
amplitudes scales with the number of edges of the graph, the reason for which
we needed to focus on a few atoms. However the calculation of the scattering
amplitude can be scaled, and study larger graphs and larger regions. In fact,
the technique scales most quadratically in the number of nodes or atoms, and
while challenging, it is feasible using a proper computing architecture. Also,
because the method requires the inversion of a large matrix, it is also
possible to use GPUs for this purpose.

Nonetheless, some of the findings of this paper are particularly interesting.
For instance, the non-monotonicity of the ON state on the dressing parameter
is an interesting feature that can be studied experimentally. Moreover, this
approach can be extended also beyond the one-body Hamiltonian, which can be
generalized to a local one-dimensional quantum field [146, 147, 148, 149].

These and other generalizations and improvements will be considered in future
papers.

###### Acknowledgements.

This work was carried out under the auspices of the NNSA of the U.S. DoE at
LANL under Contract No. DE-AC52-06NA25396. FC was also financed via an LDRD
IMS-Rapid Response grant, and this research was supported in part by grant NSF
PHY-1748958 to the Kavli Institute for Theoretical Physics (KITP). The work of
F.M.A. was partially supported by Coordenação de Aperfeiçoamento de Pessoal de
Nível Superior (CAPES, Finance Code 001). It was also supported by Conselho
Nacional de Desenvolvimento Científico e Tecnológico (CNPq), Instituto
Nacional de Ciência e Tecnologia de Informação Quântica (INCT-IQ). FMA
acknowledges financial support by CNPq Grant 313124/2023-0. The authors are
indebted to András Halbritter for several comments on the manuscript.
Author contributions. F. Caravelli conceived the study. F. Andrade, F.
Caravelli and A. A. Silva contributed to the theoretical analysis in this
order. F. Caravelli and A. A. Silva wrote the code for the analysis. F.
Caravelli wrote the first draft. All authors contributed to the writing of the
draft.
Data availability. A Mathematica code has been assigned a FigShare identifier
[150].

## References

  * [1] L. Chua, IEEE Trans. on Circ. Theor. 18, 507–519 (1971)
  * [2] D. B. Strukov, G. S. Snider, D. R. Stewart, R. S. Williams, Nature , 453, 80–83 (2008)
  * [3] D. S. Jeong et al, Adv. Electron. Mater., 2: 1600090. (2016)
  * [4] T. Serrano-Gotarredona et al., Front. Neurosci. 7: 2 (2013) https://doi.org/10.3389/fnins.2013.00002
  * [5] F. Caravelli, J. P. Carbajal, Technologies 6(4), 118 (2018)
  * [6] A. Chung, J. Deen, J.-S. Lee, M. Meyyappan, Nanotechnology 2010, 21, 412001.
  * [7] J. J. Yang, D. B. Strukov, D. R. Stewart, Nat. Nano. 8, 13–24 (2013)
  * [8] M. A. Zidan, J. P. Strachan, W. D. Lu,Nat. Electron. 1, 22–29 (2018)
  * [9] Q. Xia, J. J. Yang, Nature Mat. , 18, 309–323 (2019)
  * [10] C. Li, et al. Nat. Elect. 1, 52–59 (2018)
  * [11] C. Li, et al. Nature Mach. Int. 1, 49–57 ( 2019)
  * [12] C. Li, et al. Nature Comm. 9, 2385 (2018)
  * [13] M. Hu, et. al.Advanced Materials 30, 1705914 (2018)
  * [14] M. Prezioso, et al. Nature 521, 61–64 (2014)
  * [15] H. Jiang, et al. Nature Electronics 1, 548–554 (2018)
  * [16] A. Sawa, Mater. Today 11, 28–36 (2008)
  * [17] J. J. T. Wagenaar, et al.J. Appl. Phys. 111, 014302 (2012)
  * [18] J. M. van Ruitenbeek, M. Morales-Masis, E. Miranda, In Resistive Switching: From Fundamentals of Nanoionic Redox Processes to Memristive Device Applications Ch. 7\. Quantum Point Contact Conduction; Ielmini, D., Waser, R., Eds.; Wiley-VCH, Weinheim, pp 197– 224\. (2016)
  * [19] M. D. Pickett, R. S. Williams, Nanotechnology 23, 215202 (2012)
  * [20] H. Mähne, et al. Semicond. Sci. Technol. 29, 104002 (2014)
  * [21] S. Kim, et al. IEEE VLSI Symp. Tech. Dig. 155–156 (2012)
  * [22] F. Caravelli, A. Hamma, M. Di Ventra, EPL, 109, 2 (2015)
  * [23] A. Zegarac, F. Caravelli, EPL Perspectives 125 10001 (2019)
  * [24] F. Caravelli, F. L. Traversa, M. Di Ventra, Phys. Rev. E 95, 2 (2017)
  * [25] F. Caravelli, F. C. Sheldon, F. L. Traversa, Sci. Adv. 7, eabh1542 (2021)
  * [26] J. Hochstetter, R. Zhu, A. Loeffler et al. Nat Commun 12, 4008 (2021).
  * [27] F. Caravelli, G. Milano, C. Ricciardi, Z. Kuncic, Ann. Phys. 535, 8 (2023)
  * [28] M. Aono, T. Hasegawa, IEEE Proc. 98, 2228–2236 (2010)
  * [29] B. Cheng, et al Comm. Physics 2, 28 (2019)
  * [30] N. Agrait, A. L. Yeyati, J. M. van Ruitenbeek, Physics Reports , 377, 81–279 (2003)
  * [31] P. Makk, Sz. Csonka, A.Halbritter, Phys. Rev. B 78, 045414 (2008)
  * [32] E. Scheer,Nature 394, 154–157 (1998)
  * [33] S. Gamage et al, Perovskite Nickelate Memristive Devices ,ACS Nano 18, 2105-2116 (2024)
  * [34] R. Kubo, J. Phys. Soc. Jpn. 12 570 (1957)
  * [35] R. Landauer, IBM J. Res. Dev. 1 223 (1957)
  * [36] D.S. Fisher, P.A. Lee, Phys. Rev. B 23 6851 (1981)
  * [37] E.N. Economou, C.M. Soukoulis, Phys. Rev. Lett. 46 618 (1981)
  * [38] Y. Imry in Directions in Condensed Matter Physics, (Ed. by G. Grinstein and G. Mazenko) World Scientific, Singapore (1986)
  * [39] M. Büttiker, et alPhys. Rev. B 31 6207 (1985)
  * [40] M. Büttiker, Phys. Rev. Lett. 57 1761 (1986)
  * [41] M. Büttiker, IBM J. Res. Dev. 32 63 (1988)
  * [42] M. Büttiker, IBM J. Res. Dev. 32 317 (1988)
  * [43] M. Büttiker, Phys. Rev. B 38 9375 (1988)
  * [44] A.F. Andreev, ZhETF 46 1823 (1964)
  * [45] A.F. Andreev, ZhETF 47 2222 (1964)
  * [46] A.F. Andreev, ZhETF 49 655 (1965)
  * [47] G.E. Blonder, M. Tinkham, T.M. Klapwijk, Phys. Rev. B 25 4515 (1982),
  * [48] B. Ludoph, et al.Phys. Rev. B , 61, 8561 (2000)
  * [49] J. Halbritter, Appl Phys A: Materials Science & Processing , 43, 1–28 (1987)
  * [50] G. Milano, G. Pedretti, M. Fretto, L. Boarino, F. Benfenati, D. Ielmini, I. Valov, C. Ricciardi, and Rewiring in Multi-Terminal Self-Organizing Memristive Nanowire Networks, Adv. Intell. Syst. 2, 2000096 (2020)
  * [51] G. Milano et al,Michael Kozicki, Sayani Majumdar, Mariela Menghini, Enrique Miranda, Carlo Ricciardi, Stefan Tappertzhofen, Kazuya Terabe, and Ilia Valov*, Quantum Conductance in Memristive Devices: Fundamentals, Developments, and Applications, Adv. Mater. 34, 2201248 (2022)
  * [52] Y. Slipko, Y. V. Pershin, Phys. Rev. E 107, 064117 (2023)
  * [53] G. Milano, F. Michieletti, E. Miranda, C. Ricciardi, ResearchSquare preprint (2024)
  * [54] D. Molnár et al.Nanoscale 10, 19290–19296 (2018)
  * [55] X. Liu, et al. IEEE Electron Device Lett. 33, 236–238 (2012)
  * [56] H. Wylezich et al,J. Vac. Sci. Technol. B 33, 01A105 (2015)
  * [57] H. Wylezich, et al., ACS Appl. Mater. Interfaces 6, 17474–17480 (2014)
  * [58] A. Marchenkov, et al.,Phys. Rev. Lett., 98, 046802 ( 2007)
  * [59] W. Schulz, et. al. Phys. Rev. B , 46, 14001–14004 (1992)
  * [60] N. Moussy, H. Courtois, B. Pannetier, Eur. Phys. Lett., 55, 861–867 (2001)
  * [61] M. Vinet, C. Chapelier, F. Leoch, Phys. Rev. B , 63, 165420 (2000)
  * [62] L. Y. L. Shen, AIP Conference Proceedings , 4, 31–44 (1972)
  * [63] H. le Sueur, et al.Lett. 100, 197002 (2008)
  * [64] E. Scheer, et al.Phys. Rev.Lett. 86, 284–287 (2001)
  * [65] W. Belzig, C. Bruder, G. Sch¨on, Phys. Rev. B 54, 9443–9448. (1996)
  * [66] T. N. Török, M. Csontos, P. Makk, A. Halbritter, Nano Letters 20 2 (2020) 1192-1200
  * [67] C. W. J. Beenakker, Random-matrix theory of quantum transport, Rev. Mod. Phys. 69, 731 (1997)
  * [68] J.-L. Pichard, in Quantum Coherence in Mesoscopic Systems, edited by B. Kramer, NATO ASI Series B254 (Plenum, New York), p. 369 (1991)
  * [69] J. C. Cuevas, A. L. Yeyati, A. Martin-Rodero, Phys. Rev. Lett. , 80, 1066–1069 (1998)
  * [70] A.Gubicza, et al. Sci. Rep. 6, 30775 (2016)
  * [71] A. Gubicza, et al.Nanoscale 7, 11248–11254 (2015)
  * [72] J. C. Cuevas, E. Scheer, Molecular Electronics: An Introduction to Theory and Experiment; World Scientific, Singapore, 2010.
  * [73] S. Datta, Electronic transport in mesoscopic systems (Cambridge University Press, 1995).
  * [74] M. Di Ventra, Electrical transport in nanoscale system, (Cambridge University Press, 2008)
  * [75] R. Cron, et al.Contacts. Phys. Rev. Lett. , 86, 4104–4107 (2001)
  * [76] C. Schirm, et al.Nanotechnology , 8, 645–648 (2013)
  * [77] L. Olesen, Phys. Rev. Lett. 72 2251 (1994)
  * [78] J.M. Krans, Nature 375 767 (1995)
  * [79] E. Scheer, P. Joyez, D. Esteve, C. Urbina, and M.H. Devoret, Phys. Rev. Lett. 78 3535 (1997)
  * [80] V. Rodrigues, T. Fuhrer, D. Ugarte, Phys. Rev. Lett. 85 4124 (2000)
  * [81] S. Frank, P. Poncharal, Z.L. Wang, W.A. de Heer, Science 280 1744 (1998)
  * [82] R. Martel, et al. Appl. Phys. Lett. 73 2447 (1998)
  * [83] M. Titov, C.W.J. Beenakker, Phys. Rev. B 74 041401(R) (2006)
  * [84] J.C. Cuevas and A.L. Yeyati, Phys. Rev. B 74 180501 (2006)
  * [85] J. Gonzalez, E. Perfetto, J. Phys.: Condens. Matter 20 145218 (2008)
  * [86] N.M.R. Peres, A.H. Castro Neto, F. Guinea, Phys. Rev. B 73 195411 (2006); Phys. Rev. B 73 239902(E) (2006)
  * [87] N. Tombros, et al.Nature Phys. 7 697 (2011)
  * [88] N. Gong, et al. Nat. Comm. 9, 2102 (2018)
  * [89] A. C. Torrezan, J. W. Strachan, G. Medeiros-Ribeiro, R. S. Williams, Nanotechnology 22, 485203 (2011)
  * [90] V. Kostrykin, R. Schrader J. Phys. A: Math. Gen. 32 595-630 (1999)
  * [91] G. Giusiano,et al.Int. J. Mod. Phys. B 18, 691 (2004)
  * [92] L. Pauling, J. Chem. Phys. 4, 673 (1936).
  * [93] K. Ruedenberg, C. W. Scherr, J. Chem. Phys. 21, 1565 (1953).
  * [94] E. W. Montroll,J. Math. Phys. 11, 635-648 (1970)
  * [95] R. G. J. Mills, E. W. Montroll , J. Math. Phys. 11, 2525-2538 (1970)
  * [96] J.-P. Roth, C. R. Acad. Sci. Paris Ser. I Math. 296, 793 (1983).
  * [97] T. Kottos and U. Smilansky, Phys. Rev. Lett. 79, 4794 (1997).
  * [98] T. Kottos and U. Smilansky, Ann. Phys. (NY) 274, 76 (1999).
  * [99] H. Schanz and U. Smilansky, Phys. Rev. Lett. 84, 1427 (2000).
  * [100] A. G. M. Schmidt, B. K. Cheng, and M. G. E. da Luz, J. Phys. A 36, L545 (2003).
  * [101] F. M. Andrade, A. G. M. Schmidt, E. Vicentini, B. K. Cheng, and M. G. E. da Luz, Phys. Rep. 647, 1 (2016).
  * [102] G. Berkolaiko and P. Kuchment, Introduction to Quantum Graphs (American Mathematical Society, 2012).
  * [103] G. Berkolaiko, arXiv (2016), 1603.07356.
  * [104] J. Bolte and S. Endres, Ann. Henri Poincaré 10, 189 (2009)
  * [105] V. Kostrykin, J. Potthoff, and R. Schrader, Contemp. Math. 447, 175 (2007).
  * [106] V. Kostrykin and R. Schrader, J. Phys. A 32, 595 (1999).
  * [107] K. Chadan and P. C. Sabatier, Inverse Problems in Quantum Scattering Theory, 2nd ed. (Springer, 1989).
  * [108] V. Kostrykin and R. Schrader, Contemp. Math. 415, 201 (2006).
  * [109] E. N. Economou, Green’s Functions In Quantum Physics, 3rd ed. (Springer-Verlag, 2006).
  * [110] R. Diestel, Graph Theory, 4th ed., Graduate Texts in Mathemat- ics Vol. 173 (Springer, 2010).
  * [111] P. Kuchment, Analysis on Graphs and Its Applications, Proc. Symp. Pure. Math. 77, 291 (2008)
  * [112] R. Landauer, Phys. D38 226 (1989)
  * [113] A. Kawabata, J. Phys. Soc. Japan 58 372 (1989)
  * [114] R. Landauer and T. Martin, Phys. B 175 167 (1991)
  * [115] R. Landauer, T. Martin, Phys. B 182 288 (1992)
  * [116] M. Büttiker, Phys. Rev. B 46 12485 (1992)
  * [117] D. Waltner, U. Smilansky, Acta. Phys. Pol. A 124, 1087 (2013)
  * [118] E. Ragoucy, J. Phys. A 42, 295205 (2009)
  * [119] V. Caudrelier, E. Ragoucy, Nucl. Phys. B 828, 515 (2010).
  * [120] S. Gnutzmann, U. Smilansky, Adv. Phys. 55, 527 (2006).
  * [121] F. M. Andrade, M. G. E. da Luz, Phys. Rev. A 84, 042343 (2011).
  * [122] R. Blümel, Y. Dabaghian, and R. V. Jensen, Phys. Rev. Lett. 88, 044101 (2002).
  * [123] R. Blümell, Y. Dabaghian, and R. V. Jensen, Phys. Rev. E 65, 046222 (2002).
  * [124] M. G. E. da Luz, B. K. Cheng, and M. W. Beims, J. Phys. A 34, 5041 (2001).
  * [125] F. M. Andrade, B. K. Cheng, M. W. Beims, and M. G. E. da Luz, J. Phys. A 36, 227 (2002).
  * [126] F. M. Andrade, Phys. Lett. A 378, 1461 (2014).
  * [127] F. M. Andrade, S. Severini, Phys. Rev. A 98, 062107 (2018).
  * [128] F. M. Andrade, A. G. M. Schmidt, E.Vicentini, B.K.Cheng, M.G.E.da Luz, Phys. Rep. 647, 1 (2016).
  * [129] M. Büttiker, Phys. Rev. B 41 7906 (1990
  * [130] C.W.J. Beenakker, Phys. Rev. Lett. 67 3836 (1991)
  * [131] G.B. Lesovik, Pis’ma v ZhETF 49 513 (1989)
  * [132] M. Büttiker, Phys. Rev. Lett. 65 2901 (1990)
  * [133] B.A. Lippmann and J. Schwinger, Phys. Rev. 79, 469 (1950)
  * [134] L.D. Landau and E.M. Lifshitz Quantum Mechanics: Non-Relativistic Theory (Oxford: Pergamon Press, 1977).
  * [135] M. A. Novotny, Phys. Rev. B 90 (16), 165103 (2014)
  * [136] G. Inkoom and M. A. Novotny, J. Phys. Communications 2(11), 115018 (2018)
  * [137] M. A. Novotny, G. Inkoom, T. Novotný, Eur. Phys. Lett. 143(2), 26005 (2023)
  * [138] Ashcroft, N. W.; Mermin, N. D., Solid State Physics; Saunders College, Philadelphia, 1976.
  * [139] C. Kittel, Quantum theory of solids (New York, Wiley 1987).
  * [140] R. de L. Kronig,W. G. Penney, Proc. of the Roy. Soc. A130 (814): 499–513.
  * [141] Baradaran, M. ; Exner, P., arXiv:2106.16019
  * [142] J. Frenkel, Phys. Rev. 54 (8): 647–648 (1938).
  * [143] J. G. Simmons, J. Appl. Phys. 34, 1793–1803 (1963).
  * [144] L. Niemeyer, L. Pietronero, H. Wiesmann, Phys. Rev. Lett. 52 1033-1036 (1984)
  * [145] S. Gnutzmann, U. Smilansky, Adv. in Phys. 55, 527 (2006)
  * [146] B. Bellazzini, M. Burrello, M. Mintchev, P. Sorba, Proc.Symp.Pure Math.77:639,2008
  * [147] B. Bellazzini, M. Mintchev, P. Sorba, J.Math.Phys.51:032302 (2010)
  * [148] B. Bellazzini, M. Mintchev, J.Phys.A39:11101-11118 (2006)
  * [149] V. Caudrelier, M. Mintchev, E. Ragoucy, J. Math. Phys. 55 (2014) 083524
  * [150] A.A. Silva, A. Halbritter, F. Andrade, F.Caravelli, Mathematica code for simulation of filament formation. FigShare DOI:10.6084/m9.figshare.25487914.v1 Link