# A Transfer Hamiltonian Model for Devices Based on Quantum Dot Arrays

Authors: S. Illera; J. D. Prades; A. Cirera; A. Cornet
Journal: The Scientific World Journal
 (2015)

---

## Abstract

We present a model of electron transport through a random distribution of interacting quantum dots embedded in a dielectric matrix to simulate realistic devices. The method underlying the model depends only on fundamental parameters of the system and it is based on the Transfer Hamiltonian approach. A set of noncoherent rate equations can be written and the interaction between the quantum dots and between the quantum dots and the electrodes is introduced by transition rates and capacitive couplings. A realistic modelization of the capacitive couplings, the transmission coefficients, the electron/hole tunneling currents, and the density of states of each quantum dot have been taken into account. The effects of the local potential are computed within the self-consistent field regime. While the description of the theoretical framework is kept as general as possible, two specific prototypical devices, an arbitrary array of quantum dots embedded in a matrix insulator and a transistor device based on quantum dots, are used to illustrate the kind of unique insight that numerical simulations based on the theory are able to provide.


---

## Body

## 1. Introduction

The demand for increasing the integrated density devices has led to the emergency of a whole generation devices based on confined structures. The MOS (metal-oxide-semiconductor) transistor is the archetype of a confined two-dimensional system [1]. Nevertheless, the possibility to enhance this confinement by embedding low-dimensional structures in an insulating matrix has opened new way for further downscaling. Compared to the standard bulk technology, the corresponding devices based on these structures have increased the structural and conceptual complexity. These structures (quantum dots, wires, or layers) can be used in single-electron devices [2], new memory concepts [3], and photo- or electroluminescent devices [4]. Concerning single-electron devices, they are currently conceived to take advantage of tunnel current between quantum states belonging to nanoscale particles [5, 6]. The single-electron devices based on quantum dots (Qds) appear to be potential candidates to improve, to complete, or even to replace the current MOS technology with which they may remain compatible. In order to be able to asses the potentials and capabilities of the various novel devices, a realistic theoretical estimation of the specific device performance is thus highly desirable. Within this context, the simulations of such devices must be performed not only to understand but also to predict experimental behaviors. Moreover, from a physical point of view, we will learn a lot from these simulations if they are independent of high-level experimental parameters (as tunneling rates, defective interfaces,…, etc.) and are based on low-level concrete ones (geometrical data, barrier height…).Concerning Qds, they are particularly attractive because they possess discrete energy levels and quantum properties similar to natural atoms or molecules due to the strong confinement in all three directions. This fact affects dramatically the electronic transport properties. Until now, research has mostly concentrated on single Qds and many novel transport phenomena have been discovered, such as the staircaselike current-voltage (I-V) characteristic [7], Coulomb blockade oscillation [8], negative differential capacitance [9], and the Kondo effect [10].From experimental point of view, rapid progress in microfabrication technology has made possible coupling Qds system with aligned levels [11–13]. In fact, the use of the Coulomb blockade phenomenon in systems made up of combinations of tunnel junctions and semiconductor Qds seems to offer promising perspectives, in particular in non-volatile memory applications and also for single-electron transistor [14]. Moreover, the concept of multi-dot memory using semiconductor nanocrystals embedded in an insulator matrix as floating-gate has already been demonstrated [15] and the quantization effects have been used in self-aligned double-stacked memory to improve the retention time [16].From a theoretical perspective, researchers have recently paid much attention to electron transport through several Qds, since multiple Qd provides more Feynman paths for the electron transmission [17]. The complexity of structure and physical mechanism and the prominent role of dimensional and quantum effects characterizing the operation of these novel Qds devices preclude the use of standard macroscopic bulk semiconductor transport theory. Many authors have studied the electron transport using NEGFF (nonequilibrium Green functions formalism) [18, 19], taking into account the potential due to the self-charge. However, up to now nobody has done a computation of transport in an extended arbitrary array of Qds using this framework since this approach is usually unfeasible to implement for large systems. On the other hand, rate equation type models used for lasers or light-emitting diodes often offer a satisfactory description of the charge transport. Moreover, this approach presents a more transparent vision of the electron transport. Thus, this model is easier to tinker with, in order to deal with more complicated nanostructures based on Qds.In this work, we present in full a model based on noncoherent rate equations [20], which is suitable to study the electron transport in Qd arrays. In a previous work [21], a preliminary version of this methodology was presented and used to obtain analytical solutions for electron transport in simple cases. The methodology was also compared with nonequilibrium Green’s function calculations, obtaining favorable results [22]. Despite the simplicity of the model, it provided good results and it was also easily scalable.Now, a complete model to simulate devices based on Qd arrays is presented. This model creates a compact modeling tool to study and simulate the electronic transport as a function of geometrical and basic material parameters, making possible to use it as a device design tool. The theoretical formalism and the assumptions made in the model are thoroughly described. The interaction between Qds and between these and the leads has been introduced by transition rates and capacitive couplings. The local potential effects were computed self-consistently. Inelastic and backscattering effects were neglected. Concerning the transition rates, the use ofab initio calculations is shown to be the best way to fully understand the underlaying tunneling physics in nanostructures. However, if first-principles calculations for single tunnel events were implemented, the huge effort required would make the simulation time increase in an unacceptable manner. This impractical computational time forced us to write a compact model with some assumptions and relax the expectation of accuracy when treating with few-electron devices operating through quantum features. Specifically, we used the one-dimensional WKB approximation, which neglects spatial variation of the wavefunction over the nontransport directions to describe the tunnel processes. The hole transport was also introduced obtaining new current terms and realistic expressions for the capacity in bipolar conduction. All of these have been implemented in a computational code conforming a powerful transport simulation tool that allows reproducing, explaining, and predicting the behavior of multiple devices based on Qds. Furthermore, we also present details about the computational implementation. Finally, in order to show the capabilities of the presented methodology, two examples of practical implementations of Qd-based devices were simulated: one single Qd and a multilayer structure that conform the basic building block of future devices based on Qds.


## 2. The Model

As in any device simulation, the ultimate goal of the presented approach is the prediction of the device response of one specific architecture (geometry, material,…, etc.) to a given variation in the external conditions (bias voltage) via the solution of the dynamical equations. First of all, we are going to describe the device architecture and, after that, we write the underlying equations.


### 2.1. The Structure

Figure1(a) shows the basic building blocks of our device which is in essence an insulator layer sandwiched between two metallic electrodes. Inside the insulator layer, a random distribution of N Qds can be inserted. This is the classical structure that is obtained due to the fabrication processes, a superlattice of insulator-semiconductor bilayers. Although a single connected Qd to the leads has been obtained creating a so-called single electron transistor (SET) [29], the research mainstream is focused on the properties of structures with many Qds to create nonvolatile memories [30], light-emitting devices [31], or solar cells devices [32].Figure 1
(a) The basic building block that forms the devices based on Qds: an array of Qds (they could be ordered or randomly dispersed) embedded in an insulator matrix sandwiched between two electrodes (leads). (b) Energy band scheme of the structure. The electron starts in the left lead and crosses through the dielectric matrix by tunneling processes assisted by the intermediate Qds. (c) Scheme of the equivalent electrical circuit of the device presented in (a). The Qds (circles) are connected between them and the leads (color blocks) by tunneling junctions. The tunneling junction is described as a capacitor in parallel with a current source.
(a)
(b)
(c)Due to the fabrication processes, the insulator thickness is large enough to avoid direct tunneling between the electrodes (or leads). Therefore, the electron current needs to pass through the Qds. In Figure1(b) the energy band diagram of the system is shown. The Qds are presented as potential wells in the dielectric energy band and they assist the tunneling processes. Thus, a correct description of the tunneling processes among the Qds and Qds-electrodes is needed. These tunneling processes can be described by tunnel junctions. From the electrical point of view, a single tunnel junction is described by a capacitance and a current source that depends on the voltage drop in the junction. The equivalent electrical circuit is shown in Figure 1(c) for an arbitrary Qd array.From the electrical scheme, two coupling equations govern the response of the system: on one hand we can write the charge conservation in each Qd. In the steady state, the charge conservation for this system is analogous to the electronic Kirchhoff’s current law. Thus, using the decomposition of the tunneling junctions in current sources, we can write(1)
0
=
∑
j
‍
I
i
j
i
∈
N
,for the ith Qd (equivalent to a node), where N is the number of Qds and the summation takes into account all the linked elements to this Qd. On the other hand, the current through a junction will depend on the voltage drop in the tunneling junction. Therefore, a second equation is necessary to describe the voltage in each Qd, the same as the Kirchhoff’s voltage law. First of all, we are going to show analytical expressions for the tunnel currents through a junctions and, after that, we will calculate the potential in each Qd using the Poisson equation.


### 2.2. Current through a Tunneling Junction

A usual method employed to describe the tunneling processes in devices is the tunneling Hamiltonian approach (also called Transfer Hamiltonian approach). This theory, thoroughly studied by many authors [33–35], treats the tunnelling events as a perturbation. The matrix coefficient |
T
L
R
|
2 quantifies the probability for a particle to transfer from a state of the left side of the barrier to a state of the right side by a tunnel process.The tunneling rate is determined using time-dependent perturbation theory by considering the electron from one side of the barrier as initial state and the electron on the other side as a final state. The tunneling rate from the left to the right states (both are considered as a part of a continuum) can be calculated using the Fermi’s golden rule [36]:(2)
d
2
W
k
→
L
→
k
→
R
=
2
π
ħ
T
L
R
2
ρ
R
(
E
R
)
ρ
L
(
E
L
)
δ
(
E
R
-
E
L
)
d
E
R
d
E
L
,where ρ
L and ρ
R are the density of states of the left and right side. From this expression, we can see that we only consider ballistic transport. This means that the electron does not suffer energy loss scattering processes when it moves through the barrier. Introducing the energy distribution function in each part of the barrier we can evaluate the total tunneling rate from all occupied states on the left to all unoccupied states on the right [37]: (3)
Γ
L
→
R
=
∫
-
∞
+
∞
‍
ρ
L
E
L
f
L
E
L
·
∫
-
∞
+
∞
‍
2
π
ħ
T
L
R
2
1
-
f
R
E
R
h
h
l
h
∫
-
∞
+
∞
‍
2
π
ħ
·
ρ
R
(
E
R
)
δ
(
E
R
-
E
L
)
d
E
R
d
E
L
.The opposite tunneling rate can be calculated in a similar way. Thus, the net tunneling current I
=
-
q
[
Γ
L
→
R
-
Γ
R
→
L
] assuming symmetry in the transmission coefficient T
L
R
=
T
R
L [38] can be written as (4)
I
R
L
=
4
π
q
ħ
∫
-
∞
+
∞
‍
T
R
L
E
2
ρ
R
(
E
)
ρ
L
(
E
)
[
f
R
(
E
)
-
f
L
(
E
)
]
d
E
,where we introduce a factor of 2 to take into account the spin. Therefore, using this approach, we can describe the whole system as independent subsystems (the Qds) connected between them by a transmission probability through the dielectric media. Thus, this methodology allows us to write the currents between the different parts of the system.In the above expressions,f
R and f
L are the nonequilibrium energy distribution functions in each side of the barrier. These distributions functions take into account how the energy levels are filled (f
R) or emptied (1
-
f
R); as it is expected the electron transport only occurs if the initial state is filled and the final state is empty (see (3)). However, the distribution function of each part of the system is unknown. Assuming that the distribution functions of the electrodes (left L and right R electrodes) are well described by the equilibrium Fermi Dirac statistics using modified electrochemical potentials, μ
L
-
μ
R
=
q
V where V is the applied bias voltage, the problem is reduced to find the nonequilibrium distributions functions of each Qd (n
i).From the definition of the total charge,N
i inside the ith Qd can be expressed as (5)
N
i
=
∫
‍
ρ
i
(
E
)
n
i
(
E
)
d
E
.Here, we are going to redefine the notation used to describe the distribution functions. The distribution function for each Qd is n
i while we reserve f
L and f
R for the distribution function of the left and right leads, respectively. We can write the evolution charge in time for each Qd as N
i
=
∑
j
‍
∫
‍
I
j
i
d
t, where the subscript j takes into account all the elements that are linked to the ith Qd. Thus, from (5) we can write the evolution of charge in time as a function of the total net current flux for each Qd. This set of integrodifferential equations has a similar form as a usual rate equation. In the following, the different current terms and the elements that appear in (4) are going to be discussed.


#### 2.2.1. Electron and Holes Current Terms

Since the evolution charge in time of each Qd can be written as a function of the net current flux, it is needed to determine all the current contributions. The current contributions can be of two types; the Qds have leads contributions and also neighbors Qd current contributions. These two types of current have the same form in (4), but in each case we need to use the correct distribution function.From the point of view of the nature of these contributions, we have three different processes [39]. In Figure 2(a) we show the scheme of the different tunneling processes. The first term corresponds to electron tunneling from conduction band to conduction band (ECB). The second one is an electron tunneling from valence band to conduction band (EVB). Since the transmission coefficients are symmetric, this process also involves the inverse case, tunneling from conduction band to valence band. The last process is related to the holes: hole tunneling from valence band to valence band (HVB).Figure 2
(a) Schematics of the different tunneling processes. Electron from conduction band to conduction band (ECB), electron from valence band to conduction band (EVB), and tunneling from valence band to valence band (HVB) processes. Energy band diagram for the tunneling processes under polarization. (b) If the incident electron energy (E) is less than the modified energy barrier (q
ϕ
0
-
q
V
L
R) we use a direct tunnel expression. (c) If the incident electron energy (E) is greater than the modified energy barrier (q
ϕ
0
-
q
V
L
R), we use a Fowler-Nordheim expression. These expressions depend on the incident electron energy but also depend on the polarization voltage.
(a)
(b)
(c)One important point is how we treat the distinction between electrons and holes. From a physical point of view, the hole conduction can be viewed as electron conduction restricted to the valence band. Therefore, we can consider only electron transport but taking into account the conduction and valence band contributions to the current. Thus, we only need to consider the changing number of electrons in these two bands. Therefore, the time charge evolution equations for each Qd can be written as(6)
d
N
i
d
t
=
∫
-
∞
+
∞
‍
T
ECB
2
ρ
L
ρ
i
CB
(
f
L
-
n
i
)
d
E
+
∫
-
∞
+
∞
‍
T
HVB
2
ρ
L
ρ
i
VB
(
f
L
-
n
i
)
d
E
︸
Left
lead
contribution
+
∫
-
∞
+
∞
‍
|
T
ECB
|
2
ρ
R
ρ
i
CB
(
f
R
-
n
i
)
d
E
+
∫
-
∞
+
∞
‍
T
HVB
2
ρ
R
ρ
i
VB
(
f
R
-
n
i
)
d
E
︸
Right
lead
contribution
,
Neighboring
Qds
contribution
+
∑
j
,
j
≠
i
N
‍
∫
-
∞
+
∞
‍
T
ECB
2
ρ
i
BC
ρ
j
BC
(
n
j
-
n
i
)
d
E
+
∫
-
∞
+
∞
‍
T
HVB
2
ρ
i
VB
ρ
j
VB
(
n
j
-
n
i
)
d
E
+
∑
j
,
j
≠
i
N
‍
∫
-
∞
+
∞
‍
T
EVB
2
ρ
i
BC
ρ
j
VB
(
n
j
-
n
i
)
d
E
+
∫
-
∞
+
∞
‍
T
EVB
2
ρ
j
BC
ρ
i
VB
(
n
j
-
n
i
)
d
E
,wherei
=
1
⋯
N and we take into account all the contributions for an arbitrary ith Qd. The first pair of elements is related to the left lead contribution and the electron and the hole contributions. For simplicity, we assume infinite metallic leads; therefore, we only write the continuum DOS of the leads (ρ
L); meanwhile, in the Qd we write the DOS in separate terms, conduction (ρ
i
BC) and valence (ρ
i
VB) bands. In the next section, we will show that the Qds have discrete electron/hole energy levels instead of conduction or valence bands but this expression will not change. Similar contribution is obtained for the right lead. In these two contributions we use the Fermi Dirac distribution function to describe the leads with μ
L
-
μ
R
=
-
q
V
d
s electrochemical potentials. In each current term, we use the appropriate transmission coefficient. The last two pairs of current terms represent the current from the neighbor Qds. The subscript “j” runs over all the Qds except the Qd that we are considering. In these terms we take into account the different processes, tunneling from the conduction band (CB) to conduction band (CB) and tunneling from valence band (VB) to valence band (VB). We also need to describe the tunneling that mixes the bands (EVB processes) in two ways, from CB to VB and from VB to CB. As it is easy to see, these processes cannot occur at the same time but it is important to take both into account.The set of equations, (6), can be solved for the steady state. Under our assumption that there is no inelastic scattering, the system can be written in a matrix form and solved for each energy step to obtain the nonequilibrium distribution function for each Qd (n
i).


#### 2.2.2. Transmission Coefficients

From (2) the transmission coefficient |
T
L
R
|
2 is defined as the tunnel probability of the electrons crossing through the dielectric media. The tunneling probability is a strong function of the parallel component to the junction interface energy E
∥. At each particular total energy E, the DOS with a zero E
∥ component is heavily weighted by the tunneling probability in (4). Therefore the DOS only takes into account states with k
∥
≈
0. This approximation may in part be a justification for ignoring the tunneling electron momentum in (3) [40].Under our assumptions, we use the semiclassical and one-dimensional WKB approximation [36] for the transmission coefficient |
T
(
E
)
|
2: (7)
T
E
2
=
exp
⁡
-
2
ħ
∫
x
1
x
2
‍
2
m
diel
*
V
x
′
-
E
d
x
′
,where x
1 and x
2 are the classical turning points, m
diel
* is the effective dielectric mass, and V
(
x
) is the potential barrier of the dielectric material. This barrier is the difference between the bands of the Qd and the dielectric matrix. The transmission coefficients have been derived taking into account the effect of the electric field in the interface, E
diel. The transmission coefficient can be separated in three regions; for incident electrons with less energy than the modified height of the barrier we use a direct tunnel expression. This expression considers that the electrons see a trapezoidal potential barrier. When the incident electrons have energies between the modified height and the total barrier height we use the Fowler-Nordheim expression in which the electrons see a triangular potential barrier. Finally, for incident electrons with energy greater than the barrier we do not assume scattering; therefore, we assign |
T
(
E
)
|
2
=
1. This last case corresponds to elastic transport through the conduction band of the dielectric matrix and only occurs for large bias voltages. The transmission coefficient can be written as (8)
T
E
2
=
exp
⁡
-
4
2
m
diel
*
3
ħ
q
E
diel
q
ϕ
1
-
E
3
/
2
-
q
ϕ
0
-
E
3
/
2
hhhhhhhhhhhhhhhhhhhlhhhhhhh
for
q
ϕ
0
≥
E
,
exp
⁡
-
4
2
m
diel
*
3
ħ
q
E
diel
q
ϕ
1
-
E
-
E
c
,
1
3
/
2
hhhhhhhhhhhhhhhhhlhhhhhhhhh
for
q
ϕ
1
≥
E
hhhhhhhhhhhhhhhhhlhhhhhhhfor
q
ϕ
1
hh
≥
q
ϕ
0
1
hhhhhhhhhhhhhhhhlhhhhhhhhh
for
E
≥
q
ϕ
1
.The electric field is defined as E
diel
=
(
E
c
,
1
-
E
c
,
2
+
q
Δ
ϕ
)
/
q
d, where q
ϕ
1 is the potential barrier height, q
ϕ
0 is the modified potential barrier height, d is the tunneling distance, E
c
,
1
-
E
c
,
2 is q times the electrostatic potential between the two elements, and q
Δ
ϕ is the work function difference. In Figure 2 a scheme of the barrier is shown under external polarization and the two different tunneling mechanism are also shown.In a similar way the transmission coefficients for the holes can be derived. In that case, the potential barrier is the difference between the valence bands of the Qd and the dielectric matrix.


#### 2.2.3. Density of States

Another parameter that appears in the expression of the tunneling currents (4) is the density of states (DOS) of both sides of the barrier. As a first order approximation, we propose a simplified model to represent the discrete energy levels in the Qds. We treat each Qd as a finite spherical potential well. The height of the well is the difference between the conduction band energy level of the dielectric matrix and the one of the material that forms the Qd.Solving the spherical Schrödinger equation inside the well forl
=
0 we obtain the typical binding states [41]. The number of binding states and their energetically position depend on the height of the well V
0, the radius R, and the electron effective mass m
Qd
* and m
diel
* (mass inside the Qd and in the dielectric media, resp.). Imposing continuity of the wave function and its first derivate in r
=
R the equation that determines the binding states is (9)
cot
x
=
-
σ
0
x
2
-
m
diel
*
m
Qd
*
,where σ
0
=
(
2
m
diel
*
V
0
/
ħ
2
)
(
R
)
2 and x
=
(
2
m
Qd
*
/
ħ
2
)
(
R
)
2
E. This equation can be solved using a Newton-Raphson algorithm that gives us a discrete energy levels ϵ
i
′. Since in the Schrödinger equation the zero energy origin is located at the bottom of the well, we need to shift the energy ϵ
i
′ in order to have a common Fermi level. We obtain (10)
ρ
E
i
CB
=
∑
i
n
‍
δ
(
E
-
E
displ
-
ϵ
i
′
)
,where n is the number of bounding states in the ith Qd. The value of E
displ is half the size of the bulk material gap where we assume the Fermi level is placed. Now, we define ϵ
i
=
ϵ
i
′
+
E
displ. Similar treatment is done for the hole binding states of the Qd.Up to now we treated each Qd as an independent part of the system, but the Qds are coupled between them. This effect is introduced assuming a broadening of the discrete energy levels of the Qds. The standard way to introduce the broadening of the energy levels as a consequence of contacts is to assign a Lorentzian shape to each discrete energy level [42]: (11)
δ
(
E
-
ϵ
)
⟶
γ
/
2
π
E
-
ϵ
2
+
γ
/
2
2
,where γ is the broadening of the level and it is related with the tunnel probabilities. Therefore, the total DOS for each Qd is the total sum of the energy levels taking into account the electron and hole binding states: (12)
ρ
i
=
∑
i
n
‍
γ
/
2
π
E
-
ϵ
i
2
+
γ
/
2
2
.We have used a simplified model in order to describe the DOS structure of the Qds but the proposed approach allows to use more complicated DOS obtained usingab initio models. Therefore, this model is suitable to study the transport properties for several Qds materials and taking advantage of the atomistic theories as we have demonstrated in [43].


### 2.3. Potential Profile

Up to now we have only computed the nonequilibrium distribution function of electrons inside each Qd. Therefore, the Qds can be charged (or loose their charge) and these charge variations will affect the local potential of each Qd (V
i). As we can see in Figure 1(c) each junction is modeled as a current tunnel junction in parallel with a capacity. These capacities represent the electrostatic influence between the different parts of the system. Therefore, each Qd has a local potential due to the applied bias voltage. Since each Qd can be charged, we need to solve the Poisson equation: (13)
∇
→
·
(
ε
r
∇
→
V
i
)
=
-
q
Δ
N
i
Ω
,where ε
r is the relative permittivity of the dielectric media and Ω is the Qd volume. Δ
N
i is the change in the number of electrons, calculated with respect to the number of electrons N
0 initially in the ith Qd. The potential energy of each Qd is U
i
=
-
q
V
i. The inclusion of the charge term takes into account the carrier interaction inside each Qd. The result of this approach corresponds to the Hartree-term, that is, first approximation of the carrier-carrier interaction using a mean-field level treatment [44].The general solution of the potential energy of theith Qd involves the different capacitive coupling between it and its surrounding and its charge increasing [45] and it can be written as (14)
U
i
=
∑
j
≠
i
‍
C
i
j
C
tot
,
i
(
-
q
V
j
)
+
q
2
C
tot
,
i
Δ
N
i
,where the subscript j runs over all the components of the system, C
i
j is the capacitive coupling between the different components, and C
tot
,
i
=
∑
j
,
j
≠
i
‍
C
i
j is the total capacitive coupling of ith Qd. The charge energy constant U
0
i
=
q
2
/
C
tot
,
i is the potential increase as a consequence of the injection of one electron into the Qd. Equation (14) is a set of equations (one equation per dot) and the first term of (14), the Laplace term U
i
L, can be written in a matrix form as (15)
U
1
L
⋮
U
N
L
=
1
C
tot
,
1
0
0
0
⋮
0
0
⋯
1
C
tot
,
N
×
C
1
Lead
⋮
C
N
Lead
-
q
V
d
s
0
C
1,2
⋯
C
1
,
N
C
2,1
0
⋯
C
2
,
N
⋮
⋮
⋮
⋮
C
N
,
1
C
N
,
2
⋯
0
h
h
h
-
0
C
1,2
⋯
C
1
,
N
C
2,1
0
⋯
C
2
,
N
⋮
⋮
⋮
⋮
C
N
,
1
C
N
,
2
⋯
0
q
V
1
⋮
q
V
N
.The first term of the previous equation is the electrostatic influence of the lead in which the bias voltage (V
d
s) is applied; meanwhile, the second term is the electrostatic coupling with the neighbor Qd. The neighbors capacitive matrix is defined as N
×
N symmetric matrix with zero in the diagonal terms. Both terms are multiplied by the inverse of the total Qd capacity.The effects of the local potential on each Qd should be computed in the Qd DOSρ
i
(
E
)
→
ρ
i
(
E
-
U
i
) shifting the position of the energy levels. This fact modifies the Qd charge and the currents. In (14) we observe that the local potential depends on the increasing charge density but at the same time the charge depends on the DOS which is modified by the local potential. These considerations impose a self-consistent solution of (5) and (14). The self-consitent solution of the system can be summarized as follows.(1)
For a given external bias voltage, (6) are solved and the nonequilibrium distribution functions (n
i) of each Qd are obtained. The Qds charge are evaluated as N
i
=
∫
‍
ρ
i
(
E
)
n
i
(
E
)
d
E.
(2)
The local potential in each Qd is obtained using (14).
(3)
The DOS of the Qds is shifted according to their respective local potentialρ
i
(
E
)
→
ρ
i
(
E
-
U
i
). The transmission coefficients also change.
(4)
Repeat until the potential of the Qds converge.Before we present in detail the code implementation, the capacitive couplings among the different parts are described.


#### 2.3.1. Capacitive Elements

A realistic modelization of the capacitive coupling between the different parts of the systems [46] is needed, since the electron needs available states in the Qds in order to have transport and the DOS of each Qd depends on the local potential. Therefore, the position of the energy levels with the applied bias voltage plays an important role in the determination of the I-V curve.We use the analytical relationship for a sphere to conduct plane capacitance to model the capacitance between the leads and the Qd, which is(16)
C
i
Lead
=
4
π
ε
r
r
2
-
R
2
∑
n
=
1
∞
‍
1
sinh
⁡
(
n
arccosh
(
r
/
R
)
)
,where ε
r is the permittivity of the dielectric media, R is the Qd radius, and r is the distance between the plane and the center of the Qd.For the case of interdot capacitances (C
i
,
j) there is no analytical expression for the capacitance that takes into account Qds of different radii. We use the numerical method of image charges to calculate interdot capacitance between Qds of different sizes. In Figure 3 we show the two capacitive terms as a function of the distance and for different Qd radii.Figure 3
(a) Electrode-Qd capacity (C
i
Lead) for different Qds radii as a function of the distance. (b) Qd-Qd capacity (C
i
j) for different R
2 radii; the radius of one Qd is hold at R
1
=
1 nm. In both plots we use ε
r
=
3.9
ε
0, where ε
0 is the vacuum permittivity.
(a)
(b)


#### 2.3.2. Code Implementation

Since the methodology was presented and realistic parameters have been used, putting all together, we can describe the ballistic transport through an array of Qds embedded in an insulator matrix. We have to note that the inputs for the transport code are only material parameters and the geometry of the device. A summary of the different inputs is shown in Figure4(a).Figure 4
(a) List of the input parameters that describe the device. (b) Scheme flowchart of the self-consistent transport methodology. (c) Computational time versus number of simulated Qds. The time is referred for a single voltage point.
(a)
(b)
(c)Now, we are going to describe briefly the code implementation and some computational strategies that we have used in order to create the computational tool. The code is divided in 3 main parts.(1)
Input parameters: define the geometrical and the material parameters that form the device. The number of Qds is also defined. Calculate all the voltage independent parameters.
(2)
SCF process: start the voltage loop. For each voltage point, the SCF process is repeated until it converges to the desired error.
(3)
Output: calculate the output values.The scheme flowchart of the code is presented in Figure 4(b). The first part of the code generates the system, the Qd arrangement, and takes into account the material parameters needed to describe the Qds and the insulator matrix. Moreover, the capacitive couplings, the binding electron/hole energy levels, and the tunneling distances are also calculated.For each external bias voltage point, the self-consistent methodology is used to obtain a simultaneous convergence of the local potential and the charge in each Qd. In order to accelerate the convergence of this set of equations, we have implemented Anderson’s method [47, 48]. This procedure consists of the minimization of the “distance” between the input (U
in) and the output (U
out) potentials. The distance is defined as (17)
D
[
U
out
,
U
in
]
=
U
out
-
U
in
∣
U
out
-
U
in
1
/
2
=
F
∣
F
1
/
2
,where 〈
U
out
-
U
in
| is a vector of the differences for each Qd. Two “average” mixed potentials are defined for each iteration step as (18)
U
-
in
out
=
(
1
-
β
)
U
in
out
m
+
β
U
in
(
out
)
(
m
-
1
)
,where m is the current iteration. The aim is to obtain the “best” β value for the current iteration which minimize the distance between these two average quantities. It reads as (19)
β
=
F
(
m
)
∣
F
(
m
)
-
F
(
m
-
1
)
D
2
[
F
m
,
F
m
]
.Finally, to obtain the new guess for the next iteration, we simply mix the average U
-
in and U
-
out potentials: (20)
U
(
m
+
1
)
=
(
1
-
α
)
U
-
out
m
+
α
U
-
in
m
,where α is chosen empirically. This scheme is implemented in the code as follows.(1)
Initialize the variables with the last value of the previous voltage point. For the first voltage point the variables start with zero.
(2)
Calculate the solution of the Poisson equationU
out for given U
in.
(3)
Calculate the value ofβ.
(4)
Calculate the average mixing for the inputU
-
in.
(5)
Calculate the average value for the outputU
-
out.
(6)
Do the simple mixing between the input and output potentials.
(7)
Save values for the next iteration.
(8)
Repeat steps (2)–(7) until the desired convergence is achieved.Once the convergence has been achieved, the outputs can been obtained for this voltage step. In a final staged part, the outputs (Qd occupancy versus voltage, current versus voltage, and local potential versus voltage) are calculated for each voltage point and saved in a matrix structure. The process repeats until all the bias voltage steps have been done using the previous potential results as the initial guess for the next bias voltage iteration. In Figure 4(c) we show the computational time needed to obtain results for one voltage point as a function of the number of Qds. The computational time grows with the number of Qds but it is still reasonable and allows simulating large Qd arrays.The methodology has been explained in depth to enable the interested reader to create his/her own code and reproduce the following results. However, implementation of the code for specific devices is available (The developed code is available under agreement in contact with the authorsillera@el.ub.edu).


## 2.1. The Structure

Figure1(a) shows the basic building blocks of our device which is in essence an insulator layer sandwiched between two metallic electrodes. Inside the insulator layer, a random distribution of N Qds can be inserted. This is the classical structure that is obtained due to the fabrication processes, a superlattice of insulator-semiconductor bilayers. Although a single connected Qd to the leads has been obtained creating a so-called single electron transistor (SET) [29], the research mainstream is focused on the properties of structures with many Qds to create nonvolatile memories [30], light-emitting devices [31], or solar cells devices [32].Figure 1
(a) The basic building block that forms the devices based on Qds: an array of Qds (they could be ordered or randomly dispersed) embedded in an insulator matrix sandwiched between two electrodes (leads). (b) Energy band scheme of the structure. The electron starts in the left lead and crosses through the dielectric matrix by tunneling processes assisted by the intermediate Qds. (c) Scheme of the equivalent electrical circuit of the device presented in (a). The Qds (circles) are connected between them and the leads (color blocks) by tunneling junctions. The tunneling junction is described as a capacitor in parallel with a current source.
(a)
(b)
(c)Due to the fabrication processes, the insulator thickness is large enough to avoid direct tunneling between the electrodes (or leads). Therefore, the electron current needs to pass through the Qds. In Figure1(b) the energy band diagram of the system is shown. The Qds are presented as potential wells in the dielectric energy band and they assist the tunneling processes. Thus, a correct description of the tunneling processes among the Qds and Qds-electrodes is needed. These tunneling processes can be described by tunnel junctions. From the electrical point of view, a single tunnel junction is described by a capacitance and a current source that depends on the voltage drop in the junction. The equivalent electrical circuit is shown in Figure 1(c) for an arbitrary Qd array.From the electrical scheme, two coupling equations govern the response of the system: on one hand we can write the charge conservation in each Qd. In the steady state, the charge conservation for this system is analogous to the electronic Kirchhoff’s current law. Thus, using the decomposition of the tunneling junctions in current sources, we can write(1)
0
=
∑
j
‍
I
i
j
i
∈
N
,for the ith Qd (equivalent to a node), where N is the number of Qds and the summation takes into account all the linked elements to this Qd. On the other hand, the current through a junction will depend on the voltage drop in the tunneling junction. Therefore, a second equation is necessary to describe the voltage in each Qd, the same as the Kirchhoff’s voltage law. First of all, we are going to show analytical expressions for the tunnel currents through a junctions and, after that, we will calculate the potential in each Qd using the Poisson equation.


## 2.2. Current through a Tunneling Junction

A usual method employed to describe the tunneling processes in devices is the tunneling Hamiltonian approach (also called Transfer Hamiltonian approach). This theory, thoroughly studied by many authors [33–35], treats the tunnelling events as a perturbation. The matrix coefficient |
T
L
R
|
2 quantifies the probability for a particle to transfer from a state of the left side of the barrier to a state of the right side by a tunnel process.The tunneling rate is determined using time-dependent perturbation theory by considering the electron from one side of the barrier as initial state and the electron on the other side as a final state. The tunneling rate from the left to the right states (both are considered as a part of a continuum) can be calculated using the Fermi’s golden rule [36]:(2)
d
2
W
k
→
L
→
k
→
R
=
2
π
ħ
T
L
R
2
ρ
R
(
E
R
)
ρ
L
(
E
L
)
δ
(
E
R
-
E
L
)
d
E
R
d
E
L
,where ρ
L and ρ
R are the density of states of the left and right side. From this expression, we can see that we only consider ballistic transport. This means that the electron does not suffer energy loss scattering processes when it moves through the barrier. Introducing the energy distribution function in each part of the barrier we can evaluate the total tunneling rate from all occupied states on the left to all unoccupied states on the right [37]: (3)
Γ
L
→
R
=
∫
-
∞
+
∞
‍
ρ
L
E
L
f
L
E
L
·
∫
-
∞
+
∞
‍
2
π
ħ
T
L
R
2
1
-
f
R
E
R
h
h
l
h
∫
-
∞
+
∞
‍
2
π
ħ
·
ρ
R
(
E
R
)
δ
(
E
R
-
E
L
)
d
E
R
d
E
L
.The opposite tunneling rate can be calculated in a similar way. Thus, the net tunneling current I
=
-
q
[
Γ
L
→
R
-
Γ
R
→
L
] assuming symmetry in the transmission coefficient T
L
R
=
T
R
L [38] can be written as (4)
I
R
L
=
4
π
q
ħ
∫
-
∞
+
∞
‍
T
R
L
E
2
ρ
R
(
E
)
ρ
L
(
E
)
[
f
R
(
E
)
-
f
L
(
E
)
]
d
E
,where we introduce a factor of 2 to take into account the spin. Therefore, using this approach, we can describe the whole system as independent subsystems (the Qds) connected between them by a transmission probability through the dielectric media. Thus, this methodology allows us to write the currents between the different parts of the system.In the above expressions,f
R and f
L are the nonequilibrium energy distribution functions in each side of the barrier. These distributions functions take into account how the energy levels are filled (f
R) or emptied (1
-
f
R); as it is expected the electron transport only occurs if the initial state is filled and the final state is empty (see (3)). However, the distribution function of each part of the system is unknown. Assuming that the distribution functions of the electrodes (left L and right R electrodes) are well described by the equilibrium Fermi Dirac statistics using modified electrochemical potentials, μ
L
-
μ
R
=
q
V where V is the applied bias voltage, the problem is reduced to find the nonequilibrium distributions functions of each Qd (n
i).From the definition of the total charge,N
i inside the ith Qd can be expressed as (5)
N
i
=
∫
‍
ρ
i
(
E
)
n
i
(
E
)
d
E
.Here, we are going to redefine the notation used to describe the distribution functions. The distribution function for each Qd is n
i while we reserve f
L and f
R for the distribution function of the left and right leads, respectively. We can write the evolution charge in time for each Qd as N
i
=
∑
j
‍
∫
‍
I
j
i
d
t, where the subscript j takes into account all the elements that are linked to the ith Qd. Thus, from (5) we can write the evolution of charge in time as a function of the total net current flux for each Qd. This set of integrodifferential equations has a similar form as a usual rate equation. In the following, the different current terms and the elements that appear in (4) are going to be discussed.


### 2.2.1. Electron and Holes Current Terms

Since the evolution charge in time of each Qd can be written as a function of the net current flux, it is needed to determine all the current contributions. The current contributions can be of two types; the Qds have leads contributions and also neighbors Qd current contributions. These two types of current have the same form in (4), but in each case we need to use the correct distribution function.From the point of view of the nature of these contributions, we have three different processes [39]. In Figure 2(a) we show the scheme of the different tunneling processes. The first term corresponds to electron tunneling from conduction band to conduction band (ECB). The second one is an electron tunneling from valence band to conduction band (EVB). Since the transmission coefficients are symmetric, this process also involves the inverse case, tunneling from conduction band to valence band. The last process is related to the holes: hole tunneling from valence band to valence band (HVB).Figure 2
(a) Schematics of the different tunneling processes. Electron from conduction band to conduction band (ECB), electron from valence band to conduction band (EVB), and tunneling from valence band to valence band (HVB) processes. Energy band diagram for the tunneling processes under polarization. (b) If the incident electron energy (E) is less than the modified energy barrier (q
ϕ
0
-
q
V
L
R) we use a direct tunnel expression. (c) If the incident electron energy (E) is greater than the modified energy barrier (q
ϕ
0
-
q
V
L
R), we use a Fowler-Nordheim expression. These expressions depend on the incident electron energy but also depend on the polarization voltage.
(a)
(b)
(c)One important point is how we treat the distinction between electrons and holes. From a physical point of view, the hole conduction can be viewed as electron conduction restricted to the valence band. Therefore, we can consider only electron transport but taking into account the conduction and valence band contributions to the current. Thus, we only need to consider the changing number of electrons in these two bands. Therefore, the time charge evolution equations for each Qd can be written as(6)
d
N
i
d
t
=
∫
-
∞
+
∞
‍
T
ECB
2
ρ
L
ρ
i
CB
(
f
L
-
n
i
)
d
E
+
∫
-
∞
+
∞
‍
T
HVB
2
ρ
L
ρ
i
VB
(
f
L
-
n
i
)
d
E
︸
Left
lead
contribution
+
∫
-
∞
+
∞
‍
|
T
ECB
|
2
ρ
R
ρ
i
CB
(
f
R
-
n
i
)
d
E
+
∫
-
∞
+
∞
‍
T
HVB
2
ρ
R
ρ
i
VB
(
f
R
-
n
i
)
d
E
︸
Right
lead
contribution
,
Neighboring
Qds
contribution
+
∑
j
,
j
≠
i
N
‍
∫
-
∞
+
∞
‍
T
ECB
2
ρ
i
BC
ρ
j
BC
(
n
j
-
n
i
)
d
E
+
∫
-
∞
+
∞
‍
T
HVB
2
ρ
i
VB
ρ
j
VB
(
n
j
-
n
i
)
d
E
+
∑
j
,
j
≠
i
N
‍
∫
-
∞
+
∞
‍
T
EVB
2
ρ
i
BC
ρ
j
VB
(
n
j
-
n
i
)
d
E
+
∫
-
∞
+
∞
‍
T
EVB
2
ρ
j
BC
ρ
i
VB
(
n
j
-
n
i
)
d
E
,wherei
=
1
⋯
N and we take into account all the contributions for an arbitrary ith Qd. The first pair of elements is related to the left lead contribution and the electron and the hole contributions. For simplicity, we assume infinite metallic leads; therefore, we only write the continuum DOS of the leads (ρ
L); meanwhile, in the Qd we write the DOS in separate terms, conduction (ρ
i
BC) and valence (ρ
i
VB) bands. In the next section, we will show that the Qds have discrete electron/hole energy levels instead of conduction or valence bands but this expression will not change. Similar contribution is obtained for the right lead. In these two contributions we use the Fermi Dirac distribution function to describe the leads with μ
L
-
μ
R
=
-
q
V
d
s electrochemical potentials. In each current term, we use the appropriate transmission coefficient. The last two pairs of current terms represent the current from the neighbor Qds. The subscript “j” runs over all the Qds except the Qd that we are considering. In these terms we take into account the different processes, tunneling from the conduction band (CB) to conduction band (CB) and tunneling from valence band (VB) to valence band (VB). We also need to describe the tunneling that mixes the bands (EVB processes) in two ways, from CB to VB and from VB to CB. As it is easy to see, these processes cannot occur at the same time but it is important to take both into account.The set of equations, (6), can be solved for the steady state. Under our assumption that there is no inelastic scattering, the system can be written in a matrix form and solved for each energy step to obtain the nonequilibrium distribution function for each Qd (n
i).


### 2.2.2. Transmission Coefficients

From (2) the transmission coefficient |
T
L
R
|
2 is defined as the tunnel probability of the electrons crossing through the dielectric media. The tunneling probability is a strong function of the parallel component to the junction interface energy E
∥. At each particular total energy E, the DOS with a zero E
∥ component is heavily weighted by the tunneling probability in (4). Therefore the DOS only takes into account states with k
∥
≈
0. This approximation may in part be a justification for ignoring the tunneling electron momentum in (3) [40].Under our assumptions, we use the semiclassical and one-dimensional WKB approximation [36] for the transmission coefficient |
T
(
E
)
|
2: (7)
T
E
2
=
exp
⁡
-
2
ħ
∫
x
1
x
2
‍
2
m
diel
*
V
x
′
-
E
d
x
′
,where x
1 and x
2 are the classical turning points, m
diel
* is the effective dielectric mass, and V
(
x
) is the potential barrier of the dielectric material. This barrier is the difference between the bands of the Qd and the dielectric matrix. The transmission coefficients have been derived taking into account the effect of the electric field in the interface, E
diel. The transmission coefficient can be separated in three regions; for incident electrons with less energy than the modified height of the barrier we use a direct tunnel expression. This expression considers that the electrons see a trapezoidal potential barrier. When the incident electrons have energies between the modified height and the total barrier height we use the Fowler-Nordheim expression in which the electrons see a triangular potential barrier. Finally, for incident electrons with energy greater than the barrier we do not assume scattering; therefore, we assign |
T
(
E
)
|
2
=
1. This last case corresponds to elastic transport through the conduction band of the dielectric matrix and only occurs for large bias voltages. The transmission coefficient can be written as (8)
T
E
2
=
exp
⁡
-
4
2
m
diel
*
3
ħ
q
E
diel
q
ϕ
1
-
E
3
/
2
-
q
ϕ
0
-
E
3
/
2
hhhhhhhhhhhhhhhhhhhlhhhhhhh
for
q
ϕ
0
≥
E
,
exp
⁡
-
4
2
m
diel
*
3
ħ
q
E
diel
q
ϕ
1
-
E
-
E
c
,
1
3
/
2
hhhhhhhhhhhhhhhhhlhhhhhhhhh
for
q
ϕ
1
≥
E
hhhhhhhhhhhhhhhhhlhhhhhhhfor
q
ϕ
1
hh
≥
q
ϕ
0
1
hhhhhhhhhhhhhhhhlhhhhhhhhh
for
E
≥
q
ϕ
1
.The electric field is defined as E
diel
=
(
E
c
,
1
-
E
c
,
2
+
q
Δ
ϕ
)
/
q
d, where q
ϕ
1 is the potential barrier height, q
ϕ
0 is the modified potential barrier height, d is the tunneling distance, E
c
,
1
-
E
c
,
2 is q times the electrostatic potential between the two elements, and q
Δ
ϕ is the work function difference. In Figure 2 a scheme of the barrier is shown under external polarization and the two different tunneling mechanism are also shown.In a similar way the transmission coefficients for the holes can be derived. In that case, the potential barrier is the difference between the valence bands of the Qd and the dielectric matrix.


### 2.2.3. Density of States

Another parameter that appears in the expression of the tunneling currents (4) is the density of states (DOS) of both sides of the barrier. As a first order approximation, we propose a simplified model to represent the discrete energy levels in the Qds. We treat each Qd as a finite spherical potential well. The height of the well is the difference between the conduction band energy level of the dielectric matrix and the one of the material that forms the Qd.Solving the spherical Schrödinger equation inside the well forl
=
0 we obtain the typical binding states [41]. The number of binding states and their energetically position depend on the height of the well V
0, the radius R, and the electron effective mass m
Qd
* and m
diel
* (mass inside the Qd and in the dielectric media, resp.). Imposing continuity of the wave function and its first derivate in r
=
R the equation that determines the binding states is (9)
cot
x
=
-
σ
0
x
2
-
m
diel
*
m
Qd
*
,where σ
0
=
(
2
m
diel
*
V
0
/
ħ
2
)
(
R
)
2 and x
=
(
2
m
Qd
*
/
ħ
2
)
(
R
)
2
E. This equation can be solved using a Newton-Raphson algorithm that gives us a discrete energy levels ϵ
i
′. Since in the Schrödinger equation the zero energy origin is located at the bottom of the well, we need to shift the energy ϵ
i
′ in order to have a common Fermi level. We obtain (10)
ρ
E
i
CB
=
∑
i
n
‍
δ
(
E
-
E
displ
-
ϵ
i
′
)
,where n is the number of bounding states in the ith Qd. The value of E
displ is half the size of the bulk material gap where we assume the Fermi level is placed. Now, we define ϵ
i
=
ϵ
i
′
+
E
displ. Similar treatment is done for the hole binding states of the Qd.Up to now we treated each Qd as an independent part of the system, but the Qds are coupled between them. This effect is introduced assuming a broadening of the discrete energy levels of the Qds. The standard way to introduce the broadening of the energy levels as a consequence of contacts is to assign a Lorentzian shape to each discrete energy level [42]: (11)
δ
(
E
-
ϵ
)
⟶
γ
/
2
π
E
-
ϵ
2
+
γ
/
2
2
,where γ is the broadening of the level and it is related with the tunnel probabilities. Therefore, the total DOS for each Qd is the total sum of the energy levels taking into account the electron and hole binding states: (12)
ρ
i
=
∑
i
n
‍
γ
/
2
π
E
-
ϵ
i
2
+
γ
/
2
2
.We have used a simplified model in order to describe the DOS structure of the Qds but the proposed approach allows to use more complicated DOS obtained usingab initio models. Therefore, this model is suitable to study the transport properties for several Qds materials and taking advantage of the atomistic theories as we have demonstrated in [43].


## 2.2.1. Electron and Holes Current Terms

Since the evolution charge in time of each Qd can be written as a function of the net current flux, it is needed to determine all the current contributions. The current contributions can be of two types; the Qds have leads contributions and also neighbors Qd current contributions. These two types of current have the same form in (4), but in each case we need to use the correct distribution function.From the point of view of the nature of these contributions, we have three different processes [39]. In Figure 2(a) we show the scheme of the different tunneling processes. The first term corresponds to electron tunneling from conduction band to conduction band (ECB). The second one is an electron tunneling from valence band to conduction band (EVB). Since the transmission coefficients are symmetric, this process also involves the inverse case, tunneling from conduction band to valence band. The last process is related to the holes: hole tunneling from valence band to valence band (HVB).Figure 2
(a) Schematics of the different tunneling processes. Electron from conduction band to conduction band (ECB), electron from valence band to conduction band (EVB), and tunneling from valence band to valence band (HVB) processes. Energy band diagram for the tunneling processes under polarization. (b) If the incident electron energy (E) is less than the modified energy barrier (q
ϕ
0
-
q
V
L
R) we use a direct tunnel expression. (c) If the incident electron energy (E) is greater than the modified energy barrier (q
ϕ
0
-
q
V
L
R), we use a Fowler-Nordheim expression. These expressions depend on the incident electron energy but also depend on the polarization voltage.
(a)
(b)
(c)One important point is how we treat the distinction between electrons and holes. From a physical point of view, the hole conduction can be viewed as electron conduction restricted to the valence band. Therefore, we can consider only electron transport but taking into account the conduction and valence band contributions to the current. Thus, we only need to consider the changing number of electrons in these two bands. Therefore, the time charge evolution equations for each Qd can be written as(6)
d
N
i
d
t
=
∫
-
∞
+
∞
‍
T
ECB
2
ρ
L
ρ
i
CB
(
f
L
-
n
i
)
d
E
+
∫
-
∞
+
∞
‍
T
HVB
2
ρ
L
ρ
i
VB
(
f
L
-
n
i
)
d
E
︸
Left
lead
contribution
+
∫
-
∞
+
∞
‍
|
T
ECB
|
2
ρ
R
ρ
i
CB
(
f
R
-
n
i
)
d
E
+
∫
-
∞
+
∞
‍
T
HVB
2
ρ
R
ρ
i
VB
(
f
R
-
n
i
)
d
E
︸
Right
lead
contribution
,
Neighboring
Qds
contribution
+
∑
j
,
j
≠
i
N
‍
∫
-
∞
+
∞
‍
T
ECB
2
ρ
i
BC
ρ
j
BC
(
n
j
-
n
i
)
d
E
+
∫
-
∞
+
∞
‍
T
HVB
2
ρ
i
VB
ρ
j
VB
(
n
j
-
n
i
)
d
E
+
∑
j
,
j
≠
i
N
‍
∫
-
∞
+
∞
‍
T
EVB
2
ρ
i
BC
ρ
j
VB
(
n
j
-
n
i
)
d
E
+
∫
-
∞
+
∞
‍
T
EVB
2
ρ
j
BC
ρ
i
VB
(
n
j
-
n
i
)
d
E
,wherei
=
1
⋯
N and we take into account all the contributions for an arbitrary ith Qd. The first pair of elements is related to the left lead contribution and the electron and the hole contributions. For simplicity, we assume infinite metallic leads; therefore, we only write the continuum DOS of the leads (ρ
L); meanwhile, in the Qd we write the DOS in separate terms, conduction (ρ
i
BC) and valence (ρ
i
VB) bands. In the next section, we will show that the Qds have discrete electron/hole energy levels instead of conduction or valence bands but this expression will not change. Similar contribution is obtained for the right lead. In these two contributions we use the Fermi Dirac distribution function to describe the leads with μ
L
-
μ
R
=
-
q
V
d
s electrochemical potentials. In each current term, we use the appropriate transmission coefficient. The last two pairs of current terms represent the current from the neighbor Qds. The subscript “j” runs over all the Qds except the Qd that we are considering. In these terms we take into account the different processes, tunneling from the conduction band (CB) to conduction band (CB) and tunneling from valence band (VB) to valence band (VB). We also need to describe the tunneling that mixes the bands (EVB processes) in two ways, from CB to VB and from VB to CB. As it is easy to see, these processes cannot occur at the same time but it is important to take both into account.The set of equations, (6), can be solved for the steady state. Under our assumption that there is no inelastic scattering, the system can be written in a matrix form and solved for each energy step to obtain the nonequilibrium distribution function for each Qd (n
i).


## 2.2.2. Transmission Coefficients

From (2) the transmission coefficient |
T
L
R
|
2 is defined as the tunnel probability of the electrons crossing through the dielectric media. The tunneling probability is a strong function of the parallel component to the junction interface energy E
∥. At each particular total energy E, the DOS with a zero E
∥ component is heavily weighted by the tunneling probability in (4). Therefore the DOS only takes into account states with k
∥
≈
0. This approximation may in part be a justification for ignoring the tunneling electron momentum in (3) [40].Under our assumptions, we use the semiclassical and one-dimensional WKB approximation [36] for the transmission coefficient |
T
(
E
)
|
2: (7)
T
E
2
=
exp
⁡
-
2
ħ
∫
x
1
x
2
‍
2
m
diel
*
V
x
′
-
E
d
x
′
,where x
1 and x
2 are the classical turning points, m
diel
* is the effective dielectric mass, and V
(
x
) is the potential barrier of the dielectric material. This barrier is the difference between the bands of the Qd and the dielectric matrix. The transmission coefficients have been derived taking into account the effect of the electric field in the interface, E
diel. The transmission coefficient can be separated in three regions; for incident electrons with less energy than the modified height of the barrier we use a direct tunnel expression. This expression considers that the electrons see a trapezoidal potential barrier. When the incident electrons have energies between the modified height and the total barrier height we use the Fowler-Nordheim expression in which the electrons see a triangular potential barrier. Finally, for incident electrons with energy greater than the barrier we do not assume scattering; therefore, we assign |
T
(
E
)
|
2
=
1. This last case corresponds to elastic transport through the conduction band of the dielectric matrix and only occurs for large bias voltages. The transmission coefficient can be written as (8)
T
E
2
=
exp
⁡
-
4
2
m
diel
*
3
ħ
q
E
diel
q
ϕ
1
-
E
3
/
2
-
q
ϕ
0
-
E
3
/
2
hhhhhhhhhhhhhhhhhhhlhhhhhhh
for
q
ϕ
0
≥
E
,
exp
⁡
-
4
2
m
diel
*
3
ħ
q
E
diel
q
ϕ
1
-
E
-
E
c
,
1
3
/
2
hhhhhhhhhhhhhhhhhlhhhhhhhhh
for
q
ϕ
1
≥
E
hhhhhhhhhhhhhhhhhlhhhhhhhfor
q
ϕ
1
hh
≥
q
ϕ
0
1
hhhhhhhhhhhhhhhhlhhhhhhhhh
for
E
≥
q
ϕ
1
.The electric field is defined as E
diel
=
(
E
c
,
1
-
E
c
,
2
+
q
Δ
ϕ
)
/
q
d, where q
ϕ
1 is the potential barrier height, q
ϕ
0 is the modified potential barrier height, d is the tunneling distance, E
c
,
1
-
E
c
,
2 is q times the electrostatic potential between the two elements, and q
Δ
ϕ is the work function difference. In Figure 2 a scheme of the barrier is shown under external polarization and the two different tunneling mechanism are also shown.In a similar way the transmission coefficients for the holes can be derived. In that case, the potential barrier is the difference between the valence bands of the Qd and the dielectric matrix.


## 2.2.3. Density of States

Another parameter that appears in the expression of the tunneling currents (4) is the density of states (DOS) of both sides of the barrier. As a first order approximation, we propose a simplified model to represent the discrete energy levels in the Qds. We treat each Qd as a finite spherical potential well. The height of the well is the difference between the conduction band energy level of the dielectric matrix and the one of the material that forms the Qd.Solving the spherical Schrödinger equation inside the well forl
=
0 we obtain the typical binding states [41]. The number of binding states and their energetically position depend on the height of the well V
0, the radius R, and the electron effective mass m
Qd
* and m
diel
* (mass inside the Qd and in the dielectric media, resp.). Imposing continuity of the wave function and its first derivate in r
=
R the equation that determines the binding states is (9)
cot
x
=
-
σ
0
x
2
-
m
diel
*
m
Qd
*
,where σ
0
=
(
2
m
diel
*
V
0
/
ħ
2
)
(
R
)
2 and x
=
(
2
m
Qd
*
/
ħ
2
)
(
R
)
2
E. This equation can be solved using a Newton-Raphson algorithm that gives us a discrete energy levels ϵ
i
′. Since in the Schrödinger equation the zero energy origin is located at the bottom of the well, we need to shift the energy ϵ
i
′ in order to have a common Fermi level. We obtain (10)
ρ
E
i
CB
=
∑
i
n
‍
δ
(
E
-
E
displ
-
ϵ
i
′
)
,where n is the number of bounding states in the ith Qd. The value of E
displ is half the size of the bulk material gap where we assume the Fermi level is placed. Now, we define ϵ
i
=
ϵ
i
′
+
E
displ. Similar treatment is done for the hole binding states of the Qd.Up to now we treated each Qd as an independent part of the system, but the Qds are coupled between them. This effect is introduced assuming a broadening of the discrete energy levels of the Qds. The standard way to introduce the broadening of the energy levels as a consequence of contacts is to assign a Lorentzian shape to each discrete energy level [42]: (11)
δ
(
E
-
ϵ
)
⟶
γ
/
2
π
E
-
ϵ
2
+
γ
/
2
2
,where γ is the broadening of the level and it is related with the tunnel probabilities. Therefore, the total DOS for each Qd is the total sum of the energy levels taking into account the electron and hole binding states: (12)
ρ
i
=
∑
i
n
‍
γ
/
2
π
E
-
ϵ
i
2
+
γ
/
2
2
.We have used a simplified model in order to describe the DOS structure of the Qds but the proposed approach allows to use more complicated DOS obtained usingab initio models. Therefore, this model is suitable to study the transport properties for several Qds materials and taking advantage of the atomistic theories as we have demonstrated in [43].


## 2.3. Potential Profile

Up to now we have only computed the nonequilibrium distribution function of electrons inside each Qd. Therefore, the Qds can be charged (or loose their charge) and these charge variations will affect the local potential of each Qd (V
i). As we can see in Figure 1(c) each junction is modeled as a current tunnel junction in parallel with a capacity. These capacities represent the electrostatic influence between the different parts of the system. Therefore, each Qd has a local potential due to the applied bias voltage. Since each Qd can be charged, we need to solve the Poisson equation: (13)
∇
→
·
(
ε
r
∇
→
V
i
)
=
-
q
Δ
N
i
Ω
,where ε
r is the relative permittivity of the dielectric media and Ω is the Qd volume. Δ
N
i is the change in the number of electrons, calculated with respect to the number of electrons N
0 initially in the ith Qd. The potential energy of each Qd is U
i
=
-
q
V
i. The inclusion of the charge term takes into account the carrier interaction inside each Qd. The result of this approach corresponds to the Hartree-term, that is, first approximation of the carrier-carrier interaction using a mean-field level treatment [44].The general solution of the potential energy of theith Qd involves the different capacitive coupling between it and its surrounding and its charge increasing [45] and it can be written as (14)
U
i
=
∑
j
≠
i
‍
C
i
j
C
tot
,
i
(
-
q
V
j
)
+
q
2
C
tot
,
i
Δ
N
i
,where the subscript j runs over all the components of the system, C
i
j is the capacitive coupling between the different components, and C
tot
,
i
=
∑
j
,
j
≠
i
‍
C
i
j is the total capacitive coupling of ith Qd. The charge energy constant U
0
i
=
q
2
/
C
tot
,
i is the potential increase as a consequence of the injection of one electron into the Qd. Equation (14) is a set of equations (one equation per dot) and the first term of (14), the Laplace term U
i
L, can be written in a matrix form as (15)
U
1
L
⋮
U
N
L
=
1
C
tot
,
1
0
0
0
⋮
0
0
⋯
1
C
tot
,
N
×
C
1
Lead
⋮
C
N
Lead
-
q
V
d
s
0
C
1,2
⋯
C
1
,
N
C
2,1
0
⋯
C
2
,
N
⋮
⋮
⋮
⋮
C
N
,
1
C
N
,
2
⋯
0
h
h
h
-
0
C
1,2
⋯
C
1
,
N
C
2,1
0
⋯
C
2
,
N
⋮
⋮
⋮
⋮
C
N
,
1
C
N
,
2
⋯
0
q
V
1
⋮
q
V
N
.The first term of the previous equation is the electrostatic influence of the lead in which the bias voltage (V
d
s) is applied; meanwhile, the second term is the electrostatic coupling with the neighbor Qd. The neighbors capacitive matrix is defined as N
×
N symmetric matrix with zero in the diagonal terms. Both terms are multiplied by the inverse of the total Qd capacity.The effects of the local potential on each Qd should be computed in the Qd DOSρ
i
(
E
)
→
ρ
i
(
E
-
U
i
) shifting the position of the energy levels. This fact modifies the Qd charge and the currents. In (14) we observe that the local potential depends on the increasing charge density but at the same time the charge depends on the DOS which is modified by the local potential. These considerations impose a self-consistent solution of (5) and (14). The self-consitent solution of the system can be summarized as follows.(1)
For a given external bias voltage, (6) are solved and the nonequilibrium distribution functions (n
i) of each Qd are obtained. The Qds charge are evaluated as N
i
=
∫
‍
ρ
i
(
E
)
n
i
(
E
)
d
E.
(2)
The local potential in each Qd is obtained using (14).
(3)
The DOS of the Qds is shifted according to their respective local potentialρ
i
(
E
)
→
ρ
i
(
E
-
U
i
). The transmission coefficients also change.
(4)
Repeat until the potential of the Qds converge.Before we present in detail the code implementation, the capacitive couplings among the different parts are described.


### 2.3.1. Capacitive Elements

A realistic modelization of the capacitive coupling between the different parts of the systems [46] is needed, since the electron needs available states in the Qds in order to have transport and the DOS of each Qd depends on the local potential. Therefore, the position of the energy levels with the applied bias voltage plays an important role in the determination of the I-V curve.We use the analytical relationship for a sphere to conduct plane capacitance to model the capacitance between the leads and the Qd, which is(16)
C
i
Lead
=
4
π
ε
r
r
2
-
R
2
∑
n
=
1
∞
‍
1
sinh
⁡
(
n
arccosh
(
r
/
R
)
)
,where ε
r is the permittivity of the dielectric media, R is the Qd radius, and r is the distance between the plane and the center of the Qd.For the case of interdot capacitances (C
i
,
j) there is no analytical expression for the capacitance that takes into account Qds of different radii. We use the numerical method of image charges to calculate interdot capacitance between Qds of different sizes. In Figure 3 we show the two capacitive terms as a function of the distance and for different Qd radii.Figure 3
(a) Electrode-Qd capacity (C
i
Lead) for different Qds radii as a function of the distance. (b) Qd-Qd capacity (C
i
j) for different R
2 radii; the radius of one Qd is hold at R
1
=
1 nm. In both plots we use ε
r
=
3.9
ε
0, where ε
0 is the vacuum permittivity.
(a)
(b)


### 2.3.2. Code Implementation

Since the methodology was presented and realistic parameters have been used, putting all together, we can describe the ballistic transport through an array of Qds embedded in an insulator matrix. We have to note that the inputs for the transport code are only material parameters and the geometry of the device. A summary of the different inputs is shown in Figure4(a).Figure 4
(a) List of the input parameters that describe the device. (b) Scheme flowchart of the self-consistent transport methodology. (c) Computational time versus number of simulated Qds. The time is referred for a single voltage point.
(a)
(b)
(c)Now, we are going to describe briefly the code implementation and some computational strategies that we have used in order to create the computational tool. The code is divided in 3 main parts.(1)
Input parameters: define the geometrical and the material parameters that form the device. The number of Qds is also defined. Calculate all the voltage independent parameters.
(2)
SCF process: start the voltage loop. For each voltage point, the SCF process is repeated until it converges to the desired error.
(3)
Output: calculate the output values.The scheme flowchart of the code is presented in Figure 4(b). The first part of the code generates the system, the Qd arrangement, and takes into account the material parameters needed to describe the Qds and the insulator matrix. Moreover, the capacitive couplings, the binding electron/hole energy levels, and the tunneling distances are also calculated.For each external bias voltage point, the self-consistent methodology is used to obtain a simultaneous convergence of the local potential and the charge in each Qd. In order to accelerate the convergence of this set of equations, we have implemented Anderson’s method [47, 48]. This procedure consists of the minimization of the “distance” between the input (U
in) and the output (U
out) potentials. The distance is defined as (17)
D
[
U
out
,
U
in
]
=
U
out
-
U
in
∣
U
out
-
U
in
1
/
2
=
F
∣
F
1
/
2
,where 〈
U
out
-
U
in
| is a vector of the differences for each Qd. Two “average” mixed potentials are defined for each iteration step as (18)
U
-
in
out
=
(
1
-
β
)
U
in
out
m
+
β
U
in
(
out
)
(
m
-
1
)
,where m is the current iteration. The aim is to obtain the “best” β value for the current iteration which minimize the distance between these two average quantities. It reads as (19)
β
=
F
(
m
)
∣
F
(
m
)
-
F
(
m
-
1
)
D
2
[
F
m
,
F
m
]
.Finally, to obtain the new guess for the next iteration, we simply mix the average U
-
in and U
-
out potentials: (20)
U
(
m
+
1
)
=
(
1
-
α
)
U
-
out
m
+
α
U
-
in
m
,where α is chosen empirically. This scheme is implemented in the code as follows.(1)
Initialize the variables with the last value of the previous voltage point. For the first voltage point the variables start with zero.
(2)
Calculate the solution of the Poisson equationU
out for given U
in.
(3)
Calculate the value ofβ.
(4)
Calculate the average mixing for the inputU
-
in.
(5)
Calculate the average value for the outputU
-
out.
(6)
Do the simple mixing between the input and output potentials.
(7)
Save values for the next iteration.
(8)
Repeat steps (2)–(7) until the desired convergence is achieved.Once the convergence has been achieved, the outputs can been obtained for this voltage step. In a final staged part, the outputs (Qd occupancy versus voltage, current versus voltage, and local potential versus voltage) are calculated for each voltage point and saved in a matrix structure. The process repeats until all the bias voltage steps have been done using the previous potential results as the initial guess for the next bias voltage iteration. In Figure 4(c) we show the computational time needed to obtain results for one voltage point as a function of the number of Qds. The computational time grows with the number of Qds but it is still reasonable and allows simulating large Qd arrays.The methodology has been explained in depth to enable the interested reader to create his/her own code and reproduce the following results. However, implementation of the code for specific devices is available (The developed code is available under agreement in contact with the authorsillera@el.ub.edu).


## 2.3.1. Capacitive Elements

A realistic modelization of the capacitive coupling between the different parts of the systems [46] is needed, since the electron needs available states in the Qds in order to have transport and the DOS of each Qd depends on the local potential. Therefore, the position of the energy levels with the applied bias voltage plays an important role in the determination of the I-V curve.We use the analytical relationship for a sphere to conduct plane capacitance to model the capacitance between the leads and the Qd, which is(16)
C
i
Lead
=
4
π
ε
r
r
2
-
R
2
∑
n
=
1
∞
‍
1
sinh
⁡
(
n
arccosh
(
r
/
R
)
)
,where ε
r is the permittivity of the dielectric media, R is the Qd radius, and r is the distance between the plane and the center of the Qd.For the case of interdot capacitances (C
i
,
j) there is no analytical expression for the capacitance that takes into account Qds of different radii. We use the numerical method of image charges to calculate interdot capacitance between Qds of different sizes. In Figure 3 we show the two capacitive terms as a function of the distance and for different Qd radii.Figure 3
(a) Electrode-Qd capacity (C
i
Lead) for different Qds radii as a function of the distance. (b) Qd-Qd capacity (C
i
j) for different R
2 radii; the radius of one Qd is hold at R
1
=
1 nm. In both plots we use ε
r
=
3.9
ε
0, where ε
0 is the vacuum permittivity.
(a)
(b)


## 2.3.2. Code Implementation

Since the methodology was presented and realistic parameters have been used, putting all together, we can describe the ballistic transport through an array of Qds embedded in an insulator matrix. We have to note that the inputs for the transport code are only material parameters and the geometry of the device. A summary of the different inputs is shown in Figure4(a).Figure 4
(a) List of the input parameters that describe the device. (b) Scheme flowchart of the self-consistent transport methodology. (c) Computational time versus number of simulated Qds. The time is referred for a single voltage point.
(a)
(b)
(c)Now, we are going to describe briefly the code implementation and some computational strategies that we have used in order to create the computational tool. The code is divided in 3 main parts.(1)
Input parameters: define the geometrical and the material parameters that form the device. The number of Qds is also defined. Calculate all the voltage independent parameters.
(2)
SCF process: start the voltage loop. For each voltage point, the SCF process is repeated until it converges to the desired error.
(3)
Output: calculate the output values.The scheme flowchart of the code is presented in Figure 4(b). The first part of the code generates the system, the Qd arrangement, and takes into account the material parameters needed to describe the Qds and the insulator matrix. Moreover, the capacitive couplings, the binding electron/hole energy levels, and the tunneling distances are also calculated.For each external bias voltage point, the self-consistent methodology is used to obtain a simultaneous convergence of the local potential and the charge in each Qd. In order to accelerate the convergence of this set of equations, we have implemented Anderson’s method [47, 48]. This procedure consists of the minimization of the “distance” between the input (U
in) and the output (U
out) potentials. The distance is defined as (17)
D
[
U
out
,
U
in
]
=
U
out
-
U
in
∣
U
out
-
U
in
1
/
2
=
F
∣
F
1
/
2
,where 〈
U
out
-
U
in
| is a vector of the differences for each Qd. Two “average” mixed potentials are defined for each iteration step as (18)
U
-
in
out
=
(
1
-
β
)
U
in
out
m
+
β
U
in
(
out
)
(
m
-
1
)
,where m is the current iteration. The aim is to obtain the “best” β value for the current iteration which minimize the distance between these two average quantities. It reads as (19)
β
=
F
(
m
)
∣
F
(
m
)
-
F
(
m
-
1
)
D
2
[
F
m
,
F
m
]
.Finally, to obtain the new guess for the next iteration, we simply mix the average U
-
in and U
-
out potentials: (20)
U
(
m
+
1
)
=
(
1
-
α
)
U
-
out
m
+
α
U
-
in
m
,where α is chosen empirically. This scheme is implemented in the code as follows.(1)
Initialize the variables with the last value of the previous voltage point. For the first voltage point the variables start with zero.
(2)
Calculate the solution of the Poisson equationU
out for given U
in.
(3)
Calculate the value ofβ.
(4)
Calculate the average mixing for the inputU
-
in.
(5)
Calculate the average value for the outputU
-
out.
(6)
Do the simple mixing between the input and output potentials.
(7)
Save values for the next iteration.
(8)
Repeat steps (2)–(7) until the desired convergence is achieved.Once the convergence has been achieved, the outputs can been obtained for this voltage step. In a final staged part, the outputs (Qd occupancy versus voltage, current versus voltage, and local potential versus voltage) are calculated for each voltage point and saved in a matrix structure. The process repeats until all the bias voltage steps have been done using the previous potential results as the initial guess for the next bias voltage iteration. In Figure 4(c) we show the computational time needed to obtain results for one voltage point as a function of the number of Qds. The computational time grows with the number of Qds but it is still reasonable and allows simulating large Qd arrays.The methodology has been explained in depth to enable the interested reader to create his/her own code and reproduce the following results. However, implementation of the code for specific devices is available (The developed code is available under agreement in contact with the authorsillera@el.ub.edu).


## 3. Application and Discussion

The above-described self-consistent transport simulator has been applied to study the ballistic electronic transport in different generic devices based on silicon Qds embedded in SiO2 matrix (Si/SiO2 Qds). The parameters used to describe the materials are listed in Table 1. The m
ECB
*, m
EVB
*, and m
HVB
* are the oxide effective masses for the different tunneling processes; the values are extracted from Lee and Hu [39]. m
Qd
,
CB
* and m
Qd
,
VB
* are the electron and hole Si bulk effective masses used to obtain the binding states in the Qd. We consider that the equilibrium Fermi level is in the middle of Si E
gap
⁡; confinement potentials of the Qd are ϕ
1
,
ECB for electrons and ϕ
1
,
HVB for holes, respectively. Finally, ε
r is the relative permittivity of the SiO2 matrix.Table 1
Parameters used in the simulation in order to describe Si Qds embedded inSi
O
2 insulator matrix. m
0 and ε
0 are the free electron mass and the vacuum permittivity, respectively.
m
ECB
* (m
0)
0.40 [23, 24]
ϕ
1
,
ECB (eV)
3.1 [25, 26]
m
EVB
* (m
0)
0.30 [23]
ϕ
1
,
HVB (eV)
-
4.5 [25, 26]
m
HVB
* (m
0)
0.32 [23]
E
gap
⁡ (eV)
1.12 [25, 26]
m
Qd
,
CB
* (m
0)
0.33 [27]
ε
r
Si
O
2 (ε
0)
3.9 [26, 28]
m
Qd
,
VB
* (m
0)
0.28 [27]
ε
r
Si (ε
0)
11.7 [26, 28]The electron transport takes place from the left electrode to the right electrode through the Qds. From a physical point of view, the most general transport condition is that the energy levels of the Qds (ϵ
i) must lie between the electrochemical potentials of the leads (μ
L
-
μ
R
=
-
q
V
d
s); this condition can be summarized as μ
L
≥
ϵ
i
≥
μ
R. The type of transport will depend on the nature of these energy levels and electron or hole binding energy levels. Moreover, in order to have transport between the Qds, overlapping of the Qd DOS is necessary. Free states in the arriving Qd are also needed. These two conditions are directly related to the expression of the tunneling current described by (4). The transmission coefficients are also strongly dependent on the tunneling distance; therefore, some processes are more favored than others. Thus, the electronic transport plays with the transmission probabilities between the different processes and the available states.The total net current will be the sum of the partial tunnel currents among the Qds and the right lead (electron and hole terms). This net current is going to be dependent on the position of the Qds, the tunneling distances, and the alignment of the energy levels (the local potential and the DOS of each Qd).


### 3.1. Single Si/SiO2 Qd

Before we simulate a complete device based on large arrays of Qds, we are going to describe a single Si/SiO2 Qd in different configurations. From this small system, we will show the basis of the ballistic transport and the effects of the geometrical arrangement of the Qds in the final electrical response.In Figure5, we present the obtained I
(
V
) curves and the accumulated charge for a system composed of two leads separated by 5 nm and a single Si/SiO2 Qd connected to them placed at different positions. The position of the Qd, x, is measured from the left lead. An external bias voltage is applied to the right lead whereas the left one is kept as a reference; this means μ
L
=
0 and μ
R
=
-
q
V. Concerning the accumulate charge, it reflects the variation of the number of electrons with respect to the initial number. Therefore, if the charge is negative it implies that the system looses charge (hole accumulation). However, if the charge is positive this implies that the system increases its charge (electron accumulation).Figure 5
A single Si/SiO2 Qd of R = 1.5 nm placed in different positions between the two leads. x is the distance from the left lead to the center of the Qd. The separation among the leads is 5 nm. ((a)-(b)) I
(
V
) curve and accumulated charge for a centered Qd. The hole and electron currents are also shown. ((c)–(e)) I
(
V
) curves for different Qd positions and ((e)-(f)) accumulated charge in the same cases.
(a)
(b)
(c)
(d)
(e)
(f)Figure5(a) shows the total I
(
V
) curve and the hole and electron currents for a Qd connected symmetrically to both leads, which is x = 2.5 nm. Concerning the partial currents (electron and hole contributions), the electron current is the dominating term since the electron barrier (3.1 eV) is smaller than the hole barrier (4.5 eV). Besides, the opening of the discrete electron/hole conductive channels is clearly visible in the current steps at different voltages due to the position of the discrete electron/hole energy levels. Since the system is symmetrically coupled to the leads, the total current is symmetric in both polarization directions.Regarding the accumulated charge, in Figure5(b), the Qd remains practically uncharged. The obtained trend reflects the position of the electron and hole energy states respect to the equilibrium Fermi level and the effect of the self-charge. An electron state is the first energy level that starts to be conductive but a hole conductive channel is opened immediately being the Qd accumulated charge the difference between the electron and hole fluxes.TheI
(
V
) curves and the accumulated charge are shown in Figures 5(c)-5(d) and 5(e)-5(f) for the nearest case to the left lead (x
<
2.5) and for the right lead (x
>
2.5), respectively. As can be seen, the symmetry in the total I
(
V
) curve has been broken due to the different capacitive coupling among the leads [49]. Moreover, the Qd gains/loses net charge. Besides, we must note that the obtained trends for the x
<
2.5 cases are complementary to the x
>
2.5 cases.When an external polarization is applied, different incoming/outgoing fluxes are created and the occupation of the states differs from the equilibrium case. Besides, the transport takes place on the energy region betweenμ
L and μ
R. Thus, only the energy states placed in this energy region can gain or lose charge. From the rate equation model, the Qd nonequilibrium distribution function in the steady state can be viewed as a balance between the two leads that strongly depends on the transmission coefficients and the occupation of the leads. Since the leads occupation is well described by the electrochemical potentials μ
L
/
R, the Qd energy level occupation is dominated by the highest transmission coefficient. Therefore, the lead connected to the Qd with the highest transmission coefficient dominates the Qd occupation. Since the transmission coefficients are strongly dependent on the tunneling distance, for lower voltages, the closer lead will dominate the final response. However, for larger voltages, the band bending of the wider oxide increases giving an smaller effective tunneling distance and the Qd recovers its initial charge.


### 3.2. Electron Transport through Si/SiO2 Qds Arrays

In order to deal with devices based on Qd arrays, we are going to study the electron transport in a multilayer structure. This part is going to be one of the most important building blocks of the devices and, therefore, it is important to have a good characterization of it. From the experimental point of view, the superlaticce approach (SL) [50–52] allows creating these kind of structures: layers of Qds separated by insulator layers and size-controlled Qds.Here, we present theI
(
V
) and the total accumulated charge for a multilayer structure. The structure is formed by 2 Si Qds layers. The leads are placed perpendicular to them; therefore, the electronic transport takes place laterally. Both layers are spaced by 5 nm and the layer size is 20 nm length and 20 nm width. We simulate 10 randomly distributed Qds per layer and for the Qds radii we use a normal distribution with 1.75 nm mean value and 0.3 nm standard deviation. A top view scheme of the system is shown in the inset of Figure 6(b).Figure 6
(a) ObtainedI
(
V
) curves (electron/hole and total currents) for two layers of 10 Si/SiO2 Qds. The system arrangement is described in the text. The Qds radius distribution is also shown in the inset. (b) Total accumulated charge (∑
i
‍
Δ
N
i where i
=
1
⋯
20) of the structure as a function of the external bias voltage. In the inset, a top view of the system is presented.
(a)
(b)In Figure6(a), the total I
(
V
) curve and the electron and hole currents terms are shown. From the I
(
V
) curves, two different regimes are obtained: for low and medium voltage ranges, the current reflects the discrete nature of the Qds energy spectra. However, for the largest voltages, the continuous part of the Qd DOS begins to be conductive and the current increases in a continuous form loosing its step-like shape. Concerning the electron/hole current components, the electron term dominates since its potential barrier is lower than the hole one. Moreover, the stepping behavior is still present in the current reflecting the aperture of the different conductive channels. In this case, the transport involves several tunneling processes between the Qds. Therefore, the energy level alignment is necessary. However, the negative differential resistance (NDR) is not clearly visible since there are many electron/hole pathways and the sum of the different conduction channels masks this effect [22].Regarding the accumulated charge plotted in Figure6(b), it is the total sum of the accumulated charge in each Qd (∑
i
‍
Δ
N
i where i
=
1
⋯
20). The geometrical disposition of the Qds strongly affects the final response of the system as we have shown in the previous section. For this system, since there are many electronic conduction channels and a complex capacitive interaction between the Qds, it is hard to explain clearly the accumulated charge trend. Reflecting that, the geometrical disposition of the Qds influences directly the final response of the system.


## 3.1. Single Si/SiO2 Qd

Before we simulate a complete device based on large arrays of Qds, we are going to describe a single Si/SiO2 Qd in different configurations. From this small system, we will show the basis of the ballistic transport and the effects of the geometrical arrangement of the Qds in the final electrical response.In Figure5, we present the obtained I
(
V
) curves and the accumulated charge for a system composed of two leads separated by 5 nm and a single Si/SiO2 Qd connected to them placed at different positions. The position of the Qd, x, is measured from the left lead. An external bias voltage is applied to the right lead whereas the left one is kept as a reference; this means μ
L
=
0 and μ
R
=
-
q
V. Concerning the accumulate charge, it reflects the variation of the number of electrons with respect to the initial number. Therefore, if the charge is negative it implies that the system looses charge (hole accumulation). However, if the charge is positive this implies that the system increases its charge (electron accumulation).Figure 5
A single Si/SiO2 Qd of R = 1.5 nm placed in different positions between the two leads. x is the distance from the left lead to the center of the Qd. The separation among the leads is 5 nm. ((a)-(b)) I
(
V
) curve and accumulated charge for a centered Qd. The hole and electron currents are also shown. ((c)–(e)) I
(
V
) curves for different Qd positions and ((e)-(f)) accumulated charge in the same cases.
(a)
(b)
(c)
(d)
(e)
(f)Figure5(a) shows the total I
(
V
) curve and the hole and electron currents for a Qd connected symmetrically to both leads, which is x = 2.5 nm. Concerning the partial currents (electron and hole contributions), the electron current is the dominating term since the electron barrier (3.1 eV) is smaller than the hole barrier (4.5 eV). Besides, the opening of the discrete electron/hole conductive channels is clearly visible in the current steps at different voltages due to the position of the discrete electron/hole energy levels. Since the system is symmetrically coupled to the leads, the total current is symmetric in both polarization directions.Regarding the accumulated charge, in Figure5(b), the Qd remains practically uncharged. The obtained trend reflects the position of the electron and hole energy states respect to the equilibrium Fermi level and the effect of the self-charge. An electron state is the first energy level that starts to be conductive but a hole conductive channel is opened immediately being the Qd accumulated charge the difference between the electron and hole fluxes.TheI
(
V
) curves and the accumulated charge are shown in Figures 5(c)-5(d) and 5(e)-5(f) for the nearest case to the left lead (x
<
2.5) and for the right lead (x
>
2.5), respectively. As can be seen, the symmetry in the total I
(
V
) curve has been broken due to the different capacitive coupling among the leads [49]. Moreover, the Qd gains/loses net charge. Besides, we must note that the obtained trends for the x
<
2.5 cases are complementary to the x
>
2.5 cases.When an external polarization is applied, different incoming/outgoing fluxes are created and the occupation of the states differs from the equilibrium case. Besides, the transport takes place on the energy region betweenμ
L and μ
R. Thus, only the energy states placed in this energy region can gain or lose charge. From the rate equation model, the Qd nonequilibrium distribution function in the steady state can be viewed as a balance between the two leads that strongly depends on the transmission coefficients and the occupation of the leads. Since the leads occupation is well described by the electrochemical potentials μ
L
/
R, the Qd energy level occupation is dominated by the highest transmission coefficient. Therefore, the lead connected to the Qd with the highest transmission coefficient dominates the Qd occupation. Since the transmission coefficients are strongly dependent on the tunneling distance, for lower voltages, the closer lead will dominate the final response. However, for larger voltages, the band bending of the wider oxide increases giving an smaller effective tunneling distance and the Qd recovers its initial charge.


## 3.2. Electron Transport through Si/SiO2 Qds Arrays

In order to deal with devices based on Qd arrays, we are going to study the electron transport in a multilayer structure. This part is going to be one of the most important building blocks of the devices and, therefore, it is important to have a good characterization of it. From the experimental point of view, the superlaticce approach (SL) [50–52] allows creating these kind of structures: layers of Qds separated by insulator layers and size-controlled Qds.Here, we present theI
(
V
) and the total accumulated charge for a multilayer structure. The structure is formed by 2 Si Qds layers. The leads are placed perpendicular to them; therefore, the electronic transport takes place laterally. Both layers are spaced by 5 nm and the layer size is 20 nm length and 20 nm width. We simulate 10 randomly distributed Qds per layer and for the Qds radii we use a normal distribution with 1.75 nm mean value and 0.3 nm standard deviation. A top view scheme of the system is shown in the inset of Figure 6(b).Figure 6
(a) ObtainedI
(
V
) curves (electron/hole and total currents) for two layers of 10 Si/SiO2 Qds. The system arrangement is described in the text. The Qds radius distribution is also shown in the inset. (b) Total accumulated charge (∑
i
‍
Δ
N
i where i
=
1
⋯
20) of the structure as a function of the external bias voltage. In the inset, a top view of the system is presented.
(a)
(b)In Figure6(a), the total I
(
V
) curve and the electron and hole currents terms are shown. From the I
(
V
) curves, two different regimes are obtained: for low and medium voltage ranges, the current reflects the discrete nature of the Qds energy spectra. However, for the largest voltages, the continuous part of the Qd DOS begins to be conductive and the current increases in a continuous form loosing its step-like shape. Concerning the electron/hole current components, the electron term dominates since its potential barrier is lower than the hole one. Moreover, the stepping behavior is still present in the current reflecting the aperture of the different conductive channels. In this case, the transport involves several tunneling processes between the Qds. Therefore, the energy level alignment is necessary. However, the negative differential resistance (NDR) is not clearly visible since there are many electron/hole pathways and the sum of the different conduction channels masks this effect [22].Regarding the accumulated charge plotted in Figure6(b), it is the total sum of the accumulated charge in each Qd (∑
i
‍
Δ
N
i where i
=
1
⋯
20). The geometrical disposition of the Qds strongly affects the final response of the system as we have shown in the previous section. For this system, since there are many electronic conduction channels and a complex capacitive interaction between the Qds, it is hard to explain clearly the accumulated charge trend. Reflecting that, the geometrical disposition of the Qds influences directly the final response of the system.


## 4. Conclusions

The high efficiency concepts of the next generation of Qds based devices pose new requirements on models for the theoretical description of their transport properties. An intuitive theoretical framework suitable for this purpose is available in the noncoherent rate equations. This approach provides a simple and transparent method to describe the electron transport. Using the Transfer Hamiltonian approach to describe the tunneling current terms in ballistic regime, the rate equations can be used in order to obtain the nonequilibrium distribution functions in each Qd. The effect of self-charge has been taken into account solving the Poisson equation with the appropriate boundary conditions for each Qd that involve the capacity coupling between the different parts of the system and the accumulated charge in each Qd. As expected, the calculation of the local potential inside each Qd is one of the most critical points, since theI-V curves depend on the position of the energy level. Due to the simplicity of the model, this can be easily extended to analyze arbitrary large arrays of Qds of interest in technological applications creating a powerful and useful tool that enables to design new concept devices based on Qd properties.In order to simulate devices as realistic as possible, suitable expressions for the transmission coefficients, the energy level positions and the capacitive coupling have been used. These parameters can be described using basic material properties and geometrical representations of the system. Moreover, the hole currents have been taken into account, obtaining a complete description of the electron transport in the structure.Finally, the proposed formalism has been used to describe, using realistic material parameters, the electrical response of a single Si/SiO2 Qd. Moreover, a bilayer structure has also been simulated. This structure composes the basic building block for future devices based on Qds and demonstrates the practicability of the here presented approach.



---

*Source: 426541-2015-03-24.xml*